孫士平 張 冰 胡 政

南昌航空大學(xué)航空制造工程學(xué)院，南昌，330063

0 引言

航空結(jié)構(gòu)中復(fù)合材料層合板作為組成薄壁結(jié)構(gòu)的板殼類構(gòu)件承受著復(fù)雜載荷工況，其抗屈曲穩(wěn)定性能倍受重視，精確預(yù)測其屈曲性能成為高品質(zhì)結(jié)構(gòu)設(shè)計(jì)的必要前提。研究者采用各種數(shù)值方法開展了復(fù)合材料層合板的屈曲性能分析，對軸壓載荷屈曲問題進(jìn)行了深入研究，取得了顯著成果，但對剪切屈曲問題研究不多[1-6]。近期有TARJAN等[7]采用Ritz法計(jì)算了均布、線性變化軸壓及剪切載荷作用下正交各向異性復(fù)合材料層合長板的屈曲性能；YANG等[8]結(jié)合配點(diǎn)法推導(dǎo)出四邊固支正交各向異性方板剪切屈曲問題的一般解析解，簡化了計(jì)算過程；SELYUGIN等[9]基于有限元法分析了單/雙軸剪切載荷工況正交各向異性長板的剪切屈曲性能及彎扭耦合對剪切屈曲性能的影響；JUNG等[10]提出改進(jìn)8-ANS殼單元研究了剪切載荷和鋪層順序?qū)?fù)合材料層合殼屈曲行為的影響；劉慶輝等[11-12]和CHEN等[13]則分別采用Ritz法、Galerkin法分析了剪切、剪切與面內(nèi)線性變化載荷組合工況下，帶扭轉(zhuǎn)約束正交各向異性長板的屈曲問題，計(jì)算結(jié)果與有限元數(shù)值解相吻合。這些剪切屈曲研究主要以正交各向異性板為對象，研究層合板尺寸、彎扭耦合對剪切屈曲性能的影響規(guī)律，較少關(guān)注各向異性對稱復(fù)合材料層合板的剪切屈曲分析及鋪層順序優(yōu)化。

不同數(shù)值方法求解復(fù)合材料層合板屈曲問題的計(jì)算效率存在較大差異，應(yīng)用智能搜索算法開展復(fù)合材料層合板屈曲優(yōu)化特別是含耦合項(xiàng)而無封閉解的各向異性復(fù)合材料層合板屈曲優(yōu)化，需要耗費(fèi)大量的屈曲響應(yīng)分析次數(shù)，高效的數(shù)值方法會顯著減少優(yōu)化時(shí)間。廣義微分求積法[14](generalized differential quadrature,GDQ)是從微分求積法[15](DQ)發(fā)展而來的一種更穩(wěn)定高效的數(shù)值方法，已應(yīng)用于許多力學(xué)問題求解。KADKHODAYAN等[16]采用GDQ法分析了均質(zhì)矩形板的彈/塑性屈曲問題；李威等[17]采用GDQ法分析了受切向力作用的梁的穩(wěn)定性問題；DARVIZEH等[18]通過均布軸壓作用復(fù)合材料層合板屈曲問題的計(jì)算比較證明，GDQ法能比Ritz法更高效精確地求解復(fù)合材料層合板屈曲問題；MAAREFDOUST等[19]采用GDQ法分析了中厚斜型均質(zhì)板的剪切屈曲問題。采用GDQ法開展各向異性復(fù)合材料層合板的剪切屈曲分析與優(yōu)化研究尚未見報(bào)道。

