式中:K为整体热传导矩阵;T为温度列阵;P为热流率列阵.根据SIMP插值理论,单元热传导矩阵[15]可写为
式中:ki为单元热传导矩阵;ki0为初始单元热传导矩阵;xi为单元i的相对密度;p为惩罚因子.建立多工况热传导结构的拓扑优化数学模型为
式中:C为结构散热弱度;S为工况数;j为工况号;N为设计域内单元数;wj为工况权系数;tij为单元i的温度列阵;V(x)为优化后结构体积;V0为初始结构体积;Vi为单元i的体积;f为体积比;xmin为密度下限,为防止数值奇异,通常取0.001.假设热流率与单元相对密度无关,则通过伴随法易推导得
式中:Ci为单元散热弱度.由式(4)可知,敏度可以表达为单元散热弱度的代数运算形式.常见的拓扑优化求解算法有优化准则(Optimality Criteria,OC)法和数学规划算法两类.优化准则法具有收敛速度快、计算规模和变量数目无关的优点.根据式(4)求解得到敏度后,标准OC法更新设计变量表达式[21]为
式中:k为第k轮优化迭代过程;m为移动极限；η为阻尼系数，m及η用于确保数值计算稳定性和收敛性,分别取值为0.2和0.5.
式中:
为体积敏度;λ为拉格朗日乘子,根据体积约束方程采用二分法求解.1.2 克服拓扑优化中的数值不稳定现象棋盘格现象及网格依赖性问题是连续体拓扑优化中普遍存在的数值不稳定现象.目前常见的解决办法包括高阶等参元法、周长约束法、全局或局部梯度法、过滤法等[22].敏度过滤法是目前应用最为广泛的方法之一[23].敏度过滤通常可以通过卷积和求解PDE等方式实现.通过卷积进行敏度过滤与求解PDE进行敏度过滤本质上是一致的,但采用有限元方式离散PDE可以使用与结构分析相同的网格信息,相对于搜索式的过滤方式,可以避免占用大量内存,提高过滤效率,适用于大规模优化求解问题[24].由式(4)可知,敏度取决于单元散热弱度,因此可以通过求解PDE对单元散热弱度过滤来实现敏度过滤.单元散热弱度过滤后形成的场函数
隐含在: 
并满足Neumann边界条件：
式中:C为单元散热弱度过滤前形成的场函数;Rmin为PDE参数,其作用类似于卷积敏度过滤中的过滤半径rmin,两者对应关系[24]为
1.3 周期性结构建模及结果后处理周期性结构的本质就是在各个子区域内具有相同的拓扑形式.对于一个二维设计域问题,如图 1所示,将设计域划分为mx×my个子区域,mx和my分别为沿x轴和y轴的子区域个数.周期性约束的数学模型表达式为
 |
| 图 1 周期性结构示意图Fig. 1 Illustration of periodic structure |
| 图选项 |
式中:上标(1),(2),…,(mx×my)为子区域编号;q=1,2,…,n,其中n=
为单个子区域内的有限单元数.为实现式(10),不同子区域内,相同位置的单元敏度需相等;由式(4)可知,只需进一步满足: 
与文献[25]类似,在第 k 轮优化迭代过程中,可以通过重新分配过滤之后的单元散热弱度来实现周期性约束.为了保证当轮优化过程中的目标函数 C 不因周期性约束发生变化,要求分配前后的散热弱度总和不变.考虑周期性约束的单元散热弱度表达式为
在重新分配单元散热弱度之后,根据式(4)、式(5)来分别计算单元敏度和更新单元相对密度.由于此时每个子区域内同一位置的单元散热弱度相同,即敏度相同,就可以实现周期性结构.在每轮优化迭代中,采用上述OC算法进行优化求解,直至满足收敛准则:
式中:x为由设计变量xi组成的向量.由于过滤的平均作用,优化边界存在着大量的中间密度单元,从而导致优化边界模糊,为了得到清晰的拓扑优化结果,优化过程分两个阶段：第一阶段采用过滤以消除数值不稳定性现象,ε取值为0.5%;在此之后继续采用无过滤方式进行优化直至收敛精度达到0.1%.2 数值算例及讨论本节算例均采用平面模型,算例均在MATLAB软件中编程实现[26].用双线性四边形单元离散结构,忽略结构尺寸及材料单位,单元大小为0.25×0.25,在以下算例中,如不特殊说明,材料为各向同性导热材料,p=3,rmin=1.5.随着子区域个数增多,每个子区域内的单元数相应减少,为了严格满足体积比约束,不可避免地,个别单元要取中间值,采用文献[27]中提及的后处理方法显示宏观、微观拓扑构型.算例1 如图 2所示,平面结构尺寸为180×180,4个角点温度恒为0,设平面结构的左上角为坐标原点,加热点坐标位置为(82,80),热流率为10,导热系数为10.令mn=mx=my,限于篇幅,仅列出mn为1～6、9、10、12、15、36、40、48和60的宏观、微观最优拓扑构型及优化结果.最优拓扑构型如图 3所示.结构散热弱度随mn的变化趋势如图 4 所示,计算得到的规律同样适用于其他算例.
 |
