学习中心

使用偏微分方程（PDE）接口建模：泊松方程和拉普拉斯方程

---

本文为“使用偏微分方程（PDE）建模”系列课程的第 1 部分。首先，我们将快速介绍如何使用 COMSOL Multiphysics^® 软件中的通用偏微分方程接口。您将了解如何使用基于有限元法（FEM）和边界元法（BEM）的泊松方程和拉普拉斯方程接口建模。COMSOL 软件中内置的数学接口，可用于模拟这些特定的偏微分方程（即 拉普拉斯方程 和 泊松方程 ）。但是，为了演示如何求解您自己的方程，我们将不使用这些预置接口。您还将看到关于如何求解地月系统中牛顿重力的示例，这在软件中并没有相应的内置选项。文中演示的建模技术也适用于更广泛的环境，在后续部分的课程中，我们将进一步了解这些数学接口，以及如何对更通用的偏微分方程和偏微分方程系统进行建模。

COMSOL Multiphysics^® 中的数学接口

如下图所示，您可以在模型向导的 数学 分支下找到数学接口

模型向导中的 数学 接口。

通过这些接口，您可以使用不同的公式和数值方法求解各种类型的偏微分方程，还可以使用这些接口求解常微分方程（ODE）和其他特殊类型的方程，并将不同类型的方程耦合起来，形成非线性方程系统，模拟各种数学和物理过程。除了 优化 和 灵敏度 以及 移动界面 分支下的部分接口外，所有这些接口都是软件核心功能的一部分。

有限元法和边界元法

在 COMSOL Multiphysics^® 中，求解偏微分方程的主要方法是有限元法，软件中的大部分数学和物理场接口都基于这种方法。不过，在某些情况下，使用边界元法可能更加方便，而且边界元法经常会与有限元法法耦合使用。边界元法仅需要使用表面上的网格，即不在计算体积内使用网格，这在计算体积较大的模型时非常有利。这种情况通常发生在相关对象相距较远时，举个有趣的例子，我们来看看如何求解地月系统中引力场的拉普拉斯方程和泊松方程（目前软件中还没有相应的内置接口）。地球和月球之间的距离相对于它们的半径较远，因此适合使用边界元法对它们周围的空间建模。为了简单起见，我们将假设使用稳定的经典牛顿引力。

对于静电、腐蚀、静磁和声学，软件中有专门的耦合了有限元法和边界元法公式的物理场接口。由于重力不属于这种情况，我们将通过把基于有限元法的 系数形式偏微分方程 接口与基于边界元的 偏微分方程，边界元 接口相耦合来创建自定义方程和方程耦合。偏微分方程边界元 接口的另一个好处是，它可以自动覆盖一个无限计算体积。请注意，在本文的示例中，我们将忽略其他所有天体，但可以通过添加太阳或其他天体对示例模型进行扩展。

牛顿万有引力

先介绍一些符号。我们将使用偏微分方程来模拟重力，为此，我们需要引力矢量场 ，其 SI 加速度单位为 ，在 x、y 和 z 方向上有三个分量：

为了使计算更简单，并能计算出偏微分方程的标量值，我们还需要引力势，即，标量场 , 其单位是并不怎么常见的 。稍后，您将看到如何使用 数学 接口指定自定义单位，并将其与软件中的自动单位处理功能相集成。

我们将使用下列向量微积分中的符号：

for the gradient of the scalar field.
for the divergence of the vector field.

用于计算万有引力的高斯定律 为：,
其中 是引力常数， 是质量密度。如果需要的话，我们也可以将由于地球内部分层所造成的密度变化考虑在内。不过，在此示例中，我们将使用地球和月球的平均质量密度值：5515 和 3340 。 在软件中，我们定义这些密度值为 变量，即可以将它们更改为坐标变量 x, y, 和 z 的函数, 或者，如果您更喜欢使用球坐标系，也可以将它们更改为 r， theta，和 phi 的函数。

引力场与引力势之间的关系是

将其引入方程，我们就得到了引力的经典标量方程，即如下形式的泊松方程：

引力常数的单位是 密度的单位是 。由于这两个单位相乘，所以等式右边的单位是 或 。在这个示例中，我们将使用单位 ，其实它也可以换算成 。

由于自由空间中没有质量密度， 我们可以使用拉普拉斯方程：