本文采用GDQ法開展對稱復(fù)合材料層合板的剪切屈曲分析與鋪層優(yōu)化。針對GDQ法求解層合板屈曲問題存在計(jì)算精度差、不穩(wěn)定收斂現(xiàn)象，提出權(quán)重系數(shù)擾動策略來改進(jìn)GDQ法的計(jì)算收斂性，實(shí)現(xiàn)擾動GDQ法對層合板剪切屈曲的穩(wěn)定有效計(jì)算。在此基礎(chǔ)上，結(jié)合直接搜索模擬退火算法[20]開展復(fù)合材料層合板的鋪層順序優(yōu)化，比較邊界條件、長寬比、鋪層數(shù)和鋪層形式對剪切屈曲性能的影響，研究軸壓與剪切載荷組合工況下，不同長寬比和載荷比時(shí)層合板屈曲性能的變化規(guī)律。

1 復(fù)合材料層合板描述

1.1 基本方程

(a)層合板載荷邊界 (b)層合板鋪層結(jié)構(gòu)圖1 對稱層合板示意圖Fig.1 A schematic of symmetrical composite laminate

承受剪切和均勻軸壓載荷作用的對稱復(fù)合材料層合板如圖1所示，長a、寬b、厚h，總鋪層數(shù)為2L。層合板第l層鋪層的纖維角度為θl，第l層鋪層下表面到中面的距離為zl。

基于經(jīng)典層板理論的對稱復(fù)合材料層合板屈曲微分方程為

(1)

其中，w即w(x,y,t)為中面位移函數(shù)；t為時(shí)間；Nx和Ny分別為層合板在x和y方向承受的軸向載荷;Nxy為層合板在xy平面內(nèi)承受的剪切載荷；Dij(i,j=1, 2, 6)為層合板剛度矩陣元素, 其計(jì)算公式為

(2)

層合板的邊界條件常用S(simply supported，簡支)、C(clamped，固支)以及F(free，自由)4位字母組合表示，如邊界條件SCSC中第1和第3位字母S分別表示x=0和x=1邊為簡支邊界，第2和第4位字母C分別表示y=0和y=1邊為固支邊界。以層合板x=0邊為例，固支邊界為

(3)

簡支邊界為

(4)

1.2 控制方程的GDQ求解

GDQ法是一種求解常(偏)微分方程的數(shù)值計(jì)算方法。若一維函數(shù)f(x)在區(qū)間[0, 1]上連續(xù)可微，則其n階導(dǎo)數(shù)可用定義域內(nèi)任意N點(diǎn)的函數(shù)值線性加權(quán)組合表達(dá)：

(5)

i=1,2,…,Nn=1,2,…,N-1

(6)

i,j=1,2,…,Ni≠jn=1,2,…,N-1

對于二維函數(shù)g(x,y)，其導(dǎo)數(shù)可按下式計(jì)算：

(7)

i=1,2,…,Nj=1,2,…,Ms,q≥0

s+q=1,2,…,N+M-2

令層合板長寬比λ=a/b，取位移函數(shù)w(x,y,t)=W(x,y)eiωt，并設(shè)X=x/a、Y=y/b，將式(1)變換處理為

(8)

將式(7)代入式(8)得

(9)

將式(7)代入式(3)和式(4)并進(jìn)行歸一化處理，得到以X=0邊為例的邊界表達(dá)式，固支邊界：

(10)

簡支邊界：

(11)

記格柵點(diǎn)位移向量W=(Wb,Wd)T，其中下標(biāo)b表示邊界上及與其相鄰的(4N+4M-16)個格柵點(diǎn)，下標(biāo)d表示余下的(N-4)×(M-4)個內(nèi)部格柵點(diǎn)。利用式(10)、式(11)分別置換式(9)的第1、2、N×(M-1)、N×M方程實(shí)現(xiàn)邊界條件的施加，得到矩陣表達(dá)形式：

(12)

消除變量Wb，式(12)整理可得特征方程：

(13)

其中，Add和Bdd為(N-4)×(M-4)階方陣；Abb為4N+4(M-4) 階方陣；Adb和Abd分別為(N×M-4N-4M+16)×(4N+4M-16)和(4N+4M-16)×(N×M-4N-4M+16)階矩陣；η為屈曲載荷系數(shù)。求解特征方程(式(13))可得到最小特征值所對應(yīng)的臨界屈曲載荷系數(shù)ηcr，將其歸一化處理為量綱一臨界屈曲載荷系數(shù)k*：

