基于連續(xù)階躍參考應(yīng)力的SKO結(jié)構(gòu)拓?fù)鋬?yōu)化方法

侯劍云 丁曉紅

上海理工大學(xué)，上海，200093

0 引言

結(jié)構(gòu)拓?fù)鋬?yōu)化是指尋求材料分布的最優(yōu)拓?fù)湫螒B(tài)，其目的是在設(shè)計域空間內(nèi)尋求結(jié)構(gòu)最佳的傳力路線形式，以優(yōu)化結(jié)構(gòu)的某些性能或減輕結(jié)構(gòu)的重量。結(jié)構(gòu)拓?fù)鋬?yōu)化已發(fā)展100多年，大體分為解析方法和數(shù)值方法［1］。解析方法的理論數(shù)學(xué)求解困難，不便于在工程實(shí)際中直接應(yīng)用。近幾十年來，計算機(jī)技術(shù)在結(jié)構(gòu)分析中的普遍應(yīng)用使結(jié)構(gòu)優(yōu)化數(shù)值方法得到了迅猛發(fā)展，它主要分為兩大類：一類是可解決各類結(jié)構(gòu)的尺寸、形狀及拓?fù)鋬?yōu)化，但計算效率和通用性并不理想的均勻化方法（homogenization method）［2］；另一類是以漸進(jìn)結(jié)構(gòu)優(yōu)化方法（evolutionary structural optimization）［3－4］為代表的優(yōu)化算法，它的特點(diǎn)是在宏觀的角度上對結(jié)構(gòu)材料分布進(jìn)行設(shè)計，計算效率高、通用性好，具有啟發(fā)式特征，能得到近似最優(yōu)解。在工程實(shí)際中，近似最優(yōu)解通常被廣泛采納，因此這一類算法具有良好的發(fā)展前景。

SKO（soft kill option）方法是一種啟發(fā)式的算法，其基本原理是逐漸“軟化”低應(yīng)力的材料，“硬化”高應(yīng)力的材料，使經(jīng)過優(yōu)化后的結(jié)構(gòu)應(yīng)力水平變得更均勻。它最初是德國Karlsruhe研究中心提出的［5－6］。目前，國內(nèi)外對SKO方法的研究大體上分為兩種，一是采用以體積率作為刪除準(zhǔn)則的方法［7］，二是以參考應(yīng)力作為材料刪除準(zhǔn)則的方法［8］。前者通過多次刪除一定的體積使優(yōu)化結(jié)果收斂于目標(biāo)體積，每次刪除過程中將體積率對應(yīng)的應(yīng)力值作為刪除標(biāo)準(zhǔn)逐漸軟化和硬化材料。不同的目標(biāo)體積可得到不同的拓?fù)浣Y(jié)構(gòu)，但體積率對應(yīng)的應(yīng)力值求取困難，運(yùn)算耗時長。后者通過結(jié)構(gòu)的應(yīng)力分布選擇參考應(yīng)力逐漸軟化和硬化材料，這種方法參考應(yīng)力選取簡單，運(yùn)算耗時短，但不同的參考應(yīng)力對優(yōu)化結(jié)果影響較大。若參考應(yīng)力取值過大會導(dǎo)致過刪除，過小會導(dǎo)致優(yōu)化的前后體積改變不大。丁曉紅等［8］提出了以平均應(yīng)力、過濾應(yīng)力等為參考應(yīng)力的刪除準(zhǔn)則。在迭代過程中，參考應(yīng)力與結(jié)構(gòu)的應(yīng)力分布相關(guān)聯(lián)，由于參考應(yīng)力隨著迭代次數(shù)的增加而增大，不同設(shè)計問題的變化規(guī)律是不同的，材料刪除的規(guī)律具有不確定性，有可能存在大片區(qū)域的應(yīng)力高于參考應(yīng)力水平，材料無法繼續(xù)刪除，導(dǎo)致迭代過早收斂，僅達(dá)到局部最優(yōu)。

局部最優(yōu)是一種尋優(yōu)過程中的早熟現(xiàn)象，需要避免。為此，本文提出一種基于連續(xù)階躍參考應(yīng)力的尋優(yōu)迭代策略，可有效地改善早熟現(xiàn)象，得到最優(yōu)結(jié)構(gòu)。

