張耐民，裴金亮，趙 陽，尤天慶，劉元清

?

有限水層對(duì)物體運(yùn)動(dòng)影響的數(shù)值計(jì)算研究

張耐民1,2，裴金亮2，趙 陽1，尤天慶2，劉元清2

（1. 哈爾濱工業(yè)大學(xué)，哈爾濱，150001；2. 北京宇航系統(tǒng)工程研究所，北京，100076）

穿越自由液面過程中航行體表面攜帶的附著水層會(huì)對(duì)分離物體受力及分離運(yùn)動(dòng)產(chǎn)生重要影響。首先采用流體體積方法（Volume of Fluid，VOF）建立了圓柱運(yùn)動(dòng)的剛體動(dòng)力學(xué)方程和氣液兩相流動(dòng)N-S方程的耦合求解模型，基于試驗(yàn)數(shù)據(jù)驗(yàn)證了算法的準(zhǔn)確性。在此基礎(chǔ)上，開展了不同外形物體穿越有限水層的計(jì)算分析。研究了不同形狀、結(jié)構(gòu)質(zhì)量、運(yùn)動(dòng)形式、水層厚度等因素對(duì)物體穿越有限水層的受力情況和運(yùn)動(dòng)變化規(guī)律。

附著水層；流體體積方法；N-S方程；航行體

0 引 言

水下高速運(yùn)動(dòng)的航行體出水過程涉及眾多物理現(xiàn)象，穿越自由液面過程中航行體表面攜帶附著水是眾多復(fù)雜物理現(xiàn)象之一，與航行體的運(yùn)動(dòng)狀態(tài)、空泡狀態(tài)等密切相關(guān)。航行體出水后進(jìn)行的尾罩分離等運(yùn)動(dòng)即在附著水層的覆蓋下進(jìn)行，附著水層會(huì)對(duì)分離物體的受力及分離運(yùn)動(dòng)產(chǎn)生重要影響。針對(duì)此問題的物理及數(shù)學(xué)模型都是非常復(fù)雜的。隨著計(jì)算流體動(dòng)力學(xué)（Computational Fluid Dynamics，CFD）的發(fā)展，如流體體積法（Volume of Fluid，VOF）等能比較好地解決自由液面問題。Hirt等[1]提出VOF方法以來，帶自由表面的兩相粘性不可壓縮流體運(yùn)動(dòng)的數(shù)值計(jì)算技術(shù)得到了迅速的發(fā)展。近30年來，VOF方法在物體入水的數(shù)值計(jì)算領(lǐng)域得到了廣泛的應(yīng)用。Arai等[2]假設(shè)流動(dòng)為無粘不可壓流動(dòng)，應(yīng)用VOF方法對(duì)二維楔、圓柱、船艏的砰擊入水問題進(jìn)行了模擬；Schumann[3]用類似的方法模擬了船艏的砰擊入水；Xing等[4]基于VOF方法，應(yīng)用Comet軟件對(duì)二維圓柱的粘性不可壓縮出入水運(yùn)動(dòng)進(jìn)行了模擬，取得了與試驗(yàn)數(shù)據(jù)符合較好的結(jié)果。

本文首先應(yīng)用VOF方法實(shí)現(xiàn)圓柱運(yùn)動(dòng)的剛體動(dòng)力學(xué)方程和氣液兩相流動(dòng)N-S方程的耦合求解，基于試驗(yàn)數(shù)據(jù)驗(yàn)證了算法的準(zhǔn)確性。在此基礎(chǔ)上，分析了不同外形物體穿越有限水層的受力和運(yùn)動(dòng)變化規(guī)律。

1 VOF模型及基本控制方程組

連續(xù)性方程：

混合物的動(dòng)量方程：

混合物的能量方程：

對(duì)任意向量，發(fā)射坐標(biāo)系和體坐標(biāo)系之間的轉(zhuǎn)換關(guān)系為

剛體的平動(dòng)和轉(zhuǎn)動(dòng)的控制方程為

2 數(shù)值計(jì)算方法有效性驗(yàn)證

為驗(yàn)證此數(shù)值計(jì)算方法在計(jì)算物體出水過程中的有效性，以直徑=11 cm的中性浮力圓柱體為研究對(duì)象，當(dāng)其漂浮于水下一定深度時(shí)額外施加大小等于自重的外力，使其向上運(yùn)動(dòng)直至沖出水面。計(jì)算域?yàn)閷?0、高8的方形區(qū)域。

