應(yīng)用擾動(dòng)廣義微分求積法的復(fù)合材料層合板剪切屈曲分析與優(yōu)化

孫士平 張 冰 胡 政

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

0 引言

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

不同數(shù)值方法求解復(fù)合材料層合板屈曲問題的計(jì)算效率存在較大差異，應(yīng)用智能搜索算法開展復(fù)合材料層合板屈曲優(yōu)化特別是含耦合項(xiàng)而無封閉解的各向異性復(fù)合材料層合板屈曲優(yōu)化，需要耗費(fèi)大量的屈曲響應(yīng)分析次數(shù)，高效的數(shù)值方法會(huì)顯著減少優(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法開展對(duì)稱復(fù)合材料層合板的剪切屈曲分析與鋪層優(yōu)化。針對(duì)GDQ法求解層合板屈曲問題存在計(jì)算精度差、不穩(wěn)定收斂現(xiàn)象，提出權(quán)重系數(shù)擾動(dòng)策略來改進(jìn)GDQ法的計(jì)算收斂性，實(shí)現(xiàn)擾動(dòng)GDQ法對(duì)層合板剪切屈曲的穩(wěn)定有效計(jì)算。在此基礎(chǔ)上，結(jié)合直接搜索模擬退火算法[20]開展復(fù)合材料層合板的鋪層順序優(yōu)化，比較邊界條件、長寬比、鋪層數(shù)和鋪層形式對(duì)剪切屈曲性能的影響，研究軸壓與剪切載荷組合工況下，不同長寬比和載荷比時(shí)層合板屈曲性能的變化規(guī)律。

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

1.1 基本方程

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

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

基于經(jīng)典層板理論的對(duì)稱復(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

對(duì)于二維函數(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)個(gè)格柵點(diǎn)，下標(biāo)d表示余下的(N-4)×(M-4)個(gè)內(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))可得到最小特征值所對(duì)應(yīng)的臨界屈曲載荷系數(shù)ηcr，將其歸一化處理為量綱一臨界屈曲載荷系數(shù)k*：

(14)

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

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

2 擾動(dòng)GDQ法

基于MATLAB平臺(tái)編程實(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ì)算獲得的剪切屈曲對(duì)比指標(biāo)φ均未穩(wěn)定趨于1，都存在較大波動(dòng)。φ的波動(dòng)反映了GDQ法求解剪切屈曲問題存在精度差、計(jì)算振蕩不收斂問題。要利用GDQ法進(jìn)行剪切屈曲問題求解就必需采取有效措施改進(jìn)GDQ法的計(jì)算穩(wěn)定性和精度。

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

(16)

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

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

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

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

采用擾動(dòng)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。可以看出，擾動(dòng)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所示，分別采用擾動(dòng)GDQ法和FEM進(jìn)行計(jì)算獲得的結(jié)果如圖6所示。由圖6可知，不同組合工況下擾動(dòng)GDQ法計(jì)算結(jié)果與FEM結(jié)果的最大相對(duì)差異僅1.29%，吻合度較好。另外，給定θ時(shí)工況1的k*要明顯大于工況2的k*，雖然兩種工況均在θ=0°時(shí)取得最小k*，但分別在50°和70°附近獲得最大k*且k*值差異高達(dá)2.3倍；兩種工況下通過變化θ均能有效提高k*，相對(duì)于對(duì)應(yīng)工況的最小k*，提高幅度高達(dá)300%和500%，說明鋪層角度優(yōu)化能顯著提高復(fù)合材料層合板的抗剪切屈曲能力。

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

圖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í)θ變化對(duì)k*值的影響小于負(fù)剪切工況時(shí)的影響；另外，圖7中兩種剪切載荷下，SCSC和CSCS邊界的k*值均位于CCCC和SSSS的k*值之間，且在θ從0°趨于90°的過程中，SCSC的k*逼近CCCC的k*，而CSCS的k*逼近SSSS的k*。由此可見，鋪層角θ對(duì)剪切屈曲性能影響明顯，邊界條件變化不會(huì)改變鋪層角對(duì)剪切屈曲的影響規(guī)律，但正剪切載荷會(huì)使層合板屈曲性能的鋪層順序優(yōu)化求解復(fù)雜化。

(a)Nxy=-1N

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

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

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

總層數(shù)為2L的對(duì)稱復(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°，每個(gè)問題進(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*，說明對(duì)于復(fù)合材料層合長板，長邊的邊界約束顯著影響屈曲性能，而短邊的邊界約束對(duì)屈曲性能影響較小，可以忽略。同時(shí)，當(dāng)λ從0.5增大到3，3種邊界對(duì)應(yīng)最優(yōu)角度從略大于30°逐漸趨于60°，除SSSS為單調(diào)遞增外，SCSC在λ<1時(shí)最優(yōu)角度增大較快，在λ=1達(dá)到54°后隨λ增大而平緩增長。另外，因邊界條件和結(jié)構(gòu)的對(duì)稱性，SSSS和CCCC邊界下λ與1/λ的層合板最優(yōu)鋪層角互為補(bǔ)角。

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

考慮任意鋪設(shè)對(duì)稱復(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層對(duì)稱復(fù)合材料層合板，獲得最優(yōu)鋪層順序近似為[(57)L]s，對(duì)應(yīng)k*=93.20。以上結(jié)果說明，剪切載荷作用的對(duì)稱復(fù)合材料層合板鋪層順序優(yōu)化中，鋪設(shè)形式與鋪層數(shù)的影響較小，可簡化為單變量優(yōu)化問題來提高優(yōu)化效率；層合板長寬比對(duì)剪切屈曲的影響隨長寬比增大而減小，當(dāng)長寬比大于2后，優(yōu)化鋪層角接近60°。計(jì)算結(jié)論與復(fù)合材料結(jié)構(gòu)設(shè)計(jì)手冊(cè)[23]的相關(guān)論述一致，說明擾動(dòng)GDQ法求解剪切屈曲是有效可行的。

表3 剪切載荷作用對(duì)稱復(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向軸壓組合工況，僅為對(duì)應(yīng)值的0.4～0.8，說明更復(fù)雜載荷會(huì)弱化層合板的抗屈曲能力。從最優(yōu)鋪層看，組合工況的優(yōu)化鋪層比剪切載荷工況更復(fù)雜，僅X向軸壓組合工況λ<1時(shí)各邊界的優(yōu)化鋪層角為正值且相近，而兩種組合工況在SSSS和CCCC邊界獲得的最優(yōu)鋪層角有正有負(fù)但絕對(duì)值相近；另外，雙軸壓組合工況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)變化，存在一個(gè)最佳的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)波動(dòng)；而λ≥2后,λ對(duì)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ù)擾動(dòng)策略能有效改善載荷矩陣的奇異性，實(shí)現(xiàn)了擾動(dòng)GDQ法對(duì)層合板剪切屈曲問題的穩(wěn)定快速求解。

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

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