CN102043883A - 基于材料破损约束的连续体结构拓扑设计建模及优化设计方法 - Google Patents

Abstract

本发明公开了一种基于材料破损约束的连续体结构拓扑设计建模及优化设计方法，用于解决现有技术不同程度存在着的：(a)应力集中和奇异问题；(b)最大应力波动问题、分析和计算量大问题。该方法采用一种有效应力约束松弛式处理应力奇异现象。优化问题近似建模时，将结构应力的q₁范数度量函数作为罚函数，用结构应力的q₂范数度量函数约束、几个最潜在的主动单元应力约束和引入的体积约束替代所有单元应力约束，并结合变约束限措施控制局部应力水平。方法的优化过程分为二个优化阶段和一个阶段转化步，自动地实现设计空间的扩展和减缩。本发明方法利用二次规划法进行各阶段的求解，可获得纯的黑/白分布的最优结构拓扑，且具有良好的优化设计效率。

Description

基于材料破损约束的连续体结构拓扑设计建模及优化设计方法

技术领域

本发明涉及基于材料破损约束的连续体结构拓扑设计建模及优化设计方法，适用于复杂连续体结构系统的设计建模、分析和优化设计。属于结构优化设计技术领域。

技术背景

在现代结构设计理念中，要求结构尽可能轻量化的同时，也要求结构具备高度可靠性和足够的安全性。这两者之间常常出现矛盾。发展复杂结构和零部件的创新设计技术是解决这一矛盾、提高结构性能、降低结构重量和缩短结构设计周期的有效途径与发展新方向。大部分基于灵敏度分析的连续体结构拓扑优化方法所涉及的优化模型中，通常采用柔顺度最小和材料用量约束来表述，较少考虑材料破损约束的工程要求。虽然柔顺度最小设计能获得一个近似的满应力设计结构，然而理论结果表明：对相同的初始优化结构，最强硬的结构布局可能完全不同于最刚硬的结构布局。

近年来，研究领域提出了一些基于应力约束的连续体结构拓扑优化方法，可归纳如下：

1、基于应力的渐进结构优化法(Evolutionary Structural OptimizationMethod)是一种满应力设计准则法，其主要缺点是它仅仅在某特定的情况下获得满足应力约束的最小重量结构拓扑优化设计。

2、ICM(independent，continuous，mapping)方法引入对偶规划方法大大减少了设计变量的数目，提高了优化的效率。该方法采用结构应变能约束替代所有局部应力约束，使得优化后的结构会出现单元应力超限。

3、序列整型规划法能获得结构优化解，但求解的结构规模很小，不适合工程应用。

4、现有的基于均匀化法或实体各向同性材料惩罚法(简称SIMP法)的含应力要求的拓扑优化方法，采用每一种载荷工况的一个或二个总体应力聚集约束函数替代该载荷工况的所有单元应力约束函数来减小应力约束处理的工作量，并辅之以ε-松弛处理来解决应力奇异问题，这种方法放松了局部应力控制。采用该方法得到的最佳拓扑与局部应力约束设计获得的最佳拓扑有差异，且比较接近于柔顺度最小设计的结果。另有一些方法虽能解决应力集中问题，但最后仍需要解释得到的材料分布，很难准确获得原问题的最优解。另外ICM法和SIMP法的迭代优化过程中，被删除的大量材料单元仍以较小的刚度和质量存在于结构中，易引起结构总刚度矩阵病态以及结构出现虚假的局部模态，结构分析和灵敏度分析的规模始终保持为最大设计域的规模，其分析计算量大。考虑与质量分布关联的载荷作用时，被删除的大量材料单元的存在使得结构应力及其灵敏度分析有较大误差。

基于材料破损约束的连续体结构拓扑优化问题是未能彻底解决的一个难题，需要解决如下几个主要问题：(a)应力集中问题；(b)应力奇异问题；(c)大量应力约束归并后的几个聚集应力约束函数应能近似反映原局部应力约束的主要特征，以及最大应力波动问题；(d)结构分析和灵敏度分析的规模不变问题。

发明内容

本发明的目的在于克服现有技术之不足而提供一种基于材料破损约束的连续体结构拓扑设计建模及优化设计方法。本方法的优化过程分为二个优化阶段和一个阶段转化步，建立了二个优化阶段和阶段转化步的近似等效优化数学模型，并提出了相应优化阶段的求解算法。主要解决以结构体积为目标函数，结构单元应力小于或等于材料破坏应力为约束条件的连续体结构拓扑优化设计问题

本发明基于材料破损约束的连续体结构拓扑设计建模及优化设计方法，包括下述步骤：

第一步：有限元初始建模

根据需要设计的连续体结构的服役要求，确定一个拓扑结构的足够大设计区域，利用规则网格对该设计区域进行有限元划分，确定边界条件、载荷工况、结构拓扑变量初值和单元材料特性值，形成最大设计区域固定有限单元网格模型，形成设计用的初始结构并形成其保留单元数据文件；根据需要设计的连续体结构的用材，确定材料破损应力阀值；设置初始迭代步数l＝1；则设计用的初始结构为当前迭代步初始结构；

第二步：有限元模型更新

读取固定有限单元网格模型文件和当前迭代步初始结构保留单元数据文件，在结构周边和孔洞周围设置一层人工材料单元，形成当前迭代步初始有效结构并形成其有限单元网格模型，建立最大设计区域有限单元网格模型与当前迭代步初始有效结构有限单元网格模型的映射关系，得到映射关系数据文件；此时的当前迭代步初始有效结构也为上一迭代步优化获得的有效结构；

第三步：当前迭代步初始有效结构有限元分析

根据当前迭代步初始有效结构有限单元网格模型数据进行有限元分析，获得实际载荷作用下当前迭代步初始有效结构有限单元网格模型的结构位移、应力，生成基于伴随方程的应力灵敏度计算所需的虚载荷，计算虚载荷作用下当前迭代步初始有效结构位移、应力；读取所述映射关系数据文件，获得最大设计区域固定有限单元网格模型的结构位移、应力数据；

第四步：形成近似优化模型

在该迭代步中，当前迭代步初始有效结构随着结构拓扑变量变化将发生变化，其网格模型不变，其拓扑变量、位移和应力都将发生变化；则变化了拓扑变量的当前迭代步初始有效结构称为当前迭代步有效结构；

优化过程分为二个优化阶段和一个阶段转化步；刚开始优化迭代时，进行第一个优化阶段的近似优化模型建模；当式(1)成立时，进行阶段转化步的近似优化模型建模；后继的所有迭代步都进行第二个优化阶段的近似优化模型建模；

l≥l_ap，和|(V^(l-1)-V^(l-2))/V^(l-1)|＜ε_f    (1)

式中l是当前的迭代步数，l_ap是预先指定的限值，选取l_ap＝45；V^(l-1)和V^(l-2)分别是第(l-1)迭代步和第(l-2)迭代步优化获得的有效结构体积；ε_f是一预先指定的经验值，选取ε_f＝0.001；

4.1)第一个优化阶段和阶段转化步的近似优化模型建模

在结构拓扑优化中，始终存在应力奇异问题；这里根据实际载荷和虚载荷作用下当前迭代步初始有效结构位移、应力数据，根据式(2)建立处理应力奇异问题的单元材料破损约束松弛关系；

e＝1，2，Λ，N；r＝1，2，Λ，L    (2)

式中ρ_e是当前迭代步有效结构第e单元的拓扑变量，它是第e单元的相对密度；ρ_e与第e单元的弹性模量E_e的关系为

是当前迭代步有效结构第e单元满材料时的弹性模量；υ是预先给定的惩罚参数经验值，取υ＝3；是当前迭代步有效结构第e单元满材料时的允许von Mises破坏应力；

为第r组载荷作用下当前迭代步有效结构第e单元的von Mises应力，L为载荷组数，N为当前迭代步初始有效结构单元个数；w为一预先给定的经验参数，w取为0.7；取当前迭代步初始有效结构中可设计部件的单元数为Q，其单元编号设为n_q(q＝1，2，Λ，Q)；不设计部件的单元数为P，其单元编号设为i_p(p＝1，2，Λ，P)，当前迭代步初始有效结构所有保留单元编号为s_i(i＝1，2，Λ，P+R)，当前迭代步初始有效结构所有保留单元个数为P+R；

