凹形運(yùn)動(dòng)體三自由度入水運(yùn)動(dòng)特性數(shù)值研究

路麗睿，魏英杰，王 聰，王海洋，宋武超

(1. 哈爾濱工業(yè)大學(xué) 航天學(xué)院，黑龍江 哈爾濱 150001；2. 中國科學(xué)院蘇州生物醫(yī)學(xué)工程技術(shù)研究所，江蘇 蘇州 215000)

0 引 言

入水是指運(yùn)動(dòng)體由空氣介質(zhì)穿越自由液面進(jìn)入到水介質(zhì)的過程。該過程涉及到空泡流動(dòng)、湍動(dòng)渦等復(fù)雜現(xiàn)象，具有強(qiáng)非線性和非定常特性。同時(shí)多相流場與運(yùn)動(dòng)體運(yùn)動(dòng)存在耦合影響，運(yùn)動(dòng)體運(yùn)動(dòng)特性較之單一水下運(yùn)動(dòng)較為復(fù)雜。

19世紀(jì)初，Gekle等[1-2]對(duì)回轉(zhuǎn)體垂直入水開展了數(shù)值模擬研究，獲得了空泡發(fā)展演化規(guī)律和深閉合位置與弗汝德數(shù)間的關(guān)系，建立了空泡閉合射流模型。Quan等[3]基于Navier-Stokes方程對(duì)入水空泡生成、發(fā)展、頸縮和閉合過程進(jìn)行模擬。Reinhard等[4]開展了二維對(duì)稱模型入水過程空泡形成和擴(kuò)張的數(shù)值研究，建立了基于Wagner理論和水動(dòng)力占優(yōu)假設(shè)的入水空泡預(yù)測數(shù)值模型，分別研究了拋物線型和對(duì)稱楔形對(duì)入水空泡形成和擴(kuò)張的影響規(guī)律。Iranmanesh等[5]基于VOF自由液面追蹤技術(shù)開展了低弗汝德數(shù)下水平圓柱入水的數(shù)值仿真研究，研究了圓柱直徑、入水速度對(duì)入水空泡流的影響。葉取源等[6]分別采用非線性自由面理論和歐拉-拉格朗日混合邊界元方法模擬了錐體和圓盤入水空泡發(fā)生、發(fā)展、面閉合和深閉合等流動(dòng)現(xiàn)象，定性地分析了運(yùn)動(dòng)體質(zhì)量和速度對(duì)入水空泡演化過程的影響。Gong等[7]在SPH方法基礎(chǔ)上，引入非反射邊界條件（non-reflection boundary），開展了二維平面楔形體入水過程的數(shù)值模擬研究。張珂等[8 – 9]采用單相格子Boltzmann流動(dòng)模型開展了圓柱體垂直入水二維和三維數(shù)值模擬，結(jié)果表明單相格子Boltzmann流動(dòng)模型可以較為精確地捕捉到入水過程中的部分流動(dòng)特征。王聰?shù)萚10]對(duì)球體高速入水開展數(shù)值仿真，研究了空氣壓力對(duì)入水空演化的影響。馬慶鵬等[11]分別對(duì)球體和錐頭圓柱體入水空泡及流場環(huán)境影響開展了模擬研究。王平等[12]開展了楔形體波浪中入水?dāng)?shù)值研究，研究了波浪對(duì)楔形體受力及入水姿態(tài)的影響。

凹形運(yùn)動(dòng)體作為一種復(fù)雜結(jié)構(gòu)運(yùn)動(dòng)體，具有廣泛的工程應(yīng)用背景。目前僅有路中磊等[13 – 14]以火箭助推器入水為背景開展過“凹體”入水相關(guān)研究，其研究主要關(guān)注入水過程流場特性，對(duì)運(yùn)動(dòng)特性關(guān)注較少。除火箭助推器外，導(dǎo)彈尾罩亦屬于凹形運(yùn)動(dòng)體。在潛射導(dǎo)彈冷發(fā)射過程中，導(dǎo)彈點(diǎn)火時(shí)將尾罩從一側(cè)彈開，其運(yùn)動(dòng)具有典型的多自由度特性，研究尾罩入水運(yùn)動(dòng)特性對(duì)避免尾罩砸艇具有重要意義。本文以此為背景，開展基于重疊網(wǎng)格技術(shù)的凹形運(yùn)動(dòng)體三自由度入水過程運(yùn)動(dòng)特性的數(shù)值研究，分析初始垂向速度、水平遷移速度對(duì)凹形運(yùn)動(dòng)體三自由度入水過程運(yùn)動(dòng)特性及水動(dòng)力特性的影響。

