基于自適應隨機優化的連續陣風關鍵載荷預測

肖宇

中國商飛上海飛機設計研究院 結構強度研究所，上海 201210

對于民用客機來說，連續陣風載荷是設計載荷的關鍵組成部分，對民機設計與取證都具有極其重要的意義[1]。適航條款針對陣風載荷，提出在設計范圍內，須滿足相應的強度要求，以確保整個設計包線內飛機結構的完整性。為此，設計人員需要在包線邊界上以及內部選擇足夠多的載荷工況進行評估，確認飛機結構所有部件對應的關鍵受載工況都可被準確識別，用于下游專業進行強度校核。

因此，一方面，對于一輪載荷分析，理論上需要計算的工況會達到100萬的量級，其中對于常規構型的民用飛機，設計員可以通過現有經驗篩選，進而縮減至10萬的量級[2]。然而，對于目前的載荷分析模型，其結構上相臨2站位的關鍵工況極有可能一致，實際上供下游專業強度校核的嚴重工況較少，因此大量的工況無計算的意義。

另一方面，對于非常規布局的民用飛機[3]，如采用大展弦比機翼等設計[4]，可能涉及到各種非線性效應[5]，導致其飛行包線以及關鍵載荷工況將會與現有型號截然不同，因此將無法從現有經驗來實現載荷工況的縮減，進而導致設計周期及費用的大幅增加，而考慮到其中大量載荷工況的不必要性，有必要開展關鍵載荷的主動識別技術研究，通過對這些關鍵工況的搜尋，避免大部分不構成包線工況的計算，從而顯著提高陣風載荷計算的效率。

近些年，國外在此領域進行了較多的研究[6]，如歐盟在FFAST(Future Fast Aeroelastic Simulation Technologies)項目中第1項專題便是關鍵載荷的識別，在該項目資助下，Khodaparast等[2]采用Kriging建立代理模型，并基于基因優化算法實現了對1-cos型離散陣風的關鍵載荷識別，隨后Khodaparast和Cooper[7]又應用于模型更改后的快速載荷計算中。Cavagna等[8]則采用神經網絡構建了代理模型，并與陣風氣彈求解器耦合，開發了軟件Fast-GLP(Fast tool for worst case Gust Load Prediction)，通過數值模擬，表明了關鍵載荷預測的準確性。Castellani等[9]則采用參數化模型降階的方法進行了動載荷的快速預測，效率同樣獲得了較為顯著的提升。

但上述研究中大都針對連續的變量進行識別，同時商載、油載等重量分布情況也未予以考慮，而國內民機陣風載荷計算中[10]，在各特征重量、各飛行高度下，必須真實考慮燃油消耗邏輯次序進行重量分布計算，因此在工程中只能事先對關鍵特征參數進行離散，從而形成大量的工況文件。對于這些真實考慮燃油消耗邏輯的陣風載荷工況，對其進行關鍵載荷的準確識別，具有較大的挑戰。

鑒于此，本文針對國內民機型號工作，創新性地提出了一種采用自適應隨機優化算法的陣風關鍵工況分析流程。首先，通過DoE(Design of Experiment)[11]方法，計算初始工況并建立初始代理模型，同時采用自適應隨機優化算法，在極值載荷工況附近尋找新的可能極值工況，進行計算并更新代理模型，在迭代收斂后即可確定所關注載荷的全局最大值。然后，以某型民機側向連續陣風模型作為算例，結果顯示，在現有陣風模型計算過程中，本文方法由于避免大量不必要工況的計算，可以有效地提高連續陣風關鍵載荷的預測效率。

1 計算方法

1.1 初始值設定

采用DoE方法生成計算模型的初始樣本點，以盡量少的計算次數，獲取到盡量多的模型特性。DoE方法可以分為3類：一是因子類設計(factorial designs)，二是響應面類設計(response-surface designs)，三是隨機類設計(randomized designs)。圖1依次給出了三維空間內，3種模型產生的樣本點對比，其中黑色點即為樣本點。從圖1可以看出，二水平全因子(2LFF)采樣與Box-Behnken采樣的樣本點都在邊界上，而拉丁超立方采樣(Latin hypercube sampling)則大部分樣本點存在于內部、且具有隨機性，考慮到目前工程中陣風載荷計算工況出現在邊界的概率較大，本文基于前2種方法開展對比分析研究。

