張　旭　劉沛清　屈秋林 /

(北京航空航天大學(xué)，北京100191)

[基金項(xiàng)目]本文受國(guó)家自然科學(xué)基金(11502012、11772033)支持。

0　引言

大型飛機(jī)水上迫降性能的研究[1-4]離不開(kāi)對(duì)物體入水過(guò)程的研究。在飛機(jī)入水的過(guò)程中，會(huì)有很大的沖擊載荷和沖擊壓強(qiáng)出現(xiàn)在結(jié)構(gòu)上，同時(shí)后機(jī)身出現(xiàn)很大的吸力導(dǎo)致飛機(jī)出現(xiàn)大幅度的抬頭。后機(jī)身的吸力主要是由于飛機(jī)水平方向的運(yùn)動(dòng)引起的，而沖擊載荷和沖擊壓強(qiáng)主要是由于飛機(jī)的下降速度導(dǎo)致的。

影響物體入水沖擊過(guò)程的因素有很多，譬如物體的幾何形狀、姿態(tài)和沖擊速度等。大部分研究主要集中于軸對(duì)稱或二維物體進(jìn)入平靜水面的過(guò)程。首先研究物體入水沖擊問(wèn)題的是Von Krmn[5]，他在研究水上飛機(jī)的入水性能時(shí)使用動(dòng)量守恒和附加質(zhì)量法求解了二維楔形體入水的垂向載荷。隨后，Wagner[6]考慮了水面抬升并通過(guò)擴(kuò)張平板假設(shè)得出了更好的結(jié)果。Cointe、Armand[7]、Oliver[8]和Korobkin[9]等使用了匹配漸進(jìn)展開(kāi)法求解二維圓柱入水過(guò)程的載荷。Zhao等[10]使用邊界元方法求解了二維楔形體和艏外飄截面入水過(guò)程中的水面形狀、垂向載荷和壓強(qiáng)分布。Mei等[11]采用了一種解析方法求解二維楔形體、二維圓柱和艏外飄截面入水過(guò)程中的載荷。Greenhow[12]、Sun和Faltinsen[13]采用邊界元方法求解二維圓柱入水過(guò)程中的載荷，而Vandamme等[14]則使用了SPH方法來(lái)求解二維圓柱入水問(wèn)題。王明振等[15]采用實(shí)驗(yàn)方法研究三種水陸兩棲飛機(jī)典型橫截面在不同投放高度和質(zhì)量下的入水沖擊中壓強(qiáng)和載荷的變化。上述研究主要集中于楔形體、圓柱、艏外飄截面和飛機(jī)截面等二維物體的入水沖擊過(guò)程。

三維物體入水問(wèn)題的研究主要集中于航天器返回艙入水，而對(duì)于帶中央翼盒的機(jī)身形狀入水的研究很少。Mcgehee等[16]對(duì)Mercury返回艙的1/12比例模型以姿態(tài)角范圍在-30°～30°之間的入水沖擊載荷進(jìn)行了實(shí)驗(yàn)和理論分析。Stubbs[17]采用Apollo返回艙的1/4比例模型實(shí)驗(yàn)研究了不同速度和姿態(tài)角的著水沖擊載荷。Wang等[18]使用LS-DYNA研究了Apollo返回艙以不同姿態(tài)角入水過(guò)程中的加速度變化，結(jié)果發(fā)現(xiàn)Von Kármán方法低估了初始的沖擊加速度而Wagner方法則出現(xiàn)了高估。以上的研究表明姿態(tài)角、速度和物體形狀都會(huì)對(duì)入水過(guò)程產(chǎn)生影響，同時(shí)三維物體的形狀復(fù)雜，將二維的結(jié)果直接應(yīng)用到三維誤差較大。

本文利用數(shù)值模擬手段，對(duì)帶中央翼盒的細(xì)長(zhǎng)體機(jī)身在不同姿態(tài)角、入水速度和機(jī)身尾翹角情況下的入水沖擊過(guò)程進(jìn)行數(shù)值模擬，通過(guò)研究各參數(shù)變化的影響規(guī)律，達(dá)到揭示機(jī)理之目的。

技 術(shù) 研 究

1　物理模型

