代理模型在機(jī)翼油箱晃動適航審定中的應(yīng)用

孫 濱， 韓 松， 王東林， 李 明*， 繆圣亮

(1.中國民用航空沈陽航空器適航審定中心，沈陽 110043;2.大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室，大連 116024; 3.中航通飛研究院有限公司，珠海 519040)

1 引 言

隨著航空航天技術(shù)發(fā)展，液體晃動引起的飛行器結(jié)構(gòu)破壞逐漸得到重視和關(guān)注[1-3]。在各類航天器飛行過程中，受不穩(wěn)定氣流、發(fā)動機(jī)轉(zhuǎn)子失衡以及飛行軌跡和姿態(tài)的變換等諸多因素影響，液體燃料和冷卻劑在儲液箱內(nèi)會發(fā)生劇烈晃動，導(dǎo)致儲液箱壁面應(yīng)力集中，引發(fā)壁面破壞導(dǎo)致燃料泄漏或爆炸，造成巨大的損失。因此，在飛行器適航審定中，考慮儲油箱低頻晃動引起的壁面結(jié)構(gòu)破壞非常必要。

儲油箱晃動涉及到多因素流固耦合，問題復(fù)雜，全尺寸試驗(yàn)代價(jià)高。翟冰冰等[4]通過對大型運(yùn)輸類飛機(jī)油箱進(jìn)行適航驗(yàn)證，確定了晃動頻率、角度和試驗(yàn)載液量等試驗(yàn)參數(shù)。李寶方[5]對晃動試驗(yàn)的美國軍用規(guī)范進(jìn)行總結(jié)，并介紹了主要試驗(yàn)設(shè)備和方法。葉休乃[6]對民航適航性條例和美國軍用規(guī)范中有關(guān)油箱晃動和晃振試驗(yàn)標(biāo)準(zhǔn)進(jìn)行了總結(jié)?；蝿游锢碓囼?yàn)實(shí)施過程復(fù)雜，試驗(yàn)周期長，需要進(jìn)行大批量試驗(yàn)來涵蓋多種工況，代價(jià)較高，不能經(jīng)濟(jì)有效地反映出工況條件與響應(yīng)間的關(guān)系，不利于飛行器的快速高效適航審定。

隨著計(jì)算機(jī)和有限元的發(fā)展，數(shù)值模擬開始進(jìn)入油箱晃動領(lǐng)域。韋杰創(chuàng)[7]利用任意拉格朗日-歐拉法對航空彈性薄壁油箱開展流固耦合計(jì)算，驗(yàn)證了有限元計(jì)算的準(zhǔn)確性，并對不同的晃動幅值、晃動頻率和載液比例進(jìn)行了結(jié)果分析。李威鋒[8]建立了油箱流固耦合模型，分析了影響油箱晃動的因素，并結(jié)合動力學(xué)分析給出了油箱壽命預(yù)估。張雄等[9]提出了一種新型流固耦合不可壓物質(zhì)點(diǎn)法,并研究了運(yùn)動剛體容器中的液體晃動問題，解決了完全不可壓物質(zhì)點(diǎn)法無法處理不規(guī)則邊界和移動邊界的問題。與物理試驗(yàn)相比，數(shù)值模擬有效降低了經(jīng)濟(jì)成本，可以方便計(jì)算任意外部激勵下油箱的晃動響應(yīng)。但流固耦合計(jì)算復(fù)雜，單個(gè)算例需耗時(shí)若干小時(shí)，且基于離散工況的有限元無法給出外部激勵與結(jié)構(gòu)響應(yīng)規(guī)律，難以有效應(yīng)用于飛行器快速概念設(shè)計(jì)和適航審定。

