求解代数场方程

2014年 1月 14日

许多用户都知道,COMSOL Multiphysics可以求解偏微分方程(PDEs)和常微分方程(ODEs)以及初始值问题。但大家可能不太了解,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 环境中,这种方程有几种定义方法。例如,通过使用系数形式偏微分方程 接口或域常微分和微分代数方程接口定义。使用系数形式偏微分方程 接口,我们只需将空间和时间导数的所有系数清零,然后用稳态研究求解。使用域常微分和微分代数方程 接口也许更容易一些,因为开始时没有空间导数项。

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

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

从模型向导中选择域常微分方程和微分代数方程接口

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

稳态研究

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

显示分布式常微分方程1稳态研究的模型树

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

分布式常微分方程的设置窗口

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

非理想气体定律

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

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

在全局定义节点下定义参数

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

密度场结果

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

多重根的处理

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

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

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

(u-2)^2-p=0

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

u=2 \pm \sqrt{p}

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

统一初始值为1.9的解
统一初始值为 1.9 的解。

统一初始值为 2.1 的解
统一初始值为 2.1 的解。

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

具有非均匀分布初始值的解

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

结语

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

扩展阅读

博客分类


评论 (2)

正在加载...
康宇 邓
康宇 邓
2023-07-18

你好,我定义参数,p=x*y,红色显示报错,怎么办?

Anran Wei
Anran Wei
2023-07-19 COMSOL 员工

您好,当在全局“参数”的定义中使用坐标变量 (如 x, y, z) 时,新版本的COMSOL软件会给出一个错误的颜色指示,但模型还是可以执行计算的。本篇博客发表时(2014年)所使用的软件版本是允许在全局参数的定义中进行该操作的,但在后续的更新版本中,建议您在“变量”中定义和坐标有关的表达式。

浏览 COMSOL 博客