CN103049606B - 一种柔顺机构0-1拓扑图提取方法 - Google Patents

Abstract

本发明公开了一种柔顺机构0-1拓扑图提取方法，通过这种不断缩小拓扑优化模型可利用材料域的方式逐渐逼近最佳解，最终提取出理想的柔顺机构0-1拓扑图。本发明可有效克服基于SIMP方法松弛设计变量时所产生的中间单元问题，提取出理想的柔顺机构0-1拓扑图；它完全不同于传统渐进拓扑优化方法，具有连续设计变量，灵敏度连续性好，而且概念简单，方法易于执行，不易导致棋盘格现象；克服了传统方法的参数敏感性和网格依赖性的缺陷，具有较好的收敛稳定性。

Description

一种柔顺机构0-1拓扑图提取方法

技术领域

本发明涉及柔顺机构拓扑优化设计中的拓扑图提取技术领域，特别涉及一种柔顺机构0-1拓扑图提取方法。

背景技术

随着微机电***、微加工和微操作以及新材料等领域的迅速发展，柔顺机构的设计已经成为目前国内外机构学领域的研究热点。

采用拓扑优化方法研究柔顺机构的设计只需给定设计域和指定输入输出位置，无须从一个已知的刚性机构出发，且所得的机构具有分布式柔性的优越性能，因而引起了极大的重视。

这种方法通常是以有限元分析为基础，在拓扑优化的初始阶段，首先将设计域离散成一定数量的有限元网格，再利用优化方法确定单元材料的保留与删除，以满足预定的目标和约束条件。

在优化结果中，单元密度的理想取值为0或1，当单元密度取值为0时，表示该单元被删除，单元密度取值为1时，该单元被保留，于是，由高密度单元连接构成机构的拓扑图。但是，这类整数规划问题往往很难求解。因此，人们常采用松弛法，如均匀化方法和相对密度法，用连续变量的优化模型逼近原离散变量的优化模型，于是，拓扑设计变量可以在区间[0,1]内连续取值，使优化结果中出现中间密度单元。这一现象与工程实际不符。因此，在柔顺机构拓扑优化设计中，必须发展一种有效的方法，以克服中间单元问题，提取出清晰的柔顺机构0-1拓扑图。目前，在包括柔顺机构在内的拓扑优化领域中，主要包括以下几种方法。

一是采用形如ρ^P,(P＞1)的惩罚形式，以减小中间密度单元的刚度影响，然而，这种方法难以彻底去除中间单元。另一种方法是阈值法，但由于其阈值选择是启发式的，因而必然影响到方法的执行效果。此外，又发展了在优化模型的目标函数或约束条件中追加显式的惩罚函数方法，但这种额外的惩罚项也可能给优化软件寻找可行解造成困难。类似的惩罚函数法还有混合的SINH方法和凹约束方法，不过，这两种方法也依然不能提取出满意的0-1拓扑图。另外，基于图像的过滤技术也被尝试来克服中间单元问题，然而，由于其灵敏度计算复杂，使得计算效率大大降低，且所得的拓扑图具有网格依赖性。

近来，又提出了一种修改最佳准则表达形式的启发式方法，但这种方法的最佳拓扑图具有参数敏感性。除了上述方法之外，渐进方法也是一种常用的去除中间单元的方法。它是通过逐渐地去除无效材料来获得最佳解，其优点在于不需引入太多的数学理论，不过由于该方法属于整数规划范畴，灵敏度连续性差，且被去除的单元不能再被利用，因而计算效率低，并易于导致结果不收敛。在此基础上，又发展了双向渐进方法，这种方法旨在去除无效材料的同时，也能添加有效材料，以改进收敛效果，但是这一方法本质上仍是整数规划，且无法精确评估所添加的材料单元的灵敏度。另外，还有模拟退火和SIMP相结合的方法以及单元连结参数化方法等，从它们的优化结果来看，也依然无法提取出理想的0-1拓扑图。并且，上述方法有着易于导致棋盘格的共同缺点。

此外，基于水平集的拓扑优化方法也是解决中间单元问题的一种有效方法，但是这种方法也具有初始敏感性、不能生成新孔、计算效率低和难以收敛到不光滑的角点等缺陷，虽然目前已提出一些改进的方法，但这些问题尚未完全得到很好的解决。

发明内容

本发明的发明目的是针对现有柔顺机构拓扑图提取方法的技术不足，提供一种柔顺机构0-1拓扑图提取方法。

本发明首先建立柔顺机构第一层拓扑优化模型；然后，重复有限元分析、灵敏度过滤及基于最佳准则的设计变量更新过程，直至获得优化模型的拓扑优化最佳结果；之后，采用渐进方法，删除一定数量低密度单元，所有保留单元的集合构成下一层优化模型的可利用材料域，其中的高密度单元在下一层优化中恒定不变，同时，采用体积约束延拓方法，得到下一层优化模型的体积比，进而建立下一层柔顺机构拓扑优化模型，并对其进行优化。通过这种不断缩小拓扑优化模型可利用材料域的方式逐渐逼近最佳解，最终提取出理想的柔顺机构0-1拓扑图。

为实现上述发明目的，本发明采用的技术方案为：

提供一种柔顺机构0-1拓扑图提取方法，包括如下步骤：

1)：建立柔顺机构第一层拓扑优化模型；

2)：对柔顺机构拓扑优化模型进行优化，步骤如下:

2-1)：对柔顺机构拓扑优化模型进行有限元分析，并得出体积约束和优化目标的灵敏度；

2-2)：对优化目标的灵敏度进行过滤，以消除拓扑图中的棋盘格；

2-3)：基于最佳准则更新设计变量；

2-4)：重复步骤2-1)至2-3)，直到优化迭代收敛；

3)：判断柔顺机构拓扑优化最大迭代数或前后两次迭代的单元密度变化最大值是否小于一个阈值，若满足则终止循环并输出结果，否则继续执行以下步骤；

4)：建立下一层柔顺机构拓扑优化模型,

5)：重复步骤2-1)至4)，如此不断缩小拓扑优化模型的可利用材料域，直至提取出满足目标体积V^*的柔顺机构0-1拓扑图。

优选地，以Ω₁为设计域，其为柔顺机构初始拓扑优化模型可利用材料域；P_i和P_o分别为柔顺机构载荷输入点和位移输出点；F_in和F_d分别为柔顺机构输入载荷和沿输出位移方向的虚拟单位载荷；k_in和k_out分别为柔顺机构输入和输出弹簧刚度；该柔顺机构的目标体积比为θ^*，初始阶段的拓扑优化层数为L＝1，所述初始阶段的拓扑优化体积比为将设计域离散为N_x行N_y列，柔顺机构的应变能和互应变能如下：

