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

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

1 引 言

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

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

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

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

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

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

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

表1 晃動試驗工況

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

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

圖1 機翼油箱結構

3 機翼油箱晃動響應面代理模型

(1)

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

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

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

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

β=(XTX)-1XTy

(2)

且本文有

4 結果與討論

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

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

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

(3)

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

圖2 上蒙皮貼片點

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

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

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

表3 復材燃油箱上蒙皮代理模型物理試驗驗證

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

圖3 上蒙皮最大主應變峰值變化趨勢

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

表4 上蒙皮各工況最危險載液比

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

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

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

0.1A2-206.7L2

(4)

圖4 下蒙皮貼片點

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

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

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

表6 復材燃油箱下蒙皮代理模型物理試驗驗證

分析下蒙皮響應面代理函數(shù)最大主應變峰值隨晃動幅值、晃動頻率和載液比三個自變量的變化趨勢。給定晃動頻率和載液量，存在最危險載液比使得響應函數(shù)取極大值。當頻率為20 次/min，幅值為15°時，最大主應變峰值隨載液比的變化趨勢如 圖5(a)所示，最危險載液比為0.93。各典型工況及其最危險載液比列入表7。同時，需要指出的是，相同工況下，下蒙皮最大主應變峰值比上蒙皮最大主應變峰值大。這是因為除自由液面晃動的影響外，下蒙皮還受到液體重力的影響，導致機翼下蒙皮變形增加，因此下蒙皮最大主應變峰值會較大。

圖5 下蒙皮最大主應變峰值變化趨勢

表7 下蒙皮各工況最危險載液比

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

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

5 結 論

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