基于類桁架的多工況下平面Prager結構拓撲優化

謝汩 周克民

摘要：為開展多工況下曲面結構的拓撲優化，采用兩向正交類桁架連續體材料模型和有限元分析方法，以桿件在結點位置的密度和方向為優化設計變量進行結構優化。根據有限元分析結果，采用滿應力準則法優化各單一工況下的材料分布。按照多工況與各單工況下材料的方向剛度最大值的差值最小為原則，優化多工況下的桿件方向和密度分布。將桿件豎直位置的質心連線作為拓撲優化的平面Prager結構，以3個算例表明該方法的有效性。

關鍵詞：拓撲優化; 類桁架結構; 多工況; 殼體結構; Prager結構

中圖分類號：TU323.4; TP319.99

文獻標志碼：A

Topology optimization of Prager structure under

multiple cases based on truss-like structure

XIE Gu1， ZHOU Kemin2

（

1. Xiamen Xiangan Airport Investment Construction Co.， Ltd.， Xiamen 361006， Fujian， China;

2. College of Civil Engineering， Huaqiao University， Xiamen 361021，Fujian， China）

Abstract：

To promote the topology optimization of curved surface structures under multiple cases， the material model of two-direction orthogonality truss-like continua and the finite element analysis method are selected， and the structure is optimized by taking rod density and rod direction at node position as optimization design variables. According to the results of finite element analysis， the material distribution under

each case is optimized using the full stress criterion method. On the basis of minimum difference of material maximum directional stiffness between multiplecases and each case， the rod direction and density distributions under multiple cases are optimized. The mass centres of vertical position rods are connected to form the planar Prager structure for topology optimization. The effectiveness of the proposed method is indicated by three examples.

Key words：

topology optimization; truss-like structure; multiple cases; shell structure; Prager structure

收稿日期：2019-03-21

修回日期：2019-05-18

基金項目：

國家自然科學基金（11172106， 10872072）

作者簡介：

謝汩（1981—），男，湖南寧鄉人，研究方向為工程技術與結構優化，（E-mail）xiemi88@foxmail.com

0?引?言

空間殼體結構一直備受關注，其以受力合理、質量輕、結構剛度大等特點在現代建筑中得到廣泛應用。在綠色環保成為全球熱點的背景下，對空間殼體結構進行優化設計成為關注的焦點，其重點是確定合理的結構外形和合理的剛度分布。[1]與截面尺?寸和形狀優化比較，殼體結構的拓撲優化相對較困難。[2-4]Prager結構是一種特殊的拓撲優化殼結構，其結構外形為覆蓋在設計平面上一定區域的曲面。[5-7]殼體結構在該區域邊界支撐，載荷沿同一豎直方向，受力點平面位置不變，豎直位置隨結構變化而變化。這符合一般大跨空間殼體結構的受力特點，因此可以通過Prager結構進行殼體結構的拓撲優化。

ROZVANY等[5]在給定的豎向載荷和邊界支撐條件下，建立由垂直拱梁組成的網格式曲殼模型，推導出Prager結構拓撲優化的差分迭代方程。采用差分方法時，假設桿件位于相互平行或垂直的平面內，因此桿件方向不能優化，只能分析簡單載荷情況。Prager的結構形式隨載荷變化而變化，因而與常規結構優化問題[8-9]有較大區別。本文采用類桁架材料模型[10-12]和有限元分析方法，解決具有復雜載荷條件和邊界條件等更一般的優化問題，且為計算簡單僅研究多工況豎向載荷下的平面Prager結構。

1?類桁架連續體材料模型

假設設計域內任一點都有2組正交的桿件分布，并假設這2個正交方向為材料主軸方向，按這2個材料主軸方向建立局部坐標系。

1.1?彈性矩陣

假定這2個材料主軸方向的桿件分布密度分別為t1和t2，彈性模量為E。為避免剛度矩陣奇異，并使材料在t1=t2時各向同性，設剪切等效密度為（t1+?t2）/4。由于平行桿件之間沒有作用力，所以假設橫向變形因數為0。材料主軸方向的彈性矩陣為

