簡成文 吳遠(yuǎn)飛

（中國直升機(jī)設(shè)計(jì)研究所，江西 景德鎮(zhèn) 333001）

0 引言

飛機(jī)在高空高速飛行，當(dāng)遇到緊急情況時(shí)必須迫降。飛機(jī)的迫降形式主要有3 種：場內(nèi)迫降、場外迫降和水面迫降。其中，水面迫降是指著陸場在海洋、湖泊以及河流等水面上，水面迫降要求盡可能地靠近陸地，主要采用3 種方法（試驗(yàn)法、解析法以及數(shù)值仿真法）對飛機(jī)水上迫降性能進(jìn)行研究。20 世紀(jì)80 年代以前，研究人員主要通過動(dòng)力相似縮比模型試驗(yàn)來研究飛機(jī)的水上迫降性能。然而，隨著計(jì)算機(jī)科學(xué)技術(shù)的發(fā)展，飛機(jī)的全機(jī)模擬仿真以其精度高、成本低以及效率高的優(yōu)點(diǎn)逐漸取代了試驗(yàn)的研究方法。Anghileri M 等[1]對近似剛體結(jié)構(gòu)進(jìn)行入水試驗(yàn)，并用SPH 方法進(jìn)行仿真，對無量綱的加速度和壓力結(jié)果進(jìn)行對比，數(shù)據(jù)吻合度較高。金屬板結(jié)構(gòu)是直升機(jī)結(jié)構(gòu)中最常見的結(jié)構(gòu)，對直升機(jī)典型金屬板結(jié)構(gòu)入水沖擊響應(yīng)進(jìn)行研究，可以為后續(xù)直升機(jī)整機(jī)入水仿真分析奠定基礎(chǔ)。該文采用SPH-FEM 流固耦合計(jì)算方法對金屬板結(jié)構(gòu)入水沖擊響應(yīng)進(jìn)行研究。

1 SPH-FEM 流固耦合計(jì)算方法

SPH粒子與有限單元通過定義接觸的方式進(jìn)行流固耦合計(jì)算。將有限單元的節(jié)點(diǎn)視為背景粒子，背景粒子的變量與相應(yīng)的有限單元節(jié)點(diǎn)一致，例如粒子質(zhì)量、位置、速度和應(yīng)力等[2-3]。背景粒子只能被其他SPH 粒子搜索，在每個(gè)時(shí)間步內(nèi)，接觸部分的相關(guān)信息會(huì)從背景粒子傳遞到有限元數(shù)據(jù)中。對有限元部分來說，這種傳遞相當(dāng)于施加邊界條件，對SPH 部分來說，有限單元節(jié)點(diǎn)轉(zhuǎn)化為鄰近粒子，從而避免了邊界影響，進(jìn)而使SPH 粒子與有限元之間具有很好的連續(xù)性。

位于有限元節(jié)點(diǎn)支持域內(nèi)的SPH 粒子會(huì)對該節(jié)點(diǎn)產(chǎn)生接觸力，同時(shí)位于SPH 粒子支持域內(nèi)的有限元節(jié)點(diǎn)會(huì)對該粒子產(chǎn)生接觸力，接觸勢能如公式（1）所示。式中：NCONT為粒子i的鄰近粒子數(shù)；m為粒子質(zhì)量；ρ為密度；K為接觸罰剛度；rij為粒子間距；W為光滑函數(shù)（當(dāng)xA、xB屬于相同個(gè)體時(shí)，W（xA-xB）=0）；?pavg為粒子間光滑長度的平均值；n為粒子總數(shù)；x為粒子位置。

接觸力如公式（2）所示。

式中：?為拉普拉斯算子。

對SPH 粒子來說，是將接觸力施加到動(dòng)量方程的，如公式（3）所示。

式中：ν為粒子的速度矢量；t為時(shí)間；N為粒子數(shù)量；σ為應(yīng)力；∏ij為人工黏度項(xiàng)；x為密度、速度和能量等變量；α、β為常數(shù)，取值分別為0.04 和0.01；xi是粒子i的位置。

在SPH 方法中，Monoghan 型的人工黏性將動(dòng)能轉(zhuǎn)換為熱能，同時(shí)防止粒子相互接近時(shí)出現(xiàn)非物理穿透現(xiàn)象，如公式（4）所示。

式中：c為聲速；?為粒子的形函數(shù)。

對大多數(shù)流體仿真來說，常數(shù)α∏和β∏的取值分別為0.04 和0.01。

對有限單元來說，接觸力被當(dāng)作外力施加到動(dòng)力學(xué)方程，如公式（5）所示。

式中：、和u分別為加速度、速度和位移；M、C和K分別為系統(tǒng)的質(zhì)量矩陣、阻尼矩陣和剛度矩陣。