圖1為機(jī)身垂直入水的示意圖。笛卡爾坐標(biāo)系原點(diǎn)定義在初始的平靜水面上，x方向平行于平靜水面指向機(jī)身后方，y方向鉛垂向上，z方向由右手螺旋定則確定；姿態(tài)角a是機(jī)身軸線與平靜水面的夾角；x1和x2是機(jī)身的兩個(gè)典型橫截面；機(jī)身以恒定速度V向下沖擊平靜水面，入水后機(jī)身受到的垂向力Fy，機(jī)身浸沒(méi)在平靜水面下的浮力Fb，兩者之差定義為機(jī)身受到的沖擊力Fi。

圖2展示了計(jì)算中的機(jī)身模型(參照某型飛機(jī))。計(jì)算中采用1/10的縮比模型，計(jì)算模型的機(jī)身總長(zhǎng)L=3.85m，后體長(zhǎng)La=36.3%L，當(dāng)量直徑為D=0.31m。尾翹角分別為β=3°,5°7°，如圖2(b)所示。圖2(c)展示了不同姿態(tài)角(α=8°,12°)下的機(jī)身橫截面形狀，其中x1=2.1m處為中央翼盒的典型截面，x2=3.3m處為機(jī)身尾部的典型截面，可以發(fā)現(xiàn)它們的形狀基本相似。

(a) 機(jī)身側(cè)視圖

(b) 后機(jī)身

(c) 機(jī)身橫截面形狀圖2　機(jī)身模型

2　計(jì)算方法和驗(yàn)證

2.1 流場(chǎng)求解

計(jì)算中采用ANSYS FLUENT 14.0雙精度求解器來(lái)求解非定?？蓧嚎s的RANS方程，考慮了空氣和水的重力。湍流模擬選取可實(shí)現(xiàn)的k-ε模型和加強(qiáng)壁面函數(shù)處理。壓強(qiáng)速度耦合采用SIMPLE方法。對(duì)流項(xiàng)采用三階MUSCL格式，擴(kuò)散項(xiàng)采用二階中心格式，非定常項(xiàng)采用一階隱式格式。

2.2 自由面模型

自由面的捕捉采用VOF模型，模型通過(guò)對(duì)每一個(gè)相定義體積分?jǐn)?shù)來(lái)模擬該相，每一個(gè)網(wǎng)格單元中所有相體積分?jǐn)?shù)之和等于1。對(duì)于一個(gè)網(wǎng)格單元，q相的體積分?jǐn)?shù)定義為:γq=0意味著該單元中沒(méi)有q相；γq=1意味著該單元中充滿q相；而0<γq<1表示該單元是相間的交界面。

第q相體積分?jǐn)?shù)的連續(xù)方程為：

(1)

2.3 整體動(dòng)網(wǎng)格技術(shù)

整體動(dòng)網(wǎng)格技術(shù)(Qu等[19])用來(lái)模擬機(jī)身和水體之間的相對(duì)運(yùn)動(dòng)。計(jì)算過(guò)程中，整個(gè)計(jì)算區(qū)域及其內(nèi)部的網(wǎng)格隨著機(jī)身做剛體運(yùn)動(dòng)，這種方式不會(huì)出現(xiàn)網(wǎng)格的變形和重構(gòu)，保證了網(wǎng)格的質(zhì)量并節(jié)約了計(jì)算資源。

機(jī)身表面采用無(wú)滑移壁面邊界條件；計(jì)算區(qū)域的邊界為速度入口邊界條件，邊界上的體積分?jǐn)?shù)設(shè)定可以保證在網(wǎng)格的運(yùn)動(dòng)中自由面保持不動(dòng)。

圖3是計(jì)算區(qū)域和機(jī)身附近的網(wǎng)格劃分情況。計(jì)算中采用半模，計(jì)算區(qū)域大小為3L×2L×L，如圖3(a)所示，其中xy平面為對(duì)稱面。結(jié)構(gòu)網(wǎng)格可以更好的捕捉水面，因而計(jì)算中采用了結(jié)構(gòu)網(wǎng)格并對(duì)機(jī)身物面附近進(jìn)行加密。圖3(b)展示了機(jī)身附近的網(wǎng)格劃分情況，最終采用的網(wǎng)格數(shù)量大約140萬(wàn)。

(a) 計(jì)算區(qū)域

(b) 機(jī)身附近網(wǎng)格圖3　計(jì)算區(qū)域和機(jī)身附近網(wǎng)格

2.4 網(wǎng)格數(shù)量和時(shí)間步長(zhǎng)的依賴性驗(yàn)證