D（t1，t2，0）=E

t1t2（t1+t2）/4

（1）

假設材料主軸坐標系與整體坐標系xOy的夾角為α，則整體坐標系下的彈性矩陣為

D（t1，t2，α）=

TT（α）D（t1，t2，0）

T（α）

（2）

式中：T（α）為應變轉換矩陣。式（2）還可以寫為

D（t1，t2，α）=

Etm

1+Rtcos 2α00.5Rtsin 2α

01-Rtcos 2α0.5Rtsin 2α

0.5Rtsin 2α0.5Rtsin 2α0.5

（3）

式中：tm=（t1+t2）/2，Rt=（t1-t2）/（t1+t2）。

1.2?多工況下模型參數的確定

在多工況優化中，材料各結點處可能不止有2組桿件，也可能不互相垂直。為使問題簡化，仍然采用正交材料模型。將式（3）中第1行第1列元素展開可得

D11?（t1，t2，α）=E（（t1+t2）+（t1-t2）cos 2α）/2

（4）

式（4）為材料在x方向產生單位變形時的應力，反映該方向的剛度。假設在任意單工況l（l=1，2，…，L）下的拓撲優化結果中，材料沿αl和αl+π/2方向分別布置有密度為t1，l?和t2，l?的正交桿件，此時材料在任意θ方向的剛度為

D11?（θ;t1，l?，t2，l?，αl）=

E（（t1，l?+t2，l?）+（t1，l?-t2，l?）cos 2（αl-θ））/2

（5）

在所有單工況下，方向剛度的最大值

Sm（θ）=max?l=1，2，…，LD11?（θ;t1，l?，t2，l?，αl）

（6）

式（6）為所有工況下各方向剛度的包絡線。如果在多工況下φ和φ+π/2方向布置有密度為x1和x2的梁，則方向剛度D11?（θ;x1，x2，φ）可以寫為

D11?（θ;x1，x2，φ）=

E（（x1+x2）+（x1-x2）cos 2（φ-θ））/2

（7）

用2-范數定義式（6）和（7）的差值為

δ2≡∫π?0?（D11?（θ;x1，x2，φ）-Sm（θ））2dθ=

∫π?0?（D11?（θ;x1，x2，φ））2dθ-

2∫π?0?D11?（θ;x1，x2，φ）×

Sm（θ）dθ+∫π?0?（Sm（θ））2dθ

（8）

其中：∫π?0?（Sm（θ））2dθ只受各單工況下優化材料的分布參數t1，l?、t2，l?和αl的控制，與多工況下的設計變量x1、x2和φ無關。將式（8）前兩項記為δ21，積分可得

δ21=

πE2（2（x1+x2）2+（x1-x2）2）/8-πE2（（x1+

x2）I0+（x1-x2）（I1cos 2φ+I2sin 2φ））

（9）

其中：

（I0，I1，I2）=

1πE

∫π?0?Sm（θ）（1，cos 2θ，sin 2θ）dθ

（10）

為優化結構材料，使多工況下與各單工況下梁剛度的最大值接近，即使式（9）取極小值。由于

∫π?0?（Sm（θ））2dθ

與多工況下的優化變量無關，因此只需要求式（10）的極小值。對δ21關于（x1+x2）、（x1-x2）和φ分?別求偏導，得

δ21（x1+x2）=（4（x1+x2）-8I0）πE28

（11）

δ21（x1-x2）=（2（x1-x2）-8（I1cos 2φ+

I2sin 2φ））πE28（12）

δ21φ=-16（x1-x2）（-I1sin 2φ+I2cos 2φ）πE28

（13）

令式（11）和（12）等于0，聯立求解得

x1，x2=I0±2（I1cos 2φ+I2sin 2φ）

（14）

令式（13）等于0，求解得

φ=φ0