E_s＝∫_Ωε(u)^TDε(u)dΩ＝U^TKU；

$E_{\mathrm{ms}}={\∫}_{\mathrm{\Ω}}{\mathrm{\ϵ}\left(u_d\right)}^T\mathrm{D\ϵ}\left(u\right)\mathrm{d\Ω}=U_d^T\mathrm{KU};$

式中，E_ms是***的互应变能，互应变能越大则表明***的柔性越大；E_s是***的应变能，应变能越小则表明***的刚度越大；D为弹性矩阵，K是***整体刚度矩阵，U是F_in作用下的节点位移向量，U_d是F_d作用下的节点位移向量，ε(u)和u是设计域内任一点在载荷F作用下的应变和弹性变形，ε(u_d)和u_d是设计域内任一点在载荷F_d作用下的应变和弹性变形；

为使柔顺机构既有足够大的刚度又有足够大的柔性，通过多目标优化而得到柔顺机构的应变能和互应变能的关系如下：

式中，Min代表最小值；

首先，基于相对密度法松弛设计变量，使柔顺机构单元密度在0-1范围内取值，即：

0＜ρ_min≤ρ_i≤ρ_max＝1,i＝1,2,…,N₁；

式中，ρ_i是单元i的密度，ρ_min是单元密度下限，ρ_max是单元密度上限，N₁是自然数，表明初始阶段可利用材料域的单元数，且N₁＝N_xN_y；

柔顺机构拓扑优化模型的整体刚度矩阵是：

$\begin{array}{c}K=\underset{i=1}{\overset{N_1}{\mathrm{\Σ}}}{K}_i\\ =\underset{i=1}{\overset{N_1}{\mathrm{\Σ}}}{\∫}_{V_e}{\mathrm{\ρ}}_i^P{B}^T\mathrm{DBdV}\\ =\underset{i=1}{\overset{N_1}{\mathrm{\Σ}}}{\mathrm{\ρ}}_i^P{K}^0\end{array};$

式中，K_i是单元i的刚度矩阵，V_e是任一实心单元的材料体积，P为密度ρ_i的指数，且P∈Z,P＞1，K⁰为任一实心单元的单元刚度矩阵，且B是任一实心单元的应变矩阵；该柔顺机构第一层拓扑优化模型的体积约束如下：

$V\left(\mathrm{\ρ}\right)=\underset{i=1}{\overset{N_1}{\mathrm{\Σ}}}{V}_e{\mathrm{\ρ}}_i\≤{\mathrm{\θ}}_1^{*}{V}_0;$

式中，V₀是柔顺机构初始拓扑优化模型可利用材料的体积，ρ是由ρ_i,i＝1,2,…,N所构成的列向量；

根据上述目标和约束条件，得到柔顺机构第一层拓扑优化模型如下：

$\left\{\begin{array}{c}\begin{array}{cc}\mathrm{Min}:& f\left(\mathrm{\&rho;}\right)=-\frac{E_{\mathrm{ms}}}{E_s}\end{array}\\ \begin{array}{cc}s.t.& \mathrm{KU}=F_{\mathrm{in}}\end{array}\\ {\mathrm{KU}}_d=F_d\\ \underset{i=1}{\overset{N_1}{\mathrm{\&Sigma;}}}{V}_e{\mathrm{\&rho;}}_i\&le;{\mathrm{\&theta;}}_1^{*}{V}_0\\ 0<{\mathrm{\&rho;}}_{\min }\&le;{\mathrm{\&rho;}}_i\&le;{\mathrm{\&rho;}}_{\max }=1,i=\mathrm{1,2},...,N_1\end{array}\right..$

优选地，步骤2-1)中，优化目标的灵敏度为：

$\frac{\&PartialD;f}{{\&PartialD;\mathrm{\&rho;}}_i}=\frac{E_s\left(P{\left({\mathrm{\&rho;}}_i\right)}^{P-1}{\left(u_i\right)}^T{K}^0{u}_{\mathrm{di}}\right)-E_{\mathrm{ms}}\left(P{\left({\mathrm{\&rho;}}_i\right)}^{P-1}{\left(u_i\right)}^T{K}^0{u}_i\right)}{{\left(E_s\right)}^2},i=1,...,N_1;$

体积约束的灵敏度为： $\frac{\&PartialD;V}{{\&PartialD;\mathrm{\&rho;}}_i}=V_e,i=1,...,N_1.$

优选地，步骤2-2)中，过滤后的灵敏度如下：

$\frac{\&PartialD;\hat{f}}{{\&PartialD;\mathrm{\&rho;}}_i}=\frac{1}{{\mathrm{\&rho;}}_i\underset{j\&Element;N_e}{\mathrm{\&Sigma;}}{\hat{H}}_j}\underset{j\&Element;N_e}{\mathrm{\&Sigma;}}{\hat{H}}_j{\mathrm{\&rho;}}_j\frac{\&PartialD;f}{{\&PartialD;\mathrm{\&rho;}}_j},i=1,...,N_1;$

式中，N_e为单元i的邻域,该邻域内的各单元中心到单元i的中心的距离小于或等于过滤半径r，是卷积因子，如下：

${\hat{H}}_j=r-\left|\right|x_j-x_i\left|\right|;$

式中，x_i是单元i的中心坐标。

优选地，步骤2-3)中，最佳准则如下：

