拓扑优化是一种在指定设计域内寻找最优材料分布的结构优化方法. 自Bendsoe和Kikuchi首次提出均匀化方法以来^[1]，拓扑优化，尤其是连续体拓扑优化发展十分迅速. 此后，又逐渐推出了变密度法^[2-3]、水平集法^[4-5]、独立连续映射方法^[6-8]和渐进结构优化方法^[9-10]等拓扑优化方法，并在结构概念设计阶段得到了广泛应用. 目前关于结构拓扑优化的研究大多注重于减小结构的柔度^[11-13]. 但工程实践中更关注结构是否符合强度要求，即结构中的最大应力要低于材料的许用应力，以保证结构安全. 因此，近年来应力相关的拓扑优化问题备受关注. 应力优化问题主要有3大挑战^[14]：应力的非线性行为、应力奇异、应力的局部性质.

应力的非线性行为是指应力水平对结构的变化极为敏感，这种现象在拐角和锯齿形结构等应力集中部位更加明显. 应力的非线性行为使得优化过程不收敛或收敛至局部最优解. Le等^[14]使用变量过滤和灵敏度过滤方法显著改善了应力的非线性特性；Liu等^[15]使用固定网格技术得到了具有光滑边界的拓扑构型，消除了锯齿形边界引起的应力集中.

应力奇异现象主要出现在采用变密度法的优化问题中. 优化时由于低密度单元具有与其密度不相符的高应力，使得优化算法不能有效更新设计变量. 采用应力插值模型^[16-17]可以有效克服应力奇异现象. 离散拓扑优化方法，如水平集法和渐进结构优化方法，因为不存在中间密度单元，所以不会出现应力奇异现象，且拓扑优化结果边界清晰. 水平集法的缺点是较为依赖初始设计. Huang等^[18-19]提出的双向渐进结构优化（bi-directional evolutionary structural optimization，BESO）方法采用启发式算法更新设计变量，避免使用数学规划方法求解优化问题，降低了计算代价，优化效率较高，因此BESO方法是一种极具潜力的拓扑优化方法.

应力的局部性质是指应力仅能表示单元自身的力学性能，不能反映结构的整体性能. 当采用传统方法进行应力优化时，需要构建与设计域中单元数目相当的约束，此举显著增加了计算时长，难以应用于实际问题. 目前主要使用应力全局化方法解决应力局部性质带来的计算量大的问题. 其中一种应力全局化方法采用凝聚函数逼近结构中的最大应力，建立应力的全局函数. 典型的凝聚函数有Kreisselmeier-Steihauser（K-S）函数和P范数. 王选等^[20-21]采用K-S函数，Liu等^{[15, 22]}采用P范数研究应力优化问题，优化效果显著. 另外一种应力全局化方法是隋允康等^[23-24]提出的畸变比能法，该方法基于第四强度理论，将应力优化等效为畸变比能优化，并进一步转化为应变能优化. 畸变比能法极大提高了应力优化的计算效率，可应用于设计变量较多的应力优化问题. 但该方法得到的结果偏于保守，无法消除结构中的应力集中区域.

针对体积约束下的应力极小化连续体结构拓扑优化问题，提出一种改进的BESO方法. 该方法用一对离散设计变量0（或接近0的极小值）和1来表示设计域中单元的有无，以避免应力奇异现象，得到边界清晰的拓扑构型. 在Liu等^{[15, 24-25]}的启发下，提出使用K-S函数逼近结构中的最大畸变比能，建立应力拓扑优化模型. 推导单元灵敏度计算方法，并基于商业有限元软件ABAQUS二次开发拓扑优化程序，提高了应力优化的计算效率. 利用BESO方法逐步更新设计变量，克服应力优化的非线性问题，使结构稳定进化至最优构型. 最后使用3个典型的数值算例验证该方法的有效性.

