凡健 易振峰 姚興智 謝錦鵬 譚文超 王昱

摘要：為探究溫度變化對溫室大棚骨架結(jié)構(gòu)安全性的影響，建立考慮機械載荷和熱載荷的大棚骨架熱彈性結(jié)構(gòu)優(yōu)化設計模型。以機械載荷和熱載荷下結(jié)構(gòu)的最大應力值最小化為目標，以結(jié)構(gòu)的總材料用量為約束，考慮機械載荷與熱載荷聯(lián)合作用下的應力分布，實現(xiàn)對大棚骨架的連續(xù)體結(jié)構(gòu)優(yōu)化設計，使結(jié)構(gòu)在滿足支撐剛度的前提下最大程度緩和結(jié)構(gòu)的應力集中問題?？紤]到熱應力優(yōu)化問題中的設計依賴性和中間變量問題，采用密度過濾函數(shù)獲得清晰的最優(yōu)拓撲構(gòu)型。通過兩種典型溫室大棚骨架優(yōu)化算例，對比不同溫度變化幅度和不同材料用量下優(yōu)化結(jié)構(gòu)的最大等效應力，結(jié)果表明：該模型在相同體積比下結(jié)構(gòu)最大等效應力優(yōu)化效率最高可達到約15%，同種載荷條件下體積比增加0.1可實現(xiàn)結(jié)構(gòu)最大等效應力優(yōu)化效果提高近1%。該研究為考慮機械載荷和熱載荷的大棚骨架熱彈性結(jié)構(gòu)優(yōu)化設計提供有效的設計方法，對工程應用中的溫室大棚骨架結(jié)構(gòu)設計具有指導意義。

關鍵詞：溫室結(jié)構(gòu)；拓撲優(yōu)化；熱彈性結(jié)構(gòu)；應力

中圖分類號：S625.1

文獻標識碼：A

文章編號：2095-5553 （2024） 06-0077-06

收稿日期：2022年10月28日

修回日期：2022年12月15日

*基金項目：國家自然科學基金資助項目（51705161）；廣州市基礎與應用基礎研究項目（202102020870）

第一作者：凡健，男，1996年生，安徽滁州人，碩士研究生；研究方向為結(jié)構(gòu)拓撲優(yōu)化、結(jié)構(gòu)設計。E-mail： fanjian@stu.scau.edu.cn

通訊作者：王昱，女，1987年生，長沙人，博士，副教授；研究方向為結(jié)構(gòu)拓撲優(yōu)化、輕量化。E-mail： yu-wang@scau.edu.cn

Lightweight design of greenhouse frame structure considering mechanical load and thermal load

Fan Jian^{1， 2}， Yi Zhenfeng^{1， 2}， Yao Xingzhi^{1， 2}， Xie Jinpeng¹， Tan Wenchao^{1， 2}， Wang Yu^{1， 2}

（1. College of Engineering， South China Agricultural University， Guangzhou， 510642， China;

2. Key Laboratory of Key Technology on Agricultural Machinery and Equipment， Ministry of Education，

South China Agricultural University， Guangzhou， 510642， China）

Abstract： In order to explore the influence of temperature change on the safety of greenhouse skeleton structure， an optimal design model of greenhouse skeleton thermoelastic structure considering mechanical load and thermal load was established. With the goal of minimizing the maximum stress value of the structure under mechanical load and thermal load， and choosing the total material consumption of the structure as the constraint， the stress distribution under the combined action of mechanical load and thermal load was considered to achieve the optimal design of the continuum structure of the greenhouse framework， so that the structure could minimize the stress concentration of the structure under the premise of meeting the support stiffness. Considering the design dependency and intermediate variables in the thermal stress optimization problem， the density filter function was used to obtain a clear optimal topology. Through two typical greenhouse frame optimization examples， the maximum equivalent stress of the optimized structure under different temperature changes and material consumption was compared. The results showed that the maximum equivalent stress optimization efficiency of the model could reach about 15% under the same volume fraction. Under the same load conditions， increasing the volume ratio by 0.1 could achieve an increase of nearly 1% in the maximum equivalent stress optimization effect of the structure. The obtained conceptual design scheme of the greenhouse skeleton structure has guiding significance for the design of the greenhouse skeleton structure in engineering applications.

