怎么用python进行有限元分析_用python实现简单的有限元方法（二）

华中师范大学 hahakity

第一部分介绍了加权残差法(Weighted Residual Method)并引出有限元算法中经常出现的伽辽金(Galerkin)法。简单说明了有限元与有限差分法的区别：有限差分是对微分算子做差分近似，有限元是对待求函数做基函数展开近似。这一节继续介绍伽辽金有限元算法。

这里特意强调伽辽金有限元算法，是因为还有另一种基于泛函和变分的 Ritz 有限元算法。跳过 Ritz 变分法，直接讲解伽辽金法入门比较简单。

学习目标Galerkin 有限元法

分片基底函数展开

二阶微分算子的降阶

一维泊松方程：已知电荷密度求电势

预备知识加权残差法 (Weighted Residual Method)

Galerkin 有限元法

假设要求解的微分方程有统一的格式，

其中

是微分算子，

是要求解的场，

是源。

将

做函数基底展开 (使用爱因斯坦求和规则，省略

求和符号)，

注意此时

为基底函数，人工设定，为已知量，

为待定系数。

选取基底函数

作为伽辽金法中的检验函数，则对第 i 个检验函数，加权残差为，

其中

是全域，

,

可以将

写成矩阵乘积形式，此处用

表示列向量，

是

构成的矩阵，

是

构成的列向量，

是

构成的列向量。

Galerkin 有限元法令残差

为0，通过解

，得到待定系数

。

分片基底函数展开

将全域离散化为 n 个子域，每个子域称作一个单元(“element”)。连续的函数

近似为有限个单元上的分片函数。如果单元很小，每个单元内

近似线性变化。

先考虑简单的一维问题，将区间分成 n 等份，每份长

。

对第 e 个单元区域

，线性插值函数为

其中

，

当

时，

,

当

时，

,

根据

的定义，得到第 j 个基底函数(注意它横跨两个单元)，

的函数图像就是下图中间那个大三角形，覆盖第 j 和第 j+1 个单元。它跟

在第 j 个单元有重叠，跟

在第 j+1 个单元有重叠。此时内积都是局部进行，不重叠区域积分为0。

根据

可知第 j 行只有

,

,

三个非零矩阵元，K 矩阵是稀疏的 3 对角阵。

二阶微分算子降阶

算

时，很快遇到第一个问题，线性插值函数一阶微分在每个单元内为常数，

二阶或更高阶导数等于0。如果微分算子里有二阶导数，可以用分部积分法降阶，

二阶微分算子降阶之后变为，

只与全域边界条件有关，对于大部分不在边界上的点，这项为0。刚好在边界上的那些点，

或它的导数已知，因此可以将这一项单独列出来，令

在新的定义下，方程变为，

多出来的那项

是分部积分项在边界处的贡献

构成的列向量。

此时矩阵

第 j 行的非零矩阵元为,

每个矩阵元的积分区间由前面的示意图给出。

注意此时

或者等于 1/h，或者等于 -1/h，积分很好算。

简单的泊松方程求解

我们已经有了足够的知识，从头构造一个简单的一维 Galerkin 有限元算法。这里先演示一维有限元算法。步骤总结如下，将全域离散化为有限个子域(又称单元)

选择基底函数，对每个子域插值

用 Galerkin 法改写微分方程，函数内积计算

的矩阵元

解

这里拿一个简单的一维泊松方程演示伽辽金有限元算法。

物理问题：根据 Maxwell 方程，电场的散度正比于电荷密度

，而电场本身又可写成标量势

的负梯度

, 其中

。将第二个方程代入第一个得到给定电荷分布下静电势满足的柏松方程，

，二维情况下方程为，

，

一维情况下

其中

是真空电极化常数。

选择求解区域

，边界条件

。

选择一个简单的电荷密度分布，

此时，泊松方程有简单的解析解：

, 我们先假装不知道。