SPH-FEM 流固耦合計(jì)算步驟如圖1 所示。

圖1 SPH-FEM 流固耦合方法計(jì)算示意圖

2 金屬板入水沖擊動(dòng)力學(xué)建模

2.1 SPH 水體模型

水體單元模擬采用Langrange實(shí)體單元和SPH粒子單元，當(dāng)采用實(shí)體單元時(shí)，將實(shí)體單元賦予水動(dòng)力（Hydrodynamic）材料，狀態(tài)方程采用多項(xiàng)式狀態(tài)方程，如公式（6）、公式（7）所示。

式中：Ci為材料常數(shù)；μ為無量綱參數(shù)；Ei為內(nèi)能；ρ0為水體的初始密度；ρ為當(dāng)前水體的密度。

當(dāng)水體采用SPH 粒子單元時(shí)，將SPH 粒子單元賦予Monaghan 狀態(tài)方程材料，狀態(tài)方程如公式（8）所示。

式中：p0為水體的初始壓強(qiáng)；ρ0為水體的初始密度；γ為常量，一般來說γ=7；B為狀態(tài)方程中的系數(shù)。

當(dāng)流體最大流速為vmax時(shí)，狀態(tài)方程中的系數(shù)B如公式（9）所示。

考慮試驗(yàn)水池足夠大，不考慮水池大小對仿真分析結(jié)果的影響，建立1 200 mm×1 200 mm×405 mm 的水域，水域由SPH粒子和實(shí)體單元耦合而成。其中，試驗(yàn)件與水體接觸部分采用SPH 粒子模擬，該部分水域尺寸為800 mm×800 mm×100 mm，其余水域采用實(shí)體單元模擬。SPH粒子單元采用立方體形式，粒子間距為15 mm。

2.2 流固耦合接觸模型

SPH 粒子與FEM 單元之間的流固耦合通過接觸來設(shè)置。接觸是邊界條件高度非線性的復(fù)雜問題，需要追蹤接觸前物體的運(yùn)動(dòng)以及接觸后物體之間的相互作用。在接觸過程中，載荷隨時(shí)間、結(jié)構(gòu)的變形而變化，即結(jié)構(gòu)與載荷互相耦合。

接觸計(jì)算采用罰函數(shù)法，在每個(gè)時(shí)間步先檢查從節(jié)點(diǎn)是否穿透主面，如果無穿透，就不做處理；如果穿透，就在該從節(jié)點(diǎn)與被穿透主面之間引入一個(gè)接觸力，其大小與穿透深度、主面剛度成正比，從節(jié)點(diǎn)i受到的接觸力Fi如公式（10）所示。

式中：SLFACM為接觸剛度縮放因子；STF（SNODE）為主/從節(jié)點(diǎn)的接觸剛度；gi為穿透厚度。

同時(shí)，非線性懲罰方式能夠增加罰剛度，從而避免完全穿透。非線性罰剛度如公式（11）所示。

式中：ki為線性剛度；FSVNL為非線性剛度參數(shù)；Hcont為接觸厚度；為非線性罰剛度。

其接觸力如公式（12）所示。

該文對結(jié)構(gòu)與水體之間的接觸進(jìn)行設(shè)置，定義SPH 粒子為從節(jié)點(diǎn)，結(jié)構(gòu)為主面，建立相應(yīng)的接觸模型。

2.3 金屬材料模型

在大變形狀態(tài)下，金屬材料結(jié)構(gòu)的薄弱區(qū)域?qū)⒊霈F(xiàn)應(yīng)力集中的現(xiàn)象，當(dāng)應(yīng)力過大并超過材料彈性極限時(shí)，薄弱區(qū)域的材料將進(jìn)入塑性狀態(tài)，線彈性的本構(gòu)關(guān)系不再適用，需要滿足相應(yīng)的屈服準(zhǔn)則。而屈服準(zhǔn)則是在外部載荷作用下，物體內(nèi)某一點(diǎn)開始產(chǎn)生塑性變形時(shí)應(yīng)力所必須滿足的條件。

金屬材料屬于各向同性材料，該文中金屬的強(qiáng)度準(zhǔn)則采用最大線應(yīng)變準(zhǔn)則，定義如下：當(dāng)材料最大線應(yīng)變ε1（或|ε3|）達(dá)到極限值時(shí)，材料發(fā)生破壞。其表達(dá)式如公式（13）所示。

式中：εT、εC分別為材料拉伸、壓縮時(shí)的極限應(yīng)變。

在仿真計(jì)算中，當(dāng)殼單元應(yīng)變達(dá)到預(yù)期最大應(yīng)變值時(shí)單元失效并刪除單元。

2.4 結(jié)構(gòu)有限元模型

