基于方程的时间空间离散建模

Author Image

作者Walter Frei

2020年 9月 24日

COMSOL Multiphysics®软件的核心优势之一是可以修改计算模型中的几乎所有表达式。我们必须谨慎使用这项功能,但是可以借助它实现其他很强大的功能。在这篇博客中,我们将以一个带有平移的二维(2D)热模型为例,先将某些材料属性设置为零,然后发现此问题变为类似于求解一维(1D)瞬态模型;最后思考这项功能是如何简便、快速地实现一些类型的优化问题。

将 2 个空间维度变为 1 个空间维度和 1 个时间维度

A schematic showing a 2D stationary thermal model that is analogous to a 1D transient model, which illustrates a space-time discretization problem.
类似于 1D 瞬态模型的 2D 稳态热模型示意图。

考虑一个非常简单的矩形2D传热模型。 在没有体积热但有对流项的情况下,我们将求解温度的稳态(时间不变)传热控制方程,其表达式如下:

\rho C_p \mathbf{u} \nabla T + \nabla \cdot \mathbf{k} \nabla T = 0

这里,\rho是材料密度;C_p是比热;\mathbf{k}是热导率,用下面的对角矩阵形式表示:

\mathbf{k}= \begin{bmatrix}k_{xx} & 0\\
0 & k_{yy}\end{bmatrix}

通过扩展所有项,可以更详细地写出控制方程:

\rho C_p \left( u_x \frac{\partial T}
{\partial x}+ u_y \frac{\partial T}{\partial y}\right) + \frac{\partial}{\partial x}\left(k_{xx} \frac{\partial T}{\partial x}\right)+ \frac{\partial}{\partial y} \left(k_{yy} \frac{\partial T}{\partial y}\right)= 0

现在,我们来做一件有意思的事情。假定速度矢量完全沿着+y方向上,因此u_x = 0。将对角导热系数张量的y分量设置为零,k_{yy} = 0,以上等式简化为:

\rho C_p u_y \frac{\partial T}{\partial y} + \frac{\partial}{\partial x}\left(k_{xx} \frac{\partial T}{\partial x}\right) = 0

将导热系数张量的一项设置为零可能不常见,但这带来了一个好处:当我们写出 1D 瞬态传热控制方程时,没有体积加热或对流项,这一点就显而易见了:

\rho C_p \frac{\partial T}{\partial t} + \frac{\partial}{\partial x}\left(k_{xx} \frac{\partial T}{\partial x}\right) = 0

请注意,以上两个方程中的第一项看起来几乎相同,如果我们正确选择速度项,那么上述两个方程将相同。现在,我们看到了一些非常有趣的现象:我们的 2D 稳态模型似乎与 1D 瞬态模型相同。

实际上,这些方程是相同的,但是当我们通过数值求解时,需要注意一些事情。首先,在 2D 域的底部边界应用固定温度条件类似于瞬态 1D 情况的初始条件。其次,对于这样的 2D 模型,具有映射的网格(与y轴对齐)是有意义的,并且可以将y方向上的单元数量视为固定的时间步。这与时域模型大不相同,后者默认情况下使用自适应时间步进。而 2D 模型还将(默认情况下)包含数值稳定项,这将使二者的计算结果略有不同,除非其中一个模型达到了较高的网格划分和容差细化水平。

那么,下一个问题是:我们该怎么使用这项功能?为了了解它的功能,让我们来看一个经典的 1D 传热模型,其中一块材料被加热,如下图所示。尽管该模型看起来很简单,但是工程上关注的许多问题都可以简化到这种经典情况。当从左侧表面加热材料时,整个材料逐渐变热但不均匀。同时,左侧还存在辐射冷却和一些自由对流,近似为恒定的传热系数。

An illustration showing how transient heating problem can be reduced to a 1D model.
可以将材料板的 2D 瞬态加热模型降低为 1D 模型。

在 COMSOL Multiphysics 中,可以通过设置 2D 矩形域并使用固体和流体传热物理场接口来解决这个问题。在这个接口中,我们可以将平移运动子特征应用于固体特征,对底部边界应用温度边界条件,左侧应用热通量表面对环境辐射特征来模拟对流冷却和辐射冷却,以及添加热通量条件施加热负荷。在时间轴上具有均匀热负荷的一个模拟结果如下所示。

2D simulation results for the transient heating in a 1D material slab model, as well as the mesh.
代表 1D 平板瞬态加热的 2D 稳态模型的模拟结果和网格。

现在,假设我们要优化此加热过程。让我们尝试跟随时间来改变施加的热负荷,以使整个平板的温度在模拟结束时尽可能地接近目标温度,并设置约束使峰值温度永远不会超过设定值。通过 1D 平板的 2D 模型,我们发现这个优化问题非常容易设置。下面,让我们看看如何实现吧!

定义优化问题

COMSOL 软件定义节点下的拓扑优化分支中的密度模型特征用于定义变量控制的时间热通量。这个特征在边界上定义了一个场,该场被限制在 0 到 1 之间,表示没有加热和最大加热。该场的离散化是线性的,这意味着在每个单元上热负荷随时间可以线性变化,并且热负荷可以在ÿ方向的网格大小确定的时间跨度上变化。这并不理想——这意味着热负荷的变化与我们的时间步一样快。相反,我们希望热负荷随时间的变化慢于离散化。因此,我们需要采取下一步骤,那就是引入亥姆霍兹过滤器(Helmholtz filter)。这在拓扑优化领域中是很常见的作法,并且可以在设置节点的过滤部分实现。过滤半径需要稍大于网格尺寸,并表示一个平滑的时间。该过滤后的控制变量场的名称为dtopo1.theta,并乘以入射热通量。