${\mathrm{\&rho;}}_i^{t+1}=\left\{\begin{array}{cc}{\mathrm{\&rho;}}_{\max }& ,i\&Element;{\mathrm{\&Omega;}}_{L,2}\\ \left\{\begin{array}{ccc}{\mathrm{\&rho;}}_i^t{\left(M_i^t\right)}^{\mathrm{\&eta;}}& \mathrm{if}& \max \{(1-\mathrm{\&zeta;}){\mathrm{\&rho;}}_i^t,{\mathrm{\&rho;}}_{\min }\}<{\mathrm{\&rho;}}_i^t{\left(M_i^t\right)}^{\mathrm{\&eta;}}<\min \{(1+\mathrm{\&zeta;}){\mathrm{\&rho;}}_i^t,{\mathrm{\&rho;}}_{\max }\}\\ \min \{(1+\mathrm{\&zeta;}){\mathrm{\&rho;}}_i^t,{\mathrm{\&rho;}}_{\max }\}& \mathrm{if}& {\mathrm{\&rho;}}_i^t{\left(M_i^t\right)}^{\mathrm{\&eta;}}\&GreaterEqual;\min \{(1+\mathrm{\&zeta;}){\mathrm{\&rho;}}_i^t,{\mathrm{\&rho;}}_{\max }\}\\ \max \{(1-\mathrm{\&zeta;}){\mathrm{\&rho;}}_i^t,{\mathrm{\&rho;}}_{\min }\}& \mathrm{if}& {\mathrm{\&rho;}}_i^t{\left(M_i^t\right)}^{\mathrm{\&eta;}}<\max \{(1-\mathrm{\&zeta;}){\mathrm{\&rho;}}_i^t,{\mathrm{\&rho;}}_{\min }\}\end{array}\right.& ,i\&Element;{\mathrm{\&Omega;}}_L\backslash {\mathrm{\&Omega;}}_{L,2},L\&GreaterEqual;1\\ {\mathrm{\&rho;}}_{\min }& ,i\&Element;{\mathrm{\&Omega;}}_1\backslash {\mathrm{\&Omega;}}_L\end{array}\right.;$

式中，Ω_L是步骤2中优化模型的可利用材料域，Ω_L,2是步骤2的可利用材料域中密度为ρ_max且恒定不变的单元集合，且当L＝1时，Ω_L,2为空集，Ω_L\Ω_L,2是步骤2的可利用材料域中密度可变的单元的集合，Ω₁\Ω_L是步骤2的优化模型中所有已删除单元的集合，其中各单元密度为ρ_min且恒定不变，当L＝1时，Ω₁\Ω_L为空集，t为迭代数,η为松弛因子,0＜η＜1,ζ为一个较小的移动极限,为一个非负数，即：

式中，可由Kuhn–Tucker必要条件导出，如下：

$Q_i^t={(-\frac{\&PartialD;f}{{\&PartialD;\mathrm{\&rho;}}_i}/\mathrm{\&lambda;}\frac{\&PartialD;V}{{\&PartialD;\mathrm{\&rho;}}_i})}^t;$

式中，λ为拉格朗日乘子，在设计变量更新过程中采用二分法得到，以使更新后的密度满足体积约束，即：

$\begin{array}{c}{\mathrm{\&theta;}}_L^{*}{V}_0-\{\underset{\mathrm{icas}1}{\mathrm{\&Sigma;}}{V}_e{\mathrm{\&rho;}}_i^t{\left(M_i^t\right)}^{\mathrm{\&eta;}}+\underset{\mathrm{icas}2}{\mathrm{\&Sigma;}}{V}_e\min \{(1+\mathrm{\&zeta;}){\mathrm{\&rho;}}_i^t,{\mathrm{\&rho;}}_{\max }\}+\underset{\mathrm{icas}3}{\mathrm{\&Sigma;}}{V}_e\max \{(1-\mathrm{\&zeta;}){\mathrm{\&rho;}}_i^t,{\mathrm{\&rho;}}_{\min }\}\\ +\underset{i\&Element;{\mathrm{\&Omega;}}_1\backslash {\mathrm{\&Omega;}}_L}{\mathrm{\&Sigma;}}{V}_e{\mathrm{\&rho;}}_{\min }+\underset{i\&Element;{\mathrm{\&Omega;}}_{L,2}}{\mathrm{\&Sigma;}}{V}_e{\mathrm{\&rho;}}_{\max }\}=0,L\&GreaterEqual;1\end{array};$

式中，是步骤2的体积比，和分别是集合Ω_L\Ω_L,2中所更新的三类设计变量的和。

优选地，步骤4)中，首先修改柔顺机构拓扑优化的层数L，使L＝L+1；采用渐进方法，按单元密度从小到大顺序，删除一定数量低密度单元，使其密度值均为ρ_min，直至设计域中所保留的单元数N_L为：

$N_L=\max \left(\right[\frac{N_1(V_{L-1}^{*}-{\mathrm{\&rho;}}_{\min })}{{\mathrm{\&rho;}}_{\max }-{\mathrm{\&rho;}}_{\min }}+1],N_s),L\&GreaterEqual;2;$

式中，max代表最大值，符号[]表示取整运算，N_s是步骤2的优化结果中密度为ρ_max的单元数；所有N_L个保留单元的集合构成下一层拓扑优化模型的可利用材料域Ω_L，所述下一层即第L层，L≥2，其中N_s个密度为ρ_max的单元构成下一层优化模型的可利用材料域中恒定不变的集合Ω_L,2；

接着，修改体积比：采用体积约束延拓方法，对下一层优化模型，体积比θ_L ^*取为：

${\mathrm{\&theta;}}_L^{*}=\max ((\mathrm{\&alpha;}+\mathrm{\&Delta;\&alpha;}(L-2)){\mathrm{\&theta;}}_{L-1}^{*},{\mathrm{\&theta;}}^{*}),L\&GreaterEqual;2;$

式中，max代表最大值，α是体积约束延拓因子初值，且0＜α＜1，Δα是体积约束延拓因子增量，且0＜Δα＜1，0＜α+Δα(L-2)<1；

最后，得到所述下一层柔顺机构拓扑优化模型为：

$\left\{\begin{array}{c}\begin{array}{cc}\mathrm{Min}:& f\left(\mathrm{\&rho;}\right)=-\frac{E_{\mathrm{ms}}}{E_s}\end{array}\\ \begin{array}{cc}s.t.& \mathrm{KU}=F_{\mathrm{in}}\end{array}\\ {\mathrm{KU}}_d=F_d\\ \underset{i=1}{\overset{N_1}{\mathrm{\&Sigma;}}}{V}_e{\mathrm{\&rho;}}_i\&le;{\mathrm{\&theta;}}_L^{*}{V}_0\\ 0<{\mathrm{\&rho;}}_{\min }\&le;{\mathrm{\&rho;}}_i\&le;{\mathrm{\&rho;}}_{\max }=1,i={\mathrm{\&Omega;}}_L\backslash {\mathrm{\&Omega;}}_{L,2}\\ {\mathrm{\&rho;}}_i={\mathrm{\&rho;}}_{\max },i\&Element;{\mathrm{\&Omega;}}_{L,2}\\ {\mathrm{\&rho;}}_i={\mathrm{\&rho;}}_{\min },i\&Element;{\mathrm{\&Omega;}}_1\backslash {\mathrm{\&Omega;}}_L\end{array}\right.;$

式中，Ω_L\Ω_L,2是第L层优化模型的可利用材料域中密度可变的单元集合，Ω₁\Ω_L是第L层优化模型中所有已删除的单元集合，L≥2。

本发明相对于现有技术，具有以下有益效果：

1.本发明可有效克服基于SIMP方法松弛设计变量时所产生的中间单元问题，提取出理想的柔顺机构0-1拓扑图。

2.本发明完全不同于传统渐进拓扑优化方法，它具有连续设计变量，灵敏度连续性好，而且概念简单，方法易于执行，不易导致棋盘格现象。

3.本发明将体积约束延拓方法和渐进方法相结合，能最大限度地保留所有可能单元，避免了添加材料的灵敏度评估困扰。

4.本发明克服了传统方法的参数敏感性和网格依赖性的缺陷，具有较好的收敛稳定性。

附图说明

图1为本发明柔顺机构0-1拓扑图提取方法流程图；

图2为实施例柔顺力-位移反向机构设计域和边界条件示意图；

图3为α＝0.8,Δα＝0.05以及具有120×60网格离散方式时，实施例柔顺力-位移反向机构0-1拓扑图提取结果；

图4为α＝0.8,Δα＝0.05以及具有180×90网格离散方式时，实施例柔顺力-位移反向机构0-1拓扑图提取结果；

图5为α＝0.8,Δα＝0.05以及具有120×60网格离散方式时，实施例柔顺力-位移反向机构0-1拓扑图提取结果；

图6为α＝0.8,Δα＝0.05以及具有180×90网格离散方式时，实施例柔顺力-位移反向机构0-1拓扑图提取结果；

图7为α＝0.85,Δα＝0.025以及具有120×60网格离散方式时，实施例柔顺力-位移反向机构0-1拓扑图提取结果；

图8为α＝0.85,Δα＝0.025以及具有180×90网格离散方式时，实施例柔顺力-位移反向机构0-1拓扑图提取结果；

图9为α＝0.85,Δα＝0.025以及具有120×60网格离散方式时，实施例柔顺力-位移反向机构0-1拓扑图提取结果；

图10为α＝0.85,Δα＝0.025以及具有180×90网格离散方式时，实施例柔顺力-位移反向机构0-1拓扑图提取结果。

具体实施方式

下面结合附图和具体实施例对本发明的发明目的作进一步详细地描述，实施例不能在此一一赘述，但本发明的实施方式并不因此限定于以下实施例。除非特别说明，本发明采用的材料和加工方法为本技术领域常规材料和加工方法。

如图1所示，为本发明所提出的一种柔顺机构0-1拓扑图提取方法流程图。

本发明实施例是一种典型的柔顺力-位移反向机构，其设计域和边界条件如图2所示。其中，设计域大小为120×120，材料的弹性模量和泊松比分别为E＝1和v＝0.3，P_i和P_o分别为载荷输入点和位移输出点，输入载荷F_in＝1，同时，沿输出位移方向作用有一虚拟单位载荷F_d，输入和输出弹簧刚度分别为k_in＝1和k_out＝0.001。

柔顺机构的目标体积比为θ^*＝0.25，灵敏度过滤半径为r＝2.5。由于设计域和边界条件具有对称性，因此仅取设计域的下半区域进行优化。采用SIMP方法松弛密度变量，以柔顺机构的最大柔度和最大刚度为优化目标，得到柔顺机构第一层拓扑优化模型为：

$\left\{\begin{array}{c}\begin{array}{cc}\mathrm{Min}:& f\left(\mathrm{\&rho;}\right)=-\frac{E_{\mathrm{ms}}}{E_s}\end{array}\\ \begin{array}{cc}s.t.& \mathrm{KU}=F_{\mathrm{in}}\end{array}\\ {\mathrm{KU}}_d=F_d\\ \underset{i=1}{\overset{N_1}{\mathrm{\&Sigma;}}}{V}_e{\mathrm{\&rho;}}_i\&le;{\mathrm{\&theta;}}_1^{*}{V}_0\\ 0<{\mathrm{\&rho;}}_{\min }\&le;{\mathrm{\&rho;}}_i\&le;{\mathrm{\&rho;}}_{\max }=1,i=\mathrm{1,2},...,N_1\end{array}\right.;$

式中，ρ_i是单元i的密度，ρ_min是单元密度下限，ρ_max是单元密度上限，N是设计域内的单元数，V_e是任一实心单元的材料体积，V₀是柔顺机构初始拓扑优化模型可利用材料的体积，ρ_i是单元i的密度，U是载荷F_in作用下的节点位移矢量，U_d是载荷F_d作用下的节点位移矢量，N₁是自然数，表明初始阶段可利用材料域单元数，E_ms和E_s分别是***的互应变能和应变能，公式为：

E_s＝∫_Ωε(u)^TDε(u)dΩ＝U^TKU；

$E_{\mathrm{ms}}={\&Integral;}_{\mathrm{\&Omega;}}{\mathrm{\&epsiv;}\left(u_d\right)}^T\mathrm{D\&epsiv;}\left(u\right)\mathrm{d\&Omega;}=U_d^T\mathrm{KU};$

式中，D为弹性矩阵，ε(u)和u是设计域内任一点在载荷F作用下的应变和弹性变形，ε(u_d)和u_d是设计域内任一点在载荷F_d作用下的应变和弹性变形,K是***的整体刚度矩阵，公式为；

$\begin{array}{c}K=\underset{i=1}{\overset{N_1}{\mathrm{\&Sigma;}}}{K}_i\\ =\underset{i=1}{\overset{N_1}{\mathrm{\&Sigma;}}}{\&Integral;}_{V_e}{\mathrm{\&rho;}}_i^P{B}^T\mathrm{DBdV}\\ =\underset{i=1}{\overset{N_1}{\mathrm{\&Sigma;}}}{\mathrm{\&rho;}}_i^P{K}^0\end{array}$

式中，K_i是单元i的刚度矩阵，K⁰为任一实心单元的单元刚度矩阵，且B是任一实心单元的应变矩阵,P为密度ρ_i的指数，且P∈Z,P＞1，本实施例中，为了加快收敛，取P＝3；

为了获得优化模型的拓扑优化最佳结果，按以下公式对其进行有限元分析，求出节点位移向量U和U_d：

KU＝F；

KU_d＝F_d。

进一步得出优化目标和体积约束的灵敏度，公式如下:

$\frac{\&PartialD;f}{{\&PartialD;\mathrm{\&rho;}}_i}=\frac{E_s\left(P{\left({\mathrm{\&rho;}}_i\right)}^{P-1}{\left(u_i\right)}^T{K}^0{u}_{\mathrm{di}}\right)-E_{\mathrm{ms}}\left(P{\left({\mathrm{\&rho;}}_i\right)}^{P-1}{\left(u_i\right)}^T{K}^0{u}_i\right)}{{\left(E_s\right)}^2},i=1,...,N_1;$

$\frac{\&PartialD;V}{{\&PartialD;\mathrm{\&rho;}}_i}=V_e,i=1,...,N_1.$

再对优化目标的灵敏度进行过滤，公式如下：

$\frac{\&PartialD;\hat{f}}{{\&PartialD;\mathrm{\&rho;}}_i}=\frac{1}{{\mathrm{\&rho;}}_i\underset{j\&Element;N_e}{\mathrm{\&Sigma;}}{\hat{H}}_j}\underset{j\&Element;N_e}{\mathrm{\&Sigma;}}{\hat{H}}_j{\mathrm{\&rho;}}_j\frac{\&PartialD;f}{{\&PartialD;\mathrm{\&rho;}}_j},i=1,...,N_1;$

式中，N_e为单元i的邻域,该邻域内的各单元中心到单元i的中心的距离小于或等于过滤半径r，是卷积因子，即：式中，x_i是单元i的中心坐标。

之后，再利用最佳准则更新设计变量,公式如下:

${\mathrm{\&rho;}}_i^{t+1}=\left\{\begin{array}{cc}{\mathrm{\&rho;}}_{\max }& ,i\&Element;{\mathrm{\&Omega;}}_{L,2}\\ \left\{\begin{array}{ccc}{\mathrm{\&rho;}}_i^t{\left(M_i^t\right)}^{\mathrm{\&eta;}}& \mathrm{if}& \max \{(1-\mathrm{\&zeta;}){\mathrm{\&rho;}}_i^t,{\mathrm{\&rho;}}_{\min }\}<{\mathrm{\&rho;}}_i^t{\left(M_i^t\right)}^{\mathrm{\&eta;}}<\min \{(1+\mathrm{\&zeta;}){\mathrm{\&rho;}}_i^t,{\mathrm{\&rho;}}_{\max }\}\\ \min \{(1+\mathrm{\&zeta;}){\mathrm{\&rho;}}_i^t,{\mathrm{\&rho;}}_{\max }\}& \mathrm{if}& {\mathrm{\&rho;}}_i^t{\left(M_i^t\right)}^{\mathrm{\&eta;}}\&GreaterEqual;\min \{(1+\mathrm{\&zeta;}){\mathrm{\&rho;}}_i^t,{\mathrm{\&rho;}}_{\max }\}\\ \max \{(1-\mathrm{\&zeta;}){\mathrm{\&rho;}}_i^t,{\mathrm{\&rho;}}_{\min }\}& \mathrm{if}& {\mathrm{\&rho;}}_i^t{\left(M_i^t\right)}^{\mathrm{\&eta;}}<\max \{(1-\mathrm{\&zeta;}){\mathrm{\&rho;}}_i^t,{\mathrm{\&rho;}}_{\min }\}\end{array}\right.& ,i\&Element;{\mathrm{\&Omega;}}_L\backslash {\mathrm{\&Omega;}}_{L,2},L\&GreaterEqual;1\\ {\mathrm{\&rho;}}_{\min }& ,i\&Element;{\mathrm{\&Omega;}}_1\backslash {\mathrm{\&Omega;}}_L\end{array}\right.;$

式中，Ω_L是步骤2中优化模型的可利用材料域，Ω_L,2是步骤2的可利用材料域中密度为ρ_max且恒定不变的单元集合，当L＝1时，Ω_L,2为空集，Ω_L\Ω_L,2是步骤2的可利用材料域中密度可变的单元集合，Ω₁\Ω_L是步骤2的优化模型中所有已删除单元的集合，其中各单元密度为ρ_min且恒定不变，当L＝1时，Ω₁\Ω_L为空集，t为迭代数,η为松弛因子,0＜η＜1,ζ为一个较小的移动极限，在本实施例中，为了保证密度更新过程的稳定性，η和ζ分别取值为0.3和0.05，为一个非负数，即：

式中，可由Kuhn–Tucker必要条件导出，公式是：式中，λ为拉格朗日乘子，在设计变量更新过程中采用二分法得到，以使更新后的密度满足体积约束，即：

$\begin{array}{c}{\mathrm{\&theta;}}_L^{*}{V}_0-\{\underset{\mathrm{icas}1}{\mathrm{\&Sigma;}}{V}_e{\mathrm{\&rho;}}_i^t{\left(M_i^t\right)}^{\mathrm{\&eta;}}+\underset{\mathrm{icas}2}{\mathrm{\&Sigma;}}{V}_e\min \{(1+\mathrm{\&zeta;}){\mathrm{\&rho;}}_i^t,{\mathrm{\&rho;}}_{\max }\}+\underset{\mathrm{icas}3}{\mathrm{\&Sigma;}}{V}_e\max \{(1-\mathrm{\&zeta;}){\mathrm{\&rho;}}_i^t,{\mathrm{\&rho;}}_{\min }\}\\ +\underset{i\&Element;{\mathrm{\&Omega;}}_1\backslash {\mathrm{\&Omega;}}_L}{\mathrm{\&Sigma;}}{V}_e{\mathrm{\&rho;}}_{\min }+\underset{i\&Element;{\mathrm{\&Omega;}}_{L,2}}{\mathrm{\&Sigma;}}{V}_e{\mathrm{\&rho;}}_{\max }\}=0,L\&GreaterEqual;1\end{array};$

式中，是步骤2的体积比，和分别是集合Ω_L\Ω_L,2中所更新的三类设计变量的和。

利用最佳准则更新设计变量的具体过程是：

首先，给定λ的最小值和最大值，分别是λ₁＝0和λ₂＝100000；

然后，计算它们的中点值：λ＝0.5(λ₁+λ₂)，并代入密度变量更新公式，得到更新后的单元密度；

之后，计算体积约束余量再根据▽V值修改λ₁或λ₂的值：如果▽V＞0，则使λ₁＝λ，否则如果▽V＜0，则使λ₂＝λ；

重复λ₁和λ₂的中点值的计算、单元密度更新以及λ₁或λ₂的修改过程，直至λ₂与λ₁的差小于一个阈值为止。

重复有限元分析、灵敏度计算和过滤以及设计变量更新过程，直到优化迭代收敛。

判断最大迭代数或前后两次迭代的单元密度变化最大值小于一个阈值的条件是否满足，若满足则终止循环并输出结果；否则修改柔顺机构拓扑优化的层数L，使L＝L+1；然后，采用渐进方法，按单元密度从小到大顺序，删除一定数量低密度单元，使其密度值均为ρ_min，直至设计域中所保留的单元数N_L为：

$N_L=\max \left(\right[\frac{N_1(V_{L-1}^{*}-{\mathrm{\&rho;}}_{\min })}{{\mathrm{\&rho;}}_{\max }-{\mathrm{\&rho;}}_{\min }}+1],N_s),L\&GreaterEqual;2.$

式中，max代表最大值，符号[]表示取整运算，N_s是步骤2的优化结果中密度为ρ_max的单元数；所有N_L个保留单元的集合构成下一层(即第L,(L≥2)层)拓扑优化模型的可利用材料域Ω_L，其中N_s个密度为ρ_max的单元构成下一层优化模型的可利用材料域中恒定不变的集合Ω_L,2；

接着，修改体积比：采用体积约束延拓方法，对下一层优化模型，体积比取为：

${\mathrm{\&theta;}}_L^{*}=\max ((\mathrm{\&alpha;}+\mathrm{\&Delta;\&alpha;}(L-2)){\mathrm{\&theta;}}_{L-1}^{*},{\mathrm{\&theta;}}^{*}),L\&GreaterEqual;2.$

式中，max代表最大值，α是体积约束延拓因子初值，0＜α＜1，Δα是体积约束延拓因子增量，0＜Δα＜1，且α+Δα(L-2)≤1。为了验证不同的参数设置对本发明所提出的方法的性能的影响，在本实施例中，对参数α和Δα选取了四组数值，分别是：

(1)α＝0.8,Δα＝0.05；

(2)α＝0.8,Δα＝0.05；

(3)α＝0.85,Δα＝0.025；

(4)α＝0.85,Δα＝0.025；

并且，给定α+Δα(L-2)的上限为0.99，以使逐步逼近θ^*。

同时，对每一组参数值情形，分别采用了两种不同的网格离散方式，即：120×60和180×90。

最后，得到下一层(即第L,(L≥2)层)柔顺机构拓扑优化模型为：

$\left\{\begin{array}{c}\begin{array}{cc}\mathrm{Min}:& f=-\frac{E_{\mathrm{ms}}}{E_s}\end{array}\\ \begin{array}{cc}s.t.& \mathrm{KU}=F_{\mathrm{in}}\end{array}\\ {\mathrm{KU}}_d=F_d\\ \underset{i=1}{\overset{N_1}{\mathrm{\&Sigma;}}}{V}_e{\mathrm{\&rho;}}_i\&le;{\mathrm{\&theta;}}_L^{*}{V}_0\\ 0<{\mathrm{\&rho;}}_{\min }\&le;{\mathrm{\&rho;}}_i\&le;{\mathrm{\&rho;}}_{\max }=1,i={\mathrm{\&Omega;}}_L\backslash {\mathrm{\&Omega;}}_{L,2}\\ {\mathrm{\&rho;}}_i={\mathrm{\&rho;}}_{\max },i\&Element;{\mathrm{\&Omega;}}_{L,2}\\ {\mathrm{\&rho;}}_i={\mathrm{\&rho;}}_{\min },i\&Element;{\mathrm{\&Omega;}}_1\backslash {\mathrm{\&Omega;}}_L\end{array}\right..$

式中，Ω_L\Ω_L,2是第L,(L≥2)层优化模型的可利用材料域中密度可变的单元集合，Ω₁\Ω_L是第L,(L≥2)层优化模型中所有已删除的单元集合。

再对该拓扑优化模型进行求解，如此不断缩小可利用材料域，直至提取出满足目标体积比θ^*的理想的柔顺机构0-1拓扑图。

图3、图5、图7和图9分别是在120×60网格离散方式下，对应于参数α和Δα的四组取值情形时，本发明采用的一种典型的柔顺力-位移反向机构最终0-1拓扑图提取结果。

图4、图6、图8和图10分别是在180×90网格离散方式下，对应于参数α和Δα的四组取值情形时，本发明采用的一种典型的柔顺力-位移反向机构最终0-1拓扑图提取结果。

图3至图10表明，本发明所提出的柔顺机构0-1拓扑图提取方法能获得理想的黑白拓扑结果，同时能避免传统方法的参数敏感性和网格依赖性的缺陷，具有较好的收敛稳定性。

上述实施例仅为本发明的较佳实施例，并非用来限定本发明的实施范围。即凡依本发明内容所作的均等变化与修饰，都为本发明权利要求所要求保护的范围所涵盖。

Claims (5)

1.一种柔顺机构0-1拓扑图提取方法，其特征在于包括如下步骤：

1)：建立柔顺机构第一层拓扑优化模型；

2)：对柔顺机构拓扑优化模型进行优化，步骤如下:

2-1)：对柔顺机构拓扑优化模型进行有限元分析，并得出体积约束和优化目标的灵敏度；

2-2)：对优化目标的灵敏度进行过滤，以消除拓扑图中的棋盘格；

2-3)：基于最佳准则更新设计变量；

2-4)：重复步骤2-1)至2-3)，直到优化迭代收敛；

