连续体结构拓扑优化的一种新变密度法_毛虎平

+
q （ k） + r0 x （j k） － L （j k）
（ k） 0j
)
值 / （ m / N） 值 / （ m / N） 比 / % 23． 645 4 38． 490 3 7． 390 1 7． 582 1 7 ． 219 6 7 ． 263 68． 75 80． 3 69 ． 58 81 ． 16
图1
式（ 1 ） 惩罚函数图 80 ˑ 20
300 200 100 80
MMA 算法在优化模型中引入人工变量 yi ， z， 以改善 每一个子优化问题的性态。 优化的数学模型表示如下：
min f槇0 （ x） + a0 z +
∑( c y
m i i =1
i
+
1 d y2 2 i i
)
（ 2）
120 ˑ 30
（ k） ） （ k） q0 q （i k — — 原优化模型中， 约束函数在设计点 x j 处关 j ， j —
于设计变量的一阶 Taylor 展开式；
（ k） r0 ， r （i k） — — — 目标函数和约束的余项 。
槇 （ k） （ k） r （i k） 很 在当前迭代点处 f0 （ x j ） = f0 （ x j ） ， 因此， 容易得到。 由设计 点 开 始，基 于 对 偶 法 解 子 优 化 问 题 （ 式 （ 2） ） ， 通过共轭梯度类算法解一系列线性方程来求得 子优化问题的， 反复迭代直至收敛。
［2 － 4 ］
s． t．
∑V x
j j =1
j
－ V ≤ 0 （ 0 ≤ xj ≤ 1）
（ 1）
— — 结构柔度； 式中： C— x— — — 设计变量， 即单元密度； U— — — 结构位移向量； Vj — — — 单元 j 的体积； V— — — 结构总体积；
n
K— — — 结构刚度矩阵， K =
3
数值算例与结果分析
为验证文中方法的有效性， 对两个问题进行了计
40
机
械
设
计
第 30 卷第 5 期
而目标函数的初始值相差很大， 单 得最优值非常接近， 元数 80 ˑ 20 时， 文中方法迭代次数比 SIMP 法都小， 单 文中方法迭代次数比 SIMP 法都大， 元数 120 ˑ 30 时， 而且增加的幅度也大。
转化成一系列显式的更为简单的严格凸的近似子优化 问题， 在每一步迭代中， 通过求解一个近似的凸子问
*
收稿日期： 2012 － 04 － 10 ； 修订日期： 2012 － 11 － 15 基金项目： 国家自然科学基金资助项目（ 51275489 ） ； 中北大学校基金资助项目（ 2012 ） ； 博士科研启动基金资助项目（ 2011 ） 作者简介： 毛虎平（ 1974 —） ， 男， 山西临汾人， 讲师， 博士， 研究方向： 机械优化设计。
（ a） p = 100
中方法的 CPU 耗时和迭代次数急剧增加。 下一步将深 入研究惩罚系数 p 对拓扑优化敏感度和效率的影响 。
参考文献
（ b） p = 300 图3 算例 1 的优化结果 ［ 1］ Bendsoe M， Kikuchi N． Generation optimal topologies in J］ ． Int J structural design using a homogenization method ［ Computer Methods in Applied Mechanics and Engineering， 1988 ， 71 ： 197 － 224． ［ 2］ ［ 3］ 罗震， 陈立平， 黄玉盈， 等． 连续体结构的拓扑优化设计 ［ J］ ． 力学进展， 2004 ， 34 （ 4 ） ： 463 － 476． 文中方法和 SIMP 法比较（ 惩罚系数 p = 200 ） 单元 数 柔度初 柔度终 下降 CPU 耗 时 /s 935 1 169 159 109 迭代 次数 88 64 102 127 ［ 6］ ［ 5］ ［ 4］ 罗震， 陈立平， 黄玉盈， 等． 基于 RAMP 密度 － 刚度插值 J］ ． 计算力学学报， 2005 ， 22 （ 5 ） ： 格式的结构拓扑优化［ 585 － 591． Bendsoe M P， Sigmund O． ， Topology optimization： theory［ M］ ． New York： Methods and Applications， 2003． 程耿东， 张东旭． 受应力约束的平面弹性体的拓扑优化 ［ J］ ． 大连理工大学学报， 1995 ， 35 （ 1 ） ： 1 － 9． Xie Y M， Stecen G P． Evolutionary structural optimization J］ ． Int J Computers and Structures， for dynamic problems ［ 1996 ， 58 （ 6 ） ： 1067 － 1073．
（ k） U （j k） — — — 移动渐近线方程； 式中： L j ，
30
38 ． 541 8
p
（ k） 0j
， p
（ k） ij
— — — 原优化模型中， 目标函数在设计点 x 于设计变量的一阶 Taylor 展开式；
处关
分析表 1 得出， 目标函数最优值随着 p 的变小而 p 从 500 到 50 ， 变小， 但是变化不大。 单元数 80 ˑ 20 ， 最 优值从 7 ． 669 4 到 7 ． 174 2 ， 其变化为 6 ． 46 % ； 单元数 120 ˑ 30 ， p 从 500 到 50 ， 最优值从 7 ． 335 8 到 7 ． 001 9 ， 其变化为 4 ． 55 % 。 而迭代次数和 CPU 耗时是从 p = 100 p从 和 200 向两边扩散。 分析迭代次数， 单元数 80 ˑ 20 ， 100 到 50 ， p从 迭代次数从 69 到 113 ， 其增加 63 ． 77 % ， 100 到 500 ， 迭代次数从 69 到 161 ， 其增加 133. 33 % ； 单 p 从 100 到 50 ， 元数 120 ˑ 30 ， 迭代次数从 69 到 113 ， 其 p 从 100 到 500 ， 增加 63 ． 77 % ， 迭代次数从 69 到 161 ， 其 增加 133 ． 33 % ； 从图 2 可以看出， 简支梁的拓扑变化甚 微， 变化最大的是单元数 120 ˑ 30 中 p = 300 时， 如图 3 所示。
*
摘要： 从变密度法 SIMP 和 RAMP 存在的问题出发， 提出一种指数形式的惩罚函数， 研究了惩罚系数对拓扑优化的影 响， 解决了 SIMP 和 RAMP 刚度矩阵奇异、 惩罚系数太大得不到合理拓扑 的问题， 得出 惩罚 系数 的 取值 范围， 两个 算 例分 析验证了文中提出方法的正确性和优势， 并和 SIMP 法进行了比较。 关键词： 拓扑优化; 变密度法; 指数函数 中图分类号： O342 文献标识码： A 文章编号： 1001 － 2354 （ 2013 ） 05 － 0038 － 04 min C（ x） = U T KU
1
刚度拓扑优化模型与新惩罚函数
以柔度最小化（ 刚度最大化） 作为优化的目标函数， 以结构整体的体积约束作为优化的约束条件， 在给定的 载荷和位移边界条件下， 基于新密度函数模型建立了线 弹性结构拓扑优化设计的在静力状态下的数学模型 如下：
2
基于移动渐近线法（ MMA） 优化策略
MMA 法通过引入移动渐近线， 将隐式的优化问题
∑p
j =1
Leabharlann Baidu
x j －1
E0 K j ；
Kj — — — 第 j 个单元刚度矩阵； E0 — — — 固体材料的弹性模量； p— — — 惩罚系数； n— — — 单元数目。
新密度函数插值模型：
E p （ x j ） = p x j － 1 E0
p — — 插值以后的弹性模量。 式中： E —
模型的柔度函数和结构柔度的敏度为 ：
DOI:10.13841/j.cnki.jxsj.2013.05.002
第 30 卷第 5 期 2013 年5 月
机
械
设
计
JOURNAL OF MACHINE DESIGN
Vol． 30 No． 5 May 2013
连续体结构拓扑优化的一种新变密度法
1 1 2 3 毛虎平 ， 苏铁熊 ， 李建军 ， 王应红 （ 1． 中北大学 机电工程学院，山西 太原 030051 ； 2． 中北大学 理学院， 山西 太原 030051 ； 3． 上海汽车集团股份有限公司乘用车公司， 上海 201804 ）
（ a） 单元数 80 ˑ 20
图4
单元数的影响
4
结语
从 SIMP 法和 RAMP 法存在的缺点出发， 提出一种
指数函数形式的惩罚函数， 消除了刚度矩阵奇异的问
（ b） 单元数 120 ˑ 30 图2 惩罚系数 p 的影响
单元密度为 0 ， 刚度矩 题（ 对于 SIMP 法和 RAMP 法中， 阵奇异， 因此单元密度 ≥ 0． 001 ） ， 去掉了单元密度最 小值的限制。 并验证了文中方法的正确性。 并研究了指 数形式惩罚函数的参数对拓扑优化的影响 ， 初步得出 参数 p 取 100 200 较好。 同时， 随着单元数的增加， 文
表1 单元 数 惩罚 文中方法惩罚系数的影响 柔度终 下降 CPU 耗 迭代 时 /s 487 341 282 180 460 320 6 908 5 908 5 390 2 275 2 552 次数 161 118 108 69 165 113 452 380 327 139 155
柔度初
系数 p 值/ （ m/ N） 值/ （ m/ N） 比 / % 500 58 ． 412 6 7 ． 669 4 46 ． 226 3 7 ． 577 2 38 ． 490 3 7 ． 582 1 28 ． 291 1 7 ． 473 1 25 ． 657 7 ． 293 2 20 ． 929 7 ． 174 2 58 ． 450 3 7 ． 335 8 46 ． 271 8 7 ． 340 1 38 ． 541 8 7 ． 263 28 ． 351 7 7 ． 174 9 20 ． 997 7 ． 001 9 86 ． 87 83 ． 61 80 ． 3 73 ． 58 71 ． 57 65 ． 72 87 ． 45 84 ． 14 81 ． 16 74 ． 69 66 ． 65
2013 年 5 月
毛虎平， 等： 连续体结构拓扑优化的一种新变密度法
39
［10 －11 ］ 。 题， 然后由对偶法或原始对偶 － 内点法求解
算分析。 例 1 ： 24 ˑ 8 平面简支梁， 两端中间固定， 弹性模量 E = 2 ． 0 ˑ 10 11 Pa， 泊松比 ν = 0 ． 3 ， 体积比 f = 0 ． 5 ， 梁 体中间下边受一个集中竖直向下载荷 F = 20 kN。 计算 结果如表 1 、 表 2 所示。
n
1988 年 Bendson 和 Kikuchi 提 出 了 均 匀 化 方 ［1 ］ 法 ， 此后连续体结构的拓扑优化设计新的方法不断 、 同时应用也不断被扩展， 如变密度法 变厚 出现， ［5 ］ ［6 － 7 ］ ［8 ］ 度法 、 进化法 （ ESO ） 及 ICM 方法等。 变密度 法的关 键 问 题 是 如 何 构 造 函 数。 代 表 性 的 有 SIMP （ Solid Isotropic Microstructure with Penalization ） 和 RAMP （ Rational Approximation of Material Properties ） 。 前者构造的函数为一指数函数， 惩罚因子 p 取值越大， 效果越好， 但太大又容易引起棋盘格问题； RAMP 方法 构造的函数为有理函数。这两种方法都是使小于 1 的 大部分中间密度单元趋向于 0 ， 并且相对密度小于 0． 8 的单元都迅速趋向于 0 ， 非常接近 1 的少部分单元趋 ［9 ］ 向于 1 ， 这对于拓扑优化过程是不利的。 而且这两 一般认为单元密 种变密度 法 都 存 在 奇 异 矩 阵 问 题， 度≤0． 001 时材料为空。 文中的出发点是构造惩罚函数天然不存在奇异矩 而且使单元密度小于 0． 8 的相对密度缓慢趋 阵问题， 向于 0 ， 增加接近 1 的单元数量趋向于 1 。 从这一思想 出发进行研究， 克服 SIMP 和 RAMP 方法的不足。
50 500 300 200 100 50 表2 材料 模型 单元 数
s． t．
{
f槇i （ x） － a i z － y i ≤ f槇i x min j ≤ x j ≤ x max j
（ z ≥ 0， yi ≥ 0）
— — 约束的个数； 式中： i— ai ， ci ， di — — — 第 i 个约束时给定的大于等于零的常数， ci + di ＞ 0； a0 — — — 给定的大于等于零的常数； f槇0 （ x） ， f槇i （ x） — — — 连续可微的函数。
文中方法和 SIMP 法比较（ p = 200 ） 柔度初 柔度终 下降 CPU 耗 迭代 时 /s 407 282 2 240 5 390 次数 133 108 135 327
原优化模型中目标函数的 MMA 近似展开式为：
f槇0 （ x） =
∑( U
n j =1
p
（ k） j
（ k） 0j
－ x （j k）
n
C（ x） = C' （ x） =
∑p
j =1
x j －1
E0 U T j Kj Uj
x j －1
∑lg pp
E0 U T j Kj Uj
30 ， 50 ， 100 ， 200 ， 500 ， 5 000 时的惩 图 1 为 p 取 10 ， 由图中可以看出 p 取 100 和 200 较为理想。 罚曲线。
原优化模型中约束方程的 MMA 近似展开式为：
f槇i （ x） =
SIMP 80 ˑ 本文 20
∑(
n j =1
） p （ij k） q （i k j + r （i k） + U （j k） － x （j k） x （j k） － L （j k）
)
SIMP 120 ˑ 23 ． 735 3 本文
（ k） j