根据第四强度理论：当结构内一点处的畸变比能（也称畸变能密度）达到了材料的极限值时，该点处的材料就发生塑性屈服. 单元i的畸变比能定义为

(1) $\begin{split} {s_i} =& \frac{{1 + \upsilon }}{{6E}}\left[ {{{\left( {{\sigma _1} - {\sigma _2}} \right)}^2} + {{\left( {{\sigma _2} - {\sigma _3}} \right)}^2}} \right. + \\ &\left. {{\rm{ }}{{\left( {{\sigma _3} - {\sigma _1}} \right)}^2}} \right] = \frac{1}{{6G}}\sigma _{{\rm{vm}},\;i}^2{\rm{ }}. \end{split} $

式中：υ为泊松比，E为材料的弹性模量，G为材料的剪切弹性模量，σ₁、σ₂、σ₃分别表示单元i在3个方向上的主应力，σ_{vm, i}为单元i的von Mises应力. 由式（1）可知，最小化畸变比能等效于最小化von Mises应力^{[8, 25-26]}.

应力优化时随着材料分布的变化，结构中具有最大畸变比能的单元不断改变，若以结构中的最大畸变比能最小为目标，会导致优化过程不稳定. 为了克服畸变比能的局部性质，采用K-S函数对结构中的最大畸变比能进行全局化近似^{[8, 24-25]}：

(2) ${s^{{\rm{KS}}}} = \frac{1}{\mu }\ln\; \left( {\sum\limits_{i = 1}^n {{{\rm{e}}^{\mu {s_i}}}} } \right).$

式中：s^KS为畸变比能的全局函数，μ为常数，n为设计域中的单元个数. 当μ 趋于无穷大时，全局函数s^KS退化为结构中的最大畸变比能，即s^KS→max （s_i），i=1，2，···，n，s^KS失去全局性.

本研究讨论的是提高结构强度的拓扑优化方法，为了保证结构的安全，须最小化结构中的最大畸变比能. 1.2节建立了畸变比能的全局函数，结合BESO方法，拓扑优化模型为

(3) $\left. {\begin{split} &{\min :{s^{{\rm{KS}}}}};\\ &{{\rm{s}}.{\rm{t}}.\;V = {V^*} = \displaystyle\sum\limits_{i = 1}^n {{x_i}{v_i}} },\\ &\;\;\;\;\;\;\;\;{{x_i} = {x_{\min }}\;{\rm{or}}\;{\rm{1}}.} \end{split}} \right\}$

式中：V为当前结构的体积分数；V^*为优化前设定的约束体积分数；v_i为单元i在初始结构中的体积分数；x_i为拓扑设计变量，代表单元i的状态. 当x_i=1时，表示单元i为实单元；当x_i=x_min时，表示单元i为空单元. 取x_min=0.001，以避免结构整体刚度矩阵奇异.

使用离散设计变量，目标函数连续性不足，无法推导单元灵敏度. 为了增加目标函数的连续性，并避免重新划分网格，采用文献[14]介绍的方法，对弹性模量和单元应力进行插值：

(4) ${E_i} = {x_i}^p{E_0},$

(5) ${{{\sigma }}_i} = {{{D}}_i}{{{B}}_i}{{{u}}_i} = {x_i}^q{{{D}}_0}{{{B}}_i}{{{u}}_i}. $

式中：E_i、E₀分别为单元i和实体材料的弹性模量；σ_i为单元i的应力，对于平面应力单元，σ_i=[σ_x，σ_y，τ_xy]^T，σ_x、σ_y、τ_xy分别为单元中心处沿x向和y向的正应力以及切应力；D_i、D₀分别为单元i和实体材料的本构矩阵；B_i为单元i的应变矩阵；u_i为单元i所属节点的位移矢量；p、q为常数，本研究取p=3，q=0.5.

由式（2）可知，畸变比能全局函数s^KS对设计变量x_j的导数可表示为^[24-25]

