一种基于SOLA-VOF的溃坝流动界面追踪方法

一种基于sola-vof的溃坝流动界面追踪方法
技术领域
1.本发明涉及自然灾害模拟技术领域，具体为一种基于sola-vof的溃坝流动界面追踪方法。


背景技术：

2.溃坝是一种严重的灾害，决堤后湍急的洪水会危及到百姓的人身安全和财产安全。采用模拟方法预测溃坝的流动趋势和特性，有利于及早规避风险、减少损失。
3.自由界面追踪问题是多相流领域研究的重要课题。vof(volume of fluid)法通过追踪网格中某一相流体占体积分数来确定界面位置，具有质量缺失少、计算量小、精度高等优点，应用广泛。在vof界面重构方案中，施主—受主(donor-acceptor)法(根据单元网格内流速方向将单元分为施主和受主单元)，同时考虑单元上
4.游和下游的流量信息来获得界面形状，减小了界面模糊。sola(solution algorithm)法是一种求解速度压力耦合的高效算法，它在mac法和simple法基础上加以改进，采用有限差分法离散控制方程。在求解速度、压力时，与simple算法采用的双方程迭代方法不同，只迭代动量方程，将连续性方程作为收敛判据却不参与迭代，故而大幅缩减速度压力迭代求解时间。在提出sola算法不久，nichols等首次将sola算法与自由表面追踪问题结合求解，提出sola-vof耦合算法完成了气液两相流界面捕捉。2018年fahime修正sola-vof算法模拟了非牛顿流体中气泡的生长过程，hamid则模拟了纳米气泡生成过程；2019年hao采用sola-vof模拟了压铸中内部孔隙的产生，将sola方法扩展到不可压缩流体的计算中。除了在铸造和气泡模拟上的应用外，该算法在其他领域也有众多应用。但在使用sola-vof方法求解两相交界面时，界面附近流体琐碎会随时间逐渐累积，造成误差。
5.为满足工程上对该算法的计算精度要求，本文自编fortran程序改进了sola-vof算法，以减少界面杂质，并根据改进算法模拟了溃坝的界面演化过程，同时分析了界面张力和壁面附着力对相界面模拟的影响。


技术实现要素：