请注意，重力矢量场 和重力势 的关系类似于静电中的电场 和电势 ：

事实上，在模型向导或 AC/DC 的 电场和电流 分支下的 添加物理场 窗口中有专门的物理场接口。原则上，我们可以使用 静电 接口来建立重力模型，但这样做会有一些不合适的地方，例如边界条件的名称不匹配，单位也不对。

系数形式偏微分方程接口

COMSOL Multiphysics^® 中用于基于方程建模的最重要工具之一是 系数形式偏微分方程 接口，其中的模板提供了一个强大的通用接口，用于设定线性和非线性方程，包括经典的泊松方程和拉普拉斯方程。请注意，还有专门针对泊松方程和拉普拉斯方程的 数学 接口，您可以在模型向导中的 数学 > 经典偏微分方程 下找到，也可以使用这个接口。不过，本文的目的是学习如何使用通用偏微分方程接口，如 系数形式偏微分方程 接口。

假设标量场 , 用户界面中的 系数形式偏微分方程 接口的方程表示为：

这个接口有两个重要的边界条件。首先是广义 诺伊曼 边界条件：

,

是向外表面法线矢量，用于表征边界通量。第二个是 狄利克雷 边界条件：

用于对边界赋以因变量的固定值。

如果假设只有一个因变量，并且所有材料都是各向同性的，那么大部分系数都是标量，但有些系数是矢量，如下所示：

标量: 
矢量: 

我们将在后续的文章中进一步介绍这些系数。请注意，对于更常见的情况，场量和系数都可以是向量或矩阵（张量）。

以牛顿万有引力为例，您将看到如何通过简单的系数匹配，使用 系数形式偏微分方程 接口来定义一个偏微分方程。

泊松方程

在本文的示例中，只需要系数形式偏微分方程的稳态形式。换句话说，对于所有的关注时间，我们假设：

即，忽略系数  和  的影响。

在本文的示例中，我们也不需要很多其他系数，即，

最后得出等式（移动前面的负号）：

为了模拟特定引力的泊松方程形式，我们只使用系数 和 其他系数均保持默认值（零）。具体建模方程为：

于是得到

因变量：
扩散系数：
源项：

稍后，我们将学习如何使用模型向导从头开始进行设置。为了让您了解这个偏微分方程在用户界面中的显示效果，下图显示了相应的 系数形式偏微分方程 接口的截图。

使用 系数形式偏微分方程 接口模拟的重力方程。

在这个接口中，分别定义 G 和 rho 为模型中引力常数和密度的参数和变量。两个球体内 rho 的值将是地球  5515 ；和月球 3340 , 我们将用这两个球体来表示这些天体。

模型向导

首先创建一个新模型，选择 文件 > 新建。在模型向导中，选择 三维 作为空间维度，如下图所示。

在模型向导中为空间维度选择 三维。  

下一步，在模型向导中选择 数学 > 偏微分方程接口 > 系数形式偏微分方程。将 场名称 和 因变量名称 都改为 V，如下图所示。对于更通用的模型，场名称 可以用来引用因变量的数组或向量。

The Select Physics page with the Added physics interfaces field on the left and the Review Physics Interface window on the right showing the dependent variables and units.
模型向导中的 选择物理场 界面。

如上图所示，我们需要将因变量的物理量单位改为自定义单位，即重力势能的单位为 , 源项的单位为 。

单击 单位 栏中 因变量物理量 表格右侧的 定义因变量单位 按钮，然后输入  J/kg 作为单位，如图所示：

因变量物理量 的 自定义单位 设置。

单击 单位 栏 源项物理量 表格右侧的 定义源项单位 按钮。如下图所示，输入 J/(m^2*kg) 作为单位。

源项物理量 的 自定义单位 设置。

下一步，在模型向导中，选择 稳态 作为研究类型。

模型向导的 选择研究 界面。

单击 完成 进入模型开发器的三维工作区。

使用系数形式偏微分方程计算泊松方程 

我们需要为地球半径  r_earth、月球半径 r_moon，以及地球到月球的距离 d_moon 定义几个参数。有一个内置参数引力常数 G_const，为方便起见，我们将其重新定义为参数 G，如下图所示。

示例模型中使用的参数。

