interFoam解析


如果打不开图像,请右键在新标签页打开图像后刷新几次 如果打不开图像,请右键在新标签页打开图像后刷新几次

1. 引言

在有重力以及源项的情况下,有不可溶多相体系动量方程:

\begin{equation}\label{mom} \frac{\partial \rho \mathbf{U}}{\partial t}+\nabla \cdot (\rho \mathbf{U} \mathbf{U} ) - \nabla \cdot \tau =-\nabla p+\rho\bfg+\bfF \end{equation}

其中各项的含义可参考icoFoam解析,$\bfF$表示表面张力。连续性方程可以表示为

\begin{equation}\label{cont} \frac{\partial \rho}{\partial t}+\nabla \cdot (\rho \mathbf{U})=0 \end{equation}

若考虑两种流体均为不可压缩的,也即考虑某种流体运动的流体微元,其物质导数为$0$,即

\begin{equation}\label{D} \frac{\rD \rho}{\rD t}=\frac{\partial \rho}{\partial t}+\bfU\cdot\nabla\rho=0 \end{equation}

将方程\eqref{D}代入到\eqref{cont}有:

\begin{equation}\label{cont2} \nabla\cdot\bfU=0 \end{equation}

在这种情况下,方程\eqref{mom}和\eqref{cont2}组成了不可压缩多相界面类模型中关于速度和压力的方程。其中不同的相具有不同的密度。在速度$\bfU$已知的情况下,方程\eqref{D}中的密度为纯对流方程可以进行求解。在Front-tracking方法中,方程\eqref{D}并没有直接进行求解而是通过某种界面重组的方法进行组建。在Volume of fluid(VOF)方法中则定义了一个变量$\alpha$来表示流体的相分数。考虑某一个网格单元的气液两相系统,如果此网格单元内充满了流体,则$\alpha=1$;如果此网格单元内充满了气体,则则$\alpha=0$。如果$\alpha$的值介于$0$和$1$之间,则此网格单元内为气液混合。可见,CFD中VOF的界面是通过跟踪相分数来获得的,并不是直接算出来的,是通过计算每个网格内的相场值而后处理出来的。网格越密,VOF后处理出来的界面越薄。VOF与Front-tracking这些模型均为界面重构的微观多相流模型。可以看作为多相计算流体力学领域的直接模拟。OpenFOAM中的VOF模型求解器,interFoam,可用于求解多相流的界面问题,如崩塌的堤坝、极少量的气泡界面、液体晃荡等。界面类存在两个关键问题:1)界面尽可能要薄,2)在界面失稳的情况下算法要尽可能的稳定。本文从最基本的液膜压力降开始分析,一步一步推导VOF模型。对于大量气泡/液滴存在的情况,VOF模型对网格分辨率要求较高。在这种情况下可以使用多流体模型来计算。


Front-tracking方法(左)与VOF方法(右)示意图。
在front-tracking中,通过跟踪粒子的位置获得界面形状。在VOF中,界面由一系列相分数不为0和1的网格构成。

2. 控制方程与算法
2.1 表面张力与VOF模型

表面张力最重要的特征即其会导致界面处存在一个界面压力降$\Delta p$。这个效应需要在动量方程中有所体现。考虑下图的小界面微元,定义表面张力为$\bfF_\sigma$,其是一种为了保持界面平衡而作用在每单位长度上的力,是一种和界面相切的张力。表面张力的大小主要和流体的物理特性有关。在曲面中,均衡的流体表面张力和压力降$\Delta p$均衡。这个界面压力降主要取决于表面张力和曲面弧度。


图1. 两相界面处的曲面微元受力示意图

如图(1)所示的一个曲面微元。其中$p_1$为曲面微元上方液体向下施加的压强,$p_2$为曲面微元下方液体向上施加的压强,$s_1$和$s_2$分别表示曲面微元的俩个边长。曲面微元上的总压力的值为:

\begin{equation}\label{Fp} F_p=(p_1 - p_2) \mathrm{d}s_1 \mathrm{d}s_2 \end{equation}

基于表面张力的定义,其作用在曲面微元的四个边上,且合力向上和曲面微元上的总压力$F_p$均衡。如果表面张力系数定义为$\sigma$。那么作用在$\mathrm{d}s_1$上的表面张力为$\sigma\mathrm{d}s_1$,其在$y$方向上的分量为:$\sigma\mathrm{d}s_1 \mathrm{sin} \frac{\mathrm{d}\alpha}{2} \approx \sigma\mathrm{d}s_1 \frac{\mathrm{d}\alpha}{2} $。以此类推,我们对四个边加和之后,有作用在曲面微元上的表面张力之和的大小为:

