力学学报  2018 , 50 (2): 385-394 https://doi.org/10.6052/0459-1879-17-286

生物、工程及交叉力学

基于改进的双向渐进结构优化法的应力约束拓扑优化

王选1, 刘宏亮1, 龙凯2, 杨迪雄1*, 胡平1

1 大连理工大学 工业装备结构分析国家重点实验室,大连 116023
2 华北电力大学 新能源电力系统国家重点实验室, 北京 102206

STRESS-CONSTRAINED TOPOLOGY OPTIMIZATION BASED ON IMPROVED BI-DIRECTIONAL EVOLUTIONARY OPTIMIZATION METHOD

Wang Xuan1, Liu Hongliang1, Long Kai2, Yang Dixiong1*, Hu Ping1

1 State Key Laboratory of Structural Analysis for Industrial Equipment, Dalian University of Technology, Dalian 116023, China
2 State Key Laboratory for Alternate Electrical Power System with Renewable Energy Sources, North China Electric Power University,Beijing 102206, China

中图分类号:  O342,O346

通讯作者:  *通讯作者:杨迪雄,教授,主要研究方向:结构优化与建筑抗震减震. E-mail:yangdx@dlut.edu.cn*通讯作者:杨迪雄,教授,主要研究方向:结构优化与建筑抗震减震. E-mail:yangdx@dlut.edu.cn

收稿日期: 2017-08-22

接受日期:  2018-03-1

网络出版日期:  2018-02-10

版权声明:  2018 《力学学报》编辑部 《力学学报》编辑部 所有

基金资助:  国家自然科学基金(11272075, 11772079), 北京市自然科学基金(2182067)和中央高校基本科研业务费专项(2017MS077, 2018ZD09)资助项目.

展开

摘要

工程结构设计时经常需要限制最大名义应力,以避免发生断裂或疲劳破坏,一个有效的策略是采用拓扑优化方法. 常规的双向渐进结构优化法(bi-evolutionary structural optimization, BESO)不能有效求解应力约束拓扑优化问题,为此本文提出一种改进的双向渐进结构优化方法,处理体积和应力约束下的最小柔顺性问题. 引入基于K-S函数的全局应力度量,以减小大量局部应力约束引起的计算代价. 采用拉格朗日乘子法将应力约束函数引入到目标函数,然后由二分法确定合适的拉格朗日乘子的值使得应力约束得到满足. 而且,详细推导了基于BESO方法的应力约束拓扑优化模型及其灵敏度列式,最后通过三个典型拓扑优化算例验证改进方法的有效性. 为展示考虑应力约束的优点,将应力约束设计与传统的基于刚度的设计进行了比较. 结果表明, 改进的BESO方法优化迭代过程稳健,获得了边界灰度单元很少的清晰的拓扑构型,并实现了有效降低应力集中效应的设计.

关键词: 结构拓扑优化 ; 应力约束 ; 双向渐进结构优化法 ; 拉格朗日乘子法 ; 二分法

Abstract

It is necessary to limit maximum nominal stress for engineering structural design generally, so as to avoid that the failure of fracture or fatigue occurs. Topology optimization approach is a feasible strategy. The conventional bi-evolutionary structural optimization (BESO) method cannot effectively address the topology optimization problem with stress constraint. To overcome this limitation, this paper suggests an improved BESO method for stress-constrained topology optimization, in which the minimum compliance design problem with volume and stress constraints is considered. A global stress measure based on the K-S function is introduced to reduce the computational cost associated with the local stress constraint. Meanwhile, the stress constraint function is added to the objective function by using the Lagrange multiplier method. Moreover, the appropriate value of the Lagrangian multiplier is then determined by a bisection method so that the stress constraint is satisfied. The model of BESO method for solving stress-constrained topology optimization and its sensitivity analysis are detailed. Finally, three typical topology optimization examples are performed to demonstrate the validity of the present method, in which the stress constrained designs are compared with the traditional stiffness based designs to illustrate the merit of considering stress constraint. The optimized results indicate that the improved BESO method, as a robust algorithm with stable iterative history, achieves an ambiguous topology with clearly defined boundaries, and realizes a design that effectively reduces stress concentration effect at the critical stress areas.

Keywords: structural topology optimization ; stress constraint ; BESO method ; Lagrange multiplier method ; bisection method

0

PDF (5742KB) 元数据 多维度评价 相关文章 收藏文章

本文引用格式 导出 EndNote Ris Bibtex

王选, 刘宏亮, 龙凯, 杨迪雄, 胡平. 基于改进的双向渐进结构优化法的应力约束拓扑优化[J]. 力学学报, 2018, 50(2): 385-394 https://doi.org/10.6052/0459-1879-17-286

Wang Xuan, Liu Hongliang, Long Kai, Yang Dixiong, Hu Ping. STRESS-CONSTRAINED TOPOLOGY OPTIMIZATION BASED ON IMPROVED BI-DIRECTIONAL EVOLUTIONARY OPTIMIZATION METHOD[J]. Chinese Journal of Theoretical and Applied Mechanics, 2018, 50(2): 385-394 https://doi.org/10.6052/0459-1879-17-286

引 言

近20年来连续体结构拓扑优化获得了快速发展,成为工程创新设计的强有力工具. 各种拓扑优化方法,例如基于材料分布的均匀化方法[1]、密度法/各向同性实体材料惩罚法(solid isotropic material with penalization, SIMP)[2]、渐进结构优化法(evolutionary structural optimization, ESO)[3]、独立连续映射法(independent continuous mapping method, ICM)[4],基于几何边界描述的水平集方法(level set method, LSM) [5],以及最近发展的移动可变形组件法(moving morphable components, MMC) [6],都被用来解决各种结构拓扑优化设计问题[7,8,9,10,11,12,13,14].

