深度内容 | 北太天元PDE工具箱的有限元方法及算例详解

一、前言

偏微分方程（Partial Differential Equation, PDE）是数学中的一个重要分支，主要用于描述物理量随时间和空间变化的规律。偏微分方程的应用非常广泛，涵盖了物理学、工程学、生物学和化学等多个领域。例如，在物理学中，它们用于描述热传导、流体力学和电磁场理论等现象；在工程学中，用于描述结构力学、振动分析和传热等问题；在生物学和化学中，用于研究种群动态、化学反应和扩散过程等。此外，偏微分方程的应用也已扩展至经济学、金融预测、图像处理等现代领域。

大多数偏微分方程的解析解难以求得，数值方法成为求解偏微分方程的常用途径。数值解法旨在通过离散化的方式将偏微分方程转化为有限维的代数方程，从而可以在计算机上进行求解。常用的数值解法包括有限差分方法（Finite Difference Method, FDM）、有限元方法（Finite Element Method, FEM）和有限体积方法（Finite Volume Method, FVM）等。

北太天元PDE工具箱提供了偏微分方程数值解的有限差分法和有限元方法两类。本文主要介绍有限元方法及应用算例。

二、有限元方法基本原理

有限元方法的基本原理可以概括为将复杂问题分解为简单的小单元，通过局部近似来构建全局解u。首先将求解域划分为有限个简单的子区域，称为单元（elements）。这些单元可以是线段、三角形、四边形，四面体、六面体等。有限元法的关键是将问题从强形式转化为弱形式（积分形式），使得解的空间从连续可微函数扩展到更大的函数空间（通常是 Sobolev 空间）。在弱形式的基础上，结合单元的离散化和基函数，将问题的全局解u转化为一个矩阵形式的代数方程组$\mathbf{AU}=\mathbf{b}$。利用数值代数方法求解全局代数方程组得到节点上的近似解$\mathbf{U}=[u_1,u_2,\cdots ,u_N{]}^T$，最终利用基函数在整个区域上得到有限元逼近解：
$$u\approx u_h=\sum _{i=1}^N{u}_i{\phi }_i$$

1、弱形式

在有限元方法中，方程的弱形式化是一个核心原理。目的是将问题的解从要求在每个点上满足微分方程的强条件，转换为在一定函数空间内满足积分形式的条件。这种转化不仅可以简化问题，还可以容忍某些解的不连续性或低光滑性。

假设有如下问题：

$$\nabla \cdot (-k\nabla u)=g(u,x)\text{in}\Omega ，$$

该方程两端都乘以一个测试函数$\phi \in V_h$，并在$\Omega$上积分得

$${\int }_{\Omega }k\nabla \cdot (-k\nabla u)\phi dV={\int }_{\Omega }g\phi dV,\forall \phi \in V_h$$
然后通过分部积分可得：

$${\int }_{\Omega }k\nabla u\cdot \nabla \phi dV+{\int }_{\partial \Omega }(-k\nabla u)\cdot \mathbf{n}\phi dS={\int }_{\Omega }g\phi dV,\forall \phi \in V_h$$
该方程放宽了原方程的要求。原方程要求每一个点上都必须被明确定义，而弱形式公式只要求积分相等，可以允许解u不连续或不具备足够的光滑性。

2、伽辽金（Galerkin）法

伽辽金法用于求解微分方程的弱形式。它是有限元方法的基础，具体思想是在弱形式中采用一种特定的离散化策略：选择一组基函数 { φ 1 , φ 2 , ⋯   , φ n } \{\varphi_1,\varphi_2,\cdots,\varphi_n\} {φ1​,φ2​,⋯,φn​}，这些函数通常是多项式或其他合适的基函数，并且它们在问题的定义域上有某些特性（如满足边界条件）。这些函数组成了一个有限元逼近空间$V_h$，对于任意的函数$u_h\in V_h$可以表示为：
$$u_h=\sum _i{u}_i{\phi }_i(x)$$

代入到弱形式，可得：