1 連續(xù)階躍參考應(yīng)力SKO方法介紹

1.1 基本原理介紹

假定結(jié)構(gòu)由不同的材料組成，將每個單元的彈性模量作為參數(shù)來改變結(jié)構(gòu)的拓?fù)湫螒B(tài)。材料的彈性模量被定義為溫度的函數(shù)，即隨著溫度的升高，彈性模量變小，材料被“軟化”，當(dāng)材料“軟”到一定程度，可認(rèn)為材料被刪除；同時隨著溫度的降低，材料的彈性模量變大，材料被“硬化”。溫度T和彈性模量E間的關(guān)系可假設(shè)為線性關(guān)系，如圖1所示，此處的溫度沒有物理意義，僅是單元彈性模量改變的控制器。在應(yīng)用連續(xù)階躍參考應(yīng)力的SKO方法的過程中，會經(jīng)過多階的局部尋優(yōu)過程，每階尋優(yōu)都經(jīng)過數(shù)次迭代。溫度的迭代公式為

式中，上標(biāo)i、j表示第i階尋優(yōu)過程中第j次迭代；Tn（i，j）為第i階尋優(yōu)過程中第j次迭代時n節(jié)點(diǎn)的溫度；σn（i，j－1）為第i階尋優(yōu)過程中第j－1次迭代時n節(jié)點(diǎn)的等效節(jié)點(diǎn)應(yīng)力；σ（i，j）為第i階局部尋優(yōu)過程中第j次迭代的連續(xù)階躍參考r應(yīng)力；s為步長因子；Tmax、Tmin為迭代時的溫度上下限，上限為100，下限為0。

圖1 溫度與彈性模量關(guān)系圖

由式（1）可知，如果節(jié)點(diǎn)n的應(yīng)力大于參考應(yīng)力，則節(jié)點(diǎn)溫度降低，材料的彈性模量增大，材料被“硬化”；否則節(jié)點(diǎn)溫度升高，材料被“軟化”。Emax為所選用材料的彈性模量，用以模擬固體區(qū)域 。 當(dāng) Tn（i，j－1）在［0，100］之間時，彈性模量和溫度之間的關(guān)系滿足如圖1所示的線性關(guān)系。如果Tn（i ，j－1）≤0，則E＝Emax；當(dāng) Tn（i，j－1）≥100時，則E＝Emin。Emin非常小，接近結(jié)構(gòu)空隙率的數(shù)值，經(jīng)驗證明，取Emin＝Emax／1000即可以確保數(shù)值穩(wěn)定［7］。

1.2 局部尋優(yōu)過程和連續(xù)階躍參考應(yīng)力

在SKO方法的尋優(yōu)過程中，有可能出現(xiàn)局部最優(yōu)結(jié)構(gòu)，它存在著一些材料堆積區(qū)域，此時若提高參考應(yīng)力，則尋優(yōu)過程又可繼續(xù)進(jìn)行。本文提出一種連續(xù)階躍參考應(yīng)力的概念，使結(jié)構(gòu)能在尋優(yōu)過程中自動調(diào)節(jié)參考應(yīng)力的取值，從而逐步使結(jié)構(gòu)趨優(yōu)。其計算公式為

其中，ζ（i，j）為第i階尋優(yōu)過程中第j次迭代的權(quán)重系數(shù)；連續(xù)階躍參考應(yīng)力σr（i，j）是以第i階尋優(yōu)過程中初次迭代的平均應(yīng)力σa（i）為基數(shù)乘以一個權(quán)重系數(shù)ζ（i，j）得 到的。 權(quán) 重 系 數(shù)ζ（i，j）計 算公 式如下：

式 中 ，ζ0（i）為第i階尋優(yōu)權(quán)重系數(shù)初值；α為常數(shù)，取0.1；Δζ為權(quán)重系數(shù)增量。

第i階權(quán)重系數(shù)初值ζ0（i）為第i－1階尋優(yōu)過程中最后一次迭代的權(quán)重系數(shù)，即