Keywords： greenhouse structure; topology optimization; thermoelastic structure; stress

0 引言

溫室大棚作為一種作物栽培設施，因其建造和運行成本低，在我國北方迅速發(fā)展起來，成為中國設施農(nóng)業(yè)的主體^[1]。按照不同的功能需求，可將常見溫室大棚分為兩大類：連棟溫室和日光溫室。長期以來，溫室結(jié)構(gòu)的有關研究主要集中在溫室結(jié)構(gòu)受到各種機械載荷工況下結(jié)構(gòu)的安全性研究^[2-4]，忽略了晝夜溫差變化產(chǎn)生的熱應力對大棚結(jié)構(gòu)性能的影響。當熱應力的存在使結(jié)構(gòu)內(nèi)部應力超過許用應力時，則導致結(jié)構(gòu)破壞，給現(xiàn)實農(nóng)業(yè)生產(chǎn)帶來安全隱患^[5]。

除了在典型溫室結(jié)構(gòu)的基礎上進行尺寸優(yōu)化外，拓撲優(yōu)化作為力學工程領域的先進設計方法^[6]，近年來發(fā)展了一系列的優(yōu)化模型，被應用于解決土木、車輛、儀器設備等^[7-9]工程結(jié)構(gòu)的機械載荷和熱載荷作用下的應力問題。桁架結(jié)構(gòu)是溫室大棚骨架結(jié)構(gòu)中最為常見的形式，通過拓撲優(yōu)化方法可合理布置桁架單元的空間分布，可以提升其承載能力^[10]。然而，鮮有文獻在針對溫室大棚骨架結(jié)構(gòu)研究的同時，綜合考慮機械載荷和熱載荷的結(jié)構(gòu)應力問題。

本文將基于連續(xù)體拓撲優(yōu)化方法，針對設施農(nóng)業(yè)對大棚骨架結(jié)構(gòu)的剛度和熱變形以及輕量化的工程需求，建立考慮機械載荷和熱載荷作用下結(jié)構(gòu)應力分布的設計模型，在輕量化的同時使得結(jié)構(gòu)具有較高的剛度和較小的熱變形。為驗證該優(yōu)化模型的可行性，對兩種典型溫室大棚骨架優(yōu)化設計算例，針對常見溫室主要承重結(jié)構(gòu)和工作環(huán)境特點給出優(yōu)化設計方案，獲得滿足實際農(nóng)業(yè)生產(chǎn)需要的結(jié)構(gòu)。

1 考慮機械載荷和熱載荷的結(jié)構(gòu)優(yōu)化模型

1.1 連續(xù)體結(jié)構(gòu)拓撲優(yōu)化方法

當前連續(xù)體結(jié)構(gòu)優(yōu)化領域的主要方法有：尺寸優(yōu)化、形狀優(yōu)化和拓撲優(yōu)化。拓撲優(yōu)化作為一種啟發(fā)式的結(jié)構(gòu)設計方法，指在滿足已知的約束、載荷和邊界條件的情況下，通過將所求結(jié)構(gòu)的某種力學性能轉(zhuǎn)化成目標函數(shù)，建立其與結(jié)構(gòu)參數(shù)的數(shù)學模型，進而采用數(shù)值優(yōu)化算法優(yōu)計算出最優(yōu)的結(jié)構(gòu)參數(shù)，最終獲得最優(yōu)結(jié)構(gòu)的過程。主流的拓撲優(yōu)化方法可分為以均勻化為主的材料插值類方法^[11]和描述結(jié)構(gòu)形狀函數(shù)類方法，如水平集法^[12]。結(jié)構(gòu)的拓撲優(yōu)化旨在設計區(qū)域內(nèi)尋找一個給定體積的子區(qū)域，使得該區(qū)域?qū)哪繕撕瘮?shù)（如結(jié)構(gòu)柔順度、結(jié)構(gòu)位移等）取得極值。引入離散變量的材料密度函數(shù)，如式（1）所示。

