劉　哲，　金達(dá)鋒,　范志瑞

(1.清華大學(xué) 汽車安全與節(jié)能國家重點(diǎn)實(shí)驗(yàn)室，北京100084；2.中北大學(xué) 機(jī)械與動(dòng)力工程學(xué)院，山西 太原030051)

?

基于雙值參數(shù)化方法且考慮制造約束的炭纖維復(fù)合材料鋪層優(yōu)化

劉哲1，金達(dá)鋒1,范志瑞2

(1.清華大學(xué) 汽車安全與節(jié)能國家重點(diǎn)實(shí)驗(yàn)室，北京100084；2.中北大學(xué) 機(jī)械與動(dòng)力工程學(xué)院，山西 太原030051)

基于梯度的優(yōu)化方法對(duì)炭纖維復(fù)合材料層合板的鋪層數(shù)量和順序進(jìn)行優(yōu)化。優(yōu)化問題中以鋪層質(zhì)量為目標(biāo)，并以剛度和制造約束為約束。采用改進(jìn)雙值參數(shù)化方法對(duì)鋪層的材料性能進(jìn)行插值，并基于凸規(guī)劃對(duì)偶算法對(duì)優(yōu)化問題進(jìn)行求解。為了適應(yīng)凸規(guī)劃對(duì)偶算法的特點(diǎn)，將關(guān)于鋪層角度的制造約束表述為少量非線性約束。同時(shí)引入離散度約束和懲罰指數(shù)以消除優(yōu)化結(jié)果中的中間變量。算例結(jié)果驗(yàn)證了該優(yōu)化方法的有效性。

炭纖維復(fù)合材料； 制造約束； 優(yōu)化； 層合板； 凸規(guī)劃對(duì)偶算法

1　前言

由于炭纖維復(fù)合材料較鋼、鋁等材料具有較高的比強(qiáng)度、比剛度以及良好的抗腐蝕性，其在航空、航天等領(lǐng)域得到了廣泛應(yīng)用。炭纖維復(fù)合材料在力學(xué)性能上體現(xiàn)出較強(qiáng)的各向異性使得該種材料具有良好的可設(shè)計(jì)性。多數(shù)情況下為了使結(jié)構(gòu)在滿足一定剛度要求的情況下實(shí)現(xiàn)結(jié)構(gòu)輕量化，對(duì)炭纖維復(fù)合材料的鋪層層數(shù)和鋪層角度的優(yōu)化十分必要。

為了使鋪層設(shè)計(jì)更具有工程意義，在鋪層優(yōu)化時(shí)還需考慮優(yōu)化結(jié)果的可制造性、制造成本以及鋪層在使用過程中隨機(jī)載荷對(duì)鋪層剛度和強(qiáng)度的影響。為此鋪層優(yōu)化的結(jié)果需滿足相應(yīng)的制造約束。在優(yōu)化問題中引入制造約束雖然可能會(huì)影響優(yōu)化結(jié)果的最優(yōu)性，但是可以保證優(yōu)化結(jié)果的可制造性、低成本以及在隨機(jī)載荷下合理的剛度和強(qiáng)度。

Ghiasi[1]對(duì)近年來復(fù)合材料鋪層優(yōu)化方法做了系統(tǒng)性的綜述。其中能夠處理制造約束的優(yōu)化方法主要包括遺傳算法、分級(jí)優(yōu)化方法和多材料拓?fù)鋬?yōu)化方法。

由于遺傳算法在處理整數(shù)型設(shè)計(jì)變量的優(yōu)勢，其在鋪層優(yōu)化中得到了廣泛的關(guān)注[2-4]。研究者通過改進(jìn)基因編碼方式和遺傳操作以及向優(yōu)化問題中引入約束函數(shù)使得遺傳算法具備處理制造約束的能力。但由于遺傳算法的優(yōu)化效率較低，其很難應(yīng)用于鋪層數(shù)量較多且鋪層數(shù)量不確定的優(yōu)化問題中。

Herencia[5]和Liu[6]等利用分級(jí)優(yōu)化方法對(duì)鋪層進(jìn)行優(yōu)化。在第一級(jí)優(yōu)化中利用基于梯度的優(yōu)化算法和鋪層參數(shù)，找到最優(yōu)的鋪層參數(shù)。第二級(jí)優(yōu)化利用遺傳算法找到與最優(yōu)鋪層參數(shù)相對(duì)應(yīng)的并且滿足制造約束的實(shí)際鋪層方案。此種方法在第一級(jí)優(yōu)化中通過利用基于梯度的優(yōu)化算法可以有效減少對(duì)目標(biāo)函數(shù)和約束函數(shù)的計(jì)算次數(shù)并提高優(yōu)化效率。在第二級(jí)優(yōu)化中，雖然遺傳算法的效率較低，但由于鋪層角度與鋪層參數(shù)的函數(shù)關(guān)系簡單，故對(duì)優(yōu)化效率影響不大。但是由于目前研究僅得出某些特定類型鋪層對(duì)應(yīng)的鋪層參數(shù)可行域，此種方法具有較強(qiáng)的局限性。