由于N相当大，一般几千，甚至几万，如果采用需要式(2)约束函数导数的数学优化方法，其工作量太大以致无法完成优化求解；故需采用少量约束函数替代式(2)以便大量减少约束函数数目，在此情况下常常出现结构应力集中问题、局部应力难以控制和最大应力波动问题；为了解决结构应力集中问题和使得优化拓扑有较好的黑/白分布特征，根据所有实际载荷作用下当前迭代步初始有效结构应力和单元材料破损应力阀值

，将式(3)中设计拓扑变量的离散条件和式(4)中应力的q₁范数度量函数

作为二个罚函数，构造新的目标函数式(5)；

$\sqrt{{\mathrm{\ρ}}_{n_v}}(1-\sqrt{{\mathrm{\ρ}}_{n_v}})=0$ (v＝1，2，Λ，Q)    (3)

$f_{\mathrm{so}}^r={\left(\underset{i=1}{\overset{P+R}{\mathrm{\Σ}}}{\left(\frac{{\mathrm{\σ}}_{\mathrm{vm}}^{s_i,r}}{g^w\left({\mathrm{\ρ}}_{s_i}\right){\mathrm{\σ}}_Y^{s_i}}\right)}^{q_1}{\mathrm{\ρ}}_{s_i}{V}_0^{s_i}\right)}^{1/q_1}$ (r＝1，2，Λ，L)    (4)

min： $V+\left(\mathrm{\γ}\underset{v=1}{\overset{Q}{\mathrm{\Σ}}}(1-\sqrt{{\mathrm{\ρ}}_{n_v}})\sqrt{{\mathrm{\ρ}}_{n_v}}{V}_0^{n_v}\right)+{\mathrm{\γ}}^t\left(\underset{r=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\α}}_r{f}_{\mathrm{so}}^r\right)V^{(l-1)}/\left(\underset{r=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\α}}_r{f}_{\mathrm{so}}^{r,(l-1)}\right)---\left(5\right)$

式中

和

分别是当前迭代步初始有效结构第s_i单元和第n_v单元满材料时的体积，V^(l-1)是当前迭代步初始有效结构的体积，V是当前迭代步有效结构体积；

是当前迭代步初始有效结构的；q₁是预先给定的经验值，取q₁＝18.；α_r是第r组实际载荷下结构应力的加权参数，且

选取α_r＝1/L(r＝1，2，Λ，L)；γ^t和γ是另外两个加权经验参数，分别取γ^t＝0.1和γ＝0.55；

为了大量减少应力约束数，首先确定所有载荷作用下少量最潜在的主动应力约束所对应的区域与单元信息；根据所有实际载荷作用下当前迭代步初始有效结构的应力、单元材料破损应力阀值

和其拓扑变量，按下述四个步骤选取第r组实际载荷下当前迭代步有效结构的N_c个子区域与相应的最潜在的主动应力约束单元信息；N_c是一个预先指定的当前迭代步有效结构子区域的个数，选取N_c＝5；

第一步骤：选取

作为第一个子区域所含的当前迭代步初始有效结构的所有保留单元集合；

第二步骤：根据第r组实际载荷作用下当前迭代步初始有效结构的应力、单元材料破损应力阀值

和其拓扑变量，按式(6)确定

中最潜在的主动应力约束单元序号

，并由式(7)确定当前迭代步有效结构的第二个子区域

$\underset{s_j\∈U_1^r}{\max }\left(\frac{{\mathrm{\σ}}_{\mathrm{vm}}^{s_j,r}}{g^w\left({\mathrm{\ρ}}_{s_j}\right){\mathrm{\σ}}_Y^{s_j}}{|}_{{\mathrm{\ρ}}^{(l-1)}}\right)=\frac{{\mathrm{\σ}}_{\mathrm{vm},(l-1)}^{s_1^r,r}}{g^w\left({\mathrm{\ρ}}_{s_1^r}^{(l-1)}\right){\mathrm{\σ}}_Y^{s_{\mathrm{ma}}}}$ and $s_1^r\∈U_1^r---\left(6\right)$

$U_2^r=\left\{s_i\right|\left|\right|x\left(s_i\right)-x\left(s_1^r\right)\left|\right|>\mathrm{dand}{s}_i\∈U_1^r\}---\left(7\right)$

式中

为第r组载荷作用下当前迭代步初始有效结构第单元的vonMises应力；x(s_i)是当前迭代步初始有效结构第s_i单元几何中心点的坐标矢量，

是当前迭代步初始有效结构第s_i单元和第

单元的几何中心点距离；d是一预先指定的参数，这里采用d＝9.d_mesh，而d_mesh是结构单元网格线间的最小空间距离；

第三步骤：按第二步骤的方法，根据式(8)和式(9)，确定

中最潜在的主动应力约束单元序号

和第r组实际载荷下当前迭代步有效结构的第k个子区域

第四步骤：按第二步骤的方法，根据递推法按式(8)和式(9)确定第r组实际载荷下当前迭代步有效结构的所有的

以及

$s_2^r,...,s_{N_c}^r;$

$\underset{s_j\∈U_{(k-1)}^r}{\max }\left(\frac{{\mathrm{\σ}}_{\mathrm{vm}}^{s_j,r}}{g^w\left({\mathrm{\ρ}}_{s_j}\right){\mathrm{\σ}}_Y^{s_j}}{|}_{{\mathrm{\ρ}}^{(l-1)}}\right)=\frac{{\mathrm{\σ}}_{\mathrm{vm},(l-1)}^{s_{(k-1)}^r,r}}{g^w\left({\mathrm{\ρ}}_{s_{(k-1)}^r}^{(l-1)}\right){\mathrm{\σ}}_Y^{s_{(k-1)}^r}}$ and $s_{(k-1)}^r\∈U_{(k-1)}^r,$ k＝3，4，Λ，N_c  (8)

$U_k^r=\left\{s_i\right|\left|\right|x\left(s_i\right)-x\left(s_{(k-1)}^r\right)\left|\right|>\mathrm{dand}{s}_i\∈U_{(k-1)}^r\},$ k＝3，4，Λ，N_c；r＝1，2，Λ，L  (9)

式中

是当前迭代步初始有效结构第s_i单元和第

单元的几何中心点距离；

为了减少结构最大应力波动及控制局部应力水平，需引进结构体积约束和相应的变上、下限以便使结构拓扑变量在一小范围变化；根据V^(l-1)和最大设计区域满材料时的结构体积V₀，采用式(10)计算当前迭代步有效结构体积的变上限V^U，l、变下限V^L，l；

V^U，l＝min(V^(l-1)(1.+β₁)，V₀)，V^U，l＝max(V^(l-1)(1.-β₁)，0)    (10)

式中β₁为一预先给定的经验值，选取β₁＝0.01；

在控制结构最大应力波动及控制局部应力水平方面，还需在替代原应力约束的等效应力约束中引入变约束限

和

，以使得相关的单元应力小量变化；根据实际载荷作用下当前迭代步初始有效结构的应力、单元材料破损应力阀值

、拓扑变量，按式(11a)、(11b)、(11c)确定γ_r，1、γ_r，k(k＝1，2，Λ，N_c；r＝1，2，Λ，L)和

(r＝1，2，Λ，L)，并按式(12a、12b)计算后继步所需的变应力约束限

(k＝1，2，Λ，N_c；r＝1，2，Λ，L)和

(r＝1，2，Λ，L)；

${\mathrm{\γ}}_{r,1}=\underset{i=\mathrm{1,2},\mathrm{\Λ},P+R}{\max }\left(\frac{{\mathrm{\σ}}_{\mathrm{vm}}^{s_i,r}}{g^w\left({\mathrm{\ρ}}_{s_i}\right){\mathrm{\σ}}_Y^{s_i}}{|}_{{\mathrm{\ρ}}^{(l-1)}}\right)=\frac{{\mathrm{\σ}}_{\mathrm{vm},(l-1)}^{s_1^r,r}}{g^w\left({\mathrm{\ρ}}_{s_1^r}^{(l-1)}\right){\mathrm{\σ}}_Y^{s_1^r}}---\left(11a\right)$