$$\sum _i^n{u}_i{\int }_{\Omega }k\nabla {\phi }_i\cdot \nabla {\phi }_{j}dV+\sum _i^n{\int }_{\partial \Omega }(-k{u}_i\nabla {\phi }_i)\cdot \mathbf{n}{\phi }_{j}dS={\int }_{\Omega }g\left(\sum _i^n{u}_i\nabla {\phi }_i\right){\phi }_{j}dV,{\phi }_j\in V_h$$
然后就得到了一个线性方程组：

$$\mathbf{AU}=\mathbf{b},\mathbf{U}=[u_1,u_2,\cdots ,u_N{]}^T$$

$\mathbf{A}$是一个$n\times n$的矩阵，由上面的积分表达式决定。求解这些方程组，就给出了偏微分方程的数值近似解。

三、已实现的有限元方法求解PDE类型

北太天元PDE工具箱已实现的有限元方法求解PDE类型包括：

• 椭圆方程(组)

• 抛物方程(组)

• 双曲方程(组)

• 特征方程(组)

• 结构力学模型

• 静电学模型

• 静磁学模型

• 交流电力电磁学模型

• 导电介质直流模型

• 热传导模型

• 扩散模型

• 相场模型

• …

四、有限元方法求解PDE算例

下面从方程模块、方程组模块、应用模块对工具箱展开介绍。

1、方程模块

• 椭圆方程求解

考虑以下形式的椭圆型偏微分方程

$$-\nabla \cdot (c\nabla u)+au=f\text{, 在 }\Omega$$
该方程用于描述稳态的物理现象，定义在有界二维区域$\Omega$上，$\Omega$通常表示物理对象的一个截面或空间域。$u=u(x,y)$是在区域$\Omega$上的未知函数，$c=c(x,y)$、$a=a(x,y)$和$f=f(x,y)$是在区域$\Omega$上的已知函数，它们可以是实数或复数值函数。

考虑区域$\Omega$为正方形区域$[-1,1]\times [-1,1]$，系数为：
$$c=e^{x+y};a=0;f=\sin (\pi x)\cos (\pi y)$$

在边界$\partial \Omega$上，设定狄利克雷边界条件$u=0$。

网格模型和求解结果如下。

• 抛物方程求解

抛物型偏微分方程一般形式为：
$$d\frac{\partial u}{\partial t}-\nabla \cdot (c\nabla u)+au=f,\text{ 在 }\Omega$$
其中，$u=u(x,y,t)$表示物理量随时间$t$和空间位置$(x,y)$的变化，d是时间上的扩散系数，$c、a、f$分别为空间上的扩散系数、反应系数和源项。$f$可以依赖于时间和位置，$c、a$只依赖于位置。

考虑区域$\Omega$为正方形区域 $[-1,1]\times [-1,1]$。

对于方程取以下系数和初边值边界条件：
$$d=1;c=0.01;a=0;f=\sin (\pi x)\cos (\pi y)$$
在时间$t=0$时，初始条件设为
$$u(x,y,0)=15e^{-\frac{(x-1{)}^2+(y-1{)}^2}{0.02}}$$ 在边界$\partial \Omega$上，设定狄利克雷边界条件$u=0$。

求解后可视化如下：

• 双曲方程求解

双曲型偏微分方程一般形式为：

$$d\frac{{\partial }^{2}u}{\partial t^2}-\nabla \cdot (c\nabla u)+au=f,\text{ 在 }\Omega$$
其中，$u=u(x,y,t)$表示物理量随时间$t$和空间位置$(x,y)$的变化，$t$是时间上的扩散系数，$c、a、f$分别为空间上的扩散系数、反应系数和源项。$f$可以依赖于时间和位置，$c、a$只依赖于位置。

考虑区域$\Omega$为正方形区域$[-1,1]\times [-1,1]$。对于方程取以下系数和初边值边界条件：
$$d=1;c=0.01;a=0;f(x,y,t)=e^{-\frac{(x-1{)}^2+(y-1{)}^2}{2{\sigma }^2}}\ast \sin (\omega t)$$
初始条件设为：
$$u(x,y,0)=e^{-\frac{(x-1{)}^2+(y-1{)}^2}{2{\sigma }^2}};\frac{\partial u}{\partial t}(x,y,0)=0$$
在边界$\partial \Omega$上，设定狄利克雷边界条件 $u=0$。