3)：判断柔顺机构拓扑优化最大迭代数或前后两次迭代的单元密度变化最大值小于一个阈值的条件是否满足，若满足则终止循环并输出结果，否则继续执行以下步骤；

4)：建立下一层柔顺机构拓扑优化模型,

5)：重复步骤2-1)至4)，如此不断缩小拓扑优化模型的可利用材料域，直至提取出满足目标体积V^*的柔顺机构0-1拓扑图；

步骤1)中，以Ω₁为设计域，其为柔顺机构初始拓扑优化模型可利用材料域；P_i和P_o分别为柔顺机构载荷输入点和位移输出点；F_in和F_d分别为柔顺机构输入载荷和沿输出位移方向的虚拟单位载荷；k_in和k_out分别为柔顺机构输入和输出弹簧刚度；该柔顺机构的目标体积比为θ^*，初始阶段的拓扑优化层数为L＝1，所述初始阶段的拓扑优化体积比为将设计域离散为N_x行N_y列；柔顺机构的应变能和互应变能如下：

E_s＝∫_Ωε(u)^TDε(u)dΩ＝U^TKU；

$E_{\mathrm{ms}}={\&Integral;}_{\mathrm{\&Omega;}}\mathrm{\&epsiv;}{\left(u_d\right)}^T\mathrm{D\&epsiv;}\left(u\right)\mathrm{d\&Omega;}=U_d^T\mathrm{KU};$

