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)
其中，zl为第l层下表面到中面的距离；
为第l层铺层的变换刚度系数，由铺层折减刚度按与铺层角θl相关的坐标变换获得。
层合板的边界条件常用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,…,N n=1,2,…,N-1
其中，N为整个定义域栅格点的数目；xk(k=1,2,…,N)为任意选取的栅格点；权重系数
为拉格朗日多项式的n阶导数，其表达式为
(6)
i,j=1,2,…,N i≠j n=1,2,…,N-1
其中，
为拉格朗日多项式，其一阶导数
对于二维函数g(x,y)，其导数可按下式计算：
(7)
i=1,2,…,N j=1,2,…,M s,q≥0
s+q=1,2,…,N+M-2
其中，N、M分别为定义域内x、y方向栅格点的数目；权重系数
分别为对应x、y方向拉格朗日多项式的s阶导数和q阶导数。
令层合板长宽比λ=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法的收敛性改进
考虑四边固支(CCCC)矩形板，a=500 mm，λ=2，h=1 mm，单位长度上承受剪切力Nxy=-1N，Nx=Ny=0，采用GDQ法和有限元法(FEM)进行矩形板剪切屈曲计算。均质板材料弹性模量E=210 GPa，泊松比υ=0.3；T300/5208环氧树脂基碳纤维复合材料层合板[30/-45/60/-75]s的材料参数见表1，密度为ρ，单铺层厚度t0=0.125 mm。FEM采用ANSYS软件，均质板选用shell63单元，复合材料层合板选用shell99单元，单元大小均为25 mm×25 mm，计算获得均质板和复合材料层合板的临界剪切屈曲载荷系数
分别为10.235和78.531。记
为对比指标，
为GDQ法计算获得的矩形板临界剪切屈曲载荷系数，如果φ趋于1则说明GDQ法与FEM具有相同计算精度。图3给出了四边固支复合材料和均质矩形板取不同栅格点数(N×M)时φ的变化曲线。
图2 归一化设计域栅格点划分示意图
Fig.2 Normalized design domain grid point division diagram
表1 复合材料性能参数
Tab.1 Compositeperformance parameters
E1E2E3ν12ν23181 GPa10.3 GPa10.3 GPa0.280.3G12G23G13v13ρ7.17 GPa3.78 GPa7.17 GPa0.281 600 kg/m3
图3 基于GDQ法的矩形板φ曲线
Fig.3 Rectangular plate φ curve based on GDQ method
由图3可知，随着栅格点数的增加，无论复合材料层合板还是均质板，GDQ法计算获得的剪切屈曲对比指标φ均未稳定趋于1，都存在较大波动。φ的波动反映了GDQ法求解剪切屈曲问题存在精度差、计算振荡不收敛问题。要利用GDQ法进行剪切屈曲问题求解就必需采取有效措施改进GDQ法的计算稳定性和精度。
分析四边简支(SSSS)均质板在栅格点数取7×7时的一阶权重系数矩阵c(1)和载荷矩阵
发现，c(1)的主对角线上存在非常接近于零的元素(d(1)存在同样现象)，从而造成B矩阵接近于奇异而不能准确求逆，导致计算结果误差大而不收敛。
为此，提出权重系数扰动策略，对GDQ的n阶导数权重系数矩阵的主对角线元素进行扰动：
(16)
其中，δ为扰动参数。取δ=10-4，扰动后的
及Br如下：
可看到的主对角线上原接近于零的元素增大为扰动量，显著减小了及Br中元素的数量级差，有效地改善了Br矩阵的奇异性，通过提高计算精度实现了特征方程的稳定求解。
采用扰动策略的矩形板屈曲性能计算结果如图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)
考虑CCCC均质板，扰动参数δ取不同值对剪切屈曲计算收敛的影响如图5所示。扰动量为10-6时，因对中元素的数量级差改变较小而不能有效提高计算精度导致不收敛；扰动量大于10-6后，随着栅格点数增加计算结果均能稳定收敛；当扰动量大于10-2时，虽然计算能快速稳定收敛，但因扰动量偏大而影响Br矩阵数据的准确性，导致计算偏差。扰动量为10-4、10-3时能较好地协调计算稳定性和计算精度，本文取10-4。
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*的计算结果对比
Tab.2 Comparison of calculation results of homogeneous
rectangular plate k*
SSSSCCCC扰动GDQRitz[21]扰动GDQRitz[21]λ1.09.3259.32514.6414.641.27.9847.98312.6412.641.47.2887.28711.7311.731.66.9086.90711.2511.251.86.6896.68210.6410.642.06.5466.54610.2510.252.26.2946.29310.0310.032.46.1046.1039.909.902.65.9775.9769.779.772.85.8945.8939.629.623.05.8415.8409.549.54
图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*均呈单调递减，且降速随λ增大而逐渐减小；其中SCSC边界下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 剪切载荷作用对称复合材料层合板屈曲优化结果
Tab.3 Buckling optimization results of symmetrical composite laminate under shear loading
SSSSSCSCCCCCk*最优铺层之一优化铺层[22]k*最优铺层之一k*最优铺层之一λ0.251 356.92[31/323]s1 350.78[314]s2 170.56[304]s0.5371.96[(32/33)2]s376.42[332/32/33]s572.91[31/323]s1.0131.83[454]s[454]s169.93[542/532]s186.47[454]s2.093.20[(57/58)2]s[57.54]s142.48[594]s143.23[59/583]s3.086.48[594]s[58.8/58.9/58.82]s137.29[592/60/59]s137.73[60/59/602]s
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所示。可以看到，两种工况下，优化铺层随β增大而趋于剪切载荷工况优化结果；在β≥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)
Tab.4 Buckling optimization results of composite laminate under X axial compression combined loading
(Nx∶Ny∶Nxy=1∶0∶1)
SSSSSCSCCCCCk*最优铺层之一k*最优铺层之一k*最优铺层之一λ0.25281.23[24]s281.46[2/32/2]s985.58[15/14/152]s0.572.43[8/7/82]s74.16[9/10/92]s249.23[132/14/16]s1.035.27[452/-452]s58.61[45/44/-25/-23]s72.92[353/34]s2.034.64[472/-48/47]s53.59[493/-34]s59.20[44/43/44/40]s3.034.28[48/47/-49/48]s54.20[46/45/41/-30]s56.62[47/462/47]s
表5 双轴压组合工况复合材料层合板屈曲优化结果(Nx∶Ny∶Nxy=1∶1∶1)
Tab.5 Buckling optimization results of composite laminate under biaxial compression combined loading
(Nx∶Ny∶Nxy=1∶1∶1)
SSSSSCSCCCCCk*最优铺层之一k*最优铺层之一k*最优铺层之一λ0.25219.74[20/-21/20/19]s221.18[21/-16/22/21]s569.82[32/31/-31/-32]s0.555.32[20/-19/19/22]s56.02[20/-15/22/21]s145.70[302/-312]s1.018.35[(45/-45)2]s36.16[622/-46/-45]s39.20[21/752/74]s2.013.85[69/-75/69/65]s35.32[592/-49/-46]s36.42[602/-592]s3.013.82[68/-75/662]s35.20[59/60/-582]s35.87[59/58/-59/-56]s
表6 λ=2时SSSS边界复合材料层合板屈曲优化结果
Tab.6 Buckling optimization results of composite laminate
with SSSS boundary and λ=2
Nx∶Ny∶Nxy=1∶0∶ β Nx∶Ny∶Nxy=1∶1∶ βk*最优铺层之一k*最优铺层之一β036.801[-43/473]s13.908[-74/66/67/73]s0.1037.067[44/-463]s13.918[72/-70/-66/-67]s0.2037.039[45/-462/45]s13.920[71/-69/-74/-66]s0.3036.834[45/-47/44/-47]s13.918[70/-71/-72/-79]s0.5036.472[46/-47/462]s13.898[69/-73/-75/-62]s0.7535.760[462/-492]s13.878[70/-71/71/86]s1.034.460[472/-48/47]s13.847[69/-75/69/65]s1.531.667[494]s13.756[70/71/-70/77]s2.028.164[504]s13.615[702/72/-65]s5.015.356[532/522]s10.932[654]s108.489[(53/54)2]s7.166[613/60]s
图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后，剪切载荷对屈曲性能的影响可以忽略。
[1] 唐俊. 剪切载荷作用下复合材料层合板优化分析[J]. 复合材料学报, 1980, 6(4): 48-56.
TANG Jun. On Optimum Analysis of Laminated Plates under Shear[J]. Acta Materiae Compositae Sinica, 1980, 6(4): 48-56.
[2] GRENESTEDT J L. Layup Optimization against Buckling of Shear Panels[J]. Structural & Multidisciplinary Optimization, 1991, 3(2):115-120.
[3] FUKUNAGA H, SEKINE H. A Laminate Design for Elastic Properties of Symmetric Laminates with Extension Shear or Bending-Twisting Coupling[J]. Journal of Composite Materials, 1994, 28(8): 708-731.
[4] KOSTELETOS S. Shear Buckling Response of Laminated Plates[J]. Composite Structures, 1992, 20(3): 147-154.
[5] LOUGHLAN J. The Shear Buckling Behaviour of Thin Composite Plates with Particular Reference to the Effects of Bend-Twist Coupling[J]. International Journal of Mechanical Sciences, 2001, 43(3): 771-792.
[6] LOPATIN A V, KORBUT Y B. Buckling of Clamped Orthotropic Plate in Shear[J]. Composite Structures, 2006, 76(1/2): 94-98.
[7] TARJAN G, SAPKAS A, KOLLAR L P. Stability Analysis of Long Composite Plates with Restrained Edges Subjected to Shear and Linearly Varying Loads[J]. Journal of Reinforced Plastics and Composites, 2010, 29(9):1386-1398.
[8] 杨端生, 黄炎, 李广利. 各向异性矩形板的剪切屈曲分析[J]. 应用力学学报, 2012, 29(2):220-225.
YANG Duansheng, HUANG Yan, LI Guangli. Shear Buckling Analysis of Anisotropic Rectangular Plates[J]. Chinese Journal of Applied Mechanics, 2012, 29(2):220-225.
[9] SELYUGIN S. On Choice of Optimal Anisotropy of Composite Plates against Buckling, with Special Attention to Bending-Twisting Coupling[J]. Structural and Multidisciplinary Optimization, 2013, 48(2): 279-294.
[10] JUNG W Y, HAN S C. Shear Buckling Responses of Laminated Composite Shells Using a Modified 8-node ANS Shell Element[J]. Composite Structures, 2014, 109(1): 119-129.
[11] 刘庆辉, 乔丕忠, 郭兴文, 等. 正交各向异性板剪切屈曲载荷的近似表达式研究[J]. 工程力学, 2014, 31(3): 32-38.
LIU Qinghui, QIAO Pizhong, GUO Xingwen, et al. Approximate Formulas for Shear Buckling of Orthotropic Plates[J]. Engineering Mechanics, 2014, 31(3): 32-38.
[12] LIU Q H, QIAO P Z, GUO X W. Buckling Analysis of Restrained Orthotropic Plates under Combined In-plane Shear and Axial Loads and Its Application to Web Local Buckling[J]. Composite Structures, 2014,111(11): 540-552.
[13] CHEN Q Y, QIAO P Z. Shear Buckling of Rotationally-restrained Composite Laminated Plates[J]. Thin-Walled Structures, 2015, 94: 147-154.
[14] SHU C, DU H J. A Generalized Approach for Implementing General Boundary Condition in the GDQ Free Vibration Analysis of Plates[J]. International Journal of Solids and Structures, 1997, 34(7): 837-846.
[15] SHU C. Differential Quadrature and Its Application in Enginerring[M]. London: Springer, 2000:95-120.
[16] KADKHODAYAN M, MAAREFDOUST M. Elastic/Plastic Buckling of Isotropic Thin Plates Subjected to Uniform and Linearly Varying In-plane Loading Using Incremental and Deformation Theories[J]. Aerospace Science and Technology, 2014, 32(1):66-83.
[17] 李威，曾志松, 韩旭. 弹性地基上受切向力作用梁稳定性研究[J]. 振动与冲击, 2014, 33(8): 192-195.
LI Wei, ZENG Zhisong, HAN Xu. Stability of a Beam on Elastic Foundation Subjected to a Tangential Force[J]. Journal of Vibration and Shock, 2014, 33(8): 192-195.
[18] DARVIZEH M, DARVIZEH A, ANSARI R, et al. Buckling Analysis of Generally Laminated Composite Plates[J]. Composite Structures, 2004, 63(1): 69-74.
[19] MAAREFDOUST M, KADKHODAYAN M. Elastic/Plastic Buckling Analysis of Skew Plates under In-plane Shear Loading with Incremental and Deformation Theories of Plasticity by GDQ Method[J]. Journal of the Brazilian Society of Mechanical Sciences and Engineering, 2015, 37(2):761-776.
[20] 孙士平, 曾庆龙, 吴建军. 基于改进模拟退火算法的复合材料层合板屈曲优化[J].中国机械工程, 2015, 26(12): 1676-1683.
SUN Shiping, ZENG Qinglong, WU Jianjun. Buckling Optimization of Composite Laminates Based on Improved Simulated Annealing[J]. China Mechanical Engineering, 2015, 26(12): 1676-1683.
[21] SMITH S T, BRADFORD M A, OEHLERS D J. Elastic Buckling of Unilaterally Constrained Rectangular Plates in Pure Shear[J]. Engineering Structures, 1999, 21(5): 443-453.
[22] KIM C W, LEE J S. Optimal Design of Laminated Composite Plates for Maximum Buckling Load Using Genetic Algorithm[J]. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 2005, 219(8):869-878.
[23] 中国航空研究院.复合材料结构设计手册[M]. 北京：航空工业出版社, 2001:348-353.
Chinese Aeronautics Research Institute. Handbook of Composites Structures Design[M]. Beijing: Aviation Industry Press, 2001:348-353.
