基于均匀化理论的拓扑优化算法研究_王书亭

第32卷第10期　　　华　中　科　技　大　学　学　报(自然科学版) 　　　　Vol. 32　No. 102004年　10月　　J. Huazhong Univ. of Sci. &Tech. (Nature Science Edition ) 　　　Oct. 　2004

基于均匀化理论的拓扑优化算法研究

王书亭　左孔天

(华中科技大学机械科学与工程学院, 湖北武汉430074)

摘要:在研究均匀化理论和拓扑优化理论基础上, 推导了复合材料的均匀化求解方程, 并将均匀化理论应用于拓扑优化中, 推导了基于均匀化理论的二维拓扑优化求解算法. 关　键　词:均匀化理论; 拓扑优化; 优化算法

中图分类号:TH122　　文献标识码:A 　　文章编号:167124512(2004) 1020025203

Algorithm of topology optimization theory

W ang S Abstract :The homogenization and analyzed. The equations used for solving the presented. The homogenization theory was intro 2duced to the a 22D topology optimization algorithm based on homogenization theory was K ey w ords :theory ; topology optimization ; optimization algorithm W ang Shuting 　Dr. ; College of Mech. Sci.

430074, China.

&Tech. , Huazhong Univ. of Sci. &Tech. , Wuhan

　　拓扑优化问题是结构优化中相对复杂的一类

问题, 也是结构优化领域研究的热点[1, 2]. 本文在研究均匀化理论和拓扑优化理论基础上, 推导了二维拓扑优化问题在直角坐标系下的求解公式.

1　复合材料的均匀化理论

1. 1　均匀化方法的引入

图1　微观胞元结构

ε

为宏观度量. 微元体受体力f εi , 面力p i , 外部作用力t i 作用. 微元体在平衡条件下满足虚功方程

复合材料模型可看作由不同材料组成的微结构胞元在空间的组合与重复, 其微观胞元计算域

Y 的结构如图1所示. 这种带空洞的胞元在空间

的重复造成了材料宏观特性在空间的各向异性.

在均匀化方法中, 通过在材料中引入带空洞的胞元模型来求解材料的宏观特性参数[3～5], 从工程角度看, 本文采用的均匀化方法是寻求一种用宏观量来表达材料微观量变化特征的方法. 1. 2　复合材料问题中的均匀化求解方法

设y 为复合材料微观胞元度量, y =x/ε,

x

收稿日期:2004203217.

∫E

f v d Ω+t v d Γ+p v d S. ∫∫∫

　

εΩ

ε

ijkl

εd Ω=

x l x j

Γt

i i

εε

Ω

εi i

S

εi i

(1)

εεε

应力应变关系为σij =E ijkl e kl , 应变位移关系为εεε

e kl =(5u k /5x l +5u l /5x k ) /2. 由弹性常量对称性εεεεε有E ijkl =E jikl =E ijlk =E klij . u 依赖于宏观度量x

和微观度量y =x/ε, 因此可将u ε展开为

1

　u ε(x ) =u 0(x , y ) +εu (x , y ) +

作者简介:王书亭(19682) , 男, 博士; 武汉, 华中科技大学机械科学与工程学院(430074) .

E 2m ail :wangst @hust.edu. cn

基金项目:国家航天部载人航天计划“921”工程重点子项目; 国家高技术研究发展计划资助项目(2003AA001031) .

　　　　　　　　　　　　　　华　中　科　技　大　学　学　报(自然科学版) 　　　　　　　　　　　第32卷26

22

εu (x , y ) +…,

012

u , u , u , …可用Φ统一表达为:

(Φ(x , y ) ) =+.

x i x i εy i

(2)

0i

d Y d Ω=

y l x p v d S

|Y |S

i i

d Ω.

(3)

　　对于一个Y 周期函数Ψ(y ) 有如下两式存在:

Ψd Ω=Ψ(y ) d Y d Ω;

lim ε

+ε|Y |ΩγΩε→0

　　由于上式对任意v 均满足, 因此可以选择

v =v (y ) , :

0() 1()

　E ijkl +d Y =

5x l 5y l 5y j γ

∫

S

∫

ΩS

∫p v (y ) d S. ∫

S

i i

(9)

(4)

ε→0

lim +εε

∫

ε

d S =

|Y |

Ψ(y ) d S d Ω,

∫

(5)

式中|Y |代表胞元体积.

将式(2) 和(3) 代入式(1) , 得

01ε

　εE ijkl 2++・

y l Ωεy l y i εx l

式(9) 即为微观平衡方程.

同样地, 由式(8) , 并利用式(4) 和(5) 的性质, 并使v =v (x ) 时, 可得宏观平衡方程

01　E ijkl +d Y ・

x l l Ω|Y |γ

∫

0++y i y l x 01+x l l j

() =

|Y y