式中，E_s是***的应变能，应变能越小则表明***的刚度越大；E_ms是***的互应变能，互应变能越大则表明***的柔性越大；D为弹性矩阵，K是***整体刚度矩阵，U是F_in作用下的节点位移向量，U_d是F_d作用下的节点位移向量，ε(u)和u是设计域内任一点在载荷F作用下的应变和弹性变形，ε(u_d)和u_d是设计域内任一点在载荷F_d作用下的应变和弹性变形；

为使柔顺机构既有足够大的刚度又有足够大的柔性，通过多目标优化而得到柔顺机构的应变能和互应变能的关系如下：

Min:式中，Min代表最小值；

首先，基于相对密度法松弛设计变量，使柔顺机构单元密度在0-1范围内取值，即：

0＜ρ_min≤ρ_i≤ρ_max＝1,i＝1,2,…,N₁；

式中，ρ_i是单元i的密度，ρ_min是单元密度下限，ρ_max是单元密度上限，N₁是自然数，表明初始阶段可利用材料域的单元数，且N₁＝N_xN_y；

柔顺机构拓扑优化模型的整体刚度矩阵是：

$\begin{array}{c}K=\underset{i=1}{\overset{N_1}{\mathrm{\&Sigma;}}}{L}_i\\ =\underset{i=1}{\overset{N_1}{\mathrm{\&Sigma;}}}{\&Integral;}_{V_e}{\mathrm{\&rho;}}_i^P{B}^T\mathrm{DBdV}\\ =\underset{i=1}{\overset{N_1}{\mathrm{\&Sigma;}}}{\mathrm{\&rho;}}_i^P{K}^0\end{array};$