求解后可视化如下：

2、方程组模块

• 椭圆方程组求解

一般椭圆方程组表达式为：

$$-\nabla \cdot (c\nabla u)+au=f$$
其中，

$c=\left[\begin{array}{cc}{c}_{11}& c_{12}\\ c_{21}& c_{22}\end{array}\right],$$a=\left[\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right],$$f=\left[\begin{array}{c}{f}_1\\ f_2\end{array}\right]$

$c_{11},c_{12},c_{21},c_{22}$：控制二维系统中的扩散率，$a_{11},a_{12},a_{21},a_{22}$： 控制二维系统中的对流作用，$f_1,f_2$：体力源项或外部作用源。

• 抛物方程组求解

一般抛物方程组表达式为：
$$d\frac{\partial u}{\partial t}-\nabla \cdot (c\nabla u)+au=f$$其中，

$d=\left[\begin{array}{cc}{d}_{11}& d_{12}\\ d_{21}& d_{22}\end{array}\right],$ $c=\left[\begin{array}{cc}{c}_{11}& c_{12}\\ c_{21}& c_{22}\end{array}\right],$$a=\left[\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right],$$f=\left[\begin{array}{c}{f}_1\\ f_2\end{array}\right]$

$c_{11},c_{12},c_{21},c_{22}$:控制二维系统中的扩散率，$a_{11},a_{12},a_{21},a_{22}$：控制二维系统中的对流作用，$f_1,f_2$：体力源项或外部作用源$d_{11},d_{12},d_{21},d_{22}$： 控制时间演化的质量系数。

• 双曲方程组求解

一般双曲方程组表达式为：

$$d\frac{{\partial }^{2}u}{\partial t^2}-\nabla \cdot (c\nabla u)+au=f$$
其中，

$d=\left[\begin{array}{cc}{d}_{11}& d_{12}\\ d_{21}& d_{22}\end{array}\right],$$c=\left[\begin{array}{cc}{c}_{11}& c_{12}\\ c_{21}& c_{22}\end{array}\right],$$a=\left[\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right],$$f=\left[\begin{array}{c}{f}_1\\ f_2\end{array}\right]$

$c_{11},c_{12},c_{21},c_{22}$: 控制二维系统中的扩散率，$a_{11},a_{12},a_{21},a_{22}$：控制二维系统中的对流作用，$f_1,f_2$：体力源项或外部作用源，$d_{11},d_{12},d_{21},d_{22}$：控制时间演化的质量系数。

• 特征方程组求解

一般特征值方程组表达式为：
$$-\nabla \cdot (c\nabla u)+au=\lambda du$$其中，

$d=\left[\begin{array}{cc}{d}_{11}& d_{12}\\ d_{21}& d_{22}\end{array}\right],$$c=\left[\begin{array}{cc}{c}_{11}& c_{12}\\ c_{21}& c_{22}\end{array}\right],$ $a=\left[\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right]$

$c_{11},c_{12},c_{21},c_{22}$: 控制二维系统中的控制扩散和模态振动的传播，$a_{11},a_{12},a_{21},a_{22}$：控制二维系统中的对流作用，$d_{11},d_{12},d_{21},d_{22}$：控制振动加速度的相关系数。

3、应用模块

• 结构力学模型（应力）

结构力学模型用于描述材料在力作用下的应力和应变行为。平面应力问题位移方程为：
$$-\nabla \cdot (c\otimes \nabla u)=K$$其中， $u$是二维向量，$c$是秩为 4 的张量，$K$为体积力，张量$c$可以写成 4 个 2×2 的矩阵：

$c_{11}=\left[\begin{array}{cc}2G+\mu & 0\\ 0& G\end{array}\right]$, $c_{12}=\left[\begin{array}{ll}0& \mu \\ G& 0\end{array}\right]$,
$c_{21}=\left[\begin{array}{ll}0& \mathrm{G}\\ \mu & 0\end{array}\right]$，$c_{22}=\left[\begin{array}{cc}\mathrm{G}& 0\\ 0& 2G+\mu \end{array}\right]$

