基于改进模拟退火算法的复合材料层合板屈曲优化

孙士平　曾庆龙　吴建军

南昌航空大学,南昌,330063

基于改进模拟退火算法的复合材料层合板屈曲优化

孙士平曾庆龙吴建军

南昌航空大学,南昌,330063

针对复合材料层合板的屈曲优化问题，提出一种改进的直接搜索模拟退火算法来求解最大化临界屈曲载荷系数的铺层顺序设计问题。改进算法引入搜索范围动态调整的新点产生方式实现全局搜索和局部搜索的协调，提高了算法的计算稳定性和计算效率。以离散铺层角为设计变量，采用里兹法进行屈曲响应分析以考虑弯扭耦合影响，通过典型多极值屈曲优化问题分析比较了算法改进措施的有效性。不同角度增量、铺层数、长宽比和载荷比的层合板屈曲优化结果表明，改进算法能有效地进行层合板铺层顺序优化。

屈曲优化;模拟退火算法;层合板;里兹法;优化设计

0　引言

复合材料层合板在航空、汽车、船舶等工程领域被广泛用于板、壳等结构中,屈曲失稳成为其不容忽视的失效形式之一,层合板的屈曲性能表现因此倍受关注。然而,与传统金属板壳相比,各向异性特性、铺层层数和铺层角度的离散化特征,使得层合板的屈曲分析与优化更加复杂困难。较多的层合板屈曲优化设计以正交各向异性板和0°、±30°、±45°、±60°、90°等几种离散铺层角为研究对象，从优化算法和屈曲响应分析两方面研究问题的解决途径。通过铺层角度连续化，利用序列线性规划[1]、可行方向法[2]等数学规划法实现层合板屈曲优化,但连续解的离散化会导致设计结果非最优化或不满足约束,且优化易陷入局部极值。遗传算法(GA)[3-6]、模拟退火(SA)[7-9]、粒子群(PSO)[10]、蚁群(ACO)[11]等随机算法能较好地求解离散型层合板屈曲优化问题,其目标函数评价成本远高于数学规划法。Erdal等[9]采用具有记忆功能的直接搜索模拟退火(DSA)[12]算法开展了层合板屈曲优化研究，考虑了角度增量为10°、15°、30°时不同离散角度对优化结果的影响，但DSA算法易陷入局部解。Karakay等[13]则比较了GA、SA和ACO算法求解层合板屈曲和频率优化问题时的性能表现。其他如分层优化方法[14]、分散搜索算法[15]、进化算法[16]及两层次优化方法[10]等也相继被用于层合板屈曲问题的铺层顺序设计。

另一方面，准确求解层合板屈曲响应，实现目标函数的精确计算，可以提高优化效率、减小误差,但因存在弯扭耦合故很难获得屈曲控制方程的精确解。一些研究通过忽略弯扭耦合以简化屈曲响应分析来优化铺层顺序,并结合限定弯扭耦合相对大小[5,9,10,13,15,16]、铺层铺设形式[3,6,8,9,11,16]等附加约束来减小计算偏差,这必然会对最优解的获得造成阻碍。屈曲问题的近似求解方法如里兹法[14]和有限元法[1-2,4,7]可以考虑弯扭耦合，提供满足精度要求的计算结果。有限元法能求解复杂结构的屈曲响应,但与随机算法结合的计算成本不容忽视,研究者常采用近似模型替代原问题模型来提高计算效率[5],而里兹法能方便有效地求解规则形状和常规边界条件的层合板屈曲响应问题。

针对层合板屈曲优化的铺层顺序设计问题,本文以规则简支板为研究对象,以离散铺层角度为设计变量，通过改进DSA算法的新解产生方式提高算法全局收敛性和稳定性，采用里兹法进行层合板屈曲响应分析来考虑弯扭耦合对优化结果的影响,从而实现了层合板屈曲载荷系数最大化的铺层顺序优化。此外，本文还研究了不同角度增量下具有更大离散设计空间时改进算法的性能表现和优化解的变化规律，并比较层合板不同长宽比、载荷比以及不同铺层数对优化结果的影响。

1　优化问题描述

图1所示对称层合板由N=2n层铺层组成,长宽为a×b,每层铺层具有相同厚度t0,第k层的纤维铺层角度为θk(k=1,2,…,n),则层合板铺层顺序可表示为[θ1/θ2/…/θk/…/θn]s。图1中,Nx、Ny分别为x和y方向作用在板中面的压力载荷。