1 數(shù)值計(jì)算理論及方法

1.1 基本控制方程

入水運(yùn)動(dòng)是一個(gè)非定常問題，流動(dòng)參數(shù)均隨時(shí)空變化，主要表現(xiàn)在質(zhì)量擴(kuò)散、動(dòng)量黏滯耗散等方面。在流體流動(dòng)過程中，混合相滿足質(zhì)量守恒和動(dòng)量守恒方程。

連續(xù)性方程為：

動(dòng)量守恒方程為：

1.2 多相流模型

低速入水空泡氣液兩相介質(zhì)互不摻混，空泡壁面處兩相流速基本一致，因此可采用均質(zhì)平衡多相流理論[15]。本文選用的VOF多相流模型是一種適用于固定歐拉網(wǎng)格的界面捕捉模型，該模型可以較好描述流動(dòng)過程中多相流體的界面位置。VOF模型將兩相流體視為一種密度可變的混合流體，通過逐一求解計(jì)算流域內(nèi)每一控制單元的控制方程得到混合流體參數(shù)，并根據(jù)各相流體所占的體積分?jǐn)?shù)去實(shí)現(xiàn)對(duì)流體間界面的追蹤?；旌舷嗝芏群蛣?dòng)力粘度分別如下：

1.3 湍流模型

湍流是一種具有強(qiáng)非線性的流動(dòng)狀態(tài)，同時(shí)也是最為普遍的流動(dòng)狀態(tài)。本文采用Realizable，該模型將渦粘系數(shù)（）和湍動(dòng)能（）與耗散率（）聯(lián)系起來，分別建立湍動(dòng)能（）與耗散率（）的輸運(yùn)方程，如式（5）和式（6）所示，最終使整個(gè)控制方程封閉。

1.4 重疊網(wǎng)格技術(shù)

重疊網(wǎng)格的基本思想是對(duì)復(fù)雜流域的分解與網(wǎng)格的重新組合[16]。重疊網(wǎng)格技術(shù)將流域分解為若干個(gè)子區(qū)域，各子區(qū)域計(jì)算網(wǎng)格在計(jì)算中相互獨(dú)立，網(wǎng)格之間存在重疊、嵌套和覆蓋關(guān)系。在重疊區(qū)域中，各子區(qū)域計(jì)算網(wǎng)格通過流場信息的插值與映射進(jìn)行數(shù)據(jù)交換，從而建立起了各子區(qū)域網(wǎng)格之間的耦合關(guān)系，為子區(qū)域流場計(jì)算提供邊界條件。重疊網(wǎng)格的主要過程包含尋點(diǎn)、挖洞、建立網(wǎng)格重疊關(guān)系3個(gè)部分。在尋點(diǎn)過程中，網(wǎng)格A落入網(wǎng)格B中的節(jié)點(diǎn)稱為洞內(nèi)點(diǎn)，洞內(nèi)點(diǎn)在數(shù)值求解過程中不參加計(jì)算，與洞內(nèi)點(diǎn)緊鄰的點(diǎn)稱為邊界插值點(diǎn)，邊界插值點(diǎn)通過插值與映射接收來自其他子區(qū)域網(wǎng)格的數(shù)據(jù)信息，并將所接受的數(shù)據(jù)作為該區(qū)域求解的邊界條件。如圖1所示，網(wǎng)格A通過其邊界插值點(diǎn)接收B網(wǎng)格內(nèi)求解得到的流場信息，同樣，網(wǎng)格B也從其外邊界插值點(diǎn)獲得A網(wǎng)格的流場數(shù)據(jù)，并將該數(shù)據(jù)作為求解邊界條件。

1.5 數(shù)值計(jì)算方法

