基于CEL算法的水陸兩棲飛機水上降落動力特性分析*

姚小虎 黃愉太 歐智成 劉曉明 胡伯仁

(1.華南理工大學(xué) 土木與交通學(xué)院，廣東 廣州510640;2.成都飛機工業(yè)(集團)有限責(zé)任公司 技術(shù)中心，四川成都610092)

水陸兩棲飛機是一種能在水面上起飛和降落的特殊多用途飛機［1］，如水上轟炸機、水上偵查機、水上運輸機、水上戰(zhàn)斗機、海上救護機、滅火用水上飛機等，隨著我國航空事業(yè)的快速發(fā)展，水陸兩棲飛機也在飛機大家族中占據(jù)一席之地.

水陸兩棲飛機的水上降落屬于典型的流固耦合問題［2］，入水時水對機身船體的沖擊使船體發(fā)生變形，而船體變形反過來又將影響水的流動.入水問題的研究早期以實驗［3-4］和理論［5］為主，隨著計算機的高速發(fā)展，數(shù)值模擬也越來越多地應(yīng)用于入水問題.Brooks 等［6］最先使用LS-DYNA 模擬Apollo 指令艙撞擊水面的過程. 盧鳴飛等［7］采用FEM-SPH 算法對飛機的入水進行了分析. 劉翔［8］建立了A320飛機模型，采用ALE/Penalty Coupling 方法成功模擬出飛機迫降的流固耦合作用，并將計算結(jié)果與實驗進行了對比.張韜等［9］運用Dytran、采用一般耦合算法計算了飛機水上迫降的沖擊荷載，并將結(jié)果與剛體飛機模型進行了對比.Chandra 等［10］對2009年在美國發(fā)生的一起飛機成功水上迫降的例子采用FEM-SPH 算法進行了分析，分析中將飛機視為剛體，結(jié)果證明采用該方法進行水上迫降計算是合理的.但是目前對于水上降落問題研究的大部分文獻都是將機身處理成剛體［11］，且只給出定性結(jié)果而沒有做定量的研究，而在實際入水過程中，考慮機身的沖擊變形對提高結(jié)果的精確度是非常有必要的［12］.

文中采用耦合的歐拉－拉格朗日(CEL)算法對水上飛機降落的動力過程進行了研究，分析了某型水陸兩棲飛機機頭水上降落過程中不同姿態(tài)角和垂直速度對飛機結(jié)構(gòu)動力響應(yīng)的影響，并得到適合該機型水上降落的參考角度和速度，為該型水陸兩棲飛機水上降落提供理論指導(dǎo)依據(jù).

1 耦合歐拉－拉格朗日法

拉格朗日(Lagrangian)方法主要應(yīng)用于固體力學(xué)中的受力和變形分析，在計算過程中，所分析物體的形狀變化與有限元網(wǎng)格的變化完全保持一致，其優(yōu)點在于它可以大大簡化控制方程的求解過程，并可跟蹤質(zhì)點的運動軌跡，從而能準確地描述物體邊界的運動;但是在涉及特大變形問題時，拉格朗日算法采用的單元網(wǎng)格將會出現(xiàn)嚴重的畸變現(xiàn)象，以致無法繼續(xù)計算.歐拉方法(Eulerian)主要應(yīng)用于流體力學(xué)的計算分析，其采用的有限元網(wǎng)格形狀、大小和空間位置保持不變，在計算過程中材料可在固定的有限元網(wǎng)格中流動，因此歐拉算法能夠處理物質(zhì)的扭曲及一些大變形問題，但歐拉方法在捕捉物體邊界信息上較為困難，不能夠精確描述物質(zhì)的邊界［13］.

基于以上方法特點，Noh［14］首先提出了耦合的歐拉－拉格朗日算法(CEL)，它有機地結(jié)合拉格朗日方法和歐拉方法的優(yōu)點，克服其中的缺陷，充分發(fā)揮各自的優(yōu)勢. 在CEL 方法中，固體變形小則采用拉格朗日單元，而流體變形大時采用歐拉單元. 因此，CEL 算法在解決物體的大位移時有較大優(yōu)勢，比如碰撞、流體動力學(xué)及流固耦合作用等.

