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

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

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

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

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

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

1 數(shù)值模型及數(shù)值過程

1.1 數(shù)值模型

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

式中：U為流場(chǎng)速度；t為時(shí)間；p為壓力；ρ為流體密度；μ為黏度；α1為氣體的體積分?jǐn)?shù)（基本相），是隨時(shí)間變化的函數(shù).在數(shù)值計(jì)算中，當(dāng)α1＝1時(shí)表示該網(wǎng)格內(nèi)為空氣，當(dāng)α1＝0時(shí)表示網(wǎng)格內(nèi)為水，當(dāng)0＜α1＜1時(shí)，表示網(wǎng)格內(nèi)為水和空氣的交界面.為保證相分?jǐn)?shù)嚴(yán)格介于[0，1]，可采用高分辨率的離散模式如 HRIC或CICSAM求解式（3），但這些離散模式在時(shí)間步進(jìn)過程中會(huì)出現(xiàn)數(shù)值耗散，導(dǎo)致相界面的模糊性.為此，可通過添加人工對(duì)流項(xiàng)對(duì)相界面附近的體積分?jǐn)?shù)進(jìn)行擠壓來緩解界面的模糊性.根據(jù)兩相流理論，假設(shè)兩相的相對(duì)速度Ur＝Ul－U2，添加對(duì)流項(xiàng)將式（3）轉(zhuǎn)化為

式中：Δ·（α1α2Ur）為人工對(duì)流項(xiàng)；α2為水的體積分?jǐn)?shù)，α2＝1－α1.通過體積分?jǐn)?shù)，將流體密度ρ和黏度μ更新為

式中：ρ1和ρ2分別為空氣和水的密度；μ1和μ2分別為空氣和水的黏度.由于VOF法中并不存在所謂的多個(gè)相速度，故將相對(duì)速度采用單向流的速度進(jìn)行?；?，即

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

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

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

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

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

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

VOF法中由于界面不連續(xù)性的影響難以準(zhǔn)確計(jì)算界面的法向和曲率，為了精確捕捉自由面，采用Level－set函數(shù)來保持界面的光滑性.通過Level－set函數(shù)φ（x，y，t）定義 VOF的體積函數(shù)F（α1，t）來實(shí)現(xiàn) VOF和Level－set的耦合（CLSVOF），則有

式中：Ω為網(wǎng)格域；H為階躍函數(shù)，記為

這種界面體積分?jǐn)?shù)的光順可在Fluent軟件中通過VOF模型求解設(shè)置中勾選耦合CLSVOF模型直接實(shí)現(xiàn).

1.2 數(shù)值過程

當(dāng)物體以給定的入水速度下落時(shí)，在每個(gè)時(shí)刻t，具體的數(shù)值求解流程如下：① 取cα為0.5，對(duì)相方程采用有限差分法求解相體積分?jǐn)?shù)α1，更新流體密度ρ和黏度μ；② 對(duì)式（11）進(jìn)行有限體積離散，修正項(xiàng)通過添加表面張力系數(shù)來實(shí)現(xiàn)源項(xiàng)的添加，時(shí)間項(xiàng)采用歐拉隱式離散，壓力項(xiàng)采用PRESTO！，速度項(xiàng)采用二階迎風(fēng)格式，迭代求解離散后的非線性方程組，從而給出速度的預(yù)測(cè)值；③ 結(jié)合式（1）采用PISO算法對(duì)流域壓力場(chǎng)和速度場(chǎng)進(jìn)行修正[8].應(yīng)用動(dòng)網(wǎng)格技術(shù)，根據(jù)給定的入水速度，預(yù)測(cè)下一時(shí)刻物體的位置，并采用局部重構(gòu)方法對(duì)流域網(wǎng)格進(jìn)行重新劃分.

2 收斂性研究

