使用布辛涅斯克近似模拟自然对流

2015年 4月 7日

今天,我们将比较的 布辛涅斯克近似 与完整 纳维-斯托克斯方程 在自然对流问题中的应用。本文介绍了如何在 COMSOL Multiphysics 软件中实现布辛涅斯克近似,以及使用布辛涅斯克近似的潜在优势。

应用示例:方形空腔中的自然对流

在下面的示例中,我们将使用一个耦合了纳维-斯托克斯方程和传热方程的模型来模拟带有加热壁的方形空腔中的自然对流。空腔左壁和右壁的温度分别为 293K 和 294K;顶壁和底壁是隔热的;流体是空气,侧面的长度为 10cm。

方腔的示意图

我们将使用此模型比较三种不同建模方法的计算成本:

  1. 求解完整的纳维-斯托克斯方程(方法1)
  2.  

    \rho \left( \frac{\partial \mathbf{u}}
    {\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right) = -\nabla p + \nabla \cdot ( \mu (\nabla \mathbf{u} + (\nabla \mathbf{u})^{T}) -\frac{2}{3} \mu (\nabla \cdot \mathbf{u})\mathbf{I}) + \rho \mathbf{g}

     

  3. 用压力变换求解完整的纳维-斯托克斯方程(方法2)
  4.  

    \rho \left( \frac{\partial \mathbf{u}}{\partial t}
    + \mathbf{u} \cdot \nabla \mathbf{u}\right) = -\nabla P + \nabla \cdot ( \mu (\nabla \mathbf{u} + (\nabla \mathbf{u})^{T})- \frac{2} {3}\mu (\nabla \cdot \mathbf{u})\mathbf{I}) + (\rho-\rho_0)\mathbf{g}

     

  5. 使用带压力变换的布辛涅斯克近似(方法3)
  6.  

    \rho_{0} \left( \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u}\cdot \nabla \mathbf{u} \right) = -\nabla P + \mu \nabla^2 \mathbf{u}-\rho_0\frac{T-T_0}{T_0}\,\mathbf{g}

     

COMSOL 多物理场仿真百科中,我们介绍了这三种方法及其变量的定义。

我们使用 COMSOL Multiphysics 中的层流流体传热 接口、非等温流多物理场耦合进行静态研究来求解模型:

在模型开发器中设置模型。

在建立模型时,检查流动是层流还是湍流很重要。对于自然对流问题,我们是通过计算格拉斯霍夫数来实现的。对于理想气体,我们可以定义为
 

\mathrm{Gr} = \frac{g \left(T_\mathrm{hot}-T_\mathrm{cold}\right) L^3}{\nu^2\,T_\mathrm{cold}}

 
格拉斯霍夫数是浮力与黏性力之比。当格拉斯霍夫数低于 10^8 时表示流动为层流,而高于 10^9 时表示流动为湍流。在本文中,格拉斯霍夫数为 1.5 \times \hspace{1pt} 10^5 左右,说明流动是层流。

方法 1

使用完整的纳维-斯托克斯方程时,将浮力设置为 \rho \mathbf{g} :

在 COMSOL 中设置浮力。

使用体积力功能添加浮力项。术语 nitf1.rhog_const 分别代表与温度和压力相关的密度 \rho,以及重力加速度  \mathbf{g}

方法 2

当使用带压力变换的纳维-斯托克斯方程时,我们必须进行三个更改。

首先,需要将体积力的定义更改为 (\rho-\rho_0)\mathbf{g},例如:

定义体积力的屏幕截图。

术语 rho0 是指参考密度 \rho_0

接下来,用变量表中的材料属性计算参考密度 \rho_0 和参考黏度 \mu :

变量表格。

其中,pA 和 T0 分别代表参考温度和压力。

空气黏度被设定为常数 \mu_{0} :

将空气黏度设定为常数。

方法 3

最后,当使用布辛涅斯克近似时,需要将浮力设置为 -\rho_0\frac{T-T_0} {T_0}\,\mathbf{g}