图1　双轴向载荷作用下的对称层合板示意图

轴压载荷作用下对称层合板的屈曲控制方程为[17]

(1)

i,j=1,2,6

(2)

c=cosθs=sinθ

Q11=E1/(1-ν12ν21)

Q12=ν12E2/(1-ν12ν21)Q66=G12

式中,E1、E2分别为单层材料在1、2主轴方向的弹性模量；G12为1-2平面内的剪切模量;ν12、ν21分别为1-2平面内的纵向和横向泊松比。

从式(2)可以看到,由于铺层偏轴刚度与铺层角θ具有复杂的函数关系,造成屈曲优化问题具有多极值特征,故使得以铺层角θ为设计变量的铺层顺序设计较为困难。

四边简支板的边界条件为

(3)

式中，Mx、My分别为x轴和y轴方向的力矩。

上述边界条件中,对于弯扭耦合项D16、D26为零或其值相对较小的几种层合板,可采用封闭形式求解屈曲载荷系数λb[3]:

λb(m,l)=

(4)

其中,m、l为x、y轴向的半波数(m,l=1,2,…),不同m、l组合对应不同λb(m,l),相应屈曲载荷为λbNx和λbNy。

而D16和D26不为零时,因法向位移不能通过变量分离技术获得封闭形式解析解,故常采用里兹法近似求解。根据里兹法,取层合板中面z向位移函数：

(5)

其中,cml为待定系数,M=L为事先给定的正整数,一般在[1,10]内取值。

层合板的弯曲应变能V为

(6)

中面载荷做功为

(7)

基于最小势能原理有：

(8)

由此建立关于cml的M×L个线性方程组,问题表现为方程组系数矩阵的特征值问题,通过子空间迭代法等标准特征值求解方法获得一系列特征值,即屈曲载荷系数λb(m,l)。而临界屈曲载荷系数λcb对应层合板能抵抗屈曲的最小临界屈曲载荷,是一系列λb(m,l)值中的最小值，即

λcb=minλb(m,l)

(9)

最小临界屈曲载荷则通过公式Nx,cr=λcbNx和Ny,cr=λcbNy求得。

层合板屈曲优化问题就是要最大化临界屈曲载荷系数λcb,以层合板铺层角θ=(θ1,θ2,…,θk,…,θn)T(k=1,2,…,n)为设计变量,临界屈曲载荷系数最大化问题的数学表达为

(10)

2　模拟退火算法的改进

SA算法是一种模拟高温金属冷却过程的全局随机搜索算法[18],基于单点串行搜索并以一定概率接受差点,因此会丢失优化过程中获得的优化点。DSA算法[12]区别于SA算法之处在于DSA算法拥有一个点集合,按两种机制产生的新点被接受后仅替换该集合中的最差点,实现了优化过程最佳点的记忆保存。DSA算法的新点产生混合机制覆盖了全局粗搜索与局部细搜索，但降低了算法搜索效率和优化稳定性，易导致算法陷于局部极值[9]。本文针对层合板屈曲优化问题,通过改进DSA算法的初始温度确定方法和新点修正方法，增加动态调整的新点产生方式，提高算法求解层合板铺层顺序优化的稳定性和计算效率。

2.1初始集合和初始温度确定

首先在搜索空间S内按均匀分布随机产生20n(n为铺层角度个数)个点的集合,若集合中包含可行解则初始温度T0=500n,否则T0=2000n,使得T0与变量规模和问题难易程度相关,从而确定一个相对高且计算经济的初始温度值。因式(10)为无约束问题,因此T0=500n。继续从该集合中选出目标函数值较小的7(n+1)个点构成初始点集合A,并标记A中对应的最大目标函数值fH和最小目标函数值fL。这样避免了DSA算法以概率1抽样确定高初始温度所耗费的目标评价时间和因初始温度高而增加的迭代收敛过程，从而有助于提高算法计算效率。

2.2新点产生方式

在DSA算法中新点θ按混合机制产生:

(11)

式中,P为给定概率，本文取0.5;ρ为(0,1)区间均匀随机数;θ1为A中最佳点,θj(j=2,…,n+1)为A中随机选取的n个点;U(S)为搜索区间S内按均匀分布产生的随机点。

当新点θ在搜索区间之外时,重复式(11)直到产生满足要求的新点。

(12)