由应力集中、或高应力值所引起的结构断裂、疲劳破坏严重影响其使用寿命,因此研究应力约束拓扑优化问题具有重要意义. 与体积约束下的刚度最大化模型相比,考虑应力约束的拓扑优化有其自身的难点. 第一个难点是所谓的“奇异解现象” [15,16]. 奇异解现象起源于桁架优化问题,Cheng和Jiang对拓扑优化问题中的奇异解给出了数学解释[15],即它主要是由应力约束的不连续导致的. Duysinx和Bendsoe[17]指出在应力约束的连续体拓扑优化中也存在奇异解现象. 目前常用的处理奇异解现象的方法包括 ε松弛技术[16,17,18]qp松弛技术[19]等.

应力约束拓扑优化设计的第二个难点在于应力作为一个局部的物理量所引起的大量局部性约束. 一般情况下,为了较为准确 地计算结构的应力场,需要加密有限元网格,导致局部应力约束的数目增加,从而显著增加了局部应力约束的敏度分析的计算 代价[20,21]. 为了处理这个问题,通常采用逼近最大局部函数值的凝聚函数将局部应力约束转化为一个全局的应力约束[18,20-26]或者“分块”的应力约 束[21,23,27]. 目前较为常用的凝聚函数有P 范数[20,21,22,23]和Kieisselmeier-Steinhauser (K-S)函数[18,20,24-28],其中Yang和Chen[24]早在1996年采用K-S函数将诸多应力约束化为一个应力约束. 另外一种有效方式是隋允康等利用von Mises强度理论,借助应变能将应力约束全局化,从而显著减少了应力约束个数[4].

应力约束拓扑优化问题的第三个难点在于应力的高度非线性行为. 临近区域密度值的改变对结构某些关键区域的应力水平和非线性行为产生严重影响,特别是具有较大的空间应力梯度区域,如有凹陷角、拐角处[18,23,29]. 这就要求结构拓扑优化列式具备有效降低或者消去应力集中现象的能力,而且需要求解算法与优化列式保持数值一致性以确保优化问题的稳定收敛[23,29].

大部分考虑应力相关的拓扑优化问题都是基于SIMP法[17-27,30]和水平集方法[29,31-33]. 此外,Cai等[28,34]将水平集函数与有限胞元法结合,提出了一种求解应力约束拓扑优化问题的新框架. 隋允康 等[4,35-37]基于ICM法提出了一系列有效求解应力约束拓扑优化问题的策略,更多关于应力约束和其他力学性能约束的处 理方法可以参照文献[38]. 荣见华等[39]将ICM法与ESO方法结合,通过每步迭代不断更新位移和应力约束,提出了一种新的应力优化思路.

调查发现,常规的ESO/BESO法[8,40]研究应力约束的拓扑优化问题鲜有报道,这是因为对于考虑应力约束的优化问题而 言,结构的应力值对设计变量的变化特别敏感,从而常规的仅包含0 ~1离散变量的ESO/BESO法求解应力约束拓扑优化问题 时,会导致应力约束函数发生极度振荡现象,无法有效控制结构的应力水平. 为此,本文提出一种改进的双向渐进结构优化方法,使设计变量变化量 Δxe随着迭代步数的增加而逐渐减小进而稳定优化过程. 采用 K-S 函数来凝聚局部应力约束以减少与局部应力约束相关的计算代价, 并通过拉格朗日乘子法施加应力约束. 且详细推导了基于BESO 方法的应力约束拓扑优化模型及其灵敏度列式,最后通过三个典型拓扑优化算例验证了本文方法的有效性.

1 应力约束拓扑优化列式

1.1 拓扑优化模型

结构拓扑优化的目的是寻找一个满足某些约束条件的最优的材料分布,以使结构获得某种最优的结构行为,如结构的重量最轻或刚度最大化. 本文与文献[19,28,33]一样,考虑应力和体积约束下的柔顺性最小化问题,其对应的数学列式可表示如下

Find: x

Min: $C = {\pmb F}^{\rm T}{\pmb U}$

