为了解决体积力情况下的拓扑优化设计问题，大量的****对此开展了研究。1984年，陈树勋和叶尚辉[8]为了解决天线设计问题，采用了一种考虑自重的结构优化方法，称为导重准则法。1989年，Bendsoe[9]提出了基于固体各向同性材料惩罚(Solid Isotropic Material with Penalization，SIMP)模型的变密度法，用于解决连续体结构的拓扑优化问题。Pedersen[10]通过对采用SIMP模型的变密度法的研究，发现在考虑自重载荷时，会出现结构无法支撑自身质量的情况，该情况会导致优化结果中材料分布模糊，称为低密度区域的“附属”效应。Stolpe和Svanberg[11]提出了材料属性有理近似函数插值(Rational Approximation of Material Properties，RAMP)模型，该模型增加了拓扑优化过程中获得松弛问题0-1解的可能性，在一些情况下能够获得更加明确的材料分布。Bruyneel和Duysinx[12]通过对SIMP模型进行修正，选择基于梯度的运动渐进线逼近方法，提出了一种具有密度依赖性的拓扑优化方法，证明了体积力对结构拓扑优化的巨大影响。Lopes等[13]在SIMP模型的基础上，创立了一种刚度和密度特性参数化的材料模型，提出了一种适用于惯性力的三维结构拓扑优化设计公式。Liu和Li等[14-15]将陈树勋等提出的导重准则法引入到了连续体的拓扑优化问题中，并且对其做了改进，使其兼具优化准则法和数学规划法的特点。随后，陈祥和刘辛军[16]比较了SIMP模型和RAMP模型各自的特点，基于RAMP模型结构导重法对拓扑优化问题进行了求解。张晖等[17]针对自重载荷作用下的连续体拓扑优化问题提出了采用RAMP模型和平均敏度过滤技术相结合的求解策略。许华旸等[18-19]继续基于导重法对考虑自重载荷作用的拓扑优化做了相关研究。Jain和Saxena[20]比较自重和不同比值关系的外载荷作用下的优化结果，说明了自重对结构拓扑优化的影响及在优化设计中考虑结构自重的必要性。
自重载荷下结构拓扑优化会出现中间密度等问题，为此大量****开展了研究。Bruyneel和Duysinx[12]针对低密度区域的“附属”效应问题，采用了在合适位置添加非结构质量的方式解决问题。许华旸等[19]同样也通过添加质量点的方式获得明确的材料分布。Holmberg等[21]比较了附加非结构质量情况下结构受各种最糟糕加速度时得到的不同优化结果，强调了结构自重的影响。虽然非结构质量在考虑体积力等问题的拓扑优化研究中已有大量的研究，但是鲜有非结构质量对拓扑优化的影响研究尤其是悬臂梁结构拓扑优化的研究。鉴于此，本文对结合导重准则法和RAMP模型对自重载荷作用下的结构拓扑优化进行了研究，在该研究的基础上针对悬臂梁结构在自重载荷作用下末端容易出现中间密度单元模糊区域的情况，采用在梁末端添加非结构质量的求解策略，讨论了非结构质量对优化结果的影响规律。
1 理论和方法 1.1 变密度法 变密度法是结构拓扑优化应用最为广泛的方法，其核心思想是：对材料的伪密度进行惩罚使其归于0和1，从而完成材料的去除和保留。对于变密度法的惩罚模型，目前较为常用的有SIMP模型和RAMP模型。基于SIMP模型的变密度法又叫SIMP法，由Bendsoe[9]首次引入变密度法中，其表达式为
 | (1) |
式中：Ei为单元的弹性模量；E0和Emin分别为满材料单元和孔洞材料单元的弹性模量；xi为单元的伪密度；p为SIMP模型的惩罚因子，一般取值为3。
RAMP模型由Stolpe和Svanberg[11]提出，该模型可以让施加有设计相关载荷时的优化设计获得更加稳定的优化结果，其表达式为
 | (2) |
式中：q为RAMP模型的惩罚因子，一般取值为5~120。在结构拓扑优化中，自重载荷为设计相关载荷，故而采用RAMP模型。
变密度法拓扑优化常用的求解方法有数学规划法、优化准则法及导重准则法，导重准则法结合了数学规划法和优化准则法的优势，具有更高的求解效率。
1.2 体积约束下最大刚度的自重载荷结构拓扑优化
1.2.1 结构等效节点重力的计算 重力载荷作为一种设计相关载荷，其大小会随着设计变量的更新而变化。拓扑优化依赖于有限元计算，因此需要将单元的重力等效到单元节点上。通常，将每个单元重力按照虚功相等的原则分配到单元节点之上。
有限元分析采用四节点矩形单元，根据四节点矩形单元的形函数及单元体积密度和质量的关系，可以得到四节点矩形单元的第i个节点在重力加速度方向的等效节点力为[18]
 | (3) |