式中,γ为搜索半径调节系数;I为n维自适应变换列向量;“∘”表示Hadamard乘积;y为n维单位列向量;ωk为(0,5)内均匀分布随机数;Tk为第k次循环迭代的温度;ai为[0,1]内随机数;β为自适应概率,本文取0.1;rand(n)是按标准正态分布随机得到的n维列向量。

当产生的新点θ落在搜索空间S之外时,按下式修正θ的各分量θi:

(13)

测试数据表明式(13)可以提高计算效率且不影响优化结果。

2.3新点接受概率

新点θ的接受概率Pa计算式为

(14)

式中,f(θ)为新点θ的目标函数值;T为当前温度。

如果新点θ被接受,则替代A中最差点(fH所对应点),并更新fH和fL。

2.4降温策略

DSA算法采用指数形式降温策略：

Tk+1=αk+1Tkk≥0

(15)

降温系数αk+1∈[αmin,αmax](k≥1)按下式计算：

(16)

Lk=10n+10n(1-e(fL-fH))

2.5收敛准则

DSA算法采取双迭代收敛条件：

(17)

其中,ε1和ε2为较小实数,本文取ε1=ε2=0.001,即迭代温度足够低和A中的点足够集中时算法终止。

3　算例分析

基于MATLAB平台编程实现了DSA算法和改进DSA算法,并分别采用这两种算法进行临界屈曲载荷系数最大化的铺层顺序优化设计,比较算法改进效果。铺层材料为Graphite/Epoxy[6],E1=127.5575 GPa,E2=13.0316 GPa,G12=6.4124 GPa,ν12=0.30,厚度t0=0.127 mm。每个问题独立运行50次。

算例1考虑四边简支、两轴向承压的64层对称层合板,长a=50.8 mm,宽b=a/2,厚度t=8.128 mm,层合板承受面内压力Nx=0.175 N/mm,Ny=l1Nx,l1为两轴向载荷大小比例。假设该层合板为对称均衡层合板,由铺层角为0°2、±45°或90°2(即Δθ=45°)的双层铺层组铺设而成。l1=1时DSA算法和改进DSA算法基于式(4)和里兹法获得的优化结果分别见表1和表2,表中Nfa为平均目标函数评估次数,Nit为平均迭代次数,成功率为所获优化解与理论最优值的绝对误差小于设定精度(本文设为0.5)的次数与50次统计次数之比。

表1中DSA算法、改进DSA算法以及SA算法[9]的优化结果为独立运行50次统计所得,而GA算法[6]的优化数据是独立运行200次统计结果。结果显示改进DSA算法获得了与SA算法[9]、GA算法[6]相同的7个全局最优解，而DSA算法仅获得了其中3个,且GA算法[6]获得的优化解③和⑥实为局部极值，说明本文改进措施提高了DSA算法的全局搜索能力,能有效求解层合板屈曲优化问题。

表1　基于式(4)的λcb优化结果(l1=1)

表2数据显示DSA算法与改进DSA算法虽然均获得3个优化铺层，但改进DSA仅用了DSA所需Nfa和Nit的约1%的计算成本，获得了高于DSA的收敛成功率和优于DSA算法的λcb优化结果,说明本文改进措施有效避免了DSA算法易陷入局部极值和迭代后期收敛缓慢的不足,提高了DSA算法的计算效率和稳定性。另外,表1中基于式(4)所获优化铺层按里兹法计算获得的λcb基本各不相同,且均小于表2中基于里兹法获得的λcb,说明层合板屈曲优化结果受响应分析方法计算精度的影响明显,考虑弯扭耦合的里兹法提高了计算精度,能获得更好的优化结果。

假设上述64层对称层合板由单层铺层铺设而成,且θk∈[0°,90°],Δθ分别为30°、15°、10°、5°,l1分别为0.25、0.5、1.0、2.0时改进DSA算法求解式(10)的优化结果见表3和表4。

表3　基于式(4)的λcb优化结果(l1=0.25)

表4　不同l1时基于里兹法的λcb优化结果

算例2四边简支两轴向承压的对称层合板,长a=50.8 mm,长宽比为a/b,Nx=0.175 N/mm,Ny=l1Nx。铺层角度取值范围为[-90°,90°],总铺层数为N。基于里兹法和改进DSA算法求解式(10),获得单轴(l1=0)和双轴(l1=1)向压力作用下,不同N、a/b及Δθ时的优化结果(见表5～表8),表中Nopt、Iopt分别为50次独立运算中获得全局优化解的个数和次数。