计算

的矩阵元，

这个积分用手积挺容易出错，写一小段 python 代码，将积分区域划分为

与

。在 jupyter-notebook 里运行，

import sympy

xjm1, xj, xjp1 = sympy.symbols(['x_{j-1}', 'x_{j}', 'x_{j+1}'])

x, h, eps = sympy.symbols(['x', 'h', '\epsilon'])

rho = lambda x: sympy.sin(sympy.pi * x)

A = - sympy.integrate(rho(x) * (x - xjm1)/h, (x, xjm1, xj))

B = - sympy.integrate(rho(x) * (xjp1 - x)/h, (x, xj, xjp1))

sympy.simplify(A + B)

输出结果：

Jackson 电动力学书本中介绍了一种简化方法，假设

在每个单元里是常数，只对

积分，则积分结果为，

先不管

, 矩阵方程的解为，

。写一段简单的 python 代码，构造矩阵

和列向量

，使用 scipy.sparse.linalg.inv 函数求

的逆，进而得到结果。

from scipy.sparse import dia_matrix

from scipy.sparse.linalg import inv

from numpy import pi

class FEM:

def __init__(self, nodes, xmin=0, xmax=1):

self.nodes = nodes

x = np.linspace(xmin, xmax, nodes)

self.x = x

self.h = x[1] - x[0]

def Kmatrix(self):

n = self.nodes

m = 1/self.h * np.ones(n)

data = [m, -2*m, m]

offsets = [-1, 0, 1]

# 使用 scipy.sparse 稀疏矩阵库，构造 3 对角稀疏矩阵 K

K = dia_matrix((data, offsets), shape=(n, n)).tocsc()

return K

def bvec(self):

'''假设 rho(x) 在每个单元内值为常数，仅对 u_j(x) 做积分'''

return - np.sin(pi * self.x) * self.h

def solve(self):

K = self.Kmatrix()

b = self.bvec()

return inv(K) * b

def compare(self):

ground_truth = 1/(pi**2) * np.sin(pi * self.x)

fem_res = self.solve()

plt.plot(self.x, ground_truth, label="ground truth")

plt.plot(self.x, fem_res, label="finite element")

plt.title("number of nodes =%s"%self.nodes)

plt.xlabel("x")

plt.ylabel(r"$\phi(x)$")

plt.legend(loc='best')

plt.show()

如果用 1001 个节点，计算结果比较精确，101 或 11 个节点，计算结果误差较大。

fem = FEM(1001)

fem.compare()

二维或高维时的网格生成

上面的一维问题有限元与有限差分算法区别很小，将方程右边的

除到左边，则 K 矩阵变成了二阶导数的差分矩阵。但对于高维且有复杂边界的偏微分问题，有限元算法就显示出它的优势。这一篇变得太长，下一篇尝试构造通用的二维和三维有限元求解器，以及如何在二维圆形区域求解本文例子中的静电势。

为了简化任务，使用 gmsh 软件的 python 库生成有限元分析需要的二维和三维网格。安装方法如下，

pip install gmsh

pip install pygmsh

使用 pygmsh，可以很快生成一个二维的圆形区域网格，

import pygmsh

with pygmsh.geo.Geometry() as geom:

geom.add_circle([0.0, 0.0], 1.0, mesh_size=0.2)

mesh = geom.generate_mesh()

# 将网格文件保存为 vtk 格式

mesh.write("test.vtk")

安装 paraview，打开 test.vtk, 则可以看到如下二维圆形区域的网格图，

总结

介绍了伽辽金有限元算法，一维分片函数展开，二阶微分算子降阶，数值求解了简单的一维泊松方程，介绍了网格生成软件 gmsh。将二维和三维有限元算法的通用求解器推迟到下一节。

警告：本文的求解器非通用求解器，仅作原理展示使用。

参考文献：

【1】Tao Pang，An Introduction to computational physics (second edition)