選取文獻(xiàn)[9]中的30°楔形體進(jìn)行網(wǎng)格收斂性的研究，初始入水速度為2.42m／s.以3種不同的網(wǎng)格尺寸對(duì)流域進(jìn)行劃分，記為Mesh1、Mesh2和Mesh3，具體網(wǎng)格尺寸信息和計(jì)算時(shí)間見表1.圖2給出了不同網(wǎng)格尺寸下數(shù)值結(jié)果與文獻(xiàn)[9－10]的BEM的對(duì)比情況.當(dāng)網(wǎng)格尺寸為Mesh1時(shí)，最大垂向力比實(shí)驗(yàn)值小12%；當(dāng)網(wǎng)格尺寸減小到Mesh2時(shí)，結(jié)果偏大7%，主要原因是由于實(shí)驗(yàn)的三維效應(yīng)，在文獻(xiàn)[9]中，二維預(yù)測(cè)的最大垂向力比實(shí)驗(yàn)值大8%.當(dāng)網(wǎng)格尺寸進(jìn)一步減至Mesh3時(shí)，砰擊力趨于穩(wěn)定，這表明當(dāng)網(wǎng)格尺寸足夠小時(shí)，數(shù)值解趨于一致.對(duì)于砰擊壓力的結(jié)果（壓力位置距低端尖角的垂向距離為45mm），雖然實(shí)驗(yàn)中并未給出，但從文獻(xiàn)[10]結(jié)果的對(duì)比可以看出，數(shù)值的精度滿足要求.根據(jù)垂向力和砰擊壓力的對(duì)比結(jié)果，Mesh2這一網(wǎng)格尺寸可以保證預(yù)測(cè)的精度.為了節(jié)省計(jì)算時(shí)間，本文的流域網(wǎng)格生成均采用Mesh2中的網(wǎng)格尺寸參數(shù).

表1 用于收斂性分析的不同網(wǎng)格尺寸Tab.1 Different mesh sizes in mesh convergence analysis

圖2 網(wǎng)格的收斂性分析Fig.2 Mesh convergence analysis

3 數(shù)值結(jié)果

3.1 研究對(duì)象

本研究以楔形體為研究對(duì)象，入水砰擊模型如圖3所示.為了便于對(duì)比分析，選取文獻(xiàn)[11]中楔形體進(jìn)行建模.底升角β＝20°，傾斜角為θ，計(jì)算中分別選為5°、10°和15°.R為迎風(fēng)側(cè)，底升角為β－θ；L為背風(fēng)側(cè)，底升角為β＋θ；F為氣－液交界面.模型總質(zhì)量為44kg，單邊長度為0.205m，初始入水速度v＝3.13m／s.實(shí)驗(yàn)過程中每側(cè)布置4個(gè)壓力測(cè)點(diǎn)，具體位置如圖4所示.整個(gè)流域尺寸為4m×3m，初始水深為2m，近楔形體處局部網(wǎng)格如圖5所示.入水過程中忽略重力的影響，t＝0表示剖面底部剛好接觸靜水面.

圖3 非對(duì)稱入水砰擊模型示意圖Fig.3 Schematic diagram of asymmetrical water entry model

圖4 實(shí)驗(yàn)中壓力測(cè)點(diǎn)布置示意圖[11]Fig.4 Schematic diagram of pressure measured positions in the experiment[11]

圖5 近物體表面網(wǎng)格示意圖Fig.5 Schematic diagram of mesh near the body

3.2 改進(jìn)方法的影響分析

為探討改進(jìn)后的VOF法對(duì)自由面重構(gòu)的分辨率及壓力分布的影響，首先需要驗(yàn)證VOF法對(duì)氣－液交界面捕捉的精確度.由于文獻(xiàn)[11]實(shí)驗(yàn)中并未給出液面抬升的實(shí)驗(yàn)圖片，本研究選取文獻(xiàn)[12]中給出的自由面進(jìn)行對(duì)比分析.數(shù)值模擬中垂向入水速度為3.81m／s，橫向速度為7.62m／s，楔形體底升角為37°，傾斜角為5°.圖6給出了不同時(shí)刻下實(shí)驗(yàn)和數(shù)值模擬的自由面的對(duì)比情況，可以看出VOF法能夠較好地反映入水過程中的氣－液交界面.