$S.t.: \left\{ \begin{array}{l} {\pmb K}{\pmb U} ={\pmb F} \\ \sum_{e = 1}^{N_e } V_e x_e = V_0 f_{\rm v} \\ \sigma_{\max}^{\rm vM} = {\max } \left( \sigma _e^{\rm vM}\right) \leqslant \sigma ^\ast , \ \ e = 1,2, \cdots, N_e \end{array} \right. \ \ (1)$其中, ${\pmb x} = \left[ {x_1 ,x_2 , \cdots ,x_{N_e } } \right]$为单元设计变量的集合,$N_e$为总的设计变量的数目;$C$表示结构的柔顺性目标函数; ${\pmb K}$为结构的整体刚度阵,由单元刚度阵${\pmb K}_e$组装形成;${\pmb U},{\pmb F}$分别为结构的位移向量和载荷列阵;$V_e ,V_0$分别为第$e$个单元的体积和整个设计域的体积,$f_v $为准许的材料体积分数;$\sigma _e^{\rm vM} $为每个单元的vonMises应力,$\sigma _{\max}^{\rm vM} ,\sigma^\ast $为结构的最大von Mises应力及其约束限值.

1.2 插值格式

如上文所述,与柔顺性最小化问题不同的是,应力约束拓扑优化问题存在其自身的难点,如奇异最优解现象、大量局部应力约束导致的巨大计算量及应力约束的高度非线性问题等.为了克服应力约束的奇异性问题,获得清晰的拓扑构型,这里采用文献[23]的方法,对结构刚度和应力采用不同的幂指数进行惩罚.弹性模量、应力与单元密度之间的关系可定义为

$ E_e^\ast = x_e^{p_1 } E_0\,, \ \ {\pmb\sigma}_e = x_e^{p_2 } {\pmb\sigma}_e^0 (2) $

式中, ${\pmb\sigma}_e^0 ={\pmb D}_0{\pmb B \pmb u}_e = [\sigma _{e,x} ,\sigma _{e,y} ,\tau _{e,xy} ]^{ T}$是在 第$e$个单元中心点处计算得到的包含3个应力分量的应力列阵;$p_1 ,p_2$分别为弹性模量和应力的惩罚参数,本文取$p_1 =3$, $p_2= 0.5$;$E_0 , {\pmb D}_0$分别为实体材料的弹性模量、弹性矩阵,${\pmb B}$为应变位移矩阵. 式(1)中每个单元的von Mises应力可以通过下式计算

$ \sigma _e^{\rm vM} = \left( {\sigma _{e,x}^2 + \sigma _{e,y}^2 - \sigma _{e,x} \sigma _{e,y} + 3\tau _{e,xy}^2 } \right)^{1 / 2} (3) $

1.3 全局应力度量

为了减少局部应力约束导致的计算代价,通常采用凝聚函数来形成一个全局的应力约束. 较为常用的凝聚函数有P 范数和K-S函数. 不失一般性,本文使用K-S函数作为凝聚函数

$ \sigma^{\rm KS}=\dfrac 1{\mu} \ln \Big [ \sum^{N_e}_{e=1} \exp \Big ( \mu \dfrac{\sigma^{\rm vM}_e}{\sigma^*}\Big ) \Big ] (4) $

理论上当参数$\mu $的值趋向于无穷大时,上式中的全局函数$\sigma^{\rm KS}$趋向于${\sigma _{\max}^{\rm vM} }/{\sigma ^{\ast }}$,因此,式(1)中的全局应力约束变为

$ \sigma ^{\rm KS} \leqslant 1 (5) $

需要注意的是,由于式(4)中参数$\mu $不可能取无穷大,这样会导致全局函数$\sigma^{\rm KS}$不能很好地逼近${\sigma _{\max}^{\rm vM} } / {\sigma ^{\ast }}$,为此, 利用式(6)来修正式(5)中的约束

$ \bar {\sigma }^{\rm KS} = c\sigma ^{\rm KS} \leqslant 1 (6) $

其中,修正参数 $c = {\sigma _{\max}^{\rm vM} }/ {\left( {\sigma ^{\ast } \cdot \sigma ^{\rm KS}} \right)}$.

2 灵敏度分析

2.1 目标函数和应力约束的敏度

利用伴随法可推导优化列式目标函数的敏度为

$ \dfrac{\partial C}{\partial x_e } = - {\pmb u}_e^{\rm T} \dfrac{\partial {\pmb K}_e }{\partial x_e }{\pmb u}_e (7) $

根据链式法则很容易得到式(4)中$\sigma ^{\rm KS}$关于设计变量的敏度为

$ \dfrac{\partial \sigma^{\rm KS} }{\partial x_e } = \dfrac{\partial \sigma^{\rm KS} }{\partial \sigma _e^{\rm vM} }\left( {\dfrac{\partial \sigma _e^{\rm vM} }{\partial {\pmb\sigma}_e }} \right)^{\rm T}\dfrac{\partial {\pmb\sigma}_e }{\partial x_e } (8) $

从上式可以看出,要计算全局应力约束函数的敏度必须先确定K-S函数对von Mises应力的导数、von Mises应力对应力分量的导数、应力分量对设计变量的导数,下面分别计算这三项.

2.2 K-S函数对von Mises应力的导数

对式(4)中的K-S函数关于每个单元的von Mises应力求导可得

$ \dfrac{\partial \sigma^{\rm KS} }{\partial \sigma _e^{\rm vM} } = \dfrac{\exp\left({\mu \dfrac{\sigma _e^{\rm vM} }{\sigma^\ast }} \right)}{\sigma^\ast \left( {\sum_{e = 1}^{N_e } {\exp\left({\mu \dfrac{\sigma _e^{\rm vM} }{\sigma^\ast }} \right)} } \right)} (9) $

2.3 von Mises应力对应力分量的导数

式(3)中von Mises应力关于3个应力分量的导数可表示为

$ \left.\begin{array}{l} \dfrac{\partial \sigma _e^{\rm vM} }{\partial \sigma _{e,x} } = \dfrac{1}{2\sigma _e^{\rm vM} }\left( {2\sigma _{e,x} - \sigma _{e,y} } \right) \\ \dfrac{\partial \sigma _e^{\rm vM} }{\partial \sigma _{e,y} } = \dfrac{1}{2\sigma _e^{\rm vM} }\left( {2\sigma _{e,y} - \sigma _{e,x} } \right) \\ \dfrac{\partial \sigma _e^{\rm vM} }{\partial \tau _{e,xy} } = \dfrac{3\tau _{e,xy} }{\sigma _e^{\rm vM} } \end{array} \right\} (10) $

2.4 应力分量对设计变量的导数

根据式(2)可推导应力分量关于设计变量的导数

$ \dfrac{\partial {\pmb \sigma}_e }{\partial x_e } = p_2 x_e^{\left( {p_2 - 1} \right)} {\pmb D}_0 {\pmb B}{\pmb u}_e + x_e^{p_2 } {\pmb D}_0 {\pmb B}\dfrac{\partial {\pmb u}_e }{\partial x_e } (11) $

