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

文章目录

前言
问题描述
- 问题区域
- 偏微分方程
变分问题(弱形式)
有限元离散
二维线性有限元 : 参考基函数
2D linear finite element : affine mapping
有限元计算
- 刚度矩阵
- 荷载(右端项)
算法(伪代码)
- Stiffness matrix Calculation
- Load vector Calculation
- Deal with the boundary conditions
程序
示例结果

前言

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

问题描述

问题区域

偏微分方程

∇ 2 u ( x , y ) = 0 i n Ω ( 1 ) u = g o n Γ u ( 2 ) q = ∂ u ∂ n ⃗ = 0 o n Γ q ( 3 ) \nabla^2u(x, y) = 0 \ \ \ \ in \ \ \Omega \ \ \ (1) \\ u = g \ \ \ \ \ \ \ \ \ \ on \ \ \ \Gamma_u \ \ \ (2) \\ q = \frac{\partial u}{\partial \vec{n}} = 0 \ \ \ on \ \ \ \Gamma_q \ \ \ (3) \\ ∇2u(x,y)=0 in Ω (1)u=g on Γu (2)q=∂n ∂u=0 on Γq (3)
a). 其中 Ω \Omega Ω 为矩形 ABCD 减去内部矩形 EFGH 的区域。
b). 其中 Γ u \Gamma_u Γu 为矩形 EFGH 的边界。
c). 其中 Γ q \Gamma_q Γq 为矩形 ABCD 的边界。

变分问题(弱形式)

任取函数 v ( x , y ) ∈ V v(x, y) \in V v(x,y)∈V 乘在方程 (1) 两端并在区域 Ω \Omega Ω 上积分得到如下方程：
∫ Ω v ( x , y ) ∗ ∇ 2 u ( x , y ) d s = 0 ( 4 ) \int _{\Omega} v(x,y)*\nabla^2u(x, y) ds= 0 \ \ \ \ \ \ \ (4) ∫Ωv(x,y)∗∇2u(x,y)ds=0       (4)
对方程 (4) 运用 Green 公式(可以参考博客边界元方法(一) 中 二重积分的分步积分公式 章节)可得到如下方程：
∫ Ω ∇ u ⋅ ∇ v d s − ∫ ∂ Ω ( ∇ u ⋅ n ) v d w = 0 ( 5 ) \int_{\Omega} \nabla u \cdot \nabla vds - \int_{\partial \Omega} (\nabla u \cdot \mathbf{n})vdw= 0 \ \ \ \ \ \ \ (5) ∫Ω∇u⋅∇vds−∫∂Ω(∇u⋅n)vdw=0       (5)
因 ∂ Ω = Γ u + Γ q \partial \Omega = \Gamma_u + \Gamma_q ∂Ω=Γu+Γq 方程 (5) 等价于下面的方程：
∫ Ω ∇ u ⋅ ∇ v d s − ∫ Γ q ( ∇ u ⋅ n ) v d w − ∫ Γ u ( ∇ u ⋅ n ) v d w = 0 ( 6 ) \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) ∫Ω∇u⋅∇vds−∫Γq(∇u⋅n)vdw−∫Γu(∇u⋅n)vdw=0       (6)
因等式 (2) 我们可以选择最初的 v ( x , y ) = 0 o n Γ u v(x, y) = 0 \ \ on \ \ \ \Gamma_u v(x,y)=0  on   Γu 则方程 (6) 边成如下形式：
∫ Ω ∇ u ⋅ ∇ v d s − ∫ Γ q ( ∇ u ⋅ n ) v d w = 0. ( 7 ) \int_{\Omega} \nabla u \cdot \nabla vds - \int_{\Gamma_q} (\nabla u \cdot \mathbf{n})vdw = 0. \ \ \ \ \ \ \ (7) ∫Ω∇u⋅∇vds−∫Γq(∇u⋅n)vdw=0.       (7)
再将等式 (3) 带入到等式 (7) 中可以得到如下的最终弱形式：
T h e n t h e w e a k f o r m u l a t i o n i s t o f i n d u ∈ U s u c h t h a t Then \ \ the \ \ weak \ \ formulation \ \ is \ \ to \ \ find \ \ u \in U \ \ such \ \ that Then  the  weak  formulation  is  to  find  u∈U  such  that
∫ Ω ∇ u ⋅ ∇ v d s = 0 ( 8 ) \int_{\Omega} \nabla u \cdot \nabla vds = 0 \ \ \ \ \ \ \ \ \ \ (8) ∫Ω∇u⋅∇vds=0          (8)
f o r a n y v ∈ V for \ \ any \ \ v \in V for  any  v∈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 ∈ U h Γ u_h \in U_h^{\Gamma} uh∈UhΓ such that
a ( u h , v h ) = ∫ Ω ∇ u h ⋅ ∇ v h d s = 0 ( 11 ) a(u_h, v_h) = \int_{\Omega} \nabla u_h \cdot \nabla v_hds = 0 \ \ \ \ \ \ (11) a(uh,vh)=∫Ω∇uh⋅∇vhds=0 (11)
for any v h ∈ V h Γ v_h \in V_h^{\Gamma} vh∈VhΓ.

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 P_1 P1 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 E_i Ei 为离散的网格单元，以下用三角形单元介绍。

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