x=1x∈Ω_solid

0x∈Ω_void（1）

式中： Ω_solid——實體材料域；

Ω_void——孔洞材料域。

由于整數(shù)模型的計算求解非常困難，通常采用變量連續(xù)化方法將0～1整數(shù)變量問題變?yōu)?、1間的連續(xù)變量優(yōu)化模型。上述變量連續(xù)化后的模型是一個病態(tài)問題，優(yōu)化解表現(xiàn)為中間密度值、棋盤格和網(wǎng)格依賴性，可采用中間密度懲罰模型和棋盤格控制措施解決上述問題。常用的中間材料懲罰模型有SIMP（Solid Isotropic Material with Penalization Model）^[13]和RAMP（Rational Approximation of Material Properties）^[14]模型。由于熱載荷問題的復雜性，本文的拓撲優(yōu)化模型基于RAMP模型建立，將在材料插值模型部分詳細闡述。

1.2 熱彈性結(jié)構(gòu)優(yōu)化模型

考慮機械載荷和熱載荷的的熱彈性結(jié)構(gòu)輕量化設計問題，結(jié)構(gòu)變形為機械載荷與熱載荷聯(lián)合作用下的總體變形，以結(jié)構(gòu)的最大應力最小化作為目標函數(shù)，以結(jié)構(gòu)的體積比（即總材料用量）為約束，建立相應的熱彈性結(jié)構(gòu)優(yōu)化模型如式（2）所示，示意圖見圖1。

本文模型采用優(yōu)化領域廣泛應用的漸進優(yōu)化法MMA（Method of moving asymptotes）^[15]進行求解。

minσ_max=σ_pn

s.t.KU=F_m+F_th

V=∑Ni=1x_iv₀≤V₀

0min≤x_i≤1 i=1，2，3，…，N（2）

式中： σ_max——結(jié)構(gòu)最大應力；

σ_pn——結(jié)構(gòu)等效最大應力；

K——總剛度矩陣；

U——總節(jié)點位移矩陣；

F_m——作用在結(jié)構(gòu)上的機械載荷；

F_th——變化溫度場產(chǎn)生的熱載荷；

V——結(jié)構(gòu)優(yōu)化后的體積；

x_i——第i個單元的單元密度；

v₀——單元的體積；

V₀——給定材料體積上限；

x_min——單元的最小密度值，取0.001以避免總剛度矩陣的奇異；

N——設計變量數(shù)目。

2 優(yōu)化模型的實施

2.1 應力計算

針對基于密度的最大應力最小化結(jié)構(gòu)優(yōu)化問題，拓撲優(yōu)化迭代過程中出現(xiàn)的奇異現(xiàn)象存在于基于密度的拓撲優(yōu)化中，即當單元密度為很低的值時依然存在應變導致局部出現(xiàn)人造大應力的現(xiàn)象。應力松弛方法^[16]以減輕低密度單元的應力，可使優(yōu)化穩(wěn)定進行。同時，在衡量結(jié)構(gòu)的應力時，只能得到每個單元的局部應力大小，無法建立結(jié)構(gòu)整體應力水平和局部單元應力大小的聯(lián)系，但考察每一個單元則大大加重了計算成本。因此，使用P范數(shù)近似應力最大值^[17]，將各個單元應力分量整合成一個衡量值，提高計算效率。

對于熱彈性結(jié)構(gòu)而言，第i個單元的應力矢量σ_i如式（3）所示。

