有限元分析简介及伽辽金法

1.1 什么是有限元方法

有限元方法(Finite Element Method, FEM)是一种求解由偏微分方程描述或可表示为泛函极小化问题的数值方法。感兴趣的域被表示为有限单元(finite elements)的集合。有限元中的逼近函数是根据所求物理场的节点值确定的。FEM将一个连续的物理问题转化为节点值未知的离散化有限元问题，并得到一个线性方程组，求解该方程组就可以获得待求的物理量。有限元内部的值可以使用节点值恢复。

值得一提的是，FEM的两个特点：

在有限元上的物理场的分段近似提供了很好的精度，即使用简单的近似函数(增加单元的数量，我们可以达到任意精度)。
局部逼近将导致离散问题的方程组稀疏，这有助于解决具有大量节点未知数的问题。

1.2 FEM是如何工作的

为了说明有限元方法是如何工作的，我们列出如下几个有限元的主要求解步骤

离散化。第一步是将求解域划分为有限单元。有限元网格通常由预处理程序生成。网格的描述由几个数组组成，其中主要是节点坐标和单元的连通性。
选择插值函数。插值函数用于对单元上的场变量进行插值，这里通常选择多项式作为插值函数。多项式的次数取决于单元上的节点数。
寻找单元属性。建立待求函数的结点值与其他参数联系起来的矩阵方程，此步骤有多种方法可以使用，但最方便的方法是伽辽金方法(Galerkin method)。
组装单元方程。要找到整个求解区域的全局方程组，必须将所有的单元方程组合在一起。在组装过程中需要用到单元之间的连通性。同时，在求解之前还需要施加边界条件（边界条件不在单元方程中考虑）。
求解全局方程组。有限元整体方程组具有典型的稀疏、对称和正定性质，可采用直接法和迭代法来求解该方程组。待求函数在结点处的值被当做最终的解。
计算其他量。对于很多情况，我们需要计算其他的量，例如在力学问题中，除了位移外，还需要考虑应变和应力，这些位移是在求解整体方程组后得到的。

1.3 有限元方程的表述

有很多方法可以将问题的物理形式转换成有限元离散后的形式。如果物理问题是由偏微分方程描述的，则经常使用的是伽辽金方法(Galerkin method)。如果文理问题是由泛函最小化描述的，则经常使用变分公式(variational formulation)将其转换为有限元方程。

1.3.1 伽辽金方法(Galerkin method)

让我们先从简单的一维问题开始，介绍使用伽辽金方法的有限元公式。假设我们需要求解下面方程的数值解
ad2udx2+b=0,0≤x≤1(1)a\frac{d^2u}{dx^2}+b=0,\quad 0 \leq x \leq 1 \tag{1} adx2d2u+b=0,0≤x≤1(1)
边界条件为
u∣x=0=0adudx∣x=2L=R\begin{array}{l} \left.u\right|_{x=0}=0 \\ \left.a \frac{d u}{d x}\right|_{x=2 L}=R \end{array} u∣x=0=0adxdu∣∣x=2L=R
其中uuu是待求解量。我们将使用下图所示的两个线性一维有限元来解决这个问题。

首先我们考虑图1.1右半部分的有限元，该单元拥有两个节点，则近似函数u(x)u(x)u(x)可以表示为
u=N1u1+N2u2=[N]{u}[N]=[N1N2]{u}={u1u2}\begin{array}{l} u=N_1u_1+N_2u_2=[N]\{u\} \\ [N]=\left[\begin{array}{c}N_1 & N_2\end{array}\right] \\ \{u\}=\left\{\begin{array}{c}u_1 & u_2\end{array}\right\} \end{array} u=N1u1+N2u2=[N]{u}[N]=[N1N2]{u}={u1u2}
这里的NiN_iNi也称为形函数(shape functions)
N1=1−x−x1x2−x1N2=x−x1x2−x1N_1=1-\frac{x-x_1}{x_2-x_1} \\ N_2=\frac{x-x_1}{x_2-x_1} N1=1−x2−x1x−x1N2=x2−x1x−x1
使用这些形函数在结点出的值来插值u(x)u(x)u(x)。我们需要从全局离散方程中求解出未知的结点值u1u_1u1和u2u_2u2。

将节点处的值和形函数带入到式(1)的微分方程中，得到如下的近似形式
ad2dx2[N]{u}+b=ψ(2)a\frac{d^2}{dx^2}[N]\{u\}+b=\psi \tag{2} adx2d2[N]{u}+b=ψ(2)
由于式(2)是有限元中的近似形式，所以残差ψ\psiψ是一个非零的数。伽辽金方法通过将式(2)的项乘以形函数，在单元上积分并令其等于0最小化残差：
∫x1x2[N]Tad2dx2[N]{u}dx+∫x1x2[N]Tbdx=0\int_{x_{1}}^{x_{2}}[N]^{T} a \frac{d^{2}}{d x^{2}}[N]\{u\} d x+\int_{x_{1}}^{x_{2}}[N]^{T} b d x=0 ∫x1x2[N]Tadx2d2[N]{u}dx+∫x1x2[N]Tbdx=0
对其进行分部积分，得到该单元上微分方程的一种离散形式
∫x1x2[dNdx]Ta[dNdx]dx{u}−∫x1x2[N]Tbdx−{01}adudx∣x=x2+{10}adudx∣x=x1=0(3)\int_{x_{1}}^{x_{2}}\left[ \frac{dN}{dx} \right]^{T} a \left[ \frac{dN}{dx}\right] d x\{u\}-\int_{x_{1}}^{x_{2}}[N]^{T} b d x - \left\{ \begin{array}{c} 0 \\ 1 \end{array} \right\}a \left. \frac{du}{dx} \right|_{x=x_2}+ \left\{ \begin{array}{c} 1 \\ 0 \end{array} \right\}a \left. \frac{du}{dx} \right|_{x=x_1} =0 \tag{3} ∫x1x2[dxdN]Ta[dxdN]dx{u}−∫x1x2[N]Tbdx−{01}adxdu∣∣∣∣x=x2+{10}adxdu∣∣∣∣x=x1=0(3)
为了表示方便，我们使用如下几个记号