将式(9)~式(11)代入式(8),再结合式(6),可得

$ \dfrac{\partial \bar {\sigma }^{\rm KS} }{\partial x_e } = C \dfrac{\partial \sigma^{\rm KS} }{\partial \sigma _e^{\rm vM} }\left( {\dfrac{\partial {\pmb\sigma}_e^{\rm vM} }{\partial {\pmb\sigma}_e }} \right)^{\rm T} \cdot \Big( p_2 x_e^{\left( {p_2 - 1} \right)} {\pmb D}_0 {\pmb B}{\pmb u}_e + x_e^{p_2 } {\pmb D}_0 {\pmb B}\dfrac{\partial{\pmb u}_e }{\partial x_e } \Big) (12) $

3 拉格朗日乘子法及优化算法

3.1 拉格朗日乘子的确定

在渐进结构优化方法中,由于结构体积作为进化参数所以体积约束很容易满足. 其他的约束均通过拉格朗日乘子法将约束函数引入到目标函数中,因此,修改的目标函数可表示为
$ f = {\pmb U}^{\rm T}{\pmb K}{\pmb U }+ \lambda \left( {\bar {\sigma }^{\rm KS} - 1} \right) (13) $

其中$\lambda $为拉格朗日乘子. 可以看出:如果$\bar {\sigma }^{\rm KS} =1 $,则修改的目标函数等于原来的目标函数;否则,如果$\bar {\sigma }^{\rm KS}\leqslant 1$,取$\lambda = 0$,此时说明应力约束已经满足了;如果$\bar {\sigma }^{\rm KS} > 1$,$\lambda $趋向于无穷大,这说明需要在随后的迭代中减小结构的应力来满足应力约束.

为了计算式(13)中修改的目标函数关于设计变量的灵敏度,需要先确定拉格朗日乘子的值. 为了在程序中方便地实施计算,可将拉格朗日乘子重新定义为

$ \lambda = \dfrac{1}{\omega } - 1 (14) $

其中参数$\omega $的取值范围为[10-10,1]. 确定合适的参数$\omega$将使得应力约束条件满足. 因此,根据当前的应力值和敏度值估算下一次迭代的应力值

$ \bar {\sigma }_{k + 1}^{\rm KS} \approx \bar {\sigma }_k^{\rm KS} + \sum_{e = 1}^{N_e } {\dfrac{\partial \bar {\sigma }_k^{\rm KS} }{\partial x_e }} \Delta x_e (15)$

其中$k$为当前的迭代数. 根据下一次迭代的应力值及应力约束的限值采用二分法确定参数$\omega $,具体实施过程详见文献[41]. 一旦获得参数$\omega$的值,即可通过式(14)确定拉格朗日乘子的值,进而确定修改的目标函数关于设计变量的敏度

$ \dfrac{\partial f}{\partial x_e } = - {\pmb u}_e^{\rm T} \dfrac{\partial {\pmb K}_e }{\partial x_e }{\pmb u}_e + \lambda \dfrac{\partial \bar {\sigma }^{\rm KS} }{\partial x_e } (16) $

在渐进结构优化法中,设计变量的变化是根据敏度数的相对排序确定的[41,42,43],第$e$个单元所对应的敏度数可定义为

$ \alpha _e = - \dfrac{\partial f}{\partial x_e } (17) $

3.2 数值处理

为避免棋盘格和网格依赖性等数值不稳定问题,采用式(18)对式(17)中定义的单元敏度数过滤

$ \tilde {\alpha }_e^ = \Bigg (\sum_{i = 1}^{Nr} {w_{e,i} \alpha _i } \Bigg)\Bigg / \,\Bigg(\sum_{i = 1}^{Nr} {w_{e,i} } \Bigg) (18)$

其中权重$w_{e,i} $定义为

(19)

式中$r_{e,i} $为单元$e$和单元$i$之间的距离,$N_r$是以单元$e$为中心、半径为$r_{\min} $的圆形邻域内的单元个数.

另外,数值研究表明,考虑敏度数的历史信息可以有效地改善优化过程的稳定性,即采用前后两轮迭代的敏度数的均值作为当前敏度数的修正值来参与迭代,其表达式为

$ \hat {\alpha }_e^k = \dfrac{1}{2}(\tilde {\alpha }_e^k + \tilde {\alpha }_e^{k - 1} ) (20)$

BESO方法从完整的设计域开始,在优化过程中结构的体积不断地减小直到获得目标体积. 在每一步迭代中均以$\Delta x_e $的步长来更新设计变量,实现设计变量的增加或减少. 在常规的BESO方法[41,42]中,$\Delta x_e =1$意味着设计变量只能取0或1两个值,在文献[43]中,取$\Delta x_e = 0.02$,这意味着设计变量除了取0$\sim$1值之外,还可以取其他一些中间值. 本文考虑应力约束的结构拓扑优化问题,在研究中发现最大应力值$\bar {\sigma }^{\rm KS} $对设计变量的变化特别敏感,因此通过使设计变量变化量$\Delta x_e $随着迭代步数的增加而逐渐减小的策略来稳定优化过程.

3.3 优化算法

改进的渐进结构优化方法求解应力约束拓扑优化问题的算法流程可归纳如下.

(1) 定义构型参数:定义容许体积分数$f_{\rm v}$、体积进化率$ER$及应力约束限值$\sigma^\ast $等构型参数.

(2) 有限元分析:求解位移场和应力场,进而根据式(4)计算$\sigma ^{\rm KS}$的值.