σ_i=D₀（BU_i-α₀ΔT?^T）（3）

式中： D₀——固體材料的彈性矩陣；

B——應變—位移矩陣；

α₀——固體材料的熱膨脹系數(shù)；

ΔT——結(jié)構(gòu)受到的均勻溫升場；

?^T——列向量，?=[1 1 0]。

對于二維平面問題，單元應力σ_i是一個3×1的列向量，不便于優(yōu)化的開展。采用von·mises應力準則，將σ_i整合成一個等效應力σ^vm_i。

σ^vm_i=（σ_i1²+σ_i2²+3σ_i3²－σ_i1σ_i2）¹²（4）

式中： σ_i1、σ_i2、σ_i3——σ_i的三個分量，表示為水平方向應力、豎直方向應力和切應力。

為解決應力奇異問題，這里引入q松弛法對單元等效應力進行懲罰

σ^vm_i=x^q_iσ^vm_i（5）

式中： q——懲罰因子，這里取0.5。

為了實現(xiàn)目標函數(shù)連續(xù)可導，采用P范數(shù)近似全局最大應力

σ_pn=∑Ni=1σ^vm_i^P1P（6）

式中： P——范數(shù)系數(shù)，本文P取8。

2.2 材料插值

熱彈性結(jié)構(gòu)相對傳統(tǒng)只受機械載荷的優(yōu)化問題更加復雜，因其敏度與設計變量有關，在優(yōu)化過程中存在符號改變，給優(yōu)化過程帶來不穩(wěn)定性。傳統(tǒng)的SIMP插值方法在表示單元偽密度的設計變量x趨近于0時，其導數(shù)為0，這阻礙了熱應力敏度符號的轉(zhuǎn)變，影響優(yōu)化過程的順利進行。因此，本文引入RAMP插值方法分別懲罰剛度和熱應力。

對于第i個單元，其單元剛度矩陣和節(jié)點熱應力矢量表達如式（7）所示。

K_i=∫ΩB_i^TD_iB_idΩ（7）

F^th_i=∫ΩB_i^TD_iε^th_idΩ（8）

式中： D_i——密度為x_i的單元等效彈性矩陣；

ε^th_i——單元i的熱應變矩陣。

D_i=E（x_i）D₀（9）

ε^th_i=α（x_i）ΔT?^T（10）

式中： E（x_i）——材料的插值函數(shù)；

α（x_i）——單元密度為x_i時的等效熱膨脹系數(shù)。

因此，式（8）可被改寫為

F^th_i=E（x_i）α（x_i）ΔT∫ΩB_i^TD_i?^TdΩ（11）

其中只有α_i和D_i與變量x_i有關，引入表征材料固有特性的熱應力系數(shù)（TSC）^[18]概念來表達三者關系。

β（x_i）=E（x_i）α（x_i）（12）

根據(jù)RAMP插值方法，等效彈性矩陣D_i和熱應力系數(shù)β（x_i）可以表達成

E（x_i）=x_i1+S_E（1－x_i）E₀（13）

β（x_i）=x_i1+S_β（1－x_i）E₀α₀（14）

式中： S_E——楊氏模量插值系數(shù)；

S_β——TSC插值系數(shù)；

E₀——固體材料的楊氏模量；

α₀——固體材料的熱膨脹系數(shù)。

S_E和S_β值的選取影響迭代的穩(wěn)定性，并且當S_E＜S_β時迭代產(chǎn)生震蕩^[19]。根據(jù)數(shù)值試驗經(jīng)驗，這里選取S_E=8，S_β=0。

2.3 密度過濾

