如何求解多物理场问题

作者 Walter Frei

2013年 12月 16日

在此我们将介绍 COMSOL Multiphysics 中求解多物理场有限元问题的两类算法。到目前为止,我们已经学会了如何进行网格剖分,以及求解线性和非线性单物理场有限元问题,但是还未曾考虑过同一个域内存在多个相互依赖的不同物理场的问题。

一个简单的稳态多物理场问题

让我们从一个非常简单的稳态多物理场问题开始:耦合场中一个稳态电流穿过一根金属母线板,引起母线板内的传热和结构形变。电流通过使得电阻发热,使得母线板受热后膨胀。因为升温较为显著,所以不能忽略母线板的电、热以及结构材料性能变化。我们希望求出稳态条件下的电流、温度场、形变以及应力。下图显示了待求问题的示意图:

稳态电流经过金属母线板的多物理场问题示意图
待求的多物理场问题。

涉及的方程

在此需求解三个本构偏微分方程。首先,描述域内电压分布的方程为:

\nabla \cdot[- \sigma(T) \nabla V] = 0

通过有限元法离散后,可得到一组等式:

\mathbf{f}_V = \mathbf{K}_V(\mathbf{u}_T)\mathbf{u}_V-\mathbf{b}_V

其中下标 _{V} 表示未知电压,系数矩阵 \mathbf{K}_V 随未知温度 \mathbf{u}_T 变化而变化。假设电压分布已知,则体电阻热可从以下等式推导:

Q=\sigma(T)\bf{E \cdot E}

其中电场 \bf{E}-\nabla V 。这个热源显示在温度的本构方程中:

\nabla \cdot [ – k(T) \nabla T ] = Q(V)

由此给出了总的方程组:

\mathbf{f}_T = \mathbf{K}_T(\mathbf{u}_T)\mathbf{u}_T-\mathbf{b}_T(\mathbf{u}_V)

当我们求出域中的温度分布后便能求解结构位移:

\nabla \cdot [\mathbf{C}:(\epsilon-\epsilon_{\Delta T})] = \mathbf{0}

其中弹性矩阵 \bf{C} 经由随温度变化的杨氏模量 E(T) 推导得出。产生的热应变为 \epsilon_{\Delta T} = \alpha(T-T_0) ,应变为 \epsilon = 1/2 [{\nabla \mathbf{u}^\mathbf{T}_D + \nabla \mathbf{u}_D}] 。求解结构位移的方程组可写成:

\mathbf{f}_D = \mathbf{K}_D(\mathbf{u}_T)\mathbf{u}_D-\mathbf{b}_D(\mathbf{u}_T)

其中下标 {_D} 表示未知位移。将以上方程组合并后可得到:

\[ \mathbf{f} =
\left\{
\begin{array}{c} \mathbf{f}_V \\ \mathbf{f}_T \\ \mathbf{f}_d \end{array}
\right\}
=
\left[
\begin{array}{ccc}
\mathbf{K}_V(\mathbf{u}_T) %26 \mathbf{0} %26 \mathbf{0} \\
\mathbf{0} %26 \mathbf{K}_T(\mathbf{u}_T) %26 \mathbf{0} \\
\mathbf{0} %26 \mathbf{0} %26 \mathbf{K}_D(\mathbf{u}_T)
\end{array}
\right]
\left\{
\begin{array}{c} \mathbf{u}_V \\ \mathbf{u}_T \\ \mathbf{u}_D \end{array}
\right\}

\left\{
\begin{array}{c} \mathbf{b}_V \\ \mathbf{b}_T(\mathbf{u}_V) \\ \mathbf{b}_D(\mathbf{u}_T) \end{array}
\right\}
\]

通过检验我们可以发现这是一个非线性问题。通过 “求解非线性稳态有限元问题”博客,我们知道需要使用 Newton-Raphson 迭代法求解收敛得到结果:

\mathbf{u}^{i+1}=\mathbf{u}^i-[\mathbf{f’}(\mathbf{u}^{i})]^\mathbf{-1}\mathbf{f}(\mathbf{u}^{i})

这就是一个多物理场耦合问题的求解模式!在求解单物理场和多物理场的非线性问题之间并没有任何概念上的差异。所有涉及求解单物理场非线性问题的因素,包括所有关于阻尼、负载、非线性加速以及网格剖分的讨论,都同样适用于求解多物理场问题。

全耦合法

