有限元算例-线弹性方程求解#
本节将引导您使用 FEALPy 完成二维线弹性方程的有限元求解，现定义如下数学模型：
\[\begin{split}
\begin{cases}
-\nabla \cdot \boldsymbol{\sigma}(\mathbf{u}) = \mathbf{f} & \text{在 } \Omega = [0,1]^2
\\
\mathbf{u} = \mathbf{g} & \text{在 } \partial \Omega
\end{cases}
\end{split}\]
其中应力张量 \(\boldsymbol{\sigma}(\mathbf{u})\) 由广义胡克定律给出：
\[
\boldsymbol{\sigma}(\mathbf{u}) = 2\mu \boldsymbol{\varepsilon}(\mathbf{u}) + \lambda \mathrm{tr}(\boldsymbol{\varepsilon}(\mathbf{u})) \mathbf{I},
\]
也可以写成弹性矩阵的形式：
\[
\boldsymbol{\sigma}(\boldsymbol{u}) = \boldsymbol{D}\boldsymbol{\varepsilon}(\boldsymbol{u})
\]
应变张量定义为：
\[
\boldsymbol{\varepsilon}(\mathbf{u}) = \frac{1}{2} \left( \nabla \mathbf{u} + (\nabla \mathbf{u})^T \right).
\]
取材料参数为：杨氏模量 \(E = 1.0\)，泊松比 \(\nu = 0.3\)，并定义对应的拉梅常数：
\[
\mu = \frac{E}{2(1+\nu)}, \quad \lambda = \frac{E \nu}{(1+\nu)(1-2\nu)}.
\]
其中定义真解为：
\[\begin{split}
\mathbf{u}(x, y) = \begin{bmatrix}
\sin(\pi x) \sin(\pi y) \\
0
\end{bmatrix}
\end{split}\]
对应的源项为：
\[\begin{split}
\mathbf{f}(x, y) = \begin{bmatrix}
\displaystyle \frac{22.5\pi^2}{13} \sin(\pi x) \sin(\pi y) \\
\displaystyle -\frac{12.5\pi^2}{13} \cos(\pi x) \cos(\pi y)
\end{bmatrix}
\end{split}\]
变分形式为：
\[
\int_\Omega \boldsymbol{\sigma}(\boldsymbol{u_h}) :
\boldsymbol{\varepsilon}(\boldsymbol{v_h}) ~ \mathrm{d}\boldsymbol{x}
= \int_\Omega \boldsymbol{f}\cdot \boldsymbol{v_h} ~ \mathrm{d}\boldsymbol{x}
+\int_{\partial\Omega} \boldsymbol{g}\cdot\boldsymbol{v_h} ~ \mathrm{d}\boldsymbol{s},\quad\forall\boldsymbol{v_h}\in V_h
\]
1. 定义 PDE 模型#
本节首先定义 PDE 的数学模型及其相关参数、真解和边界条件。
import math
from fealpy.backend import backend_manager as bm
from fealpy.decorator import cartesian
# 定义域
def domain():
return [0, 1, 0, 1]
# 杨氏模量
def E():
return 1.0
# 泊松比
def nu():
return 0.3
# 剪切模量
def mu():
return E() / (2 * (1 + nu()))
# 拉梅常数
def lam():
return E() * nu() / ((1 + nu()) * (1 - 2 * nu()))
# 真解
@cartesian
def solution(p):
x, y = p[..., 0], p[..., 1]
u1 = bm.sin(math.pi * x) * bm.sin(math.pi * y)
u2 = bm.zeros_like(u1)
return bm.stack([u1, u2], axis=-1)
# 源项
@cartesian
def source(p):
x, y = p[..., 0], p[..., 1]
π2 = math.pi ** 2
f1 = (22.5 * π2 / 13.0) * bm.sin(math.pi * x) * bm.sin(math.pi * y)
f2 = -(12.5 * π2 / 13.0) * bm.cos(math.pi * x) * bm.cos(math.pi * y)
return bm.stack([f1, f2], axis=-1)
# dirichlet 边界条件
@cartesian
def dirichlet(p):
return solution(p)
2. 参数配置与初始化#
本节将进行后端设置、日志工具导入、网格初始化等准备工作。
2.1 设置后端#
选择数值计算后端（如 numpy、cupy 等）及计算设备。
from fealpy.backend import backend_manager as bm
backend = 'numpy'
device = 'cpu'
bm.set_backend(backend)
2.2 导入日志工具#
用于记录各阶段的耗时与调试信息。
from fealpy.utils import timer
from fealpy import logger
logger.setLevel('WARNING')
tmr = timer()
next(tmr)
2.3 设置初始网格和加密次数#
初始化计算域的网格，并设置网格加密的迭代次数。
from fealpy.mesh import TriangleMesh
mesh = TriangleMesh.from_box(domain(), nx=4, ny=4)
maxit = 4
2.4 定义误差存储矩阵#
用于记录每次网格加密下的数值误差。
errorMatrix = bm.zeros((1, maxit), dtype=bm.float64)
3. 有限元求解#
本节将依次完成有限元空间构建、矩阵组装、边界处理、线性系统求解、误差计算与网格加密。
流程包含：
构建 Lagrange 有限元空间
组装刚度矩阵 \(A\) 和载荷向量 \(F\)
处理 Dirichlet 边界条件
求解线性系统 \(A u_h = F\)
计算 \(L^2\) 误差 \(\|u - u_h\|_{L^2(\Omega)}\)
网格均匀加密
from fealpy.functionspace import functionspace
from fealpy.material import LinearElasticMaterial
from fealpy.fem import BilinearForm, LinearElasticityIntegrator
from fealpy.fem import LinearForm, ScalarSourceIntegrator
from fealpy.fem import DirichletBC
from fealpy.solver import cg
for i in range(maxit):
GD = mesh.geo_dimension()
space = functionspace(mesh, ('Lagrange', 1), shape=(GD, -1))
tmr.send(f'第{i}次空间时间')
uh = space.function()
LEM = LinearElasticMaterial(name='E1nu025',
lame_lambda=lam(), shear_modulus=mu(),
hypo='plane_strain', device=device)
LEI = LinearElasticityIntegrator(material=LEM, q=3)
bform = BilinearForm(space)
bform.add_integrator(LEI)
A = bform.assembly()
lform = LinearForm(space)
lform.add_integrator(ScalarSourceIntegrator(source, q=3))
F = lform.assembly()
tmr.send(f'第{i}次矩阵组装时间')
gdof = space.number_of_global_dofs()
A, F = DirichletBC(space, gd = dirichlet).apply(A, F)
tmr.send(f'第{i}次边界处理时间')
uh[:] = cg(A, F)
tmr.send(f'第{i}次求解器时间')
errorMatrix[0, i] = mesh.error(solution, uh)
if i < maxit-1:
mesh.uniform_refine(n=1)
tmr.send(f'第{i}次误差计算及网格加密时间')
4. 误差分析与收敛阶计算#
本节将输出每次网格加密下的误差，并计算收敛阶。
next(tmr)
print("最终误差",errorMatrix)
print("order : ", bm.log2(errorMatrix[0,:-1]/errorMatrix[0,1:]))
Timer received None and paused.
=================================================
ID Time Proportion(%) Label
-------------------------------------------------
1 574.838 [ms] 91.366 第0次空间时间
2 11.292 [ms] 1.795 第0次矩阵组装时间
3 2.021 [ms] 0.321 第0次边界处理时间
4 0.000 [us] 0.000 第0次求解器时间
5 0.000 [us] 0.000 第0次误差计算及网格加密时间
6 0.000 [us] 0.000 第1次空间时间
7 3.089 [ms] 0.491 第1次矩阵组装时间
8 0.000 [us] 0.000 第1次边界处理时间
9 994.205 [us] 0.158 第1次求解器时间
10 0.000 [us] 0.000 第1次误差计算及网格加密时间
11 0.000 [us] 0.000 第2次空间时间
12 4.987 [ms] 0.793 第2次矩阵组装时间
13 1.008 [ms] 0.160 第2次边界处理时间
14 998.259 [us] 0.159 第2次求解器时间
15 1.988 [ms] 0.316 第2次误差计算及网格加密时间
16 0.000 [us] 0.000 第3次空间时间
17 19.967 [ms] 3.174 第3次矩阵组装时间
18 980.377 [us] 0.156 第3次边界处理时间
19 5.002 [ms] 0.795 第3次求解器时间
20 1.997 [ms] 0.317 第3次误差计算及网格加密时间
=================================================
最终误差 [[0.0817831 0.02305386 0.00601597 0.00152303]]
order : [1.82679432 1.93813903 1.98185488]
5. 结果可视化#
在单元重心处计算真解和数值解，并进行可视化比较。
注：由于真解 \(y\) 分量为 0，故仅展示 \(x\) 分量的可视化结果。
from matplotlib import pyplot as plt
bc = bm.array([[1/3, 1/3, 1/3]], dtype=bm.float64)
ps = mesh.bc_to_point(bc)
u = solution(ps)
ux = u[..., 0]
uh = uh(bc)
uhx = uh[..., 0]
fig, axes = plt.subplots(1, 2)
mesh.add_plot(axes[0], cellcolor=ux, linewidths=0)
axes[0].set_title('真解', fontname='Microsoft YaHei')
mesh.add_plot(axes[1], cellcolor=uhx, linewidths=0)
axes[1].set_title('数值解', fontname='Microsoft YaHei')
plt.show()
