基于Sobol’法的軌道車輛前端吸能結(jié)構(gòu)靈敏度分析

陳秉智，汪駒暢,2

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

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

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

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

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

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

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

( 1 )

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

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

( 2 )

( 3 )

fi,j(xi,xj)=

( 4 )

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

Sobol’法的總方差V是由所有輸入?yún)?shù)產(chǎn)生。

( 5 )

由單個(gè)參數(shù)影響產(chǎn)生的偏方差Vi為

( 6 )

由參數(shù)間交互效應(yīng)產(chǎn)生的偏方差Vi1,i2,…,is可以表達(dá)為

( 7 )

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

( 8 )

式中：s為設(shè)計(jì)參數(shù)個(gè)數(shù)。 然而有時(shí)關(guān)系函數(shù)很復(fù)雜，比如軌道列車吸能結(jié)構(gòu)的碰撞過程，很難找到準(zhǔn)確的關(guān)系函數(shù)表達(dá)，無法通過上述公式求得Sobol’指數(shù)，只能根據(jù)樣本點(diǎn)近似求解。

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

1.1 Kriging模型求解Sobol’指數(shù)

Kriging代理模型最早應(yīng)用于地質(zhì)學(xué)中，確定礦物的儲(chǔ)量，而后逐漸發(fā)展形成了完善的系統(tǒng)理論。Kriging模型主要借助附近的已知點(diǎn)而對(duì)某一點(diǎn)進(jìn)行模擬，通過一定范圍內(nèi)的已知信息，基于最小化誤差的方差，選擇加權(quán)線性組合來估計(jì)[12]。Kriging模型一般由回歸模型和隨機(jī)部分組成，具體形式為

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

( 9 )

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

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

(10)

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

相比于傳統(tǒng)的差值技術(shù)，Kriging代理模型根據(jù)已有樣本進(jìn)行構(gòu)造時(shí)充分考慮了空間上的相關(guān)特征，僅依靠待估點(diǎn)周圍的已知樣本，而非一律使用所有樣本進(jìn)行擬合，這使得Kriging模型兼具全局性和局部性[13]，可以動(dòng)態(tài)的反應(yīng)參數(shù)間關(guān)系。

1.2 數(shù)例分析

為驗(yàn)證Kriging代理模型法近似求解Sobol’指數(shù)的實(shí)際效果，下面對(duì)一個(gè)三元非線性函數(shù)進(jìn)行參數(shù)t1、t2和t3的靈敏度計(jì)算。

y=3et1+2sint2+t1t2t3

t1,t2,t3∈[1,2]

(11)

Kriging模型采用Latin超立方采樣，二次回歸函數(shù)，高斯相關(guān)方程。計(jì)算的Sobol’指數(shù)見表1。

表1 估算靈敏度值與真實(shí)值對(duì)比

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

2 吸能結(jié)構(gòu)全局靈敏度分析

2.1 車輛耐撞性仿真分析

車輛碰撞是一個(gè)瞬態(tài)的復(fù)雜物理過程，它包括以大位移、大轉(zhuǎn)動(dòng)和大應(yīng)變?yōu)樘卣鞯膸缀畏蔷€性，以材料彈塑性變形為典型特征的材料非線性和以接觸摩擦為特征的邊界非線性[14]。目前，在國內(nèi)外廣泛用于大變形碰撞問題的主要應(yīng)用軟件均基于動(dòng)量方程、質(zhì)量方程、能量守恒方程、邊界約束方程等[15]，且大多采用時(shí)間域的顯式中心差分法。本文采用碰撞仿真軟件PAM-CRASH進(jìn)行計(jì)算。該軟件應(yīng)用于汽車、鐵路車輛和航空航天等領(lǐng)域內(nèi)的碰撞研究工作中，其完善的碰撞模擬方案受到廣泛認(rèn)可。

2.2 耐撞性計(jì)算模型

結(jié)合某地鐵車輛頭車結(jié)構(gòu)，建立前端司機(jī)室有限元模型。以四節(jié)點(diǎn)殼單元為主，三節(jié)點(diǎn)三角形單元為輔劃分模型。其單元總數(shù)為23 466，節(jié)點(diǎn)總數(shù)為27 087。該司機(jī)室為對(duì)稱結(jié)構(gòu)，將對(duì)耐撞性能可能產(chǎn)生影響的結(jié)構(gòu)板厚作為變量，位置標(biāo)示見圖1。