(6) $\frac{{\partial {s^{{\rm{KS}}}}}}{{\partial {x_j}}} = {{\displaystyle\sum\limits_{i = 1}^n {\frac{{\partial {s_i}}}{{\partial {x_j}}}} {{\rm{e}}^{\mu {s_i}}}}}\Bigg/{{\displaystyle\sum\limits_{i = 1}^n {{{\rm{e}}^{\mu {s_i}}}} }}.$

根据链式求导法则，单元i的畸变比能s_i对设计变量x_j的导数可表示为

(7) $ \frac{{\partial {s_i}}}{{\partial {x_j}}} = \frac{{\partial {s_i}}}{{\partial {\sigma _{{\rm{vm}},i}}}}\frac{{\partial {\sigma _{{\rm{vm}},i}}}}{{\partial {{{\sigma }}_{{i}}}}}\frac{{\partial {{{\sigma }}_{{i}}}}}{{\partial {x_j}}}. $

畸变比能对von Mises应力的导数根据式（1）可得

(8) $ \frac{{\partial {s_i}}}{{\partial {\sigma _{{\rm{vm}},\;i}}}} = \frac{{{\sigma _{{\rm{vm}},\;i}}}}{{3G}}. $

σ_{vm, i}可表示为

(9) $ {\sigma _{{\rm{vm}},\;i}} = {\left( {{{\sigma }}_i^{\rm{T}}{{V}}{{{\sigma }}_i}} \right)^{1/2}}. $

式中：V为转换矩阵. 对于平面应力单元，

(10) ${{V}} = \left[ {\begin{array}{*{20}{c}} {1.0}&{ - 0.5}&{0} \\ { - 0.5}&{1.0}&{0} \\ {0}&{0}&{3.0} \end{array}} \right].$

因此，σ_{vm, i}对σ_i的导数为

(11) $\frac{{\partial {\sigma _{{\rm{vm}},\;i}}}}{{\partial {{{\sigma }}_i}}} = \sigma _{{\rm{vm}},\;i}^{ - 1}{{\sigma }}_i^{\rm{T}}{{V}}. $

由式（5）可知，σ_i对x_j的导数为

(12) $\frac{{\partial {{{\sigma }}_i}}}{{\partial {x_j}}} = \frac{{\partial {x_i}}}{{\partial {x_j}}}{x_i}^{q - 1}{{{D}}_0}{{{B}}_i}{{{u}}_i} + {x_i}^q{{{D}}_0}{{{B}}_i}\frac{{\partial {{{u}}_i}}}{{\partial {x_j}}}.$

结合式（5）、（9）、（11）、（12）可知，σ_{vm, i}对x_j的导数为

(13) $\begin{split} \frac{{\partial {\sigma _{{\rm{vm}},\;i}}}}{{\partial {x_j}}} =& \frac{{\partial {\sigma _{{\rm{vm}},\;i}}}}{{\partial {{{\sigma }}_i}}}\frac{{\partial {{{\sigma }}_i}}}{{\partial {x_j}}} = \frac{{{\sigma _{{\rm{vm}},\;i}}}}{{{x_i}}}\frac{{\partial {x_i}}}{{\partial {x_j}}} + \\&\sigma _{{\rm{vm}},\;i}^{ - 1}{x_i}^q{{\sigma }}_i^{\rm{T}}{{V}}{{{D}}_{\rm{0}}}{{{B}}_i}\frac{{\partial {{{u}}_i}}}{{\partial {x_j}}}{\rm{ }}{\rm{.}} \end{split}$

线性结构静力平衡方程可表示为

(14) $ {{KU}} = {{F}}. $

式中：K为结构的整体刚度矩阵，U为结构中所有节点的位移矢量，F为外部施加的载荷矢量. 由于外力在优化过程中保持不变，式（14）两侧对设计变量x_j求导数，可得