φ0+π2 ，?φ0=12arctanI2I1

（15）

由此使式（8）或（9）取極小值的解為

φ=

φ0，I1cos 2φ0+I2sin 2φ0≥0

φ0+π2，I1cos 2φ0+I2sin 2φ0<0

（16）

2?有限元分析

2.1?單元剛度矩陣

為編程方便，將式（2）變換為

D（t1，t2，α）=E

2b=1tb3r=1

sbr

gr（α）

（17）

其中：

gr（α）=[cos 2αsin 2α1]

（18）

sbr?=111-1-11

（19）

以桿件在結點位置的密度t1，j?、t2，j?和方向αj作為設計變量，單元內部任意點（ξ，η）的彈性矩陣可由形函數

Nj（ξ，η）

插值得到，即

De（ξ，η）=E

j∈SeNj（ξ，η）

D（t1，j?，t2，j?，αj）=

E

j∈SeNj（ξ，η）

2b=1tbj

3r=1

sbr

gr?（αj）

Ar

（20）

A1=12

1

-10，

A2=14001001110，

A3=121

1

1/2（21）

式中：Se為屬于單元e的結點集合。

由此可以得出單元剛度矩陣為

K=

e

∫?Ve

BT

De

BdV=

Ee

j∈Se

2b=1tbj

3r=1

sbr

gr（αj）

Hejr

（22）

式中：Ve為單元e的體積;

Hejr?為常數矩陣，

Hejr?=E∫Ve?Nj

BT

Ar

BdV

（23）

對于將設計域劃分形成的矩形單元，該常數矩陣與單元無關，可以在有限元分析前求解。同樣，單元內部的桿件密度可由插值得到，因此結構體積為

V=

e

j∈Se

b

∫?Ve?Njtbj?dV=

e

b

j∈Se

tbj

∫?Ve?NjdV=

j

bvjtbj

（24）

式中：

vj=e∈Sj∫Ve?NjdV，

Sj，表示圍繞結點j的單元集合。

對于將設計域劃分形成的矩形單元，vj也是常數。

2.2?求解過程

基于以上推導，可將優化迭代過程歸納如下。

（1）采用有限元劃分設計域。初始化設計參數，即

t0，bj?=1

α0，j?=0 ， b=1，2; j=1，2，…，J

（25）

式中：下標0為迭代指數;J為結點總數。

（2）利用有限元分析，得到各結點位置的主應力方向θk，j?和主應力方向的應變εk，bj?。

（3）采用滿應力準則法優化桿件在結點位置的密度，將桿件方向調整到主應力方向，即

k+1，bj?=Etk，bj?εk，bj?/σp， αk+1，j?=θk，j

（26）

式中：下標k為迭代指數;σp為材料的許用應力。為避免密度過小導致剛度矩陣，限制

tk+1，bj?=max（ξtk+1，m?，k+1，bj?）， tk+1，m?=max?b，j（k+1，bj?）

（27）

式中：ξ為常數，此處取1×10-6?。按照優化后的桿件密度分布調整載荷在豎直方向的分布。

（4）檢查收斂條件，判斷桿件密度在2次迭代中的相對變化量是否足夠小，即

max?b，jtk+1，bj?-tk，bj?/tk+1，m?<δ （28）

式中：δ為一極小數，此處取1×10-4?。若不滿足收斂條件則返回第（2）步。

（5）將桿件密度在豎直方向的質心連線，可以表示最后的拓撲結構。

3?算?例

選擇3個平面問題算例，設計域為4×3，兩端鉸支，載荷只表示力沿水平方向的分布，不表示豎直位置。采用40×30個正方形單元，由于載荷大小、結構尺寸和材料的彈性模量E等皆與拓撲優化結果無關，所以不給出具體數據。

算例1選擇文獻[5]中的單工況算例進行對比，其力學模型見圖1a），載荷為均布力。解析解和桿件質心連線表示的拓撲優化結果見圖1b），這兩者非常接近。