計(jì)算結(jié)果能很好地反映圓柱體出水的主要現(xiàn)象。圓柱體上方的水由于圓柱的向上運(yùn)動(dòng)而被抬升，進(jìn)而形成了一層水附著在圓柱表面。此時(shí)由于圓柱逐漸脫離水的表面，浮力逐漸減小，導(dǎo)致加速度有一定幅度的減小。隨著圓柱的繼續(xù)向上運(yùn)動(dòng)，附著在圓柱表面的水逐漸下落至水面，并導(dǎo)致自由表面的波動(dòng)。中性浮力圓柱出水距離的變化及出水時(shí)自由液面變化情況如圖1和圖2所示，數(shù)值計(jì)算結(jié)果與試驗(yàn)吻合較好。

圖1 中性浮力圓柱出水距離的變化曲線

圖2 中性浮力圓柱出水時(shí)自由液面變化情況

注： 右側(cè)為VOF數(shù)值模擬結(jié)果，左側(cè)為Greenhow & Lin的試驗(yàn)照片[9]

3 算例及分析

采用上述方法，計(jì)算物體穿越有限水層問題，如圖3所示。計(jì)算過程考慮水的粘性，同時(shí)為簡(jiǎn)化問題，忽略水的重力影響。

圖3 物體穿越有限水層示意

物體最大截面直徑；最大橫截面積；水層厚度；物體的運(yùn)動(dòng)速度

3.1 物體質(zhì)量的影響

穿越水層物體的質(zhì)量對(duì)運(yùn)動(dòng)形式有著重要影響。在無其它外力條件下，運(yùn)動(dòng)質(zhì)量塊與靜止質(zhì)量塊發(fā)生彈性碰撞時(shí)，若靜止質(zhì)量塊質(zhì)量大于運(yùn)動(dòng)質(zhì)量塊，則運(yùn)動(dòng)的物體將會(huì)被彈回。針對(duì)穿越水層的物體也是如此，所不同的是，此處發(fā)生的不完全是彈性碰撞，穿越水層物體的運(yùn)動(dòng)形式除受其自身質(zhì)量的影響外，還與其外形有著重要的關(guān)系。

采用以上方法計(jì)算了尺寸為、質(zhì)量為的平板與圓柱，以指定初始速度0穿越水層，如圖4所示。對(duì)于平板，由于其與水層接觸面積較大，在拍水過程中，接觸面產(chǎn)生相對(duì)較大的反作用力，當(dāng)=0.25時(shí)平板被反彈，水層由于拍擊力作用，呈現(xiàn)噴濺運(yùn)動(dòng)趨勢(shì)。當(dāng)平板質(zhì)量增加至=1.00時(shí)，平板侵入水層內(nèi)部，與附著水層同步運(yùn)動(dòng)。對(duì)于圓柱形物體，由于其與水層接觸面較小，當(dāng)物體質(zhì)量減小至=0.03時(shí)，發(fā)生反彈現(xiàn)象。

圖4 不同質(zhì)量物體穿越水層相圖

3.2 不同運(yùn)動(dòng)形式

本文選取兩種典型形狀的物體（平板和45°楔形體）進(jìn)行計(jì)算分析。分別針對(duì)這兩種外形物體，計(jì)算勻速及勻加速運(yùn)動(dòng)穿越有限水層的過程。兩種運(yùn)動(dòng)形式的流場(chǎng)形態(tài)大體一致，即在穿越水層過程中，水層中有一部分流體附著在物體表面跟隨物體運(yùn)動(dòng)，并隨著物體運(yùn)動(dòng)逐漸滑落，如圖5所示。

圖5 不同形狀物體穿越水層相圖

由于物體外形不同，在穿越水層過程中，流體與物面的相對(duì)運(yùn)動(dòng)趨勢(shì)存在一定差別，以致受力情況不同，如圖6所示。

圖6 不同形狀物體穿越水層的受力曲線

對(duì)于平板，勻速運(yùn)動(dòng)穿越水層初始時(shí)刻，水層在平板作用下突然加速，產(chǎn)生較大的附連水質(zhì)量力作用于平板。之后水層運(yùn)動(dòng)速度與平板相近，平板受力逐漸減小至接近零值。對(duì)勻加速運(yùn)動(dòng)過程，平板之上水層量變化較少，在勻加速運(yùn)動(dòng)作用下，平板受力接近恒定（=1.0~2.0）。

對(duì)于楔形體，勻速運(yùn)動(dòng)穿越水層過程中，由于其受力面逐漸增大，受力由零值逐漸增大，至=0.5楔面全部沾濕，阻力達(dá)到最大值。

平板在加速運(yùn)動(dòng)中，其阻力主要來源于附連水質(zhì)量力。在運(yùn)動(dòng)中，水層在楔面上逐漸滑落，產(chǎn)生附連水質(zhì)量力的同時(shí)，產(chǎn)生粘性阻力。

