基于VOF的物體入水-下潛過程數(shù)值模擬

王 平，袁 帥，鄒文峰，陳偉斌，史文奇

（1.國家海洋環(huán)境監(jiān)測中心，大連116023；2.大連中交理工交通技術(shù)研究院，大連116024）

基于VOF的物體入水-下潛過程數(shù)值模擬

王 平1，袁 帥1，鄒文峰2，陳偉斌1，史文奇1

（1.國家海洋環(huán)境監(jiān)測中心，大連116023；2.大連中交理工交通技術(shù)研究院，大連116024）

漂浮物體在取水口附近入水后可能會下潛至涵道，為數(shù)值研究其入水-下潛過程，基于計算流體力學(xué)模型中的VOF法，并結(jié)合動網(wǎng)格方法，首先構(gòu)建了不同密度圓柱體的自由入水以及出水過程，驗證了物體出入水過程造成的射流、入水深度及物體速度變化；其次結(jié)合二維明渠的取水過程模擬了不同初始動量條件下圓柱體的入水-下潛過程。結(jié)果表明：文章建立的入水-出水?dāng)?shù)值模型可很好地刻畫二維圓柱體入水及出水的具體過程，且初始動量大小是圓柱體在入水后能否下潛至取水涵道的關(guān)鍵因素，數(shù)值結(jié)果為判斷取水口附近物體入水后的下潛風(fēng)險提供了分析依據(jù)。

物體入水；下潛；VOF；取水明渠；數(shù)值模擬

取水過程是沿海附近電廠包括核電廠冷卻水循環(huán)的關(guān)鍵步驟，但在取水過程中可能會出現(xiàn)附近漂浮物（如海冰）下潛進入至取水口涵道，進而影響電廠運營的安全。隨著計算機的發(fā)展，數(shù)值模擬物體出入水等問題已越來越成熟，因而采用數(shù)值方法研究取水口附近的物體下潛問題已成為可能，其對核電廠取水及運營安全至關(guān)重要。

數(shù)值計算物體下潛過程的研究較少，與之相近的有隋覺義［1］總結(jié)并分析了冰蓋前緣處的冰塊下潛臨界條件，練繼建等［2］采用VOF模型模擬了渠道冰凌下潛過程，其主要從受力角度分析，并未計算固體的運動過程。而數(shù)值模擬物體入水的研究較多，如Zhao［3］以及Zhu［4］基于勢流理論假設(shè)采用邊界元方法研究了楔形體的入水問題，但無法刻畫入水過程造成的射流影響；Shao［5］采用了無網(wǎng)格的SPH法計算了物體入水中的自由液面及流固相互作用等問題，很好地處理了入水時液面變形較大問題，但數(shù)值計算量較大；Gu［6］采用level set法捕捉自由水面，并模擬了固體在垂向及傾斜入水時的變化過程；Arai［7］采用帶自由表面的兩相流模擬方法（VOF）模擬了二維楔、圓柱船艏的砰擊入水問題，陳宇翔［8］采用了VOF方法對圓柱的入水過程進行了模擬，胡子俊［9］采用基于CIP的VOF法模擬了不同密度圓柱體的入水過程，方舟華［10］則采用該方法模擬了圓柱體入水和出水的整個過程，結(jié)果表明：VOF法可以很好地處理自由面重構(gòu)等強非線性現(xiàn)象。

密度小于水的漂浮物受取水過程的影響一般會聚集在取水口附近水面，但當(dāng)具有一定的初始動量的情況下，則會發(fā)生入水過程，在其入水深度達到一定位置時，則會在取水流速的影響下進入取水涵道，進而造成涵道的堵塞及取水效率的降低。本文基于計算流體力學(xué)模型，首先構(gòu)建了不同密度圓柱體的自由入水及出水過程，研究了物體入水過程造成的射流、入水深度及物體速度變化等；其次模擬了二維明渠的取水過程；最后結(jié)合二維明渠取水及圓柱體入水過程模擬了不同密度圓柱體在不同初始動量條件下的下潛過程。該計算過程為研究取水口物體下潛的可能性提供了一種數(shù)值方法，同時數(shù)值結(jié)果也為判斷取水口附近物體下潛風(fēng)險提供了分析依據(jù)。

1 數(shù)值模型

