[发明专利]一种基于SOLA-VOF的溃坝流动界面追踪方法

权利要求书

1.一种基于SOLA-VOF的溃坝流动界面追踪方法，包括SOLA-VOF算法，其特征在于：所述SOLA-VOF算法包括：

①、控制方程

设定溃坝流动为二维不可压缩层流流动，采用的控制方程如下：

连续性方程：

动量方程：

VOF方程：

对动量方程中粘度和密度采用两相混合值来修正计算：

υ＝F·υ_F+(1-F)·υ_A，ρ＝F·ρ_F+(1-F)·ρ_A (5)

式中，ρ_F、υ_F、ρ_A、υ_A分别为流体和气相的密度、运动粘度；

SOLA法采用有限差分显式格式来离散方程，利用散度定义离散连续性方程：

式中，D_i，j定义为网格单元(i,j)的散度，若D_i，j小于某一极小值，则认为连续性方程已经收敛，这种处理省去了每次迭代连续性方程的麻烦，只需要在每次迭代中判断速度是否符合精度要求；

②、VOF界面追踪

VOF法将单元内相流体的体积分数定义如下：

流体体积函数满足流体连续性方程，也即VOF方程：

VOF界面追踪包含两个部分：自由界面重构和单元间流量输运，界面重构是根据给定的F值确定每个单元内流体所占的区域，流量输运是按照重构好的界面确定下一时刻单元内流体的体积分数；

③、改进的单元流量运输算法

Hirt首先提出将单元内交界面简化为水平或竖直直线来重构界面，并采用施主-受主法来计算单元间流量输运值，为了确定单元内流体部分所占位置，需要构造界面单元函数Y(x_i)和X(y_i)：

Y(x_i)＝F_i，j-1δy_j-1+F_i，jδy_j+F_i，j+1δy_j+1 (8)

X(y_i)＝F_i-1，jδx_i-1+F_i，jδx_i+F_i+1，jδx_i+1 (9)

求出单元内流体界面的切线斜率：

式中，下标(i,j)代表单元中心位置，δx_i和δy_j为单元网格间距，i+1/2和j+1/2代表相邻两单元i和i+1、j和j+1的交界面；

如果单元界面近似为水平，反之，则近似为竖直，假设单元界面近似为水平线，此时，若则F代表的相流体处于重构水平线的下方，因此根据当前单元和四周单元的体积分数可以重构出整个界面形状；

获得当前界面重构形态后，为得到下一时刻的界面，还需要计算网格间的流量输运，Hirt的施主受主法采用如下公式计算单位时间dt内流过界面的流量dF：

dF＝min(F_AD|V_x|+CF，F_Ddx) (12)

CF＝max[(1-F_AD)|V_x|-(1-F_AD)dx，0] (13)

式中，下标D表示施主单元，即上游单元；A表示受主单元，即下游单元；AD表示A或D，取决于流体流动的方向和界面形态，CF为为修正施主单元多输入的流量，考虑到计算单元的上游单元不一定为满相流体单元，本研究将式(13)中1修正为上游单元即时流量F，修正后流量公式为：

CF＝max[(F-F_AD)|V_x|-(F-F_AD)dx，0] (14)

F＝max(F_D，F_DM) (15)

式中，DM表示施主单元的上游单元；

④、界面张力的处理

气液界面间存在界面张力，界面张力影响着相界面的变形、合并和破碎等过程，界面张力表达为：

p_s＝-σK (16)

式中，σ为界面张力系数，设为常数；K为曲率，与界面的斜率有关，当单元内界面近似为水平，反之当时，

对应SOLA显示格式离散，本研究先求出速度后，再引入界面张力对流动影响的修正，以增加计算精度：

⑤、壁面附着力的处理

壁面附着力一般是指液体分子与壁面固体分子之间的作用力，它会影响液体在固壁上的延展能力，并体现为“浸润”和“不浸润”，受温度、流体性质和壁面材质等影响，在数值计算中，附着力用接触角(壁面粘附角度)进行度量，接触角小于90°代表流体较易润湿其表面，度数越小，润湿性能越好；大于90°则认为流体不易润湿；

壁面附着力可以认为是界面张力在液体靠近固壁处的特殊情况，此时，对于紧靠壁面的边界单元，界面曲率的计算与其他位置稍有不同，要在边界外设置一层与边界单元相同步长的虚拟单元层，如对于右竖直壁面，附加层内每个单元的界面单元函数为：

Y(x_i+1)＝Y(x_i)+1/2(δx_i+δx_i+1)/tanθ (18)

X(y_i+1)＝X(y_j)-1/2(δy_j+δy_j+1)/tanθ (19)

θ是液体与固体壁面之间的接触角，该式是对竖直壁面进行处理，壁面上速度的修正与式(17)相同；

⑥、SOLA-VOF算法计算步骤

(1)根据初始条件或上一时刻的速度和压力，带入动量方程中，求出试算速度；

(2)将试算速度带入离散后的连续性方程中，求出D_i，j；

(3)判断每个单元的D_i，j是否满足D_i，jε，ε为收敛精度，此处取ε＝10^-3，若满足则进行步骤(6)；

(4)当D_i，jε，不满足收敛精度，需要修正压力；

(5)由修正的压力值来修正速度，执行步骤(2)；

(6)本时间步速度、压力已收敛，继续求解VOF方程计算单元体积分数，转到步骤(1)进行下一时间步计算。