響應(yīng)面代理模型由Box等[10]提出，核心是在試驗(yàn)設(shè)計(jì)基礎(chǔ)上建立響應(yīng)值與自變量之間的顯式函數(shù)來代理真實(shí)模型，對多變量非線性問題尤為適用[11-13]。文獻(xiàn)[14,15]對響應(yīng)面代理模型進(jìn)行改進(jìn)，不僅提高了模型精度，而且成功解決了多種復(fù)雜的工程問題，可潛在應(yīng)用于復(fù)合材料燃油箱晃動分析。響應(yīng)代理模型的精度和求解復(fù)雜程度依賴于樣本點(diǎn)的選取方法，如基于拉丁方設(shè)計(jì)開展等概率分層抽樣具有很好空間填充性和均衡性[16]的拉丁超立方采樣法。

本文使用拉丁超立方合理布置樣本點(diǎn)分布，基于復(fù)合材料燃油箱液體晃動有限元計(jì)算結(jié)果建立響應(yīng)面代理模型，給出晃動頻率、晃動幅值和載液比與油箱壁面峰值主應(yīng)變間的響應(yīng)面函數(shù)。通過數(shù)值模擬和物理試驗(yàn)進(jìn)行響應(yīng)面代理模型的校驗(yàn)?；诖砟Ｐ蛠眍A(yù)測復(fù)材油箱晃動最危工況，指導(dǎo)工程實(shí)際和飛機(jī)適航審定。

2 機(jī)翼油箱晃動物理試驗(yàn)和數(shù)值仿真

基于民航規(guī)章開展飛機(jī)機(jī)翼全尺寸復(fù)材燃油箱晃動試驗(yàn)。晃動試驗(yàn)的激勵方式根據(jù)《HB 6757-1993飛機(jī)燃油箱晃動和振動試驗(yàn)要求》施加，使油箱能夠繞晃動軸線(油箱晃動時(shí)所圍繞的固定軸線)作俯仰或橫滾方向晃動。試驗(yàn)件選取全尺寸復(fù)材飛機(jī)右機(jī)翼整體油箱，試驗(yàn)前對油箱進(jìn)行沖壓至內(nèi)部壓力比周圍大氣高10.5 kPa±0.5 kPa，并保壓半小時(shí)檢查其氣密性，試驗(yàn)后排空油箱內(nèi)液體，再次沖壓進(jìn)行氣密性檢查。試驗(yàn)夾具保證具有足夠的支持剛度來模擬機(jī)翼燃油箱與機(jī)身的實(shí)際連接。根據(jù)有限元計(jì)算結(jié)果，在油箱上下蒙皮主應(yīng)變較大且應(yīng)變梯度較小位置貼45°直角應(yīng)變花，測量壁面應(yīng)變分布和峰值主應(yīng)變。整個(gè)試驗(yàn)過程在室溫大氣環(huán)境下進(jìn)行。各晃動試驗(yàn)工況參數(shù)列入表1。由于全尺寸模型晃動試驗(yàn)經(jīng)濟(jì)成本和時(shí)間代價(jià)昂貴，本文僅進(jìn)行了三次物理試驗(yàn)用于數(shù)值仿真對比和響應(yīng)面模型驗(yàn)證。后續(xù)工作中將適當(dāng)增加試驗(yàn)次數(shù)增加有效實(shí)驗(yàn)數(shù)據(jù)。

表1 晃動試驗(yàn)工況

基于商用有限元軟件ANSYS Workbench開展機(jī)翼復(fù)材燃油箱晃動有限元流固耦合分析?；跈C(jī)型設(shè)計(jì)，將如圖 1所示的機(jī)翼結(jié)構(gòu)沿隔板分成六個(gè)油箱，1、2和3為主油箱，4，5和6為副油箱，分別考慮不同油箱的載液比。由于1號主油箱為研究晃動規(guī)律的重點(diǎn)，因此本文僅考慮1號主油箱的載液比影響。使用三維不可壓雷諾時(shí)均Navier-Stokes方程(RANS方程)和標(biāo)準(zhǔn)k-ε湍流方程描述液體流動，使用流體體積法VOF(Volume of fluid)模型進(jìn)行流體晃動自由液面追蹤。引入SIMPLE(Semi-Implicit Method for Pressure -Linked Equation)提高流體計(jì)算收斂速度。由于液體晃動不會導(dǎo)致油箱壁面的大變形，因此通過單向流固耦合將液面變化導(dǎo)致的動態(tài)變化壓力場傳給固體域結(jié)構(gòu)分析，最終得到機(jī)翼復(fù)材燃油箱壁面應(yīng)變分布。