由以上公式可知，連續(xù)階躍參考應(yīng)力的變化趨勢與權(quán)重系數(shù)基本一致。權(quán)重系數(shù)ζ（i，j）隨著尋優(yōu)階數(shù)的增加而逐階增大。以ζ0（0）＝0.85、α＝0.1、Δζ＝0.05為參數(shù)，分別經(jīng)過4階局部尋優(yōu)過程，來說明權(quán)重系數(shù)的變化，由于每階迭代次數(shù)不確定，為便于說明，假定每階尋優(yōu)的迭代次數(shù)均為10。如圖2所示，曲線劃分為5個區(qū)域。j在進(jìn)入下一階局部尋優(yōu)過程時重新計數(shù)，i在原來的基礎(chǔ)上增加一次，每階尋優(yōu)之初，Δζαj與Δζ剛好相等，式（4）、式（5）式變?yōu)棣疲╥，0）＝ζ0（i），ζ0（i）＝，故局部尋優(yōu)過程的第一步迭代與上一階局部尋優(yōu)過程的最后一步迭代的權(quán)重系數(shù)相等，函數(shù)不發(fā)生突變。而隨著迭代次數(shù)的增加，Δζαj很快衰減為零，最終權(quán)重系數(shù)ζ（i，j）等于第i階權(quán)重系數(shù)初值ζ（i）與權(quán)重系數(shù)增量Δζ之和。區(qū)域1、3、4、5的權(quán)

0重系數(shù)均滿足以上所述的提升方式，而區(qū)域2經(jīng)過的兩次迭代都達(dá)到局部最優(yōu)，說明權(quán)重系數(shù)增大得不夠，需再次提升，才能進(jìn)入下一階尋優(yōu)過程，這是權(quán)重系數(shù)的自動調(diào)節(jié)過程。

圖2 權(quán)重系數(shù)ζ變化圖

1.3 局部最優(yōu)判斷標(biāo)準(zhǔn)和收斂條件

在尋優(yōu)過程中，當(dāng)結(jié)構(gòu)前后兩次迭代的體積變化小于給定的局部體積容差ε1時，認(rèn)為結(jié)構(gòu)達(dá)到一個局部最優(yōu)狀態(tài)；當(dāng)小于給定的總體體積容差ε時，認(rèn)為得到最終的優(yōu)化結(jié)果。顯然，ε1要大于ε。實(shí)驗證明，ε1取100～1000倍的ε時可確保數(shù)值穩(wěn)定。這里的體積指的是量綱一化的名義體積。SKO方法中結(jié)構(gòu)的設(shè)計變量是彈性模量E，定義第i階尋優(yōu)過程中第j次迭代的名義體積為

式 中 ，En（i，j）為第i階尋優(yōu)過程中第j次迭代時n單元的彈性模量；hn和An分別為n單元的厚度和面積。

將名義體積量綱一化得到第i階尋優(yōu)過程中第j次迭代時量綱一體積指標(biāo)，其公式為

2 算法流程

本文的算法流程如圖3所示。具體步驟如下：

（1）建立初始物理模型，并確定設(shè)計空間，實(shí)際結(jié)構(gòu)中某些具有特殊要求的邊界或面需保留而不作為設(shè)計區(qū)域。

（2）賦予結(jié)構(gòu)具有如圖1所示的材料屬性，并用有限元離散網(wǎng)格施加約束和載荷，圖3中的T0為環(huán)境 溫 度，和 T（i，j）n一樣，T0沒有實(shí)際的物理意義。

（3）進(jìn)行結(jié)構(gòu)線性靜力有限元分析，提取節(jié)點(diǎn)應(yīng)力，通過式（3）～式（5）計算連續(xù)階躍參考應(yīng)力，權(quán)重系數(shù)ζ（i，j）不得超過規(guī)定的權(quán)重系數(shù)上限ζmax。

圖3 算法流程

（4）按式（1）、式（2）實(shí)現(xiàn)“軟硬化”單元的操作。

（5）未達(dá)到局部最優(yōu)時，更新j，重復(fù)步驟（3），當(dāng)達(dá)到局部最優(yōu)結(jié)構(gòu)后，更新i、j，再重復(fù)步驟（3），直到滿足上文所述的迭代終止條件，終止迭代。

3 設(shè)計參數(shù)選擇與算例

