用于预测雾化过程的大涡模拟方法

摘要

本发明公开了一种用于预测雾化过程的大涡模拟方法，基于真实速度场构造虚拟的液相速度场，并将构造的液相速度场应用于流动控制方程和界面输运方程的求解，模拟气液两相流的实时动态过程，以准确预测液滴破碎过程和液体射流的雾化过程。本发明提高了两相流模拟的计算精度和稳定性，可以准确计算和预测液滴的破碎过程和液柱射流的雾化过程。

说明书

技术领域

本发明涉及流体控制领域，特别地，涉及一种用于预测雾化过程的大涡模拟方法。

背景技术

在发动机燃烧室内，液体燃料的雾化决定了燃料与空气的混合效果，进而影响燃烧性能。雾化过程非常复杂，多种不稳定性(Kelvin-Helmholtzinstability,Rayleigh-Taylorinstability,Plateau-Rayleighinstability)同时存在，并伴有强烈的湍流，使得理论分析不具可行性。关于雾化已开展了大量的试验研究，但是由于雾化过程形成的液雾遮挡了液柱初始破碎过程，在观察和测量上造成了很大的困难。自从二十世纪七十年代，两相流的数值仿真取得了很大的进步，加深了对雾化机理的认识。

流体计算力学方法分为三种：雷诺平均方法，大涡模拟，直接数值模拟。雷诺平均方法只求解平均速度场，模化湍流运动对流场的影响。大涡模拟求解大尺度涡结构，模化小尺度涡结构对流动的影响。直接数值模拟求解所有尺度的涡结构。雾化过程中，湍流中的大尺度涡可扰动两相流界面，显著地影响液体射流的破碎过程，限制了雷诺平均方法的应用。由于大涡模拟的计算量比直接数值模拟小很多，大涡模拟更适合工程应用。

为了精确求解液柱和液滴的破碎过程，需跟踪液气界面。流行的界面跟踪方法有：Volumeoffluid(VOF)，LevelSet(LS),CoupledLSandVOF(CLSVOF)。其中，VOF方法(流体体积法)在20世纪70年代末由Hirt和Nichols等最先提出，基本思想是在欧拉网格系统上定义一个函数，根据每个网格内所含某种物质的体积量来定义在此网格上的值，然后用体积跟踪的方法求解方程，VOF方法可以精确地保证质量守恒，但是VOF函数的不连续性导致界面的构造很复杂且容易破碎。LS方法(水平集方法)可以很容易地构造界面，但是所得界面包围的液体质量不守恒。CLSVOF方法(水平集复合流体体积方法)可以很好地结合VOF和LS方法的优点，得到了广泛的应用。

由于液气界面两侧密度和流体粘性系数的不连续性，在求解控制方程时，常规的数值离散方法误差大，并造成算法的不稳定，液气密度比越大，算法越不稳定。为了得到收敛的结果，很多已发表的文献在数值仿真中采用较低的液气密度比，但是大多数雾化试验是采用高密度液体(如水、煤油、酒精)在大气环境中进行的，具有较高的液气密度比，现有技术无法将数值仿真结果与实验结果进行比较。因此，现有的两相流大涡模拟存在算法复杂、且外延液体速度不满足连续性方程，离散数值误差大等问题导致无法准确预测液滴破碎过程和液体射流的雾化过程的缺陷。

发明内容

本发明提供了一种用于预测雾化过程的大涡模拟方法，以解决现有两相流模拟方法导致的两相流液滴破碎及液体射流的雾化过程难以准确预测的技术问题。

本发明采用的技术方案如下：

一种用于预测雾化过程的大涡模拟方法，本发明方法基于真实速度场构造虚拟的液相速度场，并将构造的液相速度场应用于流动控制方程和界面输运方程的求解，模拟气液两相流的实时动态过程，以准确预测液滴破碎过程和液体射流的雾化过程。

进一步地，本发明用于预测雾化过程的大涡模拟方法包括：

步骤S10，根据水平集LS函数Φⁿ表示两相流界面，通过真实速度场Uⁿ和液相速度场U^Ln求解两相流控制方程，获得下一时间步对应的真实速度场Uⁿ⁺¹；

步骤S20，通过外延方法构造n+1时间步的液相速度场U^Ln+1；