（1）造型前準(zhǔn)備 流涂質(zhì)量、含水率、變形量等關(guān)鍵項(xiàng)經(jīng)過檢測合格后的黑模，可用于造型。首先，進(jìn)行墊料、填料的設(shè)計(jì)與裝配，除工藝設(shè)計(jì)床身外圍墊料外，在內(nèi)腔必須設(shè)計(jì)填料、襯料，以防止在流砂緊實(shí)過程中受壓力導(dǎo)致黑模變形，同時(shí)襯料要求布置合理，既能保證支撐又要能夠順利取出；另外，還需要在外側(cè)方門設(shè)置堵頭，目的是按工藝設(shè)定順序流砂，堵住的位置在下型反轉(zhuǎn)后再流砂。

圖1 機(jī)翼油箱結(jié)構(gòu)

3 機(jī)翼油箱晃動響應(yīng)面代理模型

(1)

式中βi,βi i和βi j分別為一階項(xiàng)、二階項(xiàng)和交叉項(xiàng)的待定系數(shù)，xi為第i個(gè)自變量，n為自變量個(gè)數(shù)，本文取n=3。

樣本點(diǎn)抽樣采用拉丁超立方等概率分層采樣法，最大化樣本點(diǎn)間距和最小化相關(guān)系數(shù)，使用盡可能小的樣本空間表征更多種隨機(jī)可能性。在晃動頻率F(15 次/min～45 次/min)，晃動幅值A(chǔ)(5°～20°)和油箱載液比L(0.3～1)范圍內(nèi)進(jìn)行采樣，步驟如下。

(1) 把自變量xi按區(qū)間等概率分為30個(gè)子區(qū)間，xi 0

(3) 多次重復(fù)步驟(2)，并利用樣本點(diǎn)之間最大距離最小化得到最優(yōu)樣本空間。

β=(XTX)-1XTy

(2)

且本文有

4 結(jié)果與討論

4.1 油箱上蒙皮代理模型與最危工況預(yù)測

針對1號主油箱，取有限元流固耦合分析中上蒙皮對應(yīng)試驗(yàn)應(yīng)變花15～19號貼片測試點(diǎn)中的最大主應(yīng)變峰值為響應(yīng)值，得到上蒙皮響應(yīng)面函數(shù)為

FL+0.3AL+0.03F2+0.1A2-83.4L2

(3)

式中F為晃動頻率，A為晃動幅值，L為載液比。

圖2 上蒙皮貼片點(diǎn)

另隨機(jī)選取6個(gè)工況，開展流固耦合數(shù)值計(jì)算得到主應(yīng)變響應(yīng)值。對比響應(yīng)面函數(shù)預(yù)測值與工況點(diǎn)計(jì)算值，驗(yàn)證函數(shù)模型準(zhǔn)確度。隨機(jī)工況點(diǎn)數(shù)據(jù)及其對比結(jié)果列入表2。由表2可知，響應(yīng)面函數(shù)預(yù)測結(jié)果與有限元計(jì)算結(jié)果最大誤差為9.6%。值得注意的是，最大誤差發(fā)生在響應(yīng)數(shù)值較小的工況，響應(yīng)數(shù)值大的工況對應(yīng)較小的誤差。出于安全性考慮，本文更多關(guān)注應(yīng)變響應(yīng)數(shù)值較大的危險(xiǎn)工況，因此可認(rèn)為響應(yīng)面函數(shù)可以有效預(yù)測響應(yīng)數(shù)值較大的危險(xiǎn)工況。

表2 復(fù)材燃油箱上蒙皮代理模型數(shù)值分析驗(yàn)證