(15) $ \frac{{\partial {{K}}}}{{\partial {x_j}}}{{U}} + {{K}}\frac{{\partial {{U}}}}{{\partial {x_j}}} = {{0}}. $

引入维度转换矩阵C_i，该矩阵可将结构整体位移矢量U转变为式（12）中的单元所属节点的位移矢量u_i，即：

(16) $ {{{u}}_i} = {{{C}}_i}{{U}}. $

将式（15）和（16）代入式（13），有

(17) $ \frac{{\partial {\sigma _{{\rm{vm}},\;i}}}}{{\partial {x_j}}} = \frac{{{\sigma _{{\rm{vm}},\;i}}}}{{{x_i}}}\frac{{\partial {x_i}}}{{\partial {x_j}}} - \sigma _{{\rm{vm}},\;i}^{ - 1}{x_i}^q{{\sigma }}_i^{\rm{T}}{{V}}{{{D}}_0}{{{B}}_i}{{{C}}_i}{{{K}}^{ - 1}}\frac{{\partial {{K}}}}{{\partial {x_j}}}{{U}}. $

将式（8）和（17）代入式（7），得

(18) $ \frac{{\partial {s_i}}}{{\partial {x_j}}} = \frac{{\sigma _{{\rm{vm}},\;i}^2}}{{3G{x_i}}}\frac{{\partial {x_i}}}{{\partial {x_j}}} - \frac{1}{{3G}}{x_i}^q{{\sigma }}_i^{\rm{T}}{{V}}{{{D}}_0}{{{B}}_i}{{{C}}_i}{{{K}}^{ - 1}}\frac{{\partial {{K}}}}{{\partial {x_j}}}{{U}}. $

将式（18）代入式（6），畸变比能全局函数s^KS对x_j的导数为

(19) $ \frac{{\partial {s^{{\rm{KS}}}}}}{{\partial {x_j}}} = \frac{{\dfrac{{\sigma _{{\rm{vm}},\;j}^2}}{{{x_j}}}{{\rm{e}}^{\mu {s_j}}} - \displaystyle\sum\limits_{i = 1}^n {{x_i}^q} {{\sigma }}_i^{\rm{T}}{{V}}{{{D}}_0}{{{B}}_i}{{{C}}_i}{{{K}}^{ - 1}}\frac{{\partial {{K}}}}{{\partial {x_j}}}{{U}}{{\rm{e}}^{\mu {s_i}}}}}{{3G\displaystyle\sum\limits_{i = 1}^n {{{\rm{e}}^{\mu {s_i}}}} }}. $

为了方便求解式（19），引入伴随变量λ，使

(20) $ {{K\lambda }} = \sum\limits_{i = 1}^n {{x_i}^q} {\left( {{{{D}}_0}{{{B}}_{{i}}}{{{C}}_{{i}}}} \right)^{\rm{T}}}{{V}}{{{\sigma }}_i}{{\rm{e}}^{\mu {s_i}}}. $

则式（19）可简化为

(21) $ \frac{{\partial {s^{{\rm{KS}}}}}}{{\partial {x_j}}} = \frac{{\dfrac{{\sigma _{{\rm{vm}},\;j}^2}}{{{x_j}}}{{\rm{e}}^{\mu {s_j}}} - {{{\lambda}} ^{\rm{T}}}\dfrac{{\partial {{K}}}}{{\partial {x_j}}}{{U}}}}{{3G\displaystyle\sum\limits_{i = 1}^n {{{\rm{e}}^{{\mu s _{{i}}}}}} }}. $

整体刚度矩阵K是由所有单元刚度矩阵经维度转换，并按照直接刚度法原理组装而成，即有

(22) $ {{K}} = \sum\limits_{i = 1}^n {{{C}}_i^{\rm{T}}} {{{k}}_i}{{{C}}_i}. $

单元j的刚度矩阵仅与本单元对应的设计变量x_j有关，因此K对x_j的导数为