| 图 2 数值算例1中受热平面结构示意图Fig. 2 Illustration of heated plane structure in numerical example 1 |
| 图选项 |
由图 3可知,当子区域个数为1时,即为非周期性结构拓扑优化构型,材料从加热点至冷却点分布,体现了最短传热路径这一特点.由图 4可知,相比较其他子区域个数的最优结构,其散热弱度值最小,说明了非周期性结构拓扑优化结果最优.优化目标总体上随着子区域个数增加而变差,这主要是随着子区域个数增加,单个子区域内的优化空间相应减少,导致宏观周期性结构热传导能力变差.但就本文算例而言,与文献[26]结论有所不同,优化结果呈非严格单调递增性.
 |
| 图 3 数值算例1中不同子区域个数的拓扑优化构型和微结构构型Fig. 3 Optimal topological configurations and microstructure configurations with different numbers of unit cells in numerical example 1 |
| 图选项 |
 |
| 图 4 数值算例1中散热弱度随子区域个数的变化曲线Fig. 4 Curve of optimal results changing with different numbers of unit cells in numerical example 1 |
| 图选项 |
当子区域个数不同时,最优拓扑呈现不同的宏观和微观结构.由于宏观结构大小一定,不同的子区域数实际上反映了材料微结构的尺寸,即最优拓扑构型表征了尺寸效应对材料微结构构型的影响.但随着子区域数目不断增加,尽管微观的单胞构型不同,但由其排列而成的宏观结构却逐渐呈现一致性.如图 5所示的体胞,结构尺寸为1×1,以尺寸a和b为变量,采用均匀化方法得到如图 6所示的宏观、微观结构.
 |
| 图 5 均匀化方法表征体胞Fig. 5 Unit cell configuration obtained by homogenization method |
| 图选项 |
 |
| 图 6 均匀化方法下的宏观、微观结构Fig. 6 Optimal topological configuration and microstructure configuration obtained by homogenization method |
| 图选项 |
以mn=60为例,材料在x和y方向分布比例为1.89,均匀化结果为1.61,材料分配比例接近.采用均匀化方法得到的散热弱度为62.989 4,由图 4可知,周期性结构散热弱度趋于收敛于均匀化方法计算得到的极限值.部分情况下,优化结构中出现的细小分支结构会对材料微结构加工制造造成影响,因此该方法在最优结构的可加工性方面需进一步改进.算例2 基本参数如算例1,加热点位置为(30,106),为简单起见,仅列出mn为1～6、9、10、12和15的宏观优化构型,如图 7所示.
 |
| 图 7 数值算例2中不同子区域数个数的拓扑优化构型Fig. 7 Optimal topological configurations with different numbers of unit cells in numerical example 2 |
| 图选项 |
由图 3、图 7可知,受热方式不同,材料微结构构型有所不同.因此在材料微结构设计时有必要考虑宏观结构受热方式.对于mn=9的情况,从加热点到右下角点没有材料分布,此外,表 1中的数据显示,比较前后情况,其柔顺度值增大,说明了宏观周期性约束有可能会导致拓扑优化结果失去最短传热路径.表 1 数值算例2中不同子区域个数的拓扑优化结果 Table 1 Optimal topological results with different numbers of unit cells in numerical example 2
| 子区域个数 | 迭代次数 | 最优散热弱度 |
| 6×6 | 68 | 41.032 8 |
| 9×9 | 60 | 47.849 2 |
| 10×10 | 55 | 45.106 5 |
表选项
算例3 基本参数同算例1,考虑两个工况且每个工况的加权系数为0.5,加热位置分别在(82,80)和(30,106).仅列出mn为1～6、9、10、12和15的宏观优化构型,如图 8所示.
 |
