2 一维有限元方法
2.1 背景
本部分以加权残差法的Galerkin形式为基础介绍有限元方法。有限元方法的第一步是将感兴趣的区域分割成包含节点的子区域，或单元，并在每个单元上定义要使用的“形状”类型。一个二阶单元有三个节点，它被优化地放置在单元的中点；一个三阶单元由四个节点组成，最好在单元的三分之一间隔处放置。单元和节点的集合称为“有限元网格”。在程序中使用的形状函数的类型直接与选择的网格类型相关。
2.2 形状函数
在第一篇文章中，我们使用分段多项式形状函数\(\phi_i(x)\)来处理一维（1-D）热传导问题。现在，我们将正式定义分段多项式插值的过程，并展示对特定的单元使用每个单元的局部函数更方便。
2.2.1 线性单元
让我们从定义一个网格开始，在\(0 \leq x \leq L\)的范围内，单元不一定是等大小的。然后我们在每个单元上定义线性形状函数，如图1所示。
这将是全局有限元网格和形状函数的表示，它允许我们描述我们问题的几何形状。如果网格由\(n\)个单元组成，它将有\(n+1\)个节点，其坐标为\(x_1, x_2, \cdots, x_\{n+1\}\)。单元域由下式给出： \[ e_i = \{x | x_{i-1} \leq x \leq x_{i+1}\} \tag{2.1}\]
与每个节点相关联的形状函数将被记为\(N_i(x)\)，并且满足\(N_i(x_i) = 1\)，\(N_j(x_i) = 0\)（对于所有\(j \neq i\)），由下式给出：
\[ N_1(x) = \begin{cases} \dfrac{x_2-x}{h_1} & \text{if } x \in [x_1, x_2] \\ 0 & \text{otherwise} \end{cases}\tag{2.2} \] \[ N_i(x) = \begin{cases} \dfrac{x-x_{i-1}}{h_{i-1}} & \text{if } x \in [x_{i-1}, x_i] \\ \dfrac{x_{i+1}-x}{h_i} & \text{if } x \in [x_i, x_{i+1}]\\ 0 & \text{otherwise} \end{cases}\tag{2.3} \] \[ N_{n+1}(x) = \begin{cases} \dfrac{x-x_n}{h_n} & \text{if } x \in [x_n, x_{n+1}]\\ 0 & \text{otherwise} \end{cases}\tag{2.4} \]
其中\(h_i = x_{i+1} - x_{i-1}\)，对于\(i = 1, 2, \ldots, n\)。 我们现在定义\(T(x)\)如下： \[T(x)=N_1(x)T_1+N_2(x)T_2+\cdots+N_{n+1}(x)T_{n+1}\\ =\sum^{n+1}_{i=1}N_i(x)T_i\tag{2.5}\] 式中\(T_i\)表示节点\(i\)，即位置\(x_i\)处的\(T\)值；\(T(x)\)为两个节点间的线性函数。我们考虑仅有一个单元的网格，也就是\(n=1\)。如果\(T(x)\)在两个节点之间为线性，其形式将为\(T(x)=\alpha_1+\alpha_2x\)。同时\(T(x_1)=T_1\)，\(T(x_2)=T_2\)，由此 \[T_1=\alpha_1+\alpha_2x_1\tag{2.6} \] \[T_2=\alpha_1+\alpha_2x_2\tag{2.7} \] 求解\(\alpha_1\)和\(\alpha_2\)， \[\alpha_1=\frac{T_1x_2-T_2x_1}{h_1}\tag{2.8} \] \[\alpha_2=\frac{T_2-T_1}{h_1}\tag{2.9} \] 因此 \[T(x)=\frac{T_1x_2-T_2x_1}{h_1}+\frac{T_2-T_1}{h_1}x\tag{2.10}\] 整理得到 \[T(x)=\Big [\frac{x_2-x}{h_1}\Big]T_1+\Big [\frac{x-x_1}{h_1}\Big]T_2\tag{2.11}\] 这与方程2.5的表达式相同，其中\(n = 1\)。在实践中，如果我们在局部坐标系统中需要描述一个长度为\(h^{(e)}\)的单元，我们将使用以下形式： \[T^{(e)}(x)=N^{(e)}_1(x)T^{(e)}_1+N^{(e)}_2(x)T^{(e)}_2\tag{2.12}\]
式中的\((e)\)表示的一个单元， \[N^{(e)}_1(x)=\Big [1-\frac{x}{h^{(e)}}\Big ]\tag{2.13}\] \[N^{(e)}_2(x)=\frac{x}{h^{(e)}\tag{2.14}}\] 如图2所示，有 \[\frac{dT^{(e)}}{dx}=\left[ \begin{matrix} \dfrac{dN^{(e)}_1}{dx}&\dfrac{dN^{(e)}_2}{dx} \end{matrix} \right] \left[ \begin{matrix}T^{(e)}_1\\T^{(e)}_2 \end{matrix}\right] =\left[ \begin{matrix} -\dfrac{1}{h^{(e)}}&\dfrac{1}{h^{(e)}} \end{matrix}\right] \left[ \begin{matrix}T^{(e)}_1\\T^{(e)}_2 \end{matrix}\right]\tag{2.15} \]
2.2.2 二次单元
从插值函数的角度来看，我们可以从图4中看到，如果我们在每个单元上使用抛物线弧而不是线性段，效果会更好。在每个单元上，函数\(T^{(e)}(x)\)将是二次的，因此形式为：
\[ T^{(e)}(x) = \alpha_1 + \alpha_2x + \alpha_3x^2\tag{2.16} \]
其中我们需要三个参数来确定；因此，仅仅要求在单元的两端插值是不够的。为了得到第三个关系，我们在单元的中间引入了另一个节点，并且我们还要获得该节点处的插值函数。图5为局部坐标系下的结果。
函数\(N^{(e)}_i\)由下式得到 \[T^{(e)}(0)=\alpha^{(e)}_1=T^{(e)}_1\tag{2.17a}\]
\[T^{(e)}(h/2)=\alpha_1+\frac{\alpha_2h^{(e)}}{2}+\frac{\alpha_3(h^{(e)})^2}{4}=T^{(e)}_2\tag{2.17b}\]
\[T^{(e)}(h)=\alpha_1+\alpha_2
h^{(e)}+\alpha_3(h^{(e)})^2=T^{(e)}_3\tag{2.17c}\]
单元的长度由下式给出： \[h^{(i)}=x_{2i+1}-x_{2i-1}\tag{2.20}\] 函数\(T^{(e)}(x)\)的导数分别由方程（2.20）和方程（2.21）得到: \[T^{(e)}(x)=N^{(e)}_1T^{(e)}_1+N^{(e)}_2(x)T^{(e)}_2+N^{(e)}_3(x)T^{(e)}_3= \left[ \begin{matrix} N^{(e)}_1&N^{(e)}_2&N^{(e)}_3 \end{matrix} \right] \left[ \begin{matrix} T^{(e)}_1\\ T^{(e)}_2\\ T^{(e)}_3 \end{matrix} \right]\tag{2.21} \] 以及 \[\frac{dT^{(e)}}{dx}=\left[ \begin{matrix} \dfrac{1}{h^{(e)}}\Bigg (4\dfrac{x}{h^{(e)}}-3\Bigg)& \dfrac{4}{h^{(e)}}\Bigg (1-2\dfrac{x}{h^{(e)}}\Bigg)& \dfrac{1}{h^{(e)}}\Bigg (4\dfrac{x}{h^{(e)}}-1\Bigg) \end{matrix}\right] \left[ \begin{matrix} T^{(e)}_1\\ T^{(e)}_2\\ T^{(e)}_3 \end{matrix} \right]\tag{2.22} \] 需要注意的是每个单元不再是常数。
2.2.3 三次单元
我们可以继续提高近似的阶数。在这种情况下，我们使用三次函数，即 \[ T^{(e)}(x)=\alpha_1+\alpha_2x+\alpha_3x^2+\alpha_4x^3\tag{2.23} \] 在这种情况下，每个单元需要四个节点。对于最佳的近似属性，这些节点最好位于\(x = 0, h^{(e)}/3, 2h^{(e)}/3, h^{(e)}\)的位置。按照与线性和二次单元相同的过程，我们可以找到形状函数。 \[N^{(e)}_1(x)=\Bigg(1-\frac{3x}{h^{(e)}}\Bigg)\Bigg(1-\frac{3x}{2h^{(e)}}\Bigg)\Bigg(1-\frac{x}{h^{(e)}}\Bigg)\tag{2.24a}\] \[N^{(e)}_2(x)=\frac{9x}{h^{(e)}}\Bigg(1-\frac{3x}{2h^{(e)}}\Bigg)\Bigg(1-\frac{x}{h^{(e)}}\Bigg)\tag{2.24b}\] \[N^{(e)}_3(x)=-\frac{9x}{2h^{(e)}}\Bigg(1-\frac{3x}{h^{(e)}}\Bigg)\Bigg(1-\frac{x}{h^{(e)}}\Bigg)\tag{2.24c}\] \[ N^{(e)}_4(x)=\frac{x}{h^{(e)}}\Bigg(1-\frac{3x}{h^{(e)}}\Bigg)\Bigg(1-\frac{3x}{2h^{(e)}}\Bigg)\tag{2.24c} \] 关于之前定义的单元，需要注意以下几点： 1. 尽管二次和三次单元的导数是关于独立变量\(x\)的函数，但它们在单元间的节点处不会连续。这里使用的插值类型被称为拉格朗日插值，它只保证函数在单元边界处的连续性。这些单元被称为\(C^0\)单元，其中零上标意味着只有零阶导数是连续的，即函数本身。 2. 显然，甚至更高阶的单元，例如四次、五次等，可以通过在单元中添加更多的插值节点来构建。实际上，我们可以构建在节点处也插值导数的单元。最简单的这类单元是立方埃尔米特单元，它在位于单元两端的节点处插值函数及其一阶导数。这些是\(C^1\)单元，因为现在一阶导数将在整个域中连续。甚至可以构建更复杂的单元。事实上，几乎没有限制可以达到的复杂度或预设单元行为的程度。然而，必须牢记，单元越复杂，计算成本就越高。实际上，在多维计算中，立方单元已经变得过于昂贵，很少用于流体和热传递应用。 3. 之前考虑的单元插值函数具有如下性质：\(N_i(x)e_j = \delta_{ij}\)，其中 \(\delta_{ij}\) 是克罗内克 delta 函数，即， \[ \delta_{ij} = \begin{cases} 1, & \text{if } i = j \\ 0, & \text{if } i \neq j \end{cases}\tag{2.25} \]
并且 \(x_j\) 是节点坐标。 4. 所选择的插值类型定义了依赖变量在单元内可以采取的“形状”，即线性的、二次等。因此，使用“形状函数”这个名称来表示定义单元的函数 \(N_i\)。还应注意，高阶函数总是精确地简化为低阶函数；即，二次单元精确地表示线性函数和常数函数，立方单元表示二次、线性和常数函数等。这一事实保证了如果使用高阶单元，可以获得更好的近似。
2.3 稳态传导方程
2.3.1 Galerkin表述
考虑在区间\(0 \leq x \leq L\)上找到温度分布\(T = T(x)\)的问题，它满足具有内部热源的稳态1-D方程：
\[ -K\frac{d^2T}{dx^2}=Q, 0<x<L\tag{2.26} \] 边界条件为 \[-K\frac{dT}{dx}=q,x=0\tag{2.27} \] \[T=T_L,x=L\tag{2.28}\] 方程（2.26）的加权残差表达式为 \[\int^L_0W\Bigg(-K\frac{d^2T}{dx^2}-Q\Bigg)dx=0\tag{2.29} \] 进一步写成 \[ \int^L_0W\Bigg(-K\frac{d^2T}{dx^2}-Q\Bigg)dx=\sum^2_{i=1}\int^{x_{i+1}}_{x_i}W\Bigg(-K\frac{d^2T}{dx^2}-Q\Bigg)dx =\int^{L/2}_0W\Bigg(-K\frac{d^2T}{dx^2}-Q\Bigg)dx+ \int^{L}_{L/2}W\Bigg(-K\frac{d^2T}{dx^2}-Q\Bigg)dx=0\tag{2.30} \] 全局函数\(T(x)\)由下式近似得到： \[ T(x)=\sum^3_{j=1}N_j(x)T_j=N_1(x)T_1+N_2(x)T_2+N_3(x)T_3 \tag{2.31}\] 式中形状函数\(N_i(x),i=1,2,3\)，由方程（2.2）至（2.4）给出，其中\(n=2\)，\(x_1=0,x_2=L/2,x_3=L\)。 令\(W_i(x)=N_i(x)\)，式（2.30）的伽辽金形式变为 \[ \int^{L/2}_0\left[K\frac{dN_i}{dx}\Bigg(\sum^3_{j=1}\frac{dN_j}{dx}T_j\Bigg)-N_iQ\right]dx+\left[N_i\Bigg(-K\frac{dT}{dx}\Bigg)\right]^{L/2}_0 +\int^L_{L/2}\left[K\frac{dN_i}{dx}\Bigg(\sum^3_{j=1}\frac{dN_j}{dx}T_j\Bigg)-N_iQ\right]dx+\left[N_i\Bigg(-K\frac{dT}{dx}\Bigg)\right]^L_{L/2} =0,i=1,2,3\tag{2.32} \] 方程（2.32）的前两项分别对应单元\(e_1\)和\(e_2\)。现在一次考虑其中的每一个，使用局部单元坐标。对于单元\(e_1\)，\(N_3(x) = 0\)；因此，我们可以将来自\(i = 1, 2\)的方程以矩阵形式写成 \[ \int^{L/2}_0\Bigg\{K\Bigg(\left[ \begin{matrix} \dfrac{dN^{(e_1)}_1}{dx}\\ \dfrac{dN^{(e_2)}_2}{dx}\\ \end{matrix} \right] \left[ \begin{matrix} \dfrac{dN^{(e_1)}_1}{dx}& \dfrac{dN^{(e_2)}_2}{dx} \end{matrix} \right]\Bigg) \left[ \begin{matrix} T^{(e_1)}_1\\ T^{(e_2)}_2 \end{matrix} \right]- \left[ \begin{matrix} N^{(e_1)}_1\\ N^{(e_2)}_2 \end{matrix} \right]Q \Bigg\}dx\\ +\Bigg\{ \left[ \begin{matrix} N^{(e_1)}_1\\ N^{(e_2)}_2 \end{matrix} \right]\Bigg(-K\frac{dT}{dx}\Bigg)\Bigg\}^{L/2}_0= \left[ \begin{matrix} 0\\0 \end{matrix} \right]\tag{2.33} \] 根据方程（2.14）和（2.15）以及\(h^{(e_1)}=L/2\)，方程（2.33）变为 \[ \int^{L/2}_0\Bigg\{K\left[ \begin{matrix} \dfrac{4}{L^2}&-\dfrac{4}{L^2}\\ -\dfrac{4}{L^2}&\dfrac{4}{L^2}\\ \end{matrix} \right] \left[ \begin{matrix} T^{(e_1)}_1\\ T^{(e_1)}_2 \end{matrix} \right]- \left[ \begin{matrix} 1-\dfrac{2x}{L}\\ \dfrac{2x}{L} \end{matrix} \right]Q \Bigg\}dx +\left[ \begin{matrix} -\Bigg( -K\dfrac{dT}{dx} \Bigg)_{x=0}\\ \Bigg( -K\dfrac{dT}{dx} \Bigg)_{x=L/2} \end{matrix} \right]=\left[ \begin{matrix} 0\\0 \end{matrix} \right] \] 应用方程（2.27）并积分，我们最终获得单元\(e_1\)的表达式 \[ \frac{2K}{L}\left[ \begin{matrix} 1&-1\\ -1&1 \end{matrix} \right] \left[ \begin{matrix} T_1\\T_2 \end{matrix} \right]=\frac{QL}{4}\left[ \begin{matrix} 1\\1 \end{matrix} \right]+ \left[ \begin{matrix} q\\ -\Bigg(-K\dfrac{dT}{dx}\Bigg)_{x=L/2} \end{matrix} \right] \] 在这里，我们使用了单元\(e_1\)中自由度的全局编号，即\(T_1\)和 \(T_2\)。以类似的方式，我们得到单元\(e_2\)的方程为: \[ \frac{2K}{L}\left[ \begin{matrix} 1&-1\\ -1&1 \end{matrix} \right] \left[ \begin{matrix} T_2\\T_3 \end{matrix} \right]=\frac{QL}{4}\left[ \begin{matrix} 1\\1 \end{matrix} \right]+ \left[ \begin{matrix} \Bigg(-K\dfrac{dT}{dx}\Bigg)_{x=L/2}\\ -\Bigg(-K\dfrac{dT}{dx}\Bigg)_{x=L} \end{matrix} \right] \] 现在必须将每个单元的积分加在一起，以得到完整的方程组。这个操作被称为“组装”方程，并产生最终的全局\(3\times 3\)矩阵。这是使用方程和自由度的全局编号来完成的；这个过程如图7所示，并且等同于将不同单元中由同一权重函数产生的方程贡献相加。在这种情况下，对应于 \(W_2=N_2\)的每个单元的贡献如图7所示相加。在实践中，这是通过使用单元矩阵中的全局自由度编号，并将它们的贡献添加到全局矩阵中的相应位置来实现的。例如，在单元\(2\)中，位置\((1,2)\)对应于全局系统中的\((2,3)\)，并将条目 \(-2k/L\) 添加到全局矩阵的位置\((2,3)\).在图7中可以清楚地看到，涉及内部节点通量的表达式相互抵消；实际上，这正好说明通量在内部必须是连续的。因此，在构建单元方程时，这些项通常被省略。然而，如果省略了这些项，单元方程中的等式就不成立。尽管如此，习惯上会将涉及通量的项，例如方程（3.35），设置为等于零。因此，全局方程组变为： \[ \frac{2K}{L} \left[ \begin{matrix} 1&-1&0\\ -1&2&-1\\ 0&-1&1 \end{matrix} \right]\left[ \begin{matrix} T_1\\ T_2\\ T_3 \end{matrix} \right]= \frac{QL}{4} \left[ \begin{matrix} 1\\2\\1 \end{matrix} \right] + \left[ \begin{matrix} q\\ 0\\ -\Bigg(-K\dfrac{dT}{dx}\Bigg)_{x=L} \end{matrix} \right]\tag{2.36} \] 在最后一个等式中，在\(x=L\)处出现了热通量。然而，如将其认为\(T_3\)处的公式，因为\(T_3\)已知且\(x=L\)处的热通量在前两个方程中没有出现，这个方程可以忽略并且重写为 \[ \frac{2K}{L}\left[ \begin{matrix} 1&-1\\ -1&2 \end{matrix} \right] \left[ \begin{matrix} T_1\\ T_2 \end{matrix} \right]= \frac{QL}{4} \left[ \begin{matrix} 1\\ 2 \end{matrix} \right]+ \left[ \begin{matrix} q\\ 0 \end{matrix} \right]+ \frac{2K}{L}T_L \left[ \begin{matrix} 0\\ 1 \end{matrix} \right]\tag{2.37} \] 上式可以求解出\(T_1\)和\(T_2\)。一旦求解完成后，第三个等式可用于计算右侧边界处的热传导，即 \[ \Bigg(-K\frac{dT}{dx}\Bigg)_{x=L}=\frac{2K}{L}(T_2-T_L)+\frac{QL}{4}\tag{2.38}\]
现在使用二次单元离散化区域，来近似方程（2.26）至（2.28）。我们现在使用相同数量的节点，但是是更高阶的近似。伽辽金描述与之前相似。由于局部和全局坐标系对于一个单元都是相同的，我们可以写成 \[\int^L_0\left[K\frac{dN_i}{dx}\Bigg(\sum^3_{j=1}\frac{dN_j}{dx}T_j\Bigg)-N_iQ\right]dx+\left[N_i\Bigg(-K\frac{dT}{dx}\Bigg)\right]^L_0=0\tag{2.39}\] 根据方程（3.18）、（3.21）和（3.23），得到矩阵形式： \[ \int^L_0\left\{K\left(\left[ \begin{matrix} \dfrac{1}{L}(\dfrac{4x}{L}-3)\\ \dfrac{4}{L}(1-\dfrac{2x}{L})\\ \dfrac{1}{L}(\dfrac{4x}{L}-1) \end{matrix} \right] \left[ \begin{matrix} \dfrac{1}{L}(\dfrac{4x}{L}-3)& \dfrac{4}{L}(1-\dfrac{2x}{L})& \dfrac{1}{L}(\dfrac{4x}{L}-1) \end{matrix} \right] \right) \left[ \begin{matrix} T_1\\T_2\\T3 \end{matrix} \right]\\ -\left[ \begin{matrix} 1-\dfrac{3x}{L}+2(\dfrac{x}{L})^2\\ \dfrac{4x}{L}(1-\dfrac{x}{L})\\ \dfrac{x}{L}(\dfrac{2x}{L}-1) \end{matrix} \right]Q \right\}dx= \left[ \begin{matrix} q\\0\\ -(-K\dfrac{dT}{dx})_{x=L} \end{matrix} \right]\tag{2.40} \] 积分后得到 \[ \frac{K}{6L}\left[ \begin{matrix} 14&-16&2\\ -16&32&-16\\ 2&-16&14 \end{matrix} \right] \left[ \begin{matrix} T_1\\T_2\\T_3 \end{matrix} \right]=\frac{QL}{6} \left[ \begin{matrix} 1\\4\\1 \end{matrix} \right]+ \left[ \begin{matrix} q\\0\\-(-K\dfrac{dT}{dx})_{x=L} \end{matrix} \right]\tag{2.41}\] 最后一个方程考虑\(T_3=T_L\)后消掉，最后的系统变为 \[ \frac{K}{6L}\left[ \begin{matrix} 14&-16\\ -16&32\\ \end{matrix} \right] \left[ \begin{matrix} T_1\\T_2 \end{matrix} \right]=\frac{QL}{6} \left[ \begin{matrix} 1\\4 \end{matrix} \right]+ \left[ \begin{matrix} q\\0 \end{matrix} \right]+\frac{K}{6L}T_L\left[ \begin{matrix} -2\\16 \end{matrix}\right]\tag{2.42}\] 右侧边界的热通量由下式给出 \[ \Bigg(-K\frac{dT}{dx}\Bigg)_{x=L}=\frac{K}{6L}(16T_2-2T_1-14T_L)+\frac{QL}{6}\tag{2.43}\] ### 2.3.2 可变传导和边界对流 我们现在将有限元算法扩展到解决薄杆上的稳态温度分布问题，该杆的左端受到对流热负荷，\(x = 0\)，右端保持固定温度\(T_L\)，\(x = L\)。我们还假设没有内部热源（即\(Q = 0\)），但杆的热导率随\(x\)变化（由于材料组成的变化或杆的横截面变化），即\(K = K(x)\)。描述这个问题的微分方程是：
\[ \frac{d}{dx}\left(K(x) \frac{dT}{dx}\right) = 0\tag{2.44} \] 其中 \[-K\frac{dT}{dx}+h(T-T_\infty)=0,x=0\\ T=T_L,x=L\tag{2.45}\] 式中，\(T_\infty\)为外部参考温度，\(h\)为热传递参数 加权残差形式为 \[ \int^L_0K(x)\frac{dW}{dx}\frac{dT}{dx}dx-\left[W\Bigg(-K(x)\frac{dT}{dx}\Bigg)\right]_{x=L}=0\tag{2.47} \] 从现在开始，我们将省略对应于规定温度的边界上的通量项，在这种情况下是\(x = L\)，这是有限元建模中的惯例。如果我们还用方程（2.45）替换左端边界的通量，我们有 \[ \int^L_0K(x)\frac{dW}{dx}\frac{dT}{dx}dx+Wh(T-T_\infty)|_{x=0}=0\tag{2.48}\] 请注意，从方程（2.44）到（2.48）的过程中并没有丢失任何信息；然而，我们不是寻求方程（2.48）的解析解，而是寻求利用我们定义的形状函数和网格离散化的可计算有限元近似。 与之前一样，我们在\(0\le x\le L\)上定义一个网格，并通过以下方式近似\(T(x)\)： \[ T(x)=\sum^{n+1}_{i=1}N_i(x)T_i\tag{2.49}\] 式中\(n\)为网格中的单元个数。如果\(N_i,i=1,\cdots,n+1\)为线性形状函数，则方程（2.48）的伽辽金形式为 \[ \int^L_0K(x)\frac{dN_i}{dx}\Bigg(\sum^{n+1}_{j=1}\frac{dN_j}{dx}T_j\Bigg)dx+N_ih(T-T_\infty)|_{x=0}=0\tag{2.50}\] 不代入\(K(x)\)的特定形式，我们将使用与\(T(x)\)选择的相同形状函数（在这种情况下是线性的）来插值\(K(x)\)的节点值，即在每个单元上\(K^{(e)}(x) = N^{(e)}_1(x)K^{(e)}_1 +N^{(e)}_2(x)K^{(e)}_2\)。第一个单元（包括对流边界条件）的方程为： \[ \int^{x_2}_{x_1}\left[ \begin{matrix} N^{(e_1)}_1&N^{(e_2)}_2 \end{matrix} \right] \left[ \begin{matrix} K^{(e_1)}_1\\ K^{(e_2)}_2 \end{matrix} \right] \Bigg( \left[ \begin{matrix} \dfrac{dN^{(e_1)}_1}{dx}\\ \dfrac{dN^{(e_2)}_2}{dx} \end{matrix} \right] \left[ \begin{matrix} \dfrac{dN^{(e_1)}_1}{dx}& \dfrac{dN^{(e_2)}_2}{dx} \end{matrix} \right] \Bigg)dx \left[ \begin{matrix} T^{(e_1)}_1\\ T^{(e_2)}_2 \end{matrix} \right]\\ +\left[\begin{matrix} h&0\\ 0&0 \end{matrix} \right]\left[ \begin{matrix} T^{(e_1)}_1\\ T^{(e_2)}_2 \end{matrix} \right]=hT_\infty\left[ \begin{matrix} 1\\0 \end{matrix} \right] \] 考虑到在守恒边界项有\(N^{(e_1)}_1(0)=1\)和\(N^{(e_1)}_2(0)=1\)，即 \[ N_ih(T-T_\infty)|_{x=0}=\left[ \begin{matrix} N_1\\N_2 \end{matrix} \right]h\left[ \begin{matrix} N_1&N_2 \end{matrix} \right] \left[ \begin{matrix} T_1-T_\infty\\ T_2-T_\infty \end{matrix} \right]\Bigg|_{x=0}\\ =\Bigg( \left[ \begin{matrix} N_1\\N_2 \end{matrix} \right] \left[ \begin{matrix} N_1&N_2 \end{matrix} \right] \left[ \begin{matrix} T_1\\T_2 \end{matrix} \right]- \left[ \begin{matrix} N_1\\N_2 \end{matrix} \right]hT_\infty \Bigg)\Bigg|_{x=0}\\ =h\left[ \begin{matrix} 1&0\\ 0&0 \end{matrix} \right] \left[ \begin{matrix} T_1\\T_2 \end{matrix} \right]-hT_\infty\left[ \begin{matrix} 1\\0 \end{matrix} \right] \] 积分后得到 \[ \frac{1}{2h^{(e_1)}}(K^{(e_1)}_1+K^{(e_1)}_2)\left[ \begin{matrix} 1&-1\\ -1&1 \end{matrix} \right] \left[ \begin{matrix} T^{(e_1)}_1\\ T^{(e_1)}_2 \end{matrix} \right]+ \left[ \begin{matrix} h&0\\ 0&0 \end{matrix} \right] \left[ \begin{matrix} T^{(e_1)}_1\\ T^{(e_1)}_2 \end{matrix} \right]=hT_\infty\left[ \begin{matrix} 1\\0 \end{matrix} \right]\tag{2.51} \] 对所有单元，\(e_i\neq e_1\)，方程为 \[ \frac{1}{2h^{(e_1)}}(K^{(e_1)}_1+K^{(e_1)}_2)=\left[ \begin{matrix} 1&-1\\ -1&1 \end{matrix} \right] \left[ \begin{matrix} T^{(e_1)}_1\\T^{(e_1)}_2 \end{matrix} \right]=\left[ \begin{matrix} 0\\0 \end{matrix} \right]\tag{2.52} \] 如果区间\(0\le x\le L\)被离散为2个相等长度的单元，组装后的方程为 \[ \frac{1}{L}\left[ \begin{matrix} K_1+K_2+Lh&-(K_1+K_2)&0\\ -(K_1+K_2)&K_1+2K_2+K_3&-(K_2+K_3)\\ 0&-(K_2+K_3)&K_2+K_3 \end{matrix} \right] \left[ \begin{matrix} T_1\\T_2\\T_3 \end{matrix} \right]=\left[ \begin{matrix} hT_\infty\\ 0\\0 \end{matrix} \right]\tag{2.53}\] 如果在\(x = L\)处给出了形式如方程（2.45）的对流边界条件，那么最后一个单元的方程可以从加权残差形式（2.47）得出。然后，忽略第二个单元中的第一个边界项，并在第二个项中替换对流条件，得到单元方程为 \[ \frac{1}{2h^{(e_n)}}(K^{(e_n)}_1+K^{(e_n)}_2)\left[ \begin{matrix} 1&-1\\ -1&1 \end{matrix} \right] \left[ \begin{matrix} T^{(e_n)}_1\\T^{(e_n)}_2 \end{matrix} \right]+ \left[ \begin{matrix} 0&0\\ 0&-h \end{matrix} \right] \left[ \begin{matrix} T^{(e_n)}_1\\T^{(e_n)}_2 \end{matrix} \right] =\left[ \begin{matrix} 0\\0 \end{matrix} \right]\tag{2.54} \] 对流边界条件可以应用于网格的任一或两个边界；方程（2.51）、（2.52）和（2.54）对于任何网格（均匀或非均匀）、变导热性和任一或两个边界处的对流都是有效的。在定义网格的节点坐标和各种问题数据 \(h\)、\(T_\infty\) 和 \(K(x)\) 后，可以获得特定解。 ## 2.4 轴对称热传导 许多涉及流体在管道中流动的稳态热传导问题都涉及到轴对称形式。为了开发相应的有限元算法，方程（2.44）至（2.46）的轴对称形式是：
\[ \frac{1}{r}\frac{d}{dr}\left(r\frac{dT}{dr}\right) = Q,r_1<r<r_2\tag{2.58} \]
加上边界条件：
\[ -\frac{dT}{dr} + h(T - T_{\infty}) = 0 \quad \text{at } r = r_1\tag{2.59} \]
和
\[ T = T_L \quad \text{at } r = r_2\tag{2.60} \]
其中 \(r_1\) 和 \(r_2\) 分别是内半径和外半径，\(T_{\infty}\) 是外部参考温度，\(h\) 是对流换热系数。 方程（2.58）的加权残差形式为 \[ \int^{2\pi}_0\int^{r_2}_{r_1}W\left[ -\frac{d}{dr}\Bigg(rK\frac{dT}{dr}\Bigg)\right]drd\theta=0\tag{2.61}\] 分别对\(r\)和\(\theta\)积分，我们得到 \[ 2\pi\int^{r_2}_{r_1}rK\frac{dW}{dr}\frac{dT}{dr}dr+(2\pi rWh(T-T_\infty))_{r=r_1}=0\tag{2.62}\] 比较方程（2.62）与方程（2.47），轴对称情况下的显著区别在于\(K(x)\)被替换为\(rK\)。物理上，平滑变化的热导率并不是特别重要，而突然变化（恒定）的热导率则很重要，例如，材料的变化或包裹在管道上的热绝缘可能会产生显著的变化。因此，可以通过将\(K\)移到积分之外来简化方程(2.62)，假设积分将根据\(K\)的阶跃变化分段计算。单元\(e_1=\{r|r^{(e_1)}_1\le r\le r^{(e_1)}_2\}\)轴对称的伽辽金积分可以写成 \[ \Bigg\{ \frac{K^{(e_i)}}{2h^{(e_i)}} (r^{(e_1)}_1+r^{(e_1)}_2) \left[ \begin{matrix} 1&-1\\ -1&1 \end{matrix} \right] +r_1h\left[ \begin{matrix} 1&0\\ 0&0 \end{matrix} \right]\delta_{e_1} -r_2h\left[ \begin{matrix} 0&0\\ 0&1 \end{matrix} \right]\delta_{e_n} \left[ \begin{matrix} T^{(e_i)}_1\\ T^{(e_i)}_2 \end{matrix} \right] \Bigg\}\\ =\Bigg\{ r_1hT_\infty\left[ \begin{matrix} 1\\0 \end{matrix} \right]\delta_{e_1}- r_2hT_\infty\left[ \begin{matrix} 0\\1 \end{matrix} \right]\delta_{e_n} \Bigg\}\tag{2.63} \] 式中，\(\delta_{e_i}\)为Kronecker delata函数，即在单元\(e_i\)处\(\delta_{e_i}=1\)，而对所有\(j\neq i\)的单元\(e_j\)处为0； \(n\)表示单元编号，所以\(e_n\)为右侧边界的单元。 ## 2.5 自然坐标系统 我们现在引入一个自然坐标系统来指定单元上的局部操作。我们将单元定义在区间 \(-1 \leq \xi \leq 1\)，并与变换形式相结合： \[ x=\frac{1}{2}(1-\xi)x_i+\frac{1}{2}(1+\xi)x_{i+1}\tag{2.64} \] 这个形式将区间\(-1\le\xi\le 1\)映射到了单元的区间\(x_i\le x\le x_{i+1}\)上。线性单元形状函数如图8所示，并由下式给出： \[ N_1(\xi)=\frac{1}{2}(1-\xi)\\ N_2(\xi)=\frac{1}{2}(1+\xi) \tag{2.65} \]
使用这种类型的坐标系统，如果注意到，可以通过变换（2.64）来分析形状函数及其导数的乘积： \[ \frac{d}{dx}=\frac{d}{d\xi}\frac{d\xi}{dx}=\frac{2}{h^{(e_i)}}\frac{d}{d\xi}\tag{2.66}\] 和 \[ dx=\frac{h^{(e_i)}}{2}d\xi\tag{2.67} \] 注意到方程（2.64）可以写成 \[ x=N_1(\xi)x_1+N_2(\xi)x_2 \] 因此，变换可以严格地用形状函数和节点坐标来表示。形状函数的积分和导数具有以下形式： \[ \int^{x_{i+1}}_{x_i}N_j(x)N_k(x)dx=\frac{h{(e_1)}}{2}\int^1_{-1}N_j(\xi)N_k(\xi)d\xi\tag{2.68} \] \[ \int^{x_{i+1}}_{x_i}\frac{dN_j}{dx}\frac{dN_k}{dx}dx=\frac{2}{h^{(e_i)}}\int^1_{-1}\frac{dN_j}{d\xi}\frac{dN_k}{d\xi}d\xi\tag{2.69} \] \[ \int^{x_{i+1}}_{x_i}\frac{dN_j}{dx}N_k(x)dx=\int^1_{-1}N_k(\xi)d\xi\tag{2.70} \] 此外，方程(2.68)中的函数乘积可以很容易地用指数形式进行解析评估。这个程序在结构力学中是众所周知的，并且可以很容易地用于任何一维单元网格。积分关系表达为 \[ \int^1_{-1}(N_1(\xi))^a(N_2(\xi))^bd\xi=2\frac{a!b!}{(1+a+b)!}\tag{2.71} \] 式中\(a\)和\(b\)为无符号整数。我们可以对更高阶的单元使用自然坐标系，特别是二次单元。自然坐标系下的形状函数为 \[ N_1(\xi)=\frac{1}{2}\xi(\xi-1)\tag{2.74} \] \[ N_2(\xi)=1-\xi^2\tag{2.75} \] \[ N_3(\xi)=\frac{1}{2}\xi(\xi+1)\tag{2.76} \] 在区间\(x_{2i-1}\le x\le x_{2i+1}\)的单元\(e_i\)映射到区间\(-1\le\xi\le 1\)的转化式为 \[ x=N_1(\xi)x_{2i-1}+N_2(\xi)x_{2i}+N_3(\xi)x_{2i+1}\tag{2.77} \] 注意，当 \(x_{2i} =(x_{2i-1} + x_{2i+1})/2\) 时，这个表达式简化为方程(2.64)。这种表示法比2.2中的表示法更具一般性，因为它允许内部节点位于单元中点以外的位置。 从\(N_i(x)\) 对\(x\)的导数可以通过以下方式获得： \[ \frac{dN_i}{d\xi}=\frac{dN_i}{dx}\frac{dx}{d\xi} \] \[ \frac{dN_i}{dx}=\frac{1}{dx/d\xi}\frac{dN_i}{d\xi}\tag{2.78} \] 于是有 \[ \frac{dx}{d\xi}=\frac{dN_1}{d\xi}x_{2i-1}+\frac{dN_2}{d\xi}x_{2i}+\frac{dN_3}{d\xi}x_{2i+1}\tag{2.79} \] 在（2.79）中使用了单元坐标的全局值\(x^{(e_i)}_1\)，\(x^{(e_i)}_2\)，\(x^{(e_i)}_3\)。\(dx/d\xi\)为坐标变换的雅可比行列式，通常表示为\(\boldsymbol{J}\)。
2.6 随时间变化的扩散方程
我们现在将有限元算法扩展到非稳态热扩散方程。假设 \(Q = 0\)，即没有源或汇。热传导的控制方程通常写为：
\[ \frac{\partial T}{\partial t}-\alpha\frac{\partial^2T}{\partial x^2}=0\tag{2.90} \]
其中\(\partial T/\partial
t\)表示温度随时间的变化率，\(\alpha =
k/(\rho c_p)\) 是热扩散率，\(\rho\) 是密度，\(c_p\) 是材料的比热容。
现在\(T = T(x,t)\)是关于空间（\(x\)）和时间（\(t\)）的函数。因此，除了边界条件外，我们还需要指定一个初始条件，也就是说，我们必须满足边界条件以及方程（2.90）：
\[
-K\frac{\partial T}{\partial x}+h(T-T_\infty)\Bigg|_{0,t}=0\tag{2.91}
\] \[
T(L,t)=T_L\tag{2.92}
\] 初始条件为 \[
T(x,0)=T_0(X)\tag{2.93}
\] 方程（2.91）和（2.92）也假设了广义记号 \(T_L(t) = T(t)\) 和 \(T_{\infty}(t) = T(t)\)；\(T_0(t)\) 是杆在时间\(t = 0\)时的温度分布。 在恒温假设下\(T_L\) 和\(T_{\infty}\)，方程（2.90）及其相关的边界条件和初始条件的解析解可以轻松获得。
2.6.1 空间离散
方程（2.90）的权重残差形式为 \[ \int W(x)\left[\frac{\partial T}{\partial t}-\alpha\frac{\partial^2T}{\partial x^2}\right]dx=0\tag{2.94} \] 需要注意的是，权重函数仅是关于\(x\)的函数，即严格与空间离散化相关。此外，我们将假设温度可以写成分离变量的形式，并在空间上使用与之前相同的形状函数进行近似，即： \[ T(x,t)=\sum^{n+1}_{i=1}N_i(x)T_i(t)\tag{2.95} \] 式中\(n\)为网格节点数。 时间依赖性不影响形状函数，并保持在依赖变量中。温度的对空间的导数由以下公式给出： \[\frac{\partial T}{\partial x}=\sum^{n+1}_{i=1}\frac{\partial N_i}{\partial x}T_i\equiv\left[\frac{\partial N_i}{\partial x}\right][T_i]\tag{2.96}\] 这与之前相同。我们将使用偏导数符号来表示形状函数和离散变量的导数，即使这些函数只依赖于一个独立变量，且这些导数是总导数。我们还有： \[ \frac{\partial T}{\partial t}=\sum^{n+1}_{i=1}N_i\frac{\partial T}{\partial t}\equiv\sum^{n+1}_{i=1}N_i\dot{T}_i\equiv[N_i][\dot{T}_i]\tag{2.97} \] 其中\(\dot{T}_i\)表示对时间的导数，\([N]\)和\([T]\)分别为行矩阵和列矩阵。 方程（2.94）的伽辽金形式为 \[ \int^L_0\Bigg\{N_i\Bigg(\sum^{n+1}_{j=1}N_j\dot{T}_j\Bigg)+\alpha\frac{\partial N_i}{\partial x}\Bigg(\sum^{n+1}_{j=1}\frac{\partial N_j}{\partial x}T_j\Bigg)\Bigg\}dx+\left[N_i\Bigg(-\alpha\frac{\partial T}{\partial x}\Bigg)\right]^{x=L}_{x=0}=0\tag{2.98}\] 在这个阶段，我们将引入标记符号来替换表达式中的求和符号: \[ \sum^{n+1}_{i=1}a_ib_i\equiv a_ib_i\tag{2.99}\] 于是方程（2.98）使用上述标记写为更简便的形式： \[ \left[\int^L_0N_iN_jdx\right]\dot{T}_j+\left[\alpha\int^L_0\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial x}dx\right]T_j+\left[N_i\Bigg(-\alpha\frac{\partial T}{\partial x}\Bigg)\right]^{x=L}_{x=0}=0\tag{2.100} \] 其中第一个积分中的时间项重写为 \[ \int^L_0N_i\frac{\partial T}{\partial t}dx=\Bigg(\int^L_0N_iN_jdx\Bigg)\dot{T}_j \] 方程中\(\dot{T}_j\)一般称作质量矩阵，这是由于\(N_iN_j\)表示包含节点\(i\)并与所有连接节点\(j\)相关联的单元的面积。方程（2.100）给出的形式被称为半离散伽辽金形式，因为只对空间变量进行了离散化，关于时间导数项还没有任何说明。 ### 2.6.2 时间离散化 在有限元方法中，有多种处理方程（2.100）的时间积分的方式。在这里，我们将介绍所谓的\(\theta\)-方法，它导致最常用的时间积分算法。在\(\theta\)-方法中，时间导数被简单差分所代替： \[ \frac{\partial T}{\partial t}=\frac{T^{n+1}-T^n}{\Delta t}\tag{2.102} \] 式中\(T^n=T(x,t_n)\)表示时间\(t=t_n\)时的变量值，\(\Delta t\)为时间增量，\(t_{n+1}=t_n+\Delta t\)。一般地，我们假设\(T(x,t_n)\)已知，并用作初始条件，用于推进解到时间水平\(t_{n+1}\)。我们现在引入一个弛豫参数\(q\)，并将解\(T\)写成以下形式： \[ T=\theta T^{n+1}+(1-\theta)T^{n},t_n\le t\le t_{n+1}\tag{2.102} \] 参数\(\theta\)通常指定在\(0\le\theta\le 1\)的范围内，用于控制算法的准确性和稳定性。最常用的\(\theta\)值是\(0\)、\(1/2\)和\(1\)。众所周知，当\(\theta<1/2\)时，仅达到条件稳定性。 \(\theta=1\)为向后隐式方法， \(\theta=1/2\)给出一个第二阶中心隐式方法， \(\theta=0\)给出显式欧拉前向方案。 将方程（2.101）和（2.102）代入方程（2.100），我们得到： \[ \Bigg(\int^L_0N_iN_jdx\Bigg)\Bigg(\frac{T^{n+1}-T^n_j}{\Delta t}\Bigg)+\Bigg(\alpha\int^L_0\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial x}dx\Bigg)(\theta T^{n+1}_j+(1-\theta)T^n_j)\\ +\left[N_i\Bigg(-\alpha\Bigg\{\theta\frac{\partial T^{n+1}}{\partial x}+(1-\theta)\frac{\partial T^n}{\partial x}\Bigg\}\Bigg)\right]^{x=L}_{x=0}=0 \] 可重写为 \[ \left[\int^L_0N_iN_jdx+\alpha\Delta t\theta\int^L_0\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial x}dx\right]T^{n+1}_j+\theta\Delta t\left[N_i\Bigg(-\alpha\frac{\partial T^n}{\partial x}\Bigg)\right]^{x=L}_{x=0}\\ =\left[\int^L_0N_iN_jdx+\alpha\Delta t(\theta-1)\int^L_0\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial x}dx\right]T^n_j+(\theta-1)\Delta t\left[N_i\Bigg(-\alpha\frac{\partial T^n}{\partial x}\Bigg)\right]^{x=L}_{x=0}\\ \tag{2.103} \] 在这个表达式中，可以看到当\(\theta=0\)时。等式左侧只有质量矩阵项。这个项导致相邻单元之间的“连接度”比传统有限差分方法中的更高。
2.7 矩阵公式化
有限元方法基于在特定节点位置对因变量的数值近似；产生一组必须直接或迭代求解的同时线性代数方程。对于之前讨论的示例问题，未知数的数量，因此方程的数量很少，可以通过手工推导的方式轻松解决。然而，对于大多数问题，节点的数量，由此产生的未知数的数量，将会大大增加，需要计算机来执行求解运算。
我们之前使用的系数矩阵是在每个单元上分析得到，并组装成全局数组。换句话说，从每个单元获得的局部系数矩阵被“塞进”一个大矩阵中，该矩阵包含所有局部单元的贡献。这个过程可以通过计算机程序中的“do循环”轻松执行。一旦我们为一个单元制定了有限元算法，该方法的通用性允许我们对所有单元使用相同的程序。因此，我们可以构建基于局部单元矩阵评估的全局矩阵，然后以我们希望的方式求解矩阵方程。
定义矩阵符号运算会很方便。时间导数项的局部质量矩阵为
$$M^{(e_k)}=[m^{e_k}_{ij}]=_{0}{h{(e_k)}} [ \begin{matrix}
N_i{(e_k)}&N_j{(e_k)}] \end{matrix}dx
$$ 其中 \(h\) 是单元的长度，\(N_i^{(e_k)}\) 和 \(N_j^{(e_k)}\) 是单元的形状函数。
全局矩阵 \(\boldsymbol{M}\) 表示为
\[ \boldsymbol{M}=[m_{ij}] = \sum_{k=1}^{n}
M^{(e_k)}=\sum^n_{k=1}\int_{0}^{h^{(e_k)}}
[
\begin{matrix}
N_i^{(e_k)}&N_j^{(e_k)}]
\end{matrix}dx\tag{2.109} \] 其中\(n\)是单元的数量。求和意味着每个单元矩阵是一个\((n+1)\times (n+1)\)
矩阵，通过在不包含单元中节点的所有位置填充零来从单元矩阵扩展而来。
例如，考虑使用三个等长的线性单元离散区间\(0 \leq x \leq
1\)。使用方程（2.108）评估局部质量矩阵得到 \[\boldsymbol{M}^{(e_1)}=\boldsymbol{M}^{(e_2)}=\boldsymbol{M}^{(e_3)}=\frac{1}{18}
\left[\begin{matrix}
2&1\\
1&2
\end{matrix}\right]\]
组装后的矩阵 \(\boldsymbol{M}\)为 \[ M = \begin{bmatrix} 1/9 & 1/18 & 0 & 0\\ 1/18 & 1/9 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ \end{bmatrix}+\begin{bmatrix} 0 & 0 & 0 & 0\\ 0 & 1/9 & 1/18 & 0 \\ 0 & 1/18 & 1/9 & 0 \\ 0 & 0 & 0 & 0 \\ \end{bmatrix}\\ +\begin{bmatrix} 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 \\ 0 & 0 & 1/9 & 1/18 \\ 0 & 0 & 1/18 & 1/9 \\ \end{bmatrix} =\begin{bmatrix} 1/9 & 1/18 & 0 & 0\\ 1/18 & 2/9 & 1/18 & 0 \\ 0 & 1/18 & 2/9 & 1/18 \\ 0 & 0 & 1/18 & 1/9 \\ \end{bmatrix}\]
由于矩阵\(\boldsymbol{M}\)
中许多条目是零，因此它被称为稀疏矩阵。实际上，除非节点编号方式非常人为，否则所有从有限元离散化中产生的全局矩阵都是稀疏的。在一维线性和二次单元的特殊情况中，如果节点按顺序编号，矩阵将分别为三对角线和五对角线矩阵。
类似地，扩散项表示为 \[\boldsymbol{K}=\sum^n_{k=1}\boldsymbol{K}^{(e_k)}=
\sum_{k=1}^{n} \int_{0}^{h^{(e_k)}}K^{(e_k)}
\left[\frac{dN_i^{(e_k)}}{dx}
\frac{dN_j^{(e_k)}}{dx}\right]dx\tag{2.110} \] 通常将矩阵\(\boldsymbol{K}\) 称为“刚度”或“传导”矩阵。
整体上包含包括如源项的已知函数的贡献被定义为列向量，长度等于节点数，我们将其表示为
\[
\boldsymbol{F}=\sum^n_{k=1}F^{(e_k)}=\sum_{k=1}^{n}\int_{0}^{h^{(e_k)}}[N_i^{(e_k)}]^TQ^{(e_k)}dx\tag{2.111}
\] 其中\(Q^{(e_k)}\)表示\(Q\)在单元\(e_k\)上的限制， \([N_i^{(e_k)}]^T\)是一个列矩阵
如果依赖变量出现，则从热通量边界条件的贡献被添加到刚度矩阵
(2.110)，或者如果只涉及已知数据，则添加到向量\(\boldsymbol{F}\)。通常将向量\(\boldsymbol{F}\)称为载荷向量。然后，时间依赖的传导方程可以表示为
\[\boldsymbol{M\dot{T}}+\boldsymbol{KT}=\boldsymbol{F}\tag{2.112}\]
其中 \(\boldsymbol{T}\)是节点未知数的向量，\(\boldsymbol{\dot{T}}\)是关于时间导数的向量。
使用方程(2.101)和(2.102)替换方程(2.112)中的时间导数\(\boldsymbol{\dot{T}}\)，得到完全离散化的线性代数方程组为
\[(\boldsymbol{M}+\theta\Delta
t\boldsymbol{K})\boldsymbol{T}^{n+1}=(\boldsymbol{M}+(\theta-1)\Delta
t\boldsymbol{K})\boldsymbol{T}^n+\Delta
t(\theta\boldsymbol{F}^{n+1}+(1-\theta)\boldsymbol{F}^n)\tag{2.113}\]
方程（2.113）为使用\(\theta\)-方法求解含时间参量的传导方程算法。这种表示的优点是，现在我们可以使用矩阵理论以符号方式处理方程（2.113）。例如，如果
\((\boldsymbol{M}+\theta\Delta
t\boldsymbol{K})\) 是一个可逆的非奇异矩阵，解\(\boldsymbol{T}^{n+1}\)由下式给出： \[
\boldsymbol{T}^{n+1}=(\boldsymbol{M}+\theta\Delta
t\boldsymbol{K})^{-1}[(\boldsymbol{M}+(\theta-1)\Delta
t\boldsymbol{K})+\Delta
t(\theta\boldsymbol{F}^{n+1}+(1-\theta)\boldsymbol{F}^n)]\tag{3.115}\]
在我们继续本书的其余部分时，我们将使用矩阵符号来简化我们的代数表达式。然而，我们也将说明构成计算机程序基础的矩阵的局部单元评估。
2.8 求解方法
在所有前面的部分中，将有限元方法应用于控制方程导致一组线性方程，包括与原始微分方程的各种项相关的矩阵，如方程（2.113）。我们可以进一步将这个方程表示为
\[A\phi=\boldsymbol{B}\tag{2.115}\]
其中 \[A =\boldsymbol{M}+\theta\Delta
t\boldsymbol{K}\tag{2.116}\] \[\phi=\boldsymbol{T}^{n+1}\tag{2.117}\]
\[\boldsymbol{B}=(\boldsymbol{M}+(\theta-1)\Delta
t\boldsymbol{K})\boldsymbol{T}^n+\Delta
t(\theta\boldsymbol{F}^{n+1}+(1-\theta)\boldsymbol{F}^n)\tag{2.118}\]
所有出现在\(\boldsymbol{A}\)和\(\boldsymbol{B}\)中的项都是已知的，所以我们可以很容易地求解方程（2.115）中的未知数\(\phi\)。然而，如前所述，求解大型方程组可能会变得非常耗时和昂贵，即使是使用大型计算机，更不用说使用个人计算机了。这种情况引入了解决大型线性方程组的特殊方法的需要，如方程（2.115）。处理这些方法的数学学科称为数值线性代数，具有一些背景的读者将通过高斯消元法、雅可比和高斯-塞德尔迭代、LU分解、逐步超松弛、共轭梯度等方法的名字来识别这些方法。数值线性代数中矩阵求解技术的发展成为了一个领域；我们的目的只是让读者熟悉基本思想，以便理解它们如何实现到计算机程序中。方程（2.115）的解是通过在方程两端乘以
\(\boldsymbol{A}^{-1}\)（注意\(\boldsymbol{A}^{-1}\boldsymbol{A}\equiv\boldsymbol{I}\)，\(\boldsymbol{I}\)为单位矩阵）。因此， \[\phi=\boldsymbol{A}^{-1}\boldsymbol{B}\tag{2.119}\]
如果\(\boldsymbol{A}\)很大，这很难执行。幸运的是，存在有效的数值方法，允许我们在不需要找到\(\boldsymbol{A}^{-1}\)的情况下得到\(\phi\)。矩阵代数，包括找到矩阵的逆，在附录A中有更详细的讨论。解决线性代数方程组系统基本上有两种方法。解决方程（2.119）形式的一些最简单和最有效的方案是迭代方法。迭代方法是一个近似方法，即它使用初始猜测来开始求解过程，然后迭代以获得解的估计值。方程（2.119）包括\(n\)个方程，无论是线性还是非线性，都有\(n\)个未知数\(\phi\)。矩阵\(\boldsymbol{A}\)包含\(n\times
n\)个系数（尽管在实践中许多的系数都是零）。如果我们假设对角线系数\(a_{ii}(i =
1,\cdots,n)\)不为零（这在有限元方法中是有效的），则可以相对容易地重新排列方程，从而得到一种求解未知值\(\phi\)的方法。考虑 \[\phi^{(k+1)}_1 =
-\frac{1}{a_{11}}(a_{12}\phi^{(k)}_2+a_{13}\phi^{(k)}_3+\cdots+a_{1n}\phi^{(k)}_{n+1})+\frac{b_1}{a_{11}}\\
\phi^{(k+1)}_2 =
-\frac{1}{a_{22}}(a_{21}\phi^{(k)}_1+a_{23}\phi^{(k)}_3+\cdots+a_{2n}\phi^{(k)}_{n+1})+\frac{b_2}{a_{22}}\\
\vdots\\
\phi^{(k+1)}_{n+1} =
-\frac{1}{a_{nn}}(a_{n1}\phi^{(k)}_1+a_{n2}\phi^{(k)}_2+\cdots+a_{n,n-1}\phi^{(k)}_{n})+\frac{b_n}{a_{nn}}\tag{2.120}
\] 其中\(k\)表示迭代索引。
要求解方程（2.113），我们先做一个初始猜测，\(\phi^{(0)}\)，然后将猜测值（例如，\(\phi^{(0)}_i=0\)）代入方程（2.120）的右侧。这产生了一个新的\(\phi\)d的估计值 \(\phi^{(1)}\)，希望它比之前的估计值\(\phi^{(0)}\)更好。我们继续这个过程，获得\(\phi^{(2)},\phi^{(3)},\cdots,\phi^{(k)}\)，直到解收敛。通常通过计算\((k)\)次迭代的相对或绝对误差来确定收敛性，即
\[
\max\frac{|\phi^{(k+1)}_i-\phi^{(k)}_i|}{|\phi^{(k+1)}_i|}<\epsilon\tag{2.121}\]
或 \[ \max
|\phi^{(k+1)}_i-\phi^{(k)}_i|<\epsilon\tag{2.122}\]
这种类型的迭代被称为线性迭代。虽然简单，但该方法可能需要很多次迭代才能在大型问题中实现收敛。为了在小计算机上实现有效的解决方案，需要加速收敛。
矩阵\(\boldsymbol{A}\)还必须满足某些条件才能保证收敛。然而，这不是本文讨论这些条件的范围，除了它们最简单的形式。除了一些评论，我们将满足于声明对于这里考虑的有限元离散化类型，这些条件确实得到满足，并且通常会发生收敛。
高斯-塞德尔迭代方法特别适合解决大型方程组。该方法简单易实现，计算效率高，并且比直接消元方法更不容易受到舍入误差的影响（由迭代次数决定）。然而，必须对这类方案提出警告——存在高斯-塞德尔方法不收敛的情况。这些情况通常发生在矩阵病态时。
由于矩阵求逆计算缓慢且可能需要过多的存储空间，使用矩阵乘法和标量分流来获得方程（2.1119）的解。尽管在迭代方法中有些多余，但这样的操作非常快速。高斯-塞德尔迭代算法可以表示为
\[ \phi_i^{(k+1)} = \frac{1}{a_{ii}} \left(
b_i - \sum_{j=1}^{i-1} a_{ij} \phi_j^{(k+1)} - \sum_{j=i+1}^{n} a_{ij}
\phi_j^{(k)} \right) \] 注意，最更新的 \(\phi_i^{(k+1)}\)
用于迭代过程中，这使得方案更快、更有效。一个收敛的条件要求矩阵是对角占优的，即
\[|a_{ii}| > \sum_{j=1,j\neq i}^{n}
|a_{ij}|,i=1,2,\cdots,n\tag{2.124}\]
其中严格不等式必须至少对一个方程成立。如果这个条件是满足的，无论初始向量的值如何，解都将收敛（这就是为什么许多人使用零作为初始估计的原因）。我们可以将方程（2.123）
用上三角和下三角矩阵表示为，
解决线性代数方程组的第二种方法包括直接求解系统；这些方法基于高斯消除。将系数矩阵\(\boldsymbol{A}\)分解为一个下三角矩阵和一个上三角矩阵的乘积，称为LU分解，是迭代方法的有吸引力的替代方法。可以利用系数矩阵\(\boldsymbol{A}\)的稀疏性质来编写耗时的消除步骤。使用\(LU\)分解的消除程序是直接求解方程组的最流行技术（Atkinson, 1985; Chapra和Canale, 2015）。\(LU\)分解方法是高斯消除的一个变体；通过“分解”\(\boldsymbol{A}\)矩阵为\(\boldsymbol{L}\)和\(\boldsymbol{U}\)矩阵的乘积，可以得到比原始消除更有效的算法。让我们将方程（2.115）重写为 \[\boldsymbol{A}\phi-\boldsymbol{B}=0\tag{2.131}\] 我们可以将方程（2.131）表示为上三角系统的形式，即 \[\boldsymbol{U}\phi-\boldsymbol{D}=0\tag{2.132}\] 如果我们现在用下三角矩阵\(\boldsymbol{L}\)预乘方程(2.132)，并要求得到的系统等于由方程(2.131)给出的原始方程组，我们有 \[\boldsymbol{L}(\boldsymbol{U}\phi-\boldsymbol{D})=\boldsymbol{A}\phi-\boldsymbol{B}\tag{2.133}\] 方程（2.133）有效的条件是 \[\boldsymbol{LU}=\boldsymbol{A}\tag{2.134}\] 和 \[\boldsymbol{LD}=\boldsymbol{B}\tag{2.135}\] 方程（2.134）被称为\(\boldsymbol{A}\)的LU分解。这种类型的分解可以对任何非奇异矩阵\(\boldsymbol{A}\)进行；然而，执行它没有唯一的方法。在文献中可以找到获得这种分解的不同方法。例如，如果将矩阵\(\boldsymbol{L}\)的所有对角单元都设置为1，则会得到一个唯一的分解，这被称为Doolittle约简。如果矩阵\(\boldsymbol{A}\)是对称的，可以以这样的方式进行分解 \[\boldsymbol{L}=\boldsymbol{U}^T\tag{2.136}\] 这被称为Cholesky方法或平方根方法，因为需要平方根操作来获得对角单元\(u_{ii}\)（这将在稍后更详细地讨论）。
方程（2.135）被称为前向替代，并产生向量\(\boldsymbol{D}\)。为了获得解向量\(\phi\)，还需要一个步骤。一旦得到\(\boldsymbol{D}\)，这可以通过求解方程（2.132）来提供，涉及与直接高斯消除中所需的相同类型的后向替代。大多数数值方法教科书都提供了使用这些过程的算法和子程序列表（见Atkinson, 1985; Chapra和Canale, 2015; Conte, 1965; Isaacson和Keller, 1966）。
带状矩阵是一个方阵，除了主对角线中心的一条带外，所有系数都等于零。这样的系统经常出现在微分方程的求解中，特别是在工程和科学问题中。在有限元中，这样的矩阵出现在线性结构和扩散相关问题中。当方程是非线性的，如在流体流动中，可能会出现非对称带状矩阵。在这本入门书中，我们只关心线性对称带状系统方程；然而，与FORTRAN软件一起提供的矩阵求解器将处理对称或非对称带状矩阵。像高斯消除或LU分解这样的直接方法在求解问题的几何形状不规则且方法产生带状结构较差的矩阵时，比迭代方法具有更大的优势。当同一个系统必须为许多不同的右手边求解时，它们也更方便，因为只需要进行一次分解。通过可以非常高效地执行的简单前向和后向替代来获得解。对于单个解的强带状系统方程，这样的方法通常效率较低，并且在存储和操作零值时引入了不必要的计算工作。
对于具有半带宽或半带宽\(\ell\)的对称矩阵，如图9所示，只需要在任何时候在核心中存储顺序为\(\ell\)的上三角部分。总的存储需求是\(\ell[n-(\ell - 1)/2]\)。
矩阵的带状特性以及因此的\(\ell\)的大小，直接取决于节点编号的方式。带宽由所有单元中节点编号的最大差异决定。当我们进一步研究2-D和3-D单元时，这一事实将变得非常重要，并指出了“优化”节点编号的必要性。对于1-D单元，如果节点编号是顺序的，对于线性单元，得到的全局矩阵是三对角的，对于二次单元，是五对角的。然而，读者应该注意，如果节点编号不是顺序的，情况将不是如此。解决对称方程组的最有效和最受欢迎的方法之一是前面提到的Cholesky分解。该算法依赖于这样一个事实：对称矩阵可以分解为彼此转置的上三角和下三角矩阵，即 \[\boldsymbol{A}=\boldsymbol{U}^T\boldsymbol{U}\tag{2.137}\] 其中上标\(T\)表示转置（见附录A）。方程(2.137)可以写成递归形式为 \[l_{ii} = \sqrt{q_{ii} - \sum_{k=1}^{i-1} l_{ik}^2}, \quad i = 1, 2, \ldots, n \tag{2.138}\] 和 \[l_{ij} = \frac{1}{u_{jj}} (a_{ij} - \sum_{k=1}^{j-1} l_{ik} u_{jk}), j < i \tag{2.139}\] 前提是矩阵A是正定的。这个过程显著加快了LU分解步骤。然后可以执行前向和后向替代以得到解向量。在这种方法中，不需要进行枢轴操作以避免除以零。矩阵A的正定性保证了对于所有\(i\)值，\(u_{ii} > 0\)。
2.9 小结
本章的目的是阐述并确立将有限元方法应用于线性偏微分方程。从标量量的简单稳态扩散方程开始，有限元方法被应用于越来越复杂的方程，最终导致在工程和科学问题中常见的1-D时变扩散传输方程。在几个示例问题中推导和使用各种程度的形函数后，可以看到一个通用的有限元过程的演变，该过程可以应用于广泛的范围问题。在整个单元集或网格上的组装过程被展示出来，产生了一组全局矩阵。将全局矩阵系统简化为形式为\(\boldsymbol{A}\phi
=\boldsymbol{B}\)的简单线性（矩阵）方程，其中\(\phi\)是包含未知变量的列向量，允许使用标准矩阵求解程序来求解\(\phi\)。
大多数与有限元方法相关的基本“工具”可以从1-D问题的分析中发展出来。有限元方法的基本要素可以总结如下：
- 将加权残差法应用于控制方程和相关的边界条件。
- 选择合适的形函数和加权因子（线性、二次、三次等）。
- 对\(x\)轴进行离散化，并评估一般有限元域的Galerkin近似。
- 将单元贡献组装成一组全局矩阵。
- 将必要的边界条件数据（由已知值组成）应用于载荷向量。
- 使用合适的矩阵求解程序，从初始条件数据\(\phi_{t_0} = \phi_0\)开始，求解未知值的列向量\(\phi\)，并通过时间\(t^{n+1} = t^n + \Delta t\)进行迭代，直至收敛到稳态条件。
在本文的这一点上，读者应该开始对有限元的数学基础有所了解，这些基础来自于这些简单的1-D示例。对于那些对有限元方法不太熟悉的人，我们建议在继续阅读剩余章节中的更高级材料之前，再次回顾一下最后这几节内容。在下一章中，我们将算法扩展到2-D问题。