對元組件進(jìn)行簡化，抽取相應(yīng)的框架及面板中面，對結(jié)構(gòu)焊接部分進(jìn)行共節(jié)點(diǎn)處理，對底部試驗(yàn)件的螺栓連接建立相應(yīng)的彈簧元。為了模擬在吊籃頂部施加質(zhì)量塊，建立相應(yīng)的剛體約束，并在剛體中心點(diǎn)施加質(zhì)量點(diǎn)。試驗(yàn)件及組件采用四邊形殼單元模擬，單元尺寸為10 mm×10 mm。

試驗(yàn)裝置結(jié)構(gòu)如圖2 所示，有限元模型如圖3 所示。

圖2 試驗(yàn)裝置結(jié)構(gòu)圖

圖3 有限元模型

3 金屬板入水沖擊試驗(yàn)

3.1 試驗(yàn)條件

試驗(yàn)裝置由水池、試驗(yàn)件、吊籃、導(dǎo)軌、臺架以及釋放裝置組成。試驗(yàn)件安裝在吊籃底部，試驗(yàn)件上布置應(yīng)變片。吊繩與吊籃相連，通過起吊裝置及其控制系統(tǒng)調(diào)整吊籃的起吊和投放高度，吊籃頂部安裝快速釋放裝置，通過相應(yīng)的控制系統(tǒng)實(shí)現(xiàn)快速投放吊籃和試驗(yàn)件的功能。在水池與起吊裝置之間設(shè)置導(dǎo)軌，吊籃與導(dǎo)軌接觸處采用輪子設(shè)計(jì)，吊籃投放后將沿導(dǎo)軌下墜，從而控制試驗(yàn)件的入水姿態(tài)。

吊籃具備可以承受垂向著水沖擊載荷及側(cè)向載荷的能力，吊籃上設(shè)置有配重安裝點(diǎn)，通過設(shè)置配重質(zhì)量及配重位置得到試驗(yàn)件不同質(zhì)量狀態(tài)下的試驗(yàn)條件。

吊籃四周布置加速度傳感器，試驗(yàn)件從規(guī)定高度下落，以相應(yīng)速度墜入下方水池。

3.2 試驗(yàn)參數(shù)

金屬平板為0.8 mm鋁合金薄板，尺寸為400 mm×400 mm。結(jié)構(gòu)總質(zhì)量為50 kg，入水速度為2 m/s。在距薄板中心r=50 mm 處布置應(yīng)變片，吊籃處布置加速度傳感器。

3.3 結(jié)果分析

典型金屬結(jié)構(gòu)著水后的加速度響應(yīng)、應(yīng)變響應(yīng)如圖4、圖5 所示。由圖4 可知，金屬平板入水大約5 ms 后，加速度與應(yīng)變達(dá)到峰值，隨后迅速衰減，并在零點(diǎn)附近做周期振蕩。且加速度響應(yīng)和應(yīng)變響應(yīng)的仿真曲線均與試驗(yàn)曲線基本保持一致。結(jié)果表明，采用SPH-FEM 法進(jìn)行金屬板入水沖擊分析可以真實(shí)地反映結(jié)構(gòu)的沖擊響應(yīng)特性。

圖4 典型金屬結(jié)構(gòu)加速度響應(yīng)

圖5 典型金屬結(jié)構(gòu)距平板中心r=50 mm 處的應(yīng)變響應(yīng)

由圖4 可知，金屬板加速度仿真值為17 g，試驗(yàn)值為17 g，結(jié)果的誤差為0%。由圖5 可知，距薄板中心r=50 mm處應(yīng)變仿真值為0.001 1 ε，試驗(yàn)值為0.001 15 ε，結(jié)果的誤差為5%。上述數(shù)據(jù)表明，該文所用仿真模型可以準(zhǔn)確地描述結(jié)構(gòu)的沖擊響應(yīng)特性。

4 結(jié)語

該文建立了金屬板及其支撐結(jié)構(gòu)、水體以及空氣的有限元模型，采用SPH-FEM 流固耦合方法對直升機(jī)典型金屬結(jié)構(gòu)進(jìn)行入水仿真分析，相關(guān)結(jié)論如下：1）結(jié)構(gòu)加速度仿真值及應(yīng)變仿真值與試驗(yàn)值相比，誤差小于或等于5%，表明該文所用仿真方法具有較高的準(zhǔn)確性。2）結(jié)構(gòu)加速度與應(yīng)變響應(yīng)在結(jié)構(gòu)入水沖擊初期達(dá)到最大，隨后迅速衰減，仿真曲線與試驗(yàn)曲線變化趨勢一致。3）該文所用的仿真模型與仿真方法具有較好的可行性和較高的準(zhǔn)確性，能夠?yàn)楹罄m(xù)直升機(jī)整機(jī)入水仿真分析奠定基礎(chǔ)。