(23) $ \frac{{\partial {{K}}}}{{\partial {x_j}}} = \sum\limits_{i = 1}^n {{{C}}_i^{\rm{T}}} \frac{{\partial {{{k}}_i}}}{{\partial {x_j}}}{{{C}}_{{i}}} = {{C}}_i^{\rm{T}}{{{k}}_{j,(0)}}{{{C}}_{{i}}}. $

式中：k_j为单元j的刚度矩阵，k_j,(0)是单元j为实单元时的刚度矩阵.

将式（23）代入式（21），可求出s^KS对设计变量x_j的导数为

(24) $ \frac{{\partial {s^{{\rm{KS}}}}}}{{\partial {x_j}}} = \frac{{\dfrac{{\sigma _{{\rm{vm}},\;j}^2}}{{{x_j}}}{{\rm{e}}^{\mu {s_j}}} - {{\lambda }}_j^{\rm{T}}{{{k}}_{j,(0)}}{{{u}}_j}}}{{3G\displaystyle\sum\limits_{i = 1}^n {{{\rm{e}}^{{\mu s_{{i}}}}}} }}. $

按照BESO方法删除低灵敏度单元的设计变量更新策略，将式（24）添加负号以后，可得到单元j的灵敏度为

(25) $ \alpha _j^0 = - \frac{{\partial {s^{{\rm{KS}}}}}}{{\partial {x_j}}} = \frac{{{{\lambda }}_j^{\rm{T}}{{{k}}_{j,(0)}}{{{u}}_j} - \dfrac{{\sigma _{{\rm{vm}},\;j}^2}}{{{x_j}}}{{\rm{e}}^{\mu {s_j}}}}}{{3G\displaystyle\sum\limits_{i = 1}^n {{{\rm{e}}^{\mu {s_i}}}} }}. $

BESO方法是一种在优化过程中逐步增加和删除一部分单元使结构进化至最优材料布局的拓扑优化方法. 该方法采用的灵敏度过滤、修正方法和优化设计流程，可避免拓扑优化中常见的棋盘格现象和网格依赖性，并稳定优化设计过程^[18].

应力优化时，由于应力的非线性行为，相邻单元的灵敏度数值可能相差较大. 直接使用式（25）计算得到的单元灵敏度更新设计变量，材料分布可能不连续，继而导致优化不收敛. 为了将单元灵敏度更加均匀地分布到临近单元上，引入一个无实际物理意义的概念—节点灵敏度，节点j的灵敏度定义为

(26) $ \alpha _j^n = {{\displaystyle\sum\limits_{i = 1}^M {{v_i}} \alpha _i^0}}\left/{{{\displaystyle\sum\limits_{i = 1}^M {{v_i}} }} }\right. . $

式中：M为与节点j相连的单元个数. 为了进一步平滑灵敏度，以单元i的形心为圆心，将半径r_min的子空间内的所有节点灵敏度整合，作为单元i的灵敏度：

(27) $ {\alpha _i} = {{\displaystyle\sum\limits_{j = 1}^l {{\omega _{ij}}} \alpha _j^n}} \left/{ {{\displaystyle\sum\limits_{j = 1}^l {{\omega _{ij}}} }}}\right.. $

式中：l为子空间内的节点个数，ω_ij为节点j相对于单元i的权重，可表示为

(28) $ {\omega _{{{i\!j}}}} = {r_{\min }} - {r_{{{i\!j}}}}. $

式中：r_ij为单元i的形心与节点j之间的距离. 该灵敏度过滤方法使得空单元也具有一定的灵敏度，从而有可能恢复为实单元. 显然，位于高灵敏度区域的空单元恢复概率更大.

BESO方法采用离散设计变量，因此更新设计变量（即材料的重新分布）可能会使目标函数在迭代过程中出现振荡. 在处理应力优化问题时，这种现象更加明显. 为此，将当前迭代的灵敏度与前一次迭代的灵敏度取均值，作为当前迭代的最终灵敏度，以稳定优化过程，即