算例2選擇在結構跨中和左右1/4跨分別不同時作用F1和F2，F1=F2/2，其力學模型見圖2a）。桿件豎向密度質心連線表示的拓撲優化結果見圖2b）。

算例3為在結構左、右半跨分別不同時作用均布載荷q1和q2，q1=q2，其力學模型見圖3a）。拓撲優化結果見圖3b）。

算例2和算例3的體積收斂曲線見圖4，可知程序收斂性良好。

4?結束語

采用類桁架材料模型和有限元方法研究多工況下平面Prager結構的拓撲優化方法。分別計算各單工況下結點的剛度矩陣，通過擬合各單工況下材料的最大方向剛度得出多工況下的剛度矩陣，從而得到拓撲優化結構。算例結果表明該方法可行且有效。在工程實際中，空間殼體結構的載荷條件和邊界條件遠比上述討論復雜，復雜拓撲優化問題可考慮將以上方法推廣到空間Prager結構的拓撲優化問題，為網殼結構初步設計提供參考。

參考文獻：

[1]

沈祖炎， 陳揚驥. 網架與網殼[M]. 上海： 同濟大學出版社， 1997： 173-178.

[2]周克民， 李俊峰， 李霞. 結構拓撲優化研究方法綜述[J]. 力學進展， 2005， 35（1）： 69-76. DOI： 10.3321/j.issn：1000-0992.2005.01.007.

[3]KEGL M，BRANK B. Shape optimization of truss-stiffened shell structures with variable thickness[J]. Computer Methods in Applied Mechanics and Engineering， 2006， 195（19-22）： 2611-2634. DOI： 10.1016/j.cma.2005.05.020.

[4]DAVID K.Global laminate optimization on geometrically partitioned shell structures[J]. Structural and Multidisciplinary Optimization， 2011， 43（3）：353-368. DOI： 10.1007/s00158-010-0576-9.

[5]ROZVANY G I N，PRAGER W. A new class of structural optimization problem： Optimal archgrids[J]. Computer Methods in Applied Mechanics and Engineering， 1979， 19（1）： 127-150. DOI： 10.1016/0045-7825（79）90038-0.

[6]ROZVANY G I N，NAKAMURA H， KUHNELL B T. Optimal archgrids： Allowance for selfweight[J]. Computer Methods in Applied Mechanics and Engineering， 1980， 24（3）： 287-304. DOI： 10.1016/0045-7825（80）90066-3.

[7]ROZVANY G I N，WANG C M， DOW M. Prager-structures： Archgrids and cable networks of optimal layout[J]. Computer Methods in Applied Mechanics and Engineering， 1982， 31（1）： 91-113. DOI： 10.1016/0045-7825（82）90049- 4.

[8]WANG C M.Optimization of multispan plane Prager-structures with variable support locations[J]. Engineering Structures， 1987， 9（3）： 157-161. DOI： 10.1016/0141-0296（87）90010-1.

[9]ANSOLA R，CANALES J， TRRAGO J A， et al. An integrated approach for shape and topology optimization of shell structures[J]. Computers and Structures， 2002， 80（5/6）： 449-?458. DOI： 10.1016/S0045-7949（02）00019-6.

[10]ZHOU K M，LI X. Topology optimization of structures under multiple load cases using fiber-reinforced composite material

model[J].Computational Mechanics， 2006， 38（2）： 163-170. DOI： 10.1007/s00466-005-0735-9.

[11]ZHOU K M.Optimization of least-weight grillages by finite element method[J]. Structural and Multidisciplinary Optimization， 2009， 38（5）： 525-532. DOI： 10.1007/s00158-008-0305-9.

[12]ZHOU K M，LI X. Topology optimization of truss-like continua with three families of members model under stress constraints[J]. Structural and Multidisciplinary Optimization， 2011， 43（4）： 487-?493. DOI： 10.1007/s00158-010-0584-9.

（編輯?武曉英）