基于流固耦合的實(shí)體元空心葉片鳥撞數(shù)值模擬

劉建明 ，蔣向華，武 卉 ，馬永峰 ，王 東

(1.北京航空航天大學(xué)能源與動(dòng)力工程學(xué)院，北京 100191；2.中航工業(yè)沈陽發(fā)動(dòng)機(jī)設(shè)計(jì)研究所，沈陽 110015)

1 引言

鳥撞擊葉片是航空發(fā)動(dòng)機(jī)較為普遍的1種外物損傷事故。對(duì)此，發(fā)動(dòng)機(jī)研制技術(shù)先進(jìn)國家開展了許多相關(guān)研究，并制定了相應(yīng)的抗鳥撞的設(shè)計(jì)規(guī)范和試驗(yàn)規(guī)范;中國也對(duì)飛機(jī)抗鳥撞試驗(yàn)提出了明確要求[1]。

鳥撞研究可分為理論、試驗(yàn)和數(shù)值模擬研究。理論研究包括鳥體切割模型、鳥撞物理過程、鳥體載荷模型等[2]；試驗(yàn)研究包括國外進(jìn)行的對(duì)葉片模型件進(jìn)行投鳥試驗(yàn)[3]，國內(nèi)采用高壓空氣炮及高速攝影等模擬鳥體對(duì)飛機(jī)進(jìn)行撞擊試驗(yàn);數(shù)值模擬研究主要包括各種算法及各種有限元程序的開發(fā)與應(yīng)用。目前，國內(nèi)的研究院所開展了較多鳥撞發(fā)動(dòng)機(jī)葉片的數(shù)值模擬研究，如用ANSYS/DYNA軟件，采用接觸沖擊算法進(jìn)行鳥撞葉片數(shù)值模擬[4]；用MSC.DYTRAN軟件，采用流固耦合算法進(jìn)行鳥撞葉片研究[1]等。在采用流固耦合算法進(jìn)行葉片鳥撞研究中，普遍采用殼單元?jiǎng)澐秩~片，而殼單元在模擬復(fù)雜幾何結(jié)構(gòu)方面存在困難。如果采用實(shí)體元?jiǎng)澐秩~片，由于在失效單元的處理及耦合面的推進(jìn)等方面存在諸多困難，而且目前流固耦合算法還缺乏完善的理論指導(dǎo),所以，在現(xiàn)有計(jì)算程序中一般只計(jì)算殼單元的耦合及失效。這就給具有失效模式的實(shí)體元?jiǎng)澐值目招娜~片的鳥撞數(shù)值模擬帶來了很大挑戰(zhàn)。

本文利用PATRAN軟件建立了鳥撞6面體元?jiǎng)澐值钠桨迦~片的有限元模型，研究了不考慮和考慮平板葉片失效2種情況下的數(shù)值模擬方法；建立了鳥撞6面體元?jiǎng)澐值目招娜~片轉(zhuǎn)子級(jí)瞬態(tài)動(dòng)力學(xué)有限元模型，研究鳥體密度、長度、半徑、速度的改變對(duì)葉片瞬態(tài)響應(yīng)的影響；最后進(jìn)行了具有失效模式的空心葉片鳥撞數(shù)值模擬。

2 鳥撞平板葉片方法

2.1 不考慮平板失效的方法

為了驗(yàn)證本文所采用方法的可行性，建立了鳥撞鋁板模型。鋁板尺寸為410 mm×500 mm×10 mm，4邊固支，用分段線性塑性材料模擬，材料參數(shù)見表1。

表1 鋁板材料參數(shù)

鳥體采用線性流體模擬，其質(zhì)量M=1.8 kg，密度ρ=928.15 kg/m3，體積模量 K=2200 MPa，用兩端半球中間圓柱實(shí)體來模擬，以30°角斜撞向鋁板中央。

鋁板用8節(jié)點(diǎn)6面體元?jiǎng)澐?，沿厚度方向劃?層單元；以鋁板體單元的上、下自由表面的節(jié)點(diǎn)為公共節(jié)點(diǎn)，分別建立了1層厚度為0.001 mm的4節(jié)點(diǎn)4邊形殼元；在板的厚度方向的4個(gè)面上建立了1圈啞元，如圖1所示。啞元和殼元形成包圍體單元的封閉的耦合面，如圖2所示。

為了驗(yàn)證鋪薄殼元對(duì)計(jì)算結(jié)果的影響，計(jì)算了體單元表面未鋪和鋪薄殼元時(shí)板中心點(diǎn)的位移隨時(shí)間的變化情況，如圖3所示。從圖中可見，位移響應(yīng)曲線幾乎重合，鋪薄殼元對(duì)計(jì)算結(jié)果的影響很小。

2.2 考慮平板失效的方法