在迭代過程中，不同的參考應(yīng)力對優(yōu)化結(jié)果的影響較大。本文的連續(xù)階躍參考應(yīng)力因權(quán)重系數(shù)而階躍式地由小到大變化。影響權(quán)重系數(shù)的參數(shù)有初始權(quán)重系數(shù)ζ（0）0，權(quán)重系數(shù)增量Δζ和權(quán)重系數(shù)上限ζmax。初始權(quán)重系數(shù)ζ（0）0取得過大，在開始時便使材料刪除得過多，由于結(jié)構(gòu)的后續(xù)分析基于刪除材料后的結(jié)構(gòu)，會造成某些需“硬化”的材料無法恢復(fù)，從而最終不能達(dá)到最優(yōu)拓?fù)浣Y(jié)構(gòu)。權(quán)重系數(shù)增量Δζ和權(quán)重系數(shù)上限ζmax，它們分別控制參考應(yīng)力的增長幅值和上限。本文通過受兩面內(nèi)力懸臂梁的拓?fù)鋬?yōu)化例子探討不同的權(quán)重系數(shù)參數(shù)對優(yōu)化結(jié)果的影響。

3.1 受兩面內(nèi)力懸臂梁的拓?fù)鋬?yōu)化

如圖4a所示的初始物理模型，懸臂梁的長寬比為2，板的一端固定，受兩垂直向下的力P。從圖4b所示的初始模型應(yīng)力分布圖可以看出，固定端的兩角和受力點(diǎn)存在應(yīng)力集中。在不考慮應(yīng)力集中區(qū)域的情況下，結(jié)構(gòu)的應(yīng)力值集中在［0.02，51.2］MPa之間。

圖4 長方形懸臂梁受兩面內(nèi)力模型

權(quán)重系數(shù)參數(shù)分別按表1所示的4種方案進(jìn)行選取，局部體積容差ε1和總體體積容差ε如表1所示。

表1 參數(shù)取值表

方案1中初始權(quán)重系數(shù)ζ0（0）與權(quán)重系數(shù)上限ζmax相等，即采用無連續(xù)階躍參考應(yīng)力的SKO方法，始終以平均應(yīng)力的0.8倍作為參考應(yīng)力進(jìn)行拓?fù)鋬?yōu)化。優(yōu)化后材料分布如圖5a所示。結(jié)構(gòu)的中間左上區(qū)域材料堆積，對比應(yīng)力分布圖（圖5b），該區(qū)域的應(yīng)力接近白色，說明應(yīng)力水平不高，但卻大于參考應(yīng)力而未被軟化。若不考慮應(yīng)力集中區(qū)域，存在材料的區(qū)域應(yīng)力值集中在［0，51.6］MPa之間，相比于初始模型的應(yīng)力水平均勻程度無明顯改善，優(yōu)化結(jié)果僅得到局部最優(yōu)。圖5c為量綱一體積指標(biāo)V、權(quán)重系數(shù)ζ和量綱一化后的參考應(yīng)力指標(biāo)σr的迭代進(jìn)程圖。權(quán)重系數(shù)和參考應(yīng)力不變，結(jié)構(gòu)的體積隨迭代次數(shù)的增加而降低，直到達(dá)到總體體積容差而終止，此次歷經(jīng)了221次迭代，最終體積為初始體積的0.61。

圖5 方案1的結(jié)果圖

方案2采用了連續(xù)階躍參考應(yīng)力的方案，初始權(quán) 重 系 數(shù) 為ζ0（0）為0.5，權(quán)重系數(shù)增量Δζ為0.05，權(quán)重系數(shù)上限ζmax為1.25。優(yōu)化后材料分布如圖6a所示，與圖5a相比有較多的細(xì)部結(jié)構(gòu)，其應(yīng)力云圖6b在不考慮應(yīng)力集中區(qū)域的情況下，存在材料的區(qū)域應(yīng)力值集中在［17.2，51.5］MPa之間，相比于方案1的應(yīng)力水平均勻程度有所提高，得到較為理想的近似最優(yōu)解。圖6c為性能指標(biāo)的迭代進(jìn)程圖。權(quán)重系數(shù)連續(xù)階躍增大最終達(dá)到權(quán)重系數(shù)上限ζmax。參考應(yīng)力變化與權(quán)重系數(shù)基本一致。體積隨著迭代次數(shù)的增加一直減小，最終趨于總體體積容差而終止。與方案1相比，此次歷經(jīng)了325次迭代，結(jié)構(gòu)的最終體積變?yōu)樵瓉淼?.41，比方案1多刪除一部分。