式中，K_i是单元i的刚度矩阵，V_e是任一实心单元的材料体积，P为密度ρ_i的指数，且P∈Z,P＞1，K⁰为任一实心单元的单元刚度矩阵，且B是任一实心单元的应变矩阵；该柔顺机构第一层拓扑优化模型的体积约束如下：

$V\left(\mathrm{\&rho;}\right)=\underset{i=1}{\overset{N_1}{\mathrm{\&Sigma;}}}{V}_e{\mathrm{\&rho;}}_i\&le;{\mathrm{\&theta;}}_1^{*}{V}_0;$

式中，V₀是柔顺机构初始拓扑优化模型可利用材料的体积，ρ是由ρ_i,i＝1,2,…,N所构成的列向量；

根据上述目标和约束条件，得到柔顺机构第一层拓扑优化模型如下：

$\left\{\begin{array}{cc}\mathrm{Min}:& f\left(\mathrm{\&rho;}\right)=-\frac{E_{\mathrm{ms}}}{E_s}\\ s.t.& \mathrm{KU}=F_{\mathrm{in}}\\ & {\mathrm{KU}}_d=F_d\\ & \underset{i=1}{\overset{N_1}{\mathrm{\&Sigma;}}}{V}_e{\mathrm{\&rho;}}_i\&le;{\mathrm{\&theta;}}_1^{*}{V}_0\\ & 0<{\mathrm{\&rho;}}_{\min }\&le;{\mathrm{\&rho;}}_i\&le;{\mathrm{\&rho;}}_{\max }=1,i=\mathrm{1,2},...,N_1\end{array}\right..$

2.根据权利要求1所述的柔顺机构0-1拓扑图提取方法，其特征在于：步骤2-1)中，优化目标的灵敏度为：