图示表明了选择的浮力。

与方法2 一样,我们还需要根据材料属性计算参考密度和黏度。方法3 的第三步也是最后一个步骤,是将流体密度设置为恒定的参考密度 \rho_{0} (布辛涅斯克近似认为,除了在浮力项中以外,密度是恒定的)。

屏幕截图显示了特定的密度。

注意:如果模型包含压力边界条件(开放域),则将方法1 的压力设置为静水压力 -rho0 * g_const * y,将方法2 和方法3 的压力设置为 0 [Pa]。COMSOL 案例库中有一个包含重力的模型边界条件示例。

网格是由 15000 个三角形单元和 1200 个边界层单元组成的。他们都是一阶单元

网格图。

结果

以上三种方法得出的速度大小和流线几乎相同。方法1 和方法2 之间的最大温差小于 2 \times \hspace{1pt} 10^{-6} K,方法1 和方法3 之间的最大温差约为 5 \times \hspace{1pt}10^{-4} K,唯一不同的是仿真时间。

图解速度大小和流线图
速度大小和流线。

由于此二维仿真的运行时间很短(大约 30s),我们通过比较求解器从收敛到稳态解所需的迭代次数来查看计算负载。在这种情况下,迭代次数几乎与 CPU 时间成正比。

下表比较了三种方法的迭代次数。

方法1 方法2 方法3
迭代次数 39 55 55

这些结果非常令人惊讶!

虽然布辛涅斯克近似可减少模型的非线性和收敛所需的迭代次数,但完整的纳维-斯托克斯方程(39 个迭代)的求解速度仍然快于布辛涅斯克近似(55 个迭代)。我们还注意到,使用带有压力转换的纳维-斯托克斯方程可产生与布辛涅斯克近似相同的迭代次数。

为了更好地理解这些结果,我们可以在禁用伪时间步进算法之后运行第二组仿真。伪时间步进用于稳定传递问题中向稳态的收敛。伪时间步进依赖于自适应反馈调节器来控制一个柯朗-弗里德里希斯-列维条件(CFL)数。伪时间步进通常是模型收敛的必要条件。但是,在本次特定的情况下,则不需要。

以下是 COMSOL Multiphysics 中带有伪时间步进的默认求解器设置的设置窗口:

COMSOL中的默认求解器设置。

以下快照显示了没有伪时间步进的更新的 求解器设置。除非可以轻松调整求解器设置,否则我们建议始终保持开启伪时间步进。

更新后的求解器设置。

请注意自然对流的求解器设置:

由于在自然对流建模中层流与传热物理场之间存在非常强的耦合,因此请始终使用完全耦合求解器。在层流物理场中添加体积力时,COMSOL 软件会自动切换到完全耦合求解器,这意味着我们正在模拟自然对流。

下表显示了没有伪时间步进的迭代次数:

方法1 方法2 方法3
迭代次数 9 7 7

这些结果比之前使用伪时间步进时的结果更有意义。这是因为最线性问题的方法3 的收敛速度比方法1 快。令人惊讶的是,方法2 和方法3 收敛的迭代次数相同。

将这些结果与第一组结果进行比较,发现第三种方法(布辛涅斯克近似)的速度加速了 8 倍(从 55 到 7 次迭代)。这些结果还表明,第一组结果中的迭代次数不仅取决于问题的线性程度,还取决于伪时间步进算法的调整。

我们学到了什么?

在这里,我们讨论了如何实施布辛涅斯克近似及其优点,以及如何使用压力转换方法。结果表明,对于该特定模型,无论是否启用伪时间步进,在计算时间上使用布辛涅斯克近似都没有实际用处。这通常是因为布辛涅斯克近似仅在非线性较小时才有效。相对于完整的纳维-斯托克斯方程,当布辛涅斯克近似的计算时间短得多时,这表明布辛涅斯克近似可能无效。