式中：a和b分别为单元1/2长和1/2宽；Ni表示四节点矩形单元第i个节点的形函数，1≤i≤4；ρ为材料的物理密度；g为重力加速度；A为单元面积；t为单元的厚度；Ve为一个单元的体积。
在四节点单元有限元分析中，节点力为二维矢量，与重力垂直方向的力为0，可以写出节点的体积力矢量为
 | (4) |
1.2.2 迭代公式和灵敏度 对结构进行拓扑优化设计时，以结构质量/体积为约束、最大刚度为目标是一种非常常用的优化方案。在结构的拓扑优化中，为了方便计算，通常以柔度最小来代替刚度最大的优化目标。
以体积为约束、柔度最小为目标的拓扑优化数学模型，其表达式为
 | (5) |
式中：X为N维的设计变量，即伪密度；C(X)为结构的柔度；ui为单元i的位移矢量；ki为单元i的刚度矩阵；g(X)为结构的质量/体积约束；M为结构在迭代优化过程中的实际质量/体积；M0为约束的质量/体积；xiL和xiU分别为xi的下限和上限。
根据Kuhn-Tucker条件对式(5)进行求解，得到导重准则法求解的迭代准则[14-15]为
 | (6) |
式中：λ为拉格朗日乘子；
迭代公式可写为
 | (7) |
式中：k为迭代的次数。
对目标函数关于设计变量求导数，得到其灵敏度[18]为
 | (8) |
式中：P0=[0, -1, 0, -1, 0, -1, 0, -1]为单元四节点载荷的位移方向矢量；Fi为第i个单元四节点等效载荷矢量；案例中采用RAMP模型，
分析式(8)可知，当区域应变能过小时，由于数值化整误差的影响，该区域的灵敏度会非常接近甚至相同，从而导致不收敛的情况。
综上所述，可以得到该案例的迭代公式及拉格朗日乘子的计算式[18]为
 | (9) |
 | (10) |
式中：α为步长因子，合适的取值可以保证迭代的收敛。关于α取值，陈树勋[22-23]进行了详细的数学推导和证明。在取值范围内，α值越小，越能精确地搜索最优解，但相应的计算量也越大。考虑到计算量问题，取α=0.5能够适应大多数情况。
利用式(9)和式(10)计算拉格朗日乘子的值会产生很大的计算量，并且不能保证式(9)和式(10)同时取得最优值，解决这个问题的方法是：利用二分法来求拉格朗日乘子，在求解迭代过程中，需要设定一个判断迭代是否停止的收敛条件，如下：
 | (11) |
取ε=0.01，即当设计变量矢量迭代更新值与更新前的值的差量的最大值小于ε时，迭代终止。
2 算例和分析 2.1 自重载荷拓扑优化经典案例 图 1为一个二维的简支梁，高40 mm，长120 mm，两端的下方均铰链固定，仅受该梁结构自身的重力作用。图 2为梁结构的有限元模型示意图，单元类型为四节点矩形单元，数量为4 800，厚度为1。
 |
| 图 1 自重载荷作用下的二维简支梁受力结构 Fig. 1 Structural diagram of 2D simply supported beam under self-weight load |
| 图选项 |
 |
| 图 2 待优化结构的有限元模型 Fig. 2 Finite element model of structure to be optimized |
| 图选项 |
对结构进行拓扑优化，体积约束设置为30%，经过45步迭代可以得到最终的优化结果。在迭代过程中，在第1次迭代到第9次迭代，结构发生非常大的改变；迭代到第26次时，优化的结构基本成型，趋于稳定，结构周围分布的中间密度单元逐渐收敛；从第26次迭代到第45次迭代，结构基本稳定，第45次迭代后得到了最终的优化结果，为拱形结构。图 3为二维简支梁结构在自重载荷作用下的拓扑优化过程。
 |
| 图 3 二维简支梁结构拓扑优化过程 Fig. 3 Topology optimization process of 2D simply supported beam structure |
| 图选项 |
对优化出来的结构进行有限元分析，如图 4所示，拱形结构的最大应力为9.230×10-2 MPa。图 5为优化前的简支梁结构，该结构的最大应力为29.48×10-2 MPa。在仅受自重载荷作用下，拱形结构最大应力相对该结构减少了约68.7%。
 |
