壓剪聯合載荷作用下復合材料壁板屈曲及后屈曲性能計算與優化方法研究

李 政，梁 珂

(西北工業大學航空學院, 西安 710072)

0 引言

復合材料由于具有較好的比強度和比剛度，以及極度自由的可設計性，在航天航空結構工程領域得到廣泛應用。目前飛行器結構設計不僅實現了復合材料使用比重的顯著增加，其使用范圍也逐漸從受力較小的構件向主承力結構進行轉變。復合材料壁板是組成飛行器結構的重要承力構件，此類結構通常會受到軸壓和剪切載荷的聯合作用，容易在材料失效破壞之前發生屈曲失穩，造成結構承載性能下降[1-2]。故而，研究復合材料壁板的壓剪承載性能，揭示鋪層參數對壁板屈曲載荷和后屈曲承載性能的影響機理，對于提高復合材料在飛行器結構中的應用至關重要[3-5]。

即便采用當前先進的力學數值算法技術，針對復合材料壁板開展屈曲和后屈曲性能優化仍然存在種種困難[6]。首先，屈曲載荷和后屈曲性能分別由結構的面外和面內剛度水平控制，單獨優化其中一個性能指標會引起另一個指標的明顯下降；其次，計及幾何非線性效應的結構屈曲分析計算量過大，常規非線性有限元方法已不適合用于優化迭代過程，不利于揭示屈曲及后屈曲優化的力學機理。早先，Pandey等[7]推導了計算復合材料受壓板非線性屈曲性能的解析公式并用于鋪層優化，發現單獨優化后屈曲性能會明顯降低屈曲載荷值。部分學者采用一種計算非線性屈曲響應的漸近數值方法(ANM)[8-11]來快速計算結構的屈曲及后屈曲響應。Raju等[12-13]成功將ANM方法用于軸壓復合材料層合板的后屈曲優化；并基于變分原理和Rayleigh-Ritz方法，提出了一種用于復合材料層合板后屈曲分析的高效方法，采用層壓參數作為設計變量，通過兩步優化策略實現了后屈曲性能優化。

本文基于改進的Koiter攝動理論[14]，提出一種復合材料壁板在壓剪聯合載荷作用下非線性屈曲計算的有限元降階方法，快速精確地獲得結構屈曲及后屈曲性能指標，并將其用于壁板鋪層優化，獲得滿足結構各類型性能設計要求的最優鋪層設計。

1 復合材料壁板屈曲及后屈曲性能計算

采用常規的有限元離散化方法，復合材料壁板結構在壓剪聯合載荷作用下的非線性平衡方程可以用一組含N個未知數的非線性代數方程組f(q)=λfext來表示。其中，N是有限元離散模型的總自由度，矢量q表示結構位移場，向量f和fext分別表示結構內力和外力，荷載系數λ表示外荷載的加載變化。常規非線性有限元方法采用Newton-Raphson算法直接求解上述平衡方程，即可得到壁板的非線性承載響應，即q-λ曲線。

本文基于改進的Koiter攝動理論來構造上述有限元模型的降階模型，該降階模型在壁板的未變形狀態處(q0,λ0)建立，用來近似非線性平衡方程。已知狀態(q0,λ0)附近的未知構型(q,λ)可以用場變化量(u,μ)表示，即q=u+q0和λ=μ+λ0。為了簡化平衡方程的非線性程度并求解場變化量(u,μ)，使用近似泰勒展開的方法將原先的平衡方程f(q)=λfext在已知構型(q0,λ0)處，關于位移變化量u展開至三階項，即得

(1)

式中，二維張量L、三維張量Q和四維張量C分別是壁板內力在已知構型(q0,λ0)處的一次、二次和三次展開項。方程(1)的右端考慮了用于激發屈曲分叉變形的擾動載荷，使所提出的方法能夠分析屈曲問題。在荷載矩陣F中，第一列向量f1=fext為外荷載，其他列向量fα，α=2,…,m+1表示各擾動荷載。擾動載荷是由壁板的密集屈曲模態和幾何剛度的乘積計算獲得，故而擾動載荷的數量m由壁板密集屈曲模態的數量決定。φ是載荷矩陣F中各載荷列向量系數所構成的載荷系數向量，其第一個分量為μ，其余分量均為0。

式(1)中的L，Q和C可以通過計算單元應變能相對于自由度的高階導數來獲得。本文采用基于Von Kármán運動學的四節點四邊形復合材料板殼單元來實現有限元分析[14]。復合材料單元的內力f和切線剛度L可分別通過計算單元應變能的梯度和Hessian矩陣得到

