基于Sobol’法的軌道車輛前端吸能結構靈敏度分析

陳秉智，汪駒暢,2

(1. 大連交通大學 機車車輛工程學院， 遼寧 大連 116028；2. 中車青島四方機車車輛股份有限公司 轉向架分廠， 山東 青島 266111)

軌道交通在我國交通運輸中一直扮演著重要的角色。但由于其運行速度高、質量大、盲點多的特點，如果發生交通事故所產生的后果必然十分嚴重，所以提高軌道車輛車體結構耐撞性尤為重要。而影響車體結構耐撞性的設計參數眾多，因此很有必要通過靈敏度分析辨別其中的關鍵參數，簡化優化模型[1]。吸能結構的靈敏度可定義為吸能結構各參數發生變化時，對耐撞性能的影響程度。

靈敏度分析方法大致可以分為局部靈敏度分析法和全局靈敏度分析法。其中局部靈敏度分析法依靠線性模型的基礎發展起來，因此，當研究非線性問題、或者輸入變量的不確定性在不同的數量級時，局部靈敏度分析很難提供有效的結果，并且無法考慮多個變量共同作用時各變量之間的交互作用對結果的影響[2]。而全局靈敏度分析法允許各參數同時變化，且能分析多個參數交互效應的靈敏度。全局靈敏度分析包括回歸法、方差法和Sobol’法等。其中Sobol’法基于方差分析全局靈敏度,通過計算各因素的輸入對輸出方差的影響,評估單個以及多個因素交互效應的靈敏度。Sobol’法形式簡單,計算方便,特別是當模型無確定的函數表達時,可以采用蒙特卡洛方法或Kriging代理模型近似計算[3]。于德介等[4]以線性阻尼、立方剛度的非線性被動隔振體為對象, 用Sobol’法計算得出了軟、硬彈簧和線性阻尼對傳遞率的影響規律。張小麗等[5]使用Sobol’法分析了新安江三水源模型的參數靈敏度, 結果表明Sobol’法得到的敏感度適合多目標間的對比。陳靜等[6]用通過構造雙層底板架結構強度和穩定性的Kriging代理模型，對雙層底板架結構各設計變量對板架強度和穩定性響應指標進行了全局靈敏度分析。邵永生等[7]將Sobol’法用于Sperling指標的全局靈敏度分析,得到各個參數對車輛平穩性指標的一階靈敏度及總靈敏度。聶祚興等[8]將Sobol’全局靈敏度分析法引入到汽車噪聲傳遞函數的靈敏度分析中，有效避免了局部靈敏度分析方法的缺陷。陳剛等人[9]采用Sobol’法對沖擊載荷作用下的二維管路系統材料的彈性模量、密度、管壁厚、內徑,支承剛度及位置及沖擊作用間隔時間等參數進行了全局靈敏度分析,甄別出影響管路系統沖擊動響應的關鍵參數。于德介等人[10]將基于方差的Sobol’全局靈敏度方法引入結構動力分析中，計算了傳遞率對結構參數的一階及總靈敏度 ,其結果為系統的優化設計提供依據。

但在軌道車輛領域，尤其耐撞性研究中，該方法應用較少。本文針對某地鐵車輛司機室前端吸能結構的眾多設計參數，采用Sobol’法進行靈敏度分析。同時考慮到碰撞過程參數間無準確關系函數以及計算新樣本較為耗時的特點，無法使用公式法和蒙特卡羅積分法進行計算。因此，本文通過構建具有較高擬合精度的Kriging代理模型求解Sobol’指數，分析各參數的一階主靈敏度值、總靈敏度值及二階交互效應靈敏度。結果有助于理解不同的輸出響應對各參數及參數間相互作用的敏感性，識別對耐撞性能產生主要影響的設計參數，為吸能結構優化設計提供有力參考。將復雜模型合理簡化，提高優化效率。

1 Sobol’全局靈敏度分析法原理

Sobol’全局靈敏度分析法認為目標的總方差是由單個參數產生的方差和參數間的相互作用產生的方差疊加而成的[11]。所以其核心思想就是將函數f(X)分解為如下形式

…+f1,2,…,n(x1,x2,…,xn)

( 1 )

式中：f0為常數項；xi、xj為不同的變量；n為變量總數；X為所有變量組成的向量；fi(xi)為函數中只含有單變量的部分；fi,j(xi,xj)為函數中雙變量部分，以此類推。