對于不可壓縮流體，在直角坐標(biāo)系下其運動規(guī)律可以用N-S方程來描述，包括連續(xù)性和動量方程分別為

式中：uj為(x,y,z)3個方向上的流速分量；Fi為單位質(zhì)量力沿(x,y,z)3個方向的分量；P為壓強；υt為紊流粘性系數(shù)，本文采用k-ε湍流閉合模型對其求解。

自由面追蹤采用VOF法，即通過求解單獨的動量方程和處理穿過區(qū)域的每一流體容積比來模擬兩種或多種不能混合的流體。定義體積函數(shù)αq代表第q種流體，對任一單元存在3種可能，即：αq=0表示該單元第q流體不存在；αq=1表示該單元充滿第q流體；0＜αq＜1表示該單元存在著流體之間的界面。

以第q流態(tài)為例，其體積函數(shù)的輸移擴散方程

在任一單元中控制方程的特性參數(shù)是由幾種流體的組合表示，以q1和q2兩種流態(tài)為例，那么單元體的密度可以表示為式（4），其他特性參數(shù)如粘性系數(shù)等均按此方法計算。

對由固體運動造成的動邊界問題則采用動網(wǎng)格技術(shù)，即完成每一步物體的運動狀態(tài)計算后重新劃分網(wǎng)格。為詳細刻畫物體入水時的射流現(xiàn)象在物體周圍設(shè)置移動加密網(wǎng)格，該計算域網(wǎng)格隨物體一起下移。而物體的運動則使用6自由度模型計算，并通過自定義函數(shù)定義其質(zhì)量及轉(zhuǎn)動慣量等。

數(shù)值計算中采用Piso算法耦合速度和壓力場，VOF模型采用自由面重構(gòu)方法來處理液面變化，湍流采用k-ε兩方程模型和壁面函數(shù)法；由于物體入水過程中受到較大的體積力，壓強計算采用體積力加權(quán)（Body Force Weighted）的二階差分格式；密度和動量均采用二階迎風(fēng)格式。iq=1

2 物體入水模擬

圓柱入水是研究入水問題的典型算例，在圓柱入水的瞬間，圓柱的速度突然下降，并改變了周圍壓力等流場參數(shù)，而與其接觸的水會被突然加速，并在圓柱兩側(cè)激起兩股射流，射流斜向上上升并逐漸變得不穩(wěn)定，并在重力作用下逐漸向下彎曲，最終融入水中。若圓柱密度小于水體，圓柱入水后速度逐漸減小，并在達到一定深度后在浮力的作用下逐漸上升，隨后完成一次出水過程。以下則利用數(shù)值方法刻畫圓柱入水及出水的具體過程。

根據(jù)Greenhow［11］的實驗設(shè)計圓柱體的直徑為11 cm，圓柱體的密度分別為零浮力（密度與水相同）和半浮力（密度為水的一半），實驗中圓柱由靜水上方0.5 m處自由下落，忽略空氣阻力，在圓柱接觸自由水面時速度為2.955 m/s。數(shù)值計算中初始時刻為圓柱體剛接觸自由水面，圓柱初始速度設(shè)置為2.955 m/s；計算域選取為1.5 m×1 m的水槽，水槽的上邊界為氣體壓力出口，水槽左右兩側(cè)均為水體壓力出口，并設(shè)定為固定壓力梯度，以減少入水時水位上升帶來的影響，水槽底端為固邊界；計算網(wǎng)格采用三角形網(wǎng)格，網(wǎng)格尺寸在0.35～1.5 cm間變化，圓柱體邊緣采用和圓柱體一起移動的四邊形網(wǎng)格（圖1），計算時間步長為2× 10-4s。

圖1 圓柱入水初始時刻的計算網(wǎng)格Fig.1 Numerical grid for the circular cylinder water?entry case

數(shù)值計算得到零浮力圓柱體在開始下落后0.015 s、0.03 s、0.04 s和0.2 s時以及半浮力圓柱體在開始下落后0.02 s、0.035 s、0.045 s和0.085 s時的自由面變化以及其和實驗結(jié)果的對比見圖2和圖3。從圖中可知本文基于CFD模型所建立的二維圓柱體入水模型很好地模擬了圓柱入水后所引起的射流形成和發(fā)展過程，同時數(shù)值計算得到兩種浮力圓柱體入水后的深度變化與實驗對比見圖4，除個別點異常外，整體結(jié)果吻合較好。