3.3 水層厚度的影響

物體在水中運(yùn)動(dòng)過程中，阻礙物體運(yùn)動(dòng)的外力（等效阻力）包括迎流面高壓產(chǎn)生阻礙運(yùn)動(dòng)的“壓力”以及背流面低壓產(chǎn)生的“吸力”。從受力分解的角度可知，阻力包括勻速運(yùn)動(dòng)而形成的位置力以及由于加速運(yùn)動(dòng)所產(chǎn)生的附連水質(zhì)量力。

水層厚度及在物體周圍的分布情況，同樣會(huì)對(duì)其運(yùn)動(dòng)受力產(chǎn)生重要影響。為了分析其具體影響情況，本文對(duì)比計(jì)算了物體在周圍全部是流體的水下運(yùn)動(dòng)過程、出水運(yùn)動(dòng)過程以及穿越不同厚度水層的運(yùn)動(dòng)過程，如圖7所示。加速度為10 m/s2的勻加速直線運(yùn)動(dòng)的不同運(yùn)動(dòng)過程的物體受力結(jié)果如圖8所示。

圖7 不同水層分布情況示意

圖8 不同形狀物體穿越水層受力曲線

在水下運(yùn)動(dòng)過程中，這兩部分力同時(shí)來源于物體的迎流面和背流面，阻力量值較大。對(duì)于出水過程，在初始階段阻力量值與水下運(yùn)動(dòng)過程相近；出水后，尾部阻力逐漸減小以至消失，阻力量值整體減小，與穿越水層工況量值相當(dāng)。

針對(duì)平板穿越水層過程，不同厚度水層使附連水質(zhì)量的影響存在差異。在水層厚度較大的情況下，平板上部水量較大，附連水質(zhì)量力更大。針對(duì)楔形體，由于水層由斜面向下滑落，水層厚度增大到一定程度時(shí)，物體上部流體量相差不大。因此水層厚度由=1增加至=1.5時(shí)，流體阻力變化不明顯。

3.4 旋轉(zhuǎn)運(yùn)動(dòng)

在穿越水層過程中，物體在做平動(dòng)的同時(shí)附加有旋轉(zhuǎn)運(yùn)動(dòng)，水層變化呈現(xiàn)明顯非對(duì)稱性，如圖9~11所示。

圖9 圓柱穿越水層相圖變化（5rad/s）

圖10 平板穿越水層相圖變化（5rad/s）

圖11 楔形體穿越水層相圖變化（5rad/s）

本文通過3種典型外形物體開展水層對(duì)物體做平面三自由度剛體運(yùn)動(dòng)的影響分析可知。圓柱在旋轉(zhuǎn)過程中僅在壁面存在切向速度，水層對(duì)旋轉(zhuǎn)運(yùn)動(dòng)的影響僅表現(xiàn)為切向流體粘性阻力。不同旋轉(zhuǎn)角速度情況下流體對(duì)旋轉(zhuǎn)角速度的衰減量相差較小，如圖12~14所示。

圖12 圓柱體穿越水層旋轉(zhuǎn)角速度變化曲線

圖13 平板穿越水層旋轉(zhuǎn)角速度變化曲線

圖14 楔形體穿越水層旋轉(zhuǎn)角速度變化曲線

平板和楔形體在旋轉(zhuǎn)過程中，物面與周圍流體同時(shí)存在切向和法向相對(duì)運(yùn)動(dòng)。不同旋轉(zhuǎn)角速度下，水層的變化存在較大差異，致使物體旋轉(zhuǎn)角速度的衰減規(guī)律存在差異。

由于平板受力面積較大，旋轉(zhuǎn)運(yùn)動(dòng)在水層影響下衰減較快。平板在運(yùn)動(dòng)初始階段平拍水面，旋轉(zhuǎn)角速度增大，角速度的衰減越快。

4 結(jié) 論

本文介紹了VOF方法及其在穿越有限水層問題模擬中的應(yīng)用，實(shí)現(xiàn)了剛體動(dòng)力學(xué)方程和氣、液兩相流動(dòng)N-S方程的耦合求解，基于試驗(yàn)數(shù)據(jù)驗(yàn)證了算法的準(zhǔn)確性，分析了不同質(zhì)量、外形及運(yùn)動(dòng)形式的物體穿越有限水層的受力及運(yùn)動(dòng)變化規(guī)律，得到如下結(jié)論：

