应用改进流体体积法的楔形体斜向入水研究

谢 行， 任慧龙， 陶凯东， 冯亿坤

（1.哈尔滨工程大学 船舶工程学院，哈尔滨150001；2.黄淮学院 智能制造学院，河南 驻马店463000）

典型的入水冲击包含流体爬升、流动的分离、破碎、空泡的形成与发展等阶段，涉及气－液－固三者的相互作用.在此过程中伴随着自由面的大变形，甚至会出现流体的翻卷、破碎、气穴等复杂物理现象，对这些现象进行精确模拟需要极高的数值精度.

由于入水冲击的复杂性，实际数值模拟中对其进行了一些简化，如剖面形式简化为楔形体，流体性质不考虑黏性，忽略空气效应.Zhao等[1]采用边界元法（BEM）对二维有限宽度楔形体的入水问题进行了研究，通过对射流进行合理的截断，避免了时间步进中数值的发散.Qin等[2]基于势流理论，将自由面和物面条件进行线性化，推导了求解速度势的半解析公式，通过修正伯努利方程，研究了非对称楔形体的入水问题，结果与Xu等[3]的全非线性方法相符.上述均为单向流模型的方法，对入水后期流体的破碎、气穴等难以进行合理模拟.

为了考虑空气的效应，基于两相流模型的流体体积（VOF）算法被加以应用.Krastev等[4]通过VOF模型求解两相流的雷诺平均（RANS）方程，合理地模拟了楔形体非对称入水过程中气泡的产生和演变.Xie等[5]基于Fluent软件对具有倾斜角的外飘剖面进行了数值模拟，采用VOF法追踪自由面，模拟过程中流体先发生了分离，之后撞击到物体表面，发生了二次砰击现象.王易君等[6]应用VOF模型和动网格技术实现了平底结构的入水模拟，入水过程中出现了明显的空气垫.在现有VOF法中，相方程不但在离散求解过程中存在数值耗散现象，而且通过插值方法对气－液交界面处压力的不连续性进行处理时，也会造成相应的数值误差.对于自由面的重构，通过VOF法求得的体积分数存在一个从0到1的过渡域，这导致重构的相界面具有明显的厚度层，难以准确地计算界面的法向和曲率.

本文基于 Navier－Stokes（N－S）方程建立非对称入水问题的气－液两相流模型，在Fluent软件中，通过自定义函数（UDF），引入压力连续函数和人工对流项，并耦合水平集（Level－set）函数对现有VOF法进行改进.应用改进的VOF法对楔形体的斜向入水问题进行数值模拟，并对入水过程中砰击压力和自由面特性进行分析.

1 数值模型及数值过程

1.1 数值模型

VOF法将气－液两相流视为密度可变的单项流模型，利用体积分数来确定自由面的位置和形状.假定气液不相容，其连续性方程、动量方程和相方程分别为

式中：U为流场速度；t为时间；p为压力；ρ为流体密度；μ为黏度；α1为气体的体积分数（基本相），是随时间变化的函数.在数值计算中，当α1＝1时表示该网格内为空气，当α1＝0时表示网格内为水，当0＜α1＜1时，表示网格内为水和空气的交界面.为保证相分数严格介于[0，1]，可采用高分辨率的离散模式如 HRIC或CICSAM求解式（3），但这些离散模式在时间步进过程中会出现数值耗散，导致相界面的模糊性.为此，可通过添加人工对流项对相界面附近的体积分数进行挤压来缓解界面的模糊性.根据两相流理论，假设两相的相对速度Ur＝Ul－U2，添加对流项将式（3）转化为

式中：Δ·（α1α2Ur）为人工对流项；α2为水的体积分数，α2＝1－α1.通过体积分数，将流体密度ρ和黏度μ更新为

式中：ρ1和ρ2分别为空气和水的密度；μ1和μ2分别为空气和水的黏度.由于VOF法中并不存在所谓的多个相速度，故将相对速度采用单向流的速度进行模化，即