圖7給出了楔形體斜向入水過程中分別采用改進(jìn)前（左側(cè)）和改進(jìn)后（右側(cè)）VOF法獲得自由面的對(duì)比情況.模擬過程中網(wǎng)格的總量大約為3.7×104，改進(jìn)前和改進(jìn)后的計(jì)算時(shí)間分別為2.6和3.0 h.可以看出改進(jìn)后所捕捉到的自由面厚度層明顯減小，而且輪廓更加尖銳.當(dāng)傾斜角增大時(shí)（圖7（b）），左側(cè)的自由面更加模糊，而右側(cè)仍能保持較好的清晰度.這是由于傾斜角越大，近物體處的自由面上相關(guān)變量（如密度、黏度和壓力）的梯度越大，僅采用插值方式（改進(jìn)前）對(duì)體積分?jǐn)?shù)求解存在較大誤差，而CSF模型的引入保證了交界面處的密度和黏度充分光順.在改進(jìn)前的方法中，自由面的模糊性可通過增加網(wǎng)格密度進(jìn)行改善，比如當(dāng)網(wǎng)格數(shù)量達(dá)到4.9×104時(shí)（增加大約1／3），自由面清晰度可以達(dá)到與改進(jìn)后的方法相同的效果，但這種措施不可避免地增加了計(jì)算的時(shí)間.

圖6 實(shí)驗(yàn)和數(shù)值模擬的自由面對(duì)比Fig.6 Comparison of free surface obtained by experimental and numerical results

圖7 改進(jìn)前后自由面形狀對(duì)比Fig.7 Comparison of free surface before and after improvement

圖8給出了不同數(shù)值方法和砰擊壓力沿物體表面分布的對(duì)比情況（P2峰值時(shí)刻）.圖中：無因次化量Cp＝2ps／（ρ2v2），ps為砰擊壓力；x′為楔形體的橫向長度x與浸濕半寬CL的比值，即x′＝x／CL；緊致插值曲線（CIP）法為文獻(xiàn)[13]的結(jié)果，半解析解為文獻(xiàn)[14]的結(jié)果，實(shí)驗(yàn)結(jié)果來自于文獻(xiàn)[11].改進(jìn)前的VOF法盡管在峰值上與修正的結(jié)果差異不大，但在峰值之間的部分有所偏大，這種偏差在一定程度上降低了兩側(cè)壓力的不對(duì)稱性，進(jìn)而會(huì)導(dǎo)致扭轉(zhuǎn)力矩的減小.另外在圖8（a）中半解析解的結(jié)果比實(shí)驗(yàn)值偏大18%，而這種誤差在圖8（b）中增加到40%，相比而言，CIP法和改進(jìn)VOF法的結(jié)果始終維持在10%左右.由圖7（b）中自由面可見，入水過程中包裹了一個(gè)氣泡，會(huì)導(dǎo)致氣墊效應(yīng)，進(jìn)而減緩砰擊壓力.由于半解析解是基于單向流模型，氣體效應(yīng)被忽略，而CIP法和VOF法是基于氣－液的兩相流模型，合理地計(jì)及了氣體的影響，這也正是它們之間存在較大誤差的原因所在.

圖8 改進(jìn)前后壓力分布對(duì)比Fig.8 Comparison of pressure distribution before and after improvement

3.3 砰擊載荷特性分析

圖9給出了不同θ下迎風(fēng)側(cè)的壓力時(shí)歷（文獻(xiàn)[11]實(shí)驗(yàn)中P1點(diǎn)），其中對(duì)稱結(jié)果表示兩側(cè)底升角均為β－θ的對(duì)稱楔形體的數(shù)值結(jié)果.為了清晰地對(duì)比結(jié)果，將計(jì)算的非對(duì)稱結(jié)果后移2ms，對(duì)稱結(jié)果前移2ms，實(shí)驗(yàn)結(jié)果保持不變.

