通过符号微分加速模型收敛

2020年 10月 29日

你是否知道 COMSOL Multiphysics® 软件内置了一个符号微分求解器,在设置和解决非线性问题时会自动启用?此功能可确保软件在求解此类模型时具有很高的可靠性,但也要付出一些代价。今天,我们将通过几个不同的示例演示此功能,并研究其优点和弊端。

利用符号微分求解简单的非线性问题

让我们从一个相对简单的、单自由度的非线性问题(隐式方程)开始:

T=Q/(1+T^2/100)

 
其中, T 表示系统的集总温度,Q 代表热载荷,上式当中的分母 (1+T^2/100) 代表系统热导率。

要求解参数 Q 的不同值,我们可以通过阻尼牛顿法找到非线性问题的解,该方法可以通过在 COMSOL Multiphysics 中写残差来完成:

r(T)=T-Q/(1+T^2/100)

 
然后将其输入全局方程 接口,如下截图所示。

COMSOL Multiphysics中的全局方程式界面的屏幕截图,其中包含用于求解非线性方程式的残差。
用于求解非线性方程的全局方程接口。

如下图所示,我们可以针对一系列不同的 Q 值来求解此方程。

绘制具有单自由度的非线性问题的解的图。
单自由度非线性问题的解。

为了便于说明,我们在下图中画出了 T 取值范围内和选定值 Q 的残差。零交叉处的点是上图的解。

绘制不同Q值的残差的图形,它表示非线性问题中的热负荷。
不同参数 Q 值的残差图。与零交叉的 T T值就是非线性方程的解。

我们还画出了相同范围 T 的和 Q 值下的残差的对应导数(雅可比矩阵)。

绘制不同热负荷Q值时残差导数的图。
不同 Q 值的残差的导数图。

我们看到,随着 Q 值的增加,雅可比矩函数非线性增强。我们也可以(对于这种简单情况)手动写出其导数:

\frac{\partial r(T)}
{\partial T}=1+\frac{2TQ}{100(1+T^2/100)^2}

 
请记住,软件能够自动通过符号微分获得此导数,并在牛顿步骤中使用它。

可以看出,正是残差第二部分的分母使该表达式有些复杂(回忆下,该项可以被认为是系统的热阻)。事实证明,如果我们愿意,可以忽略该项,通过 nojac() 算子包裹分母的该项:

R(T) =(T – Q/nojac(1+T^2/100))

并修改全局方程,如下截图所示。

COMSOL Multiphysics中的“全局方程式设置”窗口演示了如何使用nojac()运算符。
屏幕截图显示了如何使用 nojac() 算子。

结果,符号微分求解器将忽略 nojac() 算子中的所有内容。也就是说,该算子中的项对雅可比矩阵没有贡献。关于 T 的残差导数现在简化为:

\frac{\partial r(T)}{\partial T}=1

 
通过检查我们可以看到,这种简化的导数仅是 Q 为低值时的一个很好的近似。但是求解器不能对所有的 Q 值,都收敛到解。因此,我们通过避免获取某些项的导数(在这种情况下非常小)使问题变得更简单了,但是在某些情况下会使问题无解。

在 COMSOL Multiphysics® 中使用 nojac() 算子

我们考虑一下可以使用 nojac() 算子的其他几种方式。

如果我们将全局方程写为:

R(T) = nojac(T – Q/(1+T^2/100))

那么导数将恰好为零,并且这个问题根本无法求解,因此这种情况并无不有趣。

我们还可以尝试将 nojac() 仅应用于一个项:

R(T) = T – Q/(1+nojac(T)*T/100)

现在,我们的导数是:

\frac{\partial r(T)}

{\partial T}
=1+\frac{TQ}

{100(1+T^2/100)^2}

 
该导数非常接近于精确的残差,并且至少在这种情况下,即使采用雅可比矩阵公式的这种近似值,求解器将收敛于所有的 Q 值。尽管计算成本没有下降,但这里值得注意的一点是,我们可以将 nojac() 算子应用于残差项的任何子集。

那么,到目前为止我们学到了什么?

  1. 您可以使用 nojac() 算子来忽略对雅可比矩阵产生贡献的项
  2. 忽略项将影响非线性求解器的收敛
  3. 忽略项可减少计算模型的大小

但是,在这个小例子中,最后一点并不是很明显。

减少大型非线性模型的大小