式中：cα为相界面压缩因子，cα＝0表示无压缩效果，cα越大压缩效应越明显.另外，Ur的主要作用在于压缩相界面来保持尖锐性，应该在相界面的法向而不是切向上进行压缩，否则会引起虚假的扩散，因此Ur的方向与界面法向相Δ同.在两相流的交Δ界处，引入某一点的法向量记为α，如图1所示在水和空气区域内为0，在自由表面处不为0.此时交界面的单位法向量n可表示为

图1 交界处示意图Fig.1 Schematic diagram of interface

为了避免气－液交界面两侧压力的不连续性，引入Brackbill等[7]提出的连续表面张力（CSF）模型，在界面处定义一个连续函数.将过渡域内的压力定义为一个连续函数

式中：c为界面处的位置函数；C为表面张力系数.则整个区域的压力可修正为

式中：κ为界面曲率，κ＝－ ·n.基于式（2），结合式（10），有

式（1）、（4）和（11）构成了改进 VOF法的连续性方程、相方程和动量方程.可以看出式（11）比式（2）多了一个源项κCΔα.

VOF法中由于界面不连续性的影响难以准确计算界面的法向和曲率，为了精确捕捉自由面，采用Level－set函数来保持界面的光滑性.通过Level－set函数φ（x，y，t）定义 VOF的体积函数F（α1，t）来实现 VOF和Level－set的耦合（CLSVOF），则有

式中：Ω为网格域；H为阶跃函数，记为

这种界面体积分数的光顺可在Fluent软件中通过VOF模型求解设置中勾选耦合CLSVOF模型直接实现.

1.2 数值过程

当物体以给定的入水速度下落时，在每个时刻t，具体的数值求解流程如下：① 取cα为0.5，对相方程采用有限差分法求解相体积分数α1，更新流体密度ρ和黏度μ；② 对式（11）进行有限体积离散，修正项通过添加表面张力系数来实现源项的添加，时间项采用欧拉隐式离散，压力项采用PRESTO！，速度项采用二阶迎风格式，迭代求解离散后的非线性方程组，从而给出速度的预测值；③ 结合式（1）采用PISO算法对流域压力场和速度场进行修正[8].应用动网格技术，根据给定的入水速度，预测下一时刻物体的位置，并采用局部重构方法对流域网格进行重新划分.

2 收敛性研究

选取文献[9]中的30°楔形体进行网格收敛性的研究，初始入水速度为2.42m／s.以3种不同的网格尺寸对流域进行划分，记为Mesh1、Mesh2和Mesh3，具体网格尺寸信息和计算时间见表1.图2给出了不同网格尺寸下数值结果与文献[9－10]的BEM的对比情况.当网格尺寸为Mesh1时，最大垂向力比实验值小12%；当网格尺寸减小到Mesh2时，结果偏大7%，主要原因是由于实验的三维效应，在文献[9]中，二维预测的最大垂向力比实验值大8%.当网格尺寸进一步减至Mesh3时，砰击力趋于稳定，这表明当网格尺寸足够小时，数值解趋于一致.对于砰击压力的结果（压力位置距低端尖角的垂向距离为45mm），虽然实验中并未给出，但从文献[10]结果的对比可以看出，数值的精度满足要求.根据垂向力和砰击压力的对比结果，Mesh2这一网格尺寸可以保证预测的精度.为了节省计算时间，本文的流域网格生成均采用Mesh2中的网格尺寸参数.

表1 用于收敛性分析的不同网格尺寸Tab.1 Different mesh sizes in mesh convergence analysis

图2 网格的收敛性分析Fig.2 Mesh convergence analysis

3 数值结果

3.1 研究对象