Lund[7]基于多材料拓?fù)鋬?yōu)化方法提出了離散材料優(yōu)化方法(Discrete material optimization，DMO)對(duì)復(fù)合材料鋪層進(jìn)行優(yōu)化。在此種方法中各鋪層所采用材料的性能參數(shù)剛度是由各備選材料性能參數(shù)剛度的權(quán)重和所確定，各備選材料分別對(duì)應(yīng)各備選鋪層角度。在優(yōu)化中將各備選鋪層角度的權(quán)重值作為設(shè)計(jì)變量。同時(shí)，DMO方法采用基于SIMP法(Solid isotropic material with penalization)或RAMP法(Rational approximation of material properties)的材料插值公式對(duì)取中間值的設(shè)計(jì)變量進(jìn)行懲罰。隨后Bruyneel[8]提出了帶懲罰形函數(shù)方法(Shape functions and parameterization，SFP)，該方法通過借鑒有限元方法中形函數(shù)的構(gòu)造方式對(duì)材料插值公式進(jìn)行改進(jìn)。相對(duì)于DMO方法，SFP方法可有效減少設(shè)計(jì)變量的數(shù)量。但是當(dāng)備選角度較多時(shí)該方法難以找到已有的形函數(shù)并構(gòu)造相應(yīng)的材料插值公式。Gao等[9]對(duì)SFP方法進(jìn)行改進(jìn)并提出雙值參數(shù)化方法(Bi-valued coding parameterization，BCP)。BCP方法不再依賴于形函數(shù)構(gòu)造材料插值公式，故備選角度可為任意數(shù)量，并且可用n個(gè)設(shè)計(jì)變量表示2n個(gè)備選角度。通過DMO、SFP和BCP方法建立的鋪層優(yōu)化問題可以通過基于梯度的優(yōu)化算法進(jìn)行求解，因此此類方法可以處理大規(guī)模的設(shè)計(jì)變量。但是如何將制造約束引入優(yōu)化問題仍存在較大困難。

目前僅S?rensen[10]基于DMO方法就考慮制造約束的鋪層優(yōu)化問題進(jìn)行了研究。其中通過引入大量的線性約束以表示制造約束，并采用序列二次規(guī)劃算法對(duì)優(yōu)化問題進(jìn)行求解。由于在優(yōu)化問題中存在過多約束，優(yōu)化過程中存在設(shè)計(jì)變量的收斂速率較低的問題。

針對(duì)鋪層層數(shù)較多、鋪層層數(shù)可變且考慮制造約束的鋪層優(yōu)化問題，本文通過構(gòu)造數(shù)量較少的非線性約束和對(duì)BCP方法的改進(jìn)，提出了一種考慮制造約束的鋪層優(yōu)化問題的優(yōu)化方法。針對(duì)以鋪層質(zhì)量為優(yōu)化目標(biāo)，結(jié)構(gòu)剛度為約束并且考慮制造約束的鋪層優(yōu)化問題建立數(shù)學(xué)模型并基于凸規(guī)劃對(duì)偶算法進(jìn)行求解。

2　BCP方法介紹[9]

在BCP方法中，各鋪層的材料性能為備選角度材料性能的權(quán)重之和，如式(1)所示。

(1)

式中,下標(biāo)i為鋪層層數(shù)編號(hào)，編號(hào)越大代表鋪層越接近鋪層底部；下標(biāo)j為備選角度編號(hào)；D為材料的彈性矩陣；m為備選角度的數(shù)量；w為各備選角度的權(quán)重。

鋪層內(nèi)各備選角度的權(quán)重值如式(2)所示。

(2)

式中,mv為鋪層角度設(shè)計(jì)變量的數(shù)量；xik為鋪層角度設(shè)計(jì)變量，取值范圍為[-1,1]；ρi為鋪層密度設(shè)計(jì)變量，取值范圍為[0,1]；p、q為懲罰指數(shù)；sjk為一系數(shù)，其表達(dá)式如式(3)所示。

(3)

鋪層角度設(shè)計(jì)變量數(shù)目mv與備選角度數(shù)量m的關(guān)系如式(4)所示。

mv=[log2m]