基于變密度法的拓撲優(yōu)化方法是通過將0-1離散變量轉(zhuǎn)化成連續(xù)變量，這將導致在優(yōu)化過程中產(chǎn)生無物理意義的具有中間密度的單元。選用合理的材料插值模型參數(shù)，可懲罰中間密度往0-1方向逼近。在產(chǎn)生數(shù)學層面上局部最優(yōu)解的同時，結(jié)構(gòu)的局部區(qū)域伴隨產(chǎn)生材料密度為0和1周期分布的狀態(tài)，即棋盤格現(xiàn)象，增加了結(jié)構(gòu)加工成型的難度。針對上述問題以及2.1節(jié)提到的應力非線性問題，單元的應力極易受其相鄰元素的影響。通過施加密度過濾，可有效防止棋盤格現(xiàn)象并且降低相鄰元素對于單元應力的影響，使優(yōu)化收斂于全局解。

本研究中采用的密度過濾可定義為

x～_{（k，l）}=∑m，nH（m，n）x_{（k－m，l－n）}∑m，nH（m，n）（15）

式中： x_（i，j）——第i行第j列的單元密度；

x～_{（i，j）}——第i行第j列單元過濾后的偽密度；

H（m，n）——離散變量m，n的過濾核函數(shù)。

H（m，n）=max[0，r_min－δ（m，n）]（16）

式中： δ（m，n）——相距m行n列兩單元的中心距離；

r_min——最小過濾半徑。

r_min的大小影響著過濾后結(jié)構(gòu)細節(jié)枝干的數(shù)量，同時決定結(jié)構(gòu)的最小尺寸。

3 典型大棚骨架結(jié)構(gòu)輕量化設計

針對實際農(nóng)業(yè)應用領域兩種典型溫室結(jié)構(gòu)：連棟溫室和日光溫室，如圖2所示，分別選取其主要承重結(jié)構(gòu)作為設計域進行優(yōu)化。建立兩種溫室主體框架部分二維模型，以均勻密度場作為初始設計，設定外圍輪廓作為非設計域，代入本文的模型進行優(yōu)化。選擇現(xiàn)實生活中常用的Q235鋼材作為材料，其楊氏模量E為200GPa，泊松比μ為0.3，熱膨脹系數(shù)α為1.20×10^-5，許用應力為235MPa。

3.1 連棟溫室主體承重結(jié)構(gòu)優(yōu)化

連棟溫室主體結(jié)構(gòu)部分如圖3所示。該部分為軸對稱結(jié)構(gòu)，為減少優(yōu)化過程的計算量，截取框架右半邊作為設計區(qū)域。結(jié)構(gòu)長l=800cm，高h=400cm，厚度t取1cm，非設計域（圖3中黑色邊框）寬度為20cm，整個區(qū)域離散成80×40的平面四節(jié)點有限單元。結(jié)構(gòu)同時受到均勻溫度場ΔT和分布載荷t-作用，分布載荷隨著外輪廓傾角的變化呈1000N/m到10N/m漸變。右下方采用固定約束，選擇0.3作為體積比。為探究不同溫度變化對結(jié)構(gòu)應力的影響，ΔT分別選擇0℃、10℃、20℃和50℃，優(yōu)化后拓撲結(jié)構(gòu)和應力分布見圖4，優(yōu)化結(jié)果見表1（表中表達式上標init和opt分別表示初始構(gòu)型和優(yōu)化構(gòu)型），η表示優(yōu)化效率，計算如式（17）所示。

η=σ^init_pn－σ^opt_pnσ^init_pn（17）

由表1可知，當連棟溫室結(jié)構(gòu)只受到機械載荷時，優(yōu)化后的結(jié)構(gòu)與實際農(nóng)業(yè)生產(chǎn)應用的溫室結(jié)構(gòu)相似，其應力水平遠低于材料的許用應力。加入熱載荷后，結(jié)構(gòu)的最大等效應力根據(jù)溫升大小的變化同幅度的增加，伴隨著優(yōu)化效率的小幅提升。從圖4可知，等效最大應力值大于結(jié)構(gòu)實際應力最大值，對于P范數(shù)等效最大應力的方法造成的與真實應力的差距，可通過增大范數(shù)系數(shù)P的值縮小。