(14)

D0=E2h3/[12(1-ν12ν21)]

式中，E2為鋪層在方向2的彈性模量；ν為泊松比。

2 擾動GDQ法

基于MATLAB平臺編程實(shí)現(xiàn)GDQ方法的計(jì)算，離散格柵點(diǎn)按以下Chebyshew多項(xiàng)式確定位置：

(15)

柵格點(diǎn)劃分如圖2所示，其特點(diǎn)是臨近邊界處的柵格點(diǎn)較密，而遠(yuǎn)離邊界處的柵格點(diǎn)少，既有利于邊界條件的實(shí)施，也能保證在柵格點(diǎn)較少的情況下得到較精確的計(jì)算結(jié)果。

2.1 GDQ法的收斂性改進(jìn)

圖2 歸一化設(shè)計(jì)域柵格點(diǎn)劃分示意圖Fig.2 Normalized design domain grid point division diagram

表1 復(fù)合材料性能參數(shù)

圖3 基于GDQ法的矩形板φ曲線Fig.3 Rectangular plate φ curve based on GDQ method

由圖3可知，隨著柵格點(diǎn)數(shù)的增加，無論復(fù)合材料層合板還是均質(zhì)板，GDQ法計(jì)算獲得的剪切屈曲對比指標(biāo)φ均未穩(wěn)定趨于1，都存在較大波動。φ的波動反映了GDQ法求解剪切屈曲問題存在精度差、計(jì)算振蕩不收斂問題。要利用GDQ法進(jìn)行剪切屈曲問題求解就必需采取有效措施改進(jìn)GDQ法的計(jì)算穩(wěn)定性和精度。

為此，提出權(quán)重系數(shù)擾動策略，對GDQ的n階導(dǎo)數(shù)權(quán)重系數(shù)矩陣的主對角線元素進(jìn)行擾動：

(16)

采用擾動策略的矩形板屈曲性能計(jì)算結(jié)果如圖4所示，可以看到隨著柵格點(diǎn)數(shù)的增加，無論復(fù)合材料層合板還是均質(zhì)板，擾動后的GDQ法(簡稱擾動GDQ法)計(jì)算獲得的φ穩(wěn)定趨于1，與圖3相比，僅用較少的柵格點(diǎn)數(shù)(如10×10)就能獲得穩(wěn)定收斂的計(jì)算結(jié)果，說明擾動GDQ法求解層合板屈曲問題的計(jì)算效率和穩(wěn)定性相比于GDQ法得到了有效改進(jìn)。

圖4 基于擾動GDQ法的矩形板φ曲線Fig.4 Rectangular plate φ curve based on perturbed GDQ method

圖5 擾動量δ對計(jì)算收斂的影響(λ=2)Fig.5 Influence of disturbance δ on computational convergence(λ=2)

2.2 擾動GDQ法的計(jì)算驗(yàn)證

采用擾動GDQ法分別計(jì)算均質(zhì)板和復(fù)合材料層合板的屈曲性能來驗(yàn)證計(jì)算結(jié)果的可靠性，計(jì)算中柵格點(diǎn)數(shù)取15×15。

考慮SSSS和CCCC邊界，不同λ時(shí)承受剪切力Nxy=-1N作用的均質(zhì)矩形板臨界屈曲載荷系數(shù)k*計(jì)算結(jié)果如表2所示。材料彈性模量E=200 GPa，泊松比ν=0.3?？梢钥闯?，擾動GDQ法計(jì)算結(jié)果與文獻(xiàn)[21]Ritz法的計(jì)算結(jié)果基本一致。

