基于類桁架的多工況下平面Prager結(jié)構(gòu)拓?fù)鋬?yōu)化

謝汩 周克民

摘要：為開展多工況下曲面結(jié)構(gòu)的拓?fù)鋬?yōu)化，采用兩向正交類桁架連續(xù)體材料模型和有限元分析方法，以桿件在結(jié)點(diǎn)位置的密度和方向?yàn)閮?yōu)化設(shè)計(jì)變量進(jìn)行結(jié)構(gòu)優(yōu)化。根據(jù)有限元分析結(jié)果，采用滿應(yīng)力準(zhǔn)則法優(yōu)化各單一工況下的材料分布。按照多工況與各單工況下材料的方向剛度最大值的差值最小為原則，優(yōu)化多工況下的桿件方向和密度分布。將桿件豎直位置的質(zhì)心連線作為拓?fù)鋬?yōu)化的平面Prager結(jié)構(gòu)，以3個(gè)算例表明該方法的有效性。

關(guān)鍵詞：拓?fù)鋬?yōu)化; 類桁架結(jié)構(gòu); 多工況; 殼體結(jié)構(gòu); Prager結(jié)構(gòu)

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

文獻(xiàn)標(biāo)志碼：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，F(xiàn)ujian， 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

基金項(xiàng)目：

國家自然科學(xué)基金（11172106， 10872072）

作者簡介：

謝汩（1981—），男，湖南寧鄉(xiāng)人，研究方向?yàn)楣こ碳夹g(shù)與結(jié)構(gòu)優(yōu)化，（E-mail）xiemi88@foxmail.com

0?引?言

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

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

1?類桁架連續(xù)體材料模型

假設(shè)設(shè)計(jì)域內(nèi)任一點(diǎn)都有2組正交的桿件分布，并假設(shè)這2個(gè)正交方向?yàn)椴牧现鬏S方向，按這2個(gè)材料主軸方向建立局部坐標(biāo)系。

1.1?彈性矩陣

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

D（t1，t2，0）=E

t1t2（t1+t2）/4

（1）

假設(shè)材料主軸坐標(biāo)系與整體坐標(biāo)系xOy的夾角為α，則整體坐標(biāo)系下的彈性矩陣為

D（t1，t2，α）=

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

T（α）

（2）

式中：T（α）為應(yīng)變轉(zhuǎn)換矩陣。式（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?多工況下模型參數(shù)的確定

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

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

（4）

式（4）為材料在x方向產(chǎn)生單位變形時(shí)的應(yīng)力，反映該方向的剛度。假設(shè)在任意單工況l（l=1，2，…，L）下的拓?fù)鋬?yōu)化結(jié)果中，材料沿αl和αl+π/2方向分別布置有密度為t1，l?和t2，l?的正交桿件，此時(shí)材料在任意θ方向的剛度為

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）為所有工況下各方向剛度的包絡(luò)線。如果在多工況下φ和φ+π/2方向布置有密度為x1和x2的梁，則方向剛度D11?（θ;x1，x2，φ）可以寫為

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

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

（7）

用2-范數(shù)定義式（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θ只受各單工況下優(yōu)化材料的分布參數(shù)t1，l?、t2，l?和αl的控制，與多工況下的設(shè)計(jì)變量x1、x2和φ無關(guān)。將式（8）前兩項(xiàng)記為δ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）

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

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

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

δ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，聯(lián)立求解得

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）

以桿件在結(jié)點(diǎn)位置的密度t1，j?、t2，j?和方向αj作為設(shè)計(jì)變量，單元內(nèi)部任意點(diǎn)（ξ，η）的彈性矩陣可由形函數(shù)

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的結(jié)點(diǎn)集合。

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

K=

e

∫?Ve

BT

De

BdV=

Ee

j∈Se

2b=1tbj

3r=1

sbr

gr（αj）

Hejr

（22）

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

Hejr?為常數(shù)矩陣，

Hejr?=E∫Ve?Nj

BT

Ar

BdV

（23）

對于將設(shè)計(jì)域劃分形成的矩形單元，該常數(shù)矩陣與單元無關(guān)，可以在有限元分析前求解。同樣，單元內(nèi)部的桿件密度可由插值得到，因此結(jié)構(gòu)體積為

V=

e

j∈Se