圖2 零浮力圓柱體入水時自由面變化與實驗的對比Fig.2 Free surface deformation of a neutrally buoyant cylinder water?entry

圖3 半浮力圓柱體入水時自由面變化與實驗的對比Fig.3 Free surface deformation of a half buoyant cylinder water?entry

圖4 兩種浮力圓柱體的入水深度與實測值對比Fig.4 Comparison of cylinder underwater depth between simualtion and measurement

為研究物體自由入水和出水的整體過程，Colocchio等［12］在實驗水槽中距水面0.48 m釋放一直徑為0.3 m的不銹鋼圓柱，圓柱與水體密度之比為0.62，由于浮力作用，圓柱經(jīng)歷了入水和出水兩個過程。數(shù)值計算初始時刻圓柱體底部接觸水體，根據(jù)實驗結(jié)果初始速度設(shè)為2.55 m/s。數(shù)值網(wǎng)格在0.1～2 cm之間變化，時間步長為2×10-4s。

構(gòu)建上述類似的數(shù)值水槽，計算域選取為2 m×1.5 m，水槽上邊界為氣體壓力出口，水槽左右兩側(cè)同樣均為水體壓力出口，并設(shè)定為固定壓力梯度，水槽底端為固邊界。數(shù)值計算得到圓柱體底部豎向坐標(biāo)和圓柱豎向速度隨時間變化以及實測值（圖5），從圖可以看出本文建立的圓柱體入水模型與實驗數(shù)據(jù)吻合較好，圓柱體經(jīng)歷了一個明顯的下降、上浮和再下降的過程。

同時給出幾個典型時刻圓柱體周邊的速度場分布（圖6），0.4 s左右圓柱體開始浸沒在水中，同時伴隨著速度很大的頂部射流出現(xiàn)，0.575 s左右射流達到最大高度并開始下落，1.175 s左右柱體則完全浮出水面，此時圓柱體完成了一個完整的入水-出水過程，數(shù)值計算的結(jié)果和實驗時間基本一致。

圖5 圓柱體底部高度和垂向速度變化及和實測值對比Fig.5 Comparison of cylinder bottom height and vertical velocity between simulation and measurement

圖6 圓柱體入水時不同時間的周圍速度場分布Fig.6 Velocity distribution for cylinder water?entry at various water entry moments

3 物體下潛模擬

取水過程會導(dǎo)致取水口附近存在明顯的垂向流速，對于密度接近于水的物體在取水口附近會存在直接下潛的風(fēng)險，而對密度稍小于水的物體其入水過程，并不像在靜止水槽中那樣（圓柱在下潛到一定深度后發(fā)生上浮并完成出水），而是在下潛到一定深度后會受取水流速的影響直接潛入至取水廊道中。以下則利用入水及取水模型數(shù)值研究密度小于水的圓柱在入水后的可能下潛過程。

為獲得取水口附近物體發(fā)生下潛的可能性，首先構(gòu)建一個二維取水的數(shù)值水槽（圖7），水槽尺度19 m×8 m，其中取水涵道距水底1 m，涵道高3 m，為簡化計算涵道長度為3 m。水槽左側(cè)為壓力入水，計算時設(shè)定為固定壓力梯度，即水位保持固定在6 m；水槽右端為速度出口，計算時設(shè)定為穩(wěn)定流量2.4 m3/s；水槽內(nèi)水面高度為6 m，上邊界為氣體的壓力出口；圓柱體直徑為0.3 m，初始位置設(shè)置在（x=13.5 m，y= 6.15 m）處，圓柱體與水體密度之比為0.92。

數(shù)值計算初始階段圓柱體固定不動，只計算水槽內(nèi)的流場變化。由于明渠內(nèi)的水體由靜止開始計算，初期的自由液面會有一定的波動，隨著取水的進行，波動逐漸減小并將趨于穩(wěn)定，圖8給出了x=14 m處自由液面隨時間的變化，從圖可知在單獨計算明渠取水1 000 s左右后，水槽內(nèi)的水面基本保持不變，此時開始計算圓柱體的入水過程。