$\frac{\&PartialD;f}{{\&PartialD;\mathrm{\&rho;}}_i}=\frac{E_s\left(P{\left({\mathrm{\&rho;}}_i\right)}^{P-1}{\left(u_i\right)}^T{K}^0{u}_{\mathrm{di}}\right)-E_{\mathrm{ms}}\left(P{\left({\mathrm{\&rho;}}_i\right)}^{P-1}{\left(u_i\right)}^T{K}^0{u}_i\right)}{{\left(E_s\right)}^2},i=1,...,N_1;$

体积约束的灵敏度为：i＝1,...,N₁。

3.根据权利要求2所述的柔顺机构0-1拓扑图提取方法，其特征在于：步骤2-2)中，过滤后的灵敏度如下：

$\frac{\&PartialD;\hat{f}}{{\&PartialD;\mathrm{\&rho;}}_i}=\frac{1}{{\mathrm{\&rho;}}_i\underset{j\&Element;N_e}{\mathrm{\&Sigma;}}{\hat{H}}_j}\underset{j\&Element;N_e}{\mathrm{\&Sigma;}}{\hat{H}}_j{\mathrm{\&rho;}}_j\frac{\&PartialD;f}{{\&PartialD;\mathrm{\&rho;}}_j},i=1,...,N_1;$

式中，N_e为单元i的邻域,该邻域内的各单元中心到单元i的中心的距离小于或等于过滤半径r，是卷积因子，如下：

${\hat{H}}_j=r-\left|\right|x_j-x_i\left|\right|;$

式中，x_i是单元i的中心坐标。

4.根据权利要求3所述的柔顺机构0-1拓扑图提取方法，其特征在于：步骤2-3)中，最佳准则如下：