本研究以楔形体为研究对象，入水砰击模型如图3所示.为了便于对比分析，选取文献[11]中楔形体进行建模.底升角β＝20°，倾斜角为θ，计算中分别选为5°、10°和15°.R为迎风侧，底升角为β－θ；L为背风侧，底升角为β＋θ；F为气－液交界面.模型总质量为44kg，单边长度为0.205m，初始入水速度v＝3.13m／s.实验过程中每侧布置4个压力测点，具体位置如图4所示.整个流域尺寸为4m×3m，初始水深为2m，近楔形体处局部网格如图5所示.入水过程中忽略重力的影响，t＝0表示剖面底部刚好接触静水面.

图3 非对称入水砰击模型示意图Fig.3 Schematic diagram of asymmetrical water entry model

图4 实验中压力测点布置示意图[11]Fig.4 Schematic diagram of pressure measured positions in the experiment[11]

图5 近物体表面网格示意图Fig.5 Schematic diagram of mesh near the body

3.2 改进方法的影响分析

为探讨改进后的VOF法对自由面重构的分辨率及压力分布的影响，首先需要验证VOF法对气－液交界面捕捉的精确度.由于文献[11]实验中并未给出液面抬升的实验图片，本研究选取文献[12]中给出的自由面进行对比分析.数值模拟中垂向入水速度为3.81m／s，横向速度为7.62m／s，楔形体底升角为37°，倾斜角为5°.图6给出了不同时刻下实验和数值模拟的自由面的对比情况，可以看出VOF法能够较好地反映入水过程中的气－液交界面.

图7给出了楔形体斜向入水过程中分别采用改进前（左侧）和改进后（右侧）VOF法获得自由面的对比情况.模拟过程中网格的总量大约为3.7×104，改进前和改进后的计算时间分别为2.6和3.0 h.可以看出改进后所捕捉到的自由面厚度层明显减小，而且轮廓更加尖锐.当倾斜角增大时（图7（b）），左侧的自由面更加模糊，而右侧仍能保持较好的清晰度.这是由于倾斜角越大，近物体处的自由面上相关变量（如密度、黏度和压力）的梯度越大，仅采用插值方式（改进前）对体积分数求解存在较大误差，而CSF模型的引入保证了交界面处的密度和黏度充分光顺.在改进前的方法中，自由面的模糊性可通过增加网格密度进行改善，比如当网格数量达到4.9×104时（增加大约1／3），自由面清晰度可以达到与改进后的方法相同的效果，但这种措施不可避免地增加了计算的时间.

图6 实验和数值模拟的自由面对比Fig.6 Comparison of free surface obtained by experimental and numerical results

图7 改进前后自由面形状对比Fig.7 Comparison of free surface before and after improvement

图8给出了不同数值方法和砰击压力沿物体表面分布的对比情况（P2峰值时刻）.图中：无因次化量Cp＝2ps／（ρ2v2），ps为砰击压力；x′为楔形体的横向长度x与浸湿半宽CL的比值，即x′＝x／CL；紧致插值曲线（CIP）法为文献[13]的结果，半解析解为文献[14]的结果，实验结果来自于文献[11].改进前的VOF法尽管在峰值上与修正的结果差异不大，但在峰值之间的部分有所偏大，这种偏差在一定程度上降低了两侧压力的不对称性，进而会导致扭转力矩的减小.另外在图8（a）中半解析解的结果比实验值偏大18%，而这种误差在图8（b）中增加到40%，相比而言，CIP法和改进VOF法的结果始终维持在10%左右.由图7（b）中自由面可见，入水过程中包裹了一个气泡，会导致气垫效应，进而减缓砰击压力.由于半解析解是基于单向流模型，气体效应被忽略，而CIP法和VOF法是基于气－液的两相流模型，合理地计及了气体的影响，这也正是它们之间存在较大误差的原因所在.

图8 改进前后压力分布对比Fig.8 Comparison of pressure distribution before and after improvement

3.3 砰击载荷特性分析