步骤S30，通过外延液相速度场去散度化方法使U^Ln+1满足连续性方程：

步骤S40，利用构造的液相速度场U^Ln+1，通过水平集复合流体体积CLSVOF方法求解LS函数和VOF函数的输运方程，获得下一时刻的LS函数Φⁿ⁺¹和流体体积VOF函数Fⁿ⁺¹；

步骤S50，由气相转变为液相的控制体中，将真实速度场Uⁿ⁺¹重置为液相速度场U^Ln+1,即Uⁿ⁺¹＝U^Ln+1；

重复以上步骤S10至S50，以模拟两相流的实时动态过程。

进一步地，所述步骤S10中，对所述两相流控制方程进行空间过滤处理。

本发明具有以下有益效果：

本发明用于预测雾化过程的大涡模拟方法，基于真实速度场构造虚拟的液相速度场，并将构造的液相速度场应用于流动控制方程和界面输运方程的求解，模拟气液两相流的实时动态过程，提高了计算精度和稳定性，可以准确计算和预测液滴的破碎过程和液柱射流的雾化过程。

除了上面所描述的目的、特征和优点之外，本发明还有其它的目的、特征和优点。下面将参照图，对本发明作进一步详细的说明。

附图说明

构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1是本发明优选实施例大涡模拟方法的流程示意图；

图2是本发明优选实施例计算变量的分布示意图；

图3是本发明优选实施例液相速度场的构造示意图。

具体实施方式

以下结合附图对本发明的实施例进行详细说明，但是本发明可以由权利要求限定和覆盖的多种不同方式实施。

本发明的优选实施例提供了一种用于预测雾化过程的大涡模拟方法，本发明方法基于真实速度场构造虚拟的液相速度场，并将构造的液相速度场应用于流动控制方程和界面输运方程的求解，模拟气液两相流的实时动态过程，以准确预测液滴破碎过程和液体射流的雾化过程。

参照图1，本实施例控制方法包括：

步骤S10，根据水平集LS函数Φⁿ表示两相流界面，通过真实速度场Uⁿ和液相速度场U^Ln求解两相流控制方程，获得下一时间步对应的真实速度场Uⁿ⁺¹；

步骤S20，通过外延方法构造n+1时间步的液相速度场U^Ln+1；

步骤S30，通过外延液相速度场去散度化方法使U^Ln+1满足连续性方程：

步骤S40，利用构造的液相速度场U^Ln+1，通过水平集复合流体体积CLSVOF方法求解LS函数和VOF函数的输运方程，获得下一时刻的LS函数Φⁿ⁺¹和流体体积VOF函数Fⁿ⁺¹；

步骤S50，由气相转变为液相的控制体中，将真实速度场Uⁿ⁺¹重置为液相速度场U^Ln+1,即Uⁿ⁺¹＝U^Ln+1；

重复以上步骤S10至S50，以模拟两相流的实时动态过程。

作为一种较佳的方式，为了跟踪液气边界，引入两个函数：LS函数Φ和VOF函数F。LS函数Φ是到液气界面的变号距离函数。Φ＝0代表液气界面；在液体中Φ>0；在气体中Φ≤0。VOF函数F是每个计算单元中液体体积百分比。

优选地，对所述两相流控制方程(Navier-Stokes方程)进行空间过滤处理。经空间过滤后，连续方程变为：

$>\frac{\partial U_i}{\partial x_i}=0$>

其中，U_i为速度分量，x_i为位置坐标。

通过Smagorinsky涡粘性模型模化亚格子应力，动量方程变为：

$>\frac{\partial \left(U_i\right)}{\partial t}+\frac{\partial \left(U_i{U}_j\right)}{\partial x_j}=-\frac{1}{\rho }\frac{\partial P}{\partial x_i}+\frac{1}{\rho }\frac{\partial \left({\tau }_{ij}+{\tau }_{ij}^r\right)}{\partial x_j}+g_i+\frac{1}{\rho }{F}_i^{ST}$>

其中，P为压强，t为时间，ρ为密度，g_i为重力分量，为表面张力。τ_ij和分别为粘性应力张量和亚格子应力张量，且计算如下：

