壓剪聯(lián)合載荷作用下復(fù)合材料壁板屈曲及后屈曲性能計算與優(yōu)化方法研究

李 政，梁 珂

(西北工業(yè)大學(xué)航空學(xué)院, 西安 710072)

0 引言

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

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

本文基于改進(jìn)的Koiter攝動理論[14]，提出一種復(fù)合材料壁板在壓剪聯(lián)合載荷作用下非線性屈曲計算的有限元降階方法，快速精確地獲得結(jié)構(gòu)屈曲及后屈曲性能指標(biāo)，并將其用于壁板鋪層優(yōu)化，獲得滿足結(jié)構(gòu)各類型性能設(shè)計要求的最優(yōu)鋪層設(shè)計。

1 復(fù)合材料壁板屈曲及后屈曲性能計算

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

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

(1)

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

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

(2)

NxKxx+NyKyy+NxyKxy)

(3)

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

對復(fù)合材料鋪層結(jié)構(gòu)而言，其本構(gòu)關(guān)系為

(4)

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

繼續(xù)計算單元應(yīng)變能的三階和四階導(dǎo)數(shù)即可獲得三維張量Q和四維張量C，然后將所有單元計算得到L，Q和C對整個結(jié)構(gòu)進(jìn)行組裝，式(1)的顯示表達(dá)即可獲得。

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

(5)

(6)

截斷式(6)中的高階項即可得到結(jié)構(gòu)在未變形狀態(tài)點處的降階模型

(7)

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

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

圖1 壓剪聯(lián)合載荷作用下復(fù)合材料壁板的承載響應(yīng)曲線Fig.1 The response for the composite panel under combined compression and shear loads

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

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

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

2 復(fù)合材料壁板屈曲及后屈曲性能優(yōu)化

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

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

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

設(shè)計變量：θ=[θ1,θ2,…,θi,…]

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

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

上述優(yōu)化問題中的性能約束需要依據(jù)實際設(shè)計要求來決定是否施加以及如何施加。

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

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

3 算例分析

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

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

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

表1 兩種載荷工況

3.1 壁板屈曲及后屈曲分析

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

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

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

3.2 4鋪層壁板屈曲及后屈曲性能鋪層優(yōu)化

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

表3 不同載荷工況下以非線性屈曲載荷值為目標(biāo)的最優(yōu)鋪層信息

表4 在不同載荷工況下以后屈曲承載剛度值為目標(biāo)的最優(yōu)鋪層信息

表5 不同載荷工況下以承載剛度殘余系數(shù)值為目標(biāo)的最優(yōu)鋪層信息

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

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

圖4 4鋪層壁板各性能指標(biāo)最優(yōu)鋪層在屈曲點處的變形圖(載荷工況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鋪層壁板屈曲及后屈曲性能鋪層優(yōu)化

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

圖5 16鋪層壁板各性能指標(biāo)最優(yōu)鋪層在屈曲點處的變形圖(載荷工況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 不同載荷工況下以非線性屈曲載荷值為目標(biāo)的最優(yōu)鋪層信息

表8 不同載荷工況下以后屈曲承載剛度值為目標(biāo)的最優(yōu)鋪層信息

表9 不同載荷工況下以承載剛度殘余系數(shù)值為目標(biāo)的最優(yōu)鋪層信息

4 結(jié)論

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

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

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

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