| 图 4 二维简支梁拓扑优化结果应力分析 Fig. 4 Stress analysis of 2D simply supported beam topology optimization results |
| 图选项 |
 |
| 图 5 优化前二维简支梁结构应力分析 Fig. 5 Stress analysis of 2D simply supported beam structure before optimization |
| 图选项 |
将简支梁重力等效为集中力施加于结构重心位置得到的拓扑优化结构，如图 6所示，图中应力云图为对该结构进行实际自重作用情况下的结构应力分析所得。这个结构在自重载荷作用下最大应力为10.72×10-2 MPa，相对于拱形结构增加了约16%。
 |
| 图 6 自重载荷作用下结构应力分析 Fig. 6 Structural stress analysis under self-weight load |
| 图选项 |
2.2 二维悬臂梁在重力作用下的结构拓扑优化 图 7为一个高40 mm、长120 mm的悬臂梁，固定梁左边的端部，对自重载荷作用下的悬臂梁进行结构拓扑优化。
 |
| 图 7 自重载荷作用下的二维悬臂梁结构示意图 Fig. 7 Schematic diagram of 2D cantilever beam under self-weight load |
| 图选项 |
图 8为自重载荷作用下的悬臂梁结构经过拓扑优化得到的优化结果。在结果中，距离固定端较远的那一端相对较模糊，单元处于中间密度，没有形成清晰明确的材料分布。
 |
| 图 8 自重载荷作用下的二维悬臂梁拓扑优化结果 Fig. 8 Topology optimization results of 2D cantileverbeam under self-weight load |
| 图选项 |
图 9为不同长高比下的悬臂梁进行结构拓扑优化的结果。结果表明，在长高比大于120∶60时，优化结果没有形成完整的结构，在末端约1/3的区域存在中间密度单元，材料分布模糊；在长高比小于120∶60时，优化结果呈现出一个完整的悬臂梁结构，但是材料仍然没有趋于完全收敛，在末端靠近边线的地方分布有一定的中间密度单元。这是因为悬臂梁在长高比较大时，其左右两端的应变能差距会很大，并且右端的应变能会很小，根据式(8)的分析，右端区域的灵敏度相对于左端的小很多，并且大小会非常接近甚至相同，导致材料分布不收敛，出现中间密度单元。另外，仅仅受自重载荷下的悬臂梁结构，优化结果中材料分布应该集中于固定端附近更为合理，这是由于拓扑优化算法中采用了网格过滤方法来控制棋盘格现象和网格依赖性问题，导致优化结果中材料并没有集中分布于固定端附近。根据图 9所示的优化结果，对自重载荷下悬臂梁的拓扑优化设计，优化出来的悬臂梁结构的长度是无法保证的。在对悬臂梁长度有要求的情况下，优化设计得到的结果并不能满足设计需求。
 |
| 图 9 不同长高比下的悬臂梁拓扑优化结果 Fig. 9 Topology optimization results of cantilever beams with different aspect ratios |
| 图选项 |
综合考虑上面的问题，采取在结构末端施加非结构质量的策略。末端非结构质量的添加保证了悬臂梁优化所得到结构的长度，同时可以使悬臂梁末端材料分布收敛，得到清晰的优化结果。施加非结构质量的悬臂梁结构模型如图 10所示。
 |
| 图 10 施加非结构质量的悬臂梁结构示意图 Fig. 10 Schematic diagram of a cantilever beam with non-structural mass |
| 图选项 |
非结构质量的添加实际上改变了实际的工况，对拓扑优化得到的结果必然产生影响。但是通过研究非结构质量的合适大小，可以一定程度减少非结构质量对优化结果的影响。分别对图 10所示悬臂梁进行拓扑优化设计，图 11为在不同大小质量点情况下的拓扑优化结果，r为质量点与设计前悬臂梁结构质量的比值。
 |
| 图 11 不同非结构质量情况下悬臂梁的拓扑优化结果 Fig. 11 Topology optimization results of cantilever beams with different non-structural masses |
| 图选项 |
上面的优化结果表明，在非结构质量大于等于原悬臂梁结构质量4%的时候，悬臂梁结构的拓扑优化可以获得清晰明确的优化结果，收敛的情况有了很大的改善。分别对r=1、r=0.1和r=0.04时优化得到的结果做自重载荷下的有限元分析，其应力云图分别如图 12~图 14所示。如表 1中的3种结构的应力峰值所示，Sr=1>Sr=0.1>Sr=0.04，即非结构质量取到原悬臂梁质量的4%时优化能够获得相对更好的结构。
 |
