基于Sobol法的復合材料加筋板結構參數靈敏度

殷俊清, 趙誠誠, 陳永當, 程云飛, 顧金芋

(西安工程大學機電工程學院， 西安 710048)

復合材料加筋板結構廣泛應用于飛機的機翼、尾翼的翼面等部位[1]，其在服役過程中經常受拉伸、彎曲、剪切等外載荷作用而產生屈曲變形，并發生局部失穩，此時仍具有較強的后屈曲承載能力[2-3]。設計合理的復合材料加筋壁板結構可以顯著提高其可靠性、穩定性及承載能力。安全系數法是復合材料結構設計的經典方法之一，但由于缺少對隨機變量的定量分析，導致該方法不能準確反映出結構的安全程度[4-5]。復合材料結構中的幾何參數、材料屬性以及加工工藝參數等均存在不確定性，這些變量的不確定性對結構可靠性的影響不可忽略[6-7]。因此，對包含大量參數的復合材料加筋板進行可靠性優化設計時，需要研究各隨機因素對復合材料性能的影響程度。

目前，眾多研究學者結合數值分析、仿真模擬及實驗研究等手段，分析了隨機因素對復合材料性能的靈敏度。王佩艷等[8]分析了個各設計變量對目標響應值的靈敏度，在此基礎上，對復合材料加筋板結構進行了優化設計。Liu等[9]用總剛度矩陣和質量矩陣及其一階導數描述了靜力響應、特征值和特征向量的分析靈敏度分析方法，對復合材料層合板和殼體進行了靈敏度分析。Omairey等[10]研究了微觀尺度的幾何參數和材料性能不確定性對纖維增強復合材料的彈性和可靠性的影響及隨機參數的靈敏度。阮文斌等[11]對復合材料結構輸出位移和強度比的不確定性來源進行分析,得到輸入變量的全局靈敏度排序結果。

全局靈敏度分析方法可以定量分析出單個變量或多個變量相互作用對輸出響應的貢獻程度，并且具有全局性和穩定性，廣泛應用于可靠性分析、結構設計與優化等多個領域[12-14]。在復合材料不確定性分析方面，一些學者利用靈敏度分析法已經開展了探索性的研究，證明了該方法在復合材料結構分析方面的可行性[15]。現以復合材料加筋板為研究對象，選取屈曲載荷作為輸出響應值，以復合材料加筋板的材料屬性和幾何參數作為設計變量，基于Kriging方法的代理模型建立加筋板輸出響應與設計變量的函數關系，采用Sobol法求解復合材料加筋板材料屬性和幾何參數的靈敏度，研究結果可對復合材料加筋板可靠性設計提供指導。

1 復合材料加筋板的仿真分析

采用T形復合材料加筋板模型，如圖1所示，其中：長La=280 mm、寬Wa=160 mm、筋條間距Sa=100 mm、翼板寬度W=24 mm、腹板高度H=20 mm。

復合材料加筋板材料屬性如表1所示，蒙皮和筋條的鋪層分別為[0/90/±45]s和 [0/90/45/0/-45]s，蒙皮的總厚度為1 mm，加強筋的總厚度為1.25 mm。

基于ABAQUS軟件建立復合材料加筋壁板的仿真分析模型，為提高計算效率，模型中筋條和蒙皮均采用S4R單元進行建模，邊界條件：加載一端開放施加載荷方向自由度，另一端固定，側邊開放自由度，采用Buckle計算模塊對有限元模型進行特征值屈曲分析，其中第一階模態如圖2所示。

圖1 復合材料加筋板模型示意圖Fig.1 Schematic diagram of the composite reinforced plate model

表1 材料屬性Table 1 Material properties

圖2 T形復合材料加筋板的一階屈曲模態Fig.2 First-order flexion mode of T-shaped composite reinforced plates

通過仿真分析可得，T形復合材料加筋板的屈曲載荷為4.97 kN，與文獻[16]實驗數據4.89 kN相比，數值模擬數值與實驗數值誤差為1.6%，由此說明，研究模型是合理、準確的。

2 復合材料加筋板參數的靈敏度分析

基于有限元分析理論，在邊界條件確定時，復合材料加筋板的屈曲載荷可以表示成材料屬性和幾何參數的函數[17]為

Pcr=f1(E1,E2,G12,G13,G23,T,θ,W,H)

(1)

式(1)中：Pcr為復合材料加筋板屈曲載荷；T為復合材料單層板厚度；θ為纖維±45°的鋪層角度；W為翼板寬度；H為腹板高度。

由式(1)可知，Pcr的變化是由E1、E2、G12、G13、G23、T、θ、W、H改變引起的，為了有效地描述Pcr的變化，E1、E2、G12、G13、G23、T、θ、W、H均為隨機變量。復合材料加筋板中變量分布[17]如表2所示。

表2 變量分布Table 2 Distribution of variables