[k]=∫x1x2[dNdx]Ta[dNdx]dx{f}=∫x1x2[N]Tbdx+{01}adudx∣x=x2−{10}adudx∣x=x1\begin{aligned} &{[k]=\int_{x_{1}}^{x_{2}}\left[\frac{d N}{d x}\right]^{T} a\left[\frac{d N}{d x}\right] d x} \\ & \{f\}=\int_{x_{1}}^{x_{2}}[N]^{T} b d x+\left.\left\{\begin{array}{l} 0 \\ 1 \end{array}\right\} a \frac{d u}{d x}\right|_{x=x_{2}}-\left.\left\{\begin{array}{l} 1 \\ 0 \end{array}\right\} a \frac{d u}{d x}\right|_{x=x_{1}} \end{aligned} [k]=∫x1x2[dxdN]Ta[dxdN]dx{f}=∫x1x2[N]Tbdx+{01}adxdu∣∣∣∣x=x2−{10}adxdu∣∣∣∣x=x1

则式子(3)可以写成
[k]{u}={f}\boxed{{[k]\{u\}=\{f\}}} [k]{u}={f}
在固体力学中，[k][k][k]称为刚度矩阵(stiffness matrix)；{f}\{f\}{f}称为载荷向量(load vector)。在上述我们考虑的偏微分问题中，我们分别计算两个长度为LLL单元的刚度矩阵和载荷向量。
[k1]=[k2]=aL[1−1−11][k_1]=[k_2]=\frac{a}{L} \left[ \begin{array}{c} 1 & -1 \\ -1 & 1 \end{array} \right] [k1]=[k2]=La[1−1−11]

{f1}=bL2{11},{f2}=bL2{11}+{0R}\{f_1\}= \frac{bL}{2} \left\{ \begin{array}{c} 1 \\ 1 \end{array} \right\}, \quad \{f_2\}= \frac{bL}{2} \left\{ \begin{array}{c} 1 \\ 1 \end{array} \right\} + \left\{ \begin{array}{c} 0 \\ R \end{array} \right\} {f1}=2bL{11},{f2}=2bL{11}+{0R}

对于整个求解域上2个单元，3个结点的全局方程组可以通过组装单元方程组获得。在此实例中，第二个结点时单元共用结点。组装全局方程组为
aL[1−10−12−10−11]{u1u2u3}=bL2{121}+{00R}\frac{a}{L} \left[ \begin{array}{c} 1 & -1 & 0 \\ -1 & 2 & -1 \\ 0 & -1 & 1 \end{array} \right] \left\{ \begin{array}{c} u_1 \\ u_2 \\ u_3 \end{array} \right\} =\frac{bL}{2} \left\{ \begin{array}{c} 1 \\ 2 \\ 1 \end{array} \right\} + \left\{ \begin{array}{c} 0 \\ 0 \\ R \end{array} \right\} La⎣⎡1−10−12−10−11⎦⎤⎩⎨⎧u1u2u3⎭⎬⎫=2bL⎩⎨⎧121⎭⎬⎫+⎩⎨⎧00R⎭⎬⎫
施加边界条件u(x=0)=0u(x=0)=0u(x=0)=0后，最终的全局方程为
aL[10002−10−11]{u1u2u3}=bL2{021}+{00R}\frac{a}{L} \left[ \begin{array}{c} 1 & 0 & 0 \\ 0 & 2 & -1 \\ 0 & -1 & 1 \end{array} \right] \left\{ \begin{array}{c} u_1 \\ u_2 \\ u_3 \end{array} \right\} =\frac{bL}{2} \left\{ \begin{array}{c} 0 \\ 2 \\ 1 \end{array} \right\} + \left\{ \begin{array}{c} 0 \\ 0 \\ R \end{array} \right\} La⎣⎡10002−10−11⎦⎤⎩⎨⎧u1u2u3⎭⎬⎫=2bL⎩⎨⎧021⎭⎬⎫+⎩⎨⎧00R⎭⎬⎫
结点处的值uiu_iui可以通过求解上面的线性方程组得到。uuu在有限元内任一点处的值可以通过形函数计算。当a=1,b=1,L=1,R=1a=1,b=1,L=1,R=1a=1,b=1,L=1,R=1时有限元的解如下图所示。

在此例中，精确解是二次函数，而有限元方法中在每个单元上使用的是简单的线性函数。通过增加单元的个数或使用更加复杂的形函数可以得到更加精确的数值解。值得注意的是，在结点处有限元方法得到的解是精确的（仅针对此例中）。如果精确解是二次函数，则有限元方法使用线性形函数可在结点处获得精确值；而如果精确解是三次函数，则有限元方法使用二次形函数可在结点处获得精确值，以此类推。