表3中基于式(4)的优化结果与文献[9]的最优铺层完全一致,但其对应的λcb与按里兹法计算的λcb差异较大,与表4中l1=0.25时的优化结果相比,可发现基于里兹法获得的λcb均高于基于式(4)所获最优铺层的里兹法解,说明对于弯扭耦合不可忽视的层合板,式(4)因忽略弯扭耦合产生较大计算偏差而不能获得优化解,考虑弯扭耦合的里兹法能有效获得问题优化解。

表4数据显示，l1一定时λcb随Δθ减小而增大，说明增大设计空间能获得更好的优化结果。随着l1增大,λcb逐渐减小且受Δθ的影响减弱,层合板趋于用较大的铺层角来增强窄边承载能力以提高屈曲稳定性。

表5为l1=0的优化结果,可以看到,给定N和a/b时,Δθ分别取45°、15°、5°均获得相同的优化铺层,获得的λcb随N减小而减小，随a/b增大而增大,当a/b=0.5时0°铺层角最优,而a/b≥1时优化铺层角为45°和-45°的不同组合。

表6～表8则给出l1=1时N分别取32、64和128时的优化结果，数据显示,N一定时,λab随a/b增大而增大,铺层角随a/b增大逐渐从接近x轴向分布趋于接近y轴向分布,当a/b=1时的优化铺层角为45°和-45°的不同组合。

表5　基于里兹法的λcb优化结果(l1=0)

表7　基于里兹法的λcb优化结果(l1=1,N=64)

表8　基于里兹法的λcb优化结果(l1=1,N=32)

另外,从表5～表8中数据发现，给定N、a/b和Δθ下的优化问题50次独立运算所获得的λcb差异均较小,相对波动幅值小于1%,说明改进DSA算法具有较好的计算稳定性，能稳定地进行层合板屈曲问题的铺层顺序优化；而给定N和a/b时,随着Δθ减小,Iopt逐步减小,表6～表8中的λcb及其波动幅值逐渐增大,说明小Δθ值增大了设计空间有利于获得更好的优化解,但也因设计空间的扩大降低了算法搜索到全局优化解的概率和算法的计算稳定性。

4　结论

本文采用改进的直接搜索模拟退火算法开展了最大化层合板屈曲承载力的铺层顺序优化研究。改进的直接搜索模拟退火算法通过调整初始温度确定方法和引入搜索范围动态调整的新点产生方式,有效地改善了直接搜索模拟退火算法的计算效率和稳定性，提高了算法全局搜索能力。典型多极值屈曲优化结果验证了算法改进措施的有效性，揭示了考虑层合板弯扭耦合能更好地获得问题的优化结果,而不同的角度增量、铺层数以及长宽比层合板在不同比例轴压载荷作用下的优化结果表明，本文改进算法能稳定有效地求解层合板最优铺层顺序。

[1]HuHT,LinBH.BucklingOptimizationofSymmetricallyLaminatedPlateswithVariousGeometriesandEndConditions[J].CompositesScienceandTechnology,1995,55(3):277-285.

[2]TopalU,UzmanU.OptimumDesignofLaminatedCompositePlatestoMaximizeBucklingLoadUsingMFDMethod[J].Thin-WalledStructures,2007,45(7-8):660-669.

[3]SoremekunG,GürdalZ,HaftkaRT,etal.CompositeLaminateDesignOptimizationbyGeneticAlgorithmwithGeneralizedElitistSelection[J].ComputersStructures,2001,79(2):131-143.

[4]唐文艳,顾元宪,赵国忠.复合材料层合板铺层顺序优化遗传算法[J].大连理工大学学报,2004,44(2):186-189.

TangWenyan,GuYuanxian,ZhaoGuozhong.Stacking-sequenceOptimizationofCompositeLaminatePlatesbyGeneticAlgorithm[J].JournalofDalianUniversityofTechnology,2004,44(2):186-189.

[5]修英姝,崔德刚.复合材料层合板稳定性的铺层优化设计[J].工程力学,2005,22(6):212-216.

XiuYingshu,CuiDegang.PlyOptimizationDesignforStabilityofCompositeLaminates[J].EngineeringMechanics,2005,22(6):212-216.

[6]KarakayaS,SoykasapÖ.BucklingOptimizationofLaminatedCompositePlatesUsingGeneticAlgorithmandGeneralizedPatternSearchAlgorithm[J].StructuralandMultidisciplinaryOptimization,2009,39(5):477-486.