2 水上飛機機頭降落動力學(xué)模型

2.1 飛機機頭入水有限元模型

機頭包括雷達罩、駕駛艙、生活艙、船體4 部分，其中船體為入水首要沖擊區(qū)域.建立機頭入水模型，如圖1 所示.機頭蒙皮和腹板均采用殼單元，船體處加強件采用梁單元. 殼單元主要采用四邊形和三角形單元，船體底部單元局部加密.水域模型為橢圓柱體，其長軸為9 m，短軸為6 m，高3 m，單元類型為六面體的歐拉單元，考慮到濺水，歐拉單元的區(qū)域要略大于初始水域，因此初始水域深2.7 m，歐拉單元表面在初始液面上0.3 m.開始時船底與水域相切，在入水位置附近加密網(wǎng)格，如圖2 所示.

圖1 機頭模型Fig.1 Aircraft head model

2.2 材料本構(gòu)

機頭各部件均采用鋁板，材料分布及物理參數(shù)如表1 所示，各型號鋁材采用彈塑性雙線性本構(gòu)模型進行描述:

圖2 機頭入水模型Fig.2 Model of the head hitting the water

式中，σ 和ε 為應(yīng)力和應(yīng)變，σ0和ε0為屈服應(yīng)力和屈服應(yīng)變，E 為彈性模量，ET為強化模量.

表1 機頭各部件材料物理參數(shù)1)Table 1 The material physical parameters of head

受沖擊流體的材料行為極其復(fù)雜，可使用狀態(tài)方程來描述材料的沖擊力、體積和能量特性.各向同性沖擊力可由Mie-Gruneisen 狀態(tài)方程［15］得到，即

式中:Im為單位質(zhì)量的內(nèi)能;η =1 －ρ/ρ0，為名義體積壓縮率，ρ0表示初始材料密度;Г0為Gruneisen 常數(shù).激波速度Us與粒子速度Up之間有線性關(guān)系，Us=c+Up.在Us-UpHugoniot 模型中有PH=ρ0c2η·(1 －sη)－2，其中s 為極限的體積壓縮率，c 為材料中的聲速，當(dāng)名義應(yīng)變較小時，體積模量K=ρ0c2.

文中采用線性的Us-UpHugoniot 狀態(tài)方程描述水的特性，取參數(shù)s=Г0=0，聲速c=1483 m/s，密度ρ=1000 kg/m3，黏度μ=0.001 kg/(m·s).

3 機頭入水結(jié)果分析

文中考慮機頭在不同的入水角度和垂直入水速度下的落水情況，共5 種工況，如表2 所示.

表2 機頭入水工況Table 2 Cases of head attacking the water

飛機降落時受到豎直向上的氣動力2/3 g，與自重疊加，可近似認為飛機降落時有向下的1/3 g 加速度.水域施加一般重力場.水域邊界為自由流入和無反射流出.所有工況的機頭水平速度均為41 m/s，計算時間為0.1 s.

3.1 工況1 結(jié)果分析

圖3 是工況1(入水角為7°，垂直入水速度為3m/s)下機頭應(yīng)力和水域在入水過程中的變化.從云圖可看出，機頭入水時受力區(qū)域主要為船體，機身部分基本無影響，因此重點分析船體的結(jié)構(gòu)動力響應(yīng).機頭從0s 開始與水面接觸，在入水后的0.016 s 出現(xiàn)最大應(yīng)力，為550MPa，位于船體骨架的縱向腹板處，超過材料屈服強度且接近抗拉強度，說明縱向腹板受力明顯，危險性較高，機頭其他部位受力相對較小，均未超其屈服強度.隨著入水深度增加，機頭船底的受沖擊區(qū)域逐漸增大，入水過程中機頭逐漸抬頭.水域被劃出一道水槽，影響區(qū)域逐漸增大，并伴隨液體飛濺.

圖3 工況1 下機頭應(yīng)力和水域的變化Fig.3 Head’s stress and water variation in case 1

3.2 不同垂直速度下入水的結(jié)果分析

當(dāng)飛機入水角為7°，垂直速度分別為3、2、1m/s時，機頭應(yīng)力、船底沖擊力和質(zhì)心速度的變化情況如表3 和圖4 所示.

