4 二维四边形单元
4.1 背景
四边形单元是一个四边形。最常用的形状包含位于顶点的四个节点，并且是矩形配置中的双线性。其他常用的矩形单元涉及8节点二次、9节点双二次和12节点三次配置。最初，有限元方法仅使用三角形单元。然而，现在许多研究人员更倾向于使用四边形单元。在双线性近似中，四边形单元比三角形多一个节点，后者只有三个节点。此外，梯度是坐标方向的线性函数，而三角形单元中的梯度是常数。对于更高阶的单元，实现了更复杂的表示，从近似的角度来看越来越准确。然而，必须注意评估提高准确性的好处与更复杂单元相关的增加计算成本之间的权衡。
4.2 单元网格
与二维(2-D)三角形网格一样，四边形网格将一个区域细分为一组小域，即单元，只是现在这些是四边形。图1显示了一个矩形计算域被离散化为一组四边形单元。
请注意，如果我们连接一对对角线节点，我们得到两个三角形单元。双线性四边形单元包含四个局部节点，按逆时针顺序被指定为1、2、3和4。图2显示了一个包含四个局部节点的一般四边形单元。
四边形网格更接近标准的二维有限差分网格，而不是三角形网格。然而，四边形在离散化一个区域方面与三角形一样灵活。有限差分网格与四边形网格之间的不同之处在于，有限差分网格必须是正交的，即所有线交点（节点位置）必须彼此成直角；在有限元四边形中，每个单元都是独特的——每个单元的面可以有不同的大小和坡度。换句话说，有限差分网格是全局形成的，并且必须在所有表面和交点处垂直；有限元网格是局部形成的，不必是正交的（在物理空间中）——然而，它必须在计算域内几何上兼容以确保单元之间的连续性。
在这种情况下，连续性意味着插值函数必须在单元之间连续。如果控制方程包含二阶导数（例如，热扩散），则近似函数必须在单元之间连续。如果我们希望积分\(d^n\phi/dx^n\)，那么\(\phi\)必须具有连续的\((n-1)\)阶导数，以便只在nth导数中存在有限的跳跃不连续。Carey和Oden
(1983)，Cary (1997)，和Heinrich和Pepper
(1999)对有限元和三角形及四边形单元的这种要求及其影响进行了更详细的讨论。尽管需要额外的节点来定义四边形，但对于相同数量的节点点，单元的数量减少到三角形单元网格的一半。在离散化几何域时，通常首先将区域划分为四边形。如果需要三角形单元，四边形很容易细分。在几乎所有情况下，由四边形单元组成的网格就足够了，通常比由三角形单元组成的网格更准确。四边形网格扩展到三维是容易获得的，即一个盒子，并且比一组3-D四面体更容易在概念上可视化。
4.3 形状函数
与二维三角形单元一样，二维四边形单元基于线性、二次、三次或更高阶近似。线性三角形单元被归类为简单单元（包含常数和线性项的近似），而二次三角形为复杂单元（包含常数、线性和二次项）。简单和复杂单元的单元边界不需要单元边与坐标轴平行。然而，因为它是一个多路单元，为了可能的简单多项式表示，四边形单元要求其边与坐标系统平行。正如我们稍后将看到的，当边不平行于全局轴系统时，这个要求很容易通过局部坐标变换来放松。
4.3.1 双线性矩形单元
在其最简单的形式中，四边形成为一个矩形单元，其边界与坐标系统平行（在物理空间中）。
通过使用局部自然坐标系统进一步扩展单元，可以得到一个泛化的四边形，其中矩形是其子集。我们从矩形单元开始。四节点双线性矩形的插值函数给出为 \[ \begin{equation} f = a_1 + a_2 x + a_3 y + a_4 xy \end{equation}\] 选择\(xy\)项而不是\(x^2\)或\(y^2\)，因为\(xy\)不破坏对称性，并假设沿\(x\)或\(y\)的\(\phi\)呈线性变化。对于常数\(x\)的线，单元在\(y\)方向上是线性的，反之亦然；对于常数\(y\)的线，单元在\(x\)方向上是线性的，因此得名双线性单元。\(\phi\)的梯度被视为一个坐标方向上的线性函数，即， \[\begin{equation} \frac{\partial \phi}{\partial x} = \alpha_2 + \alpha_4 y, \frac{\partial \phi}{\partial y} = \alpha_3 + \alpha_4 x \end{equation} \] 图3显示了矩形单元和节点配置。如果到坐标原点的距离，位于单元中点，用\(a\)和\(b\)表示，节点值给出为 \[\phi=\phi_1,x=-b,y=-a\] \[\phi=\phi_2,x=b,y=-a\] \[\phi=\phi_3,x=b,y=a\] \[\begin{equation}\phi=\phi_4,x=-b,y=a\end{equation}\]
将这些值代入方程(1)中，对于\(\alpha_is，i = 1, 4\)，得到 \[\alpha_1=\frac{\phi_1+\phi_2+\phi_3+\phi_4}{4} \] \[\alpha_2=\frac{-\phi_1+\phi_2+\phi_3-\phi_4}{4b} \] \[\alpha_3=\frac{-\phi_1-\phi_2+\phi_3+\phi_4}{4a}\] \[\begin{equation}\alpha_4=\frac{\phi_1-\phi_2+\phi_3-\phi_4}{4ab}\end{equation}\] 如果我们现在用节点值和形状函数近似ϕ，就像我们在方程4.7中做的那样，我们得到 \[\begin{equation}\phi=N_1\phi_1+N_2\phi_2+N_3\phi_3+N_4\phi_4\end{equation}\] 其中形状函数给出为 \[N_1=\frac{1}{4ab}(b-x)(a-y)\] \[N_2=\frac{1}{4ab}(b+x)(a-y)\] \[N_3=\frac{1}{4ab}(b+x)(a+y)\] \[\begin{equation}N_4=\frac{1}{4ab}(b-x)(a+y)\end{equation}\] 我们在方程（6）中看到，形状函数可以用长度比，\(x/b\)和\(x/a\)来表示。例如， \[\begin{equation}N_1 =\frac{1}{4ab}(b-x)(a-y=\frac{1}{4}(\frac{a-y}{a})\frac{b-x}{b})\end{equation}\] 其中\(-1\le x/b\le 1\)且\(-1\le y/a\le 1\)。
4.3.2 二次矩形单元
我们也可以像对三角形单元那样开发二维高阶四边形单元。这样做有两种方法，可以得到两种不同的二次单元：一种在每个角落有一个节点，每个边的中点有一个节点；另一种有九个节点，第九个节点在单元的中心。
我们将在这里讨论八节点单元。四节点双线性和八节点二次四边形，以及九节点双二次单元，在工业中主要使用，并且被发现对大多数计算目的来说是足够的。八节点二次单元的插值多项式，如图4所示，由八项组成。未知变量\(\phi\)表示为 \[\begin{equation}\phi = \alpha_1 + \alpha_2 x + \alpha_3 y + \alpha_4 xy + \alpha_5 x^2 + \alpha_6 y^2 + \alpha_7 x^2y + \alpha_8 xy^2\end{equation}\]
如果节点值定义如图4所示。 \[\phi=\phi_1,x=-b,y=-a\] \[\phi=\phi_2,x=0,y=-a\] \[\phi=\phi_3,x=b,y=-a\] \[\phi=\phi_4,x=b,y=0\] \[\phi=\phi_5,x=b,y=a\] \[\phi=\phi_6,x=0,y=a\] \[\phi=\phi_7,x=-b,y=a\] \[\begin{equation}\phi=\phi_8,x=-b,y=0\end{equation}\] 从方程9代入方程8得到以下八个方程（用矩阵形式表示）： \[ \begin{equation}\begin{bmatrix} 1 & -b & -a & ab & b^2 & a^2 & -ab^2 & -a^b \\ 1 & 0 & -a & 0 & 0 & a^2 & 0 & 0 \\ 1 & b & -a & -ab & b^2 & a^2 & -ab^2 & a^2b \\ 1 & b & 0 & 0 & b^2 & 0 & 0 & 0 \\ 1 & b & a & ab & b^2 & a^2 & ab^2 & a^2b \\ 1 & 0 & a & 0 & 0 & a^2 & 0 & 0 \\ 1 & -b & a & -ab & b^2 & a^2 & ab^2 & -a^2b \\ 1 & -b & 0 & 0 & b^2 & 0 & 0 & 0 \\ \end{bmatrix} \begin{bmatrix} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \\ \alpha_5 \\ \alpha_6 \\ \alpha_7 \\ \alpha_8 \\ \end{bmatrix} = \begin{bmatrix} \phi_1 \\ \phi_2 \\ \phi_3 \\ \phi_4 \\ \phi_5 \\ \phi_6 \\ \phi_7 \\ \phi_8 \\ \end{bmatrix}\end{equation}\] 从中可以得到\(\alpha_is，i = 1, \cdots, 8\)。用节点值和形状函数近似ϕ得到 \[\begin{equation} \phi = N_1 \phi_1 + N_2 \phi_2 + N_3 \phi_3 + N_4 \phi_4 + N_5 \phi_5 + N_6 \phi_6 + N_7 \phi_7 + N_8 \phi_8 \end{equation} \] 形状函数为（参见Carey, 1997; Heinrich和Pepper, 1999; Zienkiewicz和Taylor, 1989） \[N_1=-\frac{1}{4(ab)^2}(b-x)(a-y)(ab+ax+by)\] \[N_2=\frac{1}{2ab^2}(b^2-x^2)(a-y)\] \[N_3=\frac{1}{4(ab)^2}(b-x)(a-y)(ay-by-ab)\] \[N_4=\frac{1}{2a^2b}(a^2-y^2)(b+x)\] \[N_5=\frac{1}{4(ab)^2}(b+x)(a+y)(ax+by-ab+yb)\] \[N_6=-\frac{1}{2ab^2}(b^2-x^2)(a+y)\] \[\begin{equation}N_7=-\frac{1}{4(ab)^2}(b-x)(a+y)(ab+xa-by)\end{equation}\] 于是\(\phi\)的梯度变为 \[\frac{\partial \phi}{\partial x}=\alpha_2+\alpha_4y+2\alpha_5x+2\alpha_7xy+\alpha_8y^2\] \[\begin{equation} \frac{\partial \phi}{\partial y}=\alpha_3+\alpha_4x+2\alpha_6y+\alpha_7x^2+2\alpha_8xy \end{equation}\] 注意到\(\partial\phi/\partial x\)在\(x\)方向上为线性，但在\(y\)方向上为二次，而\(\partial\phi/\partial y\)相反。就长度比例而言，\(N_1\)的形状函数可以重写为 \[\begin{equation}N_1=-\frac{1}{4}\Bigg(1-\frac{x}{b}\Bigg)\Bigg(1-\frac{y}{a}\Bigg)\Bigg(1+\frac{x}{a}+\frac{y}{a}\Bigg)\end{equation}\]
4.4 自然坐标系统
在大多数情况下，要建模的物理域并不包括直边、正交边界。具有曲线边界的区域可以使用具有曲线边的四边形单元进行离散化，以获得更准确的解决方案。通过用曲线坐标表示x，y坐标，从直线边转换为曲线边： \[ x = x(\xi, \eta) \] \[ \begin{equation}y = y(\xi, \eta) \end{equation} \] \(\xi,\eta\)的选择取决于单元几何形状。当坐标变量在\(-1\le\xi \le\eta\le 1\)范围内时，通常将\(\xi,\eta\)坐标系统称为“自然”坐标系统。对于矩形情况，我们可以找到一组无量纲坐标\(\xi,\eta\)，使得 \[ \xi = \frac{x}{b} \] \[ \begin{equation}\eta = \frac{y}{a} \end{equation}\]
其中\(\xi,\eta\)有效地用局部坐标替换了\(x\)和\(y\)全局坐标，用于单个单元；我们之前提到这些特定的比率称为长度比（图3和4）。因此，利用我们在方程16中的定义，我们将双线性形状函数（6）在自然坐标系统中重写为 \[N_1=\frac{1}{4}(1-\xi)(1-\eta)\] \[N_2=\frac{1}{4}(1+\xi)(1-\eta)\] \[N_3=\frac{1}{4}(1+\xi)(1+\eta)\] \[\begin{equation}N_4=\frac{1}{4}(1-\xi)(1+\eta)\end{equation}\] 我们看到局部坐标系统总是在-1和1之间变化，并且在制定导数项时利用这一事实。图5显示了以单元为中心的\(\xi,\eta\)坐标系统。注意，这个新坐标系统不再重要，它必须是正交的——它也可以是曲线的。然而，不能直接获得关于\(x\)和\(y\)的导数，并且需要额外的转换。我们将在处理等参单元时说明这一点；现在，我们假设一个简单的笛卡尔坐标系统。
对于导数项\(\partial N_i/\partial x\)和\(\partial N_i/\partial y\)的评估类似于第三篇文章中的三角形单元。雅可比矩阵通过 \[\begin{equation} \begin{bmatrix} \frac{\partial N_i}{\partial \xi} \\ \frac{\partial N_i}{\partial \eta} \\ \end{bmatrix} =\begin{bmatrix} \frac{\partial x}{\partial \xi}&\frac{\partial y}{\partial \xi}\\ \frac{\partial x}{\partial \eta}&\frac{\partial y}{\partial \eta} \end{bmatrix} =\boldsymbol{J} \begin{bmatrix} \frac{\partial N_i}{\partial x} \\ \frac{\partial N_i}{\partial y} \\ \end{bmatrix}\end{equation}\] 其中雅可比矩阵的逆，\(\boldsymbol{J}^{-1}\)，可以用来找到关于\(x\)和\(y\)的导数值以及变换的行列式。
4.5 使用高斯求积法进行数值积分
将积分表述为ξ和η的形式可以得到简单的积分限。这在我们处理曲线边界或不平行于坐标轴的边界时特别有利。积分定义为 \[\begin{equation}\int^a_{-a} \int^b_{-b} F(x, y) dxdy== \int^1_{-1} \int^1_{-1} f(\xi, \eta)|\boldsymbol{J}|d\xi \, d\eta \quad (5.35)\end{equation}\] 导热矩阵因此变为 \[ \begin{equation}\boldsymbol{K}^{(e)} = \int \int_{\Omega^{(e)}}\boldsymbol{K} \left[ \frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial x} + \frac{\partial N_i}{\partial y}\frac{\partial N_j}{\partial y} \right] dxdy=\\ \int^1_{-1} \int^1_{-1} \left[] \frac{\partial N_i}{\partial \xi}\frac{\partial N_j}{\partial \xi} + \frac{\partial N_i}{\partial \eta}\frac{\partial N_j}{\partial \eta} \right)|J|d\xi d\eta \end{equation}\] 其中导数\(\frac{\partial N_i}{\partial x}\)和\(\frac{\partial N_i}{\partial y}\)可以作为关于\(\xi,\eta\),\(\frac{\partial N_i}{\partial \xi}\)和\(\frac{\partial N_i}{\partial \eta}\)的函数，从方程18中计算得到。所有包含\(N_iN_j\)乘积的面积积分项变为 \[\begin{equation} \int\int_{\Omega^{(e)}} N_i N_j dxdy = \int^1_{-1} \int^1_{-1} N_i N_j |\boldsymbol{J}|d\xi d\eta \end{equation}\] 积分中的右侧项不再涉及简单多项式，并且由于导数表达式中出现的\(1/|\boldsymbol{J}|\)项，它们可能变得麻烦。最常见的做法是使用称为高斯求积法的数值积分过程。该过程易于编程并且可以为积分提供几乎精确的值。高斯求积法用一组点上的加权和近似定积分。例如，在一维中， \[\begin{equation} \int^1_{-1} f(\xi) d\xi\approxeq\sum^N_{i=1} W_i f(\xi_i) \end{equation}\] 其中\(N\)是积分点数\(W_i\)是权重因子,\(\xi_i\)是高斯点。
使用\(N\)个高斯点，我们可以精确积分\(2N-1\)阶的多项式。因此，如果\(f(\xi)\)是一个三次多项式\(f(\xi) = a + b\xi + c\xi^2 + d\xi^3\)，那么选择\(N = 2，我们可以精确积分它。\xi_i\)值定义为（Heinrich和Pepper, 1999） \[ \begin{equation} \xi_1 = -\frac{1}{\sqrt{3}}, \xi_2 = \frac{1}{\sqrt{3}} \end{equation} \] 权重是\(W_1 = W_2 = 1\)，函数将被精确积分。对于方程37中的双重积分， \[ \begin{equation} \int^1_{-1} \int^1_{-1} f(\xi, \eta)d\xi d\eta \approxeq \sum_{i=1}^N \sum_{j=1}^N W_i W_j f(\xi_i, \eta_j) \end{equation}\] 高斯点定义为 \[ \begin{equation} \xi_i = \mp \frac{1}{\sqrt{3}}, \eta_j = \mp \frac{1}{\sqrt{3}} \end{equation} \] 对于图 所示的四节点四边形单元，使用四个高斯点。当使用八节点二次或九节点双二次单元时，方程5.36中的积分数值上得到 \[ \int \int f(\xi, \eta) \, W_i W_j \, d\xi \, d\eta \approx \sum_{i=1}^3 \sum_{j=1}^3 W_i W_j f(\xi_i, \eta_j) \quad (5.42) \] 九个积分点在图7中描述。点和相应的权重在表5.2中给出。
表1 \(N = 3\)的高斯点和权重
| 高斯点 | \(\xi_i\) | \(\eta_i\) | \(W_i\) | \(W_j\) |
|---|---|---|---|---|
| 1 | -0.77460 | -0.77460 | 0.55556 | 0.55556 |
| 2 | 0.0 | -0.77460 | 0.88889 | 0.55556 |
| 3 | 0.77460 | -0.77460 | 0.55556 | 0.55556 |
| 4 | 0.77460 | 0.0 | 0.55556 | 0.88889 |
| 5 | 0.77460 | 0.77460 | 0.55556 | 0.88889 |
| 6 | 0.0 | 0.77460 | 0.88889 | 0.55556 |
| 7 | -0.77460 | 0.77460 | 0.55556 | 0.55556 |
| 8 | -0.77460 | 0.0 | 0.55556 | 0.88889 |
| 9 | 0.0 | 0.0 | 0.88889 | 0.88889 |
使用双线性四边形得到的解涉及每个单元的\(4\times 4\)矩阵乘法；而二次四边形则处理\(8\times 8\)矩阵，或者每个单元的产品计算是线性单元的四倍。此外，双线性单元需要在四个高斯点进行评估，而二次单元则涉及九个高斯点。因此，在构建四节点双线性单元的单元刚度矩阵方面，操作成本是获得八节点二次单元的单元刚度矩阵的九分之一。这在选择要使用的单元类型时可能是一个重要的考虑因素，因为它意味着CPU时间几乎相差一个数量级。
4.6 稳态传导方程
我们再次关注确定二维域内稳态热传导。此时，我们将热传导矩阵推广以允许各向异性材料，假设坐标轴与导热张量的主轴重合，即热传导在\(x-y\)坐标系统中由 \[\begin{equation}\boldsymbol{K} = \begin{bmatrix} K_{xx} & 0 \\ 0 & K_{yy} \\ \end{bmatrix}\end{equation}\] 给出。然后基本方程可以写成 \[\begin{equation} - \frac{\partial}{\partial x} \left( K_{xx} \frac{\partial T}{\partial x} \right) - \frac{\partial}{\partial y} \left( K_{yy} \frac{\partial T}{\partial y} \right) = Q \end{equation}\] 现在热流条件采取的形式 \[\begin{equation} -\left(K_{xx} \frac{\partial T}{\partial x} n_x + K_{yy} \frac{\partial T}{\partial y} n_y\right) = q \end{equation} \] 其中\(n_x\)和\(n_y\)定义如方程（3.29）。按照与方程（3.26）相同的程序，方程27的加权残差形式是 \[\begin{equation} \int_\Omega W \left[ - \frac{\partial}{\partial x} \left( K_{xx} \frac{\partial T}{\partial x} \right) - \frac{\partial}{\partial y} \left( K_{yy} \frac{\partial T}{\partial y} \right) - Q \right]d\Omega=0 \end{equation}\] 应用格林定理后， \[\begin{equation} \int_\Omega\left[ K_{xx} \frac{\partial W}{\partial x} \frac{\partial T}{\partial x} + K_{yy} \frac{\partial W}{\partial y} \frac{\partial T}{\partial y}-WQ\right]dxdy + \int_{\Gamma_B} W \left( -K_{xx} \frac{\partial T}{\partial x} n_x - K_{yy} \frac{\partial T}{\partial y} n_y \right)dΓ=0 \end{equation}\]
我们现在使用双线性形状函数近似\(T\)，即 \[\begin{equation} T = \sum_{i=1}^4 N_i T_i \end{equation} \] 设\(W=N_i\)。方程26的Galerkin公式变为 \[\begin{equation} \int_\Omega \left( K_{xx} \frac{\partial N_i}{\partial x} \frac{\partial N_j}{\partial x} + K_{yy} \frac{\partial N_i}{\partial y} \frac{\partial N_j}{\partial y} \right) T_jdxdy = \int_\Omega N_i Qdxdy-\int_{\Gamma_B}N_iqd\Gamma \end{equation}\] 其中方程29已被用来替换边界积分中的施加热流。我们现在将方程32转换为自然\(\xi，\eta\)坐标系统。我们在该系统中定义温度为 \[\begin{equation} T(\xi, \eta) = \sum_{i=1}^4 N_i(\xi, \eta) T_i \end{equation}\] 其中形状函数\(N_i(\xi，\eta)\)在方程17中给出，定义在图5所示的通用单元上。第一导数值从 \[ \begin{equation} \begin{bmatrix} \frac{\partial N_i}{\partial x}\\ \frac{\partial N_i}{\partial y} \end{bmatrix} = J \begin{bmatrix} \frac{\partial N_i}{\partial \xi} \\ \frac{\partial N_i}{\partial \eta} \\ \end{bmatrix} \end{equation}\] 其中 \[ \begin{equation} J = \begin{bmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi} \\ \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta} \\ \end{bmatrix}= \begin{bmatrix} J_{11}&J_{21}\\ J_{12}&J_{22} \end{bmatrix}=\frac{1}{2}\begin{bmatrix} x_2-x_1&0\\ 0&y_4-y_1 \end{bmatrix} \end{equation}\] 使用方程13获得。雅可比矩阵的逆如下式所示 \[\begin{equation} \boldsymbol{J}^{-1} =\frac{1}{|\boldsymbol{J}|}\begin{bmatrix} J_{22} & -J_{12} \\ -J_{21} & J_{11} \\ \end{bmatrix} \end{equation}\] 其中 \[\begin{equation} |J| = J_{11}J_{22} - J_{12}J_{21} \end{equation}\] 因此，方程32变为 \[\begin{equation} \int^1_{-1} \int^1_{-1} \left[K_{xx}\left(\frac{\partial N_i}{\partial\xi}\frac{\partial\xi}{\partial x}+\frac{\partial N_i}{\partial\eta}\frac{\partial\eta}{\partial x}\right)\left(\frac{\partial N_j}{\partial\xi}\frac{\partial\xi}{\partial x}+\frac{\partial N_j}{\partial\eta}\frac{\partial\eta}{\partial x}\right)\\ +K_{yy}\left(\frac{\partial N_i}{\partial\xi}\frac{\partial\xi}{\partial y}+\frac{\partial N_i}{\partial\eta}\frac{\partial\eta}{\partial y}\right)\left(\frac{\partial N_j}{\partial\xi}\frac{\partial\xi}{\partial y}+\frac{\partial N_j}{\partial\eta}\frac{\partial\eta}{\partial y}\right) \right]|\boldsymbol{J}|d\xi d\eta\{T_j\} = -\int^1_{-1}qN_id\Gamma+ \int^1_{-1}\int^1_{-1}QN_id\xi d\eta \end{equation}\] 其中\(\frac{\partial N_i}{\partial x}\)和\(\frac{\partial N_i}{\partial y}\)用方程34表示为\(\frac{\partial N_i}{\partial \xi}\)和\(\frac{\partial N_i}{\partial \eta}\)。 ## 4.7 稳态传导与边界对流 如上一篇文章所讨论的，边界上的对流发生会改变边界条件，即载荷向量和刚度矩阵。稳态传导方程与方程27相同，但沿\(\Gamma_B\)的边界条件由下式给出： \[ \begin{equation} -\left(K_{xx}\frac{\partial T}{\partial x}n_x+K_{yy}\frac{\partial T}{\partial y}n_y\right)=h(T-T_\infty) \end{equation} \] 其中\(h(T − T_\infty)\)项考虑了由于对流引起的表面通量；否则，方程组与刚刚在4.6节讨论的稳态传导问题类似，热通量由方程28给出。应用Galerkin过程并通过部分积分，方程27变为 \[ \begin{equation} \left[\int\int_\Omega\left(K_{xx}\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial x}+K_{yy}\frac{\partial N_i}{\partial y}\frac{\partial N_j}{\partial y}\right)dxdy\right]T_j\\ =-\int_{\Gamma_B}h(T-T_\infty)N_id\Gamma+\int\int_\Omega QN_idxdy \end{equation} \] 其中方程39已被用来定义对流梯度边界条件项。评估方程40的过程与之前完全相同；在这种情况下，我们必须考虑包含对流热通量的表面积分。当单元边是直线时，这些积分的评估可以进行解析。当单元不再是矩形（或者有曲面，如我们稍后在二次单元中将看到的），必须进行数值积分。
4.8 二次四边形单元
如4.3.2节所讨论的，最常用的二次四边形单元由八个节点组成——四个角节点加上四个中边节点，如图7所示。虽然使用二次单元需要更多的计算工作，但精度的提高，或者对稳态问题的收敛解的更接近的近似，是显著的。双线性单元的误差减少率是\(O(h^2)\)，其中h是单元大小；对于二次单元，误差减少为\(O(h^3)\)。使用二次单元获得解决方案的过程与双线性单元完全相同。
4.9 时间依赖性扩散
如果一个人希望以瞬态、时间依赖的模式解决问题，与时间导数项相关的质量矩阵\(\boldsymbol{M}\)为 \[\begin{equation} \int_\Omega\frac{\partial T}{\partial t}N_idxdy=\left[\int^1_{-1}\int^1_{-1}N_iN_j|\boldsymbol{J}|d\xi d\eta\right]\left[\frac{\partial T_j}{\partial t}\right]=\boldsymbol{M}\left[\frac{\partial T_i}{\partial t}\right] \end{equation} \] ## 4.10 小结 我们增加了第二组单元，即四边形单元。2-D四边形单元，像上一篇讨论的三角形单元一样，是简单的单元，也可以适应任何不规则几何形状，尽管这一点尚不明显，但在下一章中将清楚地展示。四边形和三角形单元之间的两个主要区别是 1.双线性四边形单元需要一个额外的第四个节点，从而增加了一个自由度来定义单元2.梯度是坐标方向的线性函数，而三节点三角形（单纯形）单元中的梯度是常数。常数梯度需要在变量快速变化的区域使用许多小单元。 由于（矩形）四边形是多重单元，单元的边界必须与坐标系平行以保持单元之间的连续性。然而，没有理由为什么自然坐标系必须与笛卡尔（或圆柱）正交坐标系对齐——它也可以是曲线坐标系统。在大多数情况下，四边形单元可能比三角形单元（对于相同数量的节点）产生更准确的结果。关于各种单元的准确性和收敛性质的更广泛讨论可以在Carey和Oden (1983)，Oden和Carey (1983)，Oden和Reddy (1976)，Zienkiewicz和Taylor (1989)，以及Carey (1997)中找到。在第6章中，我们将使用一个坐标变换来准确表示非矩形区域和具有弯曲边界的区域。该变换只需要一套关系来考虑插值函数和以x和y表示的单元几何形状。