[7]CorreiaVMF,SoaresCMM,SoaresCAM.BucklingOptimizationofCompositeLaminatedAdaptiveStructures[J].CompositeStructures,2003,62(3/4):315-321.

[8]洪厚全,李玉亮,徐超.复合材料层合板屈曲优化设计的模拟退火算法应用[J].环境与强度,2012,39(6):8-13.

HongHouquan,LiYuliang,XuChao.OptimiaztionofCompositeLaminatesforBucklingbySimulatesAnnealingAlgorithm[J].Structure&EnvironmentEngineering,2012,39(6):8-13.

[9]ErdalO,SonmezOF.OptimalDesignofCompositeLaminatesforMaximumLoadCapacityUsingSimulatedAnnealing[J].CompositeStructures,2005,71:45-52.

[10]BloomfieldMW,HerenciaJE,WeaverPM.EnhancedTwo-levelOptimizationofAnisotropicLaminatedCompositePlateswithStrengthandBucklingConstraints[J].Thin-walledStructures,2009,47(11):1161-1167.

[11]李真,陈秀华,望海.基于蚁群算法的复合材料层合板屈曲优化[J].上海交通大学学报,2012,46(5):768-773.

LiZhen,ChenXiuhua,WangHai.BucklingOptimizationofCompositeLaminateUsingModified

AntColonyOptimization[J].JournalofShanghaiJiaotongUniversity,2012,46(5):758-773.

[12]AliMM,TörnA,ViitanenS.ADirectSearchVariantoftheSimulatedAnnealingAlgorithmforOptimizationInvolvingContinuousVariables[J].Computers&OperationsResearch,2002,29(1):87-102.

[13]KarakayaS,SoykasapÖ.NaturalFrequencyandBucklingOptimizationofLaminatedHybridCompositePlatesUsingGeneticAlgorithmandSimulatedAnnealing[J].StructuralandMultidisciplinaryOptimization,2011,43(1):61-72.

[14]NaritaY,TurveyGJ.MaximizingtheBucklingLoadsofSymmetricallyLaminatedCompositeRectangularPlatesUsingaLayerwiseOptimizationApproach[J].JournalofMechanicalEngineeringScience,2004,218(7):681-691.

[15]RaoARM,ArvindN.AScatterSearchAlgorithmforStackingSequenceOptimisationofLaminateComposites[J].CompositeStructures,2005,70(4):383-402.

[16]IrisarriFX,BassirDH,CarrereN,etal.MultiobjectiveStackingSequenceOptimizationforIaminatedCompositeStructures[J].CompositesScienceandTechnology,2009,69(7/8):983-990.

[17]MarkET.StructuralAnalysisofPolymericCompositeMaterials[M].NewYork:MarcelDekker,Inc.,2004.

[18]KirkpatrickS,GelattCD,VecchiMP.OptimizationbySimulatedAnnealing[J].Science,1983,220(4598):671-680.

(编辑卢湘帆)

Buckling Optimization of Composite Laminates Based on Improved Simulated Annealing Algorithm

Sun ShipingZeng QinglongWu Jianjun

Nanchang Hangkong University,Nanchang,330063

An improved direct search SA algorithm was presented to find globally stacking sequence design of symmetrically composite laminates with maximum buckling load.A new point generation mechanism was introduced into improved SA to enhance the stability and efficiency of the algorithm,which the variational search area for generating new point was proposed to realize global and local search consistent.The design variable was as discrete ply angle and buckling response of laminate was analyzed by Ritz method to consider the bend-twist coupling.A multiple-global optimal problem known was explored to comparative investigate the performance of improved SA.Further,the numerical examples were conducted for different angle increments,layer numbers,aspect ratio and load ratios,the results show that the presented algrothm is capable to optimize the stacking sequence of composite laminates.

buckling optimization;simulated anealing(SA) algorrthm;composite laminates;Ritz method;optimal design

2014-08-07

国家自然科学基金资助项目(11362017);江西省教育厅科学技术研究项目(GJJ13520,GJJ13494)

TB33DOI：10.3969/j.issn.1004-132X.2015.12.020

孙士平,男,1972年生。南昌航空大学航空制造工程学院副教授。主要研究方向为材料/结构拓扑优化设计。发表论文20余篇。曾庆龙，男，1989年生。南昌航空大学航空制造工程学院硕士研究生。吴建军，男，1983年生。南昌航空大学航空制造工程学院硕士研究生。