如下图所示，几何模型由两个相距 d_moon  的球体组成。

The COMSOL Multiphysics UI showing the Model Builder and Sphere settings window, with the geometry for the Earth–Moon system in the Graphics window.
地月系统的几何模型。

在模型中，我们把长度单位从米改为千米。需要注意的是，由于我们将使用边界元法来计算空间重力，因此不需要用体网格来包围两个球体。（注：这个模型也可以定义为二维轴对称模型。）

将地月系统模拟为两个球体，长度单位设为千米。

地球和月球的平均密度被定义为在域 1（地球）和域 2（月球）中取值不同的两个同名域变量 rho，如下图所示。

地球平均密度的变量。

月球平均密度的变量。

请注意，如前所述，您也可以为空间密度变化自定义表达式。

下一步，按照前文所述，定义 系数形式偏微分方程 的设置系数（见下图）。

The COMSOL Multiphysics UI showing the Model Builder and the Coefficient Form PDE settings window, with the Earth–Moon system model in the Graphics window.
地月系统的 系数形式偏微分方程 设置。

基于有限元的系数形式偏微分方程设置到此结束。下一步，我们将添加一个 偏微分方程，边界元 接口。

使用偏微分方程，边界元接口求解拉普拉斯方程

由于我们将把地月系统作为一个理想的孤立系统来建模，因此用不着强加不必要的边界条件。我们希望将该系统视为位于没有外部干扰的无限自由空间中，为此我们将使用 偏微分方程，边界元 接口，它将或多或少地自动处理地月系统之外（和之间）的无限计算域。

偏微分方程，边界元 接口包括以下形式的稳态偏微分方程：

请注意，这种形式没有源项，因此不能用它来求解泊松方程。不过，我们可以在自由空间中使用它，因为那里的重力势能受拉普拉斯方程控制：

在确定系数时，我们发现

因变量：

扩散系数：

吸收系数，或亥姆霍兹系数：

数学 接口列表中有 偏微分方程，边界元 接口。要添加该接口，请右键单击 组件 1，在 添加物理场 窗口中选择 数学>偏微分方程接口> 偏微分方程，边界元，如下图所示。

The COMSOL Multiphysics UI showing the Model Builder, Component settings, the Earth–Moon system in the Graphics window, and the Add Physics window with PDE, Boundary Elements highlighted.
通过 添加物理场 窗口将 偏微分方程，边界元 接口添加到模型中。

默认情况下， 偏微分方程，边界元 会选择域 1 ，域 2 以及 无限空域，其中，无限空域 是 COMSOL 术语，指球体周围直至无限远的空间。我们需要删除域 1 和域 2，因为我们只想在 无限空域 中应用 偏微分方程，边界元，如下图所示。

从几何体选择中删除域 1 和域 2 之前的 偏微分方程，边界元 接口的域选择。

从几何体选择中删除域 1 和域 2 后的 偏微分方程，边界元 接口的域选择。

下一步是将 偏微分方程，边界元 接口的单位改为与 系数形式偏微分方程 接口的单位相同，如下图所示。

偏微分方程，边界元 接口的 单位 设置。

完成这一步后，我们可以将因变量的名称改为 V，如下图所示。

偏微分方程，边界元 接口的 因变量设置。

COMSOL Multiphysics^® 中有一项功能，即通过在两个接口中使用相同的因变量名称，可以在所有共同边界上获得两个偏微分方程之间的一致耦合。要做到这一点，两个接口中的单位必须相同。

回想一下，我们的建模策略是在以球体为模型的地球和月球内部使用 系数形式偏微分方程 接口定义泊松方程，然后定义一个 偏微分方程，边界元 接口来表示周围空间的重力场。在下图中，您可以看到当前模型树的设置。

地月引力模型中使用的数学接口。

现在，我们需要指定重力势能的渐近值，或者无穷远条件。我们将在 偏微分方程，边界元 接口的 设置 窗口中将该值设为零，如下图所示。

在 偏微分方程，边界元 接口中设置 无限空域 和 无穷远条件。

如果使用牛顿万有引力和静电之间的类比，我们可以把这个渐近条件看作是远离地月系统的球体上的“接地”边界条件。这样，我们就不需要为 系数形式偏微分方程 接口指定任何实际的边界条件。相反，这个 无穷远条件 是从 偏微分方程，边界元 接口定义的 无限空域，通过空域与球体的两个域之间的边界耦合向内传播。同样，这种耦合只需将为两个方程接口指定相同的因变量名称即可。