b

∫?Ve?Njtbj?dV=

e

b

j∈Se

tbj

∫?Ve?NjdV=

j

bvjtbj

（24）

式中：

vj=e∈Sj∫Ve?NjdV，

Sj，表示圍繞結(jié)點(diǎn)j的單元集合。

對于將設(shè)計(jì)域劃分形成的矩形單元，vj也是常數(shù)。

2.2?求解過程

基于以上推導(dǎo)，可將優(yōu)化迭代過程歸納如下。

（1）采用有限元劃分設(shè)計(jì)域。初始化設(shè)計(jì)參數(shù)，即

t0，bj?=1

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

（25）

式中：下標(biāo)0為迭代指數(shù);J為結(jié)點(diǎn)總數(shù)。

（2）利用有限元分析，得到各結(jié)點(diǎn)位置的主應(yīng)力方向θk，j?和主應(yīng)力方向的應(yīng)變εk，bj?。

（3）采用滿應(yīng)力準(zhǔn)則法優(yōu)化桿件在結(jié)點(diǎn)位置的密度，將桿件方向調(diào)整到主應(yīng)力方向，即

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

（26）

式中：下標(biāo)k為迭代指數(shù);σp為材料的許用應(yīng)力。為避免密度過小導(dǎo)致剛度矩陣，限制

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

（27）

式中：ξ為常數(shù)，此處取1×10-6?。按照優(yōu)化后的桿件密度分布調(diào)整載荷在豎直方向的分布。

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

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

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

（5）將桿件密度在豎直方向的質(zhì)心連線，可以表示最后的拓?fù)浣Y(jié)構(gòu)。

3?算?例

選擇3個(gè)平面問題算例，設(shè)計(jì)域?yàn)?×3，兩端鉸支，載荷只表示力沿水平方向的分布，不表示豎直位置。采用40×30個(gè)正方形單元，由于載荷大小、結(jié)構(gòu)尺寸和材料的彈性模量E等皆與拓?fù)鋬?yōu)化結(jié)果無關(guān)，所以不給出具體數(shù)據(jù)。

算例1選擇文獻(xiàn)[5]中的單工況算例進(jìn)行對比，其力學(xué)模型見圖1a），載荷為均布力。解析解和桿件質(zhì)心連線表示的拓?fù)鋬?yōu)化結(jié)果見圖1b），這兩者非常接近。

算例2選擇在結(jié)構(gòu)跨中和左右1/4跨分別不同時(shí)作用F1和F2，F(xiàn)1=F2/2，其力學(xué)模型見圖2a）。桿件豎向密度質(zhì)心連線表示的拓?fù)鋬?yōu)化結(jié)果見圖2b）。

算例3為在結(jié)構(gòu)左、右半跨分別不同時(shí)作用均布載荷q1和q2，q1=q2，其力學(xué)模型見圖3a）。拓?fù)鋬?yōu)化結(jié)果見圖3b）。

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

4?結(jié)束語

采用類桁架材料模型和有限元方法研究多工況下平面Prager結(jié)構(gòu)的拓?fù)鋬?yōu)化方法。分別計(jì)算各單工況下結(jié)點(diǎn)的剛度矩陣，通過擬合各單工況下材料的最大方向剛度得出多工況下的剛度矩陣，從而得到拓?fù)鋬?yōu)化結(jié)構(gòu)。算例結(jié)果表明該方法可行且有效。在工程實(shí)際中，空間殼體結(jié)構(gòu)的載荷條件和邊界條件遠(yuǎn)比上述討論復(fù)雜，復(fù)雜拓?fù)鋬?yōu)化問題可考慮將以上方法推廣到空間Prager結(jié)構(gòu)的拓?fù)鋬?yōu)化問題，為網(wǎng)殼結(jié)構(gòu)初步設(shè)計(jì)提供參考。

參考文獻(xiàn)：

[1]

沈祖炎， 陳揚(yáng)驥. 網(wǎng)架與網(wǎng)殼[M]. 上海： 同濟(jì)大學(xué)出版社， 1997： 173-178.

[2]周克民， 李俊峰， 李霞. 結(jié)構(gòu)拓?fù)鋬?yōu)化研究方法綜述[J]. 力學(xué)進(jìn)展， 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.

（編輯?武曉英）