考慮剪切與均布軸壓載荷組合作用，λ=2時(shí)SSSS邊界斜交鋪設(shè)復(fù)合材料層合板[(±θ)2]s的臨界屈曲載荷系數(shù)k*受鋪層角θ變化的影響。剪切力Nxy=-1N，組合工況1為Nx∶Ny∶Nxy=1∶0∶1；組合工況2為Nx∶Ny∶Nxy=1∶1∶1。鋪層材料性能如表1所示，分別采用擾動GDQ法和FEM進(jìn)行計(jì)算獲得的結(jié)果如圖6所示。由圖6可知，不同組合工況下擾動GDQ法計(jì)算結(jié)果與FEM結(jié)果的最大相對差異僅1.29%，吻合度較好。另外，給定θ時(shí)工況1的k*要明顯大于工況2的k*，雖然兩種工況均在θ=0°時(shí)取得最小k*，但分別在50°和70°附近獲得最大k*且k*值差異高達(dá)2.3倍；兩種工況下通過變化θ均能有效提高k*，相對于對應(yīng)工況的最小k*，提高幅度高達(dá)300%和500%，說明鋪層角度優(yōu)化能顯著提高復(fù)合材料層合板的抗剪切屈曲能力。

表2 均質(zhì)矩形板k*的計(jì)算結(jié)果對比

圖6 復(fù)合材料層合板k*不同方法計(jì)算結(jié)果比較Fig.6 Comparison of calculation results of composite laminate k* using different methods

考慮單位長度正剪切Nxy=1N、負(fù)剪切Nxy=-1N載荷作用工況，λ=2的復(fù)合材料層合板[(±θ)2]s在CCCC、SSSS、CSCS、SCSC 4種邊界時(shí)，k*隨θ變化曲線如圖7所示。

圖7a顯示，負(fù)剪切工況時(shí)k*隨θ的變化為單調(diào)凸變化，不同邊界下k*均在θ=0°最小而在接近60°區(qū)域最大；而圖7b中，正剪切工況下k*隨θ的變化為非凸變化，不同邊界下k*均在θ=5°附近最小而在接近55°區(qū)域最大，此時(shí)θ變化對k*值的影響小于負(fù)剪切工況時(shí)的影響；另外，圖7中兩種剪切載荷下，SCSC和CSCS邊界的k*值均位于CCCC和SSSS的k*值之間，且在θ從0°趨于90°的過程中，SCSC的k*逼近CCCC的k*，而CSCS的k*逼近SSSS的k*。由此可見，鋪層角θ對剪切屈曲性能影響明顯，邊界條件變化不會改變鋪層角對剪切屈曲的影響規(guī)律，但正剪切載荷會使層合板屈曲性能的鋪層順序優(yōu)化求解復(fù)雜化。

(a)Nxy=-1N

(b)Nxy=1N圖7 λ=2時(shí)不同邊界下鋪層角θ對k*的影響Fig.7 Influence of lay angle θ on k* at different boundary when λ=2

通過上述算例的計(jì)算數(shù)據(jù)比較分析，說明采用擾動策略的擾動GDQ法能有效進(jìn)行矩形板的剪切屈曲問題分析，獲得的計(jì)算結(jié)果準(zhǔn)確可靠。

3 復(fù)合材料層合板的屈曲優(yōu)化

總層數(shù)為2L的對稱復(fù)合材料層合板[θ1/θ2/…/θL-1/θL]s，鋪層材料參數(shù)如表1所示，板長a=500 mm，t0=0.125 mm。以鋪層角θl(l=1, 2,…,L)為設(shè)計(jì)變量，最大化層合板k*的優(yōu)化模型：

(17)

式(17)為多極值離散優(yōu)化問題，本文采用自適應(yīng)直接搜索模擬退火算法求解，離散變量θl的增量取Δθ=1°，每個問題進(jìn)行40次優(yōu)化統(tǒng)計(jì)計(jì)算結(jié)果。

3.1 剪切載荷工況