開展全尺寸復(fù)材燃油箱晃動物理試驗(yàn)以驗(yàn)證響應(yīng)面函數(shù)正確性。表3為三種工況下函數(shù)預(yù)測值與試驗(yàn)值結(jié)果對比，1號主油箱的上蒙皮最大主應(yīng)變函數(shù)預(yù)測值與試驗(yàn)結(jié)果誤差在10%以內(nèi)，滿足工程要求。有限元數(shù)值計(jì)算和全尺寸物理試驗(yàn)表明，響應(yīng)面代理函數(shù)可準(zhǔn)確預(yù)測上蒙皮不同工況下的主應(yīng)變峰值響應(yīng)，可用于快速適航審定和最危工況預(yù)測。

表3 復(fù)材燃油箱上蒙皮代理模型物理試驗(yàn)驗(yàn)證

進(jìn)一步分析復(fù)材燃油箱上蒙皮響應(yīng)面函數(shù)隨晃動幅值、晃動頻率和載液比三個(gè)自變量的變化趨勢，結(jié)果如圖 3所示。當(dāng)頻率固定為20 次/min，幅值為15°時(shí)，最大主應(yīng)變峰值變化趨勢如 圖3(a)所示，存在最危險(xiǎn)載液比0.91使得主應(yīng)變峰值達(dá)到最大。當(dāng)固定載液比為0.7，幅值為15°時(shí)，如 圖3(b)所示，響應(yīng)面代理函數(shù)呈單調(diào)遞增特性，主應(yīng)變峰值隨晃動頻率增大而增大。當(dāng)固定載液比為0.7，頻率為20 次/min時(shí)，最大主應(yīng)變峰值基本隨晃動幅值呈遞增趨勢，如圖3(c)所示。

圖3 上蒙皮最大主應(yīng)變峰值變化趨勢

給定晃動幅值和頻率時(shí)，存在載液比使響應(yīng)面函數(shù)存在極大值，極大值對應(yīng)的工況即是工程實(shí)際中應(yīng)盡量避免和關(guān)注的危險(xiǎn)工況。幾種典型晃動工況下，最危險(xiǎn)載液比列入表4，可知隨著晃動頻率與晃動幅值的增加，最危險(xiǎn)載液比呈增加趨勢。

表4 上蒙皮各工況最危險(xiǎn)載液比

4.2 油箱下蒙皮代理模型與最危工況預(yù)測

有限元流固耦合分析中，取1號主油箱下蒙皮1～5，7，37，38貼片點(diǎn)的最大主應(yīng)變峰值為響應(yīng)，得到下蒙皮響應(yīng)面函數(shù)為

0.1FA+0.1FL-3.8AL-0.03F2+

0.1A2-206.7L2

(4)

圖4 下蒙皮貼片點(diǎn)

與上蒙皮類似，同樣選取6個(gè)驗(yàn)證樣本點(diǎn)開展流固耦合數(shù)值分析，驗(yàn)證響應(yīng)面代理函數(shù)模型的準(zhǔn)確度。驗(yàn)證樣本點(diǎn)數(shù)據(jù)及其對比結(jié)果列入表5，下蒙皮響應(yīng)面函數(shù)預(yù)測結(jié)果與有限元計(jì)算結(jié)果基本一致，最大誤差為6.2%，相差僅為5.17 με。

表5 復(fù)材燃油箱下蒙皮代理模型數(shù)值分析驗(yàn)證

同樣選取三組全尺寸機(jī)翼油箱晃動試驗(yàn)工況數(shù)據(jù)對響應(yīng)面代理函數(shù)模型的準(zhǔn)確性進(jìn)行驗(yàn)證。由表6可知，三種工況下函數(shù)預(yù)測值與物理試驗(yàn)結(jié)果誤差均在10%以內(nèi)。數(shù)值分析和試驗(yàn)驗(yàn)證表明，響應(yīng)面代理函數(shù)可準(zhǔn)確預(yù)測下蒙皮不同工況下的主應(yīng)變峰值響應(yīng)，用于快速適航審定。

表6 復(fù)材燃油箱下蒙皮代理模型物理試驗(yàn)驗(yàn)證