${\mathrm{\γ}}_{r,k}=\underset{s_i\∈U_k^r}{\max }\left(\frac{{\mathrm{\σ}}_{\mathrm{vm}}^{s_i,r}}{g^w\left({\mathrm{\ρ}}_{s_i}\right){\mathrm{\σ}}_Y^{s_i}}{|}_{{\mathrm{\ρ}}^{(l-1)}}\right)=\frac{{\mathrm{\σ}}_{\mathrm{vm},(l-1)}^{s_k^r,r}}{g^w\left({\mathrm{\ρ}}_{s_k^r}^{(l-1)}\right){\mathrm{\σ}}_Y^{s_k^r}}---\left(11b\right)$

${\mathrm{\γ}}_r^2={\left(\frac{1}{P+R}\underset{i=1}{\overset{P+R}{\mathrm{\Σ}}}{\left(\frac{{\mathrm{\σ}}_{\mathrm{vm},(l-1)}^{s_i,r}}{g^w\left({\mathrm{\ρ}}_{s_i}^{(l-1)}\right){\mathrm{\σ}}_Y^{s_i}}\right)}^{q_2}\right)}^{1/q_2}---\left(11c\right)$

${\mathrm{\θ}}_{r,k}^{l,1}=\left\{\begin{array}{cc}\min (({\mathrm{\γ}}_{r,k}+\min ({\mathrm{\β}}_1{\mathrm{\γ}}_{r,k},1-{\mathrm{\γ}}_{r,k})),1.0),& {\mathrm{\γ}}_{r,k}\≤1\\ \max (({\mathrm{\γ}}_{r,k}-\min ({\mathrm{\β}}_1{\mathrm{\γ}}_{r,k},{\mathrm{\γ}}_{r,k}-1)),1.0),& {\mathrm{\γ}}_{r,k}>1\end{array}\right.---\left(12a\right)$

${\mathrm{\θ}}_r^{l,2}=\left\{\begin{array}{cc}\min (({\mathrm{\γ}}_r^2+\min ({\mathrm{\β}}_2{\mathrm{\γ}}_r^2,1-{\mathrm{\γ}}_r^2)),1.0),& {\mathrm{\γ}}_r^2\≤1\\ \max (({\mathrm{\γ}}_r^2-\min ({\mathrm{\β}}_2{\mathrm{\γ}}_r^2,{\mathrm{\γ}}_r^2-1)),1.0),& {\mathrm{\γ}}_r^2>1\end{array}\right.---\left(12b\right)$

式中β₂为一预先给定的经验值，选取β₂＝0.01；

综合上面的考虑，为了大量减少应力约束数，解决结构应力集中问题、减少结构最大应力波动及控制局部应力水平，根据当前迭代步初始有效结构应力、单元材料破损应力阀值和变约束限，用式(13)中第五项的应力q₂范数度量函数约束、式(13)中第四项的N_c个当前迭代步有效结构子区域上最潜在的主动单元应力函数约束和引入的体积约束V≤V^U，l及V≥V^L，l替代式(2)中所有的单元应力约束，以及结合变约束限V^U，l、V^L，l、

和

并与式(5)的目标函数相结合，形成以当前迭代步有效结构体积为目标函数，其单元应力小于或等于材料破坏应力为约束条件的连续体结构拓扑优化设计的等效近似优化模型；即形成近似优化模型(13)；

$\left\{\begin{array}{c}\min :V+\left(\mathrm{\γ}\underset{v=1}{\overset{Q}{\mathrm{\Σ}}}(1-\sqrt{{\mathrm{\ρ}}_{n_v}})\sqrt{{\mathrm{\ρ}}_{n_v}}{V}_0^{n_v}\right)+{\mathrm{\γ}}^t\left(\underset{r=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\α}}_r{f}_{\mathrm{so}}^r\right)V^{(l-1)}/\left(\underset{r=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\α}}_r{f}_{\mathrm{so}}^{r,(l-1)}\right)\\ s.t.:V\≤V^{U,l}\\ V\≥V^{L,l}\\ \underset{s_i\∈U_k^r}{\max }\left(\frac{{\mathrm{\σ}}_{\mathrm{vm}}^{s_i,r}}{g^w\left({\mathrm{\ρ}}_{s_i}\right){\mathrm{\σ}}_Y^{s_i}}\right)\≤{\mathrm{\θ}}_{r,k}^{l,1}(r=\mathrm{1,2},\mathrm{\Λ},L;k=\mathrm{1,2},\mathrm{\Λ},N_c)\\ f_{\mathrm{sc}}^r={\left(\frac{1}{P+R}\underset{i=1}{\overset{P+R}{\mathrm{\Σ}}}{\left(\frac{{\mathrm{\σ}}_{\mathrm{vm}}^{s_i,r}}{g^w\left({\mathrm{\ρ}}_{s_i}\right){\mathrm{\σ}}_Y^{s_i}}\right)}^{q_2}\right)}^{1/q_2}\\ {\stackrel{\‾}{\mathrm{\ρ}}}_{n_v}^1\≤{\mathrm{\ρ}}_{n_v}\≤{\stackrel{\‾}{\mathrm{\ρ}}}_{n_v}=1,v=\mathrm{1,2},\mathrm{\Λ},Q\end{array}\right.---\left(13\right)$

式中q₂是预先给定的经验值，选取q₂＝6.；

为确保当前迭代步有效结构非奇异的小量经验值，依据实体各向同性材料惩罚法(简称SIMP法)，选取

4.2)第二个优化阶段的近似优化模型建模

为了使得第二个优化阶段的优化拓扑有较好的黑/白分布特征，以及改进优化效率，需对取(0，1)中间值的结构拓扑变量进行较大惩罚，同时需删除具有较大拓扑变量值的单元；按照第一个优化阶段和阶段转化步的近似优化模型建模方法和公式进行第二个优化阶段的近似优化模型建模，其中：

将式(2)中的

的υ以及ρ_e与第e单元的弹性模量E_e的关系式

中的υ改为υ₁，选取υ₁＝6；将式(4)和式(13)中的参数q₁和q₂分别取为q₁＝5和q₂＝2；式(5)和近似优化模型(13)的目标函数中的一个加权经验参数γ赋零；

第五步：形成近似的二次数学规划模型及求解

对于有大量设计变量和少量约束函数的优化问题，二次数学规划求解法的计算效率是非常高的；这里先将近似优化模型(13)转化为二次数学规划模型；根据实际载荷和虚载荷作用下当前迭代步初始有效结构的位移、应力数据和其拓扑变量，将近似优化模型(13)的约束条件展开成拓扑变量倒变量的约束条件一次显式近似式；根据移动渐近线方法，更新拓扑设计变量上限、下限，并形成近似优化模型(13)中目标函数的第三项关于拓扑变量的一次显式近似式，再将该一次显式近似式对拓扑变量倒变量求一次导数和二次导数，建立近似优化模型(13)中目标函数的第三项关于拓扑变量倒变量的二次显式近似式；根据式(14)的当前迭代步有效结构的体积函数以及约束条件一次显式近似式和近似优化模型(13)中目标函数的第三项关于拓扑变量倒变量的二次显式近似式，将近似优化模型(13)转化为近似的二次数学规划模型；