本文基于Fluent 17.2平臺(tái)，采用有限體積法離散流動(dòng)控制方程；壓力場與速度場的耦合選用Coupled算法；壓力場的空間離散采用PRESTO!格式；各相體積率離散采用CICSAM格式；各求解變量的離散采用一階迎風(fēng)格式；對(duì)流項(xiàng)采用QUICK離散格式；相界面的幾何重構(gòu)采用Geo-Reconstruct格式。通過C語言程序，嵌入Fluent用戶自定義函數(shù)定義入水運(yùn)動(dòng)體慣性矩及邊界壓力，實(shí)現(xiàn)模型三自由度運(yùn)動(dòng)。

早年嫁過來時(shí)，采蓮掮著鋤頭跟福生做農(nóng)活，不惟下湖，還得去崖上。崖頭，就是高處的地坡。正午坐在地頭上喝水，偶爾一搭眼，就看到天宇下的崖坡上，一人一牛在耕地，遠(yuǎn)看兩個(gè)黑點(diǎn)，背后一抹駱駝云，一趟一趟，不厭其煩地劃著圈……此后，她的夢里始終留存著這樣一組畫面。農(nóng)夫，耕牛，在天宇下的土坡上，永遠(yuǎn)沒有盡頭地轉(zhuǎn)悠，一簇牛角上的紅布條像火苗似的燃燒著。而那一聲吆牛號(hào)子，那份游蕩與戚然交織的天籟啊，簡直就是入心入肺了。[1]12

2 計(jì)算模型及網(wǎng)格劃分

2.1 幾何模型

本文計(jì)算所用凹形運(yùn)動(dòng)體幾何模型如圖2所示，模型幾何參數(shù)及屬性見表1。

表 1 模型幾何參數(shù)及屬性Tab. 1 The parameter of geometric model

2.2 計(jì)算域及邊界條件

為簡化計(jì)算，可假設(shè)凹形運(yùn)動(dòng)體在z方向不存在速度分量，且不存在繞y軸的轉(zhuǎn)動(dòng)速度分量。故本算例的數(shù)值計(jì)算選擇在三維、三自由度情況下進(jìn)行。其計(jì)算域及邊界條件如圖4所示，xy平面為對(duì)稱平面，計(jì)算域?qū)挒?0D，高為7D，厚為10D。以觸水時(shí)刻作為時(shí)間零點(diǎn)，觸水點(diǎn)為坐標(biāo)原點(diǎn)，重力方向?yàn)閥負(fù)方向。

2.3 網(wǎng)格劃分

圖4給出了計(jì)算域網(wǎng)格和子域網(wǎng)格的劃分方式，2套網(wǎng)格均采用正六面體結(jié)構(gòu)化網(wǎng)格，以較好控制網(wǎng)格質(zhì)量和網(wǎng)格數(shù)量。為保證空泡區(qū)域計(jì)算精度，并保持較好的過度效果，對(duì)凹形運(yùn)動(dòng)體附近區(qū)域采用局部加密處理，運(yùn)動(dòng)體向遠(yuǎn)場過渡區(qū)域采用漸疏處理。

2.4 數(shù)值方法驗(yàn)證

為了驗(yàn)證數(shù)值方法的有效性，基于上述計(jì)算方法及邊界條件，對(duì)Aristoff[17]所開展的球體入水實(shí)驗(yàn)進(jìn)行計(jì)算。實(shí)驗(yàn)研究所選球體直徑D=25.4 mm，球體密度為 2.3 g/cm3，入水速度 v0=2.17 m/s。本節(jié)數(shù)值計(jì)算采用1:1的計(jì)算域模型，運(yùn)動(dòng)參數(shù)初始條件設(shè)置與Aristoff實(shí)驗(yàn)狀態(tài)完全一致，計(jì)算邊界條件及重疊網(wǎng)格劃分方法與2.2節(jié)和2.3節(jié)所述相同。圖5和圖6分別給出了計(jì)算所得空泡形態(tài)和運(yùn)動(dòng)軌跡與試驗(yàn)的對(duì)比圖，從圖中可以看出，數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)具有較好的一致性。

3 結(jié)果分析

3.1 凹形運(yùn)動(dòng)體三自由度入水運(yùn)動(dòng)過程