\begin{equation}\label{Fsigma} F_\sigma=\sigma\left(\frac{1}{R_1}+\frac{1}{R_2} \right) \mathrm{d}s_1 \mathrm{d}s_2 \end{equation}

其中$R_1$,$R_2$为$\mathrm{d}s_1$,$\mathrm{d}s_2$的主曲率半径。均衡状态下有$F_p = F_\sigma$,即:

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

我们用$\bfn$来表示曲面微元的单位法向,有:

\begin{equation}\label{pn} \frac{\p \bfn}{\p x} \mathrm{d}x= \mathrm{mag}\left(\bfn_x|_{x+\mathrm{d}x} \right) - \mathrm{mag}\left(\bfn_x|_x\right) =\mathrm{sind}\alpha \approx \mathrm{d}\alpha \end{equation}
也即:
\begin{equation} \frac{\partial \bfn}{\partial x} = \frac{\mathrm{d}\alpha}{\mathrm{d} x} = \frac{1}{R_1} \end{equation}
\begin{equation}\label{six} \nabla \cdot \bfn = \frac{\partial \bfn}{\partial x} + \frac{\partial \bfn}{\partial y} + \frac{\partial \bfn}{\partial z} =\left( \frac{1}{R_1}+ \frac{1}{R_2} \right) \end{equation}

将方程\eqref{six}代入到方程\eqref{three}有Young-laplacian方程:

\begin{equation} \Delta p =p_1 - p_2 = \sigma \nabla \cdot \bfn \end{equation}

如图(2)所示,对于CFD中的自由表面,某一点的法向量为$\nabla\alpha$,其为一个连续函数,在流体1和2区内部为0,仅仅在自由表面存在非零值。图(2)自由表面的单位法向量可表示为:

\begin{equation} \bfn=\frac{\nabla \alpha}{|\nabla \alpha|} \end{equation}

图2. 流体界面示意图

进一步,界面处的曲率可定义为:

\begin{equation} \kappa = - \nabla \cdot \bfn \end{equation}

为了将压力降$\Delta p $包含在动量方程中,需要将其转变为梯度函数。在转变的过程中存在一些数值问题。例如,压力降仅仅在界面附近存在。但动量方程却适用于整个计算域。Continuum Surface Force(CSF)模型可以处理这个问题。如图(3)所示,CSF模型将界面内的压力表示为压力的连续函数:

