有限元方法入门(2):二维算例(三角剖分)

上一篇 一维有限元方法 谈到,在一维空间上,由于网格剖分形式单一,所以FEM与FDM相比,并没有太大的优势。而在二维空间中,由于FDM对于规则网格的依赖性较高,FEM的优势便显现出来。

同样是选择椭圆方程,来应用二维空间上的有限元方法,考虑最常见的Poisson方程,如下:

其中 $\Omega\in \mathbb{R}^2$ 是区域上的一个有界开集,$\Gamma$ 是区域 $\Omega$ 的边界。

我们选取具体的二维Poisson方程,来讲一讲有限元方法的具体步骤:

1. 构造变分形式

对于任意的 $v \in H^1_0$,乘以 $-\nabla^2 u=f$ 的两侧,运用Green公式,并带入边界条件可以得到:

即得到边值的变分形式 $a(u, v)=(f,v)$。

2. 区域三角剖分

一般来说,三角区域剖分没有特殊的限制。但为构造最接近于规则化的网格,则需要符合Delaunay三角剖分准则,即:

  1. 唯一性:在Delaunay三角形网中任一三角形的外接圆范围内不会有其它点存在。
  2. 最大化最小角:在散点集可能形成的三角剖分中,Delaunay三角剖分所形成的三角形的最小角最大。从这个意义上讲,Delaunay三角网是“最接近于规则化的“的三角网。具体的说是指在两个相邻的三角形构成凸四边形的对角线,在相互交换后,六个内角的最小角不再增大。

为简单起见,在矩形区域 $[0,1]\times [0,1]$ 上,创建规则的直角三角网格,具体形式如下:

将 $x, y$ 区间均匀分割为 $n$ 份,则可得到 $N=(n+1)^2$ 个点,$3n^2+2n$ 条边,以及 $M=2n^2$ 个三角单元。

构造节点坐标矩阵 $P$,边元节点序号矩阵 $E$,单元节点序号矩阵 $U$,边界节点布尔矩阵 $PBnd$,单元边元序号矩阵 $UEdg$

Tips: 具体的剖分方式不再详述,参考 mesh_2d.py,实现方式较为简单。

3. 建立有限元空间

有限元空间

在第一步构造变分形式的时候,我们提到选取 $v\in H^1_0$,构造形如 $a(u, v)=(f, v)$ 的弱形式,而要进行有限元计算,则需要构造Sobolev空间 $H^1_0$ 的有限维子空间 $V_h$(其中 $u_h\in V_h$ 是次数不超过1的多项式,在区域 $[0, 1]\times [0, 1]$ 上属于函数空间 $H^1_0$,且有 $u_h(x, y)=0, x,y\in\partial \Omega$)。

最简单的分段线性函数空间 $V_h$ 是由分段线性函数组成的,它由节点上的一组值

按照一次Lagrange型插值公式

所确定,称为线性有限元空间。

试探函数由 $N$ 个节点的值决定,故 $u_h$ 的自由度为 $N$,也即 $V_h$ 是 $N$ 维线性空间。

PS: 这里在构造试探函数空间 $V_h$ 的时候,仅说明了 $V_h\subset H^1_0$,而未带入边界条件

基底函数

与一维情形相似的,我们在二维区域上同样创建一种“山形函数”,即节点的基函数在该节点处取值为1,在其他节点处取值为0,数学表示为:

同样的,$\phi_0,\ \phi_1,\ \cdots,\ \phi_{N-1},\ \phi_N$ 线性无关,且组成 $V_h$ 的基底,使得 $V_h=\mathrm{span}\{\phi_i\}^N_{i=1}$,且对于任意的 $u_h\in V_h$,都有 $u_h=\sum^N_{i=1}u_i \phi_i,\ u_i=u_h(x_i, y_i)$。

离散化变分问题

将试探函数 $u_h=\sum^N_{i=0}u_i \phi_i$ 和检验函数 $v_h\in \{\phi_i\}^N_{i=1}$,带入到变分形式中,得到

积分以内积形式表示,则方程表现为:

构造单元函数

由于 $\phi_i$ 和 $\phi_j$ 都是区域 $\Omega$ 上的分片线性函数,所以二者的内积 $(\phi_i, \phi_j)$ 需分区域计算,而不同节点上基函数的非零区域并不相同,如果按照节点-节点计算内积,则计算过程十分繁琐。同样的,在编写程序计算时,考虑到内存空间占用问题,以及为方便多线程处理,需要高效的算法。

可以考虑,将空间基函数转换到每个单元上,从单元基函数的角度,进行计算。