導(dǎo)彈尾罩在拋落時(shí)，為了避免尾罩垂向下落砸艇，往往會(huì)在導(dǎo)彈尾罩一側(cè)安裝拉斷式鉸鏈，使尾罩彈開后向一側(cè)飛離，如圖7所示。因此本文對(duì)凹形運(yùn)動(dòng)體三自由運(yùn)動(dòng)初始方向進(jìn)行如圖7所示定義，水平遷移速度沿x正方向，垂向速度沿y負(fù)方向，旋轉(zhuǎn)方向正負(fù)由右手定則給出，初始角速度方向?yàn)檎?/p>

圖8給出了凹形運(yùn)動(dòng)體以135°傾角、垂向初始速度 v0=7.5 m/s、水平遷移速度 u0=2 m/s、角速度 ω0=2.0 rad/s狀態(tài)入水過程氣液兩相體積分?jǐn)?shù)變化云圖。在凹形運(yùn)動(dòng)體入水過程中，下邊緣撞擊自由液面，流體質(zhì)點(diǎn)獲得來自運(yùn)動(dòng)體的動(dòng)能而產(chǎn)生向外排開的運(yùn)動(dòng)，如圖8中t=0 ms所示，由于運(yùn)動(dòng)體凸面與凹面幾何特性的差異，導(dǎo)致凹面尖端更易誘導(dǎo)空泡形成。隨著入水時(shí)間的增加，噴濺液體在凹面內(nèi)部閉合，形成隨動(dòng)空泡，如圖8中t=5 ms所示。在閉合瞬間，噴濺液體在凹面上半側(cè)產(chǎn)生抨擊作用而形成力矩，運(yùn)動(dòng)體旋轉(zhuǎn)方向發(fā)生改變，定義該運(yùn)動(dòng)為回旋運(yùn)動(dòng)。運(yùn)動(dòng)體的回旋導(dǎo)致由上緣誘導(dǎo)形成的二次空泡被橫向拉伸，空泡橫向發(fā)展，如圖8中t=10 ms至t=15 ms所示。隨著空泡的橫向發(fā)展，迎流面一側(cè)水流體在重力作用下下落，與背流面一側(cè)空泡界面相遇，二次空泡完成閉合，如圖8中t=20 ms至t=25 ms所示。二次空泡受慣性力作用繼續(xù)做橫向運(yùn)動(dòng)，逐漸與自由液面邊界分離而形成了獨(dú)立空泡，如圖8中t=30 ms至35 ms所示。隨動(dòng)空泡亦由于凹形運(yùn)動(dòng)體旋轉(zhuǎn)運(yùn)動(dòng)逐漸脫落形成獨(dú)立空泡，如圖8中t=15 ms所示。獨(dú)立空泡形態(tài)脫落后受運(yùn)動(dòng)體凹面旋轉(zhuǎn)運(yùn)動(dòng)形成的剪切力作用而逐漸卷曲拉長，如圖 8 中 t=15 ms至 t=35 ms所示。

對(duì)凹形運(yùn)動(dòng)體所受到的水動(dòng)力系數(shù)做如下定義：

