基于FEMSPH算法的柔性直升機(jī)著水?dāng)?shù)值模擬

（航空工業(yè)直升機(jī)設(shè)計(jì)研究所，江西景德鎮(zhèn) 333001）

0 引言

由于直升機(jī)的應(yīng)急救生設(shè)備決定了在發(fā)生事故進(jìn)行水上迫降時(shí)，人員能否安全逃生，近年來越來越受到重視[1]。目前國內(nèi)對(duì)直升機(jī)水上迫降的研究[2-3]還處于以模型試驗(yàn)為主，仿真分析輔助的階段，而國外早已采取仿真分析取代試驗(yàn)的手段。直升機(jī)著水是一個(gè)復(fù)雜的流固耦合問題，其本質(zhì)是結(jié)構(gòu)物入水后所受到的沖擊與滑水。目前很難模擬出柔性直升機(jī)以及應(yīng)急浮筒的運(yùn)動(dòng)歷程。本文旨在提供一種柔性直升機(jī)著水?dāng)?shù)值模擬方法，為早日實(shí)現(xiàn)仿真分析取代試驗(yàn)提供一種思路。

1 計(jì)算方法

本文采用有限元分析（FEM）以及光滑粒子流體動(dòng)力學(xué)（SPH）耦合的方法來模擬直升機(jī)著水的整個(gè)過程。

光滑粒子動(dòng)力學(xué)（SPH）的原理[4]是用一系列任意分布的粒子來表示問題域，粒子的速度、密度和位置受連續(xù)性方程、動(dòng)量守恒方程以及能量守恒方程控制。用積分表達(dá)式對(duì)場函數(shù)進(jìn)行近似，然后采用域內(nèi)的相鄰粒子對(duì)應(yīng)的值疊加求和取代場函數(shù)的積分表達(dá)式，然后將粒子近似法應(yīng)用于所有偏微分方程組的常函數(shù)相關(guān)項(xiàng)中，將偏微分方程組進(jìn)行離散得到粒子的場變量。SPH方法是一種具有無網(wǎng)格、自適應(yīng)屬性的純拉格朗日流體動(dòng)力學(xué)求解方法。

有限元分析（FEM）的原理是將一個(gè)連續(xù)體求解域離散成若干個(gè)單元并通過它們邊界上的節(jié)點(diǎn)相互連接，對(duì)每一個(gè)單元假定一個(gè)合適的近似解從而推導(dǎo)出求解域的近似解。針對(duì)傳統(tǒng)的有限元方法不能模擬大變形問題的缺陷，采用有限元分析（FEM）以及光滑粒子流體動(dòng)力學(xué)（SPH）耦合的方法[5]可以避免自由面大變形所造成的網(wǎng)格畸變問題，并且能準(zhǔn)確地模擬出自由液面的變化以及水花的飛濺。

2 計(jì)算模型的建立

2.1 網(wǎng)格模型

首先對(duì)機(jī)身及浮筒組件進(jìn)行網(wǎng)格繪制，確保所繪制的有限元網(wǎng)格沒有多余的邊界且網(wǎng)格的法向方向的一致。如圖1所示為某型機(jī)的機(jī)身及浮筒組件模型，機(jī)身與浮筒之間通過綁帶連接，機(jī)身及浮筒組件的最小網(wǎng)格尺寸很大程度上決定了整個(gè)仿真的速度，本文選用的機(jī)身及浮筒組件的最小網(wǎng)格尺寸為0.05m。

圖1 機(jī)身及浮筒組件模型

選取長35m，寬30m，深5m的水體進(jìn)行網(wǎng)格繪制，原則上水體網(wǎng)格尺寸應(yīng)略大于機(jī)體及浮筒組件網(wǎng)格尺寸，本文選用的水體網(wǎng)格最小尺寸為0.0625m，將繪制完成的網(wǎng)格轉(zhuǎn)換成SPH粒子，如圖2。

圖2 水體SPH粒子

2.2 材料模型

由于機(jī)身在著水過程中形變較小，為了減小計(jì)算量，將機(jī)體定義為不易變形的剛體材料。而浮筒以及綁帶在著水過程中形變較大，為準(zhǔn)確模擬其形變，將它們定義為膜單元進(jìn)行建模。膜單元的特性為只能受壓，不能受拉，與浮筒以及綁帶的特性相吻合。水體在直升機(jī)沖擊作用下發(fā)生劇烈的變形、移動(dòng)，故水體采用SPH粒子建模，其中，粒子的粘性和光滑半徑將關(guān)系到能否準(zhǔn)確的模擬出水體的狀態(tài)。