建立了1只鳥體撞擊2塊平板葉片的有限元模型，此計(jì)算是為模擬鳥體撞擊破壞第1塊板后撞擊第2塊板的情形。參照第2.1節(jié)，2塊板表面由殼單元和啞元分別形成1個(gè)包圍平板葉片的耦合面；在上板底端和下板上端之間建立啞元，此啞元和上板下端的殼元以及下板上端的殼元形成1個(gè)耦合面。鳥體數(shù)據(jù)與第2.1節(jié)的相同，平板葉片材料與第2.1節(jié)相比只添加了最大塑性應(yīng)變破壞準(zhǔn)則。定義2塊板的最大塑性應(yīng)變均為0.01，殼的最大塑性應(yīng)變?yōu)?.02，計(jì)算得到的鳥撞瞬間2塊板應(yīng)力云圖如圖4所示。從圖中可見，2塊板被撞壞。

由2塊平板葉片的鳥撞破壞模擬可知，要進(jìn)行具有破壞模式、實(shí)體元?jiǎng)澐值目招娜~片的鳥撞數(shù)值模擬，可以采用在空心葉片表面鋪殼元的方法。

3 不考慮失效的葉片鳥撞模擬

3.1 有限元模型

8節(jié)點(diǎn)6面體空心葉片模型采用PATRAN建模。所建空心葉片模型每個(gè)葉片沿葉高有14層節(jié)點(diǎn)，共13層單元。為了節(jié)省計(jì)算時(shí)間，模型只建立了在發(fā)動(dòng)機(jī)中與鳥體有初始碰撞行為的第1級(jí)葉片，由13片葉片組成，每片葉片有節(jié)點(diǎn)2044個(gè)、單元1122個(gè)；為了進(jìn)行流固耦合計(jì)算，以葉片體單元自由表面上的節(jié)點(diǎn)為公共節(jié)點(diǎn)，建立覆蓋在體表面的4節(jié)點(diǎn)的殼元，其厚度為0.001 mm。所建有限元模型如圖5所示。

3.2 相關(guān)數(shù)據(jù)

葉片材料的靜態(tài)應(yīng)變強(qiáng)化效應(yīng)采用雙線性應(yīng)力-應(yīng)變模型模擬，應(yīng)力循環(huán)效應(yīng)采用等向強(qiáng)化模型模擬，動(dòng)態(tài)應(yīng)變率效應(yīng)采用Cowper-Symonds模型模擬。其硬化方程為

式中：σd為動(dòng)態(tài)應(yīng)力；σy為靜態(tài)屈服應(yīng)力；ε˙為等效應(yīng)變率；D、P為表征動(dòng)態(tài)硬化的材料系數(shù)，計(jì)算時(shí)取D=100、P=10。

葉片的材料密度ρ=4430 kg/m3，泊松比μ=0.34，彈性模量E=112500MPa，靜態(tài)屈服應(yīng)力σy=900MPa，靜態(tài)硬化模量EH=14286MPa。葉片的轉(zhuǎn)速為12000 r/min。

鳥體的幾何尺寸（包括半徑和長度）和初始速度見表2?？紤]到鳥體在碰撞過程中的流變假設(shè)，即當(dāng)鳥體相對(duì)速度大于50 m/s時(shí)，會(huì)發(fā)生流變；鳥撞一般發(fā)生在飛機(jī)起降階段，飛機(jī)與鳥的相對(duì)速度不大。因此，選取了2種速度（70、100 m/s）進(jìn)行對(duì)比計(jì)算。

表2 鳥體的幾何尺寸及速度

鳥體采用線性流體模型，選取3種材料計(jì)算，見表3。

表3 鳥體材料

3.3 葉片的預(yù)應(yīng)力分析

鳥撞事故發(fā)生時(shí)，發(fā)動(dòng)機(jī)一般處于穩(wěn)定旋轉(zhuǎn)狀態(tài)，葉片由于旋轉(zhuǎn)引起的慣性力已經(jīng)具有一定應(yīng)力，為了得到更準(zhǔn)確的分析結(jié)果，必須考慮預(yù)應(yīng)力的影響，在瞬態(tài)分析開始前給葉片模型加上一定的應(yīng)力場。這部分由NASTRAN計(jì)算得到。

3.4 鳥撞數(shù)值模擬

選取第2.2節(jié)中給定材料的葉片、鳥體密度1、鳥體半徑1、鳥體長度1、鳥體速度1進(jìn)行計(jì)算，并以該計(jì)算工況作為1個(gè)標(biāo)準(zhǔn)工況。計(jì)算時(shí)間為0.002 s。

3.4.1 鳥體密度對(duì)葉片應(yīng)力的影響

選取鳥體密度2、3進(jìn)行對(duì)比計(jì)算。在3種密度鳥體撞擊下葉片的應(yīng)力（極大值應(yīng)力處，下同）及其峰值變化曲線分別如圖6、7所示。從圖6中可見，在遭受鳥體撞擊后，葉片會(huì)在較短時(shí)間內(nèi)產(chǎn)生應(yīng)力值較大的沖擊波，其峰值隨著鳥體密度的增大而增大，這是因?yàn)轼B體高速撞擊使葉片發(fā)生應(yīng)變硬化及應(yīng)變率硬化，硬化程度隨鳥體密度增大而增大，葉片的應(yīng)力隨著其硬化程度的增大而增大；從圖7中可見，鳥體密度對(duì)葉片應(yīng)力峰值的影響不是線性的，因?yàn)楫?dāng)鳥體密度增大到一定程度后，由于葉片的可變形能力和硬化能力已充分發(fā)揮，葉片的應(yīng)力峰值的增加速度就會(huì)變慢。