圖1 3種典型DoE方法對比Fig.1 Comparison of three typical DoE approaches

1.2 代理模型構建

考慮到單個工況計算耗時較多，有必要建立一個可以快速評估的代理模型來替代原有的陣風計算模型，進而降低計算量，便于后續優化工作的開展。代理模型典型的有多項式響應面[12]、Kriging[13]、徑向基函數[14]、人工神經網絡[15]、支持向量機[16]等，綜合考慮計算效率、精度以及模型透明性，本文采用多元自適應回歸樣條(MARS)模型[17]進行構建，該方法相對Kriging法，其計算效率較高，同時無需像徑向基函數法一樣事先假設復雜非線性基函數形狀，因此具有一定的優勢。

MARS方法由Friedman[17]首次提出，其采用簡單的鉸鏈函數(hinge function)及其積作為基函數，基函數形式及個數可以按精度要求自動確定，因此具有較強的自適應性，研究表明，該方法在樣本數50～1 000，以及變量維數3～20具有較明顯的優勢[17]。

MARS所采用的鉸鏈函數為

(1)

(2)

式中：下標+表示函數取正部分；t為節點坐標值。兩者可構成一個映射對，所有的映射對可構成集合C，即

C={(xj-t)+，(t-xj)+}

t∈{x1j,x2j，…，xNj},j∈{1，2，…,p}

(3)

式中：N為樣本空間維數；p為樣本點個數。MARS構建過程中采用的基函數皆來自于集合C中的任意一個或多個鉸鏈函數的乘積。

整個建模分為前向和后向2個過程，最終的代理模型為

(4)

式中：hm為基函數；βm為各基函數的權重；M為基函數總數。具體過程參見文獻[17]。

1.3 基于自適應隨機算法的優化模型

本文采用自適應隨機優化算法開展關鍵工況的挑選[18]。假設可能的未知關鍵工況位于當前關鍵工況的附近，且滿足截斷正態分布規律，圖2給出了參數x在標準差σ為0.2、0.5、0.7下的概率密度曲線對比，其中x取值范圍為[1,7]，最優點(即當前關鍵工況)為4。從圖2可以看出，概率密度在最優點達到最大值，隨著與最優點的距離增大概率逐漸趨近于0，同時在取值范圍外概率則嚴格等于0，確保只在參數的有效范圍內取值。另外，隨著標準差σ減小，搜索范圍逐漸減小。

圖2 截斷正態分布概率密度曲線Fig.2 Function of probability density for truncated normal distribution

獲取得到當前最優工況后，在其附近隨機生成的若干個樣本點，計算樣本點與最優點的距離值，同時采用所建立的代理模型，評估得到所有工況的目標值，在此基礎上，將歸一化的目標值與距離值以一定的權重系數求和，從而得到下一次計算的工況。

考慮到本文所有計算針對的是離散的工況文件，上述樣本點計算時都需采取一定的策略映射到真實的工況集，本文假定樣本點定義域為整數域，從而可快捷地實現映射的目的。

另外，為提高魯棒性及收斂性，在隨機數生成的過程中，采用自適應技術確定標準差σ，給定σmax和σmin，初始標準差設定為σmax，在進行nf次連續的失敗關鍵工況識別后將減小標準差，而在進行ns次連續的成功關鍵工況識別后增大標準差以提高收斂速度，最終在標準差值達到σmin后認為計算得到收斂，從而達到實現自適應調節計算參數的目的。

1.4 陣風關鍵載荷計算

本文陣風關鍵載荷的計算流程如圖3所示。首先，依所選用的DoE方法確定初始工況，并通過建立的陣風氣彈模型計算得到各載荷值；然后，在此基礎上通過MARS方法建立初始代理模型，在極值載荷工況附近通過隨機優化方法得到下一個可能的極值工況，通過記錄比較連續失敗次數nf和連續成功次數ns，實現自適應的調整標準差，重復進行迭代，直到計算流程收斂。

2 模型驗證

2.1 模型介紹

針對某高平尾、下單翼、尾吊型飛機的側向連續陣風關鍵載荷開展建模研究，其外形如圖4所示，主要計算參數如表1所示。