现在让我们看一个更大的三维模型。我们求解了一个边长为 1cm 的正方体上的传热问题,并在此部分上应用了一个极细 网格,从而我们获得了一个较大的模型;自由度大约为 300 万。

带有参数标记的传热问题示意图。
一个用于研究极细网格的简单传热问题的示意图。

我们还将材料特性定义为温度的函数,并使用类似于我们已经看到的材料特性,通过用户自定义的表达式来定义材料的导热系数:

(1+(T[1/K])^2/100)[W/m/K]

可以直接在物理场接口中完成此操作,如下截图所示。

固体传热界面的“设置”窗口具有手动定义的导热系数设置。
手动将导热系数定义为温度的函数。

当我们使用所有默认设置解决此问题时,该模型将需要用阻尼牛顿法进行四次迭代来求解。整个过程需要大约 135s 的计算时间,并且在一台典型的台式计算机上使用超过 12GB 的内存。在求解过程中,我们还会在日志中收到一条消息,内容为:

找到非对称矩阵。

该消息是材料属性中与温度相关项的结果。也就是说,材料属性中的非线性项使雅可比矩阵不对称,因此该模型将需要更多的内存来求解。现在让我们将整个材料属性表达式包裹在 nojac() 算子中:

nojac(1+(T[1/K])^2/100)[W/m/K]

重新运行该模型,我们看到该模型仍会收敛,需要大约 10GB 内存求解,但现在需要五次默认阻尼牛顿法迭代。也就是说,内存需求下降了,但是收敛的迭代次数却增加了。最有趣的是,执行求解时间为 125s,比以前少了一点!为什么会这样?求解器需要更多的迭代来收敛,但是每次迭代花费的时间更少,因为雅可比矩阵是对称的,需要处理的数据更少。

因此,在这种情况下,尽管需要更多的迭代才能收敛到解,但使用 nojac() 算子既可以减少运行求解时间,又可以减少所需的内存!还要注意,因为使用更新的材料特性表达式,仍然可以正确评估材料特性本身,所以我们能够收敛于相同的解,

避免非局部耦合项

再看最后一个示例:如下图所示,压碎包含可压缩气体的薄壁容器。可以使用结构力学模块中可用的公式进行建模。

A model of a tank with thin walls that contains compressible gas and is subject to an external load, demonstrating how to use symbolic differentiation to accelerate model convergence.
装有可压缩气体并施加有外部压碎载荷的薄壁容器。

通过利用对称性,我们可以将模型缩减为完整模型的1/8。我们施加近似于压碎的分布载荷,以及代表内部可压缩气体的压力载荷。为了表示可压缩气体,类似于在充气超弹性密封条压缩示例中所做的,我们使用了计算变形容器体积的方法。也就是说,通过使用积分耦合算子,我们在容器的整个表面上施加均匀的压力载荷,并且该压力是基于整个容器的变形来计算的。

这种非局部耦合的作用是,容器每个部分的变形都会影响压力,而压力会影响到各处的变形,因此我们的雅可比矩阵变得非常密集。使用默认的极细网格设置求解此问题时,该模型将花费大约 12GB 的内存和 400s 的时间来求解此具有约 22000 个自由度的模型。

The Face Load Settings window showing the modified pressure as an input to the crushed tank model.
屏幕截图显示了压碎容器模型中修改后的压力。

但是,通过简单的将 nojac() 算子包裹到我们所施加的压力上,我们就可以忽略非局部耦合,从而降低内存需求和求解时间。这样,我们只需要 3GB 的内存和 10s 即可收敛到相同的解!

结论

默认情况下,COMSOL 软件将始终采用精确的雅可比矩阵,因为这可以提供最大的收敛可能性。通常,我们宁愿牺牲一些时间和内存来获得尽可能可靠的解。但是,我们确实有能力和灵活性来忽略雅可比矩阵中的选定项。如果要使用此功能,我们必须显式修改材料属性,修改基本的控制方程或在模型内编写自己的变量表达式,因此这是一个有点高级的用户功能。

当然,使用此功能将近似于雅可比矩阵,这有时会减慢甚至阻止收敛,因此在这种情况下请不要使用它。但是,如果模型应用 nojac() 算子时确实收敛,收敛于完整雅可比矩阵计算相同的解,并且可能会使用更少的时间和内存,那么在构建复杂的多物理场模型时,这是一个很好的特性。


评论 (0)

正在加载...
浏览 COMSOL 博客