矩形区域上的扩散方程

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

问题描述

函数 $u(x,y,t)$ 满足扩散方程： $\frac{\partial u}{\partial t}=a\left(\frac{\partial^{2} u}{\partial x^{2}}+\frac{\partial^{2} u}{\partial y^{2}}\right)\ \left(\text{in } \Omega\right)$

初始条件：

\[u(x,y,0) = u_0 \ \left(\text{in } \Omega\right)\]

边界条件：

\[u(x,y,t) = 0 \ \left(\text{on } \partial \Omega\right)\]

其中，$\Omega$ 是一个单位长度的正方形，$\partial \Omega$ 是 $\Omega$ 的边界。

该偏微分方程的解析解为 $\begin{aligned} u= & \frac{16 u_{0}}{\pi^{2}}\left\{\sum_{n=0}^{\infty} \frac{1}{2 n+1} \sin \left[\frac{(2 n+1) \pi x}{l_{1}}\right] \exp \left[-\frac{\pi^{2}(2 n+1)^{2} a t}{l_{1}{ }^{2}}\right]\right\} \\ & \times\left\{\sum_{n=0}^{\infty} \frac{1}{2 m+1} \sin \left[\frac{(2 m+1) \pi y}{l_{2}}\right] \exp \left[-\frac{\pi^{2}(2 m+1)^{2} a t}{l_{2}{ }^{2}}\right]\right\} \end{aligned}$

令 $a=1\times10^{-4}$，$u_0=10$，试用有限元法求解 $t=1000$ 内 $u(x,y,t) $ 的演化过程。

脚本文件

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

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

defi
disp u
coor x y
shap %1 %2
gaus %3
mass %1 ec
mate ec ea 1 1e-4

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

load = +[u]*0

end

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

defi
stif s
mass m
load f

equation
vect u
read(s,unod) u
matrix = [s]*dt+[m]
forc = [f]*dt+[m]*[u]

solution u
write(s,unod) u
gidpost("u",time) u(1)

end

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

defi
a parb

open(1,file="time.dat")
read(1,*) dt,tmax
close(1)

#startsin a
time = 0.0
do while (time <= tmax-1d-8)
  time = time + dt
  if (time > tmax) then
    dt = dt-(time-tmax)
    time = tmax
  endif
  #solvsin a
enddo