$V\≈\underset{v=1}{\overset{Q}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{n_v}{V}_0^{n_v}+\underset{p=1}{\overset{P}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{i_p}{V}_0^{i_p}---\left(14\right)$

采用实体各向同性材料惩罚法(简称SIMP法)的解决棋盘格问题技术，修改该二次数学规划模型的系数项；而后，采用库塔克准则和对偶二次规划法求解获得该近似二次规划模型的拓扑变量解；

第六步：结构设计空间的减缩和扩展

优化过程分为二个优化阶段和一个阶段转化步；刚开始优化迭代时，进行第一个优化阶段的结构设计空间的减缩和扩展；当式(1)成立时，进行阶段转化步的结构设计空间的减缩和扩展，后继的所有迭代步都进行第二个优化阶段的结构设计空间的减缩和扩展；

6.1)第一优化阶段结构设计空间的减缩和扩展

为了使得优化拓扑有较好的黑/白分布特征，以及易于求解具有较大规模有限元网格模型的结构最大设计域的优化问题，采用在一些优化迭代步将已删除的单元从结构上去掉，而其他单元拓扑变量不变的策略；即根据给定的小量ε＝0.07和

，设置一个经验阀值

，使得

(q＝1，2，Λ，Q)，选取经验阀值根据第五步获得的拓扑变量解，挑选出最潜在的可删去单元集合

式中h_i为当前迭代步初始有效结构上保留单元编号，

为第五步获得的当前迭代步的第h_i单元拓扑变量；

另外，根据第五步获得的拓扑变量解，按式(15)计算

，并按式(16)确定人工材料单元中最潜在的需增添单元集合N_a；

${\mathrm{\ρ}}_a^m=f_k^{-1}(\underset{j=1}{\overset{m^h}{\mathrm{\Σ}}}{f}_k\left({\mathrm{\ρ}}_{m_j}^{\left(l\right)}\right)/m^h)---\left(15\right)$

$N_a=\left\{m_j\right|{\mathrm{\ρ}}_{m_j}^{\left(l\right)}\≥\max ({\mathrm{\ρ}}_a^m,{\stackrel{\‾}{\mathrm{\ρ}}}_m^1)\mathrm{and}({\mathrm{\ρ}}_{m_j}^{\left(l\right)}-{\mathrm{\ρ}}_{m_j}^{(l-1)})\≥0.0\}---\left(16\right)$

式中m^h为当前迭代步初始有效结构上人工材料单元个数，m_j为当前迭代步初始有效结构上人工材料单元编号，为第五步获得的当前迭代步的第m_j单元拓扑变量，

是函数

的逆函数；

设计空间的减缩和扩展操作公式表示为式(17)：

${\mathrm{\ρ}}_{h_i}^{*}=\left\{\begin{array}{cc}{\mathrm{\ρ}}_{h_i}^{\left(l\right)}& h_i\∉N_d\\ 0& h_i\∈N_d\end{array}\right.,$ ${\mathrm{\ρ}}_{m_j}^{*}=\left\{\begin{array}{cc}{\mathrm{\ρ}}_{m_j}^{\left(l\right)}& m_j\∈N_a\\ 0& m_j\∉N_a\end{array}\right.---\left(17\right)$

设计空间的减缩和扩展准则表示为式(18)：

N_d∪N_a≠0和k_nv≥3    (18)

式中k_nv是施加式(17)操作后的迭代步数；

当式(18)满足时，则执行方程(17)的操作；将设计空间减缩和扩展得到的结构拓扑变量矢量仍记为ρ^(l)，置迭代步数l＝l+1，转向第二步继续优化迭代；

当式(18)不满足时，不进行设计空间的减缩和扩展；将第五步获得的结构拓扑变量矢量记为ρ^(l)，置迭代步数l＝l+1，转向第二步继续优化迭代；

通过上述的技术完成设计空间的减缩和扩展，确保了结构拓扑变量的小量变化，在一定程度上维持了迭代方法的收敛特性，同时又大大减少了每一迭代步有效结构特性分析和其灵敏度分析量；同时，第四步和该步的技术，如变位移限技术、引进目标函数惩罚项和设计空间的减缩和扩展，在一定程度确保了相邻几轮迭代步的结构拓扑处于某个具有较好黑/白分布特征的结构拓扑的小邻域里，从而有利于获得有较好黑/白分布特征的优化拓扑；

6.2)阶段转化步结构设计空间的减缩和扩展

当第一优化阶段结构设计空间调整结束时，当前迭代步有效结构上可能有大量灰色区域；为了改进优化效率和减少灰色区域，建立了阶段转化步结构设计空间的减缩和扩展操作，该操作类似于SIMP法的最后拓扑解释工作；它的初始结构是第一优化阶段获得的最佳结构；根据第五步获得的拓扑变量解，按式(19)确定单元集合

${\stackrel{\&OverBar;}{N}}_d^T=\left\{h_i\right|{\mathrm{\&rho;}}_{h_i}^{\left(l\right)}<{\mathrm{\&rho;}}_d^{\mathrm{Tp}}\}---\left(19\right)$

式中

是一预先指定的阀值，选取

该阶段设计空间的减缩和扩展操作公式表示为：

${\mathrm{\&rho;}}_{h_i}^{*}=\left\{\begin{array}{cc}1& h_i\&NotElement;{\stackrel{\&OverBar;}{N}}_d^T\\ 0& h_i\&Element;{\stackrel{\&OverBar;}{N}}_d^T\end{array}\right.---\left(20\right)$

根据

，执行方程(20)的操作，并将当前迭代步有效结构所有人工材料单元的拓扑变量赋零，获得的优化结构作为第二个优化求解阶段的初始结构；

将设计空间减缩和扩展得到的结构拓扑变量矢量仍记为ρ^(l)，置迭代步数l＝l+1，转向第二步继续优化迭代。

6.3)第二个优化求解阶段的结构设计空间减缩和扩展

当第一个优化求解阶段迭代解接近结构最佳拓扑邻域时，求解收敛很慢，为了改进求解效率，提出了如下类于梯度法的拓扑变量更新技术；根据第五步获得的拓扑变量解，按式(21a)和(21b)计算两个参数

而后按式(22)确定单元集合

和

，形成的单元集合

含有较大拓扑变量值的单元；

${\mathrm{\&rho;}}_d^{\left(l\right)}=\frac{1}{2}\{f_k^{-1}(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{f}_k\left({\mathrm{\&rho;}}_i^{\left(l\right)}\right)/N)+\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&rho;}}_i^{\left(l\right)}/N\}---\left(21a\right)$

${\mathrm{\&rho;}}_a^{\left(l\right)}=\frac{1}{2}\{f_k^{-1}(\underset{j=1}{\overset{P+R}{\mathrm{\&Sigma;}}}{f}_k\left({\mathrm{\&rho;}}_{s_j}^{\left(l\right)}\right)/(P+R))+\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&rho;}}_i^{\left(l\right)}/N\}---\left(21b\right)$

${\stackrel{\&OverBar;}{N}}_d=\left\{h_i\right|{\mathrm{\&rho;}}_{h_i}^{\left(l\right)}\&le;{\mathrm{\&rho;}}_d^{\left(l\right)}\},$ ${\stackrel{\&OverBar;}{N}}_a=\left\{m_j\right|{\mathrm{\&rho;}}_{m_j}^{\left(l\right)}\&GreaterEqual;{\mathrm{\&rho;}}_a^{\left(l\right)}\}---\left(22\right)$

式中s_j(j＝1，2，Λ，(P+R))是当前迭代步初始有效结构的所有保留单元编号，P+R为当前迭代步初始有效结构所有保留单元个数；

为第五步获得的当前迭代步的第i单元拓扑变量；

是函数

的逆函数；

该阶段设计空间的减缩和扩展操作公式表示为式(23)：

