有限元方法基础-以二维拉普拉斯方程为例(附程序)

前言

本文从零介绍有限元方法，包括每一步的数学推导，同时附程序开发指南。可以方便新手入门。

问题描述

问题区域

偏微分方程

$$\begin{gathered} {\nabla }^{2}u(x,y)=0in\Omega (1) \\ u=gon{\Gamma }_u(2) \\ q=\frac{\partial u}{\partial \vec{n}}=0on{\Gamma }_q(3) \end{gathered}$$
a). 其中$\Omega$ 为矩形 ABCD 减去内部矩形 EFGH 的区域。
b). 其中 ${\Gamma }_u$ 为矩形 EFGH 的边界。
c). 其中 ${\Gamma }_q$ 为矩形 ABCD 的边界。

变分问题(弱形式)

任取函数 $v(x,y)\in V$ 乘在方程 (1) 两端并在区域 $\Omega$ 上积分得到如下方程：
$${\int }_{\Omega }v(x,y)\ast {\nabla }^{2}u(x,y)ds=0(4)$$
对方程 (4) 运用 Green 公式(可以参考博客 边界元方法(一) 中 二重积分的分步积分公式 章节)可得到如下方程：
$${\int }_{\Omega }\nabla u\cdot \nabla vds-{\int }_{\partial \Omega }(\nabla u\cdot \mathbf{n})vdw=0(5)$$
因 $\partial \Omega ={\Gamma }_u+{\Gamma }_q$ 方程 (5) 等价于下面的方程：
$${\int }_{\Omega }\nabla u\cdot \nabla vds-{\int }_{{\Gamma }_q}(\nabla u\cdot \mathbf{n})vdw-{\int }_{{\Gamma }_u}(\nabla u\cdot \mathbf{n})vdw=0(6)$$
因等式 (2) 我们可以选择最初的 $v(x,y)=0on{\Gamma }_u$ 则方程 (6) 边成如下形式：
$${\int }_{\Omega }\nabla u\cdot \nabla vds-{\int }_{{\Gamma }_q}(\nabla u\cdot \mathbf{n})vdw=0.(7)$$
再将等式 (3) 带入到等式 (7) 中可以得到如下的最终弱形式：
$Thentheweakformulationistofindu\in Usuchthat$
$${\int }_{\Omega }\nabla u\cdot \nabla vds=0(8)$$
$foranyv\in V$
其中
U : = { u ∈ C 0 ( Ω ) ∣ u = g , ( x , y ) ∈ Γ u }      ( 9 )   V : = { v ∈ C 0 ( Ω ) ∣ v = 0 , ( x , y ) ∈ Γ u }      ( 10 ) U := \{u \in C^0(\Omega) | u = g, \quad (x,y) \in \Gamma_u \} \ \ \ \ (9) \\\ V := \{v \in C^0(\Omega) | v = 0, \quad (x,y) \in \Gamma_u \} \ \ \ \ (10) U:={u∈C0(Ω)∣u=g,(x,y)∈Γu​}    (9) V:={v∈C0(Ω)∣v=0,(x,y)∈Γu​}    (10)

有限元离散

The Galerkin formulation : find $u_h\in U_h^{\Gamma }$ such that
$$a(u_h,v_h)={\int }_{\Omega }\nabla u_h\cdot \nabla v_{h}ds=0(11)$$
for any $v_h\in V_h^{\Gamma }$.

V h : = { u h ∈ C ( T h ) ∣ u h ∣ E i ∈ P 1 ∀ E i ∈ T h }     ( 12 )     U h Γ : = { ψ h ∈ V h ∣ ψ h ∣ Γ u = g }     ( 13 )     V h Γ : = { ϕ h ∈ V h ∣ ψ h ∣ Γ u = 0 }     ( 14 ) V_h := \{ u_h \in C(T_h) | \quad u_h | _{E_i} \in P_1 \quad \forall E_i \in T_h \} \ \ \ (12) \\\ \\\ U_h^{\Gamma} := \{ \psi_h \in V_h | \quad \psi_h |_{\Gamma_u} = g \} \ \ \ (13) \\\ \\\ V_h^{\Gamma} := \{ \phi_h \in V_h | \quad \psi_h |_{\Gamma_u} = 0 \} \ \ \ (14) Vh​:={uh​∈C(Th​)∣uh​∣Ei​​∈P1​∀Ei​∈Th​}   (12)  UhΓ​:={ψh​∈Vh​∣ψh​∣Γu​​=g}   (13)  VhΓ​:={ϕh​∈Vh​∣ψh​∣Γu​​=0}   (14)
Where $P_1$ is a collection of linear polynomial spaces.

T h : = { E n , n = 1 , ⋯   , n u m M e s h E l e m e n t s } T_h:=\{ E_n, \quad n = 1, \cdots, numMeshElements\} Th​:={En​,n=1,⋯,numMeshElements} ，其中 $E_i$ 为离散的网格单元，以下用三角形单元介绍。

二维线性有限元 : 参考基函数