結(jié)合設(shè)計(jì)經(jīng)驗(yàn)，設(shè)定各厚度參數(shù)下限為2 mm，上限為8 mm。參照BS EN 15227:2008標(biāo)準(zhǔn)[16]，同時(shí)考慮減少單元數(shù)量，提高計(jì)算效率，在司機(jī)室后端按照編組質(zhì)量進(jìn)行集中配重,整體以25 km/h的速度撞擊靜止剛性墻。

碰撞過程中，結(jié)構(gòu)的總能量變化和接觸力變化見圖2。

總吸能量和接觸力峰值是評(píng)價(jià)結(jié)構(gòu)耐撞吸能的重要指標(biāo)，因此分別以這兩個(gè)指標(biāo)作為目標(biāo)響應(yīng)，借助Matlab軟件中的Dace工具箱構(gòu)建Kriging模型，選擇二次回歸函數(shù)，高斯相關(guān)方程。

為檢驗(yàn)Kriging模型的準(zhǔn)確性，共選取了30個(gè)隨機(jī)點(diǎn)計(jì)算其誤差，模型擬合精度檢驗(yàn)見表2。其中以總吸能量為目標(biāo)函數(shù)時(shí)的最高誤差為2.07%，平均誤差為0.72%；以接觸力峰值為目標(biāo)函數(shù)時(shí)的最高誤差為1.92%，平均誤差為1.43%。

表2 Kriging模型擬合精度檢驗(yàn)

誤差在可接受范圍內(nèi)，模型精度滿足要求，可繼續(xù)進(jìn)行Sobol’靈敏度分析。

2.3 總吸能量為目標(biāo)響應(yīng)的靈敏度分析

主靈敏度表示的是單個(gè)參數(shù)變化對(duì)輸出結(jié)果的影響，而總靈敏度還反映了與其他參數(shù)交互作用共同產(chǎn)生的影響。以總吸能能量為目標(biāo)時(shí)，各輸入?yún)?shù)的主靈敏度和總靈敏度見圖3。

由圖3可見，設(shè)計(jì)參數(shù)x1，x2，x3的靈敏度值較大，對(duì)結(jié)構(gòu)總吸能量有較大影響。同時(shí)注意到參數(shù)x2,x3,x6的主靈敏度與總靈敏度相差較大，說明該參數(shù)與其他參數(shù)交互效應(yīng)對(duì)輸出參數(shù)產(chǎn)生了較大影響，因此針對(duì)各參數(shù)間的交互效應(yīng)進(jìn)一步分析，得到圖4所示二階交互效應(yīng)熱圖。

由圖4可見，參數(shù)x2與x3、x2與x6的交叉處色塊顏色較深，靈敏度較高，因此在結(jié)構(gòu)優(yōu)化中，在對(duì)總靈敏度較大的參數(shù)重點(diǎn)優(yōu)化時(shí)，應(yīng)同時(shí)注重這兩對(duì)參數(shù)間的合理搭配，以達(dá)到更優(yōu)的結(jié)果。

2.4 接觸力峰值為目標(biāo)響應(yīng)的靈敏度分析

同理，將接觸力峰值設(shè)為目標(biāo)響應(yīng)，得到的一階主靈敏度和總靈敏度值見圖5。

設(shè)計(jì)參數(shù)x1，x2的靈敏度值較大，對(duì)接觸力峰值有較大影響。同時(shí)注意到各參數(shù)的主靈敏度和總靈敏度相差較小，因此可推斷交互效應(yīng)所產(chǎn)生的影響會(huì)很小。各參數(shù)二階交互效應(yīng)熱圖也驗(yàn)證了這一預(yù)測(cè)，見圖6。

由圖6可看出，參數(shù)間二階交互效應(yīng)均較小，因此在優(yōu)化中可合理忽略參數(shù)間的搭配組合，著重對(duì)主靈敏度較大的x1、x2參數(shù)進(jìn)行優(yōu)化。

3 優(yōu)化結(jié)果對(duì)比

根據(jù)靈敏度分析結(jié)果，將原模型的7個(gè)設(shè)計(jì)參數(shù)簡化為僅包含靈敏度較高的4個(gè)參數(shù)的簡化模型，利用ISIGHT軟件，分別對(duì)原模型和簡化后的模型進(jìn)行尺寸優(yōu)化，優(yōu)化結(jié)果見表3。

表3 原模型與簡化模型優(yōu)化結(jié)果對(duì)比

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

4 結(jié)論

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

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

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

(3) 通過對(duì)某地鐵前端吸能結(jié)構(gòu)的全局靈敏度分析，可找出較為敏感的參數(shù)，對(duì)結(jié)構(gòu)優(yōu)化提供有力參考。