方案3與方案2相比，權(quán)重系數(shù)增量Δζ從0.05提升為0.12。優(yōu)化后的材料分布如圖7a所示。方案3所得到的優(yōu)化結(jié)果與方案2的優(yōu)化結(jié)果差別不大，其應(yīng)力分布圖（圖7b）與方案2的應(yīng)力分布差別不大，存在材料區(qū)域的應(yīng)力值集中在［17.2，51.4］MPa之間，相比于方案2的應(yīng)力水平均勻程度相當(dāng)，同樣得到較為理想的近似最優(yōu)解。由圖7c可知，權(quán)重系數(shù)的逐階增幅相應(yīng)比方案2高，直到趨于權(quán)重系數(shù)上限ζmax。參考應(yīng)力增大趨勢與權(quán)重系數(shù)基本一致。體積指標(biāo)隨著迭代次數(shù)的增加逐漸趨于水平，最終達(dá)到總體體積容差而終止。迭代次數(shù)由方案2的325次減少為方案3的229次。迭代效率明顯所提高，最終體積為初始體積的0.40。

圖6 方案2的結(jié)果圖

圖7 方案3的結(jié)果圖

方案4與方案2相比，初始權(quán)重系數(shù)ζ0（0）從0.5增大為1.4，權(quán)重系數(shù)增量Δζ不變，ζmax取2.0。圖8a為優(yōu)化后的材料分布圖，材料分布與方案2和方案3相比有很大不同，它的最終體積雖然減少了很多，為原來的0.23，但由應(yīng)力分布圖（圖8b）可知，存在材料的區(qū)域應(yīng)力值集中在［17.1，102.6］MPa之間，應(yīng)力水平的均勻程度不如方案2和方案3，得到的結(jié)果收斂于次優(yōu)解。因為當(dāng)初始權(quán)重系數(shù)ζ0（0）取得過大，迭代初始刪除材料過多，造成了后續(xù)需“硬化”的材料無法恢復(fù)。從圖8c可知，權(quán)重系數(shù)逐階增大直至趨于權(quán)重系數(shù)上限ζmax。參考應(yīng)力在第一階的增幅很大，后續(xù)的增幅與權(quán)重系數(shù)一致。體積指標(biāo)隨著迭代次數(shù)的增加逐漸趨于水平，最終達(dá)到總體體積容差而終止。

圖8 方案4的結(jié)果圖

由以上4種方案的優(yōu)化結(jié)果可知，初始權(quán)重系數(shù)ζ0（0）如果取值過低，迭代過程會被延長，但每次迭代的應(yīng)力變化會更平穩(wěn)，這樣可以避免應(yīng)力變化過大引起誤刪除，故一般在［0.5，1］之間取值。Δζ作為權(quán)重系數(shù)增量，它決定了每階尋優(yōu)過程中參考應(yīng)力增大幅度，為提高效率可適當(dāng)提高權(quán)重系數(shù)增量Δζ，但提高過大也會導(dǎo)致結(jié)構(gòu)的過刪除，一般在［0.05，0.15］之間取值。顯然權(quán)重系數(shù)上限ζmax要大于初始權(quán)重系數(shù)ζ0（0），但取得過大會導(dǎo)致過刪除，一般在［1.0，2.5］之間取值。

綜上所述，在結(jié)構(gòu)優(yōu)化過程中，利用連續(xù)階躍參考應(yīng)力的SKO方法可以更快更好地得到近似最優(yōu)解，并且得到的優(yōu)化結(jié)果邊界清晰，結(jié)構(gòu)合理，應(yīng)力水平更均勻。

3.2 拱橋的拓?fù)鋬?yōu)化

如圖9所示，初始模型為拱橋的平面模型，假設(shè)橋底到橋面頂部的距離為H，橋面的曲率半徑為R，R／H ＝5，橋底左右兩端固支，橋上端面受均布載荷P。

圖9 拱橋平面模型

