求解代数场方程| COMSOL 博客

非理想气体定律

考虑这样一种情况：我们正在运行 CFD 模拟，结果生成一个速度场 u=(u,v,w) 和一个压力场 p 。现在，假设同时我们希望求解气体中的对流和传导，以及温度场 T ：

-\nabla \cdot(k\nabla T)+\rho C \textbf{u}\cdot\nabla T=0

另外，假设密度 \rho 由理想气体定律给出：

\rho= \frac{pM}{RT}

要求解这个问题，我们可以简单地将 \rho 作为 T 的函数表达式代入传热方程，得到：

-\nabla \cdot(k\nabla T)+\frac{pM}{RT} C \textbf{u}\cdot\nabla T=0

其中，除了 T 以外，所有的量都是已知的，现在我们用 COMSOL 求解这个问题。

这很容易。理想气体定律是一种状态方程，我们可以用适合直接代换的封闭代数形式访问 \rho(T) 。

然而，对于理想气体定律不适用的情况，例如一种高压下的高分子量的气体。从原理上讲，这种状态方程可能是：

A(p+B\rho^2)(1-C\rho)-D\rho=0

如果你对这种状态方程的物理参数感兴趣，可以查看这个资源。

这是一个以 \rho 表示的三次方程，我们想求解 \rho 。原则上可以这样做，但很麻烦，而且这个方程一般会有三个根。你可以在这篇关于三次函数的学习资源中阅读如何手动求解这种方程。

关于这个任务，大家可能没有立即意识到：由于 p 是空间中的一个（假定已知的）压力场，实际上有一个三次多项式方程可以求解 CFD 域内每一个空间点。我们可以把这个方程称为 代数场方程 或 分布式代数方程 。手动求解所有这些根并不是一个非常容易管理的任务。我们可以通过 COMSOL 中的求解器来完成。该如何操作呢？

我们可以从概念上把这个代数方程看作是一个既没有空间导数也没有时间导数的偏微分方程的简并形式。在 COMSOL 环境中，这种方程有几种定义方法。例如，通过使用 系数形式偏微分方程 接口或 域常微分和微分代数方程 接口定义。使用 系数形式偏微分方程 接口，我们只需将空间和时间导数的所有系数清零，然后用稳态研究求解。使用 域常微分和微分代数方程 接口也许更容易一些，因为开始时没有空间导数项。

使用域常微分和微分代数方程接口

现在，让我们看一个基于上述非理想气体定律的一个例子。首先，在模型向导中选择 域常微分和微分代数方程 接口：

要求解一个代数方程，应该使用稳态研究。这可在模型向导的下一步中完成，并将使用默认的阻尼牛顿求解器，该求解器对于大多数求根来说是足够的。

完成模型向导后，模型树将显示如下:

默认情况下，分布式常微分方程的设置窗口看起来像这样:

代表密度的未知场变量在这里被称为 u （但是如果我们想的话，我们可以很容易地将其改为 rho )。现在可以将 d_a 系数设置为零。但是，由于我们已经选择了稳态研究，求解器将不会考虑这个方程的时间导数，所以我们也可以选择让它保持原样。在这种情况下，唯一重要的是我们在源项中输入的内容。接下来，我们将在源项中输入非理想气体定律的方程，用变量名 u 来表示密度：

剩下的就是定义系数 A ， B ， C 和 D 以及给定的压力场 p 。现在我们来做一些简化的假设。

假设我们的 CFD 域是一个简单的二维正方形单元。此外，假设压力场在空间上是变化的 p(x,y)=xy ，我们将选择用无量纲单位。在实际情况下，这个压力场当然会与这个方程一起求解，来自 CFD 仿真的结果。我们还将为这些系数赋予一些任意的数字： A=1 ， B=2 ， C=3 和 D=4 。这些都可以通过在全局定义下通过定义参数来完成，如下所示：

这样做之后，我们就可以求解这个方程，得到密度场：

这个表面图显示了整个单位矩形中每个值的三次多项式方程的解， p 。

多重根的处理

对于 p 的每个值，在每个点上，上面的图只显示了对应方程的三个根的三个可能解中的一个。基于物理上的考虑，我们可能需要对解施加额外的连续性标准。例如，如果一些根是复值的，也许可以忽略。还可能会遇到这样一种情况：在域的中间，解突然从一个物理上可实现的根变为另一个根。这可能对应一个相变，比如说从液体变成气体。

那么，如何知道我们得到的是这些解中的哪一个？一般来说，这是个很难回答的问题。如果我们在求解代数状态方程的同时运行一个瞬态模拟，那么时间历史可能会决定状态方程在空间的每个点上取哪个根。当使用稳态求解器时，找到的根可能由初始值设置决定，然后由牛顿求解器用作其迭代的起点，在求解器相对于该起点的收敛区域。此外，对于离散化方案的选择，如单元阶次和有限单元类型的选择，可能会影响求解器选择的根。

为了看看会发生什么，让我们用一个更简单的方程来做实验:

(u-2)^2-p=0

如果手动求解决这个问题，可以得到:

u=2 \pm \sqrt{p}

因此，如果我们在单位矩形上求解，我们将会在 (0,0) 得到唯一的解 u=2 ，在 (1,1) 得到解 u=1 或 u=3 。如果现在用上述同样的方法在 COMSOL 中输入这个较简单的方程，对于两个不同的初始值，我们将得到以下结果：

统一初始值为 1.9 的解。

统一初始值为 2.1 的解。

请注意，在一般情况下，我们可以使用不均匀的初始值分布，在计算域的不同部分得到不同的根：

具有一个非均匀分布初始值的解 u=1.9*(y=0.5) 。这里，网格与这个初始值分布所定义的步骤一致，并且关闭了单元间的平滑。

对于某些类型的仿真，仅仅解决控制的偏微分方程可能是不够的。COMSOL Multiphysics 能够在常规物理场接口的基础上求解代数方程，如流体流动，这对更高级的物理场仿真可能很关键。COMSOL 的求解功能并不局限于代数方程，还可以求解所谓的超越方程，包括诸如 sin() 或 cos() 等函数，甚至是涉及积分的更广义的方程。这篇博客中介绍的例子都使用了 COMSOL Multiphysics 的核心功能。

The Van der Walls state law

手动求解三次方程

求解非线性稳态有限元问题

COMSOL Multiphysics 中基于方程的建模网络研讨会视频