(4)

當(dāng)備選角度數(shù)量為4時(shí)，根據(jù)式(4)可知鋪層角度設(shè)計(jì)變量數(shù)目為2。根據(jù)式(3)系數(shù)sjk的取值如表1所示。

根據(jù)式(2)和表1中的系數(shù)取值，四種備選角度的權(quán)重分布見圖1。

表 1　系數(shù)sjk取值(mv=2,m=4)

圖 1　備選角度權(quán)重分布圖(mv=2,m=4,ρ=1,p=1)

3　制造約束介紹

為了保證鋪層在隨機(jī)載荷下的剛度、強(qiáng)度并降低制造成本，鋪層設(shè)計(jì)需滿足一定的制造約束。結(jié)合實(shí)際工程經(jīng)驗(yàn)提出優(yōu)化過程中所需滿足的制造約束如下：

(1)為了降低制造成本，備選鋪層角度需從{0°,+45°,-45°,90°}中選取，且各層厚度均相同；

(2)為了保證結(jié)構(gòu)的可制造性，鋪層的需具有連續(xù)性，即最底層鋪層必須存在且兩個(gè)確定存在的鋪層之間不允許出現(xiàn)密度為0的鋪層；

(3)為了在隨機(jī)載荷下有效抑制裂紋在基體內(nèi)的發(fā)展進(jìn)而防止結(jié)構(gòu)的失效，不允許出現(xiàn)連續(xù)n層及以上具有相同鋪層角度的鋪層，n通常取3或4。

圖2和圖3分別列舉出違反制造約束2和制造約束3的鋪層方案。圖2中鋪層4為最底層鋪層，此時(shí)鋪層4的密度為0違反了制造約束2中關(guān)于最底層鋪層必須存在的要求。同時(shí)，在第1和3層鋪層確定存在的情況下，處于二者之間的鋪層2其密度為0。此種鋪層方案違反了制造約束2中關(guān)于兩個(gè)確定存在的鋪層之間不允許出現(xiàn)密度為0鋪層的要求。圖3中四層鋪層均存在，但1至3層的鋪層角度均為同一確定角度，也即鋪層中出現(xiàn)了連續(xù)3層鋪層角度相同的情況。當(dāng)n=3時(shí)，該鋪層方案違反了制造約束3的要求。

4　優(yōu)化問題的建立

優(yōu)化問題以結(jié)構(gòu)的質(zhì)量作為優(yōu)化目標(biāo)，并以剛度和制造約束作為約束。當(dāng)各層的面積、厚度和材料密度均相等時(shí)最小化鋪層質(zhì)量等價(jià)于最小化鋪層層數(shù)。在優(yōu)化過程中優(yōu)化算法需對(duì)鋪層的層數(shù)和各層的鋪層角度進(jìn)行優(yōu)化。

圖 2　違反制造約束2的鋪層示意圖

圖 3　違反制造約束3的鋪層示意圖(n=3)

4.1改進(jìn)BCP方法的設(shè)計(jì)變量

由于在優(yōu)化中引入制造約束，BCP方法中部分設(shè)計(jì)變量并不適應(yīng)此類優(yōu)化問題，故需對(duì)BCP方法中設(shè)計(jì)變量進(jìn)行改進(jìn)。

保留BCP方法中的鋪層角度設(shè)計(jì)變量。制造約束1中要求備選鋪層角度數(shù)量為4，根據(jù)式(4)可知各層鋪層角度設(shè)計(jì)變量數(shù)量為2。

在制造約束2中最底層鋪層必須存在，因此在BCP方法中對(duì)應(yīng)最底層鋪層的密度在優(yōu)化中應(yīng)始終保持為1。如果優(yōu)化過程中各層密度可自由變化則可能會(huì)違反制造約束2，因此需對(duì)各層密度值進(jìn)行約束?？梢酝ㄟ^要求鋪層密度從底部至頂部逐層遞減來滿足制造約束2。將式(5)引入優(yōu)化問題的數(shù)學(xué)模型后，雖然可使優(yōu)化結(jié)果滿足制造約束2，但由于向優(yōu)化數(shù)學(xué)模型中引入過多的線性約束進(jìn)而降低了使得優(yōu)化問題的求解效率。

ρi-ρi+1≤0i=1,......,pn-1

(5)

式中,pn為鋪層的總層數(shù)。

在此不再以BCP方法中各層密度值作為優(yōu)化問題的設(shè)計(jì)變量，而是以相鄰兩層的密度值增量作為密度設(shè)計(jì)變量并記為xi3。為了便于區(qū)分，現(xiàn)將各層的密度值ρi稱為物理密度。各層物理密度值ρi與密度設(shè)計(jì)變量xi3的函數(shù)關(guān)系見式(6)。當(dāng)密度設(shè)計(jì)變量xi3≥0時(shí)，各鋪層物理密度值均滿足式(5)。