由于復合材料結構的復雜性，等參數與屈曲載荷之間存在著隱函數關系。在求解復合材料加筋板參數靈敏度時，需多次通過ABAQUS軟件進行數值模擬來計算屈曲載荷的數值。因此，進行復合材料加筋板參數靈敏度分析時就需要對屈曲載荷與多個參數之間的函數進行擬合。Kriging模型作為線性回歸的一種改進技術，包含了線性回歸部分和非參數部分[18]，具有很高的擬合精度同時可以簡化運算，提高計算效率。因此，采用Kriging方法來擬合屈曲載荷與各個參數之間的關系模型。

2.1 基于Kriging方法建立關系模型

(2)

式(2)中：C=[C1(X),C2(X),…,C9(X)]T為多項式函數；β=[β1,β2,…,β9]T為回歸系數；Z(X) 服從(0，δ2)的正態分布，其協方差為

Cov[Z(Xi),Z(Xj)]=δ2R(Xi,Xj)

(3)

式(3)中：δ為方差；R為m×m的相關矩陣，R(Xi,Xj) 為兩個任意試驗點Xi和Xj的相關函數，R選取Gaussian函數形式為

(4)

(5)

(6)

(7)

(8)

矩陣R可通過引入最小化函數來確定，進而求解出Kriging代理模型。

通過拉丁超立方采樣抽取500個實驗樣本，隨機挑選450組數據建立屈曲載荷與多個參數的Kriging的代理模型，將剩余的50組數據用于模型的精度檢測，實際值與預測值結果如圖3所示，最大相對誤差為7.76%，預測的平均誤差為3%。計算結果表明，使用Kriging建立的代理模型來預測復合材料加筋板的屈曲載荷是合理的。

圖3 屈曲載荷實際值與預測值結果Fig.3 Actual versus predicted flexion load results

2.2 基于Sobol法的靈敏度分析

Sobol靈敏度分析法通常用于非線性、非單調的函數模型，能夠便捷地計算出各輸入參數對響應值的一階、高階以及全局靈敏度系數，它可以表示為多個子函數相互組合的函數，即

(9)

式(9)中：x=(x1,x2,…,xn)為自變量，定義域為Kn={(x1,x2,…,xn)|0≤xi≤1,i=1,2,…,n}；f0為常數；fpj為子函數表達形式，每個子函數間均為正交，且在其對應定義域內積分為0，即

(10)

因此，子函數的計算表達式為

(11)

式(11)中：dxi和dxij分別為dx1dx2…dxn中不含dxi和dxidxj的乘積；Kn-1為不含xi的定義域；Kn-2為不含xi和xj的定義域。

通過計算，可以求解出函數的每一項偏方差和總方差為

(12)

(13)

在Sobol的方法中，輸入參數的靈敏度是由其對響應值總方差的貢獻率進行評價的。因此，一階響應指數Si和總一階響應指數STi為

(14)

(15)

式中：Si為一階靈敏度系數，反映了單一參數xi對復合材料加筋板屈曲載荷的影響程度；STi為全局靈敏度系數，反映綜合考慮輸入參數x1,x2,…,xn的變化情況時，各參數間交互作用對復合材料加筋板屈曲載荷的影響。

2.3 參數靈敏度分析結果

在得到Kriging代理模型Pcr(X)后，使用Sobol方法對復合材料加筋板的材料屬性和幾何參數進行靈敏度分析，分析流程如圖4所示，分析結果如圖5所示。

圖4 復合材料加筋板參數靈敏度分析流程Fig.4 Parameter sensitivity analysis process for composite reinforced sheet

圖5 復合材料加筋板參數靈敏度Fig.5 Parameter sensitivity of composite stiffened plates

復合材料加筋板的材料屬性和幾何參數對屈曲載荷的貢獻量通常由全局靈敏度系數進行評價，由圖5(b)可知，設計參數對屈曲載荷的全局靈敏度系數大小為ST6>ST8>ST5>ST1>ST9>ST2>ST4>ST3>ST7，X6與X8對應的一階靈敏度系數與全局靈敏度系數相差較大，說明參數復合材料單層板厚度T和翼板寬度W之間相互交互作用最為明顯。

根據圖5的復合材料加筋板參數靈敏度計算結果，可以看出，復合材料加筋板的單層板厚度和翼板的寬度對響應值屈曲載荷的影響較大，其中，單層板厚度對屈曲載荷的影響最大，剩余幾個變量的靈敏度數值較小，則表明這些參數對響應值的影響程度較小，因此，當設計目標為提高屈曲載荷時，應該重點關注復合材料加筋板的單層板厚度和翼板的寬度，盡量減少它們的隨機不確定性。

3 結論

(1)在有限元分析的基礎上，選取了9個變量為設計參數，基于Kriging方法的代理模型建立了屈曲載荷與材料屬性、幾何參數之間的函數表達式。

(2)運用Sobol靈敏度分析方法，對復合材料加筋板材料屬性和幾何參數進行靈敏度分析，結果顯示單層板厚度和翼板寬度對屈曲載荷影響較大，則說明單層板厚度、翼板寬度為關鍵參數，為復合材料加筋板的結構設計提供指導。