${\mathrm{\&rho;}}_i^{t+1}=\left\{\begin{array}{cc}{\mathrm{\&rho;}}_{\max }& ,i\&Element;{\mathrm{\&Omega;}}_{L,2}\\ \left\{\begin{array}{ccc}{\mathrm{\&rho;}}_i^t{\left(M_i^t\right)}^{\mathrm{\&eta;}}& \mathrm{i\; f}& \max \{(1-\mathrm{\&zeta;}){\mathrm{\&rho;}}_i^t,{\mathrm{\&rho;}}_{\min }\}<{\mathrm{\&rho;}}_i^t{\left(M_i^t\right)}^{\mathrm{\&eta;}}<\min \{(1+\mathrm{\&zeta;}){\mathrm{\&rho;}}_i^t,{\mathrm{\&rho;}}_{\max }\}\\ \min \{(1+\mathrm{\&zeta;}){\mathrm{\&rho;}}_i^t,{\mathrm{\&rho;}}_{\mathrm{ma\; x}}\}& \mathrm{if}& {\mathrm{\&rho;}}_i^t{\left(M_i^t\right)}^{\mathrm{\&eta;}}\text{\&GreaterEqual;min{}(1+\mathrm{\&zeta;}){\mathrm{\&rho;}}_i^t,{\mathrm{\&rho;}}_{\max }\}\\ \max \{(1-\mathrm{\&zeta;}){\mathrm{\&rho;}}_i^t,{\mathrm{\&rho;}}_{\min }& \mathrm{if}& {\mathrm{\&rho;}}_i^t{\left(M_i^t\right)}^{\mathrm{\&eta;}}\&le;\max \{(1-\mathrm{\&zeta;}){\mathrm{\&rho;}}_i^t,{\mathrm{\&rho;}}_{\min }\}\end{array}\right.& ,i\&Element;{\mathrm{\&Omega;}}_L\backslash {\mathrm{\&Omega;}}_{L,2},L\&GreaterEqual;1\\ {\mathrm{\&rho;}}_{\min }& ,i\&Element;{\mathrm{\&Omega;}}_1\backslash {\mathrm{\&Omega;}}_L\end{array};\right.$

式中，Ω_L是步骤2中优化模型的可利用材料域，Ω_L,2是步骤2的可利用材料域中密度为ρ_max且恒定不变的单元集合，当L＝1时，Ω_L,2为空集，Ω_L\Ω_L,2是步骤2的可利用材料域中密度可变的单元集合，Ω₁\Ω_L是步骤2的优化模型中所有已删除单元的集合，其中各单元密度为ρ_min且恒定不变，当L＝1时，Ω₁\Ω_L为空集，t为迭代数,η为松弛因子,0＜η＜1，ζ为一个较小的移动极限,为一个非负数，即：

式中，由Kuhn–Tucker必要条件导出，如下：

$Q_i^t={(-\frac{\&PartialD;f}{{\&PartialD;\mathrm{\&rho;}}_i}/\mathrm{\&lambda;}\frac{\&PartialD;V}{{\&PartialD;\mathrm{\&rho;}}_i})}^t;$

式中，λ为拉格朗日乘子，在设计变量更新过程中采用二分法得到，以使更新后的密度满足体积约束，即：

$\begin{array}{c}{\mathrm{\&theta;}}_L^{*}{V}_0-\{\underset{\mathrm{icas}1}{\mathrm{\&Sigma;}}{V}_e{\mathrm{\&rho;}}_i^t{\left(M_i^t\right)}^{\mathrm{\&eta;}}+\underset{\mathrm{icas}2}{\mathrm{\&Sigma;}}{V}_e\min \{(1+\mathrm{\&zeta;}){\mathrm{\&rho;}}_i^t,{\mathrm{\&rho;}}_{\max }\}+\underset{\mathrm{icas}3}{\mathrm{\&Sigma;}}{V}_e\max \{(1-\mathrm{\&zeta;}){\mathrm{\&rho;}}_i^t,{\mathrm{\&rho;}}_{\min }\}\\ +\underset{i\&Element;{\mathrm{\&Omega;}}_1\backslash {\mathrm{\&Omega;}}_L}{\mathrm{\&Sigma;}}{V}_e{\mathrm{\&rho;}}_{\min }+\underset{i\&Element;{\mathrm{\&Omega;}}_{L,2}}{\mathrm{\&Sigma;}}{V}_e{\mathrm{\&rho;}}_{\max }\}=0,L\&GreaterEqual;1\end{array};$

式中，是步骤2的体积比，和分别是集合Ω_L\Ω_L,2中所更新的三类设计变量的和。

5.根据权利要求4所述的柔顺机构0-1拓扑图提取方法，其特征在于：步骤4)中，首先修改柔顺机构拓扑优化的层数L，使L＝L+1；采用渐进方法，按单元密度从小到大顺序，删除一定数量低密度单元，使其密度值均为ρ_min，直至设计域中所保留的单元数N_L为：

$N_L=\max \left(\right[\frac{N_1(V_{L-1}^{*}-{\mathrm{\&rho;}}_{\min })}{{\mathrm{\&rho;}}_{\max }-{\mathrm{\&rho;}}_{\min }}+1],N_s),L\&GreaterEqual;2;$

式中，max代表最大值，符号[]表示取整运算，N_s是步骤2的优化结果中密度为ρ_max的单元数；所有N_L个保留单元的集合构成下一层拓扑优化模型的可利用材料域Ω_L，所述下一层即第L层，L≥2，其中N_s个密度为ρ_max的单元构成下一层优化模型的可利用材料域中恒定不变的集合Ω_L,2；

接着，修改体积比：采用体积约束延拓方法，对下一层优化模型，体积比取为：

${\mathrm{\&theta;}}_L^{*}=\max ((\mathrm{\&alpha;}+\mathrm{\&Delta;\&alpha;}(L-2)){\mathrm{\&theta;}}_{L-1}^{*},{\mathrm{\&theta;}}^{*}),L\&GreaterEqual;2;$

式中，max代表最大值，α是体积约束延拓因子初值，且0＜α＜1，Δα是体积约束延拓因子增量，且0＜Δα＜1，0＜α+Δα(L-2)<1；

最后，得到所述下一层柔顺机构拓扑优化模型为：

$\left\{\begin{array}{cc}\mathrm{Min}:& f\left(\mathrm{\&rho;}\right)=-\frac{E_{\mathrm{ms}}}{E_s}\\ s.t.& \mathrm{KU}=F_{\mathrm{in}}\\ & {\mathrm{KU}}_d=F_d\\ & \underset{i=1}{\overset{N_1}{\mathrm{\&Sigma;}}}{V}_e{\mathrm{\&rho;}}_i\&le;{\mathrm{\&theta;}}_L^{*}{V}_0\\ & <{\mathrm{\&rho;}}_{\min }\&le;{\mathrm{\&rho;}}_i\&le;{\mathrm{\&rho;}}_{\max }=1,i\&Element;{\mathrm{\&Omega;}}_L\backslash {\mathrm{\&Omega;}}_{L,2}\\ & {\mathrm{\&rho;}}_i={\mathrm{\&rho;}}_{\max },i\&Element;{\mathrm{\&Omega;}}_{L,2}\\ & {\mathrm{\&rho;}}_i={\mathrm{\&rho;}}_{\min },i\&Element;{\mathrm{\&Omega;}}_1\backslash {\mathrm{\&Omega;}}_L\end{array}\right.;$

式中，Ω_L\Ω_L,2是第L层优化模型的可利用材料域中密度可变的单元集合，Ω₁\Ω_L是第L层优化模型中所有已删除的单元集合，L≥2。