空间曲线上的面积最小曲面问题

关键词： 稳态 非线性 Laplace方程

非线性偏微分方程的求解是有限元方法的重要应用领域。本节以空间曲线上的面积最小曲面问题为例，讲解如何使用 FEtch 系统求解非线性椭圆型偏微分方程。

问题描述

给定一条空间曲线，以此空间曲线为边界，张在它上面的曲面有无数多个。那么，其中哪个曲面的面积是最小的呢？

此问题用数学语言描述为：在 $XY$ 平面区域 $\Omega$ 的边界 $\partial \Omega$ 上给定函数值 $\varphi(x, y)$，求张在这一条曲线上的曲面，使得曲面的面积最小。

即 \(A(u)=\operatorname{inf} \int_{\Omega} \sqrt{1+|\nabla u|^{2}} \mathrm{~d} \Omega=\operatorname{inf} \iint_{\Omega} \sqrt{1+\left(\frac{\partial u}{\partial x}\right)^{2}+\left(\frac{\partial u}{\partial y}\right)^{2}} \mathrm{~d} x \mathrm{~d} y\)

控制方程

微分形式

基于 Euler-Lagrange 方程，可将上述问题转化为求解非线性的 Laplace 方程 \(\frac{\partial}{\partial x}\left(k(u) \frac{\partial u}{\partial x}\right)+\frac{\partial}{\partial y}\left(k(u) \frac{\partial u}{\partial y}\right)=0\ \left(\text{in } \Omega\right) \tag{1}\label{eq1}\)

其中，系数 $k(u)$ 为 \(k(u)=1/\sqrt{1+\left(\frac{\partial u}{\partial x}\right)^{2}+\left(\frac{\partial u}{\partial y}\right)^{2}}\)

弱形式

运用变分原理，由微分方程 $\eqref{eq1}$ 两端乘以 $u$ 的变分，并在全域上积分，得 \(\int_{\Omega}\left(\frac{\partial}{\partial x}\left(k(u) \frac{\partial u}{\partial x}\right)+\frac{\partial}{\partial y}\left(k(u) \frac{\partial u}{\partial y}\right)\right) \delta u d \Omega=0\) 对于二阶导数的部分，利用分部积分将其中的一阶导数转化到 $u$ 的变分上，得 \(\int_{\Omega}\left(k(u) \frac{\partial u}{\partial x} \frac{\partial \delta u}{\partial x}+k(u) \frac{\partial u}{\partial y} \frac{\partial \delta u}{\partial y}\right) d \Omega=0 \tag{2}\label{eq2}\)

这就是最终的弱解形式。

离散化

式 $\eqref{eq2}$ 进行空间离散化后，可写成矩阵形式 \(S(U)U=0\tag{3}\label{eq3}\) 其中，$S(U)$ 为刚度矩阵，$U$​ 为待求解的未知向量。该式一个非线性方程组，无法直接求解，需要引入迭代过程。

最常用的处理方法是采用直接迭代法（也称 Picard 法），将式 $\eqref{eq3}$ 线性化，得 \(S(U^i) U^{i+1}=0 \tag{4}\label{eq4}\)

其中，上标 $i$ 表示迭代步数。上式需要反复地迭代计算，直至误差 $\Vert U^{i+1}-U^i\Vert$ 小于某一设定的阈值。直接迭代法可以通过使用松弛因子，有效地提高收敛速度。松弛迭代法已被写入 NFE 算法模板库，是 FEtch 系统默认的非线性迭代算法。

可以看出，迭代过程的引入是非线性问题较线性问题计算量明显增大的主要原因。因此，如何在保证计算稳定性的同时提高收敛速度，一直是非线性求解领域永恒的话题。

脚本编写

一个具体问题的有限元语言表述，通常需要包含 5 种类型的文件：PDE、FBC、MDI、GCN 和 NFE。由于方程 $\eqref{eq2}$ 不存在边界积分项，所以不需要编写 FBC 文件。下面依次给出其他 4 个文件的详细介绍。

PDE 文件

在这个例子中，只有一个物理场，需要给出温度场对应的 PDE 文件和 FBC 文件。根据有限元语言的语法规则，将方程 $\eqref{eq2}$ 的体积分项写入 PDE 文件如下。

这里将 shap 行和 gaus 行的参数设为待定，参数值 %1、%2 和 %3 将从 MDI 文件中取得。

MDI 文件

MDI 文件负责组织 PDE 和 FBC 文件，声明程序具体使用的坐标系统、单元类型和积分方法。这里文件命名为 nlap.mdi，具体内容如下。

MDI 文件声明的单元类型和积分方式将以参数的形式传递给 PDE 文件中的的 %1、%2 和 %3。如果后续想要使用其他的单元类型或者积分方法，只需修改 MDI 文件即可。

GCN 文件

GCN 文件给出计算流程信息，系统规定其文件名字要与 MDI 文件一致。作为线性椭圆方程，稳态热传导问题的计算流程非常简单，只需填写少量代码，具体内容如下。

NFE 文件

GCN 文件使用了 nell 算法，系统会基于模板文件自动生成以下的 nlapa.nfe 文件。如果用户有其他需要，可以在此基础上进一步地修改。

算例

问题描述

在本例中，取 $\Omega$ 为位于坐标原点的一个单位圆盘，边界上的 $\varphi(x, y)=x^2$，求此时面积最小的曲面。

前处理

填写材料参数文件

系统自动生成了 nlap.mat 文件，用于设置材料参数。参数初始值全部为 PDE 文件给出的默认值。可根据具体问题进行修改。对于本算例，文件内容如下。

赋单元材料号

赋迭代初始值

施加边界条件

其中：
u-I 为边界类型，-1 表示约束，1 表示自由
u-D 为相应的约束值或荷载值

网格剖分

根据所开发程序的计算要求，采用四边形线性单元进行网格剖分。共 1203 个单元，1292 个节点。

计算结果

可以看到，限制在单位圆区域，以 $x^2$ 为边界的面积最小的曲面是马鞍面。计算结果与理论解完全一致，充分证明了算法和程序的有效性。