圖7 不同密度鳥體撞擊下的葉片應(yīng)力峰值曲線

3.4.2 鳥體長度對(duì)葉片應(yīng)力的影響

選取鳥體長度2、3進(jìn)行對(duì)比計(jì)算。3種長度鳥體撞擊下葉片的應(yīng)力及其峰值曲線分別如圖8、9所示。從圖中可見，鳥體長度的影響與鳥體密度的影響相類似。鳥體長度對(duì)葉片應(yīng)力峰值的影響不大，這是因?yàn)槿~片應(yīng)力峰值的產(chǎn)生主要來自鳥體前端與葉片發(fā)生碰撞的部位，由于應(yīng)力波的作用，鳥體后端的速度明顯減弱，對(duì)葉片碰撞的影響也就大大減弱。

圖8 不同長度鳥體撞擊下的葉片應(yīng)力曲線

圖9 不同長度鳥體撞擊下的葉片應(yīng)力峰值曲線

3.4.3 鳥體半徑對(duì)葉片應(yīng)力的影響

選取鳥體半徑2、3進(jìn)行對(duì)比計(jì)算。3種半徑鳥撞下葉片應(yīng)力及其峰值變化曲線分別如圖10、11所示。從圖中可見，鳥體半徑對(duì)葉片應(yīng)力的影響與密度的影響相類似。

圖10 不同半徑鳥體撞擊下的葉片應(yīng)力曲線

圖11 不同直徑鳥體撞擊下的葉片應(yīng)力峰值曲線

3.4.4 鳥體速度對(duì)葉片應(yīng)力的影響

選取鳥體速度1、2進(jìn)行對(duì)比計(jì)算。2種鳥體速度下的葉片應(yīng)力及其峰值變化曲線分別如圖12、13所示。從圖中可見，鳥體速度變化對(duì)葉片應(yīng)力峰值有較大影響，因?yàn)轼B體的速度越大，鳥體撞擊時(shí)動(dòng)量也越大，并且撞擊的作用時(shí)間更短，因而產(chǎn)生的應(yīng)力峰值也越大。

圖12 不同速度鳥體撞擊下的葉片應(yīng)力曲線

圖13 不同速度鳥體撞擊下的葉片應(yīng)力峰值曲線

4 考慮失效的葉片鳥撞模擬

為了模擬帶失效模式空心葉片的鳥撞響應(yīng)，給葉片添加了最大塑性應(yīng)變破壞準(zhǔn)則。鳥撞瞬間葉片的破壞情況如圖14所示。

5 結(jié)論

（1）采用在實(shí)體元平板葉片表面鋪1層薄殼元以形成封閉耦合面的方法，驗(yàn)證了加殼單元后對(duì)鳥撞計(jì)算結(jié)果的影響。計(jì)算結(jié)果表明，鋪了1層薄殼單元對(duì)計(jì)算結(jié)果影響很小。

（2）建立了帶失效模式的實(shí)體元平板葉片的鳥撞模型，并進(jìn)行了計(jì)算。計(jì)算結(jié)果比較符合實(shí)際破壞情況，從而為帶失效模式空心葉片的鳥撞模擬提供了方法和依據(jù)。

（3）建立了鳥撞實(shí)體元?jiǎng)澐值目招娜~片轉(zhuǎn)子級(jí)瞬態(tài)動(dòng)力學(xué)有限元模型，模擬了葉片在穩(wěn)定旋轉(zhuǎn)狀態(tài)下遭受不同密度、長度、半徑及速度鳥體撞擊下的葉片的瞬態(tài)響應(yīng)。上述參數(shù)的增加均會(huì)增大葉片的應(yīng)力峰值，并且對(duì)應(yīng)力峰值的影響都不是線性的。

對(duì)鳥撞的數(shù)值模擬可為鳥撞試驗(yàn)和抗鳥撞設(shè)計(jì)提供一定參考。

[1]蔣向華,王延榮.采用流固耦合方法的整級(jí)葉片鳥撞擊數(shù)值模擬[J].航空動(dòng)力學(xué)報(bào),2008,23（2）:299-304.

[2]關(guān)玉璞,陳偉,高德平.航空發(fā)動(dòng)機(jī)葉片外物損傷研究現(xiàn)狀[J].航空學(xué)報(bào),2007,28（4）:851-857.

[3]Horsley J.The Rolls-Royce Way of Validating Fan Integrity[R].AIAA-93-2602.

[4]陳偉,關(guān)玉璞,高德平.發(fā)動(dòng)機(jī)葉片鳥撞擊瞬態(tài)響應(yīng)的數(shù)值模擬[J].航空學(xué)報(bào)，2003,24（6）:531-533.

[5]卞文杰,等.瞬態(tài)動(dòng)力學(xué)CAE解決方案MSC.Dytran基礎(chǔ)教教程[M].北京：北京大學(xué)出版社,2004.