表3 機頭不同垂直速度入水結(jié)果Table 3 The results of head attacking the water at different vertical speeds

通過對3 種工況的比較可知，機頭應(yīng)力都是在入水初期(0.02 s 左右)出現(xiàn)最大值. 對于船殼最大應(yīng)力，3 種工況相差不大，均保持在430 ～450MPa 之間，滿足材料強度要求.機頭船體骨架應(yīng)力隨垂直速度變化明顯，在3 和2 m/s 時其最大應(yīng)力甚至超過了船殼應(yīng)力，分別為550 和529 MPa，接近材料抗拉強度，1m/s 時骨架應(yīng)力迅速減小為337MPa，遠小于抗拉強度，三者出現(xiàn)位置均位于船體骨架的縱向腹板處，可見腹板受垂直速度影響明顯.隨著垂直速度的減小，船體骨架的應(yīng)力也減小，這是因為隨著入水深度增加，速度越大，船體骨架尤其縱向腹板的受擠壓程度也越大，容易發(fā)生局部屈曲和變形.當(dāng)入水角度相同時，垂直速度越大，船底沖擊力越小. 在入水初期，船底底面法向與水域垂直，此時底面受到的沖擊力出現(xiàn)峰值，而后沖擊力逐漸衰減，隨著機頭慢慢抬頭，受沖擊面法向發(fā)生變化，受水平速度影響，沖擊力再次出現(xiàn)峰值.從質(zhì)心速度減小量可看出，入水過程中機頭整體保持穩(wěn)定(水平晃動速度v2極小)，垂直速度(v3)越大，入水深度也越大，從而阻力增加，水平速度(v1)減小得也越多. 綜合對比3 種工況，當(dāng)入水速度為1 m/s時機頭入水最安全.

3.3 不同角度下入水的結(jié)果分析

飛機在入水的垂直速度為3 m/s，入水角度分別為7°、6°、5°時的結(jié)果如表4 和圖5 所示.

圖4 不同垂直速度下入水的機頭應(yīng)力、船底沖擊力和質(zhì)心速度時程曲線Fig.4 Time-history curves of head stress，bottom pressure and mass centre of the head at different vertical speeds

表4 機頭不同角度入水結(jié)果Table 4 The results of the head attacking the water at different angles

機頭入水初期應(yīng)力迅速增大，在0.02 s 左右出現(xiàn)最大值，并隨入水角度減小而減小.3 種工況下，船殼的最大應(yīng)力相當(dāng)，變化不大，均在430 ～450MPa之間波動，滿足材料強度要求，但對于船底骨架，最大應(yīng)力變化明顯，在7°和6°時超過船殼分別達到了550 和536 MPa，出現(xiàn)在縱向腹板處，接近抗拉強度，為最危險點.在入水角為5°時減小到365 MPa，滿足強度要求.入水瞬間船底沖擊力出現(xiàn)峰值，然后衰減，隨著機頭抬頭，沖擊面發(fā)生變化，沖擊力再次出現(xiàn)峰值，當(dāng)入水角度為6°時，船底的沖擊力最大，這是因為船底為曲面，6°時船殼底面法向方向更接近豎直方向，因此入水的速度對于船底法向方向的分量越大，沖擊力也越大.入水過程中，6°時船底底面與水面更接近平行，產(chǎn)生阻力相對較小，其水平速度(v1)較其他兩工況減小量更少，對于垂直速度(v3)，由于入水深度大，7°時減小明顯，3 種工況下機頭整個過程都保持穩(wěn)定(v2極小). 綜上可得出，當(dāng)飛機水上降落時入水角度為5°的安全性更高.

4 結(jié)論

文中建立了某型水陸兩棲飛機機頭水上降落流固耦合模型，采用CEL 算法分析了機頭在不同入水角度和垂直入水速度下的結(jié)構(gòu)響應(yīng). 船殼的應(yīng)力在入水初期迅速增大，之后處于波動狀態(tài).各工況下船殼最大應(yīng)力值相差不大，而船體骨架受力變化明顯，其最大值大部分情況下甚至超過船殼，且都出現(xiàn)在縱向腹板處，因此在結(jié)構(gòu)設(shè)計中應(yīng)加強骨架尤其是縱向腹板處的結(jié)構(gòu)強度. 機頭入水瞬時船底沖擊力即出現(xiàn)峰值.入水速度減小，機頭應(yīng)力和船底沖擊力都隨之減小;隨著入水角度的減小，機頭應(yīng)力也減小，船底沖擊力在6°時最大. 因此，根據(jù)計算結(jié)果，當(dāng)入水角度為5°、垂直入水速度為1 m/s 時，飛機水上降落相對最安全.