圖10a為拱橋優(yōu)化后材料分布圖，圖10b為趙州橋?qū)嵕皥D，圖10c為迭代過程中性能參數(shù)的變化圖，初始權(quán)重系數(shù)ζ（0）0為0.5，權(quán)重系數(shù)增量Δζ為0.12，經(jīng)過195次迭代結(jié)束，權(quán)重系逐階增大，達(dá)到ζmax后不再增加，體積指標(biāo)隨著迭代次數(shù)增加而逐漸趨于平穩(wěn)，最后達(dá)到總體體積容差而終止。最終體積為原來的0.41，得到的結(jié)果與實(shí)際拱橋結(jié)構(gòu)相似。

圖10 模型的優(yōu)化結(jié)果圖

3.3 受水平力的直立穿孔板的拓?fù)鋬?yōu)化

如圖11所示，初始模型為一直立的平板，中間有兩長方形孔，H／L ＝2，H1／L1＝2，H／H1＝3，平板下端固定，左側(cè)受兩水平集中載荷P。

圖11 受水平力的直立穿孔板的初始模型

圖12a為受水平力的直立穿孔板的優(yōu)化后材料分布圖。如圖12b所示，初始權(quán)重系數(shù)ζ（0）0為0.5，權(quán)重系數(shù)增量Δζ為0.12，第2、第3階局部尋優(yōu)過程中權(quán)重系數(shù)增大幅度超過權(quán)重系數(shù)增量，這是因為迭代過程中出現(xiàn)相鄰幾次迭代均處于局部最優(yōu)，權(quán)重系數(shù)自動調(diào)節(jié)至更高，以便能進(jìn)入下階尋優(yōu)。隨著迭代次數(shù)的增加，體積變化逐漸趨于平穩(wěn)，最后達(dá)到總體體積容差而終止。最終體積為原來的0.34，結(jié)構(gòu)得到近似最優(yōu)解。

圖12 模型的優(yōu)化結(jié)果圖

4 結(jié)束語

本文提出的基于連續(xù)階躍參考應(yīng)力的SKO方法，能在尋優(yōu)過程中自動調(diào)節(jié)參考應(yīng)力值，有效地解決尋優(yōu)過程中的材料堆積問題，使結(jié)構(gòu)趨于最優(yōu)。算例表明：基于連續(xù)階躍參考應(yīng)力的SKO優(yōu)化方法實(shí)現(xiàn)對平面應(yīng)力的線彈性結(jié)構(gòu)的拓?fù)鋬?yōu)化設(shè)計，方法易于實(shí)現(xiàn)，收斂性好，得到的結(jié)果結(jié)構(gòu)合理，邊界清晰，優(yōu)化后的結(jié)構(gòu)應(yīng)力水平變得更均勻。該方法也同樣適用于三維實(shí)體模型的拓?fù)鋬?yōu)化。

［1]周克民，李俊峰，李霞.結(jié)構(gòu)拓?fù)鋬?yōu)化研究方法綜述［J］.力學(xué)進(jìn)展，2005，35（1）：69－76.

［2]Bendsoe M P，Kikuchi N.Generating Optimal Tpologies in Structural Design Using Homogenization Method［J］.Computer Methods in Applied Mechanics and Engineering，1988，71（2）：197－224.

［3]Xie Y M，Steven G P.A Simple Evolutionary Procedure for Structural Optimization［J］.Computer ＆Structure，1993，49：885－896.

［4]榮見華，謝憶民，姜節(jié)勝，等.漸進(jìn)結(jié)構(gòu)優(yōu)化設(shè)計的現(xiàn)狀與進(jìn)展［J］.長沙交通學(xué)院學(xué)報，2001，17（3）：17－23.

［5]Mattheck C.Design in Nature：Learning from Trees［M］.Berlin：Springer－Verlag，1998.

［6]Harzheim L，Graf G.A Review of Optimization of Cast Parts Using Topology Optimization I－Topology Optimization without Manufacturing Constraints［J］.Structural and Multidisciplinary Optimization，2005，30：491－497.

［7]Harzheim L，Graf G.Topshape：An Attempt to Create Design Proposals including Manufacturing Constraints［J］.International Journal of Vehicle Design，2002，28（4）：389－409.

［8]丁曉紅，程莉.基于SKO方法的滿應(yīng)力結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計［J］.中國機(jī)械工程，2009，20（8）：1765－1770.