選取尾翹角為β=5°的機(jī)身模型以α=12°姿態(tài)角，V=0.5 m/s下沉速度的入水過(guò)程，開(kāi)展網(wǎng)格數(shù)量和時(shí)間步長(zhǎng)的依賴性驗(yàn)證。

圖4是70萬(wàn)(Coarse)、140萬(wàn)(Normal)和280萬(wàn)(Fine)網(wǎng)格數(shù)量計(jì)算的垂向載荷系數(shù)時(shí)間歷程對(duì)比，其中S=π(D/2)2為典型機(jī)身界面的面積?？梢钥闯?0萬(wàn)網(wǎng)格與140萬(wàn)差距較大，而140萬(wàn)與280萬(wàn)差距較小，因此下文計(jì)算中選取140萬(wàn)網(wǎng)格是合理的。

圖4　不同網(wǎng)格數(shù)量計(jì)算的垂向載荷系數(shù)時(shí)間歷程

圖5是Δt=10-3s、10-4s和10-5s時(shí)間步長(zhǎng)計(jì)算的垂向載荷系數(shù)時(shí)間歷程對(duì)比，可以看出時(shí)間步長(zhǎng)為10-3s的結(jié)果與10-4s的結(jié)果差距較大，而10-4s與10-5s的結(jié)果差距較小，因此下文計(jì)算中選取的時(shí)間步長(zhǎng)為10-4s。

圖5　不同時(shí)間步長(zhǎng)計(jì)算的垂向載荷系數(shù)時(shí)間歷程

2.5 計(jì)算方法驗(yàn)證

本文選取Lin和Shihe[20]的圓柱入水實(shí)驗(yàn)來(lái)驗(yàn)證本文計(jì)算方法的精度。實(shí)驗(yàn)中的壓強(qiáng)傳感器為Kyowa PGM-2KC(直徑5.5 mm，頻率24 KHz，量程2×105Pa)，照相系統(tǒng)為CCD-16 MHz。圓柱材料為丙烯酸樹(shù)脂，圓柱長(zhǎng)度20 cm，直徑20 cm，重量12.5 kg。自由降落的初始高度為0～20 cm，入水速度范圍0.76 m/s～1.98 m/s。壓強(qiáng)監(jiān)測(cè)點(diǎn)相對(duì)于中心線的圓心角為0°、7.5°、15°和30°。驗(yàn)證選取入水速度為0.99 m/s的實(shí)驗(yàn)結(jié)果。計(jì)算中采用上述數(shù)值模擬方法對(duì)二維圓柱的入水過(guò)程進(jìn)行模擬。

圖6　實(shí)驗(yàn)和模擬的壓強(qiáng)系數(shù)的時(shí)間歷程對(duì)比

圖6是實(shí)驗(yàn)和模擬計(jì)算的壓強(qiáng)系數(shù)的時(shí)間歷程對(duì)比。其中壓強(qiáng)系數(shù)定義為Cp=(p-p0)/(0.5ρV2)，其中p為當(dāng)?shù)貕簭?qiáng)，p0為大氣壓強(qiáng)，ρw為水的密度，V為觸水時(shí)刻的速度。從圖中可以看出模擬結(jié)果與實(shí)驗(yàn)結(jié)果的峰值和趨勢(shì)都吻合得很好。實(shí)驗(yàn)結(jié)果中7.5°和15°中壓強(qiáng)幾乎沒(méi)有波動(dòng)而模擬中有波動(dòng)，這個(gè)差異是由于實(shí)驗(yàn)中三維效應(yīng)導(dǎo)致該位置沒(méi)有空氣泡出現(xiàn)，而二維計(jì)算模型中出現(xiàn)了空氣泡。本文的數(shù)值方法可以很好地模擬入水過(guò)程。

3　計(jì)算結(jié)果與分析

為了研究入水速度、姿態(tài)角和機(jī)身尾翹角對(duì)機(jī)身垂直入水沖擊過(guò)程的影響，計(jì)算了不同尾翹角(β=3°,5°,7°)的機(jī)身以常速度進(jìn)入水體的過(guò)程，其中入水速度分別為V=0.25 m/s,0.5 m/s，姿態(tài)角分別為α=8°,12°。

3.1 典型的機(jī)身入水過(guò)程

選取尾翹角為β=5°的機(jī)身模型以姿態(tài)角α=12°，V=0.5 m/s下沉速度的垂向入水過(guò)程為典型算例。圖7為該過(guò)程中的水面位置和機(jī)身腹部壓強(qiáng)分布。