圖9給出了圓柱體開始入水時取水明渠內(nèi)的流速矢量分布，從圖中可知，受取水影響取水口附近水體存在明顯的垂向流速，其將決定圓柱體是否會潛入涵道。忽略空氣影響，考慮兩種工況，工況1：假設(shè)圓柱從距水面0.8 m的高度下落，圓柱體接觸水面時的速度約4 m/s；工況2：假設(shè)圓柱從距水面0.45 m的高度下落，圓柱體接觸水面時的速度約3 m/s，數(shù)值計算均從圓柱體接觸自由液面開始。

數(shù)值計算在圓柱體與固邊界發(fā)生碰撞溢出后停止，計算得到工況1條件下圓柱體的入水-下潛過程，見圖10-a，圖中給出了入水后每隔0.5 s圓柱體的位置變化，以及入水開始時明渠內(nèi)的液面及流線圖；工況2中圓柱體在入水5 s后與取水口前壁發(fā)生碰撞，此時數(shù)值計算停止，圖10-b給出了此刻明渠內(nèi)的液面及流線圖，以及入水后每隔0.5 s的圓柱體位置變化。

從圖10可知，工況1條件下的圓柱體在入水后下潛至取水涵道內(nèi)，而工況2由于與取水口前壁碰撞并未直接潛入取水涵道，因為初始條件的不同，兩種工況下圓柱體的入水路徑有所區(qū)別，且工況1圓柱體潛入涵道的時間要短于工況2，說明工況1條件下的圓柱體在入水后更易發(fā)生下潛。

圖7 二維取水明渠及圓柱體下潛計算示意Fig.7 Calculation sketch of two?dimensional water channel and cylinder dive

圖8 明渠取水計算時x=14 m處的水位變化Fig.8 Surface change inx=14 m of open channel

圖9 圓柱體開始入水時的取水明渠內(nèi)的流場分布Fig.9 Velocity distribution of open channel at water entry moment

圖10 兩種工況下圓柱體的入水-下潛過程Fig.10 Cylinder water entry and dive process of two cases simulated by present model

為對比取水對圓柱體下潛過程的影響，同時模擬了在沒有取水（靜水）時，工況1和2條件下的圓柱體入水-出水過程，數(shù)值得到工況1、工況2及其靜水條件下圓柱體的入水深度和y向速度變化對比如圖11。從圖可知，工況1（靜水）時，圓柱體在入水3 s后開始上浮，而在受到取水流速的影響后，圓柱體則會繼續(xù)向下移動，并在4 s后潛入涵道，其y向速度值也在3 s后不斷增加并在潛入至涵道后逐漸減小；工況2（靜水）條件下，圓柱在入水2.5 s后開始上浮，而受取水影響后，圓柱體的下潛深度及y向速度在入水4 s后均逐漸增大，并在入水后5 s時與取水口前壁發(fā)生碰撞。

圖11 不同情況下圓柱體的入水深度及y向速度變化對比Fig.11 Comparison of cylinder bottom height and vertical velocity between different cases

4 結(jié)論

為數(shù)值研究取水口附近物體下潛的可能性，基于計算流體力學(xué)模型中的VOF法，并結(jié)合動網(wǎng)格方法構(gòu)建了不同密度圓柱體的自由入水及出水過程，驗證了物體入水過程造成的射流、入水深度及物體速度變化等；其次模擬了二維明渠的取水過程；最后結(jié)合二維明渠取水及圓柱體入水模擬了不同初始動量條件下圓柱體的入水-下潛過程。

數(shù)值計算結(jié)果表明：本文建立的入水-出水模型可很好地模擬二維圓柱體入水及出水形成的射流及速度變化，且初始動量大小決定了圓柱體在入水后是否能下潛至取水涵道，該計算過程為研究取水口附近物體入水后的下潛可能性提供了一種數(shù)值方法。由于計算量的限制，本文并未考慮三維實際情況，后續(xù)研究應(yīng)予以關(guān)注。

［1］隋覺義，方達憲.冰蓋前緣處冰塊下潛臨界條件研究［J］.水利學(xué)報，1993（8）：10-18. SUI J Y，F(xiàn)ANG D X.Study on the critical condition of ice dive in ice cover leading edge［J］.Journal of Hydraulic Engineering，1993（8）：10-18.