(6)

在優(yōu)化過程中基于式(6)求得各鋪層的物理密度后，將各層的物理密度和鋪層角度設(shè)計(jì)變量代入式(2)后即可求得各鋪層內(nèi)所有備選鋪層角度的權(quán)重。隨后根據(jù)各備選鋪層角度的權(quán)重求得各鋪層內(nèi)材料的彈性矩陣并代入有限元模型中求得結(jié)構(gòu)響應(yīng)。

圖4為優(yōu)化問題的設(shè)計(jì)變量以及各設(shè)計(jì)變量與鋪層物理密度和角度的關(guān)系。由圖4可知，除最底層鋪層外每個(gè)鋪層均包含一個(gè)密度設(shè)計(jì)變量和兩個(gè)鋪層角度設(shè)計(jì)變量，最底層鋪層僅包含兩個(gè)鋪層角度設(shè)計(jì)變量。同時(shí)各鋪層的鋪層角度僅與各鋪層相應(yīng)的鋪層角度設(shè)計(jì)變量相關(guān)，但各鋪層的物理密度與其本身以及位于其底部鋪層的密度設(shè)計(jì)變量相關(guān)。

圖 4　鋪層優(yōu)化設(shè)計(jì)變量示意圖

4.2物理密度上限控制參數(shù)T

為了保證物理密度在相鄰層之間的變化較為平緩，引入?yún)?shù)T對(duì)各層物理密度的上限進(jìn)行控制。各鋪層物理密度取值的上限是其底部相鄰鋪層物理密度值的分段函數(shù),由圖5和式(7)所示。參數(shù)T控制著兩段線性函數(shù)的交點(diǎn)。在一定程度上參數(shù)T控制了鋪層總層數(shù)的上限，并可使鋪層質(zhì)量集中于鋪層的底部。式(8)為密度設(shè)計(jì)變量的下限值。

參數(shù)T的初值一般較小，并且在優(yōu)化過程中采用連續(xù)化方法逐步增加參數(shù)T的值。表2給出當(dāng)鋪層層數(shù)為5且鋪層底部兩層的物理密度固定為1和0.95時(shí)，不同的參數(shù)T取值下其余各鋪層的鋪層密度上限值。

由表2可知，隨著參數(shù)T的增加各鋪層的物理密度上限值逐漸增加。當(dāng)T=0.5時(shí)，各鋪層的物理密度上限均相等且為0.95。

圖 5　鋪層物理密度上限示意圖(T=0.1)

(7)

式中：

x13≥ρup(i+1)-ρupi

(8)

表 2　不同T下鋪層物理密度上限值(ρ4=0.950, ρ5=1.000)

4.3最頂層鋪層物理密度非負(fù)約束

由于沒有設(shè)定密度設(shè)計(jì)變量的上限值，故存在密度設(shè)計(jì)變量取值較大并造成鋪層物理密度小于0的情況出現(xiàn)。物理密度小于0的鋪層并無實(shí)際意義，因此需要在優(yōu)化問題中引入如式(9)所示的約束以保證最頂層鋪層的物理密度值非負(fù)。

(9)

4.4目標(biāo)函數(shù)及改進(jìn)BCP方法的備選角度權(quán)重插值公式

在優(yōu)化過程中，凸規(guī)劃對(duì)偶算法不能保證約束不被違反。當(dāng)鋪層的物理密度為負(fù)值時(shí)會(huì)給優(yōu)化帶來問題。首先，式(2)中ρi的q次冪值可能為虛數(shù)并造成該鋪層內(nèi)各備選角度的權(quán)重值無實(shí)際物理意義。其次，如果將鋪層的質(zhì)量用各鋪層的物理密度之和表示，則鋪層的物理密度取負(fù)值更有利于降低目標(biāo)函數(shù)的值。此種情況下不利于解的收斂，甚至使優(yōu)化問題無法找到可行解。

為此需對(duì)目標(biāo)函數(shù)和BCP方法中各備選角度權(quán)重插值公式進(jìn)行改進(jìn)并實(shí)現(xiàn)對(duì)中間密度有效的懲罰。目標(biāo)函數(shù)可表示為各鋪層物理密度絕對(duì)值之和。但由于絕對(duì)值函數(shù)連續(xù)不可導(dǎo)，則無法利用凸規(guī)劃對(duì)偶算法對(duì)優(yōu)化問題進(jìn)行求解。因此可以通過構(gòu)造一階導(dǎo)數(shù)連續(xù)的函數(shù)來近似絕對(duì)值函數(shù)，如式(10)所示。