圖7　典型入水過(guò)程(β=5°，α=12°，V=0.5 m/s)中水面形狀和機(jī)身腹部壓強(qiáng)分布

從圖7(a)中可以看出在剛?cè)胨蟛痪玫膖=0.02 s，壓強(qiáng)峰值很大且出現(xiàn)在噴濺根部。隨著入水深度增加，壓強(qiáng)峰值略有減小，且與底部區(qū)域的壓強(qiáng)差距有所減小，如圖7(b)所示。圖7(c)為t=0.2 s時(shí)刻，此時(shí)壓強(qiáng)峰值已經(jīng)轉(zhuǎn)移到了機(jī)身底部，而且水面已經(jīng)到了中央翼盒的位置，導(dǎo)致局部的壓強(qiáng)增大。從圖7(d)可以看出，當(dāng)機(jī)身浸沒(méi)深度繼續(xù)增加，壓強(qiáng)峰值一直在機(jī)身的底部，同時(shí)由于沖擊效應(yīng)中央翼盒與水面接觸的前部位置出現(xiàn)了局部的壓強(qiáng)極大值。

在船舶研究領(lǐng)域，人們通常使用切片法，用純2D物體的入水沖擊特性代替3D物體截面的入水沖擊特性。為了研究這種近似帶來(lái)的誤差，本文對(duì)比典型機(jī)身截面和形狀相同的純2D物體入水沖擊特性的差別。

圖8是入水過(guò)程中機(jī)身x2截面和2D物體的表面壓強(qiáng)系數(shù)對(duì)比，其中y=0為初始時(shí)刻平靜水面的位置。可以發(fā)現(xiàn)3D和2D的壓強(qiáng)系數(shù)分布具有一定的相似性，而它們的差別主要來(lái)自3D的縱向流動(dòng)和沖擊。在圖8(a)中的t=0.02 s時(shí)刻，壓強(qiáng)系數(shù)峰值并沒(méi)有出現(xiàn)在物面底部即y/D最小的地方，而是出現(xiàn)在水面附近的噴濺根部，由于存在縱向流動(dòng)3D截面上的壓強(qiáng)系數(shù)峰值小于2D；隨著y/D繼續(xù)增大壓強(qiáng)系數(shù)迅速減小；最后y/D>0.1的部分是由于上方?jīng)]有接觸到水因而物面壓強(qiáng)系數(shù)基本為零。隨著入水深度的增加，在圖8(b)所示的t=0.1 s時(shí)刻，壓強(qiáng)系數(shù)峰值已經(jīng)出現(xiàn)在了y/D最小的物面底部，由于縱向存在沖擊這一時(shí)刻3D的壓強(qiáng)系數(shù)峰值大于2D；隨著y/D的增大壓強(qiáng)系數(shù)迅速減小，在y/D=0.1附近出現(xiàn)了一些負(fù)壓，這是因?yàn)闆_擊后出現(xiàn)了水面抬升，高于平靜水面的水受到了重力作用；隨后壓強(qiáng)系數(shù)上升到零附近之后保持不變(y/D>0.3)。最后圖8(c)中的t=0.3 s時(shí)刻，壓強(qiáng)系數(shù)峰值也是在y/D最小的底部，此時(shí)縱向沖擊影響減弱使得3D的壓強(qiáng)系數(shù)峰值略小于2D；隨后迅速減小直到零附近保持不變。

(a) t=0.02 s

(b) t=0.1 s

(c) t=0.3 s圖8　入水過(guò)程(β=5°，α=12°，V=0.5 m/s)中機(jī)身x2截面和2D物體表面壓強(qiáng)系數(shù)對(duì)比

圖9是入水過(guò)程中機(jī)身x1截面和2D物體的表面壓強(qiáng)系數(shù)對(duì)比?？梢钥闯鲈摻孛嫔蠅簭?qiáng)系數(shù)分布的變化規(guī)律與圖8是基本一致的。在入水初期，如圖 9(a)所示，3D與2D的壓強(qiáng)系數(shù)峰值差別較大，這是由于此時(shí)的后機(jī)身大部分已經(jīng)浸沒(méi)水面，縱向流動(dòng)的影響很大使得3D的壓強(qiáng)系數(shù)峰值遠(yuǎn)小于2D。

(a) t=0.34 s

(b) t=0.4 s

(c) t=0.6 s圖9　入水過(guò)程(β=5°，α=12°，V=0.5 m/s)中機(jī)身x1截面和2D物體表面壓強(qiáng)系數(shù)對(duì)比

