• 返回主页
  • 第17章:多相流

    翻译:刘芹芹
    校对:李娇娇(北京计算科学研究中心)

    本章讨论ANSYS Fluent的通用多相流模型。首先简单介绍多相流模型,离散相一章简要讨论拉格朗日离散相模型,凝固和融化一章讨论凝固和融化模型。关于ANSYS Fluent中如何使用多相流模型的信息,可参见《ANSYS Fluent用户指南》中的多相流模拟部分。以下小节介绍多相流模型的各种理论知识:

    17.1. 介绍

    17.2. 选择多相流模型

    17.3. VOF模型

    17.4. Mixture模型

    17.5. 欧拉模型

    17.6. 湿蒸汽模型

    17.7. 多相流中的质量传输

    17.8. 多相流中的组分传输

    17.1. 介绍

    很多在自然界和工程中遇到的流动是多相混合的流动。物理学中物质的‘相’指气相、液相、固相,而多相流中‘相’的概念更加广义。在多相流中,‘相’不仅可以定义为不同类型的物质,也可定义为相同类型物质。比如,某固相颗粒中,不同尺寸的固体颗粒也可视为不同的相。

    本节分为以下两部分:

    17.1.1. 多相流模式

    17.1.2. 多相流例子

    17.1.1. 多相流模式

    多相流主要分为五类:气-液,液-液、气-固、液-固、三相流。

    17.1.1.1. 气-液或液-液流

    下面是不同气-液或液-液流动的分类:


    图17.1:多相流流态示意图

    17.1.1.2. 气-固流

    下面是不同气-固流动的分类:

    17.1.1.3.液-固流

    以下是液-固流:

    17.1.1.4.三相流

    三相流是前面列出的几种流动模式组合

    17.1.2. 多相流例子

    下面列出多相流的一些具体例子:

    17.2. 选择多相流模型

    解决任何多相流问题的第一步,就是确定多相流模型。模型比较一节中提供了一些基本原则。具体指导原则一节中给出了详细的方法:如何确定流动中(包含气泡、液滴、或颗粒)各相之间的耦合程度,以及不同耦合程度适合的多相流模型。

    本节分为以下四部分:

    17.2.1. 多相流建模

    17.2.2. 模型比较

    17.2.3. 多相流的时间离散格式

    17.2.4. 稳定性和收敛性

    17.2.1. 多相流建模

    计算流体力学的发展为进一步深入了解多相流动力学提供了基础。目前多相流的数值计算有两种方法:欧拉-拉格朗日方法和欧拉方法;李东岳:原文中为欧拉-欧拉方法,为防止其和双流体模型混淆,此处译为欧拉方法

    17.2.1.1. 欧拉方法

    在欧拉方法中,不同的相被处理成互相贯穿的连续介质。各相的体积不能被其他相占有,因此引入‘体积分数’的概念。其假定体积分数是空间和时间的连续函数,所有相的体积分数和等于$1$。每一相都有各自的控制方程,且所有相的这些方程形式相同。另外附加一些经验性的关系式来使这些方程封闭。

    ANSYS Fluent提供了三种基于欧拉方法的多相流模型:VOF模型、Mixture模型、和双流体模型。

    17.2.1.1.1. VOF模型

    VOF模型是一种网格固定的的表面跟踪技术。该模型用于观察两种及以上互不相融流体间的分界面。VOF模型中,两种流体共用一组动量方程,计算域中各流体的体积分数在每个计算单元上被跟踪。VOF模型的应用场合有:分层流、自由面流动、灌注、晃动,液体中大气泡的流动、水坝决堤时的水流、任意液-气的稳态或瞬态分界面问题。

    17.2.1.1.2. Mixture模型

    Mixture模型用于两相或多相(流体或颗粒)的混合模拟。和双流体模型一样,所有相被处理成相互贯穿的连续介质。Mixture模型求解混合物动量方程,并通过相对速度来描述离散相。Mixture模型的应用场合主要有:低负载的粒子负载流、气泡流、沉降以及旋风分离器。Mixture模型也可用于模拟离散相没有相对速度的均匀多相流。

    17.2.1.1.3. 双流体模型

    双流体模型是ANSYS Fluent中最复杂的多相流模型。该模型中的每一相都具有一组动量方程和连续性方程。各相之间通过压力和相间交换进行耦合,耦合的处理方式取决于流动中相的类型。比如颗粒流和非颗粒流的处理方式就不同。对于颗粒流,是通过运动学理论获得相间的耦合特性。相间的动量交换也取决于流动中相的类型。ANSYS Fluent的用户可以通过自定义函数(UDF)个性化定制动量交换的计算方式。双流体模型的应用场合有:鼓泡床、上浮、颗粒悬浮、以及流化床。

    17.2.2. 模型比较

    一旦用户确定使用欧拉方法下的模型处理实际多相流问题,可以基于下述原则进一步选择合适的多相流模型:

    如本节中讨论过的,VOF模型适用于分层/自由表面流动,Mixture模型和欧拉模型适用于相间混合或分离、或离散相体积分数超过10%的流动(如果离散相体积分数小于等于10%,可用离散相模型。

    到底是选择Mixture模型还是双流体模型,用户可以考虑以下方面:

    ANSYS Fluent中的多相流模型和动网格兼容,关于动网格,可详见应用动网格的流动一节。关于ANSYS Fluent中其他模型与多相流模型的兼容性,可参考用户指南中的附录A。

    17.2.2.1.详细指南

    对于分层流和活塞流,如模型比较一节中所述,毫无疑问应选择VOF模型。然而对于其他流态不是很明确的流动,需要定义一些参数来辅助。一般来说,可使用参数颗粒负载率 ,以及Stokes数来进行判断选择合适的模型(注意,此处的‘颗粒’一词适用于颗粒、液滴、气泡)。

    17.2.2.1.1颗粒负载率的影响

    颗粒负载率对相间作用有着重要的影响。颗粒负载率被定义为离散相$d$与连续相$c$的质量密度比:

    \begin{equation} \beta=\frac{\alpha_\mathrm{d}\rho_\mathrm{d}}{\alpha_\mathrm{c}\rho_\mathrm{c}} \end{equation}

    物质密度比:

    \begin{equation} \gamma=\frac{\rho_\mathrm{d}}{\rho_\mathrm{c}} \end{equation}

    对于气-固流,物质密度比通常大于1000;对于液-固流,物质密度比通常在1左右;对于气-液流,物质密度比通常小于0.001。

    利用这些参数,可以估算颗粒相中颗粒与颗粒之间的平均距离。Crowe et al.给出的一种距离估算方法:

    \begin{equation} \frac{L}{d_\mathrm{d}}=\left(\frac{\pi}{6}\frac{1+\kappa}{\kappa}\right)^{1/3} \end{equation}

    式中$\kappa=\beta/\gamma$。这些参数信息对确定离散相处理方式十分重要。比如,对于颗粒负载率为1的气-固流,相间距离$L/d_\mathrm{d}$大约为8左右;因此颗粒可彼此视为孤立的(即非常低的负载率)。

    根据计算的颗粒负载率,相间耦合程度可分为以下三类:

    17.2.2.1.2.Stokes数的影响

    对于中等负载率的系统,需要估算Stokes数来选择最合适的模型。Stokes数定义为颗粒响应时间与系统响应时间之比:

    \begin{equation} St=\frac{\tau_d}{t_s} \end{equation} 式中的$\tau_d=\rho_\mathrm{d}d_\mathrm{d}^2/(18\mu_\mathrm{c})$,$t_s=L_s/\mathbf{U}_s$。这里的$t_s$指得是基于所研究系统的特征长度$L_s$和特征速度$\mathbf{U}_s$的比值。

    当$St\ll 1.0$的情况下,颗粒紧密跟随主流,离散相模型、Mixture模型、双流体模型都适用;用户可选择计算资源消耗最小的模型(大多数情况下为Mixture模型),或者根据其他因素选择最合适的模型。当$St>1.0$,颗粒运动将独立于主流运动,需选择离散相模型或欧拉模型。当$St\approx 1.0$,三种模型也可任选其中之一。用户可以根据计算资源消耗的大小或者其他因素选择最合适的模型。

    17.2.2.1.2.1.例子

    某煤粉分离器的特征长度为$1$m,特征速度为$10$m/s,当颗粒直径为$30\mu $m时,Stokes数值为$0.04$。而当颗粒直径为$300\mu$m,Stokes数等于$4$。显然后者不能使用Mixture模型。

    某矿物处理系统的特征长度为$0.2$m,特征速度为$2$m/s,当颗粒直径为$300\mu$m时,Stokes数等于$0.005$。这种情况下,用户可以选择Mixture模型和双流体模型。由于这种情况下相体积分数太高而不能选择离散相模型(原因如下)。

    17.2.2.1.3. 其他考虑因素

    离散相模型只限于在颗粒体积分数低的情况下使用(稠密离散相模型不受此限制)。此外,和双流体模型相比,使用离散相模型模拟时,用户可以附加燃烧模型。为获得颗粒分布,用户也可以外挂PBM模型(详见群体平衡模型指南)。

    17.2.3. 多相流的时间格式

    为了精确的模拟多相流,有必要选择高阶的空间离散格式和时间离散格式。ANSYS Fluent中除了一阶时间格式外,Mixture模型、双流体模型以及隐式的VOF模型都提供了二阶时间格式。


    注意:

    显式格式的VOF模型不能使用二阶时间格式。

    二阶时间格式适用于所有的输运方程,包括动量传输方程、能量传输方程、组分传输方程、湍流模型、相体积分数方程、压力修正方程、标量传输方程等。在多相流中,通用的传输方程可写为:

    \begin{equation} \frac{\partial \alpha\rho\Phi}{\partial t}+\nabla\cdot\left(\alpha\rho\mathbf{U}\Phi\right)=\nabla\cdot\boldsymbol{\tau}+S_\Phi \label{1} \end{equation}

    式中$\Phi$为混合相(对于Mixture模型)或某一相的变量,$\alpha$为相体积分数(Mixture模型中$\alpha=1$),$\rho$为混合相(对于Mixture模型)或某一相的的密度(李东岳:原文为“$\rho$为混合相密度”,值得商榷);$\mathbf{U}$为混合相(对于Mixture模型)或某一相的速度;$\boldsymbol{\tau}$为应力项;$S_\Phi$为源项。

    二阶精度时间格式为全隐式格式,离散后的方程($\ref{1}$)为(李东岳:原公式有误):

    \begin{equation} \frac{3(\alpha_\mathrm{P}\rho_\mathrm{P}\Phi_\mathrm{P}V_\mathrm{P})^{n+1}-4(\alpha_\mathrm{P}\rho_\mathrm{P}\Phi_\mathrm{P}V_\mathrm{P})^{n}+(\alpha_\mathrm{P}\rho_\mathrm{P}\Phi_\mathrm{P}V_\mathrm{P})^{n-1}}{2\Delta t}=\sum_\mathrm{N}\left(A_\mathrm{N}(\Phi_\mathrm{N}-\Phi_\mathrm{P})\right)^{n+1}+S_U^{n+1}-S_p^{n+1}\Phi_\mathrm{P}^{n+1} \end{equation}

    其化简形式为:

    \begin{equation} A_\mathrm{P}\Phi_\mathrm{P}=\sum_\mathrm{N}A_\mathrm{N}\Phi_\mathrm{N}+S_\Phi \end{equation}

    其中

    \begin{equation} A_\mathrm{P}=\sum_\mathrm{N}A^{n+1} + S_\mathrm{P}^{n+1}+\frac{1.5(\alpha_\mathrm{P}\rho_\mathrm{P}V_\mathrm{P})^{n+1}}{\Delta t} \end{equation} \begin{equation} S_\Phi=S_U^{n+1}+\frac{2(\alpha_\mathrm{P}\rho_\mathrm{P}\Phi_\mathrm{P}V_\mathrm{P})^{n}-0.5(\alpha_\mathrm{P}\rho_\mathrm{P}\Phi_\mathrm{P}V_\mathrm{P})^{n-1}}{\Delta t} \end{equation}

    该时间格式基于ANSY Fluent已有的一阶欧拉格式,因此比较容易使用。虽然该时间格式是无条件稳定的,但由于$^{n-1}$时间步的系数可能为负,因此如果时间步长太大,容易产生震荡。

    这个问题可以通过引入有界的二阶格式消除。但通常只在可压缩流中产生震荡,所以ANSYS Fluent仅在可压缩流中应用有界的二阶格式。因此对于单相或多相的可压缩流,二阶时间格式默认是有界的二阶格式。

    17.2.4. 稳定性和收敛性

    求解多相流的过程本身具有难度,用户可能会遇到稳定性和收敛性问题。

    对于稳态求解:推荐选择多相耦合算法(Multiphase Coupled Solver),详见《ANSYS Fluent用户指南》中的欧拉多相流耦合求解器介绍。根据耦合算法的迭代性质,其要求具有一个比较好的初始场。如果因为高阶格式、或系统复杂性导致求解困难,用户可以降低库郎数,默认的库郎数为$200$,最低可降至$4$。当迭代过程顺利进行之后,再增加库郎数。此外,速度和压力采用显性松弛,其他变量采用隐性松弛。对于耦合求解器,如果将体积分数方程的松弛因子设置较低的值,可能会导致求解过程显著变慢(为提高计算速度,松弛因子必须在$0.5$以上),但对于SIMPLE类算法则相反,通常要求体积分数方程的松弛因子取较低值。

    对于非稳态求解:合适的初始场可以避免求解不稳(求解不稳通常是因为不合适的初始场而引起的);如果担心计算时长,那么最好是采用SIMPLE算法。当体积力比较重要的时候,或者需要调用高阶数值格式的情况下,推荐先采用小的时间步长,在计算少许时间步后(得到较好的压力场后)再增加时间步长。

    当使用PISO算法(李东岳:原文为Non-Iterative Time Advancement,NITA,此处理解为PISO算法)计算非稳态流动时,好的初始场十分重要。比如,网格质量差或者有大的体积力,可能会引起求解不稳。

    如果在稳态计算、或拟稳态计算中使用MRF模型遇到收敛性问题,用户可以切换成非稳态求解器,其同样可以收敛到一个稳定解。对于MRF模型,当使用PISO算法时,用户应该注意在MRF区域的边界处,网格质量差或动量方程中一个比较大的源项可导致PISO算法稳健性问题。对于MRF问题,最好使用瞬态SIMPLE算法(李东岳:原文为Iterative Time Advancement,ITA,此处理解为瞬态SIMPLE算法)迭代类时间步进法,其能对每个时间步提供了更多的迭代次数控制项。

    此外,ANSYS Fluent还提供了一个全多相耦合求解器(Full Multiphase Coupled solver),可以同时求解所有动量方程、压力修正方程和体积分数方程,只是目前该求解器稳健性还不如其他求解器。

    此外,ANSYS Fluent还有一个选项用于在双流体框架下中求解互不相融的分层流。其主要特征是具有多个速度分布,其他特征则和单一速度分布的VOF模型相似。

    17.3. VOF模型理论

    本节分为以下几部分:

    17.3.1. VOF模型综述

    17.3.2. VOF模型局限性

    17.3.3. VOF稳态和瞬态计算

    17.3.4. 体积分数方程

    17.3.5. 物质属性

    17.3.6. 动量方程

    17.3.7. 能量方程

    17.3.8. 其他标量方程

    17.3.9. 表面张力和依附

    17.3.10. 明渠流

    17.3.11. 明渠造波边界条件

    17.3.12. Level-Set和VOF耦合模型

    17.3.1. VOF模型综述

    VOF模型是通过求解一组动量方程来模拟两种及两种以上互不相融的流体,并在整个计算域跟踪各流体的体积分数。典型应用场合包括:射流破裂、液体中大气泡的流动、溃坝、任意液-气的稳态或瞬时分界面问题。

    17.3.2. VOF模型局限性

    ANSYS Fluent的VOF模型具有以下局限性:

    17.3.3. VOF稳态和瞬态计算

    通常,ANSYS Fluent中的VOF模型只用于瞬态计算,如果用户只关心求解问题的稳态解,也可执行稳态计算。但是只有在求解问题明确的不依赖于初始条件,而且各相具有明确的边界前提下,VOF稳态计算才是合理的。比如,旋转杯子中的自由表面形状严重依赖于杯子内部初始的液体高度,此时必须选择瞬态计算。另一方面,渠道中的水流,可以选择稳态计算。

    VOF模型的公式是基于各相之间互不渗入这一事实。模型中每增加一个相就引入另外一个变量,即该相的体积分数。在每个网格单元内,所有相的体积分数和等于$1$。流场中的所有变量及属性均为各相共享,只要知道每个网格单元中各相的体积分数,就能用体积平均求出这些变量和属性。因此,任何网格单元的变量及属性是仅仅代表一个相,还是代表多个相的混合,完全取决于相的体积分数值的大小。换句话说,如果某单元内第$i$相的体积分数记为$\alpha_i$,那么有以下三种情况:

    根据局部的$\alpha_i$值,计算域内的所有网格单元被赋予了合适的属性。

    17.3.4. 体积分数方程

    ANSYS Fluent通过求解一个相或多个相的体积分数的方程,实现对相界面的跟踪。对第$q^{th}$相,方程具有以下形式:

    \begin{equation} \frac{\p}{\p t}\alpha_\rq\rho_\rq+\nabla\cdot(\alpha_\rq\rho_\rq\bfU_\rq)=S_\rq+\sum_{p=1}^n(\dot{m}_\mathrm{pq}-\dot{m}_\mathrm{qp}) \end{equation}

    式中$\dot{m}_\mathrm{qp}$为从第$q$相到第$p$相的质量传递,$\dot{m}_\mathrm{pq}$为从第$p$相到第$q$相的质量传递。默认情况,$S_\rq$为0。但用户也可以对每相指定一个常数源项或使用用户自定义源项,详见通用多相流模型中多相流中的质量传递一节。

    主相的体积分数方程不用求解;它可以根据下面的约束条件公式计算:

    \begin{equation} \sum_{q=1}^n \alpha_\rq=1 \end{equation}

    体积分数方程的时间格式可以为隐式或显式。

    17.3.4.1. 隐式

    当使用隐式时间格式的时候,体积分数方程可离散如下(李东岳:原式有误):

    \begin{equation} \frac{(\alpha_\rq\rho_\rq)^{n+1}-(\alpha_\rq\rho_\rq)^{n}}{\Delta t}V+\sum_f(\rho_\rq\bfU_f\alpha_{\rq,f})^{n+1}=S_\rq V+V\sum_{p=1}^n(\dot{m}_\mathrm{pq}-\dot{m}_\mathrm{qp}) \end{equation}

    式中:

    $n+1$表示当前时间步;

    $n$表示前一个时间步;

    $\alpha_\rq^{n+1}$表示在$n+1$时间步的网格单元体心的体积分数值;

    $\alpha_\rq^{n}$表示在$n$时间步的网格单元体心的体积分数值;

    $\alpha_{\rq,f}^{n+1}$表示在$n+1$时间步的网格单元面心的体积分数值;

    $\bfU_{f}^{n+1}$表示在$n+1$时间步的网格单元面心的速度;

    $V$表示网格单元的体积;

    当前时间步的相体积分数是当前时间步其他相变量的函数,因此,求解任意时间步的某一次相的体积分数时需要对相体积分数方程进行迭代求解。

    面心的通量(速度)可通过选择空间离散格式进行插值获得。ANSYS Fluent中隐式时间格式可以使用的空间插值格式可参考《ANSYS Fluent用户指南》中体积分数的空间插值格式一节。

    隐式时间格式对瞬态求解和稳态求解都适用。详见《ANSYS Fluent用户指南》中选择体积分数公式一节。

    17.3.4.2. 显式

    显式情况下的体积分数方程可离散如下(李东岳:原式有误):

    \begin{equation} \frac{(\alpha_\rq\rho_\rq)^{n+1}-(\alpha_\rq\rho_\rq)^{n}}{\Delta t}V+\sum_f(\rho_\rq\bfU_f\alpha_{\rq,f})^{n}=S_\rq V+V\sum_{p=1}^n(\dot{m}_\mathrm{pq}-\dot{m}_\mathrm{qp}) \end{equation}

    式中:

    $n+1$表示当前时间步;

    $n$表示前一个时间步;

    $\alpha_\rq^{n+1}$表示在$n+1$时间步的网格单元体心的体积分数值;

    $\alpha_\rq^{n}$表示在$n$时间步的网格单元体心的体积分数值;

    $\alpha_{\rq,f}^{n+1}$表示在$n+1$时间步的网格单元面心的体积分数值;

    $\bfU_{f}^{n+1}$表示在$n+1$时间步的网格单元面心的速度;

    $V$表示网格单元的体积;

    当前时间步的体积分数基于前一时间步的已知量,因此可直接计算,显式公式在每个时间步不需要对标量输运方程进行迭代求解。

    面心的通量(速度)可使用相对应的插值格式进行计算,比如几何重构格式(Geo-Reconstruct)、CICSAM格式、压缩格式(Compressive)以及修正HRIC格式(见界面附近的插值一节)。ANSYS Fluent中显式公式可以使用一些插值格式详见《ANSYS Fluent用户指南》中体积分数的空间离散格式。

    ANSYS Fluent可对求解体积分数方程的时间步长进行自动优化,但是用户可以通过修改库朗数来进一步限定时间步长。用户可以选择在每个时间步内更新体积分数,也可以选择在每个迭代步后更新。更多选项详见《ANSYS Fluent用户指南》中显式体积分数方程的时间依赖性参数设置一节。


    注意:

    使用时间显式,必须选择瞬态求解。

    17.3.4.3.界面附近的插值

    ANSYS Fluent中采用守恒的离散格式。

    在ANSYS Fluent的几何重构和施主受主(donor-acceptor)方法中,对相间界面附近的单元采用一种特殊插值处理。


    图17.2:界面重构示意图

    两相流中显式格式和隐式格式对网格单元面变量的插值格式,和单相流相同,如一阶迎风格式、二阶迎风格式、QUCIK格式、修正的HRIC格式、compressive格式或者CICSAM格式等。

    17.3.4.3.1. 几何重构格式

    在几何重构格式中,对于完全由单一相充满的网格单元,采用ANSYS Fluent中的标准插值方法计算表面通量;而对于相界面附近的网格单元,采用几何重构格式计算面通量。

    几何重构格式使用分段线性的方法处理界面。ANSYS Fluent中几何重构格式是最准确的,其适用于非结构网格。Youngs的工作使得基于非结构网格的几何重构格式得以推广。该方法假定网格单元的相界面是线性的,且利用该线性函数计算通过网格单元的面对流通量(参见图17.2)

    几何重构格式第一步:根据界面处网格单元的体积分数及其导数信息,计算这些网格单元中界面(线性的)的位置;第二步:根据上步得到的线性界面和表面上法向、切向速度分布,计算通过每个网格单元面的对流通量;第三步:根据通量守恒计算每个单元的体积分数。


    注意:

    使用几何重构格式,必须选择瞬态求解。且用户需要保证计算域中没有无厚度的面。如果存在这种无厚度的面,用户需要把它们分割开,详见《ANSYS Fluent用户指南》中面区域分割一节。

    17.3.4.3.2. 施主受主格式

    在施主受主格式中,对于完全由单一相充满的网格单元,采用ANSYS Fluent中的标准插值方法计算表面通量;而对于相界面附近的网格单元,采用施主受主格式计算面通量。这种方法认为流体从网格单元(施主)中流出,毗连的网格单元(受主)接受这部分流体。流出流入的流体具有相同的量因此这种方法可避免界面处的数值耗散。通过网格单元面的流体量受限于施主网格单元中可充满的最大相体积以及受主网格单元中最大可接受的相体积。

    界面的方向也用于确定网格单元的面通量。重构的界面方向要么水平、要么垂直,这取决于该网格单元以及毗连网格单元中的第$_\rq$相的体积分数的梯度。除了界面方向以及运动可以影响网格单元的面通量计算外,其还取决于各种不同的速度插值格式如迎风格式、逆风格式或结合两种格式来求解流量。


    注意:

    使用施主受主方法,必须选择瞬态求解。且用户需要保证计算域中没有无厚度的面。如果存在这种无厚度的面,用户需要把它们分割开,详见《ANSYS Fluent用户指南》中面区域分割一节。

    17.3.4.3.3. CICSAM

    CICSAM由Ubbink and Issa提出,其是一种高分辨率的插值方法。CICSAM特别适用相间高粘度比的流动。ANSYS Fluent中CICSAM作为一种显式格式,优点是能产生和几何重构格式一样的清晰界面。

    17.3.4.3.4. 界面压缩格式以及基于界面的变量

    界面压缩格式是一种基于限制器的二阶几何重构格式。

    \begin{equation} \alpha_f=\alpha_\rP+\beta\nabla\alpha_\rP \end{equation}

    $\alpha_f$为网格单元面的相分数;

    $\alpha_\rP$为施主单元的相分数;

    $\beta$为限制器;

    $\nabla\alpha_\rP$为施主网格单元相分数的梯度值;

    不同限制器的值对应与不同程度的界面压缩,如下表所示:

    限制器的值格式
    0一阶迎风
    1二阶重构,且有界与最小/最大相分数
    2压缩格式
    0-1一阶格式和二阶格式混合的格式
    1-2二阶格式和压缩格式混合的格式

    表17.1:限制器以及对应的格式

    压缩格式可用于不同类型的界面。如果开启尖锐界面选项,压缩格式仅适用于重构尖锐界面。如果开启尖锐/混合界面选项,压缩格式可重构尖锐界面也可重构相较为分散的混合界面。

    17.3.4.3.5. 有界梯度最大化格式(BGM)

    VOF模型结合BGM可以获得清晰的重构界面(和几何重构方法一样尖锐)。目前BGM仅用于稳态计算,不能用于瞬态计算。BGM力求实现相分数的梯度最大化来实现界面的尖锐 by maximizing the degree to which the face value is weighted towards the extrapolated downwind value(李东岳:此处意义不明

    17.3.5. 物性

    传输方程的物性取决于每个网格单元内的相的组成。比如在两相流系统中用下标$_1$和$_2$分别表示各相,如果要跟踪第一相的体积分数,那么每个网格单元的密度可根据下式计算:

    \begin{equation} \rho=\alpha_1\rho_1+(1-\alpha_1)\rho_2 \end{equation}

    通常,对$n$相流系统,基于体积分数平均的密度可以通过以下形式计算:

    \begin{equation} \rho=\sum\alpha_\rq\rho_\rq \end{equation}

    其他的所有物性(比如,粘度)都可以通过这种形式计算。

    17.3.6. 动量方程

    VOF通过在整个计算域中求解一组动量方程得到速度场,且速度场被所有相共享。下式动量方程中的物性参数$\rho$和$\mu$取决于所有相的体积分数。

    \begin{equation} \frac{\p \rho\bfU}{\p t}+\nabla\cdot(\rho\bfU\bfU)=-\nabla p+\nabla\cdot\left(\mu(\nabla\bfU+\nabla\bfU^{\mathrm{T}})\right)+\rho\mathbf{g}+\mathbf{F} \end{equation}

    这种不同相共享一个场的缺点是:如果相间存在较大的速度差异,界面附近的速度计算不准。

    同时需要注意,当粘度比大于$1\times10^3$,可能会导致收敛困难。CICSAM适用于具有较高粘度比的流动,因此可用来解决这种收敛性差的问题。

    17.3.7. 能量方程

    能量方程,同样被所有相共享,如下所示:

    \begin{equation} \frac{\p \rho E}{\p t}+\nabla\cdot(\bfU(\rho E+p))=\nabla\cdot(k_{eff}\nabla T)+S_h \end{equation}

    其中$E$表示能量,$T$表示温度,其均为基于质量平均的变量:

    \begin{equation} E=\frac{\sum_{\rq=1}^N \alpha_\rq\rho_\rq E_\rq}{\sum_{\rq=1}^N \alpha_\rq\rho_\rq} \end{equation}

    式中每一相的$E_\rq$是该相的比热容和温度的函数。

    密度$\rho$和有效导热系数$k_{eff}$被所有相共享,源项$S_h$包括辐射,以及其他的体积热源项。

    与速度场一样,当各相温差较大,将限制界面附近温度的计算准确性。各相物性如果相差几个数量级时也会引起同样的问题。举例,如果计算中包含液态金属和空气,这两种材质的热导率可能会相差$4$个数量级。如此大的物性参数差别会导致离散方程组具有各向异性系数(李东岳:此处可以理解为梯度较大的情况下,网格不同方向的面具有不同的物性系数。若梯度较为均一,在网格足够密的情况下,物性系数可以认为是渐变的,并不会引起各向异性),反过来会导致收敛性和精度问题。

    17.3.8.附加标量传输方程

    在用户的求解过程中可以含有附加标量传输方程,这取决于用户的问题定义。比如在湍流问题中,需附加求解一组湍流传输方程,湍流变量(比如$k$,$\varepsilon$等)被整个计算域内的所有相共享。

    17.3.9. 表面张力和粘附

    VOF模型考虑相界面的表面张力效应。ANSYS Fluent可通过指定相与壁面、相与多孔跃阶的接触角来增强模型的精准度。表面张力系数可设为常数、或为温度的多项式函数、或为任意变量的函数(使用UDF定义)。在控制方程中若可考虑表面张力系数的变化引起的附加切向应力项,则会引起马兰哥尼对流(李东岳:表面张力会由于表面活性剂、溶液浓度产生切向表面张力梯度进而影响流动。日常生活中比较典型的马兰哥尼对流现象即为红酒的挂杯现象,在外太空中的微重力下,马兰哥尼对流现象也比较明显,可参考NASA)。只有在零或接近零重力的场合下,各种表面张力系数效应才有重要影响。

    17.3.9.1.表面张力

    表面张力是流体中分子吸引力作用的结果。比如,以水中的气泡为例,在气泡内部,任一分子受到其周围分子作用力的合力为0。但是,在气泡表面,合力不为0,方向为径向向内。在整个球形表面上的这种径向力共同作用下,表面有收缩的趋势,因此使得表面内凹一侧的压力较大。表面张力是一种仅作用在表面上的力,其使得表面在上述力的相互作用下保持平衡。即表面张力和径向向内的分子间吸引力、径向向外的表面内外压力差构成一个平衡系统。在流体分离时,如果某流体不是球形,表面张力会通过减小表面积使得表面自由能最小化。

    ANSYS Fluent的提供两种表面张力模型:连续表面张力模型(Continuum Surface Force,CSF)、连续表面应力模型(Continuum Surface Stress,CSS)。


    注意:

    使用三角形或四面体网格计算表面张力效应不如四边形或六面体网格精确。
    因此,在表面张力效应具有重要影响的场合,应使用四边形或六面体网格。

    17.3.9.1.1. 连续表面张力(CSF)模型

    CSF模型由Brackbill et al.提出,在VOF计算中,该模型通过动量方程的源项添加表面张力效应。为理解源项的起源,以球形气泡为例,即表面张力为常量,且只考虑界面上的法向力。可以证明表面内外的压力降取决于表面张力系数$\sigma$,以及两个界面法向上的表面曲率半径$R_1$和$R_2$:

    \begin{equation} p_1 - p_2=\sigma\left(\frac{1}{R_1}+\frac{1}{R_2} \right) \end{equation}

    式中,$p_1$和$p_2$为界面两侧的压力(李东岳:如何推导上述方程请参考OpenFOAM中的VOF模型)。

    ANSYS Fluent中,CSF模型中的界面法向$\bfn$由相分数的梯度计算得出:

    \begin{equation} \bfn=\nabla\alpha_\rq \end{equation}

    界面的曲率$k$,定义为界面法向$\bfn$的散度:

    \begin{equation} k=\nabla\cdot\frac{\bfn}{|\bfn|} \end{equation}

    根据散度定理,表面处的力可转换为为体积力。这个体积力就是动量方程中增加的源项,其表达式为:

    \begin{equation} \bfF_\rP=\sum_{i< j} \sigma_{ij}\frac{\alpha_i\rho_i k_j\nabla\alpha_j+\alpha_j\rho_j k_i\nabla\alpha_i}{1/2(\rho_i+\rho_j)} \end{equation}

    对于含有两相或多相的网格单元,上述表达式允许该网格单元附近的力平滑叠加。如果网格单元仅含有两相,其可以化简为:

    \begin{equation} \bfF_\rP=\sigma_{ij}\frac{\rho k_i\nabla\alpha_i}{1/2(\rho_i+\rho_j)} \label{17.23} \end{equation}

    式中,$\rho$为体积平均后的密度。方程$\eqref{17.23}$表明网格单元的表面张力源项与该单元的密度成正比。

    17.3.9.1.2. 连续表面应力(CSS)模型

    与非守恒形式CSF模型不同,另一种模型为守恒的连续表面应力(CSS)模型。CSS模型避免了曲率的显式计算,CSS可以看成是一种基于表面应力来计算表面张力源项的模型。

    CSS模型中,由表面张力引起的表面应力张量$\bfT$可表示为:

    \begin{equation} \bfT=\sigma(\bfI-\hat{\bfn}\otimes\hat{\bfn})|\bfn| \end{equation} \begin{equation} \bfn=\nabla\alpha \end{equation}\begin{equation} \hat{\bfn}=\frac{\bfn}{|\bfn|} \end{equation}

    其中$\bfI$表示单位张量。$\bfT$可以进一步化简为: \begin{equation} \bfT=\sigma\left(|\nabla\alpha|\bfI-\frac{\nabla\alpha\otimes\nabla\alpha}{\nabla\alpha}\right) \end{equation} 随后表面张力可以表示为: \begin{equation} \bfF=\nabla\cdot\bfT \end{equation}

    17.3.9.1.3. CSS和CSF的比较

    CSS模型比CSF模型具有些许优势,特别是在涉及动态表面张力的场合。还需要注意的是CSF和CSS模型都会因为压力梯度和表面张力的不平衡在界面引入一些不符合物理速度。

    CSF模型中,表面张力用以下非守恒形式表示:

    \begin{equation} \bfF_{CSF}=\sigma k\nabla\alpha \end{equation}

    式中$k$为曲率。该表达式仅对恒定表面张力有效。

    对于动态表面张力,CSF模型需要用户在界面切向方向上,引入一个基于表面张力梯度的附加项。

    CSS模型中,表面张力用以下守恒形式表示:

    \begin{equation} \bfF_{CSS}=\nabla\cdot\left(\sigma\left(|\nabla\alpha|\bfI-\frac{\nabla\alpha\otimes\nabla\alpha}{\nabla\alpha}\right)\right) \end{equation}

    CSS模型不需要对曲率进行显式计算。因此,该模型在低分辨率区域(如尖角处)计算更加准确。

    CSS模型因为本身的守恒形式,对于变表面张力的情况不需要引入附加项。

    17.3.9.1.4.什么时候表面张力效应很重要

    判断表面张力效应是否重要,可采用两个无量纲数来判断:雷诺数$\mathrm{Re}$和毛细管数$\mathrm{Ca}$,或者雷诺数$\mathrm{Re}$和韦伯数$\mathrm{We}$。当$\mathrm{Re}\ll 1$的时候,表面张力效应的重要性取决于毛细管数:

    \begin{equation} \mathrm{Ca}=\frac{\mu\bfU}\sigma \end{equation}

    当$\mathrm{Re}\gg 1$,表面张力效应的重要性取决于韦伯数:

    \begin{equation} \mathrm{We}=\frac{\rho L\bfU^2}\sigma \end{equation}

    式中,$\bfU$为自由流速度。如果$\mathrm{Ca}\gg 1$或者$\mathrm{We}\gg 1$,可忽略表面张力效应。

    要在计算中包含表面张力效应,参考《ANSYS Fluent用户指南》中包含表面张力和粘附效应一节。

    17.3.9.2. 壁面粘附

    在VOF模型中也可以指定壁面粘附模型(通过指定壁面接触角)并与表面张力模型结合使用。该模型是Brackbill et al.提出的。该模型并不是把该边界条件直接用在壁面本身,而是通过指定流体和壁面之间界面法向来指定接触角。这种动态的边界条件导致壁面附近界面曲率发生变化。

    如果壁面处的接触角记为$\theta_w$,那么壁面附近毗连网格单元的表面法向为:

    \begin{equation} \hat{\bfn}=\hat{\bfn}_w\mathrm{cos}\theta_w+\hat{\mathbf{t}}_w\mathrm{sin}\theta_w \end{equation}

    式中$\hat{\bfn}_w$,$\hat{\mathbf{t}}_w$分别为壁面的单位法向量和切向量。指定的接触角和计算的壁面毗连的网格单元的表面法向共同确定界面的局部曲率,该曲率被用来调整表面张力计算中的体积力。

    要在计算中包含壁面粘附效应,请参考《ANSYS Fluent用户指南》中包含表面张力和粘附效应一节。

    17.3.9.3. 跳跃粘附

    VOF模型也提供了一个和壁面粘附模型类似的跳跃粘附模型。跳跃粘附模型中,如果多孔跃阶边界两侧的接触角相同,那么这个接触角对任一侧都适用。

    因此,如果多孔跃阶处的接触角记为$\theta_w$,那么壁面附近毗连网格单元的表面法向为:

    \begin{equation} \hat{\bfn}=\hat{\bfn}_w\mathrm{cos}\theta_w+\hat{\mathbf{t}}_w\mathrm{sin}\theta_w \end{equation}

    式中$\hat{\bfn}_w$,$\hat{\mathbf{t}}_w$分别为壁面的单位法向量和切向量。

    要在用户模型中包含跳跃粘附,参考《ANSYS Fluent用户指南》中的设置边界条件一节。

    在多孔跃阶边界处,ANSYS Fluent提供了两种跳跃粘附处理方法:

    17.3.10. 明渠流

    ANSYS Fluent可通过VOF和明渠流边界条件模拟明渠流效应(比如,河流、大坝等自由表面流动)。这些流动的共性是在流动的流体和流体的上方(通常是大气)之间存在自由表面。在这种情况,波浪的推进和自由表面的行为很重要。流动一般由重力和惯性力控制。该模型最适合于海洋领域应用以及排水系统的流动分析。

    明渠流的特征用无量纲弗劳德数来描述,弗劳德数定义为惯性力与静水压力之比:

    \begin{equation} \mathrm{Fr}=\frac{|\bfU|}{\sqrt{|\bfg| y}} \end{equation}

    式中$\bfU$为速度,$\bfg$为重力加速度,$y$为特征长度,这里指的是渠道底部和自由表面之间的距离。$\sqrt{|\bfg| y}$为是波浪的传播速度。对于静止的观测者来说,波速定义为:

    \begin{equation} U_w=|\bfU|+\sqrt{|\bfg| y} \end{equation}

    基于弗劳德数,明渠流又可分为以下三类:

    17.3.10.1.上游边界条件

    明渠流的上游边界条件有两种:

    17.3.10.1.1. 压力进口

    进口处的总压$p_0$可给定为:

    \begin{equation} p_0=\frac{1}{2}\rho\bfU^2+(\rho-\rho_0)|\bfg|\left(\frac{\bfg}{|\bfg|}\cdot(\mathbf{b}-\mathbf{a})\right) \end{equation}

    式中$\mathbf{b}$,$\mathbf{a}$分别为网格单元面心的位置向量和自由表面在边界处的位置向量,在这里假定自由表面是水平的且垂直于重力方向。$\bfg$为重力加速度向量,$|\bfg|$为重力加速度的大小,$|\bfU|$为速度的大小,$\rho$为网格单元内流体的密度,$\rho_0$为参考密度。

    另外,动压力$p_d$可以表示为:

    \begin{equation} p_d=\frac{\rho}{2}\bfU^2 \end{equation}

    静压力$p_s$可以表示为:

    \begin{equation} p_s=(\rho-\rho_0)|\bfg|\left(\frac{\bfg}{|\bfg|}\cdot(\mathbf{b}-\mathbf{a})\right) \label{ps} \end{equation}

    其也可以展开为:

    \begin{equation} p_s=(\rho-\rho_0)|\bfg|\left(\frac{\bfg}{|\bfg|}\cdot\mathbf{b}+y_{local} \right) \end{equation}

    其中自由表面和参考位置之间的距离$y_{local} $为:

    \begin{equation} y_{local}=-\mathbf{a}\cdot\frac{\bfg}{|\bfg|} \label{ylocal} \end{equation}

    17.3.10.1.2. 质量流量

    明渠流中各相的质量流量定义为:

    \begin{equation} m_\rq=\rho_\rq\bfS_\rq\cdot\bfU_\rq \end{equation}

    17.3.10.1.3. 有关体积分数的说明

    在明渠流中,ANSYS Fluent根据边界条件对话框中指定的参数自动计算体积分数,因此体积分数选项没有被激活。

    对次临界进口流动($\mathrm{Fr}<1$),ANSYS Fluent根据边界附近网格单元的体积分数来重构边界处的体积分数:

    对于超临界进口流动($\mathrm{Fr}>1$),体积分数可根据自由表面到底部的固定高度来计算。

    17.3.10.2.下游边界条件

    17.3.10.2.1. 压力出口

    静压取决于压力设定方法。如果选择自由表面基准法,静压由方程$\eqref{ps}$和方程$\eqref{ylocal}$计算确定,否则用户必须指定静压力为表压。

    对于次临界出口流动($\mathrm{Fr}<1$),如果是两相流,压力要么是由边界指定的压力分布确定,要么是由毗连网格单元计算确定。对于超临界流动($\mathrm{Fr}>1$),压力通常是由毗连网格单元计算确定。

    17.3.10.2.2. 出流边界

    如果在求解流动问题前并不知道速度和压力等信息,可以在明渠流的出口使用出流边界。如果出流边界条件未知,ANSYS Fluent将从内部插值获得相关信息。

    然而,该边界条件具有以下局限性:

    17.3.10.2.3. 设定回流体积分数

    ANSYS Fluent内部通过毗连的网格单元值计算出口边界的体积分数,因此,该选项没有被激活。

    17.3.10.3.数值岸

    某些场合下需要抑制波浪在出口边界传播引起的数值反射。为避免波浪反射,可在压力出口边界附近单元区域的动量方程中引入一个阻尼汇项,可参考下面的文章:Park et al.Zwart et al.

    \begin{equation} S=-\left[C_1\rho V+\frac{1}{2}C_2\rho|V|V\right]f(z)f(x) \end{equation}

    式中:

    $\hat{z}$为垂直方向( 即重力方向);

    $\hat{x}$为流动方向;

    $S$为$\hat{z}$方向上的动量源项;

    $C_1$为线性阻尼阻力($1/s$)(默认值为10);

    $C_2$为二次阻尼阻力($1/m$)(默认值为10);

    $V$为沿$\hat{z}$方向的速度;

    $z$为距自由表面的距离;

    $x$为沿$\hat{x}$方向的距离;

    $f(x)$为$\hat{x}$方向上的阻尼函数;

    $f(z)$为$\hat{z}$方向上的阻尼函数;

    $\hat{x}$方向和$\hat{z}$方向的比例因子$r_x$以及$r_z$分别由下面的公式定义:

    \begin{equation} r_x=\frac{x-x_s}{x_e-x_s} \end{equation} \begin{equation} r_z=\frac{z-z_{fs}}{z_b-z_{fs}} \end{equation}

    $\hat{x}$方向和$\hat{z}$方向的阻尼函数分别如下式所示:

    \begin{equation} f(x)=(r_x)^2 \end{equation} \begin{equation} f(z)=1-r_z \end{equation}

    其中$x_s$和$x_e$分别为$\hat{x}$方向阻尼区域的起点和终点,$z_{fs}$和$z_b$分别为沿 $\hat{z}$方向上的自由表面和底部高度。


    注意:

    ANSYS Fluent的数值岸处理中,在$\hat{z}$方向采用线性阻尼,在$\hat{x}$方向采用二次阻尼。
    且只有明渠流和明渠造波边界条件中才有该选项。

    要在模拟中使用数值岸,详见明渠流的数值岸一节。

    17.3.11. 明渠造波边界条件

    明渠造波边界条件可以模拟规则波和不规则波的传播,该模型在海洋领域十分有用,比如可以用来分析波浪运动学以及波浪对海上结构体或运动物体的冲击作用。

    波陡一般定义为波高与波长之比,相对深度定义为波高与液体深度之比。小幅振波理论一般适用于较低的波陡和较低的相对深度,而有限幅振波理论更适用于大波陡或者较大的相对深度。

    ANSYS Fluent通过速度入口边界条件为表面重力波提供以下几种选项:

    每个波理论中的短重力波的表达式都是基于无限液体深度的,而浅水波或中等水深波的表达式都是基于有限液体深度的。

    波高$H$定义如下:

    \begin{equation} H=2A=A_t+A_c \end{equation}

    式中$A$为波幅,$A_t$为波谷处幅值,$A_c$为波峰处幅值。对于线性波理论,$A_c= A_t$;对于非线性波理论,$A_c\neq A_t$。

    波数$k$定义为:

    \begin{equation} k=\frac{2\pi}{\lambda} \end{equation}

    式中$\lambda$为波长。

    波数$k$的矢量形式为$\bfK$:

    \begin{equation} \bfK=k_x\hat{x}+k_y\hat{y} \end{equation}

    式中$\hat{x}$为波传播的参考方向,在这里定义$\hat{z}$为重力的反方向,则$\hat{y}$为垂直$\hat{x}$和$\hat{z}$的方向。$\hat{x}$方向和$\hat{y}$方向上的波数分别通过下式定义: \begin{equation} k_x=k \mathrm{cos}\theta \end{equation} \begin{equation} k_y=k \mathrm{sin}\theta \end{equation}

    式子中$\theta$为波向角,其定义为在$\hat{x}$和$\hat{y}$的平面上波面和参考波传播方向的夹角。

    有效波频$\omega_e$定义为:

    \begin{equation} \omega_e=\omega+\bfK\bfU \end{equation}

    式中$\omega$为波频,$\bfU$为流动的平均速度。当存在运动的物体且流动和运动物体的运动相关时,物体运动的影响可被包含在$\bfU$中。

     

    波速$c$定义为:

    \begin{equation} c=\frac{\omega}{k} \end{equation}

    通过叠加所有波的速度分量获得进口的速度矢量可表示为:

    \begin{equation} \bfV=\bfU+u\hat{x}+v\hat{y}+w\hat{z} \end{equation}

    其中$u$,$v$以及$z$分别为表面重力波在$\hat{x}$,$\hat{y}$和$\hat{z}$方向的速度分量。

    17.3.11.1. Airy 波理论

    某一位置线性波的波剖面为:

    \begin{equation} \eta(t)=A\mathrm{cos}\alpha \end{equation}

    式中$\alpha$定义如下:

    \begin{equation} \alpha=k_xx+k_yy-\omega_et+\varepsilon \end{equation}

    其中$x$,$y$分别代表$\hat{x}$和$\hat{y}$方向的空间坐标,$\varepsilon$为相位差,$t$ 为时间。

    对浅水波和中等深水波,波频$\omega$定义如下

    \begin{equation} \omega=\sqrt{\bfg k\mathrm{tanh}(kh)} \end{equation}

    对短重力波,波频$\omega$定义如下

    \begin{equation} \omega=\sqrt{\bfg k} \end{equation}

    其中$h$为波高,$k$为波数,$\bfg$为重力加速度。

    依据浅水波、中等深水波,还是短重力波的分类,入射波边界条件处的速度分量可相应定义如下:

  • 对浅水波和中等深水波,其速度分量:
  • \begin{equation} \begin{bmatrix} u\\ v \end{bmatrix}=\frac{\bfg k A}{\omega}\frac{\mathrm{cosh}\left[k(z+h)\right]}{\mathrm{cosh}(kh)} \begin{bmatrix} \mathrm{cos}\theta\\ \mathrm{sin}\theta \end{bmatrix}\mathrm{cos}\alpha \end{equation} \begin{equation} w=\frac{\bfg k A}{\omega}\frac{\mathrm{sinh}\left[k(z+h)\right]}{\mathrm{cosh}(kh)}\mathrm{sin}\alpha \end{equation}
  • 对Short Gravity Waves(李东岳:意义不详),其速度分量:
  • \begin{equation} \begin{bmatrix} u\\ v \end{bmatrix}=\frac{\bfg k A}{\omega}e^{kz} \begin{bmatrix} \mathrm{cos}\theta\\ \mathrm{sin}\theta \end{bmatrix}\mathrm{cos}\alpha \end{equation} \begin{equation} w=\frac{\bfg k A}{\omega}e^{kz}\mathrm{sin}\alpha \end{equation}

    其中$z$为自由表面水准在$\hat{z}$方向(即重力反方向)上的高度。更多关于如何使用和设置该模型的信息,详见《ANSYS Fluent用户指南》中的模拟明渠造波边界条件一节。


    待更新