关于球面中使用的有限元网格，如下图所示，地球使用 500 千米的固定网格单元尺寸，月球使用 200 千米的固定网格单元尺寸，这是通过设置相同的最大和最小单元尺寸值来实现的。球体内部将使用四面体有限元网格。边界元网格由球体表面的三角形网格组成，与 无限空域 相邻。

The COMSOL Multiphysics UI showing the Model Builder; the Mesh 1, Size 1 settings; and the Earth and Moon models in the Graphics window.
地球模型的有限元网格设置。

现在，我们可以选择 计算 来求解偏微分方程系统。计算时间为几分钟。

结果

我们可以绘制出表面重力的大小：

对应于两个行星体上的表达式 sqrt(Vx^2+Vy^2+Vz^2)，如下图所示。您可以在文本框中键入类似的数学表达式，用于将表达式的结果可视化。

表面图中用于显示表面重力大小的表达式。

显示两个行星体表面重力大小的表面图。

我们可以使用鼠标点击每个球体的表面，读出地球的数值为 9.81 ，月球的数值为 1.62 ，这与预期值相差不到 1%。

点击结果图中每个球体的表面后生成的显示表面重力的表格。

我们可以通过绘制 pdebe.normgrad_u 变量来显示地球和月球之间沿中心线的重力大小。pdebe.normgrad_u 是一个内置变量，用于高效计算使用边界元法求解得到的重力矢量场的大小。最终，我们得到如下图所示的图形。

在表达式框中使用内置变量的绘图设置。

使用对数比例显示地球和月球之间直线上的重力。

我们可以在大约 345000 千米处看到一个最低点，但请注意，这并不是所谓的 拉格朗日点。要计算拉格朗日点的位置，需要将旋转坐标系中因轨道运动而产生的反作用加速度效应计算在内。如下图所示，为了绘制这幅图，我们增加了 栅格 数据设置，其中增加了体栅格单元，以便在空域中可视化和评估边界元的物理量。

三维栅格 数据集的设置。

我们还可以用颜色、箭头和等高线来直观地显示引力场，如下图所示。

A slice plot of the Earth–Moon system compared to a blue color scale and Rainbow color scale with green arrows indicating the gravity direction.
在这个地月系统可视化图中，切面图显示的是空间重力值低于 0.05 m/s² 时的情况，表面图显示表面的重力，等值线显示重力强度，箭头显示重力场的方向。

例如，从这幅图中我们可以看到，地球引力在月球附近也占据了很大的面积。只有当非常接近月球时，月球的引力才会占主导地位。

最后，我们可以利用切割平面图的变形修改来绘制“重力井”，如下图所示。

地月系统重力井图。

平方反比定律

我们可能不太熟悉引力的偏微分方程，因此让我们来研究一下它与我们更熟悉的引力平方反比定律之间的关系。平方反比定律实际上包含在泊松方程中，为了解这一点，我们考虑一个质量为 ，位于原点的点源 ，其相应的“点密度”为 ，满足方程：

对于数学爱好者来说，这个等式只有弱意义。

这个方程的解析解是

其中，  是与点源（原点）的距离。

通过在球坐标中求取梯度，我们可以得到平方反比定律：

其中， 是径向的向外单位矢量。

有关球坐标的更多信息，请参阅相关 维基百科页面。

We get the force  on a particle in this gravity field, having mass , by the relationship
质量为 的粒子在该重力场中受到的力 为，

从而得出平方反比定律：

基于边界元法的内置接口

有关使用边界元法的内置物理场接口的更多信息，请参阅以下博客：

• 声学：“如何使用边界元法进行声学仿真”
• 腐蚀：“使用边界元法简化腐蚀仿真过程”
• 静磁学：“验证边界元法在静磁学仿真中的应用”
• 静电学：“如何使用导线、表面和固体建立静电模型”

进阶学习

本文使用的模型可供下载。您可以打开模型文件，详细查看地月系统引力场的拉普拉斯方程和泊松方程究竟是如何求解的。在学习中心的另一篇 文章 中，我们将进一步研究系数形式偏微分方程。