由圖9可以看出：非對(duì)稱結(jié)果與實(shí)驗(yàn)值相當(dāng)接近（誤差小于5%），說明了本文預(yù)測(cè)結(jié)果的合理性；在圖9（a）中對(duì)稱結(jié)果比實(shí)驗(yàn)值偏大15%，這種差異性在底升角更小的圖9（b）中達(dá)到了28%，這主要是由于斜向入水的過程中，背風(fēng)側(cè)的底升角β＋θ較大，對(duì)流體的阻礙作用要小于兩側(cè)底升角都為β－θ的情況.研究結(jié)果表明：對(duì)于底升角β的楔形體以傾斜角θ入水時(shí)，如果單純地采用底升角為β－θ情況下的對(duì)稱楔形體進(jìn)行估算，砰擊壓力會(huì)造成結(jié)果的偏大；對(duì)稱結(jié)果的壓力上升速度要緩于實(shí)驗(yàn)值，相比于對(duì)稱的情況，非對(duì)稱情況下迎風(fēng)側(cè)流體的爬升速度更快些.

圖9 迎風(fēng)側(cè)的砰擊壓力時(shí)歷Fig.9 Pressure time histories on the upwind side

圖10給出了入水過程中兩側(cè)砰擊壓力的分布情況.隨著時(shí)間的增加（t≤14.4ms），迎風(fēng)側(cè)的壓力開始高于背風(fēng)側(cè)的壓力，兩側(cè)壓力峰值的最大相對(duì)誤差達(dá)到97%（t＝10.8ms）.之后迎風(fēng)側(cè)壓力減小，背風(fēng)側(cè)壓力仍在上升，在t＝21.6ms時(shí)超過迎風(fēng)側(cè)壓力，最后兩側(cè)壓力趨于對(duì)稱化（t＝25.2ms）.事實(shí)上這種復(fù)雜的變化主要是由于流體從折角處分離造成的.在t＝14.4ms時(shí)迎風(fēng)側(cè)浸入水中，導(dǎo)致該側(cè)壓力開始減小.而此時(shí)背風(fēng)側(cè)尚未被流體完全淹沒，壓力仍在增加，直到t＝21.6ms左右，背風(fēng)側(cè)浸入水中，該側(cè)壓力也開始減小.當(dāng)物體完全浸入水中后，由于物體處流體的波動(dòng)開始減小，此時(shí)的傾斜角對(duì)兩側(cè)壓力分布的影響可以忽略（t＝25.2ms）.圖10（b）中出現(xiàn)了與圖10（a）相似的規(guī)律性，但不同的是兩側(cè)的差異性達(dá)到307%（t＝7.2ms）.需要注意的是，底升角較小時(shí)迎風(fēng)側(cè)包裹了一個(gè)氣泡，如圖11所示，此時(shí)物面上的壓力仍是連續(xù)的，這也就意味著楔形體與空氣接觸的部位的壓力并非為0.由于此時(shí)流體會(huì)對(duì)氣體域產(chǎn)生擠壓，而楔形體作為空氣域的一個(gè)接觸部位，自然會(huì)承受這種擠壓作用，從而出現(xiàn)了大于大氣壓的壓力值.這種氣穴效應(yīng)，如果采用單向流模型，由于自由面邊界條件的限制，水和空氣的交界處壓力為0（大氣壓），會(huì)導(dǎo)致與空氣接觸部位的壓力為0.這種現(xiàn)象在文獻(xiàn)[15]的數(shù)值與實(shí)驗(yàn)的對(duì)比結(jié)果中也有所體現(xiàn)，在數(shù)值結(jié)果中氣穴處始終為0，而實(shí)驗(yàn)中卻出現(xiàn)了明顯的砰擊壓力.

圖10 不同時(shí)刻下兩側(cè)壓力分布Fig.10 Pressure distribution at different instants

圖11 t＝21.6ms時(shí)的自由面形狀Fig.11 Free surface at the instant of t＝21.6ms

4 結(jié)論

應(yīng)用改進(jìn)VOF法研究了楔形體斜向入水問題，得出以下結(jié)論：

（1）改進(jìn)VOF法能夠有效保持自由面的尖銳特征，傾斜角較大時(shí)，這種優(yōu)勢(shì)比改進(jìn)前更加明顯.

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

（3）當(dāng)?shù)咨切∮?0°時(shí)，空氣對(duì)砰擊壓力的影響顯著.計(jì)及空氣效應(yīng)的結(jié)果與實(shí)驗(yàn)的誤差控制在10%以內(nèi)，而未考慮空氣影響的結(jié)果誤差高達(dá)40%.

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