基于三角单元的三个顶点,构造简单的一次Lagrange型函数 $f(x,y)=ax+by+c$,即有:

其中

即表现为三维空间中的三个平面:

Tips: 关于三角单元上内积的数值方法,在 附录 中会详细谈到。

4. 建立有限元方程组

组装载荷矩阵F

$F$ 称为总载荷列阵

基于每一单元,计算载荷列阵,记为 $F_i$,称为单元载荷列阵

单元载荷矩阵组装总载荷矩阵:

1
2
3
4
输入: 单元载荷矩阵 Fi, 总载荷矩阵 F, 当前单元节点序号矩阵 Ui
for i = 0 -> 2 do
K[Ui[i]] <- K[Ui[i]] + Fi[i]
end for

组装刚度矩阵K

$K$ 称为总刚度矩阵

基于每一单元,计算刚度矩阵,记为 $K_i$,称为单元刚度矩阵

由单元刚度矩阵组装总刚度矩阵:

1
2
3
4
5
6
输入: 单元刚度矩阵 Ki, 总刚度矩阵 K, 当前单元节点序号矩阵 Ui
for i = 0 -> 2 do
for j = 0 -> 2 do
K[Ui[j]][Ui[i]] <- K[Ui[j]][Ui[i]] + Ki[i][j]
end for
end for

位移列阵U

$U$ 称为总位移列阵

$U$ 中的每一个元素 $u_i$ 即为函数 $u_h$ 在 $x_i$ 点处的取值

边界条件处理

Dirichlet边界条件

Dirichlet边界条件可知 $u(x)=0, x\in\partial\Omega$,应用在总刚度矩阵中,则有

1
2
3
4
5
6
7
8
9
10
输入: 边界节点布尔矩阵 PBnd, 总刚度矩阵 K, 总载荷矩阵 F
for i = 0 -> N-1 do
if PBnd[i] = 1 then
for j = 0 -> N-1 do
K[i][j] <- 0
end for
K[i][i] <- 1
F[i] <- 0
end if
end for

应用边界条件后的总刚度矩阵记为 $K_0$,总载荷矩阵记为 $F_0$。

与一维情形下的情况不同,刚度矩阵不是三对角矩阵,所以不用保持三对角形式。

5. 求解有限元方程组

求解线性方程组 $K_0 U=F_0$,具体方法此处不再详述。

1
2
from scipy.linalg import solve
u_lst = solve(a_mat, f_lst)

6. 误差检验

若已知方程真解 $u(x,y)$,且计算得到数值解 $u_h(x,y)$,则可计算误差

定义 $L^2$ 空间上的误差范数

若计算得

则说明数值计算结果是有效的。

7. 计算结果

计算数值解

易知选取的边值问题真解为 $u(x, y)=\sin(\pi x)\sin(\pi y)$

n=2

取 $n=2$,剖分得到 $N=9$ 个点,$M=8$ 个单元,计算得到 $E_2=0.3336632375743195$

n=4

取 $n=4$,剖分得到 $N=25$ 个点,$M=32$ 个单元,计算得到 $E_4=0.10742836489850027$

n=8

取 $n=8$,剖分得到 $N=81$ 个点,$M=128$ 个单元,计算得到 $E_8=0.028863378482441915$

n=16

取 $n=16$,剖分得到 $N=289$ 个点,$M=512$ 个单元,计算得到 $E_{16}=0.007355120927260682$

Tips: 具体的计算结果在solution.txt中展现。

计算误差阶数

根据 $L^2$ 误差,计算误差阶数如下:

n $L^2$ Error $P$ order
2 0.3336632375743195
4 0.10742836489850027 3.1059137676494366
8 0.028863378482441915 3.7219608565175686
16 0.007355120927260682 3.9242561431538148

可见 $P=\frac{E_N}{E_{2N}}$ 随着 $n$ 的增大,逐渐趋近于 4,所以该数值方法是有效的。

PS: 文中只简单的展示了部分核心代码,具体的实现请转到 Github: fempy,代码实现较为仓促,如有问题,欢迎批评指正~

附:三角单元积分运算

积分是限制有限元计算速度的一大关键因素。

对于一维单元上的函数 $f(x)$,计算 $\iint f(x)\mathrm{d}x$,二维三角单元上的函数 $f(x,y)$,计算 $\iint f(x,y) \mathrm{d}x\mathrm{d}y$,一般的数值计算方式会极大的影响迭代速度,则需要寻找一种高效且高精度的数值方法来近似。