$>\left(\begin{array}{ccccc}{\tau }_{ij}=2{\mathrm{\mu S}}_{ij}& {\tau }_{ij}^r=2{\mu }_r{S}_{ii}& S_{ij}=\frac{1}{2}(\frac{\partial U_i}{\partial x_j}+\frac{\partial U_j}{\partial x_i})& {\mu }_r=\rho {\left(C_S\Delta \right)}^{2}S& S=\sqrt{2S_{ij}{S}_{ij}}\end{array}\right)$>

ρ＝ρ_G+(ρ_L-ρ_G)H(φ)μ＝μ_G+(μ_L-μ_G)H(φ)

μ和μ_r分别表示动力粘性系数和亚格子粘性系数，过滤宽度Δ取为当地计算单元体积的立方根，S_ij为应变张量，C_S为Smagorinsky系数。下标G和L分别表示气体和液体，H(φ)为Heaviside函数，表示如下：

$>H\left(\phi \right)=\left(\begin{array}{ccc}1& if& \phi >0\\ 0& if& \phi \le 0\end{array}\right)$>

表面张力为：

$>\left(\begin{array}{ccc}{F}_i^{ST}=\sigma \kappa \frac{\partial H}{\partial x_i}& \kappa =\frac{\partial n_i}{\partial x_i}& n_i=-\frac{1}{\sqrt{\frac{\partial \phi }{\partial x_k}\frac{\partial \phi }{\partial x_k}}}\frac{\partial \phi }{\partial x_i}\end{array}\right)$>

其中σ为表面张力系数，κ为曲率，ni为法向量分量。

VOF函数的控制方程为：

$>\frac{\partial F}{\partial t}+U_i\frac{\partial F}{\partial x_i}=0$>

LS函数的控制方程为：

$>\frac{\partial \phi }{\partial t}+U_i\frac{\partial \phi }{\partial x_i}=0$>

作为一种较佳的方式，为了计算下一时间步的速度场：

首先，通过对流项、扩散项和重力项计算中间速度场：

$>\frac{U_i^{*}-U_i^n}{\delta t}=-\frac{\partial \left(U_i^n{U}_j^n\right)}{\partial x_j}+\frac{1}{\rho }\frac{\partial ({\tau }_{ij}^n+{\tau }_{ij}^{r^n})}{\partial x_j}+g_i$>

其次，中间速度场通过压力项修正得到n+1时间步的速度场：

$>\frac{U_i^{n+1}-U_i^{*}}{\delta t}=-\frac{1}{\rho }\frac{\partial P^{n+1}}{\partial x_i}$>

因为n+1时间步的速度场满足连续性方程，通过对上一方程求散度可得如下压力泊松方程(通过此泊松方程可求解n+1时间步的的压力场)：

$>\frac{\partial }{\partial x_i}\left(\frac{1}{\rho }\frac{\partial P_i^{p+1}}{\partial x_i}\right)=\frac{1}{\delta t}\frac{\partial U_i^{*}}{\partial x_i}$>

参照图2，图2展示了计算变量的分布示意图，LS函数Φ、压强P、VOF函数F位于计算单元中心，速度以交错方式分布于相应的计算单元表面上。u，v是速度场在x和y方向上的分量。

参照图3，液相速度初始化为动量方程得到的速度场

$>u_{i-1/2,j}^L=u_{i-1/2,j}$>如果$>\left(\begin{array}{cc}{\phi }_{i-1/2,j}>0& ({\phi }_{i-1/2,j}=\frac{{\phi }_{i-1,j}+{\phi }_{i,j}}{2})\end{array}\right)$>

$>v_{i,j-1/2}^L=v_{i,j-1/2}$>如果$>\left(\begin{array}{cc}{\phi }_{i,j-1/2}>0& ({\phi }_{i,j-1/2}=\frac{{\phi }_{i,j-1}+{\phi }_{i,j}}{2})\end{array}\right)$>

气体中(φ≤0)的液相速度是通过将沿界面法向方向从液体向气体外延得到的，求解下面的外延方程到稳态：

$>\left(\begin{array}{ccc}\frac{\partial {\overrightarrow{U}}^L}{\partial \tau }+\overrightarrow{n}·▿{\overrightarrow{U}}^L=0& if& \phi \le 0\end{array}\right)$>

一阶前向欧拉方法用于时间离散，以液相速度分量u^L为例：