| 图 12 r=1悬臂梁拓扑优化结果应力分析 Fig. 12 Stress analysis of r=1 cantilever beam topology optimization results |
| 图选项 |
 |
| 图 13 r=0.1悬臂梁拓扑优化结果应力分析 Fig. 13 Stress analysis of r=0.1 cantilever beam topology optimization results |
| 图选项 |
 |
| 图 14 r=0.04悬臂梁拓扑优化结果应力分析 Fig. 14 Stress analysis of r=0.04 cantilever beam topology optimization results |
| 图选项 |
表 1 悬臂梁结构自重载荷下应力峰值 Table 1 Peak stress of cantilever beam structure under self-weight load
| 非结构质量比值 | 峰值应力/MPa |
| 1 | 10.15×10-2 |
| 0.1 | 8.677×10-2 |
| 0.04 | 6.219×10-2 |
表选项
3 结论 采用基于导重准则法和RAMP模型的变密度法对自重载荷作用下的结构进行了拓扑优化研究, 讨论了出现的问题及采取的解决措施。结合理论推导和算例的分析，可以得到以下结论：
1) 通过对灵敏度公式的分析，由于数值化整误差的存在，在结构的低应变能区域容易出现不收敛的情况，导致材料分布模糊，结合算例进行了说明。
2) 自重载荷根据材料分布均匀分配情况下进行拓扑优化相对于将自重等效集中施加于重心位置的情况会获得不同的优化结果，并且得到的结构在自重下性能更优。
3) 在悬臂梁末端添加非结构质量可以有效解决悬臂梁末端材料分布模糊的情况，保证了悬臂梁的有效工作长度，但同时也改变了工况，选择合适大小的非结构质量可以在尽可能减少工况改变的影响的同时获得稳定收敛的优化结果。
4) 通过对比分析，在选择合适大小的非结构质量的情况下优化出来的结构相对于其他大小非结构质量能够获得更好的性能，表明了方法的有效性。
参考文献
| [1] | 赵军鹏, 王春洁. 载荷不确定条件下的结构拓扑优化算法[J]. 北京航空航天大学学报, 2014, 40(7): 959-964. ZHAO J P, WANG C J. Algorithm of structural topology optimization under loading uncertainty[J]. Journal of Beijing University of Aeronautics and Astronautics, 2014, 40(7): 959-964. (in Chinese) |
| [2] | 付志方, 王春洁. 载荷不确定的周期性结构稳健拓扑优化[J]. 北京航空航天大学学报, 2017, 43(4): 747-753. FU Z F, WANG C J. Robust topology optimization of periodic structures under uncertain loading[J]. Journal of Beijing University of Aeronautics and Astronautics, 2017, 43(4): 747-753. (in Chinese) |
| [3] | 贾娇, 程伟, 龙凯. 基于SIMP法的周期性传热材料拓扑优化[J]. 北京航空航天大学学报, 2015, 41(6): 1042-1048. JIA J, CHENG W, LONG K. Topology optimization for periodic thermal conductive material using SIMP method[J]. Journal of Beijing University of Aeronautics and Astronautics, 2015, 41(6): 1042-1048. (in Chinese) |
| [4] | 李明高, 李明, 韩璐, 等. 基于多目标的高速列车隔热结构拓扑优化[J]. 北京航空航天大学学报, 2016, 42(5): 878-884. LI M G, LI M, HAN L, et al. Multi-objective topological optimization of heat insulation structures used in high-speed train[J]. Journal of Beijing University of Aeronautics and Astronautics, 2016, 42(5): 878-884. (in Chinese) |
| [5] | 郭连水, BANDI P, RENAUD J E. 一种大变形多空间域连续体结构拓扑优化方法[J]. 北京航空航天大学学报, 2009, 35(2): 227-230. GUO L S, BANDI P, RENAUD J E. Method of multi-domain topology optimization for continuum structures[J]. Journal of Beijing University of Aeronautics and Astronautics, 2009, 35(2): 227-230. (in Chinese) |
| [6] | 刘丰睿, 赵丽滨, 韩邦成, 等. 磁悬浮控制力矩陀螺框架结构的拓扑优化设计[J]. 北京航空航天大学学报, 2010, 36(4): 455-458. LIU F R, ZHAO L B, HAN B C, et al. Topology optimization design of frame structure for magnetic suspension control moment gyroscope[J]. Journal of Beijing University of Aeronautics and Astronautics, 2010, 36(4): 455-458. (in Chinese) |
| [7] | 陶扬, 王春洁, 付志方. 卫星天线复合材料框架的铺层优化设计[J]. 北京航空航天大学学报, 2017, 43(5): 1031-1037. TAO Y, WANG C J, FU Z F. Layup design optimization of composite frame for satellite antenna[J]. Journal of Beijing University of Aeronautics and Astronautics, 2017, 43(5): 1031-1037. (in Chinese) |
| [8] | 陈树勋, 叶尚辉. 天线结构优化设计的准则法[J]. 固体力学学报, 1984(4): 482-497. CHEN S X, YE S H. Guide-weight method for structural topology optimization[J]. Chinese Journal of Computational Mechanics, 1984(4): 482-497. (in Chinese) |
| [9] | BENDSOE M P. Optimal shape design as a material distribution problem[J]. Structural Optimization, 1989, 1(4): 193-202. DOI:10.1007/BF01650949 |
| [10] | PEDERSEN N L. Maximization of eigenvalues using topology optimization[J]. Structural and Multidisciplinary Optimization, 2000, 20(1): 2-11. DOI:10.1007/s001580050130 |
| [11] | STOLPE M, SVANBERG K. An alternative interpolation scheme for minimum compliance topology optimization[J]. Structural and Multidisciplinary Optimization, 2001, 22(2): 116-124. DOI:10.1007/s001580100129 |
| [12] | BRUYNEEL M, DUYSINX P. Note on topology optimization of continuum structures including self-weight[J]. Structural and Multidisciplinary Optimization, 2005, 29(4): 245-256. DOI:10.1007/s00158-004-0484-y |
| [13] | LOPES R A, STUMP F V, SILVA E C N. Topology optimization of three dimensional structures under self-weight and inertial forces[C]//6th World Congress of Structural and Multidisciplinary Optimization, 2005: 5-6. |
| [14] | LIU X J, LI Z D, WANG L P, et al. Solving topology optimization problems by the guide-weight method[J]. Frontiers of Mechanical Engineering, 2011, 6(1): 136-150. |
| [15] | LIU X J, LI Z D, CHEN X. A new solution for topology optimization problems with multiple loads: The guide-weight method[J]. Science China Technological Sciences, 2011, 54(6): 1505. DOI:10.1007/s11431-011-4334-z |
| [16] | 陈祥, 刘辛军. 基于RAMP插值模型结合导重法求解拓扑优化问题[J]. 机械工程学报, 2012, 48(1): 135-140. CHEN X, LIU X J. Solving topology optimization problems based on RAMP method combined with guide-weight method[J]. Journal of Mechanical Engineering, 2012, 48(1): 135-140. (in Chinese) |
| [17] | 张晖, 刘书田, 张雄. 考虑自重载荷作用的连续体结构拓扑优化[J]. 力学学报, 2009, 41(1): 98-104. ZHANG H, LIU S T, ZHANG X. Topology optimization of continuum structures subjected to self weight loads[J]. Chinese Journal of Theoretical and Applied Mechanics, 2009, 41(1): 98-104. (in Chinese) |
| [18] | XU H Y, GUAN L W, CHEN X, et al. Guide-weight method for topology optimization of continuum structures including body forces[J]. Finite Elements in Analysis and Design, 2013, 75: 38-49. DOI:10.1016/j.finel.2013.07.002 |
| [19] | 许华旸, 关立文, 王立平, 等. 惯性载荷下飞行模拟器大臂结构的拓扑优化[J]. 机械工程学报, 2014, 50(9): 14-23. XU H Y, GUAN L W, WANG L P, et al. Topology optimization for the arm of flight simulator under inertial loads[J]. Journal of Mechanical Engineering, 2014, 50(9): 14-23. (in Chinese) |
| [20] | JAIN N, SAXENA R. Effect of self-weight on topological optimization of static loading structures[J]. Alexandria Engineering Journal, 2018, 57(2): 527-535. DOI:10.1016/j.aej.2017.01.006 |
| [21] | HOLMBERG E, THORE C J, KLARBRING A. Worst-case topology optimization of self-weight loaded structures using semi-definite programming[J]. Structural and Multidisciplinary Optimization, 2015, 52(5): 915-928. |
| [22] | 陈树勋. 精密复杂结构的几种现代设计方法[M]. 北京: 北京航空航天大学出版社, 1992: 226-244. CHEN X S. Several modern design methods for precise and complex structures[M]. Beijing: Beihang University Press, 1992: 226-244. (in Chinese) |
| [23] | 陈树勋. 工程结构系统的分析, 综合与优化设计[M]. 北京: 中国科学文化出版社, 2008: 99-108. CHEN S X. Analysis, synthesis and optimization of engineering structural system[M]. Beijing: China Science Culture Publishing House, 2008: 99-108. (in Chinese) |