2.3 計(jì)算條件的定義

計(jì)算條件的定義主要在于接觸模型的定義。接觸模型的建立，使得多個(gè)拓?fù)洳弧跋噙B”的結(jié)構(gòu)可以傳遞載荷。直升機(jī)著水涉及多個(gè)接觸關(guān)系，包括水體與機(jī)身、浮筒與機(jī)身、浮筒綁帶與機(jī)身及浮筒、浮筒自接觸等。以浮筒與機(jī)身的接觸為例，它是用來計(jì)算機(jī)身和浮筒的運(yùn)動(dòng)關(guān)系和傳力關(guān)系。在著水過程中浮筒受水體沖擊，會(huì)向機(jī)身靠近，并擠壓機(jī)身。建立機(jī)身和浮筒的接觸，可以不斷檢測計(jì)算浮筒和機(jī)身的幾何空間關(guān)系，并通過接近程度計(jì)算接觸力，接觸力作為作用力和反作用力施加給接觸對(duì)雙方的節(jié)點(diǎn)。接觸厚度的大小設(shè)定和模型之間的初始距離關(guān)系需要嚴(yán)密控制，即需要使接觸厚度小于接觸雙方模型的初始距離，否則會(huì)產(chǎn)生初始穿透，在整個(gè)時(shí)間歷程中可能存在的接觸關(guān)系均需要考慮。

其次是載荷條件的設(shè)置，將水體粒子、機(jī)身及浮筒組件賦予相應(yīng)的初始速度、重力加速度等，然后再進(jìn)行求解器設(shè)置，設(shè)定時(shí)間、控制迭代步長等，最后根據(jù)實(shí)際需求輸出機(jī)身及浮筒組件重心處的加速度，獲得直升機(jī)著水后的動(dòng)態(tài)響應(yīng)，確定最佳入水姿態(tài)角；輸出浮筒載荷、浮筒對(duì)機(jī)體表面的壓力等參數(shù)用于部件結(jié)構(gòu)設(shè)計(jì)及強(qiáng)度校核。

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

計(jì)算模型建立后，選取6°俯仰角，15.43m/s水平速度，1.5m/s垂向速度這一工況進(jìn)行數(shù)值模擬。圖3為這一工況下直升機(jī)的橫向及垂向過載隨時(shí)間的變化。圖4為這一工況下直升機(jī)的橫向及垂向速度隨時(shí)間的變化。

圖3 直升機(jī)的橫向及垂向過載隨時(shí)間的變化

如圖3、圖4所示，直升機(jī)在觸水初期垂向過載加速度仍接近重力加速度，這主要是因?yàn)楦⊥才c機(jī)體之間的連接是通過綁帶的柔性連接，后浮筒觸水后，后浮筒向上運(yùn)動(dòng)，機(jī)體仍保持自由落體運(yùn)動(dòng)，浮筒與機(jī)體之間的相對(duì)位置發(fā)生了變化。在0.06s時(shí)，機(jī)體底部開始觸水，機(jī)體受到了水體的垂向沖擊力，垂向過載逐漸增加，垂向速度仍在增加，增加幅度減小，當(dāng)垂向過載大于1時(shí)，垂向速度開始逐漸減小。在0.17s時(shí)，垂向過載達(dá)到最大值3.56g，隨后垂向過載逐漸減小，垂向速度逐漸趨于零。在0.37s時(shí)，垂向過載小于1，垂向速度開始增加，機(jī)體在逐漸震蕩中趨于穩(wěn)定。在整個(gè)著水歷程中，直升機(jī)的橫向過載先增大后減小，在0.17s時(shí)達(dá)到最大值0.44g，直升機(jī)的橫向速度在不斷減小，直至趨于穩(wěn)定。

圖4 直升機(jī)的橫向及垂向速度隨時(shí)間的變化

4 結(jié)論

（1）直升機(jī)柔性浮筒著水?dāng)?shù)值模擬直升機(jī)入水過程中的速度、過載峰值變化趨勢與理論分析保持一致。

（2）直升機(jī)以某一垂向速度著水時(shí)，實(shí)際觸水時(shí)垂向速度大于該垂向速度。

（3）直升機(jī)以6°俯仰角，15.43m/s水平速度，1.5m/s垂向速度在靜水中著水時(shí)垂向過載最大達(dá)到3.56g。