由于使用布辛涅斯克近似法计算加速较小,并且事实上并不容易知道布辛涅斯克近似法是否有效,因此我们通常建议求解完整的纳维-斯托克斯方程。但是,实施压力变换(方法2 和 3)确实避免了舍入误差,并简化了瞬态问题以及具有开放边界的模型的实施。这将是下一篇博客文章的研究目标。

与方法2(带压力变化的纳维-斯托克斯方程)相比,使用包含更多的实施步骤的方法3(带压力变化的布辛涅斯克近似)并不会减少迭代次数。方法3 的最终仿真时间可能会稍短,因为它不需要计算与温度和压力相关的密度以及与温度相关的黏度,但这种加速可能并不明显。

通过禁用伪时间步进算法,迭代次数根据所选方法的不同会减少4到8倍。但是请记住,如果没有伪时间步进或其他负载渐变非线性渐变策略,大多数问题将不会收敛。

您可以使用CFD 模块传热模块来设置和求解该模型。如果您对我在本篇博客文章中介绍的模型有任何疑问,请联系我们的技术支持团队。如果您还不是 COMSOL Multiphysics 用户,并且想进一步了解我们的软件,请通过此表格与我们联系 – 我们将很乐意与您联系。


评论 (11)

正在加载...
辉 王
辉 王
2023-03-21

你好,我想了解三种方法的仿真对比情况,能把你的报告三种仿真方法的仿真模型分享一下吗?谢谢

Hu Chunxia
Hu Chunxia
2023-04-10

求报告中的三种模型

Haoze Wang
Haoze Wang
2023-04-26 COMSOL 员工

您好,建议您参考以下案例并结合博客中介绍的3种方法调整边界条件的设置:http://cn.comsol.com/model/buoyancy-flow-in-air-53441

Xiaochen Nie
Xiaochen Nie
2023-08-29

COMSOL可以在[达西定率]模块中添加体积力吗?我想给磁流体在多孔介质中的流动添加磁体积力(开尔文力)。
Can COMSOL add volume force to the Darcy Ratio module? I would like to add magnetic force (Kelvin force) to the flow of magnetic fluids in porous media. Looking forward to your reply.

Xiaohan Jiang
Xiaohan Jiang
2023-09-05 COMSOL 员工

达西定律中没有预定义好的添加体积力的项,如果所仿真的多孔介质域尺寸不是特别大,推荐使用布林克曼方程,改物理场中有预定义好的体积力项。或者,如果可以有确定的方程形式,也可以将包含体积力项的达西定律通过自定义PDE形式在 COMSOL 中实现。

头 石
头 石
2023-09-05

您好,达西定律和多孔介质传热模块耦合,求解密度变化浮力流(大尺寸模型,结构较复杂),考虑重力及密度变化后稳态模型不收敛,瞬态模型提示“找不到一致的初始值”。这个应该如何解决,或者应该改用什么物理场?

越 赵
越 赵
2023-09-07 COMSOL 员工

您好,自然对流问题通过较难求得一个稳态的流动状态,我们通常求解瞬态的自然对流问题。如果出现不一致的初始值问题,通常是因为边界条件和初始值的不匹配导致的,您可以参考:https://cn.comsol.com/support/knowledgebase/1172 来修改您的模型。

航青 许
航青 许
2024-09-19

请问可以通过在达西定律接口定义弱贡献的形式添加体积力项吗?

Haoze Wang
Haoze Wang
2024-09-23 COMSOL 员工

您好,如果您这里描述的达西定律中的“体积力”指的是由于孔隙流速较高后,达西速度和压降之间呈现非线性关系,则可以使用非达西流选项,参考博客:https://cn.comsol.com/blogs/modeling-darcian-and-non-darcian-flow-in-porous-media

航青 许
航青 许
2024-09-23

我想在达西定律里面添加毛细力作用,能用弱贡献的形式添加吗?或者是修改方程变量可以实现吗?

Yi Li
Yi Li
2024-09-27 COMSOL 员工

可以通过修改底层变量或弱贡献来为方程添加指定的源项,可以参考官网弱形式系列博客。

浏览 COMSOL 博客