采用梁架單元、偶極子網格法建立陣風載荷氣彈耦合分析模型，如圖5所示。系統、燃油、商載、發動機等質量以等效方式掛載至相應節點，另外，在機翼和平尾上布置若干單元用于氣動與結構的插值。

圖3 陣風關鍵載荷的計算流程圖Fig.3 Flow chart of computation of critical gust load

圖4 模型外形示意圖Fig.4 Schematic diagram of model’s shape表1 主要參數Table 1 Main parameters

參數數值全機長/m33.5全機高/m8.4機翼面積/m279.9展長/m26.6展弦比8.9

圖5 陣風載荷計算模型示意圖Fig.5 Schematic diagram of gust load computation model

在模態空間內進行陣風氣彈耦合計算[10]，控制方程為

(5)

式中：u為廣義位移解；ω為圓頻率；ρ和v分別為空氣密度和速度；g為結構阻尼系數；M、C和K分別為廣義質量矩陣、阻尼矩陣和剛度矩陣；Q為廣義氣動力系數矩陣；P為陣風引起的廣義外部激勵，其主要由陣風模型決定。本文計算采用適航條款規定的von Karman譜來定義陣風形狀：

(6)

式中：w為陣風速度均方根值；L為陣風尺度；Ω為折算頻率[10]。

計算模型中將飛行高度離散為0～12 000 m的8個特征高度(H0，H2,…,H7)，速度離散為對應最大陣風強度的設計速度(V0)、設計巡航速度(V1)和設計俯沖速度(V2)，飛行特征重量有最大設計起飛重量(C0)、最大設計著陸重量(C1)、最大設計零油重量(C2)、45 min余油重量(C3)4種，其各高度的重量分布已經充分考慮燃油的消耗邏輯，另外還有前重心(G0)、后重心(G1)2種不同的重心分布，構成192個陣風計算工況，如表2所示，該工況集的載荷包絡值作為本文關鍵載荷預測對比的基準值。

表2 計算參數Table 2 Parameters for computation

2.2 方法驗證

針對3個感興趣的站位載荷進行對比分析，分別為機頭、機身(機翼連接處)和垂尾根部，如圖6 中黑色六角形所示，主要關心的載荷為以均方根值表示的增量彎矩、剪力以及扭矩。計算中最大標準差值σmax設置為0.02，最小標準差值σmin設置為0.005，nf和ns分別為5和3。

圖7給出了機頭無量綱彎矩收斂曲線，其中以機頭處的彎矩基準值作為無量綱參數，從圖中可以看出，雖然機頭彎矩量值較小，但是本文建立的方法依舊可以找出關鍵工況，另外，采用Box-Behnken迭代次數比2LFF采樣要多，表明對現行陣風載荷計算，采用2LFF選取初始樣本點具有較大的優勢。

圖8給出了機身無量綱彎矩收斂曲線，對于機身機翼連接處，其彎矩具有較大的量級，2種方法都可以找到該站位的關鍵工況，其中Box-Behnken相對迭代次數較多，對于2LFF來說，其樣本點恰巧包括了關鍵載荷工況，但本文優化模型計算依舊進行了若干次迭代，直至標準差達到收斂閾值σmin。圖9給出了垂尾根部無量綱扭矩的收斂曲線，從圖中可以看出，與機身類似，Box-Behnken收斂需要較多迭代，2LFF相對優勢較為明顯。

圖6 站位示意圖Fig.6 Schematic diagram of stations

圖7 機頭無量綱彎矩收斂曲線Fig.7 Convergence curve of dimensionless bending moment of nose

圖8 機身無量綱彎矩收斂曲線Fig.8 Convergence curve of dimensionless bending moment of mid-airframe

圖9 垂尾根部無量綱扭矩收斂曲線Fig.9 Convergence curve of dimensionless root torque of vertical tail

3 參數評估

3.1 代理模型評估

下面對收斂后的代理模型MARS-2LFF進行評估。首先給出機身彎矩結果，如圖10所示，其中黑色網格為實際模型計算得到的基準值，綠色直線為代理模型計算得到的值，藍色點為代理模型構建所取的點，C1G0表示在最大設計著陸重量、前重心的組合下，載荷隨速度、高度的變化。從圖10可以看出，除個別邊界點外，本文采用MARS方法構建的代理模型一定程度上可以模擬真實物理模型的特性，同時本文方法可以準確地尋找確定最大載荷值，而對于邊界及外部點的預測，現有代理模型都不可避免地存在問題，本文方法同樣不例外。