此外，升溫幅度的不同，結(jié)構(gòu)拓撲也有相應的改變，但始終保持三條棱的樣式。溫升為10℃時，優(yōu)化后的結(jié)構(gòu)應力水平在許用應力之下；20℃及以上的溫升，結(jié)構(gòu)的最大應力則嚴重超過了許用應力。研究發(fā)現(xiàn)，優(yōu)化后的結(jié)構(gòu)最大應力出現(xiàn)在右下角支撐區(qū)域處，通過對初始框架模型右下角進行局部優(yōu)化后再進行拓撲優(yōu)化，可有效降低結(jié)構(gòu)應力集中。

3.2 日光溫室主體承重結(jié)構(gòu)優(yōu)化

日光溫室主體框架部分如圖5所示。結(jié)構(gòu)長l=800cm，高h=400cm，厚度t取1cm，非設計域?qū)挾葹?0cm。整個區(qū)域同樣離散成80×40的平面四節(jié)點有限單元。左側(cè)斜坡受到均布載荷800N/m，右側(cè)斜坡受到均布載荷t-為100N/m。為探究不同體積比對優(yōu)化結(jié)果的影響，設定溫升ΔT為10℃不變，分別選取體積比等于0.2、0.3、0.4，對應優(yōu)化后的拓撲結(jié)構(gòu)和應力分布見圖6，結(jié)果如表2所示。

通過表2可知，在同一溫升下，體積比增加0.1，優(yōu)化后結(jié)構(gòu)的最大等效應力減小了約1%且幅度逐漸變小。因此，在實際農(nóng)業(yè)生產(chǎn)應用中進行優(yōu)化設計時，應當權(quán)衡好應力大小和材料用量的關系，選擇最合適的參數(shù)組合。

4 結(jié)論

本文建立考慮機械載荷和熱載荷的大棚骨架熱彈性結(jié)構(gòu)優(yōu)化設計模型，將其應用于溫室大棚骨架的輕量化設計，同時探究熱載荷對溫室大棚骨架的結(jié)構(gòu)應力分布的影響。通過探究不同溫升和不同材料用量對優(yōu)化后結(jié)構(gòu)應力狀況的影響。

1） 熱載荷的加入大幅提高結(jié)構(gòu)的最大應力：在10℃溫升時，結(jié)構(gòu)最大應力小于材料許用應力，設計滿足使用要求；溫升在20℃及以上時，最大應力超出許用應力，最大應力位置不變指導設計需要對結(jié)構(gòu)初始進行改良。

2） 同等溫升條件下，總材料用量增加0.1，結(jié)構(gòu)的最大等效應力減小約1%。因此，在設計時應結(jié)合實際條件需要，選擇出最佳的結(jié)構(gòu)應力大小和總的材料用量組合。

本文不僅為考慮機械載荷和熱載荷的大棚骨架熱彈性結(jié)構(gòu)優(yōu)化設計提供有效的設計方法，還給出具有指導意義的符合工程需求的大棚骨架結(jié)構(gòu)概念設計方案。

參 考 文 獻

[1]劉志杰， 鄭文剛， 胡清華， 等. 中國日光溫室結(jié)構(gòu)優(yōu)化研究現(xiàn)狀及發(fā)展趨勢[J]. 中國農(nóng)學通報， 2007（2）： 449-453.

Liu Zhijie， Zheng Wengang， Hu Qinghua， et al. Current situation and development on structure optimization of solar greenhouse in China [J]. Chinese Agricultural Science Bulletin， 2007（2）： 449-453.

[2]肖林剛， 鄒志榮， 吳樂天， 等. 寒冷干旱地區(qū)日光溫室結(jié)構(gòu)的優(yōu)化設計[J]. 農(nóng)學學報， 2014， 4（5）： 52-55.

