謝 行， 任慧龍， 陶凱東， 馮億坤

（1.哈爾濱工程大學 船舶工程學院，哈爾濱150001；2.黃淮學院 智能制造學院，河南 駐馬店463000）

典型的入水沖擊包含流體爬升、流動的分離、破碎、空泡的形成與發(fā)展等階段，涉及氣－液－固三者的相互作用.在此過程中伴隨著自由面的大變形，甚至會出現流體的翻卷、破碎、氣穴等復雜物理現象，對這些現象進行精確模擬需要極高的數值精度.

由于入水沖擊的復雜性，實際數值模擬中對其進行了一些簡化，如剖面形式簡化為楔形體，流體性質不考慮黏性，忽略空氣效應.Zhao等[1]采用邊界元法（BEM）對二維有限寬度楔形體的入水問題進行了研究，通過對射流進行合理的截斷，避免了時間步進中數值的發(fā)散.Qin等[2]基于勢流理論，將自由面和物面條件進行線性化，推導了求解速度勢的半解析公式，通過修正伯努利方程，研究了非對稱楔形體的入水問題，結果與Xu等[3]的全非線性方法相符.上述均為單向流模型的方法，對入水后期流體的破碎、氣穴等難以進行合理模擬.

為了考慮空氣的效應，基于兩相流模型的流體體積（VOF）算法被加以應用.Krastev等[4]通過VOF模型求解兩相流的雷諾平均（RANS）方程，合理地模擬了楔形體非對稱入水過程中氣泡的產生和演變.Xie等[5]基于Fluent軟件對具有傾斜角的外飄剖面進行了數值模擬，采用VOF法追蹤自由面，模擬過程中流體先發(fā)生了分離，之后撞擊到物體表面，發(fā)生了二次砰擊現象.王易君等[6]應用VOF模型和動網格技術實現了平底結構的入水模擬，入水過程中出現了明顯的空氣墊.在現有VOF法中，相方程不但在離散求解過程中存在數值耗散現象，而且通過插值方法對氣－液交界面處壓力的不連續(xù)性進行處理時，也會造成相應的數值誤差.對于自由面的重構，通過VOF法求得的體積分數存在一個從0到1的過渡域，這導致重構的相界面具有明顯的厚度層，難以準確地計算界面的法向和曲率.

本文基于 Navier－Stokes（N－S）方程建立非對稱入水問題的氣－液兩相流模型，在Fluent軟件中，通過自定義函數（UDF），引入壓力連續(xù)函數和人工對流項，并耦合水平集（Level－set）函數對現有VOF法進行改進.應用改進的VOF法對楔形體的斜向入水問題進行數值模擬，并對入水過程中砰擊壓力和自由面特性進行分析.

1 數值模型及數值過程

1.1 數值模型

VOF法將氣－液兩相流視為密度可變的單項流模型，利用體積分數來確定自由面的位置和形狀.假定氣液不相容，其連續(xù)性方程、動量方程和相方程分別為

式中：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法中并不存在所謂的多個相速度，故將相對速度采用單向流的速度進行?；?/p>

式中：cα為相界面壓縮因子，cα＝0表示無壓縮效果，cα越大壓縮效應越明顯.另外，Ur的主要作用在于壓縮相界面來保持尖銳性，應該在相界面的法向而不是切向上進行壓縮，否則會引起虛假的擴散，因此Ur的方向與界面法向相Δ同.在兩相流的交Δ界處，引入某一點的法向量記為α，如圖1所示在水和空氣區(qū)域內為0，在自由表面處不為0.此時交界面的單位法向量n可表示為

圖1 交界處示意圖Fig.1 Schematic diagram of interface

為了避免氣－液交界面兩側壓力的不連續(xù)性，引入Brackbill等[7]提出的連續(xù)表面張力（CSF）模型，在界面處定義一個連續(xù)函數.將過渡域內的壓力定義為一個連續(xù)函數

式中：c為界面處的位置函數；C為表面張力系數.則整個區(qū)域的壓力可修正為

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

式（1）、（4）和（11）構成了改進 VOF法的連續(xù)性方程、相方程和動量方程.可以看出式（11）比式（2）多了一個源項κCΔα.

VOF法中由于界面不連續(xù)性的影響難以準確計算界面的法向和曲率，為了精確捕捉自由面，采用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時，砰擊力趨于穩(wěn)定，這表明當網格尺寸足夠小時，數值解趨于一致.對于砰擊壓力的結果（壓力位置距低端尖角的垂向距離為45mm），雖然實驗中并未給出，但從文獻[10]結果的對比可以看出，數值的精度滿足要求.根據垂向力和砰擊壓力的對比結果，Mesh2這一網格尺寸可以保證預測的精度.為了節(jié)省計算時間，本文的流域網格生成均采用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）相似的規(guī)律性，但不同的是兩側的差異性達到307%（t＝7.2ms）.需要注意的是，底升角較小時迎風側包裹了一個氣泡，如圖11所示，此時物面上的壓力仍是連續(xù)的，這也就意味著楔形體與空氣接觸的部位的壓力并非為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法能夠有效保持自由面的尖銳特征，傾斜角較大時，這種優(yōu)勢比改進前更加明顯.

（2）隨著傾斜角的增大，迎風側砰擊壓力迅速增加，兩側壓力的不對稱性更加明顯，而對稱情況下的砰擊壓力偏大20%左右.

（3）當底升角小于10°時，空氣對砰擊壓力的影響顯著.計及空氣效應的結果與實驗的誤差控制在10%以內，而未考慮空氣影響的結果誤差高達40%.

當砰擊過程中的空氣效應不明顯時（底升角大于10°），一些基于單相流模型的方法如半解析解、BEM等具有計算速度快的優(yōu)勢，而空氣作用明顯時（底升角小于10°），基于兩相流的VOF法具有較高的數值精度.對于砰擊過程中出現的氣泡并未建立合理的數學模型，尚不明確其產生、破碎、逃逸等過程，后續(xù)研究有必要進一步探討氣泡演變機理.