平面单元的积分方法
使用下标 \(i,j,k\) 表示空间维度的指标; 使用上标 \(a\) 表示单元节点的指标; 使用上标 \(i\) 表示单元积分点的指标
三角形单元
高阶三角形单元的形函数
一阶的三角形单元形函数为 \[ N^1 = L_1, \quad N^2 = L_2, \quad N^3 = L_3 \] 三角形单元随着阶次增加, 其所覆盖的完备的多项式阶数也恰好逐次增加, 不会像四边形单元, 增加节点数也会增加不完备的多项式项. 更高阶的三角形单元形函数可以通过 Lagrange 函数进行构造, 其构造方法为 \[ N^{a} = l_{I}^{I}(L_1) l_{J}^{J}(L_2) l_{K}^{K}(L_3) \] 式中, \(I+J+K = M\), \(M\) 为三角形单元的阶数. \(l_{k}^{n}(\xi)\) 为 Lagrange 函数, 定义为 \[ l_{k}^{n}(\xi) \triangleq \frac{(\xi - \xi_0)(\xi - \xi_1) \cdots (\xi - \xi_{k-1})(\xi - \xi_{k+1}) \cdots (\xi - \xi_n)} {(\xi_k - \xi_0)(\xi_k - \xi_1) \cdots (\xi_k - \xi_{k-1})(\xi_k - \xi_{k+1}) \cdots (\xi_k - \xi_n)} \] 如此构造的 \(N^a\) 是满足性质 \(N^a(L^{b}) = \delta_{ab}\) 的. 对于节点 \(a\), \[ l_{I}^{I}(L_1^a)=1,\quad l_{J}^{J}(L_2^a)=1,\quad l_{K}^{K}(L_3^a)=1 \] 对于其它节点 \(a'\), 其索引为 \((I', J', K')\), 总存在一个索引分量 (不妨选取 \(I'\)), 使得 \(I' < I\), 这样就得到 \[ l_{I}^{I}(L_1^{a'}) = 0 \Rightarrow N^a(L^{a'}) = 0 \] 对于二阶单元, 备选的 Lagrange 函数为 \[ l_0^0(\xi) = 1, \quad l_1^1(\xi) = 2\xi, ~\{0,~\frac{1}{2}\}, \quad l_2^2(\xi) = \xi(2\xi-1), ~\{0,~\frac{1}{2},~1\} \] 花括号中是 Lagrange 函数的插值点. 因此, 二阶三角形单元的形函数为
角点处的节点: \[ N^a = L_a (2L_a - 1), \quad a=1,2,3 \] 边上的节点: \[ N^4 = 4 L_1 L_2, \quad N^5 = 4 L_2 L_3, \quad N^S6 = 4 L_3 L_1 \]
三角形单元内的坐标可以使用面积坐标进行表示 \[ x(L_1, L_2, L_3) = x_1 L_1 + x_2 L_2 + x_3 L_3 \\ y(L_1, L_2, L_3) = y_1 L_1 + y_2 L_2 + y_3 L_3 \] 式中, \((x_i, y_i)\) 是三角形节点在笛卡尔全局坐标系下的坐标, 面积坐标 \(L_i\) 在 \([0,1]\) 内变化, 并满足如下关系: \[ L_{i}(x^a, y^a) = \delta_{ia}, \quad \sum_{i} L_i = 1 \] 在有限元中, 需要计算单元 \(\Omega^e\) 上形如 \[ I = \iint_{\Omega^e} f(x,y) ~\mathrm{d}x \mathrm{d}y \] 的积分, 而面积坐标提供了一种坐标变换的方式, 使得在任意形状三角形区域的积分可以转化到标准区域进行. 这对数值积分非常友好, 因为 Gauss 型数值积分的精度通常是在标准区域上进行衡量的. 进行积分变换, 第一步是要计算 Jacobian 矩阵, 衡量微元之间的仿射变换关系: \[ \begin{pmatrix} \mathrm{d}x \\ \mathrm{d}y \end{pmatrix} = \frac{\partial (x,y)}{\partial(L_1,L_2)} \begin{pmatrix} \mathrm{d}L_1 \\ \mathrm{d}L_2 \end{pmatrix} \] 对于实体单元, 微元之间的系数为 Jacobian 矩阵的行列式值: \[ \mathrm{d}x \mathrm{d}y = \det \frac{\partial (x,y)}{\partial(L_1,L_2)} \mathrm{d}L_1 \mathrm{d}L_2 \] 对于三节点三角形单元. 行列式值为常值, 等于 2 倍的单元面积: \[ J \triangleq \det \frac{\partial (x,y)}{\partial(L_1,L_2)} = 2S_{\Omega^e} \] 经过坐标变换, 就得到原积分 \(I\) 为 \[ I = J \iint_{\triangle} \tilde f(L_1, L_2) ~\mathrm{d}L_1 \mathrm{d}L_2 \tag{1} \] 式中, \(\triangle\) 表示标准单元区域, 在文中设置的面积坐标下, \(\triangle\) 为 \((0,0)\) \((1,0)\) \((0,1)\) 三节点连线围成的直角三角形区域. 在标准区域 \(\triangle\) 中, 考虑如下积分:
\[
\iint_{\triangle} \tilde f(L_1, L_2) ~\mathrm{d}L_1 \mathrm{d}L_2
\] 数值积分将积分近似为单元区域某些点的取值和权重的乘积,
并进行求和: \[
\iint_{\triangle} \tilde f(L_1, L_2) ~\mathrm{d}L_1 \mathrm{d}L_2
\approx \sum_i \tilde f(\xi^i) w^i
\] 为了方便验证求积公式的代数精度,
给出多项式函数在标准区域的积分解析公式: \[
\iint_{\triangle} L_1^a L_2^b L_3^c ~\mathrm{d}L_1 \mathrm{d}L_2
= \frac{a!b!c!}{(a+b+c+2)!}
\] 一种常用的三角形区域的数值积分方式为 Hammer 积分,
这里列举选用一个积分点和三个积分点的情景. 只选取一个积分点时 \[
\xi^{1} = (\frac{1}{3}, \frac{1}{3}), \quad w^{1} = \frac{1}{2}
\] 可以验证, 这种求积公式是满足一阶代数精度的. 选取三个积分点时,
有两种可能: 第一种求积公式为 \[
\xi^{1} = (0, \frac{1}{2}), \quad
\xi^{2} = (\frac{1}{2}, 0), \quad
\xi^{3} = (\frac{1}{2}, \frac{1}{2}) \\
w^{1} = w^{2} = w^{3} = \frac{1}{6}
\]
第二种求积公式为
现在将标准区域的求积公式代入到原积分公式 (1) 当中得到 \[ I \approx J \sum_i \tilde f(\xi^i) w^i \] 令 \(\widetilde{w}^i \triangleq J w^i\), 就得到 \[ I \approx \sum_i \tilde f(\xi^i) \widetilde{w}^i \]
[!NOTE]
一些教材 (Zienkiewicz) 在处理三角形单元的数值积分时, 将 Jacobian 行列式写成 2 倍的三角形面积, 并将系数 2 并入权重当中: \[ I \approx S_{\Omega^e} \sum_i \tilde f(\xi^i) \hat{w}^i, \quad \hat{w}^i = 2 w^i \] 由此得到的一阶 Hammer 积分: \[ \xi^{1} = (\frac{1}{3}, \frac{1}{3}), \quad w^{1} = 1 \] 和二阶 Hammer 积分 \[ \xi^{1} = (0, \frac{1}{2}), \quad \xi^{2} = (\frac{1}{2}, 0), \quad \xi^{3} = (\frac{1}{2}, \frac{1}{2}) \\ w^{1} = w^{2} = w^{3} = \frac{1}{3} \]
四边形单元
对于四边形四节点单元, 使用的形函数为 \[ N^{I}(\xi, \eta) = \frac{1}{4} (1 + \xi^I \xi) (1 + \eta^I \eta) \] 式中, \((\xi^I, \eta^I)\) 为节点 \(I\) 的参数坐标. 等参变换为 \[ x(\xi, \eta) = \sum_{I=1}^{4} x^I N^{I}(\xi, \eta), \quad y(\xi, \eta) = \sum_{I=1}^{4} y^I N^{I}(\xi, \eta) \] 将积分转换到参数坐标区域中, 对应的 Jacobian 为 \[ J(\xi, \eta) = \frac{1}{4} \begin{bmatrix} \sum\limits_{I=1}^{4} x^I \xi^I (1 + \eta^I \eta) & \sum\limits_{I=1}^{4} x^I \eta^I (1 + \xi^I \xi) \\ \sum\limits_{I=1}^{4} y^I \xi^I (1 + \eta^I \eta) & \sum\limits_{I=1}^{4} y^I \eta^I (1 + \xi^I \xi) \end{bmatrix} \] 因此, 对于原积分为 \[ I = \iint_{\Omega^e} f(x,y) ~\mathrm{d}x \mathrm{d}y \] 将积分区域转换到参数坐标区域, 得到 \[ I = \int_{ -1}^{1}\int_{-1}^{ 1} \tilde f(\xi, \eta) J(\xi, \eta) ~\mathrm{d}\xi \mathrm{d}\eta \] 对上述乘积区域的积分 (\((-1,1) \times (-1,1)\)), 可以使用 Gauss 积分方式进行计算. 首先, 一维情景下, 在区间 \((-1,1)\) 上的 Gauss 积分列式为 \[ \int_{-1}^{1} f(\xi) ~\mathrm{d}\xi \approx \sum_{i=1}^{N_{\rm ipt}} w^{i} f(\xi^{i}) \] 式中, \(N_{\rm ipt}\) 是积分点的数量. 常用的 Gauss 积分点和权重如下表所示
在给出一维情景的 Gauss 积分列式之后, 二维或更高维的情景可以通过一维的例子方便地得到. 对于维度为 \(N_{\rm d}\) 的情景, 在乘积空间 \(\Pi(-1,1)\) 中, Gauss 积分的列式为 \[ \int_{\scriptsize \Pi(-1,1)} f(\xi_1, \xi_{2}, \ldots, \xi_{N^{\rm d}}) \mathrm{d} \xi \approx \underbrace{\sum_{i=1}^{N_{\rm ipt}} \sum_{j=1}^{N_{\rm ipt}} \cdots \sum_{k=1}^{N_{\rm ipt}}}_{N^{\rm d}} f(\xi_1^i, \xi_{2}^j, \ldots, \xi_{N^{\rm d}}^k) \underbrace{w^{i} w^{j} \cdots w^{k}}_{N^{\rm d}} \] 例如, 二维情景的 Gauss 积分为 \[ I = \int_{ -1}^{1}\int_{-1}^{ 1} \tilde f(\xi, \eta) J(\xi, \eta) ~\mathrm{d}\xi \mathrm{d}\eta \approx \sum_{i=1}^{N_{\rm ipt}} \sum_{j=1}^{N_{\rm ipt}} \tilde f(\xi^{i}, \eta^{j}) J(\xi^{i}, \eta^{j}) w^i w^{j} \] 当 \(N_{\rm ipt} = 1\) 时, 积分公式为 \[ I \approx 4\tilde f(0, 0) J(0, 0) \]
[!NOTE]
四边形单元 Jacobian 在参数坐标点 \((0,0)\) 取值为 \[ J(0,0) = \frac{1}{4}S_{\Omega^e} \]