Xiao Lingang， Zou Zhirong， Wu Letian， et al. Structural optimization and design of the sunlight greenhouse in cold and arid Regions [J]. Journal of Agriculture， 2014， 4（5）： 52-55.

[3]劉麗霞， 雷法究. 兩種典型日光溫室結(jié)構(gòu)安全性能分析[J]. 林業(yè)機械與木工設備， 2016， 44（11）： 18-23.

Liu Lixia， Lei Fajiu. Safety performance analysis of two typical greenhouse structure [J]. Forestry Machinery & Woodworking Equipment， 2016， 44（11）： 18-23.

[4]Briassoulis D， Dougka G， Dimakogianni D， et al. Analysis of the collapse of a greenhouse with vaulted roof [J]. Biosystems Engineering， 2016， 151： 495-509.

[5]李世萍， 秦小楠， 馬倩， 等. 北京力學會第18屆學術年會論文集： 工程應用： 輕鋼溫室結(jié)構(gòu)溫度應力及其對策[C]. 北京： 北京力學會， 2012.

Li Shiping， Qin Xiaonan， Ma Qian， et al. Proceedings of the 18th annual academic conference of beijing mechanics society： Engineering application： Temperature stress of light steel greenhouse structure and its countermeasures [C]. Beijing： Beijing Society of Theoretical and Applied Mechanics， 2012.

[6]Bends?e M， Sigmund O. Topology optimization： Theory， method and applications [M]. New York： Springer-Verlag Berlin Heideberg， 2003.

[7]He Y， Cai K， Zhao Z L， et al. Stochastic approaches to generating diverse and competitive structural designs in topology optimization [J]. Finite Elements in Analysis and Design， 2020， 173： 103399.

[8]Szepessy Z， Zoltan I. Thermal dynamic model of precision wire-wound resistors [J]. IEEE Transactions on Instrumentation and Measurement， 2002， 51（5）： 930-934.

[9]Mankame N D， Ananthasuresh G K. Topology synthesis of electrothermal compliant mechanisms using line elements [J]. Structural and Multidisciplinary Optimization， 2004， 26（3）： 209-218.

[10]Xu L， Min H. Optimum design of cold formed steel residential roof trusses [J]. International Specialty Conference on Cold-Formed Steel Structures， 2000.

[11]Suzuki K， Kikuchi N. A homogenization method for shape and topology optimization [J]. Computer Methods in Applied Mechanics and Engineering， 1991， 93（3）： 291-318.

[12]Wang M Y， Wang X， Guo D. A level set method for structural topology optimization [J]. Computer Methods in Applied Mechanics and Engineering， 2003， 192（1）： 227-246.

[13]Rietz A. Sufficiency of a finite exponent in SIMP （power law） methods [J]. Structural and Multidisciplinary Optimization， 2001， 21（2）： 159-163.

[14]Stolpe M， Svanberg K. An alternative interpolation scheme for minimum compliance topology optimization [J]. Structural and Multidisciplinary Optimization， 2001， 22（2）： 116-124.

[15]Svanberg K. The method of moving asymptotes-a new method for structural optimization [J]. International Journal for Numerical Methods in Engineering， 1987， 24（2）： 359-373.

[16]Cheng G D， Guo X. ε-relaxed approach in structural topology optimization [J]. Structural optimization， 1997， 13（4）： 258-266.

[17]Duysinx P， Bends?e M P. Topology optimization of continuum structures with local stress constraints [J]. International Journal for Numerical Methods in Engineering， 1998， 43（8）： 1453-1478.

[18]Deaton J D， Grandhi R V. Stiffening of restrained thermal structures via topology optimization [J]. Structural and Multidisciplinary Optimization， 2013， 48（4）： 731-745.

[19]Gao T， Zhang W. Topology optimization involving thermo-elastic stress loads [J]. Structural and Multidisciplinary Optimization， 2010， 42（5）： 725-738.