(2)

NxKxx+NyKyy+NxyKxy)

(3)

式中，K0為單元的線彈性剛度矩陣，A表示單元面積，Cm為材料矩陣，Nx，Ny和Nxy表示單元內力矢量N的3個分量，Nnl為單元內力矢量N的非線性部分，Bl，Bnl，Kxx,Kyy,Kxy為單元的幾何關系矩陣。

對復合材料鋪層結構而言，其本構關系為

(4)

式中，N，M分別為單元的薄膜內力和內力彎矩，A和D分為是薄膜剛度和彎曲剛度。

繼續計算單元應變能的三階和四階導數即可獲得三維張量Q和四維張量C，然后將所有單元計算得到L，Q和C對整個結構進行組裝，式(1)的顯示表達即可獲得。

仿照非線性平衡方程的展開方式，將位移和載荷系數向量u和φ基于攝動參數ξ進行展開可得

(5)

(6)

截斷式(6)中的高階項即可得到結構在未變形狀態點處的降階模型

(7)

上述降階模型(7)本質上是一個具有(1+m)個未知數的非線性方程組，可以選擇路徑跟蹤技術(例如弧長法)來求解該降階模型，進而得到載荷系數μ與廣義位移ξ之間的關系，然后再利用關系式λ=μ+λ0，q=u+q0和位移展開式(5)即可獲得壁板非線性承載響應的非線性預測解(q-λ)。

受壓/剪壁板面內位移隨載荷的響應曲線呈現出典型的“雙線性”曲線特征(圖1)，即前屈曲階段基本為線性承載，屈曲發生后曲線出現轉折，后屈曲階段出現明顯的非線性特性。求解上述在未變形狀態處建立的降階模型，獲得的非線性預測解能夠準確跟蹤壁板的前屈曲和初始后屈曲承載響應，壁板的屈曲和后屈曲性能指標可以通過分析該非線性預測解得到。本文所關心的3個指標，即非線性屈曲載荷、后屈曲承載剛度、承載剛度殘余系數的具體計算方式如下：

圖1 壓剪聯合載荷作用下復合材料壁板的承載響應曲線Fig.1 The response for the composite panel under combined compression and shear loads

1)非線性屈曲載荷Pcr：該載荷值可通過探測壁板面內位移隨載荷變化曲線的切線拐點獲得。當非線性預測曲線在某載荷水平下的斜率值出現突變時，該載荷即為壁板的非線性屈曲載荷。

2)后屈曲承載剛度Kpost：該剛度為非線性預測曲線在屈曲載荷點處的斜率值。壁板后屈曲承載的非線性特性顯著，用初始后屈曲剛度來表征壁板在后屈曲階段的總體承載表現。

3)承載剛度殘余系數kres：該系數值體現了壁板屈曲后承載剛度的折減情況。具體計算公式為kres=Kpost/Kpre，Kpre是壁板前屈曲階段的承載剛度。Kres越接近1，表明壁板屈曲后承載能力損失的越小。

2 復合材料壁板屈曲及后屈曲性能優化

上一節提出了復合材料壁板在壓剪聯合載荷作用下屈曲和后屈曲承載性能計算的快速分析方法，即攝動有限元降階法。本節的目標是將采用該方法計算得到的屈曲/后屈曲承載性能指標應用于壁板鋪層優化，以找到復合材料壁板各種最優性能所對應的鋪層方式θ。

在鋪層性能優化中，復合材料壁板的總厚度設定為常數，各鋪層的纖維鋪向角度θi為離散整數設計變量，則屈曲及后屈曲性能優化問題設計如下

最大化：Pcr(θ)或Kpost(θ)或kres(θ)

設計變量：θ=[θ1,θ2,…,θi,…]

鋪向角約束：-π/2≤θi≤π/2

性能約束(可選)：L.B.≤Pcr(θ)或Kpost(θ)或kres(θ)≤U.B.

上述優化問題中的性能約束需要依據實際設計要求來決定是否施加以及如何施加。

本文采用MATLAB優化工具箱來實現復合材料壁板的鋪層優化，其所嵌套的遺傳優化算法(GA)是一種成熟的優化方法，因此在下面的數值算例測試中采用了GA優化器所默認的參數設置，并用增廣拉格朗日遺傳算法(ALGA)和外部罰函數法求解非線性約束問題。使用攝動有限元降階方法來計算壁板屈曲及后屈曲性能指標的速度非常快，因此本文采用獲得全局最優解性能較佳的遺傳類算法是一種合適的選擇。