图9给出了不同θ下迎风侧的压力时历（文献[11]实验中P1点），其中对称结果表示两侧底升角均为β－θ的对称楔形体的数值结果.为了清晰地对比结果，将计算的非对称结果后移2ms，对称结果前移2ms，实验结果保持不变.

由图9可以看出：非对称结果与实验值相当接近（误差小于5%），说明了本文预测结果的合理性；在图9（a）中对称结果比实验值偏大15%，这种差异性在底升角更小的图9（b）中达到了28%，这主要是由于斜向入水的过程中，背风侧的底升角β＋θ较大，对流体的阻碍作用要小于两侧底升角都为β－θ的情况.研究结果表明：对于底升角β的楔形体以倾斜角θ入水时，如果单纯地采用底升角为β－θ情况下的对称楔形体进行估算，砰击压力会造成结果的偏大；对称结果的压力上升速度要缓于实验值，相比于对称的情况，非对称情况下迎风侧流体的爬升速度更快些.

图9 迎风侧的砰击压力时历Fig.9 Pressure time histories on the upwind side

图10给出了入水过程中两侧砰击压力的分布情况.随着时间的增加（t≤14.4ms），迎风侧的压力开始高于背风侧的压力，两侧压力峰值的最大相对误差达到97%（t＝10.8ms）.之后迎风侧压力减小，背风侧压力仍在上升，在t＝21.6ms时超过迎风侧压力，最后两侧压力趋于对称化（t＝25.2ms）.事实上这种复杂的变化主要是由于流体从折角处分离造成的.在t＝14.4ms时迎风侧浸入水中，导致该侧压力开始减小.而此时背风侧尚未被流体完全淹没，压力仍在增加，直到t＝21.6ms左右，背风侧浸入水中，该侧压力也开始减小.当物体完全浸入水中后，由于物体处流体的波动开始减小，此时的倾斜角对两侧压力分布的影响可以忽略（t＝25.2ms）.图10（b）中出现了与图10（a）相似的规律性，但不同的是两侧的差异性达到307%（t＝7.2ms）.需要注意的是，底升角较小时迎风侧包裹了一个气泡，如图11所示，此时物面上的压力仍是连续的，这也就意味着楔形体与空气接触的部位的压力并非为0.由于此时流体会对气体域产生挤压，而楔形体作为空气域的一个接触部位，自然会承受这种挤压作用，从而出现了大于大气压的压力值.这种气穴效应，如果采用单向流模型，由于自由面边界条件的限制，水和空气的交界处压力为0（大气压），会导致与空气接触部位的压力为0.这种现象在文献[15]的数值与实验的对比结果中也有所体现，在数值结果中气穴处始终为0，而实验中却出现了明显的砰击压力.

图10 不同时刻下两侧压力分布Fig.10 Pressure distribution at different instants

图11 t＝21.6ms时的自由面形状Fig.11 Free surface at the instant of t＝21.6ms

4 结论

应用改进VOF法研究了楔形体斜向入水问题，得出以下结论：

（1）改进VOF法能够有效保持自由面的尖锐特征，倾斜角较大时，这种优势比改进前更加明显.

（2）随着倾斜角的增大，迎风侧砰击压力迅速增加，两侧压力的不对称性更加明显，而对称情况下的砰击压力偏大20%左右.

（3）当底升角小于10°时，空气对砰击压力的影响显著.计及空气效应的结果与实验的误差控制在10%以内，而未考虑空气影响的结果误差高达40%.

当砰击过程中的空气效应不明显时（底升角大于10°），一些基于单相流模型的方法如半解析解、BEM等具有计算速度快的优势，而空气作用明显时（底升角小于10°），基于两相流的VOF法具有较高的数值精度.对于砰击过程中出现的气泡并未建立合理的数学模型，尚不明确其产生、破碎、逃逸等过程，后续研究有必要进一步探讨气泡演变机理.