3.2 入水速度、姿態(tài)角和尾翹角的影響

圖10是入水過(guò)程(β=5°，α=12°，V=0.5 m/s)中機(jī)身沖擊力Fi系數(shù)的時(shí)間歷程。機(jī)身剛進(jìn)入水面時(shí)(Vt/D<0.07)，沖擊壓強(qiáng)很大，同時(shí)浸潤(rùn)面積迅速增加，導(dǎo)致沖擊力系數(shù)迅速增大；隨著入水深度的增加(0.070.4)，雖然浸潤(rùn)面積依然在增加，但沖擊壓強(qiáng)變小，使得沖擊力系數(shù)小幅震蕩。

可以看出，中央翼盒沖擊水面前后沖擊力系數(shù)發(fā)生了很大的變化，因而本文以中央翼盒沖擊水面的時(shí)刻為界，將機(jī)身最低點(diǎn)觸碰水面的時(shí)刻到中央翼盒觸碰水面的時(shí)刻稱為入水前期，之后為入水后期。

圖10　入水過(guò)程(β=5°，α=12°，V=0.5 m/s)中機(jī)身沖擊力系數(shù)時(shí)間歷程

3.2.1 入水速度

圖11為不同入水速度下(β=5°，α=12°，V=0.25 m/s，0.5 m/s)機(jī)身沖擊力系數(shù)的時(shí)間歷程對(duì)比。在同一無(wú)量時(shí)刻，入水速度較小時(shí)的浸潤(rùn)面積較大，導(dǎo)致沖擊力系數(shù)較大。

圖11　不同入水速度下(β=5°，α=12°，V=0.25 m/s，0.5 m/s)機(jī)身沖擊力系數(shù)的時(shí)間歷程

3.2.2 姿態(tài)角

圖12為不同姿態(tài)角下(β=5°，α=8°,12°，V=0.5 m/s)機(jī)身沖擊力系數(shù)的時(shí)間歷程對(duì)比。入水初期，姿態(tài)角α=8°的沖擊力系數(shù)比姿態(tài)角α=12°的略大，這是由于α=8°的浸潤(rùn)面積略大；在入水后期，由于機(jī)身以α=8°入水時(shí)中央翼盒沖擊水面的時(shí)刻早于以α=12°入水的時(shí)刻，而且沖擊水面時(shí)的底部抬升角較小，這導(dǎo)致以α=8°入水的沖擊力系數(shù)峰值大于以α=12°入水的峰值。

圖12　為不同姿態(tài)角下(β=5°，α=8°,12°，V=0.5 m/s)機(jī)身沖擊力系數(shù)的時(shí)間歷程

3.2.3 尾翹角

圖13為不同尾翹角下(β=3°，5°，7°，α=12°，V=0.5 m/s)機(jī)身沖擊力系數(shù)的時(shí)間歷程對(duì)比。入水初期，尾翹角大的機(jī)身浸潤(rùn)面積較大，因而沖擊力系數(shù)較大；入水后期，尾翹角大的機(jī)身的中央翼盒先沖擊水面導(dǎo)致沖擊力先變大，但由于沖擊時(shí)的底部抬升角差別不大，因而沖擊力系數(shù)峰值差別較小。

圖13　不同尾翹角下(β=3°，5°，7°，α=12°，V=0.5 m/s)機(jī)身沖擊力系數(shù)的時(shí)間歷程

4　結(jié)論

本文通過(guò)數(shù)值模擬方法研究了不同尾翹角的機(jī)身以不同速度和不同姿態(tài)角的入水過(guò)程。結(jié)果分析表明：

1) 機(jī)身入水過(guò)程中壓強(qiáng)峰值首先出現(xiàn)在噴濺根部，隨后轉(zhuǎn)移至機(jī)身底部。三維效應(yīng)的影響表現(xiàn)為入水過(guò)程中機(jī)身截面的壓強(qiáng)峰值先小于二維的峰值，隨后大于二維峰值，最后又小于二維峰值。

2) 入水初期機(jī)身沖擊力系數(shù)迅速增大，而后略有回落，入水后期由于中央翼盒沖擊水面會(huì)導(dǎo)致沖擊力系數(shù)再次迅速增大，而后小幅震蕩。

3) 速度越大、姿態(tài)角越大、尾翹角越小，機(jī)身沖擊力系數(shù)越小。