(29) $ {\hat \alpha _i} = \left( {{{{\alpha _{i,\;k}} + {\alpha _{i,\;k - 1}}}}} \right) \left/{{2} }\right. . $

式中：α_{i, k}和α_{i, k−1}分别为单元i第k次迭代及第k−1次迭代的灵敏度.

第k次迭代在增加和删除单元之前，首先须确定第k+1次迭代的体积分数：

(30) $ {V_{k + 1}} = \max \;\left( {{V^*},{V_k} (1 - R_{\rm{E}})} \right). $

式中：R_E为体积分数的进化率. 确定体积分数以后，将所有单元按灵敏度数值从大到小排列. 对于实单元，如果灵敏度α_i≤α_del，k，则转为空单元；对于空单元，如果灵敏度α_i≥α_add，k，则转为实单元. 单元增加阈值α_add，k和删除阈值α_del，k按照下列步骤确定.

1）取α_N=α_add，k=α_del，k. α_N为按从大到小次序排列的第N个单元的灵敏度，N=V_k+1×n.

2）计算增加单元比例R_A根据初步确定的单元增加阈值α_add，k，记大于α_add，k的空单元个数为t，取R_A=t/n. 若R_A小于预先设定的最大单元增加比值R_A，max，则跳过3），否则按照3）重新确定α_add，k.

3）计算α_add，k. 为了限制每次迭代增加单元的个数，每次迭代转为实单元的空单元个数T最多为T=R_A，max×n. 将所有空单元的灵敏度从大到小排列，取前T个空单元转为实单元，即将这些空单元的设计变量由x_min更新为1. 在增加单元个数确定以后，由于下一次迭代的体积分数已经确定，需要删除的单元个数亦即确定. 按照单元灵敏度从小到大的顺序删除相应个数的实单元，即将这些实单元的设计变量由1更新为x_min.

设定最大单元增加比R_A，max的意义在于稳定优化过程. 如果直接使用第1）步方法设定的增加阈值和删除阈值更新设计变量，该次迭代可能会增加较多的单元. 与此同时会删除较多的单元，以维持下一步的体积分数. 但材料分布可能因此出现剧烈变化，难以收敛至最优解.

BESO方法迭代步骤可归纳如下.

1）构建待优化结构的有限元模型，取设计域内所有单元对应的设计变量初始值为1；

2）指定优化设计参数，如进化率R_E、灵敏度过滤半径r_min、最大单元增加比R_{A, max}、参数μ、p和q等；

3）执行有限元分析；

4）按照式（25）计算初始单元灵敏度；

5）根据式（26）、（27）、（29）确定最终单元灵敏度；

6）使用3.3节介绍的方法更新设计变量；

7）根据设计变量更新有限元模型；

8）重复步骤3）~7），直至达到设定的收敛条件.

使用3个典型的应力拓扑优化算例来验证所提出方法的有效性. 所有算例的材料弹性模量E=1 MPa，泊松比υ=0.3，R_E=2%，最大单元增加比R_{A, max}=1%. 由于应力的非线性行为，即使结构的构型已经基本确定，个别单元的增加和删除仍会使得结构中的最大应力在较小范围内波动，文献[18]使用的收敛准则不再适用. 本研究在结构体积分数达到设定约束值且材料分布稳定以后，选择一个应力最小的结果作为拓扑优化的最终结果^[22].

算例1

算例1为常用于应力优化的L型梁^{[14-15, 22]}，结构尺寸及边界条件如图1所示. L型梁顶端固定，在右上角点施加一个F=4 N的垂向载荷，为避免应力集中，该力均匀地分布在5个节点上. 使用四节点矩形平面应力单元来离散设计域，单元尺寸为1 mm，灵敏度过滤半径r_min=3 mm. 约束体积分数V^*=0.5，初始结构的最大von Mises应力σ_max=2.03 MPa，出现在L型梁的拐角处.