［2］練繼建，羅昉昕，趙新.渠道冰凌下潛的數(shù)值模擬研究［J］.水利學(xué)報，2013，44（6）：673-679. LIAN J J，LUO F X，ZHAO X.Numerical simulation of ice submergence in water conveyance channel［J］.Journal of Hydraulic En?gineering，2013，44（6）：673-679.

［3］Zhao R，F(xiàn)altinsen O M.Water?entry of two?dimensional bodies［J］.Journal of Fluid Mechanics，1993，246：593-612.

［4］Zhu X Y，F(xiàn)altinsen O M，Hu C H.Water entry and exit of a horizontal circular cylinder［J］.Journal of Offshore Mechanics and Arc?tic Engineering，2007，129：253-264.

［5］Shao S D.Incompressible SPH simulation of water entry of a free?falling object［J］.Int.J.Numer.Meth.Fluids，2009，59：91-115.

［6］Gu H B，Qian L，Causon D M，et al.Numerical simulation of water impact of solid bodies with vertical and oblique entries［J］. Ocean engineering，2013，75（5）：128-137.

［7］Arai M，Chen L Y，Inoue Y.A computing method for the analysis of water impact of arbitrary shaped bodies［J］.Journal of Computa?tional Physics，1981，39：201-225.

［8］陳宇翔，郜冶，劉乾坤.應(yīng)用VOF方法的水平圓柱入水?dāng)?shù)值模擬［J］.哈爾濱工程大學(xué)學(xué)報，2011，11（32）：1 439-1 442. CHEN Y X，GAO Y，LIU Q K.Numerical simulation of water?entry in a horizontal circular cylinder using the volume of fluid（VOF）method［J］.Journal of Harbin Engineering University，2011，11（32）：1 439-1 442.

［9］胡子俊，葉洲騰，趙西增，等.基于CIP方法的水平圓柱入水?dāng)?shù)值模擬［C］//左其華，竇希萍.第17屆中國海洋（岸）工程學(xué)術(shù)討論會論文集.南寧：海洋出版社，2015：117-181.

［10］方舟華，趙西增.物體入水的CIP方法模擬［C］//左其華，竇希萍.第17屆中國海洋（岸）工程學(xué)術(shù)討論會論文集.南寧：海洋出版社，2015：289-294.

［11］Greenhow M，Lin W M.Nonlinear?free surface effects：experiments and theory［R］.Cambridge：Massachusetts Inst of Tech Cam?bridge Dept of Ocean Engineering，1983.

［12］Colicchio G，Greco M，Miozzi M，et al.Experimental and numerical investigation of the water?entry and water?exit of a circular cyl?inder［C］//Proceedings of the 24th Int.Workshop on Water Waves and Floating Bodies.Zelenogorsk，Russia，2009：19-22.

Numerical method and study of the body water?entry and dive based on VOF

WANG Ping1,YUAN Shuai1,ZOU Wen?feng2,CHEN Wei?bin1,SHI Wen?qi1
（1.National Marine Environmental Monitoring Center,Dalian 116023,China;2.CCCC&DLUT Institute of Transportation Technology Co.Ltd.,Dalian 116024,China）

The body would dive into duct near water intake when it was entering into the water.To numerically study the body water?entry and dive,combined the volume of fluid(VOF)method with dynamic mesh,the cylinder water?entry and water?out with different density was simulated by the Fluent module.The free surface deformation, bottom height and vertical velocity caused by cylinder water?entry and water?out were verified.The cases of cylinder water?entry and dive into water intake with different initial momentum conditions were researched in case of water?entry and water?out combining with two?dimensional water channel.Numerical results show that the specific pro?cess of cylinder water?entry and water?out could be simulated accurately by the module,and initial momentum is the key factor on the cylinder diving into duct after entering into the water.Numerical results provide the analysis basis for judging whether the body dives into water intake.

water?entry;dive;the volume of fluid(VOF)method;water channel;numerical simulation

TV 131；O 242.1

A

1005-8443（2017）02-0120-06

2016-09-19；

2016-10-25

海洋公益性行業(yè)科研專項（200905007；201505019）

王平（1988-），男，安徽省人，助理研究員，主要從事近岸動力學(xué)研究。

Biography：WANG Ping（1988-），male，assistant professor.