t (x ) j

i i

γ

i

v i (x ) d Ω+

(10)

Γt

12d ++(x l y f v d Ω+t v d Γ+p v d S. ∫∫

ε

ε

Ω

εi

(9) , 该等式为u 0k 与p i 的线性关系

式. 由下两式可求出特征变形量χkl 和ψk

∫

γ

E ijpm

kl ()

d Y =

y m y j

∫

γ

S

E ijkl

()

d Y ; y j

(11)

i

Γt

i i

S

εi i

此式按ε的权值可分解为以下3式:

　εE

ε

2

Ω

ε

E ijkl

ε

ijkl

ε

d Ω=0;

y l y i

∫

γ

E ijkl

()

d Y =

y l y j

p v (y ) d Y. (12) ∫

i i

(6)

Ω

ε

01++x l y l y j

由于式(9) 的线性关系, u 1可表示为

　u 1i =-χkl i (x , y ) (5u k (x ) /5x i ) -ψi (x , y ) + u i (x ) ,

u 代表材料微观力学特性; u i (x ) 为余项. 将上式

1

d Ω=

y l x 0

∫

S

εi

p v i d S ;

ε

(7)

带入式(10) 可得　

　

∫

Ω

ε

E ijkl

ε

011+++x l y l x j x l

|Y |E

Ω

ijkl

-E ijpm

kl d Y ・y m Y |

Ω

2d Ω=

y l y ∫

Ω

εf i v i d Ω+

ε

∫

Γt

t i v i d Γ. (8)

0(x ) (x )

d Ω=

x l x j

|

Γt

i i

・i

2

　　将式(6) 乘以ε并利用式(4) 性质有

E ijkl d Y d Ω=0, 并进行分部积

|Y |ΩY y l y i

∫

|Y |

ijkl

∫

γ

E ijkl

()

d Y d Ω+y l x j

v i (x ) d Ω+

f d Y

|Y |γ

・

分有:

　

t v (x ) d Γ. ∫

Ω

Y

0-E ijkl

5y i 5y l

令

v i d Y +

H

E ijkl

(x ) =

|Y |

E ∫

s

0n v d s d Ω=0. 5y l j i

E

ijkl

-E ijpm

kl d Y ; y m

(13)

由于v (y ) 为任意, 由上式得-(5/5y i ) (E ijkl ・5u 0k /5y l ) =0, E ijkl (5u 0k /5y l ) n j =0, 进而有u 0=

ε00

u (x ) . 即u 扩展的第一项u 是x 的函数. 由式

τij (x ) =

|Y ||

b i (x ) =

f d Y , Y |γ

E ijkl

d Y ; y l

(14) (15)

γ

i

(7) , (4) 和(5) 　

|Y |E

ijkl

01++5x l 5y l 5y j

则方程(10) 可表示为

0() () H

　E ijkl (x ) d Ω=

5x l 5x j Ω

∫

第10期　　　　　　　　　　王书亭等:基于均匀化理论的拓扑优化算法研究　　　　　　　　　　　　　27

∫∫t (x ) v (x ) d Γ. ∫

Ω

τij (x )

i

(x )

d Ω+x j

i

Ω

b i (x ) v i (x ) d Ω+

(16)

Γt

只需求解以下三种情况:k =l =1; k =l =2; k =1, l =2. a . 第一种情况, k =l =1. 扩展式(11) 有

1111　E 1111+E 1122+

y 1y 2y 1γ

E 1122

　　综上所述, 复合材料均匀化求解方法可归结

为三个非耦合问题:其中两个是微观层问题, 即式(13) 中的均匀化弹性张量的求解和式(14) 中残余应力的求解. 另一个是宏观层的问题, 即式(16) 中的宏观位移场的求解.

E 1212

1111+5y 25y 1++5y 25y 1

1111+E 2222d Y =5y 15y 25y +E 1122d Y. y 1y 2

|Y |

(17)

2　基于均匀化理论的二维拓扑优化

E 1111

算法

2. 1　线弹性结构拓扑优化问题的表达形式

i =j =1时扩展式(13) 有

H

(x ) =　E 1111

・

为了计算方便, 在拓扑优化过程中,

将连续体

离散为带长方形孔洞的正方形单元结构, 如图2所示, 单元外边长为1, 长方形孔洞边长为a , b ,

E

-1111-d Y. (18) y 2

T

11

1

, , +5y 15y 25y 25y 1, , +y 1y 2y 2y 1

D 11

D 12D 22

d 2

d =

D 12

;

T

;

00

D .

D =d 1

图2　带方形孔洞的单元结构

则式(17) , (18) 成为:

其取值范围为a , b ∈[0, 1].用最小柔度表达为

(L (u ) ) ; minimize e e e

a , b , θ

∫

Y H

εT (v ) D ε(

|Y |

∫

Y

εT (v ) d 1d Y ;

D 11(x ) =

s. t :a (u , v ) =L (v ) ,

(1-∫

Ω

γ

T

(D 11-d 1ε(

a b ) d Ω≤Ωm .

e e

　　定义结构柔度为

L (v ) =

这样, 通过在计算域中离散, 可求得

H

进一步求得D 11.

b . 第二种情况, k =l =2. 与前面推导类似, 扩展式(11) , 并令i =j =1和i =j =2扩展式

22

(13) , 且令χ221=

∫

Ω

f v d Ω+

∫

ΓT

tv d s ,

e

式中a e , b e , θ为单元设计变量, e =1, 2, …, N , N 为有限单元的数目. 2. 2　二维拓扑优化中均匀化弹性张量的求解

∫

Y H

εT (v ) D ε(

|Y ||

∫

Y

εT (v ) d 2d Y ;

D 12(x ) =

H

D 22(x ) =

将均匀化方法应用于拓扑优化中时, 结构胞元中一般不存在内力边界. 当不考虑体力项的影响时, 宏观平衡方程(16) 可写为

0() () H

　E ijkl (x ) d Ω=

x l x j Ω

(D Y |γγ

T

(D 12-d 1ε(

T

-d 2ε(

22

∫

t (x ) v (x ) d Γ. ∫

Γt

i

i

H H

在计算域中离散, 可求得

c . 第三种情况, k =1; l =2. 同理, 扩展式

12

(11) 并针对i =1, j =2时扩展式(13) , 且令χ1=

　　微观特征变形χkl 是Y 周期的变量, 用于求

H

解均匀化弹性张量E ijkl , 当k 和l 分别取不同值时得到不同的求解方程. 二维线弹性问题时, i , j ,

H

k , l , p , m =1, 2. 由于χkl 和E ijkl 的对称性, 因此

∫

Y

εT (v ) D ε(

|Y |

∫

Y

εT (v ) d 3d Y ;

H

D 66(x ) =

γ

T

(D 66-d 3ε(

(下转第30页)

　　　　　　　　　　　　　　华　中　科　技　大　学　学　报(自然科学版) 　　　　　　　　　　　第32卷30

　　二维数值算例表明:用均匀化理论来进行结

构拓扑优化是切实可行的, 用优化准则法对拓扑优化中的设计变量更新能够得到较好的优化结果, 计算具有较好的数值收敛性和稳定性.

参

考

文

献

[1]王书亭, 左孔天. 基于均匀化理论的拓扑优化算法研

究. 华中科技大学学报(自然科学版) , 2004, 32(10) :

图1　算例模型

25～27,30

[2]Guedes J M , K ikuchi N. Preprocessing and postpro 2

s. t. F =KU , Ω/Ω0=0. 5; 0≤a e ≤1; 0≤b e ≤1.

(8)

cessing for materials based on the homogenization method with adaptive finite element method. Comp Methods Appl Mech Engrg , , 83:143～198[3]Behrooz Hassani , Ernest Homogenization and

2Verlag Lon 2, , G in. The COC algorithm , part Ⅱ:

, geometry and generalized shape optimiza 2tion. Computer Methods in Applied Mechanics and En 2gineering , 1991, 89:197～224

经过115次迭代后目标函数的迭代收敛曲线如图

2所示

.

图(上接第27页)

华中科技大学出版社,1997.

[2]程耿东, 顾元宪, 王　健. 我国机械优化研究与应用的

从而在计算域中离散, 可求得

D 66.

H

综述和展望. 机械强度, 1995, 17:68～74

d . 考虑角度变化时的情况. 当考虑单元角度

[3]刘书田, 程耿东. 复合材料应力分析的均匀化方法. 力

变化时, 均匀化弹性张量可进一步表示为H e e e T e H e e e D (a , b , c ) =R (θ) D (a , b ) R (θ) ,

e

) =　R (θ

e cos 2θe

sin 2θe e θ-2cos sin θ

e

sin 2θe

cos 2θe e θ2cos sin θ

e e

cos θsin θe e

θ-cos sin θe e

cos 2θ-sin 2

学学报,1997,29:306～313

[4]Martin Philip Bendsoe , Noboru K ikuchi. G enerating

optimal topologies in structural design using a homoge 2nization method. Computer Methods in Applied Me 2chanics and Engineering , 1988,71:197～224

[5]Suzuki K , K ikuchi N. A homogenization method for

shape and topology optimization. Computer Methods in Applied Mechanics and Engineering , 1991, 93:291～381

.

　　由上述方法求得单元均匀化弹性张量后, 进一步可求得结构整体位移场和应力场, 完成二维

拓扑优化问题的一个迭代步计算.

参

考

文

献

[1]余　俊, 周　济. 优化方法程序库OPB 22原理及应用.