為保證分解形式的惟一性和可求性，規定右邊每一項對其所包含的任一變量的積分為零，因此方程右邊各項前三項可以表達為式( 2 )～式( 4 )。

( 2 )

( 3 )

fi,j(xi,xj)=

( 4 )

式中：i、j、k為設計變量編號。式( 1 )中右邊其他子項均可由此類推得到。

Sobol’法的總方差V是由所有輸入參數產生。

( 5 )

由單個參數影響產生的偏方差Vi為

( 6 )

由參數間交互效應產生的偏方差Vi1,i2,…,is可以表達為

( 7 )

因此，將方差的比值定義為參數的Sobol’靈敏度指數，該數值的大小表示對輸出結果的影響程度。Sobol’靈敏度指數Si1,i2,…,is為

( 8 )

式中：s為設計參數個數。 然而有時關系函數很復雜，比如軌道列車吸能結構的碰撞過程，很難找到準確的關系函數表達，無法通過上述公式求得Sobol’指數，只能根據樣本點近似求解。

較為常用的基于樣本點的求解方法是蒙特卡洛積分法和Kriging代理模型法。其中蒙特卡羅積分法通過構建2個N·s的矩陣A、B(N為規定的樣本數)及分別交換各列得到2s個新矩陣作為輸入矩陣，計算總方差與偏方差，進而求得Sobol’指數。這種方法計算較為簡便，但其精度受采樣方式及樣本數量N直接影響，且需要重新計算的樣本點個數為N·2(s+1)，在計算本就十分耗時的耐撞性分析中不太實用。因此，本文采用Kriging代理模型法求解Sobol’指數。

1.1 Kriging模型求解Sobol’指數

Kriging代理模型最早應用于地質學中，確定礦物的儲量，而后逐漸發展形成了完善的系統理論。Kriging模型主要借助附近的已知點而對某一點進行模擬，通過一定范圍內的已知信息，基于最小化誤差的方差，選擇加權線性組合來估計[12]。Kriging模型一般由回歸模型和隨機部分組成，具體形式為

y(X)=gT(X)β+z(X)

( 9 )

式中：β為回歸參數；gT(X)為回歸模型，通常是多項式函數，為模型提供全局近似；z(X)為隨機部分，服從正態分布N(0,σ2)，為模型提供局部近似。z(X)的協方差矩陣Cov|z(xi,xj)|可以表達為

Cov|z(xi,xj)|=σ2|R(xi,xj)|

(10)

式中：R(xi,xj)為任意2個樣本點空間相關性函數，對模擬的精度起決定性作用。一般可選擇線性方程、指數方程、高斯方程、三次樣條等。其中計算效果最好、被廣泛應用的是高斯方程。

相比于傳統的差值技術，Kriging代理模型根據已有樣本進行構造時充分考慮了空間上的相關特征，僅依靠待估點周圍的已知樣本，而非一律使用所有樣本進行擬合，這使得Kriging模型兼具全局性和局部性[13]，可以動態的反應參數間關系。

1.2 數例分析

為驗證Kriging代理模型法近似求解Sobol’指數的實際效果，下面對一個三元非線性函數進行參數t1、t2和t3的靈敏度計算。

y=3et1+2sint2+t1t2t3

t1,t2,t3∈[1,2]

(11)

Kriging模型采用Latin超立方采樣，二次回歸函數，高斯相關方程。計算的Sobol’指數見表1。

表1 估算靈敏度值與真實值對比

由表1可看出，Kriging代理模型法計算的Sobol’指數與公式法求得的真實值十分接近，達到了較高的精度。故該方法可用于吸能結構全局靈敏度分析中近似求解Sobol’指數。

2 吸能結構全局靈敏度分析

2.1 車輛耐撞性仿真分析

車輛碰撞是一個瞬態的復雜物理過程，它包括以大位移、大轉動和大應變為特征的幾何非線性，以材料彈塑性變形為典型特征的材料非線性和以接觸摩擦為特征的邊界非線性[14]。目前，在國內外廣泛用于大變形碰撞問題的主要應用軟件均基于動量方程、質量方程、能量守恒方程、邊界約束方程等[15]，且大多采用時間域的顯式中心差分法。本文采用碰撞仿真軟件PAM-CRASH進行計算。該軟件應用于汽車、鐵路車輛和航空航天等領域內的碰撞研究工作中，其完善的碰撞模擬方案受到廣泛認可。