\begin{equation}\label{ten} p_t=\left\{\begin{matrix} p_b & c=1\\ p_a+\sigma c \nabla \cdot \bfn & 0< c< 1\\ p_a & c=0 \end{matrix}\right. \end{equation}

其中$c$为界面处的位置函数,其可表示为$c=\alpha_\mathrm{t} - \alpha_\mathrm{a}$($t$与$a$分别表示位置)。方程\eqref{ten}将界面处间断的压力函数表示为一个连续函数,这为CSF模型的模型基础。


图3. CSF模型下连续的压力函数以及$c$的定义

考虑整个区域内,方程\eqref{ten}可以表示为:

\begin{equation} \bfF_\sigma=\sigma \kappa \nabla \alpha \end{equation}

其中$\bfF_\sigma$仅在界面处有值。同时,方程\eqref{mom}还需要做进一步的变化,以使得边界条件的定义更加简单(在双流体模型中也可以进行类似的数学操作)。首先定义

\begin{equation}\label{prgh} p_\mathrm{rgh}=p-\rho \bfg \cdot\bfh \end{equation}

其中$\bfh$表示网格单元体心的位置矢量。对方程\eqref{prgh}进行梯度操作:

\begin{equation}\label{nablaprgh} \nabla p_\mathrm{rgh}=\nabla p-\bfg\cdot \bfh \nabla \rho - \rho \bfg \end{equation}

将其代入到方程\eqref{mom}有VOF模型的动量方程:

\begin{equation}\label{momF} \frac{\partial \rho \mathbf{U}}{\partial t}+\nabla \cdot (\rho \mathbf{U} \mathbf{U} ) - \nabla \cdot \tau =-\nabla p_\mathrm{rgh}-\bfg\cdot \bfh \nabla \rho+\sigma \kappa \nabla \alpha \end{equation}

下面来推导VOF中的相方程。方程\eqref{mom}中的密度$\rho$可以表示为:

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

即:

\begin{equation}\label{alpha} \frac{\mathrm{D}\rho}{\mathrm{D}t} = \frac{\mathrm{D} \left(\alpha\left(\rho_1 - \rho_2 \right)+\rho_2 \right)}{\mathrm{D}t}= \frac{\mathrm{D}\alpha}{\mathrm{D}t}=\frac{\partial \alpha}{\partial t}+\bfU\cdot\nabla\alpha=0 \end{equation}

方程\eqref{alpha}即为不可压缩VOF模型中的相方程。

2.2. 离散
综合考虑,不可压缩VOF模型中的方程为:

\begin{equation}\label{vof1} \frac{\partial \alpha}{\partial t}+\bfU\cdot\nabla\alpha=0 \end{equation}
\begin{equation}\label{vof2} \frac{\partial \rho \mathbf{U}}{\partial t}+\nabla \cdot (\rho \mathbf{U} \mathbf{U} ) - \nabla \cdot \tau =-\nabla p_\mathrm{rgh}-\bfg\cdot \bfh \nabla \rho+\sigma \kappa \nabla \alpha \end{equation}
\begin{equation}\label{vof3} \nabla\cdot\bfU=0 \end{equation}

首先对方程\eqref{vof2},省略方程右侧的项,进行有限体积离散有(更详细的步骤可参考icoFoam解析):

\begin{equation} {A_\mathrm{P}}\mathbf{U}_\mathrm{P}^r{\rm{ + }}\sum {A_\mathrm{N}\mathbf{U}_\mathrm{N}^r} = S_\mathrm{P}^n, \label{apanmom} \end{equation}

求解后有预测速度$\bfU_{\rP}^r$。预测速度并不符合连续性方程,考虑最终收敛的情况,对连续性方程进行离散有:

\begin{equation} \sum \left(\mathbf{U}_{\mathrm{P},f}^{n+1} \cdot \bfS_f\right)=0. \label{contfv} \end{equation}

如果能用压力表示方程\eqref{contfv}中的$\mathbf{U}_{\mathrm{P},f}^{n+1}$,则压力泊松方程即可构建。首先,收敛情况下方程\eqref{apanmom}可以写为

\begin{equation} {A_\mathrm{P}}\mathbf{U}_\mathrm{P}^{n+1}{\rm{ + }}\sum {A_\mathrm{N}\mathbf{U}_\mathrm{N}^{n+1}} = S_\mathrm{P}^n -\nabla p_{\mathrm{rgh},\rP}-\bfg\cdot \bfh \nabla \rho_\rP+\sigma \kappa \nabla \alpha_\rP, \label{apanmom2} \end{equation}

定义

\begin{equation} \mathbf{HbyA}^{n+1}_\mathrm{P} = \frac{1}{{{A_\mathrm{P}}}}\left( { - \sum {{A_\mathrm{N}}\mathbf{U}_\mathrm{N}^{n+1}} + S_\mathrm{P}^n} \right). \label{hbya} \end{equation}

\begin{equation} \mathbf{U}_\mathrm{P}^{n+1} = \mathbf{HbyA}^{n+1}_\mathrm{P} -\frac{1}{A_\rP}\left(\nabla p_{\mathrm{rgh},\rP}^{n+1}+\bfg\cdot \bfh \nabla \rho_\rP^{n+1}-\sigma \kappa \nabla \alpha_\rP^{n+1}\right), \label{up} \end{equation}

依据Rhie-Chow插值,面上的速度为

\begin{equation} \mathbf{U}_\mathrm{P,f}^{n+1} = \bfHbyA_{\rP,f}^{n+1} -\frac{1}{A_{\rP,f}}\left(\nabla_f^\bot p_{\mathrm{rgh},\rP}^{n+1}+\bfg\cdot \bfh \nabla_f^\bot \rho_\rP^{n+1}-\sigma \kappa \nabla_f^\bot \alpha_\rP^{n+1}\right), \label{upf} \end{equation}

将方程\eqref{upf}代入到\eqref{contfv}有

\begin{equation} \sum\left( \bfHbyA_{\rP,f}^{n+1}-\frac{1}{A_{\rP,f}}\left(\bfg\cdot \bfh \nabla_f^\bot \rho_\rP^{n+1}+\sigma \kappa \nabla_f^\bot \alpha_\rP^{n+1} \right)\right)=\sum \frac{1}{A_{\rP,f}}\left(\nabla_f^\bot p_{\mathrm{rgh},\rP}^{n+1}\right). \label{contfv2} \end{equation}

\begin{equation} \nabla\cdot\left(\frac{1}{A}\nabla p_{\mathrm{rgh}}^{n+1}\right)=\nabla\cdot\left(\bfHbyA^{n+1}-\frac{1}{A}\left(\bfg\cdot \bfh \nabla\rho^{n+1}+\sigma \kappa \nabla \alpha^{n+1}\right)\right) . \label{po} \end{equation}

求解方程\eqref{po}既有收敛的压力。

2.3. 界面尖锐

正如在引言中所讨论的,在VOF模型中,相分数作为一个特殊的被动标量,具有明显的数值特性。首先是相分数$\alpha$在数值上应该保证在$0-1$之间,任何越界的相分数都不符合物理。其次,如果考虑网格非常细密的情况下,相分数应该只存在两个数值:0或者1。这种非常致密的网格会导致异常高的计算资源需求,因此VOF模型中不可避免的会存在$\alpha$在$0-1$之间的情况。但VOF应尽可能的使相界面足够的尖锐。在动网格算法中,界面的尖锐可以依附于网格的变化。在Front-tracking方法中,界面尖锐可以通过示踪颗粒来获得。历史上存在不同的可用于VOF模型中保证界面尖锐的方法。在OpenFOAM中采用的是Weller提出的方法,其通过一种人工的对流项来对相界面附近的相分数进行挤压来抗衡这种数值耗散带来的相界面模糊性,并且这一人工对流项在数值上需要保证非相界面处为零。依据添加人工对流项的思想,VOF模型可以表示为

\begin{equation}\label{compress} \frac{\partial \alpha}{\partial t}+\nabla\cdot\left(\alpha\bfU\right)+\nabla\cdot\left(\alpha(1-\alpha)\bfU_c\right)=0 \end{equation}

方程\eqref{compress}中的第三项为人工添加的可压缩项,其在纯相(非界面处)计算域为$0$。$\bfU_c$为需要模化的速度。其应该在界面的法向上进行压缩而不是切向,否则引起虚假的扩散。因此$\bfU_c$的方向应与$\bfn$同向。因此有

\begin{equation}\label{compress2} \bfU_c=f\left(\frac{\nabla\alpha}{|\nabla\alpha|}\right) \end{equation}

另一个问题是压缩速度的大小问题。很明显压缩速度不能过分大,这并不符合物理。因此压缩速度最大值也只不过是$\bfU$,则有:

\begin{equation}\label{compress3} \bfU_c=c|\bfU|\frac{\nabla\alpha}{|\nabla\alpha|} \end{equation}

其中的$c$表示可控的压缩因子。当$c=0$,无压缩效果。$c$越大,压缩效应越快也越明显。最终有相方程:

\begin{equation}\label{compress4} \frac{\partial \alpha}{\partial t}+\nabla\cdot\left(\alpha\bfU\right)+\nabla\cdot\left(\alpha(1-\alpha)c|\bfU|\frac{\nabla\alpha}{|\nabla\alpha|}\right)=0 \end{equation}
3. 代码分析

fvVectorMatrix UEqn
(
	fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
  + MRF.DDt(rho, U)
  + turbulence->divDevRhoReff(rho, U)
 ==
	fvOptions(rho, U)
);

UEqn.relax();

fvOptions.constrain(UEqn);

if (pimple.momentumPredictor())
{
	solve
	(
		UEqn
	 ==
		fvc::reconstruct
		(
			(
				mixture.surfaceTensionForce()
			  - ghf*fvc::snGrad(rho)
			  - fvc::snGrad(p_rgh)
			) * mesh.magSf()
		)
	);

	fvOptions.correct(U);
}
//方程\eqref{vof2}

volVectorField HbyA(constrainHbyA(rAU()*UEqn.H(), U, p_rgh));
surfaceScalarField phiHbyA
(
	"phiHbyA",
	fvc::flux(HbyA)
  + MRF.zeroFilter(fvc::interpolate(rho*rAU())*fvc::ddtCorr(U, phi, Uf))
);

surfaceScalarField phig
(
	(
		mixture.surfaceTensionForce()
	  - ghf*fvc::snGrad(rho)
	)*rAUf*mesh.magSf()
);


phiHbyA += phig;//方程\eqref{upf}

fvScalarMatrix p_rghEqn
(
	fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
);
//方程\eqref{po}

surfaceScalarField phic(mixture.cAlpha()*mag(phi/mesh.magSf()));
surfaceScalarField phir(phic*mixture.nHatf());
//方程\eqref{compress3}

tmp talphaPhi1Un
(
	fvc::flux
	(
		phiCN(),
		cnCoeff*alpha1 + (1.0 - cnCoeff)*alpha1.oldTime(),
		alphaScheme
	)
  + fvc::flux
	(
	   -fvc::flux(-phir, alpha2, alpharScheme),
		alpha1,
		alpharScheme
	)
);

MULES::explicitSolve
(
	geometricOneField(),
	alpha1,
	phiCN,
	alphaPhi10,
	Sp,
	(Su + divU*min(alpha1(), scalar(1)))(),
	oneField(),
	zeroField()
);
//方程\eqref{compress4}

4. 验证算例

interFoam自带若干算例。在OpenFOAM-6中,自带的算例主要为以下几个:

除上述算例外,自带算例还包含若干其他算例在此不一一介绍。本文介绍的验证算例为管道段塞流/泰勒泡。细长管道内的气液流动通常可分为5种流形:泡状流(bubbly flow)、段塞流(slug flow)、搅拌流(churn flow)、塞状环状流(slug-annular flow)和环装流(annular flow)。泡状流指气泡以不规则的类球形形态分散在液相中,气相可以看作为离散相,液相可以看作为连续相。在气液相表观速度相差不大或基本相等的情况下,会形成段塞流,又称泰勒流。段塞流的特点为气塞与气塞互相被各自分开,气泡与管道壁面之间有一层很薄的液膜,气塞单元或者液塞单元的长度大于或者等于管道横截面的宽度。在通过VOF模型进行计算的时候,模型需要能够捕获这种段塞的行为而不是分层流。

本验证算例来源于CFD中文网。用户红豆沙提供的网格文件,东岳流体进行的算例设置。


泰勒泡的几何与网格

4.1. 模型设置

算例采用的网格为采用ANSYS ICEM生成的2D纯六面体结构网格。气相从下方进入,液体从上方进入,算例运行一段时间后,应出现段塞流特征。0文件夹的主要设置如下:

constant文件夹的物性主要设置如下:

4.2. 算例运行

相关算例运行首先需要执行网格转换命令:

fluentMeshToFoam fluent.msh

网格生成完毕,后直接运行VOF求解器运行即可:

interFoam

本算例运行几秒钟即可结束,0.7613秒的计算结果如下图所示:


4.3. 点评

下面几个方向的透彻研究,可以用来发表SCI论文。

  1. 经测试,不同的气相流速、液相流速展现非常不同的流动行为,用户可尝试调节不同的流速,研究流速对流形的影响;
  2. 本算例采用水表示连续相,用户可尝试使用非牛顿流体作为连续相,研究连续相物性对流形的影响,类似的研究可参考Sontti, S.G., Atta, A., 2017. CEJ, 330, 245-261;
  3. interFoam采用界面压缩技术保证界面的尖锐,用户可尝试调节system文件夹下的fvSolution文件中的cAlpha关键词的值,研究不同界面压缩因子对界面尖锐的影响;同时,VOF作为网格依赖类求解器,加密网格研究网格无关性至关重要;
  4. 接触角对气泡的形状有着重要的影响,用户可研究不同的接触角对液体润湿、段塞长度的影响;
  5. 在OpenFOAM-1812中,一个新的界面模型被引入,用户可采用OpenFOAM-1812中的interFoam对比官方版的interFoam,研究界面尖锐格式对界面的影响;
  6. 有关速度U和压力p_rgh的边界条件,感兴趣的用户可以调用下更简单的zeroGradient(速度)以及fixedValue(压力),可以加深对边界条件的理解;
点击下载

2019.04.26增加验证算例 | 2019.04.14大修 | 2018.05.12张童伟:修正方程\eqref{ten}中的$0 < c < 0$为$0< c< 1$ | 2018.04.20齐雪宇:修正方程\eqref{prgh}、\eqref{nablaprgh}中的$\bfg \bfh$为$\bfg\cdot\bfh$

东岳流体 2014-2019
勘误、讨论、补充内容请前往CFD中文网