FEniCS编程：从基础到复杂问题求解

### FEniCS编程：从基础到复杂问题求解 #### 1. 矩阵与双线性形式 在FEniCS编程中，矩阵$M$常被称为“质量矩阵”，而$K$则被称为“刚度矩阵”。与之相关的双线性形式在FEniCS程序的组装过程中非常重要，具体如下： - $a_K(u, v) = \int_{\Omega} \nabla u \cdot \nabla v dx$ - $a_M(u, v) = \int_{\Omega} uv dx$ 每个时间层的线性系统$AU_k = b$，可以先计算$M$和$K$，然后在$t = 0$时形成$A = M + dtK$，而$b$在每个时间层计算为$b = MU_{k - 1} + dtMF_k$。 为了避免在每个时间层进行组装，需要对之前的`d1_d2D.py`程序进行以下修改： 1. 定义单独的形式$a_M$和$a_K$ 2. 将$a_M$组装到$M$，将$a_K$组装到$K$ 3. 计算$A = M + dtK$ 4. 将$f$定义为表达式 5. 将$f$的公式插值到有限元函数$F_k$ 6. 计算$b = MU_{k - 1} + dtMF_k$ 相关代码如下： ```python # 1. a_K = inner(nabla_grad(u), nabla_grad(v))*dx a_M = u*v*dx # 2. and 3. M = assemble(a_M) K = assemble(a_K) A = M + dt*K # 4. f = Expression("beta - 2 - 2*alpha", beta=beta, alpha=alpha) # 5. and 6. while t <= T: f_k = interpolate(f, V) F_k = fk.vector() b = M*u_1.vector() + dt*M*F_k ``` 完整程序在`d2_d2D.py`文件中。 #### 2. 物理实例：地面温度受昼夜变化的影响 考虑一个$d$维的盒形区域$\Omega$，在区域顶部$x_{d - 1} = 0$处有一个振荡温度$T_0(t) = T_R + T_A \sin(\omega t)$，其中$T_R$是参考温度，$T_A$是表面温度变化的振幅，$\omega$是温度振荡的频率。在其他边界，假设温度的法向导数为零。初始时，整个区域的温度为$T_R$。土壤的热导率$\kappa$可能随空间变化。 该问题的初边值问题如下： - $\varrho c\frac{\partial T}{\partial t} = \nabla \cdot (\kappa \nabla T)$ 在$\Omega \times (0, t_{stop}]$内 - $T = T_0(t)$ 在$\Gamma_0$上 - $\frac{\partial T}{\partial n} = 0$ 在$\partial \Omega \setminus \Gamma_0$上 - $T = T_R$ 在$t = 0$时 其中，$\varrho$是土壤密度，$c$是热容量，$\kappa$是热导率，$\Gamma_0$是表面边界$x_{d - 1} = 0$。 使用$\theta$ - 格式进行时间离散，演化方程$\frac{\partial P}{\partial t} = Q(t)$离散为$\frac{P_k - P_{k - 1}}{dt} = \theta Q_k + (1 - \theta) Q_{k - 1}$，应用到我们的PDE中得到： $\varrho c\frac{T_k - T_{k - 1}}{dt} = \theta \nabla \cdot (\kappa \nabla T_k) + (1 - \theta) \nabla \cdot (\kappa \nabla T_{k - 1})$ 将这个时间离散的PDE转化为弱形式，得到： - $a(T, v) = \int_{\Omega} (\varrho c T v + \theta dt \kappa \nabla T \cdot \nabla v) dx$ - $L(v) = \int_{\Omega} (\varrho c T_{k - 1} v - (1 - \theta) dt \kappa \nabla T_{k - 1} \cdot \nabla v) dx$ 由于Neumann边界条件，边界积分消失。 以下是创建网格和函数空间的代码： ```python degree = int(sys.argv[1]) D = float(sys.argv[2]) W = D/2.0 divisions = [int(arg) for arg in sys.argv[3:]] d = len(divisions) # no of space dimensions if d == 1: mesh = Interval(divisions[0], -D, 0) elif d == 2: mesh = Rectangle(-W/2, -D, W/2, 0, divisions[0], divisions[1]) elif d == 3: mesh = Box(-W/2, -W/2, -D, W/2, W/2, 0, divisions[0], divisions[1], divisions[2]) V = FunctionSpace(mesh, "Lagrange", degree) ``` 设置上边界的Dirichlet条件的代码如下： ```python T_R = 0; T_A = 1.0; omega = 2*pi T_0 = Expression("T_R + T_A*sin(omega*t)", T_R=T_R, T_A=T_A, omega=omega, t=0.0) def surface(x, on_boundary): return on_boundary and abs(x[d-1]) < 1E-14 bc = DirichletBC(V, T_0, surface) ``` $\kappa$函数在特定矩形区域内为常数$\kappa_1$，在其他区域为常数$\kappa_0$。可以通过以下两种方式定义： 方式一：使用`Expression`子类 ```python class Kappa(Function): def eval(self, value, x): """x: spatial point, value[0]: function value.""" d = len(x) # no of space dimensions material = 0 # 0: outside, 1: inside if d == 1: if -D/2. < x[d-1] < -D/2. + D/4.: material = 1 elif d == 2: if -D/2. < x[d-1] < -D/2. + D/4. and \ -W/4. < x[0] < W/4.: material = 1 elif d == 3: if -D/2. < x[d-1] < -D/2. + D/4. and \ -W/4. < x[0] < W/4. and -W/4. < x[1] < W/4.: material = 1 value[0] = kappa_0 if material == 0 else kappa_1 ``` 方式二：使用字符串表达式 ```python kappa_str = {} kappa_str[1] = "x[0] > -D/2 && x[0] < -D/2 + D/4 ? kappa_1 : kappa_0" kappa_str[2] = "x[0] > -W/4 && x[0] < W/4 "\ "&& x[1] > -D/2 && x[1] < -D/2 + D/4 ? "\ "kappa_1 : kappa_0" kappa_str[3] = "x[0] > -W/4 && x[0] < W/4 "\ "x[1] > -W/4 && x[1] < W/4 "\ ```