${\mathrm{\&rho;}}_{h_i}^{*}=\left\{\begin{array}{cc}1& h_i\&NotElement;{\stackrel{\&OverBar;}{N}}_d\\ 0& h_i\&Element;{\stackrel{\&OverBar;}{N}}_d\end{array}\right.,$ ${\mathrm{\&rho;}}_{m_j}^{*}=\left\{\begin{array}{cc}1& m_j\&Element;{\stackrel{\&OverBar;}{N}}_a\\ 0& m_j\&NotElement;{\stackrel{\&OverBar;}{N}}_a\end{array}\right.---\left(23\right)$

该阶段设计空间的减缩和扩展准则表示为式(24)：

$n_d+n_a>r_s^1\&times;N,$ k_nv≥3    (24)

式中n_d、n_a分别为集合

和中单元个数，N为当前迭代步初始有效结构单元个数；

为一经验参数，常取

当完成方程(23)的操作时，则该迭代步有单元删去或单元增添，此时该迭代步又称为一个大循环迭代步，用l_d表示大循环迭代步数，此时置l_d＝l_d+1；当阶段转化步结束时，设置初始大循环迭代步数l_d＝1；

根据式(11a)中的γ_r，1(r＝1，2Λ，L)，按式(25)计算d_lim；

$d_{\lim }=\underset{r=\mathrm{1,2},\mathrm{\&Lambda;},L}{\max }\{1-{\mathrm{\&gamma;}}_{r,1}\}---\left(25\right)$

优化迭代终止准则表示为式(26)：

and d_lim≥0.0或d_lim≤d_con    (26)

式中

和

分别是第(l_d-1)和第l_d大循环迭代步优化获得的有效结构；ε_T和d_con是二个预先指定的经验值，选取ε_T＝0.001和d_con＝0.001；

当式(24)不满足时，不进行设计空间的减缩和扩展，此时大循环迭代步数不计数；将第五步获得的拓扑变量矢量记为ρ^(l)，置迭代步数l＝l+1，转向第二步继续优化迭代；

当式(24)满足时，则执行方程(23)的操作；此时检查优化迭代终止条件式(26)，当式(26)成立时，则停止优化求解，获得最佳结构拓扑；当式(26)不满足时，将设计空间的减缩和扩展得到的结构拓扑变量矢量仍记为ρ^(l)，置迭代步数l＝l+1，l_d＝l_d+1，转向第二步继续优化迭代。

本发明的有益效果是：由于采用有理分式近似材料模型和引入材料破坏应力与单元拓扑变量的关系函数，以及用结构应力的q₁范数度量函数作为罚函数，有效地解决了应力奇异和结构应力集中问题；将结构应力的q₂范数度量函数约束、选取的几个子区域上最潜在的主动单元应力约束和引入的体积约束替代所有的单元应力约束，并结合变约束限措施能较精确地控制局部应力水平、大幅度减少了应力约束数及减少结构最大应力波动；利用结构设计空间优化概念和考虑密度变量的离散条件产生黑/白分布的结构拓扑设计，自动地实现了设计空间的减缩和扩展，减少了结构特性分析和其灵敏度分析量；优化结构的单元规模变化及二次规划法的采用较大地提高了优化效率。建立的基于材料破损约束的连续体结构拓扑设计近似模型能近似反映原局部应力约束的主要特征。建立的模型和提出的求解算法能解决应力集中问题和应力奇异问题，并减小结构特性分析和其灵敏度分析量，提高优化设计效率。

附图说明

附图1是本发明优化方法的流程图；

附图2是实施例一个L型梁结构的初始设计域；

附图3是采用本发明方法获得的实施例的最佳拓扑；

附图4是附图3的最佳拓扑结构的von Mises应力云图。

下面结合具体实施方式对本发明作详细说明

具体实施方式

实施例：一个L型梁结构的拓扑优化

一个长1米、高1米、厚度0.01米的L型梁示于图2，其上端面固定。一垂直载荷F＝55.kN均匀地作用于右铅垂端最上部0.04167米的区域内，且不考虑重力作用。梁的破坏应力为σ^Y＝355Mpa，其杨氏模量为E＝210GPa，泊松比v＝0.3，质量密度是ρ₀＝7800kg/m³；最大设计区域结构作为设计用的初始结构，初始结构体积为0.0064m³。

优化设计目标：体积最小，最大Von Mises应力小于破坏应力。最大设计域划分为9216个相同的四边形平面应力单元。

优化步骤如下：

1、根据需要设计的连续体结构的服役要求，将图2所示L型梁结构的初始设计域作为最大设计区域，利用尺寸相同单元的规则网格将该最大设计区域划分为9216个相同的四边形平面应力单元；将该L型梁结构最上端面的所有节点的二个自由度固定支持，设置其他边界面的所有节点的二个自由度自由；按均布载加载方式在该L型梁结构右铅垂端最上部0.04167米的区域内施加总载荷值为F＝55.kN的均布载荷；L型梁结构9216个单元都设为可设计单元；将其所有拓扑变量初值都设置为1.，将L型梁结构每个单元的杨氏模量、泊松比、质量密度和破坏应力分别设为E＝210GPa、v＝0.3、ρ0＝7800kg/m³和σ^Y＝355Mpa；形成最大设计区域固定有限单元网格模型，最大设计区域结构作为设计用的初始结构，并形成设计用的初始结构保留单元数据文件；

2、在结构周边和孔洞周围设置一层人工材料单元，形成当前迭代步初始有效结构的有限单元网格模型，建立最大设计区域有限单元网格模型与当前迭代步初始有效结构有限单元网格模型的映射关系数据文件；

3、获得实际载荷下当前迭代步初始有效结构有限单元网格模型的结构位移、应力，生成基于伴随方程的应力灵敏度计算所需的虚载荷，计算虚载荷作用下当前迭代步初始有效结构位移、应力；读取所述映射关系数据，获得最大设计区域固定有限单元网格模型的结构位移、应力数据；

4、将当前迭代步有效结构可设计拓扑变量的离散条件和应力的q₁范数度量函数作为二个罚函数，用当前迭代步有效结构应力的q₂范数度量函数约束、自动选取的N_c个当前迭代步有效结构子区域上最潜在的主动单元应力函数约束和引入的体积约束替代原优化问题中所有的单元应力约束，以及结合变约束限方案，形成优化问题的形成近似优化模型，并获得近似二次规划模型。再按避免棋盘格的算法修正二次规划模型的系数项。

5、采用库塔克准则及对偶二次规划求解算法求解二次规划模型的拓扑变量解。

6、进行设计空间的减缩和扩展，更新结构拓扑变量，形成新的结构系统有限元模型。如此重复上述第二步至第六步进行结构拓扑优化迭代直至最终收敛到优化结果。

以本实验室的Intel core i7920(2.66G CPU速度)为工作平台，以上述L梁结构拓扑优化为例，分别采用已有方法和本发明方法获得的最佳结构拓扑的体积、最大应力、所需迭代步数以及优化所需总时间对比如下表所示。本发明方法获得的最佳结构拓扑和其von Mises应力云图见图3和图4。

由本表、图3和图4可见，基于材料破损约束的连续体结构拓扑设计建模及优化设计方法获得的最佳拓扑具有较小的体积，且没有应力集中现象。最佳结构体积为初始体积的0.4203的情况下，最大von Mises应力下降了24.57％，且在破坏应力之下。尽管本发明方法需较多的迭代步，但由于结构有限元分析、应力灵敏度分析和优化求解变量规模大幅度减小，总的计算时间减少了10.32％。本发明方法解决了应力集中问题和应力奇异问题，能较精确地控制局部应力水平及减少结构最大应力波动，减少了结构分析和灵敏度分析的规模，可获得纯的黑/白分布的最佳拓扑。

Claims (1)

1.基于材料破损约束的连续体结构拓扑设计建模及优化设计方法，包括下述步骤：

第一步：有限元初始建模

根据需要设计的连续体结构的服役要求，确定一个拓扑结构的足够大设计区域，利用规则网格对该设计区域进行有限元划分，确定边界条件、载荷工况、结构拓扑变量初值和单元材料特性值，形成最大设计区域固定有限单元网格模型，形成设计用的初始结构并形成其保留单元数据文件；根据需要设计的连续体结构的用材，确定材料破损应力阀值；设置初始迭代步数l＝1；则设计用的初始结构为当前迭代步初始结构；

第二步：有限元模型更新

读取固定有限单元网格模型文件和当前迭代步初始结构保留单元数据文件，在结构周边和孔洞周围设置一层人工材料单元，形成当前迭代步初始有效结构并形成其有限单元网格模型，建立最大设计区域有限单元网格模型与当前迭代步初始有效结构有限单元网格模型的映射关系，得到映射关系数据文件；此时的当前迭代步初始有效结构也为上一迭代步优化获得的有效结构；

第三步：当前迭代步初始有效结构有限元分析

根据当前迭代步初始有效结构有限单元网格模型数据进行有限元分析，获得实际载荷作用下当前迭代步初始有效结构有限单元网格模型的结构位移、应力，生成基于伴随方程的应力灵敏度计算所需的虚载荷，计算虚载荷作用下当前迭代步初始有效结构位移、应力；读取所述映射关系数据文件，获得最大设计区域固定有限单元网格模型的结构位移、应力数据；

第四步：形成近似优化模型

在该迭代步中，当前迭代步初始有效结构随着结构拓扑变量变化将发生变化，其网格模型不变，其拓扑变量、位移和应力都将发生变化；则变化了拓扑变量的当前迭代步初始有效结构称为当前迭代步有效结构；

优化过程分为二个优化阶段和一个阶段转化步；刚开始优化迭代时，进行第一个优化阶段的近似优化模型建模；当式(1)成立时，进行阶段转化步的近似优化模型建模；后继的所有迭代步都进行第二个优化阶段的近似优化模型建模；

l≥l_ap，和|(V^(l-1)-V^(l-2))/V^(l-1)|＜ε_f    (1)

式中l是当前的迭代步数，l_ap是预先指定的限值，选取l_ap＝45；V^(l-1)和V^(l-2)分别是第(l-1)迭代步和第(l-2)迭代步优化获得的有效结构体积；ε_f是一预先指定的经验值，选取ε_f＝0.001；

4.1)第一个优化阶段和阶段转化步的近似优化模型建模

根据实际载荷和虚载荷作用下当前迭代步初始有效结构位移、应力数据，根据式(2)建立处理应力奇异问题的单元材料破损约束松弛关系；

e＝1，2，Λ，N；r＝1，2，Λ，L    (2)

式中ρ_e是当前迭代步有效结构第e单元的拓扑变量，它是第e单元的相对密度；ρ_e与第e单元的弹性模量E_e的关系为

是当前迭代步有效结构第e单元满材料时的弹性模量；υ是预先给定的惩罚参数经验值，取υ＝3；

是当前迭代步有效结构第e单元满材料时的允许von Mises破坏应力；

为第r组载荷作用下当前迭代步有效结构第e单元的von Mises应力，L为载荷组数，N为当前迭代步初始有效结构单元个数；w为一预先给定的经验参数，w取为0.7；取当前迭代步初始有效结构中可设计部件的单元数为Q，其单元编号设为n_q(q＝1，2，Λ，Q)；不设计部件的单元数为P，其单元编号设为i_p(p＝1，2，Λ，P)，当前迭代步初始有效结构所有保留单元编号为s_i(i＝1，2，Λ，P+R)，当前迭代步初始有效结构所有保留单元个数为P+R；

根据所有实际载荷作用下当前迭代步初始有效结构应力和单元材料破损应力阀值

，将式(3)中设计拓扑变量的离散条件和式(4)中应力的q₁范数度量函数

作为二个罚函数，构造新的目标函数式(5)；

$\sqrt{{\mathrm{\&rho;}}_{n_v}}(1-\sqrt{{\mathrm{\&rho;}}_{n_v}})=0$ (v＝1，2，Λ，Q)    (3)

$f_{\mathrm{so}}^r={\left(\underset{i=1}{\overset{P+R}{\mathrm{\&Sigma;}}}{\left(\frac{{\mathrm{\&sigma;}}_{\mathrm{vm}}^{s_i,r}}{g^w\left({\mathrm{\&rho;}}_{s_i}\right){\mathrm{\&sigma;}}_Y^{s_i}}\right)}^{q_1}{\mathrm{\&rho;}}_{s_i}{V}_0^{s_i}\right)}^{1/q_1}$ (r＝1，2，Λ，L)    (4)

min： $V+\left(\mathrm{\&gamma;}\underset{v=1}{\overset{Q}{\mathrm{\&Sigma;}}}(1-\sqrt{{\mathrm{\&rho;}}_{n_v}})\sqrt{{\mathrm{\&rho;}}_{n_v}}{V}_0^{n_v}\right)+{\mathrm{\&gamma;}}^t\left(\underset{r=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_r{f}_{\mathrm{so}}^r\right)V^{(l-1)}/\left(\underset{r=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_r{f}_{\mathrm{so}}^{r,(l-1)}\right)---\left(5\right)$

式中

和分别是当前迭代步初始有效结构第s_i单元和第n_v单元满材料时的体积，V^(l-1)是当前迭代步初始有效结构的体积，V是当前迭代步有效结构体积；

是当前迭代步初始有效结构的

；q₁是预先给定的经验值，取q₁＝18.；α_r是第r组实际载荷下结构应力的加权参数，且

选取α_r＝1/L(r＝1，2，Λ，L)；γ^t和γ是另外两个加权经验参数，分别取γ^t＝0.1和γ＝0.55；

根据所有实际载荷作用下当前迭代步初始有效结构的应力、单元材料破损应力阀值

和其拓扑变量，按下述四个步骤选取第r组实际载荷下当前迭代步有效结构的N_c个子区域与相应的最潜在的主动应力约束单元信息；N_c是一个预先指定的当前迭代步有效结构子区域的个数，选取N_c＝5；

第一步骤：选取

作为第一个子区域所含的当前迭代步初始有效结构的所有保留单元集合；

第二步骤：根据第r组实际载荷作用下当前迭代步初始有效结构的应力、单元材料破损应力阀值

和其拓扑变量，按式(6)确定

中最潜在的主动应力约束单元序号

，并由式(7)确定当前迭代步有效结构的第二个子区域

$\underset{s_j\&Element;U_1^r}{\max }\left(\frac{{\mathrm{\&sigma;}}_{\mathrm{vm}}^{s_j,r}}{g^w\left({\mathrm{\&rho;}}_{s_j}\right){\mathrm{\&sigma;}}_Y^{s_j}}{|}_{{\mathrm{\&rho;}}^{(l-1)}}\right)=\frac{{\mathrm{\&sigma;}}_{\mathrm{vm},(l-1)}^{s_1^r,r}}{g^w\left({\mathrm{\&rho;}}_{s_1^r}^{(l-1)}\right){\mathrm{\&sigma;}}_Y^{s_{\mathrm{ma}}}}$ and $s_1^r\&Element;U_1^r---\left(6\right)$

$U_2^r=\left\{s_i\right|\left|\right|x\left(s_i\right)-x\left(s_1^r\right)\left|\right|>\mathrm{dand}{s}_i\&Element;U_1^r\}---\left(7\right)$

式中

为第r组载荷作用下当前迭代步初始有效结构第

单元的vonMises应力；x(s_i)是当前迭代步初始有效结构第s_i单元几何中心点的坐标矢量，是当前迭代步初始有效结构第s_i单元和第单元的几何中心点距离；d是一预先指定的参数，这里采用d＝9.d_mesh，而d_mesh是结构单元网格线间的最小空间距离；

第三步骤：按第二步骤的方法，根据式(8)和式(9)，确定

中最潜在的主动应力约束单元序号

和第r组实际载荷下当前迭代步有效结构的第k个子区域

第四步骤：按第二步骤的方法，根据递推法按式(8)和式(9)确定第r组实际载荷下当前迭代步有效结构的所有的

以及 $s_2^r,..,s_{N_c}^r;$

$\underset{s_j\&Element;U_{(k-1)}^r}{\max }\left(\frac{{\mathrm{\&sigma;}}_{\mathrm{vm}}^{s_j,r}}{g^w\left({\mathrm{\&rho;}}_{s_j}\right){\mathrm{\&sigma;}}_Y^{s_j}}{|}_{{\mathrm{\&rho;}}^{(l-1)}}\right)=\frac{{\mathrm{\&sigma;}}_{\mathrm{vm},(l-1)}^{s_{(k-1)}^r,r}}{g^w\left({\mathrm{\&rho;}}_{s_{(k-1)}^r}^{(l-1)}\right){\mathrm{\&sigma;}}_Y^{s_{(k-1)}^r}}$ and $s_{(k-1)}^r\&Element;U_{(k-1)}^r,$ k＝3，4，Λ，N_c  (8)

$U_k^r=\left\{s_i\right|\left|\right|x\left(s_i\right)-x\left(s_{(k-1)}^r\right)\left|\right|>\mathrm{dand}{s}_i\&Element;U_{(k-1)}^r\},$ k＝3，4，Λ，N_c；r＝1，2，Λ，L  (9)

式中

是当前迭代步初始有效结构第s_i单元和第

单元的几何中心点距离；

根据V^(l-1)和最大设计区域满材料时的结构体积V₀，采用式(10)计算当前迭代步有效结构体积的变上限V^U，l、变下限V^L，l；

V^U，l＝min(V^(l-1)(1.+β₁)，V₀)，V^U，l＝max(V^(l-1)(1.-β₁)，0)    (10)

式中β₁为一预先给定的经验值，选取β₁＝0.01；

根据实际载荷作用下当前迭代步初始有效结构的应力、单元材料破损应力阀值、拓扑变量，按式(11a)、(11b)、(11c)确定γ_r，1、γ_r，k(k＝1，2，Λ，N_c；r＝1，2，Λ，L)和

(r＝1，2，Λ，L)，并按式(12a、12b)计算后继步所需的变应力约束限

(k＝1，2，Λ，N_c；r＝1，2，Λ，L)和

(r＝1，2，Λ，L)；

${\mathrm{\&gamma;}}_{r,1}=\underset{i=\mathrm{1,2},\mathrm{\&Lambda;},P+R}{\max }\left(\frac{{\mathrm{\&sigma;}}_{\mathrm{vm}}^{s_i,r}}{g^w\left({\mathrm{\&rho;}}_{s_i}\right){\mathrm{\&sigma;}}_Y^{s_i}}{|}_{{\mathrm{\&rho;}}^{(l-1)}}\right)=\frac{{\mathrm{\&sigma;}}_{\mathrm{vm},(l-1)}^{s_1^r,r}}{g^w\left({\mathrm{\&rho;}}_{s_1^r}^{(l-1)}\right){\mathrm{\&sigma;}}_Y^{s_1^r}}---\left(11a\right)$

${\mathrm{\&gamma;}}_{r,k}=\underset{s_i\&Element;U_k^r}{\max }\left(\frac{{\mathrm{\&sigma;}}_{\mathrm{vm}}^{s_i,r}}{g^w\left({\mathrm{\&rho;}}_{s_i}\right){\mathrm{\&sigma;}}_Y^{s_i}}{|}_{{\mathrm{\&rho;}}^{(l-1)}}\right)=\frac{{\mathrm{\&sigma;}}_{\mathrm{vm},(l-1)}^{s_k^r,r}}{g^w\left({\mathrm{\&rho;}}_{s_k^r}^{(l-1)}\right){\mathrm{\&sigma;}}_Y^{s_k^r}}---\left(11b\right)$

${\mathrm{\&gamma;}}_r^2={\left(\frac{1}{P+R}\underset{i=1}{\overset{P+R}{\mathrm{\&Sigma;}}}{\left(\frac{{\mathrm{\&sigma;}}_{\mathrm{vm},(l-1)}^{s_i,r}}{g^w\left({\mathrm{\&rho;}}_{s_i}^{(l-1)}\right){\mathrm{\&sigma;}}_Y^{s_i}}\right)}^{q_2}\right)}^{1/q_2}---\left(11c\right)$

${\mathrm{\&theta;}}_{r,k}^{l,1}=\left\{\begin{array}{cc}\min (({\mathrm{\&gamma;}}_{r,k}+\min ({\mathrm{\&beta;}}_1{\mathrm{\&gamma;}}_{r,k},1-{\mathrm{\&gamma;}}_{r,k})),1.0),& {\mathrm{\&gamma;}}_{r,k}\&le;1\\ \max (({\mathrm{\&gamma;}}_{r,k}-\min ({\mathrm{\&beta;}}_1{\mathrm{\&gamma;}}_{r,k},{\mathrm{\&gamma;}}_{r,k}-1)),1.0),& {\mathrm{\&gamma;}}_{r,k}>1\end{array}\right.---\left(12a\right)$

${\mathrm{\&theta;}}_r^{l,2}=\left\{\begin{array}{cc}\min (({\mathrm{\&gamma;}}_r^2+\min ({\mathrm{\&beta;}}_2{\mathrm{\&gamma;}}_r^2,1-{\mathrm{\&gamma;}}_r^2)),1.0),& {\mathrm{\&gamma;}}_r^2\&le;1\\ \max (({\mathrm{\&gamma;}}_r^2-\min ({\mathrm{\&beta;}}_2{\mathrm{\&gamma;}}_r^2,{\mathrm{\&gamma;}}_r^2-1)),1.0),& {\mathrm{\&gamma;}}_r^2>1\end{array}\right.---\left(12b\right)$

式中β₂为一预先给定的经验值，选取β₂＝0.01；

根据当前迭代步初始有效结构应力、单元材料破损应力阀值和变约束限，用式(13)中第五项的应力q₂范数度量函数约束、式(13)中第四项的N_c个当前迭代步有效结构子区域上最潜在的主动单元应力函数约束和引入的体积约束V≤V^U，l及V≥V^L，l替代式(2)中所有的单元应力约束，以及结合变约束限V^U，l、V^L，l、

和

，并与式(5)的目标函数相结合，形成以当前迭代步有效结构体积为目标函数，其单元应力小于或等于材料破坏应力为约束条件的连续体结构拓扑优化设计的等效近似优化模型；即形成近似优化模型(13)；

$\left\{\begin{array}{c}\min :V+\left(\mathrm{\&gamma;}\underset{v=1}{\overset{Q}{\mathrm{\&Sigma;}}}(1-\sqrt{{\mathrm{\&rho;}}_{n_v}})\sqrt{{\mathrm{\&rho;}}_{n_v}}{V}_0^{n_v}\right)+{\mathrm{\&gamma;}}^t\left(\underset{r=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_r{f}_{\mathrm{so}}^r\right)V^{(l-1)}/\left(\underset{r=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_r{f}_{\mathrm{so}}^{r,(l-1)}\right)\\ s.t.:V\&le;V^{U,l}\\ V\&GreaterEqual;V^{L,l}\\ \underset{s_i\&Element;U_k^r}{\max }\left(\frac{{\mathrm{\&sigma;}}_{\mathrm{vm}}^{s_i,r}}{g^w\left({\mathrm{\&rho;}}_{s_i}\right){\mathrm{\&sigma;}}_Y^{s_i}}\right)\&le;{\mathrm{\&theta;}}_{r,k}^{l,1}(r=\mathrm{1,2},\mathrm{\&Lambda;},L;k=\mathrm{1,2},\mathrm{\&Lambda;},N_c)\\ f_{\mathrm{sc}}^r={\left(\frac{1}{P+R}\underset{i=1}{\overset{P+R}{\mathrm{\&Sigma;}}}{\left(\frac{{\mathrm{\&sigma;}}_{\mathrm{vm}}^{s_i,r}}{g^w\left({\mathrm{\&rho;}}_{s_i}\right){\mathrm{\&sigma;}}_Y^{s_i}}\right)}^{q_2}\right)}^{1/q_2}\\ {\stackrel{\&OverBar;}{\mathrm{\&rho;}}}_{n_v}^1\&le;{\mathrm{\&rho;}}_{n_v}\&le;{\stackrel{\&OverBar;}{\mathrm{\&rho;}}}_{n_v}=1,v=\mathrm{1,2},\mathrm{\&Lambda;},Q\end{array}\right.---\left(13\right)$

式中q₂是预先给定的经验值，选取q₂＝6.；

为确保当前迭代步有效结构非奇异的小量经验值，依据实体各向同性材料惩罚法(简称SIMP法)，选取

4.2)第二个优化阶段的近似优化模型建模

按照第一个优化阶段和阶段转化步的近似优化模型建模方法和公式进行第二个优化阶段的近似优化模型建模，其中：

将式(2)中的

的υ以及ρ_e与第e单元的弹性模量E_e的关系式中的υ改为υ₁，选取υ₁＝6；将式(4)和式(13)中的参数q₁和q₂分别取为q₁＝5和q₂＝2；式(5)和近似优化模型(13)的目标函数中的一个加权经验参数γ赋零；

第五步：形成近似的二次数学规划模型及求解

根据实际载荷和虚载荷作用下当前迭代步初始有效结构的位移、应力数据和其拓扑变量，将近似优化模型(13)的约束条件展开成拓扑变量倒变量的约束条件一次显式近似式；根据移动渐近线方法，更新拓扑设计变量上限、下限，并形成近似优化模型(13)中目标函数的第三项关于拓扑变量的一次显式近似式，再将该一次显式近似式对拓扑变量倒变量求一次导数和二次导数，建立近似优化模型(13)中目标函数的第三项关于拓扑变量倒变量的二次显式近似式；根据式(14)的当前迭代步有效结构的体积函数以及约束条件一次显式近似式和近似优化模型(13)中目标函数的第三项关于拓扑变量倒变量的二次显式近似式，将近似优化模型(13)转化为近似的二次数学规划模型；

$V\&ap;\underset{v=1}{\overset{Q}{\mathrm{\&Sigma;}}}{\mathrm{\&rho;}}_{n_v}{V}_0^{n_v}+\underset{p=1}{\overset{P}{\mathrm{\&Sigma;}}}{\mathrm{\&rho;}}_{i_p}{V}_0^{i_p}---\left(14\right)$

采用实体各向同性材料惩罚法(简称SIMP法)的解决棋盘格问题技术，修改该二次数学规划模型的系数项；而后，采用库塔克准则和对偶二次规划法获得该近似二次规划模型的拓扑变量解；

第六步：结构设计空间的减缩和扩展

优化过程分为二个优化阶段和一个阶段转化步；刚开始优化迭代时，进行第一个优化阶段的结构设计空间的减缩和扩展；当式(1)成立时，进行阶段转化步的结构设计空间的减缩和扩展，后继的所有迭代步都进行第二个优化阶段的结构设计空间的减缩和扩展；

6.1)第一优化阶段结构设计空间的减缩和扩展

根据给定的小量ε＝0.07和

，设置一个经验阀值

，使得(q＝1，2，Λ，Q)，选取经验阀值

根据第五步获得的拓扑变量解，挑选出最潜在的可删去单元集合

式中h_i为当前迭代步初始有效结构上保留单元编号，

为第五步获得的当前迭代步的第h_i单元拓扑变量；

另外，根据第五步获得的拓扑变量解，按式(15)计算

，并按式(16)确定人工材料单元中最潜在的需增添单元集合N_a；

${\mathrm{\&rho;}}_a^m=f_k^{-1}(\underset{j=1}{\overset{m^h}{\mathrm{\&Sigma;}}}{f}_k\left({\mathrm{\&rho;}}_{m_j}^{\left(l\right)}\right)/m^h)---\left(15\right)$

$N_a=\left\{m_j\right|{\mathrm{\&rho;}}_{m_j}^{\left(l\right)}\&GreaterEqual;\max ({\mathrm{\&rho;}}_a^m,{\stackrel{\&OverBar;}{\mathrm{\&rho;}}}_m^1)\mathrm{and}({\mathrm{\&rho;}}_{m_j}^{\left(l\right)}-{\mathrm{\&rho;}}_{m_j}^{(l-1)})\&GreaterEqual;0.0\}---\left(16\right)$

式中m^h为当前迭代步初始有效结构上人工材料单元个数，m_j为当前迭代步初始有效结构上人工材料单元编号，

为第五步获得的当前迭代步的第m_j单元拓扑变量，

是函数

的逆函数；

设计空间的减缩和扩展操作公式表示为式(17)：

${\mathrm{\&rho;}}_{h_i}^{*}=\left\{\begin{array}{cc}{\mathrm{\&rho;}}_{h_i}^{\left(l\right)}& h_i\&NotElement;N_d\\ 0& h_i\&Element;N_d\end{array}\right.,$ ${\mathrm{\&rho;}}_{m_j}^{*}=\left\{\begin{array}{cc}{\mathrm{\&rho;}}_{m_j}^{\left(l\right)}& m_j\&Element;N_a\\ 0& m_j\&NotElement;N_a\end{array}\right.---\left(17\right)$

设计空间的减缩和扩展准则表示为式(18)：

N_d∪N_a≠0和k_nv≥3    (18)

式中k_nv是施加式(17)操作后的迭代步数；

当式(18)满足时，则执行方程(17)的操作；将设计空间减缩和扩展得到的结构拓扑变量矢量仍记为ρ^(l)，置迭代步数l＝l+1，转向第二步继续优化迭代；

当式(18)不满足时，不进行设计空间的减缩和扩展；将第五步获得的结构拓扑变量矢量记为ρ^(l)，置迭代步数l＝l+1，转向第二步继续优化迭代；

6.2)阶段转化步结构设计空间的减缩和扩展

它的初始结构是第一优化阶段获得的最佳结构；根据第五步获得的拓扑变量解，按式(19)确定单元集合

${\stackrel{\&OverBar;}{N}}_d^T=\left\{h_i\right|{\mathrm{\&rho;}}_{h_i}^{\left(l\right)}<{\mathrm{\&rho;}}_d^{\mathrm{Tp}}\}---\left(19\right)$

式中

是一预先指定的阀值，选取

该阶段设计空间的减缩和扩展操作公式表示为：

${\mathrm{\&rho;}}_{h_i}^{*}=\left\{\begin{array}{cc}1& h_i\&NotElement;{\stackrel{\&OverBar;}{N}}_d^T\\ 0& h_i\&Element;{\stackrel{\&OverBar;}{N}}_d^T\end{array}\right.---\left(20\right)$

根据

，执行方程(20)的操作，并将当前迭代步有效结构所有人工材料单元的拓扑变量赋零，获得的优化结构作为第二个优化求解阶段的初始结构；

将设计空间减缩和扩展得到的结构拓扑变量矢量仍记为ρ^(l)，置迭代步数l＝l+1，转向第二步继续优化迭代。

6.3)第二个优化求解阶段的结构设计空间减缩和扩展

根据第五步获得的拓扑变量解，按式(21a)和(21b)计算两个参数

，而后按式(22)确定单元集合

和，形成的单元集合

含有较大拓扑变量值的单元；

${\mathrm{\&rho;}}_d^{\left(l\right)}=\frac{1}{2}\{f_k^{-1}(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{f}_k\left({\mathrm{\&rho;}}_i^{\left(l\right)}\right)/N)+\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&rho;}}_i^{\left(l\right)}/N\}---\left(21a\right)$

${\mathrm{\&rho;}}_a^{\left(l\right)}=\frac{1}{2}\{f_k^{-1}(\underset{j=1}{\overset{P+R}{\mathrm{\&Sigma;}}}{f}_k\left({\mathrm{\&rho;}}_{s_j}^{\left(l\right)}\right)/(P+R))+\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&rho;}}_i^{\left(l\right)}/N\}---\left(21b\right)$

${\stackrel{\&OverBar;}{N}}_d=\left\{h_i\right|{\mathrm{\&rho;}}_{h_i}^{\left(l\right)}\&le;{\mathrm{\&rho;}}_d^{\left(l\right)}\},$ ${\stackrel{\&OverBar;}{N}}_a=\left\{m_j\right|{\mathrm{\&rho;}}_{m_j}^{\left(l\right)}\&GreaterEqual;{\mathrm{\&rho;}}_a^{\left(l\right)}\}---\left(22\right)$

式中s_j(j＝1，2，Λ，(P+R))是当前迭代步初始有效结构的所有保留单元编号，P+R为当前迭代步初始有效结构所有保留单元个数；

为第五步获得的当前迭代步的第i单元拓扑变量；是函数

的逆函数；

该阶段设计空间的减缩和扩展操作公式表示为式(23)：

${\mathrm{\&rho;}}_{h_i}^{*}=\left\{\begin{array}{cc}1& h_i\&NotElement;{\stackrel{\&OverBar;}{N}}_d\\ 0& h_i\&Element;{\stackrel{\&OverBar;}{N}}_d\end{array}\right.,$ ${\mathrm{\&rho;}}_{m_j}^{*}=\left\{\begin{array}{cc}1& m_j\&Element;{\stackrel{\&OverBar;}{N}}_a\\ 0& m_j\&NotElement;{\stackrel{\&OverBar;}{N}}_a\end{array}\right.---\left(23\right)$

该阶段设计空间的减缩和扩展准则表示为式(24)：

$n_d+n_a>r_s^1\&times;N,$ k_nv≥3    (24)

式中n_d、n_a分别为集合

和

中单元个数，N为当前迭代步初始有效结构单元个数；

为一经验参数，常取

当完成方程(23)的操作时，则该迭代步有单元删去或单元增添，此时该迭代步又称为一个大循环迭代步，用l_d表示大循环迭代步数，此时置l_d＝l_d+1；当阶段转化步结束时，设置初始大循环迭代步数l_d＝1；

根据式(11a)中的γ_r，1(r＝1，2Λ，L)，按式(25)计算d_lim；

$d_{\lim }=\underset{r=\mathrm{1,2},\mathrm{\&Lambda;},L}{\max }\{1-{\mathrm{\&gamma;}}_{r,1}\}---\left(25\right)$

优化迭代终止准则表示为式(26)：

and d_lim≥0.0或d_lim≤d_con    (26)

式中

和

分别是第(l_d-1)和第l_d大循环迭代步优化获得的有效结构；ε_T和d_con是二个预先指定的经验值，选取ε_T＝0.001和d_con＝0.001；

当式(24)不满足时，不进行设计空间的减缩和扩展，此时大循环迭代步数不计数；将第五步获得的拓扑变量矢量记为ρ^(l)，置迭代步数l＝l+1，转向第二步继续优化迭代；

当式(24)满足时，则执行方程(23)的操作；此时检查优化迭代终止条件式(26)，当式(26)成立时，则停止优化求解，获得最佳结构拓扑；当式(26)不满足时，将设计空间的减缩和扩展得到的结构拓扑变量矢量仍记为ρ^(l)，置迭代步数l＝l+1，l_d＝l_d+1，转向第二步继续优化迭代。

Images