應用擾動廣義微分求積法的復合材料層合板剪切屈曲分析與優化

孫士平 張 冰 胡 政

南昌航空大學航空制造工程學院，南昌，330063

0 引言

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

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

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

1 復合材料層合板描述

1.1 基本方程

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

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

基于經典層板理論的對稱復合材料層合板屈曲微分方程為

(1)

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

(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法是一種求解常(偏)微分方程的數值計算方法。若一維函數f(x)在區間[0, 1]上連續可微，則其n階導數可用定義域內任意N點的函數值線性加權組合表達：

(5)

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

(6)

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

對于二維函數g(x,y)，其導數可按下式計算：

(7)

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

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

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

(8)

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

(9)

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

(10)

簡支邊界：

(11)

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

(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)階矩陣；η為屈曲載荷系數。求解特征方程(式(13))可得到最小特征值所對應的臨界屈曲載荷系數ηcr，將其歸一化處理為量綱一臨界屈曲載荷系數k*：

(14)

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

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

2 擾動GDQ法

基于MATLAB平臺編程實現GDQ方法的計算，離散格柵點按以下Chebyshew多項式確定位置：

(15)

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

2.1 GDQ法的收斂性改進

圖2 歸一化設計域柵格點劃分示意圖Fig.2 Normalized design domain grid point division diagram

表1 復合材料性能參數

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

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

為此，提出權重系數擾動策略，對GDQ的n階導數權重系數矩陣的主對角線元素進行擾動：

(16)

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

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

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

2.2 擾動GDQ法的計算驗證

采用擾動GDQ法分別計算均質板和復合材料層合板的屈曲性能來驗證計算結果的可靠性，計算中柵格點數取15×15。

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

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

表2 均質矩形板k*的計算結果對比

圖6 復合材料層合板k*不同方法計算結果比較Fig.6 Comparison of calculation results of composite laminate k* using different methods

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

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

(a)Nxy=-1N

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

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

3 復合材料層合板的屈曲優化

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

(17)

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

3.1 剪切載荷工況

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

圖8 3種邊界下的優化結果比較Fig.8 Comparison of optimization results for 3 boundaries

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

表3 剪切載荷作用對稱復合材料層合板屈曲優化結果

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

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

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

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

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

表6 λ=2時SSSS邊界復合材料層合板屈曲優化結果

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

4 結論

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

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

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