基于式(17)進(jìn)行單位長度剪切力Nxy=-1N作用復(fù)合材料層合板[(±θ)2]s的優(yōu)化，獲得λ為0.5～3時(shí)，SSSS、SCSC、CCCC 3種邊界下的優(yōu)化結(jié)果，如圖8所示。圖8顯示，隨著λ的增大，3種邊界下k*均呈單調(diào)遞減，且降速隨λ增大而逐漸減??；其中SCSC邊界下k*隨λ增大，從接近SSSS的k*而逐漸趨于接近CCCC的k*，說明對于復(fù)合材料層合長板，長邊的邊界約束顯著影響屈曲性能，而短邊的邊界約束對屈曲性能影響較小，可以忽略。同時(shí)，當(dāng)λ從0.5增大到3，3種邊界對應(yīng)最優(yōu)角度從略大于30°逐漸趨于60°，除SSSS為單調(diào)遞增外，SCSC在λ<1時(shí)最優(yōu)角度增大較快，在λ=1達(dá)到54°后隨λ增大而平緩增長。另外，因邊界條件和結(jié)構(gòu)的對稱性，SSSS和CCCC邊界下λ與1/λ的層合板最優(yōu)鋪層角互為補(bǔ)角。

圖8 3種邊界下的優(yōu)化結(jié)果比較Fig.8 Comparison of optimization results for 3 boundaries

考慮任意鋪設(shè)對稱復(fù)合材料層合板[θ1/θ2/θ3/θ4]s在SSSS、SCSC、CCCC邊界和不同λ時(shí)的優(yōu)化結(jié)果如表3所示。表3數(shù)據(jù)顯示，不同λ時(shí)，3種邊界的優(yōu)化鋪層均具有相同或相近的角度值，并隨λ從0.25增大到3，最優(yōu)鋪層角從近30°遞增趨于60°，而k*先快速下降后趨于平緩減小，與圖8的規(guī)律相同；另外，表3中SSSS邊界優(yōu)化鋪層順序與文獻(xiàn)[22]優(yōu)化結(jié)果基本相同，差異在于文獻(xiàn)[22]取增量Δθ=0.1°而更精確。進(jìn)一步優(yōu)化λ=2時(shí)，SSSS邊界鋪層順序?yàn)閇θ1/θ2/…/θL-1/θL]s的2L層對稱復(fù)合材料層合板，獲得最優(yōu)鋪層順序近似為[(57)L]s，對應(yīng)k*=93.20。以上結(jié)果說明，剪切載荷作用的對稱復(fù)合材料層合板鋪層順序優(yōu)化中，鋪設(shè)形式與鋪層數(shù)的影響較小，可簡化為單變量優(yōu)化問題來提高優(yōu)化效率；層合板長寬比對剪切屈曲的影響隨長寬比增大而減小，當(dāng)長寬比大于2后，優(yōu)化鋪層角接近60°。計(jì)算結(jié)論與復(fù)合材料結(jié)構(gòu)設(shè)計(jì)手冊[23]的相關(guān)論述一致，說明擾動GDQ法求解剪切屈曲是有效可行的。

表3 剪切載荷作用對稱復(fù)合材料層合板屈曲優(yōu)化結(jié)果

3.2 剪切與軸壓載荷組合工況

考慮復(fù)合材料層合板[θ1/θ2/θ3/θ4]s承受剪切與軸壓載荷組合工況，分X向軸壓組合工況(Nx∶Ny∶Nxy=1∶0∶1)和雙軸壓組合工況(Nx∶Ny∶Nxy=1∶1∶1)，其中Nxy=-1N。在SSSS、SCSC、CCCC 3種邊界下優(yōu)化獲得不同λ時(shí)的k*和鋪層順序，如表4、表5所示。表4和表5數(shù)據(jù)顯示，不同邊界和載荷工況下，隨著λ從0.25增大到1，k*均先快速減小，在λ>1后k*趨于平緩減小，但雙軸壓組合工況的k*要小于X向軸壓組合工況，僅為對應(yīng)值的0.4～0.8，說明更復(fù)雜載荷會弱化層合板的抗屈曲能力。從最優(yōu)鋪層看，組合工況的優(yōu)化鋪層比剪切載荷工況更復(fù)雜，僅X向軸壓組合工況λ<1時(shí)各邊界的優(yōu)化鋪層角為正值且相近，而兩種組合工況在SSSS和CCCC邊界獲得的最優(yōu)鋪層角有正有負(fù)但絕對值相近；另外，雙軸壓組合工況CCCC邊界λ=1時(shí)最優(yōu)鋪層中21°鋪層角的出現(xiàn)顯得較為特殊。