优化时，μ分别取2、5、10，得到的拓扑优化结果及相应的von Mises应力云图如图2（a）~（c）所示. 在3种情况下，最终结果的σ_max分别为1.71、1.55和1.33 MPa. 相较于初始最大应力，分别下降了15.8%、23.7%和34.5%，优化效果显著. 图2（d）为柔度最小化设计的优化结果和von Mises应力云图，与μ=2得到的结果构型相近，应力集中效应未得到缓解. 由图2（a）~（c）可知，当μ 增大时，结构中的应力集中区域逐渐消失，结构中的应力分布更加均匀. 当μ=10时，结构中的大部分单元处于同一应力水平，优化效果最好. 继续增大μ，畸变比能全局函数的非线性程度增大，病态现象严重. 当μ >10时，材料分布开始出现不连续现象，优化过程不收敛. 这说明通过合理设置参数 μ，可以有效缓解结构中的应力集中效应，提高结构强度.

L型梁优化的体积分数及σ_max迭代历史如图3所示。图中，C表示柔度最小化设计，n为迭代次数. 由图3可知，当μ=2时，在应力优化初始阶段，由于应力集中未被消除，σ_max下降不明显. 第27次迭代后，随着应力集中部位逐渐被删除，σ_max逐步下降，且下降过程平稳. 当μ =5、10时，由于L型梁的拐角在优化初始阶段即被删除，并代之以圆弧形结构，应力集中被消除，σ_max显著下降. 当结构体积分数达到设定值时，尽管拓扑构型保持不变，但局部仍在增加和删除单元，而应力对拓扑变化较为敏感，因此σ_max仍在一定范围内波动. 其中，由于目标函数的非线性程度增大，μ =10时的最大应力迭代历史较μ =5时的迭代历史波动明显. 因为柔度最小化设计无法消除应力集中的影响，所以当结构体积分数下降时，σ_max还会略有上升.

变密度法和水平集法一般不设置体积约束或将体积分数作为目标函数，但BESO方法需要设置体积约束，以保证优化稳定收敛^[18]. 为研究体积分数对拓扑优化的影响，当μ=10时，V^*分别取0.4、0.2对L型梁进行优化，得到的结果如图4所示. 由图2（c）、4（a）可知，当V^*=0.5、0.4时，最优结果的材料分布基本一致. 当V^*=0.2时，尽管材料的分布趋势与图2（c）、4（a）所示结果一致，但因为体积分数较小，L型梁左下角只保留了一根斜撑，所以体积分数的取值会影响拓扑优化结果. 随着体积分数下降，结构中的应力水平逐步上升，图4（a）、（b）所示结构的最大应力分别为1.74、2.60 MPa.

如图5所示为μ =10时L型梁的进化过程. 观察发现，在优化初始阶段，L型梁拐角处的单元即被移除，σ_max从2.03 MPa下降至1.85 MPa. 随着材料的不断移除，σ_max下降幅度逐渐减小. 在第60次迭代时，σ_max下降至1.33 MPa，相比于柔度最小化设计的2.1 MPa，下降了36.7%.

Le等^[14]使用变密度法和数学规划方法处理相同问题时，需要迭代110次. 使用本研究提出的方法，在第60次迭代时，结构中的最大应力和材料分布即达到稳定，而拓扑构型与Le等^[14]的结果一致，优化效率显著提升. Xia等^[22]研究应力优化问题时，除了使用灵敏度过滤方法以外，还使用了变量过滤，以稳定优化过程. 本研究提出的方法仅使用灵敏度过滤，避免了变量过滤运算带来的额外计算量. Liu等^[15]使用代理设计变量以得到光滑的边界，但该方法同样会增加计算代价.

算例2