本文的算例考核工作采用計算機配置為8核的i7 CPU、16G內存，所提取的最優解均滿足優化收斂條件，即GA優化算法自然停止并在窗口輸出 “exitflag=1”的正常收斂信號，每個優化問題采用GA算法重復求解至少3次。

3 算例分析

受壓剪聯合載荷作用的復合材料壁板如圖2所示，其長和寬分別為105 mm和75 mm，材料常數為E1=157 362 MPa，E2=10 092 MPa，G12=G13=5 321 MPa，G23=5 321 MPa， μ12=0.277， 復合材料鋪層的單層厚度為0.125 mm。本算例采用14×10的復合材料四邊形板單元來對壁板進性網格剖分。壁板采用四邊簡支的約束條件，即四條邊均約束面外(z向)位移，右端約束面內軸向(x向)位移并在邊中點處約束面內的另外一個方向(y向)位移。

圖2 壓剪聯合載荷作用下的復合材料壁板Fig.2 The composite panel under combined compression and shear loads

如表1和圖2所示，壁板考慮兩種不同的壓剪復雜載荷條件：1)載荷工況1采用面內均勻壓縮載荷和剪切載荷，壓剪載荷的施加比例為1∶0.2；2)載荷工況2采用非均勻的壓縮載荷和均勻剪切載荷。注意，兩種載荷工況下，均勻壓載和不均勻壓載所施加的載荷總量相當。

表1 兩種載荷工況

3.1 壁板屈曲及后屈曲分析

先采用第2節提出的攝動有限元降階方法分析某給定鋪層方式的復合材料壁板。假設4鋪層壁板的纖維方向為[15,30]s，按照載荷工況2施加不均勻面內壓載和均勻剪切載荷。依照第2節中介紹的分析步驟，在未變形狀態點處建立降階模型時需基于壁板的前2階密集屈曲模態，得到的降階模型總自由度數為3。通過求解降階模型可以獲得壁板屈曲響應的非線性預測解，這里取面內壓縮位移隨加載變化曲線，如圖3所示。圖3中的兩條響應曲線為分別采用基于常規非線性有限元方法和本文介紹的攝動有限元降階方法計算獲得，壁板在屈曲點處的變形也在圖3中給出。可知，采用兩種方法所獲得的壁板承載響應曲線從未變形狀態到初始后屈曲狀態吻合得較好。兩種方法計算獲得的屈曲及后屈曲性能指標在表2中列出，兩種方法計算出來的性能指標值也較為一致。在計算效率方面，采用本文的方法開展單次非線性屈曲計算的CPU時間在2 s以內，而采用常規非線性有限元方法則需要10 s 左右，計算效率顯著提升。

圖3 復合材料壁板非線性屈曲的承載響應曲線(載荷工況2)Fig.3 The nonlinear response for the composite panel

表2 兩種方法計算的壁板屈曲/后屈曲性能指標比較

3.2 4鋪層壁板屈曲及后屈曲性能鋪層優化

針對4鋪層壁板，采用對稱鋪設方式[θ1,θ2]s，鋪層優化選取兩個鋪層的纖維角度θ1,θ2作為優化變量，在兩種不同的載荷工況下暫不考慮性能約束，分別單獨以非線性屈曲載荷、后屈曲承載剛度和承載剛度殘余系數為優化目標開展鋪層優化工作。優化計算得到的最優鋪層信息及最優鋪層對應的屈曲及后屈曲性能指標如表3～5所示。可以看出，單獨優化后屈曲承載剛度會使屈曲載荷值顯著下降，較最優屈曲載荷下降約50%。反之，最優屈曲載荷對應的后屈曲承載剛度僅為最優值的16%。面內壓載荷的不均勻施加會使屈曲載荷的最優值下降約14%，但最優后屈曲剛度卻稍有提高。可見，屈曲載荷值優化對面內壓載荷的不均勻性更為敏感。總體來說，屈曲載荷值最優鋪層的纖維鋪向角集中在45°附近，后屈曲承載剛度值最優鋪層的纖維鋪向角集中在0°附近。經過優化的承載剛度殘余系數能達到1附近，表明采用此類鋪層的壁板在屈曲后基本沒有承載剛度的損失。在計算效率方面，得益于本文提出的屈曲/后屈曲響應快速計算方法，單次優化計算僅約需要90 min左右。

表3 不同載荷工況下以非線性屈曲載荷值為目標的最優鋪層信息

