twoLiquidMixingFoam解析


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

1. 引言

twoLiquidMixingFoam可用于计算两相混合问题。两相混合问题和VOF模型非常相似。区别之处仅仅在于两相混合问题1)需要考虑相间的分子扩散,因为两种物质是互溶的。2)不需要考虑表面张力效应,因为并不存在界面力。同时,因此,twoLiquidMixingFoam也不需要保证相界面的尖锐,不需要进行界面压缩。在阅读本文内容之前,可参考interFoam解析了解更基本的描述。

2. 控制方程与算法

首先,VOF中的动量方程为

\begin{equation}\label{vof1} \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{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 \end{equation}
其即为twoLiquidMixingFoam中的动量方程。同样的,VOF中的相方程为
\begin{equation}\label{alpha} \frac{\partial \alpha}{\partial t}+\nabla\cdot(\bfU\alpha)=0 \end{equation}
在其中添加相间扩散项有:
\begin{equation}\label{alpha2} \frac{\partial \alpha}{\partial t}+\nabla\cdot(\bfU\alpha)=\nabla\cdot\left(\left(D_c+D_t\right)\nabla\alpha\right) \end{equation}
其中$D_c$表示层流扩散系数,其需要用户自行给定。$D_t$表示湍流扩散系数,湍流扩散系数可表示为
\begin{equation}\label{Dt} D_t=\frac{\nu_t}{\mathrm{St}} \end{equation}
其中$\nu_t$表示湍流粘度,$\mathrm{St}$表示施密特数。连续性方程可以写为:
\begin{equation}\label{mass} \frac{\p\rho}{\p t}+\nabla\cdot(\rho\bfU)=0 \end{equation}
方程\eqref{mass},\eqref{alpha2}以及\eqref{vof2}构成twoLiquidMixingFoam的基本控制方程。twoLiquidMixingFoam中速度和压力的耦合可参考interFoam解析。在这里讨论一下相方程\eqref{alpha2}的数值离散。在离散相方程的时候,采用了操作符分裂的方法,这种方法可以更容易的调用高分辨率格式。具体的,首先采用高分辨率格式求解对流项,然后再求解扩散项:
\begin{equation}\label{split1} \frac{\partial \alpha}{\partial t}+\nabla\cdot(\bfU\alpha)=0 \end{equation} \begin{equation}\label{split2} \frac{\partial \alpha}{\partial t}-\frac{\partial \alpha}{\partial t}-\nabla\cdot\left(\left(D_c+D_t\right)\nabla\alpha\right)=0 \end{equation}
方程\eqref{split2}中$\frac{\partial \alpha}{\partial t}$相互抵消,第一项隐性离散增加对角占优,第二项显性离散互相抵消。通过这种方式,方程\eqref{split1}可以调用显性推进的高分辨率格式。在OpenFOAM继续通过反扩散算法保证相分数的有界。方程\eqref{split2}具备椭圆的特征采用更加耦合的隐性离散保证稳定性。

3. 代码分析
for
    (
        subCycle alphaSubCycle(alpha1, nAlphaSubCycles);
        !(++alphaSubCycle).end();
    )
    {
        #include "alphaEqn.H"\\方程\eqref{split1}
        rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
    }
	
	fvScalarMatrix alpha1Eqn\\方程\eqref{split2}
    (
        fvm::ddt(alpha1)\\隐性离散
      - fvc::ddt(alpha1)\\显性离散
      - fvm::laplacian
        (
            volScalarField("Dab", Dab + alphatab*turbulence->nut()),
            alpha1
        )
    );

    alpha1Eqn.solve();

东岳流体®版权所有
勘误、讨论、补充内容请前往CFD中文网