算例2为底端带有缺口的MBB梁，结构尺寸及边界条件如图6所示. 在该梁顶端中部施加一个F=10 N的垂向载荷，为了避免应力集中，该力均匀地分布在梁上端中部的11个节点上. 使用四节点矩形平面应力单元来离散设计域，单元尺寸为1 mm，灵敏度过滤半径r_min=3 mm，约束体积分数V^*=0.5. 初始结构的最大von Mises应力为2.49 MPa，出现在梁底端的缺口处.

如图7（a）~（d）所示为μ=2、5、10以及柔度最小化设计的优化结果及von Mises应力云图. μ =2时得到的优化结果构型与柔度最小化优化结果构型相似，MBB梁缺口处的应力集中并没有消除. 随着μ逐渐增大，优化结果中的应力水平渐趋均匀，且应力集中现象得到抑制. 当μ >10时，优化过程变得不稳定，结果不连续. 原因在于当 μ 取较大值时，目标函数退化为结构中的最大畸变比能. 但应力优化时，具有最大畸变比能的单元不断改变，增加和删除单元的位置也不断改变，材料分布变化较为剧烈，因此优化过程出现振荡.

如图8所示为μ =10时MBB梁的进化过程. 优化初始阶段，梁缺口处的单元被剔除，并逐渐形成光滑的边界，缓解了应力集中现象. 在第40~70次迭代，拓扑构型稳定，结构中的最大应力变化不大. 第70次迭代时，最大von Mises应力下降至1.61 MPa，较初始结构最大应力下降了35.4%. 使用变密度法^[14]和水平集法^[5]得到相同拓扑构型时，需要迭代上百次乃至上千次，计算效率较低.

算例3

算例3为三维悬臂梁，其尺寸及边界条件如图9所示. 悬臂梁根部固定，在右侧上端施加一个F=25 N的垂向载荷. 为了避免应力集中，该力均匀地分布在悬臂梁右上角的25个节点上. 使用8节点六面体单元来离散设计域，单元尺寸为1 mm，灵敏度过滤半径r_min=3 mm，约束体积分数V^*=0.3. 初始结构中σ_max=2.53 MPa，出现在加载区域.

如图10（a）~（d）所示为μ =2、5、10以及柔度最小化设计的拓扑优化结果及von Mises应力云图. 随着μ 逐渐增大，结构中的最大应力逐渐减小，应力分布更均匀. μ =10时的最优设计相比于柔度最小化设计，σ_max下降了6.8%.

如图11所示为μ=10时三维悬臂梁的进化过程. 随着材料的不断移除，结构中的最大应力逐渐下降，在第50次迭代以后，结构中的最大应力无明显变化.

本研究提出的考虑体积约束的应力极小化渐进结构优化方法，设计域中没有中间密度单元，避免了应力奇异，并可以得到边界清晰的拓扑优化结果. 畸变比能全局化策略显著提升了应力优化问题的计算效率. BESO方法使用的灵敏度过滤和修正方法克服了应力的非线性行为，使优化过程更加稳定.

数值算例结果表明，使用不同的K-S凝聚函数参数μ 可以得到不同的优化结果. 当μ 值较小时，得到的优化结果与柔度最小化设计结果相似. 适当增大μ 值，结构中的最大应力逐渐减小，关键区域的应力集中效应得到缓解，结构中的应力分布更加均匀. 但当μ 取值过大时，畸变比能全局函数的非线性程度增大，并退化为结构中的最大畸变比能，全局函数失去了全局性. 材料分布因此变化剧烈，优化过程无法收敛. 研究表明，设定的约束体积分数不同，拓扑优化结果也不同，即优化结果依赖体积分数. 为克服体积分数的影响，下一步将研究以体积分数作为目标函数的应力约束优化问题.

所提出的方法可以解决二维和三维结构的应力拓扑优化问题，而且当单元数量较多时，仍能保持较高的计算效率，具有实际工程应用价值.