2.2 耐撞性計算模型

結合某地鐵車輛頭車結構，建立前端司機室有限元模型。以四節點殼單元為主，三節點三角形單元為輔劃分模型。其單元總數為23 466，節點總數為27 087。該司機室為對稱結構，將對耐撞性能可能產生影響的結構板厚作為變量，位置標示見圖1。

結合設計經驗，設定各厚度參數下限為2 mm，上限為8 mm。參照BS EN 15227:2008標準[16]，同時考慮減少單元數量，提高計算效率，在司機室后端按照編組質量進行集中配重,整體以25 km/h的速度撞擊靜止剛性墻。

碰撞過程中，結構的總能量變化和接觸力變化見圖2。

總吸能量和接觸力峰值是評價結構耐撞吸能的重要指標，因此分別以這兩個指標作為目標響應，借助Matlab軟件中的Dace工具箱構建Kriging模型，選擇二次回歸函數，高斯相關方程。

為檢驗Kriging模型的準確性，共選取了30個隨機點計算其誤差，模型擬合精度檢驗見表2。其中以總吸能量為目標函數時的最高誤差為2.07%，平均誤差為0.72%；以接觸力峰值為目標函數時的最高誤差為1.92%，平均誤差為1.43%。

表2 Kriging模型擬合精度檢驗

誤差在可接受范圍內，模型精度滿足要求，可繼續進行Sobol’靈敏度分析。

2.3 總吸能量為目標響應的靈敏度分析

主靈敏度表示的是單個參數變化對輸出結果的影響，而總靈敏度還反映了與其他參數交互作用共同產生的影響。以總吸能能量為目標時，各輸入參數的主靈敏度和總靈敏度見圖3。

由圖3可見，設計參數x1，x2，x3的靈敏度值較大，對結構總吸能量有較大影響。同時注意到參數x2,x3,x6的主靈敏度與總靈敏度相差較大，說明該參數與其他參數交互效應對輸出參數產生了較大影響，因此針對各參數間的交互效應進一步分析，得到圖4所示二階交互效應熱圖。

由圖4可見，參數x2與x3、x2與x6的交叉處色塊顏色較深，靈敏度較高，因此在結構優化中，在對總靈敏度較大的參數重點優化時，應同時注重這兩對參數間的合理搭配，以達到更優的結果。

2.4 接觸力峰值為目標響應的靈敏度分析

同理，將接觸力峰值設為目標響應，得到的一階主靈敏度和總靈敏度值見圖5。

設計參數x1，x2的靈敏度值較大，對接觸力峰值有較大影響。同時注意到各參數的主靈敏度和總靈敏度相差較小，因此可推斷交互效應所產生的影響會很小。各參數二階交互效應熱圖也驗證了這一預測，見圖6。

由圖6可看出，參數間二階交互效應均較小，因此在優化中可合理忽略參數間的搭配組合，著重對主靈敏度較大的x1、x2參數進行優化。

3 優化結果對比

根據靈敏度分析結果，將原模型的7個設計參數簡化為僅包含靈敏度較高的4個參數的簡化模型，利用ISIGHT軟件，分別對原模型和簡化后的模型進行尺寸優化，優化結果見表3。

表3 原模型與簡化模型優化結果對比

靈敏度分析辨別了設計參數中較為敏感的參數，簡化模型進行優化，迭代步數減少46%，極大降低了優化工作量。得到的優化結果與考慮所有設計參數的優化結果相比，總吸能量僅減少了1.5%，而接觸力峰值則相應降低了70 kN，差異極小且各有優劣。此結果反映了Sobol’靈敏度方法在軌道車輛吸能結構優化設計過程中的可行性和實用性，不僅能較大程度上減少實驗成本和計算耗時，且能得到極好的優化結果，具有現實指導意義。

4 結論

本文使用Sobol’方法，通過構建Kriging代理模型計算了軌道車輛前端吸能結構的全局靈敏度，并將基于靈敏度分析簡化模型后的優化結果與原模型優化結果進行了對比，研究結果表明：

(1) Sobol’全局靈敏度法能處理大型復雜非線性問題，并可分析各參數間交互效應對結果的影響。

(2) 構建Kriging代理模型計算Sobol’靈敏度僅基于已有樣本，同時具有較高精度。

(3) 通過對某地鐵前端吸能結構的全局靈敏度分析，可找出較為敏感的參數，對結構優化提供有力參考。