(10)

式中,ε為極小的正數(shù)。

同時(shí)，在權(quán)重插值公式內(nèi)取消懲罰指數(shù)q，改進(jìn)后的各備選角度權(quán)重插值公式如式(11)所示。

(11)

為了實(shí)現(xiàn)對(duì)中間密度的懲罰，將懲罰指數(shù)q引入鋪層質(zhì)量函數(shù)。各層質(zhì)量相加即得到鋪層的總質(zhì)量如式(12)所示。

(12)

4.5制造約束3

在制造約束3中，不允許出現(xiàn)連續(xù)n層及以上具有相同鋪層角度的鋪層。該制造約束可表示為：

(13)

(14)

該函數(shù)是符號(hào)函數(shù)的連續(xù)化近似。因此可根據(jù)近似符號(hào)函數(shù)的輸出對(duì)某一項(xiàng)約束是否被違反進(jìn)行判斷。為了使式(13)中約束滿足時(shí)近似符號(hào)函數(shù)的輸出為-1，將式(13)中的約束松弛為如式(15)所示。

i=1,2,...,pn-n-1

(15)

隨后根據(jù)式(15)中各約束對(duì)應(yīng)的近似符號(hào)函數(shù)的輸出之和對(duì)鋪層是否違反制造約束3進(jìn)行判斷。當(dāng)且僅當(dāng)式(15)中各約束均不被違反時(shí)，鋪層方案才不違反制造約束3。因此各約束對(duì)應(yīng)的近似符號(hào)函數(shù)的輸出之和只有小于等于(n+1-pn)時(shí)，才能保證制造約束3不被違反。針對(duì)制造約束3需引入優(yōu)化問題的約束如式(16)所示。

(16)

引入過濾函數(shù)分別對(duì)鋪層的物理密度以及鋪層角度設(shè)計(jì)變量進(jìn)行過濾，使過濾后的變量更接近整數(shù)。隨后用過濾后的鋪層物理密度以及角度設(shè)計(jì)變量計(jì)算各鋪層內(nèi)各角度所對(duì)應(yīng)的權(quán)重值，如式(17)所示。

(17)

基于Heaviside函數(shù)構(gòu)造鋪層物理密度的過濾函數(shù)，如式(18)所示。

(18)

式中，α1為過濾函數(shù)的閾值。

圖 6　物理密度過濾函數(shù)在不同α1下的函數(shù)圖像(β=10)

如圖6所示物理密度過濾函數(shù)連續(xù)，且當(dāng)β值較大時(shí)大部分物理密度在過濾后取值均為0或1。過濾函數(shù)的閾值α1控制著函數(shù)圖像沿橫軸的平移量。

過濾后物理密度對(duì)過濾前物理密度ρi的一階導(dǎo)數(shù)如式(19)所示。

(19)

由式(19)可知，在ρi=α1處關(guān)于ρi的左導(dǎo)數(shù)等于右導(dǎo)數(shù)，則物理密度過濾函數(shù)一階導(dǎo)數(shù)連續(xù)。同理可構(gòu)造鋪層角度設(shè)計(jì)變量過濾函數(shù)。如式(20)所示。

(20)

式中，α2為過濾函數(shù)的閾值。

圖 7　鋪層角度設(shè)計(jì)變量過濾函數(shù)在不同α2下的函數(shù)圖像(β=10)

如圖7所示鋪層角度設(shè)計(jì)變量過濾函數(shù)連續(xù)，且當(dāng)β值較大時(shí)大部分變量在過濾后取值均為-1或1。過濾函數(shù)的閾值α2控制著函數(shù)圖像沿橫軸的平移量。

過濾后角度設(shè)計(jì)變量對(duì)過濾前角度設(shè)計(jì)變量xik的一階導(dǎo)數(shù)如式(21)所示。

(21)

由式(21)可知，在xik=α2處關(guān)于xik的左導(dǎo)數(shù)等于右導(dǎo)數(shù)，則鋪層角度設(shè)計(jì)變量過濾函數(shù)一階導(dǎo)數(shù)連續(xù)。

4.6鋪層角度設(shè)計(jì)變量離散度約束

在實(shí)際優(yōu)化中僅依靠懲罰指數(shù)p很難保證鋪層角度設(shè)計(jì)變量最終收斂至1或-1。因此有必要在優(yōu)化問題中引入相關(guān)約束以消除鋪層角度設(shè)計(jì)變量中的中間變量。為了衡量鋪層角度設(shè)計(jì)變量偏離1和-1的程度，引入離散度概念并利用其構(gòu)造相關(guān)約束。對(duì)于某一鋪層的角度設(shè)計(jì)變量而言，其離散度的計(jì)算公式如式(22)所示。