$>\left(\begin{array}{ccc}\frac{{u_{i-1/2,j}^L}^{n+1}-{u_{i-1/2,j}^L}^n}{\Delta \tau }=-{\left(n_x\frac{\partial u^L}{\partial x}\right)}_{i-1/2,j}^n-{\left(n_y\frac{\partial u^L}{\partial y}\right)}_{i-1/2,j}^n& if& {\phi }_{i-1/2,j}\le 0\end{array}\right)$>

伪时间步长Δτ＝0.3min(Δx_i-1,Δx_i,Δy_j-1,Δy_j,Δy_j+1)。一阶迎风格式用于空间离散：

$>{\left(\frac{\partial u^L}{\partial x}\right)}_{i-1/2,j}=\left(\begin{array}{ccc}\frac{u_{i-1/2,j}^L-u_{i-3/2,j}^L}{{\mathrm{\Delta x}}_{i-1}}& if& {\left(n_x\right)}_{i-1/2,j}>0\\ \frac{u_{i+1/2,j}^L-u_{i-1/2,j}^L}{{\mathrm{\Delta x}}_i}& if& {\left(n_x\right)}_{i-1/2,j}\le 0\end{array}\right)$>

作为较佳的方式，外延的液体速度应满足连续性条件

首先计算单元(i,j)中的速度源项：

$>S_{i,j}=u_{i-1/2,j}^L{\mathrm{\Delta y}}_j-u_{i+1/2,j}^L{\mathrm{\Delta y}}_j+v_{i,j-1/2}^L{\mathrm{\Delta x}}_i-v_{i,j\mathrm{+1/2}}^L{\mathrm{\Delta x}}_i$>

修正气体中的液相速度以满足连续性条件：

$>\left(\begin{array}{cc}{u}_{i-1/2,j}^L=u_{i-1/2,j}^L-a_w\frac{S_{i,j}}{A}|n_x{|}_{i-1/2,j}& u_{i+1/2,j}^L=u_{i+1/2,j}^L+a_e\frac{S_{i,j}}{A}|n_x{|}_{i+1/2,j}\end{array}\right)$>

$>\left(\begin{array}{cc}{v}_{i,j-1/2}^L=v_{i,j-1/2}^L-a_s\frac{S_{i,j}}{A}|n_y{|}_{i,j-1/2}& v_{i,j+1/2}^L=v_{i,j+1/2}^L+a_n\frac{S_{i,j}}{A}|n_y{|}_{i,j+1/2}\end{array}\right)$>

$>\left(\begin{array}{cc}{a}_w=\left(\begin{array}{cc}1& \begin{array}{cc}if& ({\phi }_{i-1/2,j}<0{\mathrm{and\phi }}_{i,j}>{\phi }_{i-1,j})\end{array}\\ 0& else\end{array}\right)& a_e=\left(\begin{array}{cc}1& \begin{array}{cc}if& ({\phi }_{i+1/2,j}<0{\mathrm{and\phi }}_{i,j}>{\phi }_{i+1,j})\end{array}\\ 0& else\end{array}\right)\end{array}\right)$>

$>\left(\begin{array}{cc}{a}_s=\left(\begin{array}{cc}1& \begin{array}{cc}if& ({\phi }_{i,j-1/2}<0{\mathrm{and\phi }}_{i,j}>{\phi }_{i,j-1})\end{array}\\ 0& else\end{array}\right)& a_n=\left(\begin{array}{cc}1& \begin{array}{cc}if& ({\phi }_{i,j+1/2}<0{\mathrm{and\phi }}_{i,j}>{\phi }_{i,j+1})\end{array}\\ 0& else\end{array}\right)\end{array}\right)$>

A＝a_e|n_x|_i-1/2,jΔy_j+a_w|n_x|_i+1/2,jΔy_j+a_s|n_y|_i,j-1/2Δx_i+a_n|n_y|_i,j+1/2Δx_i。

从以上的描述可以得知，本发明基于真实速度场构造虚拟的液相速度场，并将构造的液相速度场应用于流动控制方程和界面输运方程，模拟气液两相流的实时动态过程，提高了计算精度和稳定性，可以准确计算和预测液滴的破碎过程和液柱射流的雾化过程。且本发明通过液相速度外延算法构造下一时间步的液相速度场，通过外延液相速度无散度化处理，进一步控制了两相流模拟时由于离散数值导致的误差，提高了计算精度和稳定性。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。