以 $Lagrange$ 型的节点基函数为例。考虑参考三角形单元$\hat{E}=△{\hat{A}}_1{\hat{B}}_1{\hat{C}}_1$ 上的 2D 线性基函数，其中 ${\hat{A}}_1=(0,0)$, ${\hat{B}}_1=(1,0)$, ${\hat{C}}_1=(0,1)$。
定义基函数如下：
$${\hat{\psi }}_j(\hat{x},\hat{y})=a_j\hat{x}+b_j\hat{y}+c_{j}j=1,2,3,(15)$$
满足如下条件：
$${\hat{\psi }}_j({\hat{A}}_i)={\delta }_{ij}=\{\begin{array}{cc}0& ifj\ne i\\ 1& ifj=i\end{array}(16)$$
求解方程 (15) (16) 我们可以得到如下参考单元上的 2D 线性基函数：
$$\begin{gathered} {\hat{\psi }}_1(\hat{x},\hat{y})=-\hat{x}-\hat{y}+1 \\ {\hat{\psi }}_2(\hat{x},\hat{y})=\hat{x} \\ {\hat{\psi }}_3(\hat{x},\hat{y})=\hat{y} \end{gathered}$$

2D linear finite element : affine mapping

任意三角形单元 $E=△A_1{B}_1{C}_1$ 与 参考三角形单元$\hat{E}=△{\hat{A}}_1{\hat{B}}_1{\hat{C}}_1$ 的映射。
设：
$$A_i=(x_i,y_i),i=1,2,3.$$
consider the affine mapping
$$\begin{array}{rl}\left(\begin{array}{c}x\\ y\end{array}\right)& =\left(A_2-A_1,A_3-A_1\right)\left(\begin{array}{c}\hat{x}\\ \hat{y}\end{array}\right)+A_1\\ & =\left(\begin{array}{ll}{x}_2-x_1& x_3-x_1\\ y_2-y_1& y_3-y_1\end{array}\right)\left(\begin{array}{l}\hat{x}\\ \hat{y}\end{array}\right)+\left(\begin{array}{l}{x}_1\\ y_1\end{array}\right).\end{array}$$
得到：
$$\begin{array}{rl}\hat{x}& =\frac{\left(y_3-y_1\right)\left(x-x_1\right)-\left(x_3-x_1\right)\left(y-y_1\right)}{\left(x_2-x_1\right)\left(y_3-y_1\right)-\left(x_3-x_1\right)\left(y_2-y_1\right)},\\ \hat{y}& =\frac{\left(y_2-y_1\right)\left(x-x_1\right)-\left(x_2-x_1\right)\left(y-y_1\right)}{\left(x_2-x_1\right)\left(y_3-y_1\right)-\left(x_3-x_1\right)\left(y_2-y_1\right)}.\end{array}$$

定义 $Jacobimatrix$ :
$$J=\left(\begin{array}{ll}{x}_2-x_1& x_3-x_1\\ y_2-y_1& y_3-y_1\end{array}\right)$$
则有
$$|J|=\left(x_2-x_1\right)\left(y_3-y_1\right)-\left(x_3-x_1\right)\left(y_2-y_1\right),$$
$$\begin{array}{rl}\hat{x}& =\frac{\left(y_3-y_1\right)\left(x-x_1\right)-\left(x_3-x_1\right)\left(y-y_1\right)}{|J|}\\ \hat{y}& =\frac{-\left(y_2-y_1\right)\left(x-x_1\right)+\left(x_2-x_1\right)\left(y-y_1\right)}{|J|}\end{array}$$

有限元计算

选择有限元空间 U h = s p a n { ψ i } i = 1 N b U_h =span\{ \psi_i\}_{i=1}^{N_b} Uh​=span{ψi​}i=1Nb​​ 其中 $N_b$ 为有限元节点数目(注与网格的几何节点数不一定相同)。
因此 $u_h\in U_h^{\Gamma }$ 则有
$$u_h=\sum _{i=1}^{N_b}{u}_i{\psi }_i(17)$$
且满足 $u_h=g(x,y)\in {\Gamma }_u$

综上， 选择 " $testfunction$" $v_h={\psi }_i(i=1,...,N_t)$, 则离散有限元方程为：
$$\begin{array}{rl}& {\int }_{\Omega }\nabla \left(\sum _{j=1}^{N_b}{u}_j{\psi }_j\right)\cdot \nabla {\psi }_{i}dxdy=0\\ \Rightarrow & \sum _{j=1}^{N_b}{u}_j\left[{\int }_{\Omega }\nabla {\psi }_j\cdot \nabla {\psi }_{i}dxdy\right]=0,i=1,\cdots ,N_t\end{array}(18)$$

刚度矩阵

$$A=\left[a_{ij}\right]={\left[{\int }_{\Omega }\nabla {\psi }_j\cdot \nabla {\psi }_{i}dxdy\right]}_{i,j=1}^{i=N_t,j=N_b}$$

荷载(右端项)

$$\vec{b}={\left[b_i\right]}_{i=1}^{i=N_t}=0$$

算法(伪代码)

Stiffness matrix Calculation

Load vector Calculation

Deal with the boundary conditions

程序

具体程序介绍可以参考我的另一篇博客 ：有限元方法求解二维拉普拉斯方程C++实现

示例结果