在二维三角剖分网格上,由于三角单元的形状不一定完全相同,在一般的三角单元上进行运算时,需要定义不同的运算形式。为方便起见,我们定义等参单元的概念,将一般的三角单元变换为等参三角单元,从而定义统一的运算方式。

这里引入面积坐标的概念,选取三角单元中任一点 $p=(x,y)$,连接点 $p$ 与三角单元的三个节点 $p_1=(x_1, y_1), p_2=(x_2, y_2), p_3=(x_3, y_3)$,则可将三角单元划分为三个部分,并定义:

则 $(\lambda_1, \lambda_2, \lambda_3)$ 为三角单元的面积坐标,且易知 $\lambda_1+\lambda_2+\lambda_3=1$

则可定义仿射矩阵 $L$,将全局(global)坐标系中的点 $(x, y, 1)$ 映射为面积坐标系中的点 $(\lambda_1, \lambda_2, \lambda_3)$,面积坐标是局部(local)坐标系的一种。当 $p$ 选取三角单元的三个节点时,可以得到:

假设

由于

则可知

从而对于三角单元 $K$ 中的任意一点 $p=(x,y)$ 都有:

展开为:

方法一:直接计算

由面积坐标的定义可知 $\varphi_1=\lambda_1, \varphi_2=\lambda_2, \varphi_3=\lambda_3$,且有

则可知 $\nabla \varphi_i$ 都是常数:

则可以直接进行计算。

单元载荷矩阵

1
2
3
4
5
6
7
from scipy.integrate import quad, dblquad
f_lst = np.zeros(3)
fun = lambda x, y: f(p[2, 0] + x * (p[0, 0] - p[2, 0]) + y * (p[1, 0] - p[2, 0]),
p[2, 1] + x * (p[0, 1] - p[2, 1]) + y * (p[1, 1] - p[2, 1]))
f_lst[0] = dblquad(lambda x, y: fun(x, y) * x * det_j, 0, 1, lambda x: 0, lambda x: 1 - x)
f_lst[1] = dblquad(lambda x, y: fun(x, y) * y * det_j, 0, 1, lambda x: 0, lambda x: 1 - x)
f_lst[2] = dblquad(lambda x, y: fun(x, y) * (1 - x - y) * det_j, 0, 1, lambda x: 0, lambda x: 1 - x)

单元刚度矩阵

1
2
3
4
5
6
7
8
9
10
11
12
13
14
a_lst = np.zeros(9)
sa = (p[0, 0] - p[2, 0]) * (p[1, 1] - p[2, 1]) - (p[1, 0] - p[2, 0]) * (p[0, 1] - p[2, 1])
a_lst[0] = ((p[1, 1] - p[2, 1]) ** 2 + (p[2, 0] - p[1, 0]) ** 2) / sa
a_lst[1] = ((p[2, 1] - p[0, 1]) * (p[1, 1] - p[2, 1]) +
(p[0, 0] - p[2, 0]) * (p[2, 0] - p[1, 0])) / sa
a_lst[2] = ((p[0, 1] - p[1, 1]) * (p[1, 1] - p[2, 1]) +
(p[1, 0] - p[0, 0]) * (p[2, 0] - p[1, 0])) / sa
a_lst[3] = a_lst[1]
a_lst[4] = ((p[2, 1] - p[0, 1]) ** 2 + (p[0, 0] - p[2, 0]) ** 2) / sa
a_lst[5] = ((p[0, 1] - p[1, 1]) * (p[2, 1] - p[0, 1]) +
(p[1, 0] - p[0, 0]) * (p[0, 0] - p[2, 0])) / sa
a_lst[6] = a_lst[2]
a_lst[7] = a_lst[5]
a_lst[8] = ((p[0, 1] - p[1, 1]) ** 2 + (p[1, 0] - p[0, 0]) ** 2) / sa

方法二:Gauss积分

根据Gauss求积公式,选择不同的Gauss求积点 $\epsilon_i$,近似计算单元区域上的积分:

将定义在局部坐标系中的Gauss求积点 $(\lambda_1, \lambda_2)$,通过仿射变换,从自然坐标系变换到局部坐标系中:

坐标变换的Jacobi矩阵为:

则三角单元上的Gauss求积公式为:

三角等参单元上的Gauss点及Gauss权重数据,已经打包为 .npy 格式,由numpy库可直接读取。

Tips: 具体的Gauss积分形式,参考 basic.pyfem_2d.py,实现方式较为简单。

坚持数学知识分享,您的支持将鼓励我继续创作
0%