以 L a g r a n g e Lagrange Lagrange 型的节点基函数为例。考虑参考三角形单元 E ^ = △ A ^ 1 B ^ 1 C ^ 1 \hat{E}=\bigtriangleup \hat{A}_1\hat{B}_1\hat{C}_1 E^=△A^1B^1C^1 上的 2D 线性基函数，其中 A ^ 1 = ( 0 , 0 ) \hat{A}_1=(0, 0) A^1=(0,0), B ^ 1 = ( 1 , 0 ) \hat{B}_1=(1, 0) B^1=(1,0), C ^ 1 = ( 0 , 1 ) \hat{C}_1=(0, 1) C^1=(0,1)。
定义基函数如下：
ψ ^ j ( x ^ , y ^ ) = a j x ^ + b j y ^ + c j j = 1 , 2 , 3 , ( 15 ) \hat{\psi}_j(\hat{x}, \hat{y}) = a_j\hat{x} + b_j\hat{y} + c_j \ \ \ \ \ \ j = 1, 2, 3,\ \ \ \ \ (15) ψ^j(x^,y^)=ajx^+bjy^+cj j=1,2,3, (15)
满足如下条件：
ψ ^ j ( A ^ i ) = δ i j = { 0 i f j ≠ i 1 i f j = i ( 16 ) \hat{\psi}_j(\hat{A}_i) = \delta _{ij} = \left\{\begin{matrix} 0 & if \ \ \ j \neq i \\ 1 & if \ \ \ j = i \\ \end{matrix}\right. \ \ \ \ \ (16) ψ^j(A^i)=δij={01if j=iif j=i (16)
求解方程 (15) (16) 我们可以得到如下参考单元上的 2D 线性基函数：
ψ ^ 1 ( x ^ , y ^ ) = − x ^ − y ^ + 1 ψ ^ 2 ( x ^ , y ^ ) = x ^ ψ ^ 3 ( x ^ , y ^ ) = y ^ \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} ψ^1(x^,y^)=−x^−y^+1ψ^2(x^,y^)=x^ψ^3(x^,y^)=y^

2D linear finite element : affine mapping

任意三角形单元 E = △ A 1 B 1 C 1 E=\bigtriangleup A_1B_1C_1 E=△A1B1C1 与参考三角形单元 E ^ = △ A ^ 1 B ^ 1 C ^ 1 \hat{E}=\bigtriangleup \hat{A}_1\hat{B}_1\hat{C}_1 E^=△A^1B^1C^1 的映射。
设：
A i = ( x i , y i ) , i = 1 , 2 , 3. A_i = (x_i, y_i), \ \ \ \ \ i = 1, 2, 3. Ai=(xi,yi), i=1,2,3.
consider the affine mapping
( x y ) = ( A 2 − A 1 , A 3 − A 1 ) ( x ^ y ^ ) + A 1 = ( x 2 − x 1 x 3 − x 1 y 2 − y 1 y 3 − y 1 ) ( x ^ y ^ ) + ( x 1 y 1 ) . \begin{aligned} \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{aligned} (xy)=(A2−A1,A3−A1)(x^y^)+A1=(x2−x1y2−y1x3−x1y3−y1)(x^y^)+(x1y1).
得到：
x ^ = ( y 3 − y 1 ) ( x − x 1 ) − ( x 3 − x 1 ) ( y − y 1 ) ( x 2 − x 1 ) ( y 3 − y 1 ) − ( x 3 − x 1 ) ( y 2 − y 1 ) , y ^ = ( y 2 − y 1 ) ( x − x 1 ) − ( x 2 − x 1 ) ( y − y 1 ) ( x 2 − x 1 ) ( y 3 − y 1 ) − ( x 3 − x 1 ) ( y 2 − y 1 ) . \begin{aligned} \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{aligned} x^y^=(x2−x1)(y3−y1)−(x3−x1)(y2−y1)(y3−y1)(x−x1)−(x3−x1)(y−y1),=(x2−x1)(y3−y1)−(x3−x1)(y2−y1)(y2−y1)(x−x1)−(x2−x1)(y−y1).

定义 J a c o b i m a t r i x Jacobi \ \ \ matrix Jacobi matrix :
J = ( x 2 − x 1 x 3 − x 1 y 2 − y 1 y 3 − y 1 ) 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=(x2−x1y2−y1x3−x1y3−y1)
则有
∣ J ∣ = ( x 2 − x 1 ) ( y 3 − y 1 ) − ( x 3 − x 1 ) ( y 2 − y 1 ) , |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), ∣J∣=(x2−x1)(y3−y1)−(x3−x1)(y2−y1),
x ^ = ( y 3 − y 1 ) ( x − x 1 ) − ( x 3 − x 1 ) ( y − y 1 ) ∣ J ∣ y ^ = − ( y 2 − y 1 ) ( x − x 1 ) + ( x 2 − x 1 ) ( y − y 1 ) ∣ J ∣ \begin{aligned} \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{aligned} x^y^=∣J∣(y3−y1)(x−x1)−(x3−x1)(y−y1)=∣J∣−(y2−y1)(x−x1)+(x2−x1)(y−y1)

有限元计算

选择有限元空间 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 N_b Nb 为有限元节点数目(注与网格的几何节点数不一定相同)。
因此 u h ∈ U h Γ u_h \in U_h^{\Gamma} uh∈UhΓ 则有
u h = ∑ i = 1 N b u i ψ i ( 17 ) u_h=\sum_{i=1}^{N_b}u_i\psi_i \ \ \ \ \ \ \ (17) uh=i=1∑Nbuiψi (17)
且满足 u h = g ( x , y ) ∈ Γ u u_h = g \ \ \ \ (x, y) \in \Gamma_u uh=g (x,y)∈Γu

综上，选择 " t e s t f u n c t i o n test \ \ \ function test function" v h = ψ i ( i = 1 , . . . , N t ) v_h = \psi_i(i=1,...,N_t) vh=ψi(i=1,...,Nt), 则离散有限元方程为：
∫ Ω ∇ ( ∑ j = 1 N b u j ψ j ) ⋅ ∇ ψ i d x d y = 0 ⇒ ∑ j = 1 N b u j [ ∫ Ω ∇ ψ j ⋅ ∇ ψ i d x d y ] = 0 , i = 1 , ⋯ , N t ( 18 ) \begin{aligned} & \int_{\Omega} \nabla\left(\sum_{j=1}^{N_{b}} u_{j} \psi_{j}\right) \cdot \nabla \psi_{i} d x d y=0 \\ \Rightarrow & \sum_{j=1}^{N_{b}} u_{j}\left[\int_{\Omega} \nabla \psi_{j} \cdot \nabla \psi_{i} d x d y\right]=0, \quad i=1, \cdots, N_{t} \end{aligned} \ \ \ \ \ (18) ⇒∫Ω∇(j=1∑Nbujψj)⋅∇ψidxdy=0j=1∑Nbuj[∫Ω∇ψj⋅∇ψidxdy]=0,i=1,⋯,Nt (18)

刚度矩阵

A = [ a i j ] = [ ∫ Ω ∇ ψ j ⋅ ∇ ψ i d x d y ] i , j = 1 i = N t , j = N b A=\left[a_{i j}\right]=\left[\int_{\Omega} \nabla \psi_{j} \cdot \nabla \psi_{i} d x d y\right]_{i, j=1}^{i=N_{t}, j=N_{b}} A=[aij]=[∫Ω∇ψj⋅∇ψidxdy]i,j=1i=Nt,j=Nb

荷载(右端项)

b ⃗ = [ b i ] i = 1 i = N t = 0 \vec{b}=\left[b_{i}\right]_{i=1}^{i=N_{t}}=0 b =[bi]i=1i=Nt=0

算法(伪代码)

Stiffness matrix Calculation

Load vector Calculation

Deal with the boundary conditions

程序

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

示例结果