剪切模量 $G$定义为 $G=\frac{E}{2(1+v)}$，$\mu$定义为$\mu =\frac{2Gv}{1-v}$， 外力$F=\left[\begin{array}{c}{K}_x\\ K_y\end{array}\right]$，

方程涉及杨氏模量 $E$、泊松比 $\nu$、以及作用在材料上的体积力$K$。

• 静电学模型

静电学模型用于描述静止电荷产生的电场分布。该模型的方程形式为：
$$-\nabla \cdot (\epsilon \nabla V)=\rho$$其中， $V$表示电势，$\epsilon$为介电系数，$\rho$为空间电荷密度。

电场$E$通过$V$的梯度给出，即 $E=-\nabla V$。

该方程用于分析静电场，例如电容器中的电场分布、绝缘体中的电荷行为等。

• 交流电力电磁学模型

交流电力电磁学模型用于描述交流电场和磁场的相互作用。该模型的方程形式为：

$$-\nabla \cdot ((1/\mu )\nabla E)+\left(i\omega \sigma -{\omega }^2\epsilon \right)E=0$$
其中 ，$E$表示电场强度，$\omega$是角频率，$\mu$为磁导率，$\sigma$为电导率，$\epsilon$为介电系数。

该方程描述了电场的传播、衰减及其与介质的相互作用，广泛用于天线设计、电力传输等领域。

• 热传导模型

热传导问题的抛物型偏微分方程：
$$\rho C\frac{\partial T}{\partial t}-\nabla \cdot (k\nabla T)=Q+h(Text-T)$$方程描述了二维热传导现象，以及轴对称三维问题经简化后的热传导问题，

其中，$T$：温度，$\rho$：密度，$C$：比热，$k$：导热系数，$Q$：热源，$h$：传热系数，$T_{ext}$：环境温度。

稳态情形下的方程简化为：

$$-\nabla \cdot (k\nabla T)=Q+h(Text-T)$$边界条件可以为Dirichlet条件、Neumann条件或Robin条件。

• 相场方程

Cahn-Hilliard 方程描述二相系统中的相分离过程，适用于合金、流体混合物等具有相分离特性的系统。方程如下：
$$\frac{\partial c}{\partial t}=\nabla \cdot \left(M\nabla \left(\frac{\partial f(c)}{\partial c}-\kappa {\nabla }^{2}c\right)\right)$$其中，$c$表示浓度或组分分布，$M$为系数，$f(c)$是自由能密度函数，通常采用双井势函数形式，以确保在系统中形成两相区域。 $\kappa$是一个常数，控制界面的宽度和平滑性。Cahn-Hilliard 方程广泛应用于材料科学中的二元合金相分离、流体混合物的扩散过程、以及相变和传热模拟。它为研究物质在特定条件下的分离行为提供了重要工具。该方程可以帮助我们深入理解不同物质在不同条件下的相互作用、扩散行为，以及它们如何在空间中分离、重组。

另外，北太天元工具箱还开发了其他相场方程有限元方法，可用于描述两相界面的演化，如枝晶生长、材料裂纹拓展、形核和界面演化等现象。

五、PDE工具箱特色

北太天元PDE工具箱对外部网格具有很强的通用性，可读入外部网格用于计算。对于导入的外部网格，可自动识别求解区域网格的边界，并直接应用用户指定的边界条件进行计算。同时，PDE工具箱具有后处理技术。

六、参考文献

1. 王烈衡, 许学军. 有限元方法的数学基础，科学出版社，2004.

2. 杜其奎, 陈金如. 有限元方法的数学理论，科学出版社，2012.

3. Vidar Thomée. Galerkin Finite Element Methods for Parabolic Problems. Springer-Verlag Berlin Heidelberg, 2006.

4. S.C. Brenner, L.R. Scott. The Mathematical Theory of Finite Element Methods. Springer New York, 2008.