基于改進模擬退火算法的復合材料層合板屈曲優化

孫士平　曾慶龍　吳建軍

南昌航空大學,南昌,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年生。南昌航空大學航空制造工程學院碩士研究生。