复杂区域上的 Laplace 方程

以 Laplace 方程为例，我们一起来看一下如何使用 FEtch 系统快速地求解偏微分方程问题。

问题描述

假设在二维区域 $\Omega$ 上温度 $u(x,y)$ 满足 Laplace 方程 $\frac{\partial^{2} u}{\partial x^{2}}(x, y)+\frac{\partial^{2} u}{\partial y^{2}}(x, y)=0 \quad((x, y) \in \Omega)$ 区域 $\Omega$ 如下图所示，并给定以下边界条件

\[\begin{align} u(x, y)&=4 \quad\left((x, y) \in \Gamma_{6}\cup \Gamma_{7}\right) \\ \frac{\partial u}{\partial n}(x, y)&=x \quad\left((x, y) \in \Gamma_{2} \cup \Gamma_{4}\right) \\ \frac{\partial u}{\partial n}(x, y)&=y \quad\left((x, y) \in \Gamma_{5}\right)\\ \frac{\partial u}{\partial n}(x, y)&=\frac{x+y}{\sqrt{2}} \quad\left((x, y) \in \Gamma_{1} \cup \Gamma_{3}\right) \end{align}\]

其中 $\frac{\partial u}{\partial n}$ 表示在点 $u(x,y)$ 处沿区域边界法线方向的方向导数。该偏微分方程的解析解为 $u(x,y)=xy+4$。试用有限元法求解该问题。

imag

脚本文件

这是一个比较简单的单物理场问题。用有限元语言表述，只需要填写 5 种类型的文件，即 PDE、FBC、NFE、GCN 和 MDI 文件。

PDE 文件 (laplace.pde)，用于描述微分方程的弱形式的体积分项。在有限元理论框架中，对应于单元矩阵的计算，是整个程序中最核心的部分。

defi
disp u
coor x y
shap %1 %2
gaus %3
mate ek ef 1.0 0.0

stif
dist = +[u/x;u/x]*ek+[u/y;u/y]*ek

load = +[u]*ef

end

FBC 文件 (laplace.fbc)，用于描述微分方程的弱形式的边界积分项，和 PDE 文件一样用于单元矩阵的计算。

disp u
coor x
shap %1 %2
gaus %3
coef gx gy
mate ex ey 0 0
 
stif
null

load=+[u]*(ex*gx+ey*gy)
 
end

NFE 文件 (lapa.nfe)，与 PDE 文件相配合，用于给出待求解的矩阵格式的线性方程组（对应于整体矩阵组装），以及对计算结果作进一步处理。

defi
stif s
load f

equation
matrix = [s]
forc = [f]

solution u
gidpost("u",1.0) u(1)

end

GCN 文件 (lap.gcn)，用于定义求解器类型。同时与 NFE 文件相配合，一起组织计算流程。

defi
a ell

#startsin a
#solvsin a

MDI 文件 (lap.mdi)，用于组织和管理物理场的信息。

2dxy
#a 0 1 u
pde laplace t3g2
fbc laplace l2g2
#

实际操作中，这里的 NFE 文件可以调用系统的库文件自动生成，因此用户只需填写 PDE、FBC、GCN 和 MDI 四个文件即可，统计下来，总行数不超过 30 行！

我们通过 FEtch 客户端提交脚本文件后，会自动生成并下载有限元程序文件（lap.exe 和 run.bat）、材料参数文件（lap.mat）和 GiD 软件的前处理模板文件（lap.cnd 和 lap.bas）。编程任务至此已经顺利完成。

前处理

首先编辑 lap.mat 文件，定义边界条件的参数。

1 3 aet3g2
num ek eq
1 1d0 0d0
#
3 3 all2g2
num ex ey
1 1d0 0d0
2 0d0 1d0
3 0.70710678d0 0.70710678d0
#

前处理和后处理都是通过第三方软件 GiD 来完成的。

我们首先进入 GiD 软件界面进行前处理。依次进行几何建模、施加边界条件、赋材料值和网格剖分，最终导出计算所需的数据输入文件（lap.dat)。

边界条件的施加情况如下

第一类边界条件

imag

第二类边界条件

imag

使用的网格如下

imag

计算结果

经过计算，生成了后处理文件（lap.flavia.res）。我们再次进入 GiD 软件界面进行后处理，查看云图效果如下。可以看出计算结果非常接近于解析解。

至此，Laplace 方程的求解已顺利完成。

参考文献

[1] 戴嘉尊, 邱建贤. 微分方程数值解法.第2版[M]. 东南大学出版社, 2012.