表4 在不同載荷工況下以后屈曲承載剛度值為目標的最優鋪層信息

表5 不同載荷工況下以承載剛度殘余系數值為目標的最優鋪層信息

然后，以鋪層[15,30]s作為優化前的參考鋪層，將載荷工況2下參考鋪層和各性能最優鋪層發生屈曲時的變形圖繪制在圖4中。最后，在后屈曲承載剛度的優化過程中引入屈曲載荷約束，即Pcr>250 N，得到的最優鋪層如表6所示。可以看到，當引入屈曲載荷約束后，獲得的最優解均滿足約束要求，可以得到滿足一定抗屈曲性能的最優后屈曲承載剛度鋪層結構，該鋪層的兩個纖維角度變量分別在45°和0°附近。

表6 不同載荷工況下以后屈曲承載剛度值為目標的最優鋪層信息-考慮屈曲載荷約束

圖4 4鋪層壁板各性能指標最優鋪層在屈曲點處的變形圖(載荷工況2)Fig.4 The deformation of the optimal 4-ply panel at the buckling point for different performance index under the second load condition

3.3 16鋪層壁板屈曲及后屈曲性能鋪層優化

針對16鋪層壁板同樣采取對稱鋪層方式[θ1,θ2,θ3,θ4,θ5,θ6,θ7,θ8]s,鋪層優化選取θ1,θ2,θ3,θ4,θ5,θ6,θ7,θ8作為優化變量進一步豐富了優化空間，在兩種載荷工況下分別以非線性屈曲載荷、后屈曲承載剛度和承載剛度殘余系數作為優化目標進行優化任務。優化計算得到的優化纖維鋪層信息以及對應的屈曲和后屈曲性能指標如表7～8所示。以屈曲載荷為優化目標的計算會引起后屈曲剛度大幅下降，相較于最優后屈曲剛度下降了85%。最優后屈曲剛度對應鋪層的屈曲載荷與最優屈曲載荷相比下降了64%。該部分結果與4鋪層結果對照發現，改變壁板厚度并不能緩解屈曲載荷和后屈曲剛度的矛盾。比較兩種工況下的優化結果，不均勻載荷對屈曲載荷仍有10%的影響，但是對后屈曲剛度的影響已經變得很小。最優屈曲載荷對應的鋪層集中在45°左右，最優后屈曲剛度對應的鋪層集中在0°左右，最優承載剛度殘余系數對應的鋪層內層角度較大而外層角度較小。雖然可行解數量增加為原來的1806倍，但是單次優化計算用時仍在10～20 h范圍內。載荷工況1下的各優化鋪層發生屈曲時的變形圖如圖5所示。

圖5 16鋪層壁板各性能指標最優鋪層在屈曲點處的變形圖(載荷工況1)Fig.5 The deformation of the optimal 16-ply panel at the buckling point for different performance index under the first load condition

表7 不同載荷工況下以非線性屈曲載荷值為目標的最優鋪層信息

表8 不同載荷工況下以后屈曲承載剛度值為目標的最優鋪層信息

表9 不同載荷工況下以承載剛度殘余系數值為目標的最優鋪層信息

4 結論

本文針對壓剪聯合載荷作用下的復合材料壁板，采用非線性攝動有限元降階方法高效精準地計算了壁板的屈曲和后屈曲性能指標，并將其成功嵌套在鋪層性能優化流程，獲得滿足性能設計要求的最優鋪層信息,得出以下結論：

1)與傳統非線性有限元方法比較，本文采用攝動有限元降階方法所獲得的壁板承載響應曲線從未變形狀態到初始后屈曲狀態吻合得較好，計算獲得的屈曲及后屈曲性能指標精度較高，單次非線性屈曲計算的CPU時間在2 s以內，計算效率提升約80%。

2)單獨優化后屈曲承載剛度會使屈曲載荷值顯著下降，反之亦然。屈曲載荷值優化對面內壓載荷的不均勻性更為敏感。屈曲載荷值最優鋪層的纖維鋪向角集中在45°附近，后屈曲承載剛度值最優鋪層的纖維鋪向角集中在0°附近。經過優化的承載剛度殘余系數能達到1附近，表明采用此類鋪層的壁板在屈曲后基本沒有承載剛度的損失。

3)在后屈曲承載剛度優化中引入屈曲載荷約束后，可以得到滿足一定抗屈曲性能要求的最優后屈曲承載剛度鋪層結構，該鋪層的兩個纖維角度變量分別在45°和0°附近。