圖9給出了凹形運(yùn)動(dòng)體入水過程運(yùn)動(dòng)特性曲線，圖10給出了該過程水動(dòng)力特性曲線，圖11給出了凹形結(jié)構(gòu)體位置姿態(tài)隨時(shí)間變化簡圖。從圖9中可以看出在入水之前，凹形運(yùn)動(dòng)體y方向速度以較小變化率減小，角速度亦呈現(xiàn)較為平緩的正向增加趨勢；在入水后的極短時(shí)間內(nèi)，y方向速度迅速減小，角速度也達(dá)到峰值，之后迅速減小為0并反向增大，稱該階段為沖擊階段；之后，y方向速度變化率逐漸減緩，速度持續(xù)減小并逐漸趨于0，角速度變化率在t=5～13 ms之間出現(xiàn)短暫的微幅增大，t=10 ms之后又逐漸減小，角速度達(dá)到負(fù)方向峰值后亦逐漸減小。對(duì)比圖10可知，在入水之前，凹形運(yùn)動(dòng)體靠近自由液面，其下緣與自由液面間形成氣墊，受空氣墊阻力作用，y方向速度緩慢減小，由于空氣墊阻力與凹形運(yùn)動(dòng)體所示重力形成了與角速度方向相同的力矩，導(dǎo)致角速度亦呈現(xiàn)緩慢增加趨勢。在觸水后，進(jìn)入沖擊階段，沖擊階段主要包含入水沖擊和噴濺沖擊2個(gè)部分，在入水沖擊過程中，凹形運(yùn)動(dòng)體受到入水沖擊力作用，y方向速度迅速減小。y方向入水沖擊力在質(zhì)心處產(chǎn)生與角速度方向相同的力矩，x方向沖擊力在質(zhì)心處產(chǎn)生與角速度方向相反的力矩。由于y方向沖擊力遠(yuǎn)大于x方向沖擊力，且其值遠(yuǎn)大于空氣墊阻力，因此合力矩與角速度方向相同，運(yùn)動(dòng)體角速度迅速增加。隨著入水時(shí)間的增加，運(yùn)動(dòng)體背流面所形成的液體噴濺達(dá)到凹形運(yùn)動(dòng)體內(nèi)壁面上側(cè)（圖8中t=0～5 ms之間），此時(shí)入水沖擊階段已經(jīng)完成，運(yùn)動(dòng)體主要受噴濺沖擊力作用，進(jìn)入噴濺沖擊階段。噴濺沖擊力在質(zhì)心處所產(chǎn)生的力矩方向與旋轉(zhuǎn)方向相反，致使角速度迅速減小并反向增加，形成回旋運(yùn)動(dòng)。噴濺沖擊階段結(jié)束后，凹形運(yùn)動(dòng)體下緣形成了空泡隨動(dòng)，且由于空泡內(nèi)壓力小于凸面處所受環(huán)境壓力，導(dǎo)致迎流面與背流面壓差形成力矩作用，該力矩方向與運(yùn)動(dòng)體回旋方向相同，導(dǎo)致角速度變化率短暫增大。在t=13 ms附近，隨動(dòng)空泡脫落，運(yùn)動(dòng)體處于全沾濕狀態(tài)，其旋轉(zhuǎn)受到流體阻力和轉(zhuǎn)動(dòng)慣性力共同作用，最終在t=18 ms時(shí)角速度達(dá)到峰值并開始減小。從圖11還可以看出，由于運(yùn)動(dòng)體撞水面積較小，使得入水沖擊較噴濺沖擊弱；在t=14 ms附近，運(yùn)動(dòng)體在x方向的迎流面積達(dá)到最大，x方向阻力達(dá)到峰值。結(jié)合圖8中t=10～16 ms階段，由于運(yùn)動(dòng)體上緣與空泡尾腔分離，空泡界面與運(yùn)動(dòng)體邊緣之間存在張力，張力的拖曳作用在y方向的分量使得y方向阻力開始微幅增大，并在t=16 ms時(shí)達(dá)到峰值。

3.2 垂向速度對(duì)凹形運(yùn)動(dòng)體運(yùn)動(dòng)特性的影響