(22)

式中,Mndi為第i層鋪層角度設(shè)計(jì)變量的離散度；Mndik為第i層鋪層第k個(gè)角度設(shè)計(jì)變量的離散度。

當(dāng)xik=0時(shí)，Mndik為1且為最大值；當(dāng)xik=-1或xik=1時(shí)，Mndik為0且為最小值。因此Mndi的最大值為且最小值為0。根據(jù)離散度函數(shù)的特點(diǎn)，對(duì)離散度函數(shù)的上限進(jìn)行約束，如式(23)所示。

(23)

式中,Mndup為鋪層角度設(shè)計(jì)變量離散度上限值。

在實(shí)際優(yōu)化過程中，為防止優(yōu)化結(jié)果過早地收斂至局部最優(yōu)解，需采用連續(xù)化方法逐漸減少鋪層角度設(shè)計(jì)變量離散度上限Mndup的值。

4.7優(yōu)化數(shù)學(xué)模型

對(duì)于鋪層優(yōu)化問題，優(yōu)化變量為密度設(shè)計(jì)變量和鋪層角度設(shè)計(jì)變量，優(yōu)化目標(biāo)為鋪層質(zhì)量如式(12)所示，其約束包括結(jié)構(gòu)柔順度、式(9)、式(16)和式(23)。

find{xik}(i=1,...,pn;k=1,...3)

s.t.F=Ku

C-Cmax≤0

(24)

-1≤xi1≤1-1≤xi2≤1xi3≥ρup(i+1)-ρupi

式中,F為節(jié)點(diǎn)力向量；K為結(jié)構(gòu)剛度矩陣；u為節(jié)點(diǎn)位移向量；C為結(jié)構(gòu)柔順度；Cmax為結(jié)構(gòu)柔順度上限值。

5　優(yōu)化流程

本文優(yōu)化問題的求解基于凸規(guī)劃對(duì)偶算法，該算法可對(duì)多變量、少約束、高度非線性的優(yōu)化問題進(jìn)行高效率地求解。為了防止優(yōu)化問題的解過早地收斂于局部最優(yōu),在優(yōu)化過程中采用連續(xù)化方法對(duì)參數(shù)q、p、T以及Mndup的值進(jìn)行改變。

在迭代初期鋪層總層數(shù)的范圍以及最優(yōu)鋪層角度均不確定，此時(shí)如果對(duì)鋪層物理密度以及鋪層角度設(shè)計(jì)變量的離散度提出較高的要求會(huì)造成過早收斂于局部最優(yōu)解或無法收斂的情況。因?yàn)閰?shù)T控制各層物理密度的上限，所以在此階段應(yīng)首先確定參數(shù)T的大致范圍。在大致確定了鋪層的總層數(shù)并使得結(jié)構(gòu)剛度以及各制造約束得到滿足后，逐漸提高對(duì)鋪層物理密度以及鋪層角設(shè)計(jì)變量離散度的要求，并直至參數(shù)q、p以及Mndup達(dá)到預(yù)定的上限或下限為止。需要注意的是，隨著對(duì)離散度要求的提高，結(jié)構(gòu)剛度通常會(huì)相應(yīng)降低。當(dāng)結(jié)構(gòu)剛度下降明顯時(shí)，鋪層總層數(shù)也需相應(yīng)增加，因此參數(shù)T的值也需要繼續(xù)隨著結(jié)構(gòu)剛度的變化而改變。

(25)

綜上所述，優(yōu)化流程如圖8所示。

圖 8　優(yōu)化流程示意圖

6　數(shù)值算例

算例對(duì)兩種不同問題進(jìn)行優(yōu)化，各算例的載荷、邊界條件以及幾何尺寸分別見圖9和圖10。各數(shù)值算例中所采用的材料均相同，材料性能如表3所示。

所有算例的有限元分析均基于有限元軟件ABAQUS完成，在有限元模型中選取四節(jié)點(diǎn)減縮積分殼單元S4R。單個(gè)鋪層層厚為0.125 mm。優(yōu)化流程的剩余步驟均基于MATLAB完成。

圖 9　算例1幾何尺寸、載荷及邊界條件示意圖

圖 10　算例2幾何尺寸、載荷及邊界條件示意圖

PropertiesValuesLongitudinaltensilemodulus162GPaLongitudinalcompressivemodulus145GPaTransversetensilemodulus9．2GPaTransversecompressivemodulus9．5GPaIn?planeshearmodulus5．0GPaPoisson＇sratio0．3