a）穿越水層物體質(zhì)量較小時(shí)，在接觸面產(chǎn)生的反作用力會(huì)將物體反彈，當(dāng)接觸面較小時(shí)，存在反彈現(xiàn)象時(shí)的物體質(zhì)量越??；

b）加速運(yùn)動(dòng)過程會(huì)產(chǎn)生較大的附連水質(zhì)量力作用于物體，勻速運(yùn)動(dòng)穿越水層后期，水層運(yùn)動(dòng)速度與物體相近，附連水質(zhì)量力逐漸減?。?/p>

c）針對(duì)不同水層厚度及在物體周圍的分布形式，隨著水層厚度增加，阻力逐漸增大，受液面滑落影響，水層厚度增加到一定程度時(shí)，流體力增量趨于不明顯；

d）對(duì)于軸對(duì)稱體，旋轉(zhuǎn)運(yùn)動(dòng)僅受物面與流體相對(duì)切向運(yùn)動(dòng)產(chǎn)生的摩擦力影響，不同旋轉(zhuǎn)角速度情況下流體對(duì)旋轉(zhuǎn)角速度的衰減量相差較小；對(duì)于非軸對(duì)稱體，物面與流體相對(duì)法向運(yùn)動(dòng)產(chǎn)生的正壓力對(duì)旋轉(zhuǎn)運(yùn)動(dòng)有較大的衰減作用。

[1] Hirt C W, Nichols B D. Volume of fluid (VOF) method for dynamics of free boundaries[J]. Journal of Computa-tional Physics, 1981(39): 201-225.

[2] Arai M, Cheng L Y, Inoue Y. A computing method for the analysis of water impact of arbitrary shaped bodies[J]. Journal of the Society of Naval Architects of Japan, 1995(176): 233-240.

[3] Schumann C. Volume of fluid computations of water entry of bow sections[C]. Poitiers: Euromech 374, 1998.

[4] Xing K Y, Jensen G, Peric M. Numerical simulation of water entry and water exit of a horizontal circular cylinder[C]. Perth: the 6th International Conference on Hydrodynamics, 2004.

[5] 王文華, 王言英. 圓柱在波浪中入水的數(shù)值模擬[J]. 上海交通大學(xué)學(xué)報(bào), 2010, 44(10): 1393-1399.

Wang W H, Wang Y Y. Numerical study on cylinder entering water in wave[J]. Journal of Shanghai Jiaotong University, 2010, 44(10): 1393-1399.

[6] 陳宇翔, 郜冶, 劉乾坤. 應(yīng)用VOF方法的水平圓柱入水?dāng)?shù)值模擬[J]. 哈爾濱工程大學(xué)學(xué)報(bào), 2011, 11(32): 1439-1442.

Chen Yuxiang, Gao Ye, Liu Qiankun. 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): 1439-1442.

[7] 王福軍. 計(jì)算流體動(dòng)力學(xué)分析[M]. 北京: 清華大學(xué), 2004.

Wang Fujun. Computational fluid dynamics analysis[M]. Beijing: Tsinghua University, 2004.

[8] Shao S D. Incompressible SPH simulation of water entry of a free-falling object[J]. Int. J. Numer. Meth. Fluids, 2009(59): 91-115.

[9] Greenhow M, Lin W M. Nonlinear free-surface effects: Experiments and theory. MIT report 83-19[C]. Cambridge: Massachusetts Institute of Technology, 1983.

Numerical Calculation Research on the Effect of Finite Water Layer on Moving Object

Zhang Nai-min1,2, Pei Jin-liang2, Zhao Yang1, You Tian-qing2, Liu Yuan-qing2

(1. Harbin Institute of Technology, Harbin, 150001;2. Beijing Institute of Astronautical Systems Engineering, Beijing, 100076)

While underwater vehicle crossing free surface, the deformation of attached water layer is one of complex phenomena during water exiting. After water exit, the separation of tail cap occur in the water layer, which has a great impact on the hydrodynamic force and the separation movement. In this paper, the MAC method is used simulate free water surface, the couple solving of multi-phase flow N-S equation and cylinder object moving equation has been conducted. Base on experimental data, the accuracy of present method has been proven. By the method, the different shape object traversing finite water layer are calculated and analyzed. The hydrodynamic force and movement during crossing water layer has been researched on the impact of different factors, which include object shape, weight of object, motion form and distribution of water layer.

Water layer; Volume of fluid; N-S equation; Underwater vehicle

1004-7182(2018)01-0012-06

10.7654/j.issn.1004-7182.20180103

U666.1

A

2017-09-11；

2017-10-11

張耐民（1980-），男，博士研究生，高級(jí)工程師，主要研究方向?yàn)樗潞叫畜w總體設(shè)計(jì)