考慮不同垂向速度，圖12給出了凹形運(yùn)動(dòng)體入水過程質(zhì)心平移特性變化規(guī)律，分別取初始垂向速度v0=5.0 m/s，7.5 m/s，10.0 m/s ，初始水平遷移速度 u0=1 m/s，初始角速度 ω0=2.0 rad/s。從圖 12（a）中可以看出隨著垂向速度的增加，y方向速度趨于0時(shí)，其最大位移逐漸增大，但增大幅度隨速度的增加而減小。由圖12（b）可知，在沖擊階段，隨著垂向速度的增加，速度衰減程度逐漸增大。從圖12（c）中還可發(fā)現(xiàn)，隨著垂向速度的增大，“Z”型路徑逐漸平坦，表明垂向速度增加對(duì)運(yùn)動(dòng)體回?cái)[存在抑制作用。

為分析上述現(xiàn)象原因，圖13給出了運(yùn)動(dòng)體x方向和y方向的水動(dòng)力系數(shù)變化曲線。從圖13（b）可以看出，隨著入水速度的增加，y方向入水沖擊力以及噴濺沖擊力峰值增量大幅加大，導(dǎo)致沖擊引起的速度衰減幅度增加。在入水沖擊之后，y方向受到的阻力亦隨著速度增加而加大，沖擊與阻力對(duì)y方向最大位移具有負(fù)激勵(lì)作用，隨著速度的增加，負(fù)激勵(lì)增大，導(dǎo)致位移增加幅度減小。

圖14給出了垂向速度對(duì)旋轉(zhuǎn)運(yùn)動(dòng)特性的影響規(guī)律，從圖中可以看出，隨著垂向速度的增大，入水初期角位移峰值逐漸前移，角速度峰值逐漸增大。回旋運(yùn)動(dòng)階段角速度亦隨著垂向速度的增大而增大。由于回旋運(yùn)動(dòng)主要由噴濺沖擊引起，噴濺沖擊力對(duì)回旋運(yùn)動(dòng)角速度具有正激勵(lì)作用，從圖13（a）可以看出，隨著垂向速度的增加，噴濺沖擊力大幅增加，導(dǎo)致回旋運(yùn)動(dòng)角速度增大。

3.3 水平遷移速度對(duì)凹形運(yùn)動(dòng)體運(yùn)動(dòng)特性的影響

圖15給出了水平遷移速度對(duì)凹形運(yùn)動(dòng)體平移運(yùn)動(dòng)的影響規(guī)律，初始水平遷移速度分別取u0=2.0 m/s，3.0 m/s，4.0 m/s，初始垂向速度取 v0=6.0 m/s，初始角速度ω0=2.0 rad/s。從圖15（a）中可以看出，不同水平遷移速度下，凹形運(yùn)動(dòng)體運(yùn)動(dòng)具有相同的運(yùn)動(dòng)趨勢，且隨著水平遷移速度增加，y方向速度趨于0時(shí)，其位移最大值增加。當(dāng)u0=3 m/s時(shí)，回?cái)[運(yùn)動(dòng)引起了質(zhì)心的上漂。從圖15（b）中可發(fā)現(xiàn)，隨著水平遷移速度的增加，入水沖擊階段，y方向速度衰減幅值逐漸減小。觀察圖15（c）可發(fā)現(xiàn)，隨著速度的增加，質(zhì)心軌跡回?cái)[幅值增大，質(zhì)心運(yùn)動(dòng)軌跡的“Z”型特征隨著水平遷移速度的增加而逐漸顯著。

圖16給出了水平遷移速度對(duì)x和y方向水動(dòng)力系數(shù)的影響規(guī)律，從圖中可發(fā)現(xiàn)，沖擊階段y方向受到的入水及噴濺沖擊作用與入水瞬間速度衰減幅值呈現(xiàn)相同趨勢。隨著水平速度的增加，x方向的入水沖擊力增大，但噴濺沖擊力卻減小。出現(xiàn)該現(xiàn)象是由于水平速度增加，噴濺到達(dá)運(yùn)動(dòng)體背流面時(shí)與運(yùn)動(dòng)體的相對(duì)速度減小，導(dǎo)致噴濺沖擊作用減小。在沖擊階段之后，y方向所受到的阻力以及x方向阻力均隨水平遷移速度的增加而增大。

4 結(jié) 語

本文基于重疊網(wǎng)格技術(shù)，開展了凹形運(yùn)動(dòng)體入水過程運(yùn)動(dòng)特性數(shù)值研究，分析了初始垂向速度、水平遷移速度對(duì)凹形運(yùn)動(dòng)體入水過程運(yùn)動(dòng)特性及水動(dòng)力特性的影響，得到如下結(jié)論：

1）在凹形運(yùn)動(dòng)體三自由度入水過程中，噴濺沖擊作用對(duì)凹形運(yùn)動(dòng)體角運(yùn)動(dòng)具有重要影響，運(yùn)動(dòng)體在水中旋轉(zhuǎn)方向與噴濺沖擊形成的力矩方向相同。

2）凹形運(yùn)動(dòng)體入水運(yùn)動(dòng)存在2個(gè)運(yùn)動(dòng)占優(yōu)階段，在入水前期為平移運(yùn)動(dòng)占優(yōu)，入水后期為旋轉(zhuǎn)運(yùn)動(dòng)占優(yōu)。

3）凹形運(yùn)動(dòng)體在下落過程中，質(zhì)心運(yùn)動(dòng)軌跡會(huì)產(chǎn)生回?cái)[而呈現(xiàn)“Z”型特征，垂向速度的增加對(duì)質(zhì)心回?cái)[運(yùn)動(dòng)具有抑制作用，水平遷移速度的增加對(duì)其具有促進(jìn)作用。