但是也需要理解上述方法的一个(在某些情况下非常致命的)缺点。通过代入Newton-Raphson 迭代法,我们要求解导数 \mathbf{f'(u}^{i}) ,将其写成以下的形式:

\mathbf{f}'(\mathbf{u}^i)
=
\left[
\begin{array}{ccc}
\mathbf{ K}_V%26\mathbf{K}_{V,T}%26\mathbf{0} \\
\mathbf{-b}_{T,V}%26\mathbf{ K}_T%26\mathbf{0} \\
\mathbf{0}%26(\mathbf{K}_{D,T}-\mathbf{b}_{D,T})%26\mathbf{K}_D
\end{array}
\right]

其中使用逗号求导符号,例如,\mathbf{K}_{V,T}=\partial \mathbf{K} _{V}(\mathbf{u}_T)/\partial \mathbf{u}_T

显而易见,以上矩阵为非对称矩阵,这将导致一个问题:如果系统矩阵不能确定,我们可能需要使用更加耗费内存的直接求解器求解。(尽管迭代求解器在预条件选择正确的情况下能处理范畴更为宽广的问题,但并不能保证足以应对所有特例。)使用直接求解器解决此类多物理场问题将会耗费更多的内存和时间。

然而还是存在有一种替代方法。以上的求解方式被称为全耦合法,假设需要同时考虑物理场之间的所有耦合。实际上,在求解很多类型的多物理场问题时,我们可以忽略这些非对角线的项,采取更节约内存、更具时效性的分离求解法。

分离法

分离法顺序分析每一个物理场,使用之前求解物理场得到的结果来得到一个待求物理场的载荷和材料属性。因此,对于上面这个例子,我们先通过 Newton-Raphson 迭代求解电压:

\mathbf{u}^{i+1}_V=\mathbf{u}^i_V-[\mathbf{K}_V(\mathbf{u}^i_T)]^\mathbf{-1} \mathbf{b}_V

其中对于第一次迭代,我们必须对电压和温度各做出一个初始估计值 ( \mathbf{u}_V^{i=0} , \mathbf{u}_T^{i=0} )。使用温度场的初始条件来求解得到初始的材料性能。接下来求解温度:

\mathbf{u}^{i+1}_T=\mathbf{u}^i_T-[\mathbf{K}_T(\mathbf{u}^i_T)]^\mathbf{-1} \mathbf{b}_T(\mathbf{u}_V^{i+1})

其中,在第一次迭代 i=0 中 ,使用温度场的初始条件来对材料性能 \mathbf{K}_T(\mathbf{u}_T^{i=0}) 求解,但是荷载是通过之前得到的电压的解 \mathbf{b}_T(\mathbf{u}_V^{i=1}) 计算得到。按照类似的方法,可以求解位移场:

\mathbf{u}^{i+1}_D=\mathbf{u}^i_D-[\mathbf{K}_D(\mathbf{u}^{i+1}_T)]^\mathbf{-1} \mathbf{b}_D(\mathbf{u}_T^{i+1})

其中结构的载荷以及材料性能是由先前计算出的温度场推导而来。

然后继续使用迭代:电压、温度和位移被有序地重复计算。这个算法将持续进行,直到达到收敛为止,如“求解非线性稳态有限元问题 ”博客所述

这个方法最大的优势在于可以在每一个子步中使用最优求解器。这不仅意味着您将在每一个子步求解更小的问题,还可相应的选择更节约内存的快速求解器。虽然分离法一般需要更多步骤的迭代才可达到收敛,但是与全耦合法中的迭代相比,分离法的每一次迭代在耗时上都会有显著的减少。

以下为分离式求解器处理一个具有 n 个多物理场模型的算法:

  1. 设定模型中所有物理场的初始条件
  2. 初始化记录迭代次数的计数器
  3. 以分离法的求解顺序求得第一个物理场,使用上一步骤来计算材料性能
  4. 利用已知解的来相对应地求解第二个物理场
  5. 使用第(n-1)重已知解求解第 n 重物理场
  6. 重复步骤 2-6 直到达到收敛,或迭代次数超过预期峰值

对于一般的物理场问题,用户仍需确定物理场的求解顺序,但是软件会提供为内置的多物理场接口提供合适的求解顺序的默认建议。 COMSOL Multiphysics 会按照分离法的求解顺序为每一个物理场提供默认的线性求解器设定。

当使用分离法时,收敛结果和使用全耦合法得到的解完全一致。通常分离法会经过更多次的迭代以达到收敛,但是每个子步的内存消耗和耗时会降低,从而进一步降低整体的求解耗时和内存使用率。

对求解多物理场问题的总结

在这篇博客中我们介绍了两类处理多物理场问题的算法一全耦合法以及分离法 。全耦合法与求解单物理场非线性问题的 Newton-Raphson 迭代法一致,它非常占用内存,但是十分实用,广泛适用于各物理场之间存在强相互作用的多物理场问题。在另一方面,分离法假设能对每一个物理场单独求解,且该物理场能在多重物理场间完成迭代直至模型达到收敛。

博客分类


评论 (3)

正在加载...
欣润 吕
欣润 吕
2018-08-13

请问”全耦合法”下面的”f(ui)”是否应改为”f'(ui)”?
此外f'(ui)的矩阵表达式中,第一行第二列是否应改为KV,T · uV ? (即是否漏掉了uV);
以及第三行第二列是否应改为 KD,T · uD-bD,T ?(即是否漏掉了uD)。
谢谢!

Kaixi Tang
Kaixi Tang
2023-03-20 COMSOL 员工

你好,对于你提到的问题,我理解博客中的相关公式是没有问题的。
具体来说,对于你提到的的矩阵,其表示的是f(u)求导,而非是对其求逆:其中第一列,表示的是三个控制方程依次对“uv”求导;第二列表示三个控制方程依次对“uT”求导;第三列亦然。比如你提到的第三行第二列,对应的是力学控制方程对“uT”求导,所以得到的是KD,T-bD,T(KD、bD与uT有关,遵循链式法则)。

Tengyue Gao
Tengyue Gao
2018-10-19

吕欣润,您好!
感谢您的评论。
模型相关的问题,请您联系我们的技术支持团队:
在线支持中心:www.denkrieger.com/support
Email: support@comsol.com
谢谢!

浏览 COMSOL 博客
Baidu
map