圖5 不同入水角度下機頭應(yīng)力、船底沖擊力和質(zhì)心速度時程曲線Fig.5 Time-history curves of head stress，bottom pressure and mass centre of the head at different attacking angles

［1］尹士達，林明.我國水上飛機的特點和主要性能［J］.江蘇航空，1998(3):12-13.Yin Shi-da，Lin Ming. The main characteristics and performance of seaplane in China ［J］. Jiangsu Aviation，1998(3):12-13.

［2］張阿漫，戴紹仕.流固耦合動力學(xué)［M］.北京:國防工業(yè)出版社，2011:49-85.

［3］Thomas W L.Ditching investigation of a 1/20-scale model of the space shuttle orbiter ［M］. Bethpage:Grumman Aerospace Corporation，1975.

［4］Sawan A Shah. Water impact investigations for aircraft ditching analysis［D］.Melbourne:Royal Melbourne Institute of Technology，2010.

［5］Von Karman T.The impact of seaplane floats during landing［C］∥Proceedings of NACA TN 321，National Advisory Committee for Aeronautics. Washington:［s. n.］，1929:1-8.

［6］Brooks J R，Anderson L A. Dynamics of a space module impacting water ［J］. Spacecraft and Rockets，1994，31(3):509-515.

［7］盧鳴飛，華誠.基于FEM-SPH 耦合算法的飛機水上迫降的力學(xué)建模與分析［J］. 復(fù)旦學(xué)報:自然科學(xué)版，2013，52(5):583-591.Lu Ming-fei，Hua Cheng.Model and numerical simulations of water ditching by FEM-SPH coupling method ［J］.Journal of Fudan University:Natural Science，2013，52(5):583-591.

［8］劉翔.飛機水上迫降的運動特性分析和討論［D］. 武漢:武漢理工大學(xué)理學(xué)院，2012.

［9］張韜，李書，江翔，等.民用飛機水上迫降分析模型和數(shù)值仿真［J］.南京航空航天大學(xué)學(xué)報，2010，42(3):392-394.Zhang Tao，Li Shu，Jiang Xiang，et al.Analysis model and numeral simulation for civil plane ditching［J］.Journal of Nanjing University of Aeronautics ＆ Astronautics，2010，42(3):392-394.

［10］Chandra C，Wong T Y，Bayandor J.Crash worthiness assessment in aircraft ditching incidents［C］∥Proceeding of 27th Congress of International Council of the Aeronautical Sciences.Nice:［s.n.］，2010:1-13.

［11］賀謙，陳效鵬，李磊，等.飛機水面降落的機身載荷研究［J］.強度與環(huán)境，2009，36(4):28-33.He Qian，Chen Xiao-peng，Li Lei，et al.Research on the load of aircraft landed on water［J］. Structure ＆ Environment Engineering，2009，36(4):28-33.

［12］宋長福.民機機身結(jié)構(gòu)入水沖擊問題數(shù)值仿真研究［D］.上海:上海交通大學(xué)航空航天學(xué)院，2011.

［13］徐文杰.基于CEL 算法的滑坡涌浪研究［J］.工程地質(zhì)學(xué)報，2012，20(3):350-354.Xu Wen-jie. CEL algorithm study of reservoir surge induced by landslide［J］.Journal of Engineering Geology，2012，20(3):350-354.

［14］Noh W F.CEL:a time-dependent two-space-dimensional coupled Eulerian-Lagrangian code［M］.New York:Academic Press，1964:117-179.

［15］Shyue Keh-ming.A fluid-mixture type algorithm for compressible multicomponent flow with Mie-Gruneisen equation of state ［J］. Journal of Computational Physics，2001，171(2):678-707.