第6章 三维单元
6.1 背景
到目前为止，读者应该已经很好地理解了有限元方法用于解决一维和二维问题。第2章和第3章中制定的原则同样适用于2-D单元和1-D单元。在第4章到第6章中，我们看到了有限元方法利用二维三角形或四边形单元，其准确性各不相同。将该方法扩展到三维是合乎逻辑的。与2-D问题相比，建立计算域和边界条件所需的数据量显著增加；同样，计算工作量也大幅度增加。由于与3-D单元相关的节点数量增加（以及相关的带宽），计算网格通常受到限制。在本章中，开发并应用了四面体和六面体单元到熟悉的热传导方程。解决了一个边界受到辐射热传递和对流影响的时变热传导问题。
6.2 单元网格
将问题域离散化为由单元和节点点组成的计算域没有严格的理论基础。为了产生一个好的网格；即使严格遵守有限元的求解步骤，一些不恰当的操作仍会使结果不准确。 在1-D域中，将离散化为单元阵列并不难理解——在变量快速变化的地方通常需要小的单元尺寸。在2-D区域中，生成一个可接受的网格变得更加麻烦；可能需要制作几个网格，才能创建出一个产生合理结果的合适网格。在3-D域中，网格不仅更难制作，而且难以可视化——如果问题域复杂，可能需要许多试验网格。应该认识到，生成一个好的网格是一种艺术。计算机科学领域的计算几何技术使人们能够基于所需数量的节点和单元产生一个最优网格（见Pepper和Gewali，2002）。在生成网格上可以花费相当多的时间；然而，一旦网格制作完成，有限元方法允许在不增加任何额外努力的情况下改变边界条件，具有很大的灵活性。
6.3 形状函数
6.3.1 四面体
用于定义线性、3-D四面体的四个节点。四面体如图3所示。插值函数写为 \[ \begin{equation} \phi = \alpha_1 + \alpha_2 x + \alpha_3 y + \alpha_4 z \end{equation} \] 对于位于四面体顶点的四个节点，产生一组四个方程 $$ \[\begin{equation} \begin{matrix} \phi_1 = \alpha_{1} + \alpha_{2} x_1 + \alpha_{3} y_1 + \alpha_{4} z_1 \\ \phi_2 = \alpha_{1} + \alpha_{2} x_2 + \alpha_{3} y_2 + \alpha_{4} z_2 \\ \phi_3 = \alpha_{1} + \alpha_{2} x_3 + \alpha_{3} y_3 + \alpha_{4} z_3 \\ \phi_4 = \alpha_{1} + \alpha_{2} x_4 + \alpha_{3} y_4 + \alpha_{4} z_4 \end{matrix} \end{equation}\] $$
可以写成矩阵形式为 \[ \begin{equation} \phi=\boldsymbol{C}\alpha \end{equation}\] 或者
\[\begin{equation} \begin{bmatrix} \phi_1 \\ \phi_2 \\ \phi_3 \\ \phi_4 \end{bmatrix} = \begin{bmatrix} 1 & x_1 & y_1 & z_1\\ 1 & x_2 & y_2 & z_2\\ 1 & x_3 & y_3 & z_3\\ 1 & x_4 & y_4 & z_4 \end{bmatrix} \begin{bmatrix} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \end{bmatrix}\end{equation}\]
取\(\boldsymbol{C}\)的逆，得到\(\alpha\)列向量为 \[ \begin{equation} \alpha=\boldsymbol{C}^{-1}\phi \end{equation}\] \[ \begin{equation} \phi=\begin{bmatrix} 1&x&y&z\end{bmatrix}= \begin{bmatrix} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \end{bmatrix}=\begin{bmatrix} 1&x&y&z\end{bmatrix}\boldsymbol{C}^{-1}\phi \end{equation} \]
\(\begin{bmatrix}1&x&y&z\end{bmatrix}\)实际上是形状函数矩阵，即 \[\begin{equation} \phi \equiv \boldsymbol{N} \end{equation}\] 其中 \[ \begin{equation}N = \begin{bmatrix} 1 & x & y & z \end{bmatrix}\boldsymbol{C}^{-1} \end{equation}\]
如果希望使用二次四面体，如图4所示，则需要10个节点来定义单元。单元节点位置直接与二维空间的帕斯卡三角形相关，并且在三维中它们以金字塔形式展示类似的属性（见Zienkiewicz和Taylor，1989）。
\(\phi\)的表达式变为 \[\begin{equation} \phi = \alpha_1 + \alpha_2 x + \alpha_3 y + \alpha_4 z + \alpha_5 x^2 + \alpha_6 xy + \alpha_7 y^2 + \alpha_8 yz + \alpha_9 z^2 + \alpha_{10} xz \end{equation}\] 或者 \[\begin{equation} \phi = \begin{bmatrix} 1 & x & y & z & x^2 & xy & y^2 & yz & z^2 & xz \end{bmatrix} \boldsymbol{C}^{-1}\phi\end{equation}\] 其中 \[\begin{equation}\boldsymbol{C} = \begin{bmatrix} 1 & x_1 & y_1 & z_1 & x_1^2 & x_1y_1 & y_1^2 & z_1^2 & x_1z_1\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & x_{10} & y_{10} & z_{10} & x_{10}^2 & x_{10}y_{10} & y_{10}^2 & z_{10}^2 & x_{10}z_{10} \end{bmatrix}\end{equation} \] 四面体可以建立体积坐标系统，就像二维三角形的面积坐标系统一样。在3-D情况下，使用四个距离比率，每个比率都垂直于一个边：\(L_1\)、\(L_2\)、\(L_3\)和\(L_4\)。图5中显示了体积\(V_1\)。
![图5 体积坐标系统 [L_i = V_i/V, i = 1, \cdots, 4]](/2024/09/19/%E6%9C%89%E9%99%90%E5%85%83%E7%90%86%E8%AE%BA%E6%8E%A8%E5%AF%BC%EF%BC%88%E5%85%AD%EF%BC%89/5.png)
更高阶单元的形状函数与线性四面体类似；在这种情况下，不是直线，而是常数\(L\)平面平行于相应节点的对面。例如，图6中节点1的二次四面体形状函数为 \[\begin{equation} N_1 = (2L_1-1)L_1 \end{equation}\] 其中节点\(2\)、\(3\)、\(4\)、\(6\)、\(9\)和\(10\)处\(L_1=0\)。对于\(L_1=1/2\)，平面通过节点\(5\)、\(7\)和\(8\)。在节点\(5\)处， \[ \begin{equation} N_5 = 4L_1L_2 \end{equation}\] 其中\(L_2=0\)通过节点\(1\)、\(3\)、\(4\)、\(7\)、\(8\)和\(10\)。额外的形状函数关系可以通过参考给出二维二次三角形单元形状函数的方程进行分析。
7.3.2 六面体
四面体是一个简单的单元，它直接从三角形单元扩展而来。然而，在复杂的三维网格中，四面体难以可视化，并且可能需要相当大的努力来定义单元的连接性。使用六面体单元可以获得更精确的解，虽然需要更多的节点，但所需的努力较少。这个家族中最简单的单元是八个节点的三线性六面体，如图7所示。插值函数为 \[\begin{equation} \phi = \alpha_1 + \alpha_2 x + \alpha_3 y + \alpha_4 z + \alpha_5 xy + \alpha_6 xz + \alpha_7 yz + \alpha_8 xyz \end{equation}\]
方程17可以使用方程3的矩阵形式表示为 \[\begin{equation} \begin{bmatrix} \phi_1 \\ \vdots \\ \phi_8 \end{bmatrix} = \begin{bmatrix} 1 & x_1 & y_1 & z_1 & x_1y_1 & x_1z_1 & y_1z_1 & x_1y_1z_1\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots\\ 1 & x_8 & y_8 & z_8 & x_8y_8 & x_8z_8 & y_8z_8 & x_8y_8z_8\\ \end{bmatrix} \begin{bmatrix} \alpha_1 \\ \vdots \\ \alpha_8 \end{bmatrix} \end{equation}\]
由于\(\alpha = \boldsymbol{C}^{-1}\phi\)，我们可以从表达式 \[\begin{equation} \phi = \begin{bmatrix} 1 & x & y & z & xy & xz & yz & xyz \end{bmatrix}\boldsymbol{C}^{-1}\phi \end{equation}\] 中找到\(\phi\)的值，如之前一样得到形状函数： \[\begin{equation} N = \begin{bmatrix} 1 & x & y & z & xy & xz & yz & xyz \end{bmatrix}\boldsymbol{C}^{-1} \end{equation}\]
就像二维四边形单元一样，也可以为六面体设计一个自然坐标系统，如图8所示。形状函数定义为（经过一些代数运算！） \[\begin{equation} \begin{bmatrix} N_1\\N_2\\N_3\\N_4\\N_5\\N_6\\N_7\\N_8 \end{bmatrix} =\frac{1}{8} \begin{bmatrix} (1-\xi)(1-\eta)(1-\zeta)\\ (1+\xi)(1-\eta)(1-\zeta)\\ (1+\xi)(1+\eta)(1-\zeta)\\ (1-\xi)(1+\eta)(1-\zeta)\\ (1-\xi)(1-\eta)(1+\zeta)\\ (1+\xi)(1-\eta)(1+\zeta)\\ (1+\xi)(1+\eta)(1+\zeta)\\ (1-\xi)(1+\eta)(1+\zeta) \end{bmatrix} \end{equation}\]
上式可以写成更简化的方式 \[ \begin{equation} N_i=\frac{1}{8}(1+\xi\xi_i)(1+\eta\eta_i)(1+\zeta\zeta_i),i=1,2,\cdots,8 \end{equation}\] 式中\(\xi_i\)，\(\eta_i\)及\(\zeta_i=\pm 1\)，取决于节点的位置。
将 8 节点二次单元扩展到 20 节点二次六面体可以通过在每个单元的边上添加一个中点节点来轻松实现。图9显示了三维二次六面体。角节点的形状函数由以下关系获得 \[ \begin{equation} N_i=\frac{1}{8}(1+\xi\xi_i)(1+\eta\eta_i)(1+\zeta\zeta_i)(\xi\xi_i+\eta\eta_i+\zeta\zeta_i-2),i=1,2,\cdots,8 \end{equation}\] 对于中边节点， \[ \begin{equation} N_i=\frac{1}{4}(1-\xi^2)(1+\eta\eta_i)(1+\zeta\zeta_i),i=10,12,18,20 \end{equation}\] \[ \begin{equation} N_i=\frac{1}{4}(1-\eta^2)(1+\xi\xi_i)(1+\zeta\zeta_i),i=9,11,17,19 \end{equation}\] \[ \begin{equation} N_i=\frac{1}{4}(1-\zeta^2)(1+\xi\xi_i)(1+\eta\eta_i),i=13,14,15,16 \end{equation}\]
9 节点双二次单元扩展到 27 节点三二次单元，每个六个面的中心都有一个节点，单元中心也有一个。
注意，当其中一个自然坐标为\(\pm 1\)时，三维单元简化为它们的二维对应物。同样，当两个坐标等于\(\pm 1\)时，形状函数简化为一维形状函数，允许三维单元连接到线单元。单元之间的连续性存在于三维单元的单元面（面积）上；节点值在相邻单元共有的单元面上以相同的方式描述未知变量（\(\phi\)）的变化。
6.4 数值积分
在这一点上，评估涉及多个节点的形状函数的积分的困难应该是显而易见的。使用体积坐标，四面体单元可以精确地使用方程14进行评估；然而，与内积乘法相关的代数量很快就会变得非常繁琐。
在二维三角形单元中出现的面积积分形式为（来自第4章）\(I = \int_{0}^{1} \int_{0}^{1-L_2}f(L_1, L_2, L_3)|\boldsymbol{J}dL_1dL_2\)，其中插值函数以面积坐标表示。由于雅可比是面积坐标的函数，雅可比逆的显式形式几乎是不可能的。虽然可以使用方程4.45（二维形式的方程14）来评估一些积分，但随着项数的急剧增加，出错的机会也变得相当大。最好让计算机进行积分，即 \[\begin{equation} \int_{0}^{1} \int_{0}^{1-L_2}f(L_1, L_2, L_3)|\boldsymbol{J}|dL_1dL_2= \frac{1}{2}\sum_{i=1}^{n} w_i g\left[(L_1)_i, (L_2)_i, (L_3)_i\right] \end{equation}\] 其中积分的阶数（\(n\)）是通过累加坐标\(L_1\)、\(L_2\)和\(L_3\)的幂来确定的。以类似的方式，四面体单元被评估为 \[\begin{equation} \int_{0}^{1} \int_{0}^{1-L_1}\int_{0}^{1-L_1-L_2}f(L_1, L_2, L_3,L_4)|\boldsymbol{J}|dL_1dL_2dL_3 =\frac{1}{6} \sum_{i=1}^{n} w_i g\left[(L_1)_i, (L_2)_i, (L_3)_i,(L_4)_i\right] \end{equation}\] 其中四面体的数值积分采样点在表1中显示。如果我们希望积分\(L_1L_2L_4\)，指数的总和为3，将需要一个四次积分方案（\(n=4\)个点）。评估六面体单元矩阵的过程与用于二维四边形的过程几乎相同。雅可比矩阵现在变为\(3\times 3\)，而不是四边形的\(2\times 2\)。积分的一般形式为 \[\begin{equation} \int_{-1}^1 \int_{-1}^1 \int_{-1}^1 f(\xi, \eta, \zeta)|\boldsymbol{J}|d\xi d\eta d\zeta= \sum_{i=1}^{n} \sum_{j=1}^{n} \sum_{k=1}^{n} w_i w_j w_k f(\xi_i, \eta_j, \zeta_k) |J(\xi_i, \eta_j, \zeta_k)| \end{equation} \] 其中\(w_i, w_j, w_k\)是与相应的高斯点相关的权重因子。注意方程29是通过结合最初为一维单元（第3章）建立的一维公式得到的。权重因子和积分点坐标在表2中列出。
图10显示了立方积分方案的高斯点位置。请注意，现在需要评估刚度矩阵\(\boldsymbol{K}\)的积分点数量，对于三线性单元变为8个，对于三次单元变为27个。如果一个人想使用三阶六面体（需要64个节点），将需要64个积分点。尽管手工评估矩阵的工作量巨大，但计算机可以轻松快速地执行此类操作。