所有算例中參數(shù)q、p、Mndup以及T的初值均相同且取值為：q=1，p=1，Mndup=1，T=0.1。式(18)和式(20)中過濾函數(shù)的過濾閾值α1和α2分別為0.5和0，極小正數(shù)ε取0.001，β取5。n取值為3，備選角度為{0°,45°,90°,-45°}，備選角度與鋪層角度設(shè)計(jì)變量的對(duì)應(yīng)關(guān)系見圖11。

鋪層角度設(shè)計(jì)變量初始值均取0，也即不偏重于任何一個(gè)備選角度。鋪層密度設(shè)計(jì)變量的初值取參數(shù)T初始值對(duì)應(yīng)的鋪層密度設(shè)計(jì)變量的上限值。

圖 11　備選角度編碼示意圖

算例1中，柔順度的上限Cmax分別取1×10-3J和6×10-4J，分別針對(duì)兩種柔順度上限進(jìn)行優(yōu)化，其優(yōu)化結(jié)果見表4。

表 4　算例1優(yōu)化結(jié)果

迭代終止時(shí)所有約束均滿足，且優(yōu)化結(jié)果的結(jié)構(gòu)柔順度值均等于相應(yīng)柔順度上限值。相關(guān)優(yōu)化參數(shù)的取值見表5。

表 5　算例1中優(yōu)化參數(shù)最終取值

算例2中，柔順度的上限Cmax分別取5×10-3J和2×10-3J，分別針對(duì)兩種柔順度上限進(jìn)行優(yōu)化，其優(yōu)化結(jié)果如表6所示。

表 6　算例2優(yōu)化結(jié)果

迭代終止時(shí)所有約束均滿足，且優(yōu)化結(jié)果的結(jié)構(gòu)柔順度值均等于相應(yīng)柔順度上限值。相關(guān)優(yōu)化參數(shù)的取值如表7所示。

表 7　算例2中優(yōu)化參數(shù)最終取值

根據(jù)各算例的優(yōu)化結(jié)果可知，大部分鋪層的物理密度值為0或1，也即明確了鋪層的有無。雖然目標(biāo)函數(shù)中的懲罰指數(shù)q對(duì)中間密度有懲罰作用，但是對(duì)于某一柔順度上限值其通常對(duì)應(yīng)的理論最優(yōu)鋪層層數(shù)并非整數(shù)，因此在優(yōu)化結(jié)果中通常會(huì)出現(xiàn)在鋪層物理密度值0和1的過渡區(qū)域存在物理密度為中間值的鋪層。

由于將相鄰兩層物理密度之差作為優(yōu)化設(shè)計(jì)變量，因此所有算例中鋪層的物理密度從底部到頂部呈現(xiàn)遞減趨勢且滿足式(5)。優(yōu)化結(jié)果中避免出現(xiàn)夾空層的現(xiàn)象出現(xiàn)。同時(shí)，最底層鋪層的物理密度在優(yōu)化過程中為定值1進(jìn)而保證了最底部鋪層必須存在的要求。

除算例1中針對(duì)Cmax=6×10-4J的優(yōu)化結(jié)果外，其余算例中鋪層角度設(shè)計(jì)變量離散度上限值Mndup均減少至0.02。此種情況下除鋪層物理密度為0的鋪層外，其余鋪層對(duì)應(yīng)的鋪層角度設(shè)計(jì)變量均接近-1或1，即可明確地得到相應(yīng)鋪層的鋪層角度。

對(duì)于算例1，文獻(xiàn)[10]在不考慮制造約束的前提下進(jìn)行了鋪層優(yōu)化，最終得到各層的最優(yōu)鋪層角度均為0°。通過枚舉法驗(yàn)證了解的最優(yōu)性。根據(jù)表4所示各層的鋪層角度設(shè)計(jì)變量值可知，優(yōu)化結(jié)果中大部分鋪層的鋪層角度為0°，其中還存在少量鋪層角度為90°的鋪層。90°鋪層存在的主要目的在于防止連續(xù)3層或3層以上0°鋪層的出現(xiàn)。對(duì)于算例2，由表6可知最優(yōu)鋪層角度以90°為主，優(yōu)化結(jié)果中存在的0°鋪層主要目的是防止連續(xù)3層或3層以上90°鋪層的出現(xiàn)。因此，該優(yōu)化方法在保證鋪層方案滿足制造約束的同時(shí)有效地優(yōu)化了鋪層角度。但需要注意的是解的最優(yōu)性并不能得到保證。