| 图 8 数值算例3中不同子区域个数的拓扑优化构型Fig. 8 Optimal topological configurations with different numbers of unit cells in numerical example 3 |
| 图选项 |
由图 3、图 7和图 8比较发现,单、多工况下材料微结构有所不同,表明多目标建模会影响材料微结构构型进而会影响周期性宏观结构.3 结 论基于变密度法建立的周期性传热结构可以实现满足宏观导热性能的周期性材料设计,即“结构性”的材料设计.1) 在设计域内附加周期性约束,通过重新分配单元散热弱度以及基于PDE消除棋盘格和网格依赖性的方法可以实现宏观周期性传热结构设计.2) 数值实验结果表明,当宏观结构承载方式不同时,微观材料胞元呈现不同的拓扑优化构型.子区域个数不同时,优化得到的不同材料微结构实际上反映了材料的尺寸效应对材料设计的影响,并最终导致周期性结构有所不同.就本文算例而言,当子区域个数不断增加时,优化结果呈非严格单调递增性,且逐渐趋向收敛于均匀化方法对应的极限值.3) 在优化过程中出现的细小分支结构将会对材料微结构生产加工造成影响,因此结构在可加工性方面需进一步改进.
参考文献
| [1] | Sigmund O.Materials with prescribed constitutive parameters-an inverse homogenization problem[J].International Journal of Solids and Structures, 1994, 31(17):2313-2329. |
| Click to display the text | |
| [2] | Sigmund O, Torquato S.Design of materials with extreme thermal expansion using a three-phase topology optimization method[J].Journal of the Mechanics and Physics of Solids, 1997, 45(6):1037-1067. |
| Click to display the text | |
| [3] | Zhang W H, Wang F W, Dai G M, et al.http://html.rhhz.net/BJHKHTDXXBZRB/Topology%20optimal%20design%20of%20material%20microstructures%20using%20strain%20energy-based%20method[J].Chinese Journal of Aeronautics, 2007, 20(4):321-328. |
| Click to display the text | |
| [4] | 张卫红, 孙士平.多孔材料/结构尺度关联的一体化拓扑优化技术[J].力学学报, 2006, 38(4):522-529. Zhang W H, Sun S P.Integrated design of porous materials and structures with scale-coupled effect[J].Chinese Journal of Theoretical and Applied Mechanics, 2006, 38(4):522-529(in Chinese). |
| Cited By in Cnki (24) | |
| [5] | Hunag X, Xie Y M.Optimal design of periodic structures using evolutionary topology optimization[J].Structural Multidisciplinary Optimization, 2008, 36(6):597-606. |
| Click to display the text | |
| [6] | Zuo Z Z, Xie Y M, Huang X D.Optimal topological design of periodic structures for natural frequencies[J].Journal of Structural Engineering, 2011, 137(10):1229-1240. |
| Click to display the text | |
| [7] | Xie Y M, Zuo Z H, Huang X D, et al.Convergence of topological patterns of optimal periodic structures under multiple scales[J].Structural Optimization, 2012, 46(1):41-50. |
| Click to display the text | |
| [8] | Yan J, Cheng G D, Liu L.A uniform optimum material based model for concurrent optimization of thermoelastic structures and materials[J].International Journal for Simulation and Multidisciplinary Design Optimization, 2008, 2(4):259-266. |
| Click to display the text | |
| [9] | Liu L, Yan J, Cheng G D.Optimum structure with homogeneous optimum truss-like material[J].Computers and Structures, 2008, 86(13-14):1417-1425. |
| Click to display the text | |
| [10] | Niu B, Yan J, Cheng G D.Optimum structure with homogeneous optimum cellular material for maximum fundamental frequency[J].Structural and Multidisciplinary Optimization, 2009, 39(2):115-132. |
| Click to display the text | |
| [11] | 刘远东, 尹益辉, 郭中泽.静动态力学条件下多孔金属的材料/结构多级优化设计研究[J].中国科学:技术科学, 2012, 42(10):1172-1178. Liu Y D, Yin Y H, Guo Z Z.Static and dynamic design based on hierarchical optimization for materials and structure of porous metals[J].Scientia Sinica:Technologica, 2012, 42(10):1172-1178(in Chinese). |
| Cited By in Cnki (1) | |
| [12] | Li Q, Steven G P, Querin O M, et al.Shape and topology design for heat conduction by evolutionary structural optimization[J].International Journal of Heat and Mass Transfer, 1999, 42(17):3361-3371. |
| Click to display the text | |
| [13] | Li Q, Steven G P, Xie Y M, et al.Evolutionary topology optimization for temperature reduction of heat conducting fields[J].International Journal of Heat and Mass Transfer, 2004, 47(23):5071-5083. |
| Click to display the text | |
| [14] | Bendsoe M P, Sigmund O.Topology optimization-theory, methods and applications[M].Berlin:Springer, 2003:270-272. |
| [15] | 左孔天, 陈立平, 张云清, 等.用拓扑优化方法进行热传导散热体的结构优化设计[J].机械工程学报, 2005, 41(4):13-21. Zuo K T, Chen L P, Zhang Y Q, et al.Structual optimal design of heat conductive body with topology optimization method[J].Chinese Journal of Mechanical Engineering, 2005, 41(4):13-21(in Chinese). |
| Cited By in Cnki (46) | |
| [16] | 龙凯, 左正兴.稳态热传导下的连续体结构拓扑优化[J].中国机械工程, 2007, 18(24):2939-2943. Long K, Zuo Z X.Topological optimization of continuum structure for heat conduction[J].China Mechanical Engineering, 2007, 18(24):2939-2943(in Chinese). |
| Cited By in Cnki (11) | |
| [17] | 龙凯, 左正兴.多工况稳态热传导下的连续体拓扑优化[J].中国机械工程, 2007, 18(16):1925-1929. Long K, Zuo Z X.Topological optimization of continuum structure for heat conduction under multiple loading[J].China Mechanical Engineering, 2007, 18(16):1925-1929(in Chinese). |
| Cited By in Cnki (5) | |
| [18] | 刘书田, 贺丹.渐进密度AESO方法及其在热传导结构拓扑优化中的应用[J].计算力学学报, 2009, 26(2):151-156. Liu S T, He D.Progressive AESO algorithm and application in topology optimization of heat conduction structures[J].Chinese Journal of Computational Mechanics, 2009, 26(2):151-156(in Chinese). |
| Cited By in Cnki (11) | |
| [19] | Gao T, Zhang W H, Zhu J H, et al.Topology optimization of heat conduction problem involving design-dependent heat load effect[J].Finite Elements in Analysis and Design, 2008, 44(14):805-813. |
| Click to display the text | |
| [20] | 张晖, 刘书田, 张雄.拓扑相关热载荷作用下稳态热传导结构拓扑优化[J].中国机械工程, 2009, 20(11):1339-1343. Zhang H, Liu S T, Zhang X.Topology optimization of steady-state heat conduction problems with design-dependent heat loads[J].China Mechanical Engineering, 2009, 20(11):1339-1343(in Chinese). |
| Cited By in Cnki (8) | |
| [21] | Bendsoe M P.Optimization of structural topology, shape and material[M].Berlin:Spriner, 1995:83. |
| [22] | Sigmund O, Petersson J.Numerical instabilities in topology optimization:A survey on procedures dealing with checkerboards, mesh-dependencies and local minima[J].Structural and Multidisciplinary Optimization, 1998, 16(1):68-75. |
| Click to display the text | |
| [23] | Sigmund O, Maute K.Sensitivity filtering from a continuum mechanics perspective[J].Structural and Multidisciplinary Optimization, 2012, 46(4):471-475. |
| Click to display the text | |
| [24] | Lazarov B S, Sigmund O.Filters in topology optimization based on Helmholtz-type differential equations[J].International Journal for Numerical Methods in Engineering, 2011, 86(6):765-781. |
| Click to display the text | |
| [25] | 焦洪宇, 周奇才, 李文军, 等.基于变密度法的周期性拓扑优化[J].机械工程学报, 2013, 49(13):66-72. Jiao H Y, Zhou Q C, Li W J, et al.Periodic topology optimization using variable density method[J].Journal of Mechanical Engineering, 2013, 49(13):66-72(in Chinese). |
| Cited By in Cnki (9) | |
| [26] | Andreassen E, Clausen A, Schevenels M, et al.Efficient topology optimization in MATLAB using 88 lines of code[J].Structural and Multidisciplinary Optimization, 2011, 43(1):1-16. |
| Click to display the text | |
| [27] | Sigmund O, Maute K.Topology optimization approaches[J].Structural and Multidisciplinary Optimization, 2013, 48(6):1030-1055. |