(3) 材料体积更新:通过$V_{k + 1} = V_k (1 -ER)$决定下一个迭代步的材料体积,其中为材料体积进化率. 一旦到达目标体积分数$f_{\rm v} $,后续迭代中保持目标体积分数不变 .

(4) 灵敏度分析:根据式(7)和式(12)分别计算目标函数和应力约束函数的敏度,在此基础上根据3.1节的算法确定拉格朗日乘子的值.

(5) 敏度过滤及稳定性处理:根据式(17)定义单元敏度数,然后对敏度数进行过滤,再根据式(21)修正单元敏度数.

(6) 设计变量更新:通过下式来更新设计变量

$ x_e^{k + 1} = \left\{ \!\!\begin{array}{ll} {\min\left( {x_e^k + \Delta x_e \,, \ 1} \right)} \,, & \hat {\alpha }_e^k \geqslant \alpha _{\rm th} \max\left( {x_e^k - \Delta x_e\,, \ x_{\min} } \right) \,, & \hat {\alpha }_e^k < \alpha _{\rm th} \right.(21) $

其中$\Delta x_e = \max\left( {\zeta \Delta x_e ,10^{ - 3}} \right)$,$\zeta$是小于1的正常数,其目的是使设计变量变化量$\Delta x_e$随着迭代步数的增加而逐渐减小. $\Delta x_e$初值的选取不宜太大或太小,初值太大会导致目标函数和应力约束函数震荡,太小则不利于收敛. $x_{\min} $取很小的正数,作为设计变量下限以防止刚度矩阵奇异,这里取$x_{\min}= 10^{ - 4}$. $\alpha _{\rm th} $为敏度数阈值,可通过下一次迭代的目标体积及单元敏度数的排序确定,具体过程见文献[42].

(7) 收敛性判断:重复式(2) ~式(6),直到满足下列收敛条件之一,优化过程停止迭代

$ {\Big | \sum_{i = 1}^5 {(C_{k - i + 1} } - C_{k - 5 - i + 1} )\Big|}\Bigg / {\sum_{i = 1}^5 {C_{k - i + 1} } } \leqslant \mu _1 (22)$

$ {\Big | \sum_{i = 1}^5 (\sigma _{\max, k - i + 1}^{\rm vM} - \sigma _{\max, k - 5 - i + 1}^{\rm vM} )\Big | } \Bigg / \sum_{i = 1}^5 \sigma _{\max, k - i + 1}^{\rm vM} \leqslant \mu _2 (23)$

式中,$k$为当前的迭代步数,$\mu _1 ,\mu _2 $分别为柔顺性和最大von Mises应力的容许收敛误差,本文取$\mu _1 = \mu _2= 1\times 10^{ - 6}$.

4 数值算例与讨论

本节通过3个平面应力问题拓扑优化算例来展示本文改进的BESO方法处理应力约束拓扑优化问题的有效性.除非特殊声明,所有算例中变量和几何参数均使用无量纲参数,完全实体材料的杨氏模量$E_0$和泊松比分别设置为1.0和0.3,平面应力单元的厚度设为1.所有设计变量的初始值取$x_e = 1$,取参数$\mu = 10$,$\zeta = 0.99$,体积进化率$ER =0.01$,设计变量改变量$\Delta x_e = 0.04$.

4.1 算例1

考虑经典的L形梁优化问题[20,23-26,28-31],几何区域和边界条件如图1所示,结构左边顶部固支,右上角A点受载荷$F = 10$作用,为了避免载荷作用点处的应力集中,将集中载荷平均分配到如图 1所示的邻近的5个节点上. 初始结构在载荷$F$作用下的最大的von Mises应力为7.770 6. 在本算例中,设置材料体积分数和应力约束限值分别为$f_v=0.4$, $\sigma^{\ast } = 6.0$. 结构由6400个四节点四边形单元来离散,过滤半径取$r_{\min} = 4$.

图1   L形梁设计域和边界条件

Fig. 1   Design domain and boundary condition for L-shaped beam

图2图3展示了本文方法获得的优化结果,其中图2 是仅考虑体积约束的刚度最大化优化结果,图3是同时考虑了体积和应力约束的优化结果.

图2可以看出,不考虑应力约束的L形梁在拐角处出现了明显的应力集中现象,而在考虑应力约束的情况下,L形梁的拐角处 出现了带有弧度的圆角(图3),这说明本文方法是非常有效的,可以获得边界灰度单元很少的清晰的拓扑构型,产生有效降 低应力集中效应的设计.

图2   不考虑应力约束L形梁的拓扑优化设计

Fig. 2   Optimal topology designs for L-shaped beam without stress constraint

图3   应力和体积约束下L形梁的拓扑优化设计

Fig. 3   Optimal topology designs for L-shaped beam with stress and volume constraints

图4给出了体积和应力约束下L形梁优化问题目标函数和约束函数的迭代历史. 最终的柔顺度、体积分数及最大的von Mises应力值(无量纲)列于表1中. 从图4可以看出,本文方法迭代过程稳健,最终设计很好地满足了体积约束和应力约束. 从表1可以看出,基于体积约束下的最小柔顺性优化模型获得的最终结构的最大应力值(9.649 3)大于初始结构的最大应力值(7.770 6),优化后结构虽然减重了,但是L形梁拐角处的应力集中效应会使结构容易发生断裂等破坏行 为;而基于式(1)描述的优化模型获得的结构的最大应力值(6.004 6)比初始结构的最大应力值(7.770 6)要小, 这样在对结构进行减重的同时还增加了结构的强度.

图4   体积和应力约束下L形梁目标函数和约束函数的迭代历史

Fig. 4   Iterative histories of objective and constraint functions of L-shaped beam with volume and stress constraints

表1   算例1 L形梁优化结果

Table 1   Optimal results for L-shaped beam in Ex.1

CaseComplianceVolumeMax stress
Case 120 2550.49.649 3
Case 222 2190.46.004 6

新窗口打开

4.2 算例2

第2个算例考虑L形梁右边中点B处(图1)受集中载荷作用的优化问题. 在本算例中,设置材料体积分数和应力约束限值分别为$f_v=0.4, \sigma ^{*}= 6.0$. 载荷的大小和方向以及结构的离散网格与算例1相同.为了避免应力集中效应,与算例1类似,将集中载荷均匀分布到相邻的5个网格节点上.初始结构在载荷$F$作用下的最大的von Mises应力为7.772 5,过滤半径取$r_{min} = 3.5$.

图5展示了本文方法获得的优化结果,其中图5(a)是体积约束下的刚度最大化优化结果,图5(b)是体积约束下的应力最小化优化结果,图5(c)是体积和应力约束下的刚度最大化优化结果. 为了便于辨识,其优化结果分别标记为Case 1 、Case 2和Case 3. 此处Case 3的应力约束限值是根据Case 2的优化结果选取的. 最终的柔顺度、体积分数及最大的von Mises应力值列于表2中.

图5   L形梁的拓扑优化设计

Fig. 5   Optimal topology designs for L-shaped beam

表2   算例2 L形梁优化结果

Table 2   Optimal results for L-shaped beam in Ex.2

CaseComplianceVolumeMax stress
Case 120 0080.410.173 9
Case 238 6400.45.778 5
Case 321 0990.45.999 7

新窗口打开

图5表2可知,无论是否考虑应力约束或最小化,本文改进的BESO方法均可以获得清晰的拓扑构型. 考虑应力的拓扑构型与不 考虑应力的拓扑构型不同,体积约束下刚度最大化的L形梁(Case 1的结果)在拐角处存在明显的应力集中效应,优化后结构的最大应力(10.173 9)要比初始结构的最大应力(7.772 5)要大,而考虑 应力的优化设计(Case 2和Case 3的结果)能够有效降低应力集中效应. 工程中的应力集中效应可能导致结构破坏事故,因此,对于设计者而言,从构型概念设计时就限制最大应力是非常有必要的.

4.3 算例3

第3个算例考虑T形梁优化问题,几何区域和边界条件如图6所示,结构左边和右边完全固支,顶部左边点受集中载荷$F =10$作用,为了避免载荷作用点处的应力集中,将载荷平均分配到如图6所示的邻近的5个节点上.初始结构在载荷$F$作用下的最大的von Mises应力为4.563 3. 在本算例中,设置材料体积分数和应力约束限值分别为$f_v=0.4$, $\sigma ^{\ast }= 3.8$,过滤半径取$r_{\min}= 3.5$. 结构由 8 800 个四节点四边形单元来离散.

图6   T形梁设计域和边界条件

Fig. 6   Design domain and boundary condition for T-shaped beam

图7展示了本文方法获得的优化结果,其中图7(a)是体积约束下的刚度最大化优化结果,图7(b)是体积约束下的应力最小化优化结果,图7(c)是体积和应力约束下的刚度最大化优化结果. 为了便于辨识,其优化结果分别标记为Case1 , Case 2和Case 3. 此处Case 3的应力约束限值也是根据Case 2的优化结果选取的. 最终的柔顺度、体积分数及最大的von Mises应力值列于表3中.

图7   T形梁的拓扑优化设计

Fig. 7   Optimal topology designs for T-shaped beam

表3   T形梁优化结果

Table 3   Optimal results for T-shaped beam

CaseComplianceVolumeMax stress
Case 17 107.10.46.432 3
Case 210 3850.43.226 9
Case 38 496.20.43.801 5

新窗口打开

从优化结果可以看出基于刚度优化的设计(图7(a))与应力相关设计(图7(b)、图7(c))的拓扑构型明显不同,体积约束下基于刚度优化的T形梁在两个拐角处出现了明显的应力集中现象,其最大应力(6.432 3)比初始结构的最大应力(4.563 3)要大,而考虑应力的设计能够有效降低应力集中效应. 这说明本文改进的BESO方法可以有效地获得降低应力集中效应的拓扑构型设计.

5 结 论

本文建立改进的BESO方法求解了应力约束的拓扑优化问题,采用K-S函数来凝聚局部应力约束以减少与局部应力约束相关的计算代价,利用拉格朗日乘子法施加应力约束,由二分法确定拉格朗日乘子的值. 然后,详细推导了基于BESO方法的应力约束拓扑优化模型及其灵敏度列式,三个算例表明本文方法可有效地处理应力相关的拓扑优化问题. 比较发现,考虑应力约束和不考虑应力约束的拓扑构型不同;考虑应力约束的设计能够有效降低结构关键区域的应力集中效应.

The authors have declared that no competing interests exist.


参考文献

[1] Bendsoe MP, Kikuchi N.

Generating optimal topologies in structural design using a homogenization method

.Computer Methods in Applied Mechanics and Engineering, 1988, 71(1): 197-224

[本文引用: 1]     

[2] Bendsoe MP, Sigmund O.

Topology Optimization: Theory, Methods and Applications

. Berlin: Springer, 2003

[本文引用: 1]     

[3] Xie YM, Steven GP.

A simple evolutionary procedure for structural optimization

.Computers and Structures, 1993, 49(3): 885-896

[本文引用: 1]     

[4] 隋允康, 叶红玲, 彭细荣.

连续体结构拓扑优化应力约束凝聚化的ICM 方法

. 力学学报, 2007, 23(4): 554-563

[本文引用: 3]     

(Sui Yunkang, Ye Hongliang, Peng Xirong, et al.

The ICM method for continuum structural topology optimization with condensation of stress constraints

.Chinese Journal of Theoretical Applied Mechanics, 2007, 23(4): 554-563 (in Chinese))

[本文引用: 3]     

[5] Wang MY, Wang X, Guo DM.

A level set method for structural topology optimization

.Computer Methods in Applied Mechanics and Engineering, 2003, 192(1): 227-246

[本文引用: 1]     

[6] Guo X, Zhang WS, Zhong W.

Doing topology optimization explicitly and geometrically—a new moving morphable components based framework

.Journal of Applied Mechanics, 2014, 81(6): 081009

[本文引用: 1]     

[7] Sigmund O, Maute K.

Topology optimization approaches

.Structural and Multidisciplinary Optimization, 2013, 48(4): 1031-1055

[本文引用: 1]     

[8] 谢亿民, 黄晓东, 左志豪.

渐进结构优化法 (ESO) 和双向渐进结构优化法 (BESO) 的近期发展

. 力学进展, 2011, 41(4): 462-471

[本文引用: 2]     

(Xie Yiming, Huang Xiaodong, Zuo Zhihao, et al.

Recent developments of evolutionary structural optimization (ESO) and bidirectional evolutionary structural optimization (BESO) methods

.Advances in Mechanics, 2011, 41(4): 462-471 (in Chinese))

[本文引用: 2]     

[9] 张卫红, 郭文杰, 朱继宏.

部件级多组件结构系统的整体式拓扑布局优化

. 航空学报, 2015, 36(6): 2662-2669

[本文引用: 1]     

(Zhang Weihong, Guo Wenjie, Zhu Jihong.

Integrated layout and topology optimization design of multi-component systems with assembly units

.Acta Aeronautica et Astronautica Sinica, 2015, 36(6): 2662-2669 (in Chinese))

[本文引用: 1]     

[10] 龙凯, 王选, 韩丹.

基于多相材料的稳态热传导结构轻量化设计

. 力学学报, 2017, 49(1): 359-366

[本文引用: 1]     

(Long Kai, Wang Xuan, Han Dan.

Structural light design for steady heat conduction using multi-material

.Chinese Journal of Theoretical and Applied Mechanics, 2017, 49(1): 359-366 (in Chinese))

[本文引用: 1]     

[11] 王选, 胡平, 祝雪峰.

考虑结构自重的基于NURBS插值的3D拓扑描述函数法

. 力学学报, 2016, 48(4): 1437-1445

[本文引用: 1]     

(Wang Xuan, Hu Ping, Zhu Xuefeng, et al.

Topology description function approach using NURBS interpolation for 3Dstructures with self-weight loads

.Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(4): 1437-1445 (in Chinese))

[本文引用: 1]     

[12] 隋允康, 彭细荣.

求解一类可分离凸规划的对偶显式模型DP-EM方法

. 力学学报, 2017, 49(3):1135-1144

[本文引用: 1]     

(Sui Yunkang, Peng Xirong.

A dual explicit model based DPEM method for solving a class of separable convex programming

.Chinese Journal of Theoretical and Applied Mechanics, 2017, 49(3): 1135-1144 (in Chinese))

[本文引用: 1]     

[13] 陈文炯, 刘书田, 张永存.

基于拓扑优化的自发热体冷却用植入式导热路径设计方法

. 力学学报, 2016, 48(1): 406-412

[本文引用: 1]     

(Chen Wenjiong, Liu Shutian, Zhang Yongcun.

Optimization design of conductive pathways for cooling a heat generating body with high conductive inserts

.Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(1): 406-412 (in Chinese))

[本文引用: 1]     

[14] Guo X, Zhang WS, Zhang J, et al.

Explicit structural topology optimization based on moving morphable components (MMC) with curved skeletons

.Computer Methods in Applied Mechanics and Engineering, 2016, 310: 711-748

[本文引用: 1]     

[15] Cheng GD, Jiang Z.

Study on topology optimization with stress constraints

.Engineering Optimization, 1992, 20(1): 129-148

[本文引用: 2]     

[16] Cheng GD, Guo X.

ε-relaxed approach in structural topology optimization

.Structural Optimization, 1997, 13(4): 258-266

[本文引用: 2]     

[17] Duysinx P, Bendsoe MP.

Topology optimization of continuum structures with local stress constraints

.International Journal for Numerical Methods in Engineering, 1998, 43(6): 1453-1478

[本文引用: 3]     

[18] Paris J, Navarrina F, Colominas I, et al.

Topology optimization of continuum structures with local and global stress constraints

.Structural and Multidisciplinary Optimization, 2009, 39(4): 419-437

[本文引用: 4]     

[19] Bruggi M, Venini P.

A mixed FEM approach to stress-constrained topology optimization

.International Journal for Numerical Methods in Engineering, 2008, 73(10): 1693-1714

[本文引用: 1]     

[20] Verbart A, Langelaar M, Van Keulen F.

A unified aggregation and relaxation approach for stress-constrained topology optimization

.Structural and Multidisciplinary Optimization, 2017, 55(1): 663-679

[本文引用: 5]     

[21] Holmberg E, Torstenfelt B, Klarbring A.

Stress constrained topology optimization

.Structural and Multidisciplinary Optimization, 2013, 48(1): 33-47

[本文引用: 3]     

[22] Rong JH, Xiao TT, Yu LH, et al.

Continuum structural topological optimizations with stress constraints based on an active constraint technique

.International Journal for Numerical Methods in Engineering, 2016, 108(4): 326-360

[本文引用: 1]     

[23] Le C, Norato J, Bruns T, et al.

Stress-based topology optimization for continua

.Structural and Multidisciplinary Optimization, 2010, 41(4): 605-620

[本文引用: 6]     

[24] Yang RJ, Chen CJ.

Stress-based topology optimization

.Structural Optimization, 1996, 12(2-3): 98-105

[本文引用: 2]     

[25] Luo YJ, Kang Z.

Topology optimization of continuum structures with Drucker-Prager yield stress constraints

.Computers and Structures, 2012, 90: 65-75

[26] Luo YJ, Wang MY, Kang Z.

An enhanced aggregation method for topology optimization with local stress constraints

.Computer Methods in Applied Mechanics and Engineering, 2013, 254: 31-41

[本文引用: 2]     

[27] Paris J, Navarrina F, Colominas I, et al.

Block aggregation of stress constraints in topology optimization of structures

.Advances in Engineering Software, 2010, 41(2): 433-441

[本文引用: 2]     

[28] Cai SY, Zhang WH.

Stress constrained topology optimization with free-form design domains

.Computer Methods in Applied Mechanics and Engineering, 2015, 289: 267-290

[本文引用: 3]     

[29] Wang MY, Li L.

Shape equilibrium constraint: A strategy for stress-constrained structural topology optimization

.Structural and Multidisciplinary Optimization, 2013, 47(2): 335-352

[本文引用: 3]     

[30] Zhou MD, Sigmund O.

On fully stressed design and p-norm measures in structural optimization

.Structural and Multidisciplinary Optimization, 2017, 56(2): 731-736

[本文引用: 1]     

[31] Guo X, Zhang WS, Zhong W.

Stress-related topology optimization of continuum structures involving multi-phase materials

.Computer Methods in Applied Mechanics and Engineering, 2014, 268: 632-655

[本文引用: 2]     

[32] Xia Q, Shi T, Liu S, et al.

A level set solution to the stress-based structural shape and topology optimization

.Computers and Structures, 2012, 90: 55-64

[33] Guo X, Zhang WS, Wang MY, et al.

Stress-related topology optimization via level set approach

.Computer Methods in Applied Mechanics and Engineering, 2011, 200(47): 3439-3452

[本文引用: 1]     

[34] Cai SY, Zhang WH, Zhu JH, et al.

Stress constrained shape and topology optimization with fixed mesh: A B-spline finite cell method combined with level set function

.Computer Methods in Applied Mechanics and Engineering, 2014, 278: 361-387

[本文引用: 1]     

[35] 隋允康, 张学胜, 龙连春.

应力约束处理为应变能集成的连续体结构拓扑优化

. 计算力学学报, 2007, 24(3): 602-608

[本文引用: 1]     

(Sui Yunkang, Zhang Xuesheng, Long Lianchun.

The ICM method of structural topology optimization with stress constraints approached by the integration of strain energies

. Chinese Journal of Computational Mechanics, 2007, 24(3): 602-608 (in Chinese))

[本文引用: 1]     

[36] 隋允康, 边炳传.

屈曲与应力约束下连续体结构的拓扑优化

. 工程力学, 2008, 25(6): 6-12

(Sui Yunkang, Bian Binchuan.

Topology optimization of continuum structures under bucking and stress constraints

.Engineering Mechanics, 2008, 25(6): 6-12 (in Chinese))

[37] 隋允康,铁军.

结构拓扑优化ICM显式化与抛物型凝聚函数对于应力约束的集成化

. 工程力学, 2010, 27(增刊I): 224-237

[本文引用: 1]     

(Sui Yunkang, Tie Jun.

The ICM explicitation approach to structural topology optimization and the integrating approach to stress constraints based on the parabolic aggregation function

.Engineering Mechanics, 2010, 27(Sup. issue I): 224-237 (in Chinese))

[本文引用: 1]     

[38] 隋允康, 叶红玲. 连续体结构拓扑优化的ICM方法. 北京:科学出版社, 2013

[本文引用: 1]     

(Sui Yunkang, Ye Hongling.Continuum Topology Optimization Methods ICM. Beijing: Science Press, 2013 (in Chinese))

[本文引用: 1]     

[39] 荣见华, 葛森, 邓果.

基于位移和应力灵敏度的结构拓扑优化设计

. 力学学报, 2009, 41(4): 518-529

[本文引用: 1]     

(Rong Jianhua, Ge Sen, Deng Guo, et al.

Structural topology optimization based on displacement and stress sensitivity analysis

.Chinese Journal of Theoretical Applied Mechanics, 2009, 41(4): 518-529 (in Chinese))

[本文引用: 1]     

[40] Munk DJ, Vio GA, Steven GP.

Topology and shape optimization methods using evolutionary algorithms: A review

.Structural and Multidisciplinary Optimization, 2015, 52(2): 613-631

[本文引用: 1]     

[41] Huang X, Xie YM.

Evolutionary topology optimization of continuum structures with an additional displacement constraint

.Structural and Multidisciplinary Optimization, 2010, 40(1): 409-416

[本文引用: 3]     

[42] Huang X, Xie YM.

Convergent and mesh-independent solutions for the bi-directional evolutionary structural optimization method

.Finite Elements in Analysis and Design, 2007, 43(12): 1039-1049

[本文引用: 3]     

[43] Huang X, Li Y, Zhou SW, et al.

Topology optimization of compliant mechanisms with desired structural stiffness

.Engineering Structures, 2014, 79: 13-21

[本文引用: 2]     

/