假如剪切載荷與軸壓載荷具有如下比例關(guān)系Nxy=βNx(β≥0)，Nx=-1 N，變化載荷比β，優(yōu)化獲得λ=2時(shí)兩種組合工況下SSSS邊界復(fù)合材料層合板的屈曲優(yōu)化結(jié)果如表6所示?？梢钥吹?，兩種工況下，優(yōu)化鋪層隨β增大而趨于剪切載荷工況優(yōu)化結(jié)果；在β≥0.3時(shí)，k*隨β增大而減小且在β≥1處快速下降后降速趨緩，而當(dāng)β<0.3時(shí)，k*表現(xiàn)為非單調(diào)變化，存在一個最佳的k*。因此可以在軸壓載荷工況中適當(dāng)匹配剪切載荷來改善層合板的抗屈曲能力。進(jìn)一步優(yōu)化不同λ時(shí)X向軸壓組合工況SSSS邊界復(fù)合材料層合板，獲得k*隨β變化曲線如圖9所示。圖9顯示，復(fù)合材料層合板的抗屈曲能力整體上隨著剪切載荷的增大而降低，除λ=0.5外，僅當(dāng)β∈(0,0.5)時(shí),k*存在一定的非單調(diào)波動；而λ≥2后,λ對k*的影響小到可以忽略。

表4 X向軸壓組合工況復(fù)合材料層合板屈曲優(yōu)化結(jié)果(Nx∶Ny∶Nxy=1∶0∶1)

表5 雙軸壓組合工況復(fù)合材料層合板屈曲優(yōu)化結(jié)果(Nx∶Ny∶Nxy=1∶1∶1)

表6 λ=2時(shí)SSSS邊界復(fù)合材料層合板屈曲優(yōu)化結(jié)果

圖9 X向軸壓組合工況SSSS邊界層合板k*隨β變化曲線Fig.9 Variation curve of laminate k*with β under SSSS boundary and X axial compression combined loading

4 結(jié)論

(1)GDQ法求解復(fù)合材料層合板剪切屈曲問題的振蕩不收斂的現(xiàn)象源于GDQ載荷矩陣的奇異性，本文提出的權(quán)重系數(shù)擾動策略能有效改善載荷矩陣的奇異性，實(shí)現(xiàn)了擾動GDQ法對層合板剪切屈曲問題的穩(wěn)定快速求解。

(2)剪切載荷作用時(shí)，層合板鋪層角對屈曲性能的影響在負(fù)剪切載荷時(shí)為非凸變化，而在正剪切載荷時(shí)為凸變化；層合板鋪層數(shù)、鋪設(shè)方式和邊界條件對優(yōu)化鋪層角的影響較小，可簡化為單變量優(yōu)化問題；而隨著層合板長寬比增大，剪切屈曲性能逐漸減弱，優(yōu)化鋪層角趨于60°。

(3)剪切與軸壓載荷組合作用時(shí)，除較小的剪切載荷有助于改善層合板屈曲性能外，層合板屈曲性能隨剪切載荷增大而減小，優(yōu)化鋪層趨于剪切載荷優(yōu)化結(jié)果，而當(dāng)層合板長寬比大于2后，剪切載荷對屈曲性能的影響可以忽略。