6.(一)解决的技术问题
7.针对现有技术的不足，本发明提供了一种基于sola-vof的溃坝流动界面追踪方法，解决了上述背景技术中所存在的问题。
8.(二)技术方案
9.为实现上述目的，本发明提供如下技术方案：一种基于sola-vof的溃坝流动界面追踪方法，包括sola-vof算法，所述sola-vof算法包括：
10.①
、控制方程
11.设定溃坝流动为二维不可压缩层流流动，采用的控制方程如下：
12.连续性方程：
13.动量方程：
[0014][0015]
vof方程：
[0016]
对动量方程中粘度和密度采用两相混合值来修正计算：
[0017]
υ＝f
·
υf (1-f)
·
υa，ρ＝f
·
ρf (1-f)
·
ρaꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0018]
式中，ρf、υf、ρa、υa分别为流体和气相的密度、运动粘度；
[0019]
sola法采用有限差分显式格式来离散方程。利用散度定义离散连续性方程：
[0020][0021]
式中，d
i，j
定义为网格单元(i,j)的散度，若d
i，j
小于某一极小值，则认为连续性方程已经收敛，这种处理省去了每次迭代连续性方程的麻烦，只需要在每次迭代中判断速度是否符合精度要求；
[0022]
②
、vof界面追踪
[0023]
vof法将单元内相流体的体积分数定义如下：
[0024][0025]
流体体积函数满足流体连续性方程，也即vof方程：
[0026][0027]
vof界面追踪包含两个部分：自由界面重构和单元间流量输运，界面重构是根据给定的f值确定每个单元内流体所占的区域，流量输运是按照重构好的界面确定下一时刻单元内流体的体积分数；
[0028]
③
、改进的单元流量运输算法
[0029]
hirt首先提出将单元内交界面简化为水平或竖直直线来重构界面，并采用施主-受主法来计算单元间流量输运值，为了确定单元内流体部分所占位置，需要构造界面单元函数y(xi)和x(yi)：
[0030]
y(xi)＝f
i，j-1
δy
j-1
f
i，j
δyj f
i，j 1
δy
j 1
ꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0031]
x(yi)＝f
1-1，j
δx
i-1
f
i，j
δxi f
i 1，j
δx
i 1
ꢀꢀꢀꢀꢀꢀꢀꢀ
(9)
[0032]
求出单元内流体界面的切线斜率：
[0033][0034][0035]
式中，下标(i,j)代表单元中心位置，δxi和δyj为单元网格间距，i 1/2和j 1/2代表相邻两单元i和i 1、j和j 1的交界面；
[0036]
如果单元界面近似为水平，反之，则近似为竖直。假设单元界面近似为水平线，此时，若则f代表的相流体处于重构水平线的下方。因此根据当前单元和四周单元的体积分数可以重构出整个界面形状；
[0037]
获得当前界面重构形态后，为得到下一时刻的界面，还需要计算网格间的流量输运，hirt的施主受主法采用如下公式计算单位时间dt内流过界面的流量df：
[0038]
df＝min(f
ad
|v
x
| cf，fddx)
ꢀꢀꢀꢀꢀꢀꢀꢀ
(12)
[0039]
cf＝max[(1-f
ad
)|v
x
|-(1-f
ad
)dx，0]
ꢀꢀꢀꢀꢀꢀꢀꢀ
(13)
[0040]
式中，下标d表示施主单元，即上游单元；a表示受主单元，即下游单元；ad表示a或d，取决于流体流动的方向和界面形态，cf为为修正施主单元多输入的流量，考虑到计算单元的上游单元不一定为满相流体单元，本研究将式(13)中1修正为上游单元即时流量《f》，修正后流量公式为：
[0041]
cf＝max[(《f》-f
ad
)|v
x
|-(《f》-f
ad
)dx，0]
ꢀꢀꢀꢀꢀꢀꢀꢀ
(14)
[0042]
《f》＝max(fd，f
dm
)
ꢀꢀꢀꢀꢀꢀꢀꢀ
(15)
[0043]
式中，dm表示施主单元的上游单元；
[0044]
④
、界面张力的处理
[0045]
气液界面间存在界面张力，界面张力影响着相界面的变形、合并和破碎等过程。界面张力表达为：
[0046]
ps＝-σk
ꢀꢀꢀꢀꢀꢀꢀꢀ
(16)
[0047]
式中，σ为界面张力系数，设为常数；k为曲率，与界面的斜率有关，当单元内界面近似为水平，反之当时，
[0048]
对应sola显示格式离散，本研究先求出速度后，再引入界面张力对流动影响的修正，以增加计算精度：
[0049][0050][0051]
⑤
、壁面附着力的处理
[0052]
壁面附着力一般是指液体分子与壁面固体分子之间的作用力，它会影响液体在固壁上的延展能力，并体现为“浸润”和“不浸润”，受温度、流体性质和壁面材质等影响。在数值计算中，附着力用接触角(壁面粘附角度)进行度量，接触角小于90
°
代表流体较易润湿其表面，度数越小，润湿性能越好；大于90
°
则认为流体不易润湿；
[0053]
壁面附着力可以认为是界面张力在液体靠近固壁处的特殊情况，此时，对于紧靠壁面的边界单元，界面曲率的计算与其他位置稍有不同，要在边界外设置一层与边界单元相同步长的虚拟单元层。如对于右竖直壁面，附加层内每个单元的界面单元函数为：
[0054]
y(x
i 1
)＝y(xi) 1/2(δxi δx
i 1
)/tanθ
ꢀꢀꢀꢀꢀꢀꢀꢀ
(18)
[0055]
x(y
i 1
)＝x(yj)-1/2(δyj δy
j 1
)/tanθ
ꢀꢀꢀꢀꢀꢀꢀꢀ
(19)
[0056]
θ是液体与固体壁面之间的接触角，该式是对竖直壁面进行处理，壁面上速度的修正与式(17)相同；
[0057]
⑥
、sola-vof算法计算步骤
[0058]
(1)根据初始条件或上一时刻的速度和压力，带入动量方程中，求出试算速度；
[0059]
(2)将试算速度带入离散后的连续性方程中，求出d
i，j
；
[0060]
(3)判断每个单元的d
i，j
是否满足d
i，j
《ε，ε为收敛精度，此处取ε＝10-3
，若满足则进行步骤(6)；
[0061]
(4)当d
i，j
》ε，不满足收敛精度，需要修正压力；
[0062]
(5)由修正的压力值来修正速度，执行步骤(2)；
[0063]
(6)本时间步速度、压力已收敛，继续求解vof方程计算单元体积分数，转到步骤(1)进行下一时间步计算。
[0064]
(三)有益效果
[0065]
本发明提供了一种基于sola-vof的溃坝流动界面追踪方法，具备以下有益效果：
[0066]
本发明采用改进的sola-vof法模拟溃坝流动界面演化过程，并分析界面张力和壁面附着力对界面流动的影响。对比分析可得，(1)改进的施主—受主算法能够提高计算精度，有效减少界面琐碎流体杂质，光滑界面曲线；(2)采用该算法模拟溃坝流动过程与文献结果吻合良好；(3)界面张力对流动表现出抑制作用，减缓流体延展；(4)流体与壁面润湿性能越好，流体流动延展性越好。
附图说明
[0067]
图1为本发明中施主受主法的界面重构对比示意图；
[0068]
图2为本发明中液体与容器壁间接触角示意图；
[0069]
图3为本发明中测试流场下改进公式前后流动界面对比示意图；
[0070]
图4为本发明中溃坝水流流动过程示意图；
[0071]
图5为本发明中界面张力对液面演变的影响示意图；
[0072]
图6为本发明中溃坝不同位置处单元液相体积分数f曲线示意图；
[0073]
图7为本发明中不同接触角下流体体积分数对比示意图。
具体实施方式
[0074]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0075]
本发明提供一种技术方案：一种基于sola-vof的溃坝流动界面追踪方法，包括sola-vof算法，所述sola-vof算法包括：
[0076]
①
、控制方程
[0077]
设定溃坝流动为二维不可压缩层流流动，采用的控制方程如下：
[0078]
连续性方程：
[0079]
动量方程：
[0080][0081]
vof方程：
[0082]
对动量方程中粘度和密度采用两相混合值来修正计算：
[0083]
υ＝f
·
υf (1-f)
·
υa，ρ＝f
·
ρf (1-f)
·
ρaꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0084]
式中，ρf、υf、ρa、υa分别为流体和气相的密度、运动粘度；
[0085]
sola法采用有限差分显式格式来离散方程。利用散度定义离散连续性方程：
[0086][0087]
式中，d
i，j
定义为网格单元(i,j)的散度，若d
i，j
小于某一极小值，则认为连续性方程已经收敛，这种处理省去了每次迭代连续性方程的麻烦，只需要在每次迭代中判断速度是否符合精度要求；
[0088]
②
、vof界面追踪
[0089]
vof法将单元内相流体的体积分数定义如下：
[0090][0091]
流体体积函数满足流体连续性方程，也即vof方程：
[0092][0093]
vof界面追踪包含两个部分：自由界面重构和单元间流量输运，界面重构是根据给定的f值确定每个单元内流体所占的区域，流量输运是按照重构好的界面确定下一时刻单元内流体的体积分数；
[0094]
③
、改进的单元流量运输算法
[0095]
hirt首先提出将单元内交界面简化为水平或竖直直线来重构界面，并采用施主-受主法来计算单元间流量输运值，如图1所示，为了确定单元内流体部分所占位置，需要构造界面单元函数y(xi)和x(yi)：
[0096]
y(xi)＝f
i，j-1
δy
j-1
f
i，j
δyj f
i，j 1
δy
j 1
ꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0097]
x(yi)＝f
i-1，j
δx
i-1
f
i，j
δxi f
i 1，j
δx
i 1
ꢀꢀꢀꢀꢀꢀꢀꢀ
(9)
[0098]
求出单元内流体界面的切线斜率：
[0099][0100][0101]
式中，下标(i,j)代表单元中心位置，δxi和δyj为单元网格间距，i 1/2和j 1/2代表
相邻两单元i和i 1、j和j 1的交界面；
[0102]
如果单元界面近似为水平，反之，则近似为竖直。假设单元界面近似为水平线，此时，若则f代表的相流体处于重构水平线的下方。因此根据当前单元和四周单元的体积分数可以重构出整个界面形状；
[0103]
获得当前界面重构形态后，为得到下一时刻的界面，还需要计算网格间的流量输运，hirt的施主受主法采用如下公式计算单位时间dt内流过界面的流量df：
[0104]
df＝min(f
ad
|v
x
| cf，fddx)
ꢀꢀꢀꢀꢀꢀꢀꢀ
(12)
[0105]
cf＝max[(1-f
ad
)|v
x
|-(1-f
ad
)dx，0]
ꢀꢀꢀꢀꢀꢀꢀꢀ
(13)
[0106]
式中，下标d表示施主单元，即上游单元；a表示受主单元，即下游单元；ad表示a或d，取决于流体流动的方向和界面形态，cf为为修正施主单元多输入的流量，考虑到计算单元的上游单元不一定为满相流体单元，本研究将式(13)中1修正为上游单元即时流量《f》，修正后流量公式为：
[0107]
cf＝max[(《f》-f
ad
)|v
x
|-(《f》-f
ad
)dx，0]
ꢀꢀꢀꢀꢀꢀꢀꢀ
(14)
[0108]
《f》＝max(fd，f
dm
)
ꢀꢀꢀꢀꢀꢀꢀꢀ
(15)
[0109]
式中，dm表示施主单元的上游单元；
[0110]
④
、界面张力的处理
[0111]
气液界面间存在界面张力，界面张力影响着相界面的变形、合并和破碎等过程。界面张力表达为：
[0112]
ps＝-σk
ꢀꢀꢀꢀꢀꢀꢀꢀ
(16)
[0113]
式中，σ为界面张力系数，设为常数；k为曲率，与界面的斜率有关，当单元内界面近似为水平，反之当时，
[0114]
对应sola显示格式离散，本研究先求出速度后，再引入界面张力对流动影响的修正，以增加计算精度：
[0115][0116][0117]
⑤
、壁面附着力的处理
[0118]
壁面附着力一般是指液体分子与壁面固体分子之间的作用力，它会影响液体在固壁上的延展能力，并体现为“浸润”和“不浸润”，受温度、流体性质和壁面材质等影响。在数值计算中，附着力用接触角(壁面粘附角度)进行度量，图2为液体靠近固壁的接触角示意图，θ为液体与壁面间的接触角，接触角小于90
°
代表流体较易润湿其表面，度数越小，润湿性能越好；大于90
°
则认为流体不易润湿；
[0119]
壁面附着力可以认为是界面张力在液体靠近固壁处的特殊情况，此时，对于紧靠壁面的边界单元，界面曲率的计算与其他位置稍有不同，要在边界外设置一层与边界单元
相同步长的虚拟单元层。如对于右竖直壁面，附加层内每个单元的界面单元函数为：
[0120]
y(x
i 1
)＝y(xi) 1/2(δxi δx
i 1
)/tanθ
ꢀꢀꢀꢀꢀꢀꢀꢀ
(18)
[0121]
x(y
i 1
)＝x(yj)-1/2(δyj δy
j 1
)/tanθ
ꢀꢀꢀꢀꢀꢀꢀꢀ
(19)
[0122]
θ是液体与固体壁面之间的接触角，该式是对竖直壁面进行处理，壁面上速度的修正与式(17)相同；
[0123]
⑥
、sola-vof算法计算步骤
[0124]
(1)根据初始条件或上一时刻的速度和压力，带入动量方程中，求出试算速度；
[0125]
(2)将试算速度带入离散后的连续性方程中，求出d
i，j
；
[0126]
(3)判断每个单元的d
i，j
是否满足d
i，j
《ε，ε为收敛精度，此处取ε＝10-3
，若满足则进行步骤(6)；
[0127]
(4)当d
i，j
》ε，不满足收敛精度，需要修正压力；
[0128]
(5)由修正的压力值来修正速度，执行步骤(2)；
[0129]
(6)本时间步速度、压力已收敛，继续求解vof方程计算单元体积分数，转到步骤(1)进行下一时间步计算。
[0130]
结果分析
[0131]
1、本文选用剪切流和旋转流流场来验证上述施主—受主法中流量改进算法的正确性。剪切流速度分布为：
[0132][0133]
式中，剪切流参考点取(x0，y0)＝(0.5m，0.5m)，流体圆心在点(0.5m,0.3m)，半径r＝0.2m，计算区间为[0,1m]
×
[0,1m]，时间步长为0.0001s，空间步长δx＝δy＝0.001m。
[0134]
旋转流速度分布为：u(x,y)＝-π(y-0.5)m/s,v(x,y)＝π(x-0.5)m/s，半径为r＝0.3m，其它设置与剪切流相同，改进算法前后对比结果如图3所示。
[0135]
由图3可以看出，修正后的流体边界比修正前更加清晰平滑，且能明显看到流体界面附近的琐碎减少，证实了修正算法的正确性.
[0136]
2、溃坝模拟结果及实验对比
[0137]
为了验证改进的sola-vof法追踪界面的精确性，本研究模拟了二维溃坝的界面演变过程，并与实验数据和学者张嫚嫚的模拟结果做参照对比。
[0138]
如图4所示，溃坝实验容器长1.6m
×
高1m，右侧有一液柱长0.6m
×
高0.3m。液体粘度υf＝8.9
×
10-7
m2/s，密度ρf＝997kg/m3；气体粘度为υc＝1.7894
×
10-5
m2/s，密度ρc＝1.25kg/m3，界面张力系数σ＝0.072n/m(此处为水—空气)。采用网格数320
×
200，网格单元尺寸为δx＝δy＝0.005m。左、右、下侧壁面采用无滑移边界，上壁面采用自由滑移边界。实验与模拟结果如图4。
[0139]
由图4可以看出，不同时刻模拟结果与实验结果基本吻合，图c中流体左侧上扬高度与实验结果相比较低，这是因为本模型采用层流模型，流动不够剧烈。0.3s时刻，实验观测的流体界面前沿到达总长2/3处；本模拟结果中，界面前沿到达0.53m处，与实验结果基本一致。因此由以上结果可以认为sola-vof方法能很好追踪界面流动状态。
[0140]
3、界面张力影响
[0141]
为了探究界面张力对溃坝液体流动的影响，本文选取了图5所示两时刻下是否考
虑界面张力对应的液体界面包络线，可以看出考虑界面张力下的等值线起伏更加复杂，这是因为张力作用在交界面上，对界面产生挤压，与仅考虑重力的情况下相比，流体受力更加复杂，界面波动更明显。同时，在图5中0.3s和1.0s是否考虑张力的结果图中，取y＝0.1m和y＝0.15m线上各点体积分数f值进行比较，见图6。图6中斜线段上的点代表气液混相单元的体积分数，可以看出，在混相区内，考虑界面张力的影响导致流体体积分数更小(如0.3s，y＝0.1m处，考虑张力时流体体积分数f＝0.038《未考虑张力时f＝0.108)，这是由于界面张力作用阻碍了界面处液体的流动所致。
[0142]
4、壁面附着力的影响
[0143]
流体与壁面的润湿性能与流体属性和壁面粗糙度等因素有关。而这些影响在壁面上可以直观的表示为接触角的不同。为了探究不同条件下接触角引起的壁面附着力的差异对整个流体流动的影响，分别模拟了接触角为30、60、90
°
下溃坝的流动。选取0.3s时刻直线x＝1m上体积分数进行分析，如下图，在y＝0.14～0.16m的单元为气液两相单元，90
°
接触角对应的体积分数最小，30
°
对应最大。且在x＝1m，y＝0.15m处，30、60、90
°
接触角下单元流体体积分数f分别为0.298》0.189》0.101。可见接触角越小，润湿性能越好，导致流体更易在壁面附近流动，因此单元流体体积分数更大，如图7所示。
[0144]
需要说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。
[0145]
尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由所附权利要求及其等同物限定。