A screenshot of the Topology Optimization Density Model feature Settings window.
拓扑优化密度模型特征的屏幕快照,它沿板的受热面定义,并在代表时间的整个轴上定义。

下面屏幕快照中显示了我们的目标表达式,它是一个在顶部边界(代表最终时间)的积分探针积分的表达式基于计算出的解T和想要达到的温度T_target。通过定义一个目标,即该边界上计算出的温度与目标温度之间平方差的积分,我们得到一个微分函数,当最终温度尽可能接近目标时,该函数具有最小值。我们定义的目标名称是comp1.obj,我们将从优化研究步骤中引用该目标。

A screenshot of the Settings window for the Objective Probe feature in COMSOL Multiphysics.
在顶部边界定义的目标探针特征的屏幕快照,定义了我们要最小化的目标表达式。

我们在优化接口中定义了一个约束,即:

comp1.constr < 1

如上面的屏幕快照所示,该表达式comp1.constr是通过一个全局变量探针定义的。该表达式首先取以下平均值:

\left( \left( \frac{T}{T_{max}}\right)^{p_{exp}} \right)

在加热边界上,取p_{exp}一根。这是P 范数。如p_{exp}接近无穷大时,此约束将强制温度不沿着加热边界升高到最高温度以上。现在,由于数值原因,p_{exp}不能是无穷大即不能完全满足此约束,但是如果p_{exp}足够大,则约束可以近似满足。如果p_{exp}设置过大,则我们将引入一个极度非线性的约束函数,它会减慢收敛速度,并且还会出现一些可能的数值溢出问题,因此可以将此值视为模型中的调整参数。

A screenshot of the Global Variable Probe settings for the probe used to define the constraint expression.
屏幕截图显示了定义约束表达式的探针。

作为优化模块的一部分,我们使用研究分支中的优化研究步骤引入变量、约束和控制参数。目标和约束可以是全局变量,例如上面提到的探针。下面我们来看看这种情况如何设置。

A screenshot of the Settings window for the Optimization study step.
“研究”分支中的“优化”研究步骤的屏幕快照,定义了优化求解器类型、目标、设计变量和约束。

上面的屏幕截图显示了研究分支中的优化研究步骤。在此特征下,我们定义了优化求解器、目标函数、设计变量和约束。在设置的顶部,我们看到使用了 SNOPT 求解器。该求解器期望目标函数和约束相对于设计变量是可区分的。在优化容差模型最大计算次数支配 求解器将采取多少次迭代尝试找到最优值。接下来的三个部分定义了目标函数控制变量和参数以及约束

定义了目标函数、设计变量和约束后,现在只需求解模型即可。

A graph plotting the results for design variable, filtered design variable, and temperature.
沿时间轴绘制的设计变量,过滤设计变量和温度图。

A plot comparing the temperature to the target temperature of the material slab at the final time in the simulation.
最终时间的温度(蓝色)与目标温度(黑色)比较,表明温度场在 1D 域中非常接近目标温度。

这个包括导热系数的非线性材料模型,需要约一分钟的时间优化得到上图所示的解决方案。该解决方案绘制了板加热边随时间变化的温度,还有过滤后的变量和未过滤的设计变量,以及板横截面的最终温度。我们观察到的特别有趣的现象是:优化求解器找到的解决方案为,在刚开始的时候为加热峰值,然后热载荷逐渐减小,以不超过峰值温度的限制,当整个板的温度场趋于最佳,热载荷逐渐趋于零,最后在短时间内增加了加热量,以抵消环境冷却。

这时,我们可能会问自己这种方法的优点是什么?实际上,我们也可以单纯地在时域中解决整个优化问题。那么我们使用这种方法获得了什么呢?答案是:速度和简便性。

通过重新设计为时间空间函数,我们解决了稳态优化问题,该问题比求解瞬态优化问题要快。我们还可以使用内置的拓扑优化功能,包括亥姆霍兹过滤器,可以很容易地设置随时间任意的、受约束的、平滑的强制函数。那么,除了一些概念上的复杂性之外,这种方法的缺点是会占用更多的内存。通过同时求解空间和时间,系统矩阵与时间轴的网格成比例地变大,并且该网格必须在仿真之前固定。但是,对于这种 2D 模型,即使具有非常精细的网格,内存要求也并不太高。

具有两个空间维度的模型的内存需求确实会变大,但并不是很大。下面的链接提供了一个类似的二维模型优化,其第三维度代表时间维度的案例,这个案例包括电流和传热耦合计算焦耳热,时间优化,并使用了和刚才介绍的相同的技巧,祝您建模愉快!

动手尝试

尝试使用时空模型优化随时间变化的热负荷。单击下面的按钮访问模型文件。(注意:您需要使用有效的软件许可证登录到 COMSOL Access 帐户才能下载 MPH 文件。)


评论 (0)

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