分析下蒙皮響應(yīng)面代理函數(shù)最大主應(yīng)變峰值隨晃動幅值、晃動頻率和載液比三個(gè)自變量的變化趨勢。給定晃動頻率和載液量，存在最危險(xiǎn)載液比使得響應(yīng)函數(shù)取極大值。當(dāng)頻率為20 次/min，幅值為15°時(shí)，最大主應(yīng)變峰值隨載液比的變化趨勢如 圖5(a)所示，最危險(xiǎn)載液比為0.93。各典型工況及其最危險(xiǎn)載液比列入表7。同時(shí)，需要指出的是，相同工況下，下蒙皮最大主應(yīng)變峰值比上蒙皮最大主應(yīng)變峰值大。這是因?yàn)槌杂梢好婊蝿拥挠绊懲?，下蒙皮還受到液體重力的影響，導(dǎo)致機(jī)翼下蒙皮變形增加，因此下蒙皮最大主應(yīng)變峰值會較大。

圖5 下蒙皮最大主應(yīng)變峰值變化趨勢

表7 下蒙皮各工況最危險(xiǎn)載液比

當(dāng)固定載液比和晃動幅值時(shí)，響應(yīng)面代理函數(shù)為晃動頻率的一元二次函數(shù)，在晃動頻率范圍內(nèi)呈單調(diào)遞增趨勢。以固定載液比0.7，幅值15°為例，最大主應(yīng)變峰值隨晃動頻率的變化趨勢如圖5(b)所示。當(dāng)固定載液比與晃動頻率時(shí)，下蒙皮最大主應(yīng)變峰值會隨著晃動幅值逐漸上升。以固定載液比0.7，頻率20 次/min為例，最大主應(yīng)變峰值隨晃動幅值的變化趨勢如圖5(c)所示。

對于復(fù)材燃油箱下蒙皮，在給定晃動頻率和晃動幅值時(shí)，存在一個(gè)最危險(xiǎn)載液比使得最大主應(yīng)變峰值最大?？蛰d或滿載時(shí)油箱不發(fā)生液體晃動，最危險(xiǎn)載液比時(shí)液體晃動明顯，液面與壁面之間存在相對運(yùn)動，液體對壁面有動態(tài)壓力作用，進(jìn)而造成主應(yīng)變呈非線性變化。固定載液比和幅值時(shí)，隨著晃動頻率的增大，周期晃動過程中的加速度最大值增加，油箱壁面對液體的作用力為液體晃動提供加速度，油箱內(nèi)液體對壁面的壓力增加，機(jī)翼變形增加，導(dǎo)致下蒙皮最大主應(yīng)變峰值會隨著晃動頻率的增大而增大。固定載液比和晃動頻率時(shí)，晃動幅值越大，液體慣性力對壁面作用越明顯，從而最大主應(yīng)變峰值隨著晃動幅值的增大而增大。同時(shí)，相同工況下，下蒙皮最大主應(yīng)變峰值變化趨勢與上蒙皮一致，且下蒙皮最大主應(yīng)變峰值比上蒙皮大，所以機(jī)翼油箱破壞多發(fā)生于下蒙皮位置。

5 結(jié) 論

本文通過拉丁超立方采樣方法和響應(yīng)面代理函數(shù)，以晃動頻率、晃動幅值和油箱載液比為自變量，以油箱上下蒙皮最大主應(yīng)變峰值為響應(yīng)，建立了上下蒙皮響應(yīng)面函數(shù)。流固耦合數(shù)值計(jì)算和全尺寸機(jī)翼油箱晃動試驗(yàn)表明了響應(yīng)面模型的有效性。同時(shí)，本文基于響應(yīng)面函數(shù)分析油箱晃動，發(fā)現(xiàn)隨著晃動幅值和頻率的增加，上下蒙皮最大主應(yīng)變峰值響應(yīng)呈非線性趨勢，存在最危險(xiǎn)載液比使得響應(yīng)函數(shù)取極大值。此外，下蒙皮最大主應(yīng)變峰值數(shù)值遠(yuǎn)大于相同工況下上蒙皮的主應(yīng)變峰值，油箱破壞多發(fā)生于下蒙皮位置。