7　結(jié)論

以相鄰層的密度增量作為鋪層密度設(shè)計(jì)變量可有效避免優(yōu)化結(jié)果中出現(xiàn)夾空層的現(xiàn)象；控制參數(shù)T和目標(biāo)函數(shù)的構(gòu)造保證物理密度在相鄰層之間的變化較為平緩且有效消除了中間物理密度；引入關(guān)于鋪層角度的非線性制造約束和鋪層角度設(shè)計(jì)變量離散度約束保證了優(yōu)化結(jié)果中出現(xiàn)明確且滿足制造約束的鋪層角度。

[1]Ghiasi H, Pasini D, Lessard L. Optimum stacking sequence design of composite materials Part I: Constant stiffness design[J]. Composite Structures, 2009, 90(1): 1-11.

[2]Liu B, Haftka R T, Akgün M A, et al. Permutation genetic algorithm for stacking sequence design of composite laminates[J]. Computer Methods in Applied Mechanics and Engineering, 2000, 186(2): 357-372.

[3]Hansel W, Treptow A, Becker W, et al. A heuristic and a genetic topology optimization algorithm for weight-minimal laminate structures[J]. Composite Structures, 2002, 58(2): 287-294.

[4]Toropov V V, Jones R, Willment T, et al. Weight and manufacturability optimization of composite aircraft components based on a genetic algorithm[P]. 6th World Congresses of Structural and Multidisciplinary Optimization, Rio de Janeiro, Brazil, 2005, 30.

[5]Herencia J E, Weaver P M, Friswell M I. Optimization of anisotropic composite panels with T-shaped stiffeners including transverse shear effects and out-of-plane loading[J]. Structural and Multidisciplinary Optimization, 2008, 37(2): 165-184.

[6]Liu D, Toroporov V V, Querin O M, et al. Bilevel optimization of blended composite wing panels[J]. Journal of Aircraft, 2011, 48(1): 107-118.

[7]Lund E, Stegmann J. On structural optimization of composite shell structures using a discrete constitutive parametrization[J]. Wind Energy, 2005, 8(1): 109-124.

[8]Bruyneel M. SFP a new parameterization based on shape functions for optimal material selection: Application to conventional composite plies[J]. Structural and Multidisciplinary Optimization, 2011, 43(1): 17-27.

[9]Gao T, Zhang W, Duysinx P. A bi-value coding parameterization scheme for the discrete optimal orientation design of the composite laminate[J]. International Journal for Numerical Methods in Engineering, 2012, 91(1): 98-114.

[10]S?rensen S N, Lund E. Topology and thickness optimization of laminated composites including manufacturing constraints[J]. Structural and Multidisciplinary Optimization, 2013, 48(2): 249-265.

Ply optimization of carbon fiber-reinforced plastics laminates based on a gradient-based optimization method with manufacturing constraints

LIU Zhe1,JIN Da-feng1,FAN Zhi-rui2

(1.StateKeyLaboratoryofAutomotiveSafetyandEnergy,TsinghuaUniversity,Beijing100084,China;2.SchoolofMechanicalandPowerEngineering,NorthUniversityofChina,Taiyuan030051,China)

We optimized the number of plies and the stacking sequence of carbon fiber-reinforced plastics laminates using a gradient-based optimization method. To attain industrial relevance and reduce the risk of failure and the cost, certain manufacturing constraints were considered in the optimization. Manufacturing constraints were: 1) candidate fiber orientations were limited to a finite set, 2) the bottom layer must be solid; 3) an intermediate void was not allowed, and 4) two or three identical contiguous layers were allowed at most. In the optimization, the mass was minimized with stiffness and manufacturing constraints. In order to avoid introducing linear constraints to the optimization model, a material property interpolation scheme called a bi-value coding parameterization scheme was modified and the manufacturing constraints concerning fiber orientation were expressed as nonlinear constraints. Meanwhile, a dispersion constraint and penalty index were introduced to eliminate intermediate variables. The optimization problem was solved by a convex programming dual algorithm. Numerical examples validate the optimization method.

Carbon fiber composite; Manufacturing constraint; Optimization; Laminates; Convex programming dual algorithm

JIN Da-feng, Ph. D, Associate Professor. E-mail：Jindf@tsinghua.edu.cn

introduction: LIU Zhe, Ph. D Candidate. E-mail：lz880818@hotmail.com

1007-8827(2016)01-0068-09

TB330.1

A

2015-11-29；

2016-01-14

金達(dá)鋒，博士，副教授. E-mail：Jindf@tsinghua.edu.cn

劉哲,博士生. E-mail：lz880818@hotmail.com