基于連續階躍參考應力的SKO結構拓撲優化方法

侯劍云 丁曉紅

上海理工大學，上海，200093

0 引言

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

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

局部最優是一種尋優過程中的早熟現象，需要避免。為此，本文提出一種基于連續階躍參考應力的尋優迭代策略，可有效地改善早熟現象，得到最優結構。

1 連續階躍參考應力SKO方法介紹

1.1 基本原理介紹

假定結構由不同的材料組成，將每個單元的彈性模量作為參數來改變結構的拓撲形態。材料的彈性模量被定義為溫度的函數，即隨著溫度的升高，彈性模量變小，材料被“軟化”，當材料“軟”到一定程度，可認為材料被刪除；同時隨著溫度的降低，材料的彈性模量變大，材料被“硬化”。溫度T和彈性模量E間的關系可假設為線性關系，如圖1所示，此處的溫度沒有物理意義，僅是單元彈性模量改變的控制器。在應用連續階躍參考應力的SKO方法的過程中，會經過多階的局部尋優過程，每階尋優都經過數次迭代。溫度的迭代公式為

式中，上標i、j表示第i階尋優過程中第j次迭代；Tn（i，j）為第i階尋優過程中第j次迭代時n節點的溫度；σn（i，j－1）為第i階尋優過程中第j－1次迭代時n節點的等效節點應力；σ（i，j）為第i階局部尋優過程中第j次迭代的連續階躍參考r應力；s為步長因子；Tmax、Tmin為迭代時的溫度上下限，上限為100，下限為0。

圖1 溫度與彈性模量關系圖

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

1.2 局部尋優過程和連續階躍參考應力

在SKO方法的尋優過程中，有可能出現局部最優結構，它存在著一些材料堆積區域，此時若提高參考應力，則尋優過程又可繼續進行。本文提出一種連續階躍參考應力的概念，使結構能在尋優過程中自動調節參考應力的取值，從而逐步使結構趨優。其計算公式為

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

式 中 ，ζ0（i）為第i階尋優權重系數初值；α為常數，取0.1；Δζ為權重系數增量。

第i階權重系數初值ζ0（i）為第i－1階尋優過程中最后一次迭代的權重系數，即

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

0重系數均滿足以上所述的提升方式，而區域2經過的兩次迭代都達到局部最優，說明權重系數增大得不夠，需再次提升，才能進入下一階尋優過程，這是權重系數的自動調節過程。

圖2 權重系數ζ變化圖

1.3 局部最優判斷標準和收斂條件

在尋優過程中，當結構前后兩次迭代的體積變化小于給定的局部體積容差ε1時，認為結構達到一個局部最優狀態；當小于給定的總體體積容差ε時，認為得到最終的優化結果。顯然，ε1要大于ε。實驗證明，ε1取100～1000倍的ε時可確保數值穩定。這里的體積指的是量綱一化的名義體積。SKO方法中結構的設計變量是彈性模量E，定義第i階尋優過程中第j次迭代的名義體積為

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

將名義體積量綱一化得到第i階尋優過程中第j次迭代時量綱一體積指標，其公式為

2 算法流程

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

（1）建立初始物理模型，并確定設計空間，實際結構中某些具有特殊要求的邊界或面需保留而不作為設計區域。

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

（3）進行結構線性靜力有限元分析，提取節點應力，通過式（3）～式（5）計算連續階躍參考應力，權重系數ζ（i，j）不得超過規定的權重系數上限ζmax。

圖3 算法流程

（4）按式（1）、式（2）實現“軟硬化”單元的操作。

（5）未達到局部最優時，更新j，重復步驟（3），當達到局部最優結構后，更新i、j，再重復步驟（3），直到滿足上文所述的迭代終止條件，終止迭代。

3 設計參數選擇與算例

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

3.1 受兩面內力懸臂梁的拓撲優化

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

圖4 長方形懸臂梁受兩面內力模型

權重系數參數分別按表1所示的4種方案進行選取，局部體積容差ε1和總體體積容差ε如表1所示。

表1 參數取值表

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

圖5 方案1的結果圖

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

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

圖6 方案2的結果圖

圖7 方案3的結果圖

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

圖8 方案4的結果圖

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

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

3.2 拱橋的拓撲優化

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

圖9 拱橋平面模型

圖10a為拱橋優化后材料分布圖，圖10b為趙州橋實景圖，圖10c為迭代過程中性能參數的變化圖，初始權重系數ζ（0）0為0.5，權重系數增量Δζ為0.12，經過195次迭代結束，權重系逐階增大，達到ζmax后不再增加，體積指標隨著迭代次數增加而逐漸趨于平穩，最后達到總體體積容差而終止。最終體積為原來的0.41，得到的結果與實際拱橋結構相似。

圖10 模型的優化結果圖

3.3 受水平力的直立穿孔板的拓撲優化

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

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

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

圖12 模型的優化結果圖

4 結束語

本文提出的基于連續階躍參考應力的SKO方法，能在尋優過程中自動調節參考應力值，有效地解決尋優過程中的材料堆積問題，使結構趨于最優。算例表明：基于連續階躍參考應力的SKO優化方法實現對平面應力的線彈性結構的拓撲優化設計，方法易于實現，收斂性好，得到的結果結構合理，邊界清晰，優化后的結構應力水平變得更均勻。該方法也同樣適用于三維實體模型的拓撲優化。

［1]周克民，李俊峰，李霞.結構拓撲優化研究方法綜述［J］.力學進展，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]榮見華，謝憶民，姜節勝，等.漸進結構優化設計的現狀與進展［J］.長沙交通學院學報，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方法的滿應力結構拓撲優化設計［J］.中國機械工程，2009，20（8）：1765－1770.