圖11給出了垂尾根部扭矩代理模型與基準值對比結果，從圖中可以看出，代理模型與基準值同樣吻合得較好，同時結合圖10(c)與(d)可以看出，連續陣風載荷的量值受特征重量參數的影響較弱，其很大程度上取決于飛行高度與速度。

圖10 機身彎矩代理模型評估Fig.10 Evaluation of bending moment of mid-airframe

圖11 垂尾根部扭矩代理模型評估Fig.11 Evaluation of root torque of vertical tail

圖12和圖13給出MARS-2LFF預測值與基準值的相關性對比圖，分別對應192種工況的機身及垂尾根部彎矩和剪力。相比之下，垂尾的預測結果較好，這可能是由于垂尾根部受載情況相對簡單，而機身與機翼連接處環境相對復雜。另外，由于MARS模型構建過程中為提高泛化能力進行了后向過程，因此對于關鍵工況及附近載荷預測同樣存在誤差，但是偏差總體較小。

圖12 MARS-2LFF代理模型相關性(機身)Fig.12 Correlation plot of MARS-2LFF surrogate model(mid-airframe)

圖13 MARS-2LFF代理模型相關性(垂尾)Fig.13 Correlation plot of MARS-2LFF surrogate model(vertical tail)

3.2 收斂性評估

針對自適應隨機優化算法，并結合本文建立的側向連續陣風氣彈模型，進行參數影響的研究。模型涉及到參數主要有4個，分別是影響搜索半徑的標準差最大值σmax、最小值σmin以及2個自適應參數ns、nf。

表3～表6給出了機身、垂尾根部扭矩隨參數的變化規律。由表3～表6可知，一方面，本文模型可以較準確地預測關鍵工況及其載荷，即使設置標準差最大值、最小值相同，即僅進行一次圖3 中的自適應循環，誤差最大也僅為0.8%，且分布在對角線附近，這表明模型的準確性只取決于σmax與σmin的相對倍數，其絕對值對結果影響較小。另一方面，自適應參數ns和nf的減小一定程度上可以降低總迭代次數，而對精度影響較小。同時，針對本文模型，最多56次的工況迭代即可得到準確的關鍵載荷，相對工程實際的192次計算效率至少提升70%。

表3 標準差對機身扭矩預測影響(ns=1，nf=3)Table 3 Influence of standard deviation on prediction of mid-airframe torque (ns=1，nf=3)

表4 標準差對機身扭矩預測影響(ns=3，nf=5)Table 4 Influence of standard deviation on prediction of mid-airframe torque (ns=3，nf=5)

表5 標準差對垂尾根部扭矩預測影響(ns=1，nf=3)Table 5 Influence of standard deviation on prediction of vertical tail root torque (ns=1，nf=3)

表6 標準差對垂尾根部扭矩預測影響(ns=3，nf=5)Table 6 Influence of standard deviation on prediction of vertical tail root torque (ns=3，nf=5)

表7給出了采用參數σmin=σmax=1、ns=1、nf=3預測得到的結果，雖然該參數迭代次數較少，但仍可成功實現對高度(H)、速度(V)和重心位置(G)的準確識別，僅對特征重量項(C)的預測錯誤，即便如此，扭矩誤差值也在1%以內。

表7 關鍵工況及載荷對比Table 7 Comparisons of critical cases and loads

4 結 論

本文采用2LFF采樣技術結合MARS方法，建立了一個可進行快速預測的代理模型，并基于自適應隨機優化算法，發展了一套陣風關鍵工況及載荷的主動識別預測方法，并以某型民機側向連續陣風計算模型為例，進行了驗證，同時對所建立算法各參數的影響進行了深入的研究，結果表明：

1) 2LFF采樣相對于其他DoE方法可更有效地用于目前陣風載荷代理預測模型的構建。

2) 采用MARS構建的代理模型，結合自適應隨機優化算法可有效地識別陣風關鍵工況及載荷，可有效避免大量不必要工況的計算從而顯著提高效率。

3) 自適應參數對預測精度影響較小，適當減小可以有效地降低總迭代次數。