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

摘要

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

主权项

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{{\rho }_{n_v}}(1-\sqrt{{\rho }_{n_v}})=0$(v＝1，2，Λ，Q)    (3)$f_{\mathrm{so}}^r={\left(\underset{i=1}{\overset{P+R}{\Sigma }}{\left(\frac{{\sigma }_{\mathrm{vm}}^{s_i,r}}{g^w\left({\rho }_{s_i}\right){\sigma }_Y^{s_i}}\right)}^{q_1}{\rho }_{s_i}{V}_0^{s_i}\right)}^{1/q_1}$(r＝1，2，Λ，L)    (4)min：$V+\left(\gamma \underset{v=1}{\overset{Q}{\Sigma }}(1-\sqrt{{\rho }_{n_v}})\sqrt{{\rho }_{n_v}}{V}_0^{n_v}\right)+{\gamma }^t\left(\underset{r=1}{\overset{L}{\Sigma }}{\alpha }_r{f}_{\mathrm{so}}^r\right)V^{(l-1)}/\left(\underset{r=1}{\overset{L}{\Sigma }}{\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\in U_1^r}{\max }\left(\frac{{\sigma }_{\mathrm{vm}}^{s_j,r}}{g^w\left({\rho }_{s_j}\right){\sigma }_Y^{s_j}}{|}_{{\rho }^{(l-1)}}\right)=\frac{{\sigma }_{\mathrm{vm},(l-1)}^{s_1^r,r}}{g^w\left({\rho }_{s_1^r}^{(l-1)}\right){\sigma }_Y^{s_{\mathrm{ma}}}}$and$s_1^r\in 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\in 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\in U_{(k-1)}^r}{\max }\left(\frac{{\sigma }_{\mathrm{vm}}^{s_j,r}}{g^w\left({\rho }_{s_j}\right){\sigma }_Y^{s_j}}{|}_{{\rho }^{(l-1)}}\right)=\frac{{\sigma }_{\mathrm{vm},(l-1)}^{s_{(k-1)}^r,r}}{g^w\left({\rho }_{s_{(k-1)}^r}^{(l-1)}\right){\sigma }_Y^{s_{(k-1)}^r}}$and$s_{(k-1)}^r\in 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\in 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)；${\gamma }_{r,1}=\underset{i=\mathrm{1,2},\Lambda ,P+R}{\max }\left(\frac{{\sigma }_{\mathrm{vm}}^{s_i,r}}{g^w\left({\rho }_{s_i}\right){\sigma }_Y^{s_i}}{|}_{{\rho }^{(l-1)}}\right)=\frac{{\sigma }_{\mathrm{vm},(l-1)}^{s_1^r,r}}{g^w\left({\rho }_{s_1^r}^{(l-1)}\right){\sigma }_Y^{s_1^r}}---\left(11a\right)$${\gamma }_{r,k}=\underset{s_i\in U_k^r}{\max }\left(\frac{{\sigma }_{\mathrm{vm}}^{s_i,r}}{g^w\left({\rho }_{s_i}\right){\sigma }_Y^{s_i}}{|}_{{\rho }^{(l-1)}}\right)=\frac{{\sigma }_{\mathrm{vm},(l-1)}^{s_k^r,r}}{g^w\left({\rho }_{s_k^r}^{(l-1)}\right){\sigma }_Y^{s_k^r}}---\left(11b\right)$${\gamma }_r^2={\left(\frac{1}{P+R}\underset{i=1}{\overset{P+R}{\Sigma }}{\left(\frac{{\sigma }_{\mathrm{vm},(l-1)}^{s_i,r}}{g^w\left({\rho }_{s_i}^{(l-1)}\right){\sigma }_Y^{s_i}}\right)}^{q_2}\right)}^{1/q_2}---\left(11c\right)$${\theta }_{r,k}^{l,1}=\left\{\begin{array}{cc}\min (({\gamma }_{r,k}+\min ({\beta }_1{\gamma }_{r,k},1-{\gamma }_{r,k})),1.0),& {\gamma }_{r,k}\le 1\\ \max (({\gamma }_{r,k}-\min ({\beta }_1{\gamma }_{r,k},{\gamma }_{r,k}-1)),1.0),& {\gamma }_{r,k}>1\end{array}\right.---\left(12a\right)$${\theta }_r^{l,2}=\left\{\begin{array}{cc}\min (({\gamma }_r^2+\min ({\beta }_2{\gamma }_r^2,1-{\gamma }_r^2)),1.0),& {\gamma }_r^2\le 1\\ \max (({\gamma }_r^2-\min ({\beta }_2{\gamma }_r^2,{\gamma }_r^2-1)),1.0),& {\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(\gamma \underset{v=1}{\overset{Q}{\Sigma }}(1-\sqrt{{\rho }_{n_v}})\sqrt{{\rho }_{n_v}}{V}_0^{n_v}\right)+{\gamma }^t\left(\underset{r=1}{\overset{L}{\Sigma }}{\alpha }_r{f}_{\mathrm{so}}^r\right)V^{(l-1)}/\left(\underset{r=1}{\overset{L}{\Sigma }}{\alpha }_r{f}_{\mathrm{so}}^{r,(l-1)}\right)\\ s.t.:V\le V^{U,l}\\ V\ge V^{L,l}\\ \underset{s_i\in U_k^r}{\max }\left(\frac{{\sigma }_{\mathrm{vm}}^{s_i,r}}{g^w\left({\rho }_{s_i}\right){\sigma }_Y^{s_i}}\right)\le {\theta }_{r,k}^{l,1}(r=\mathrm{1,2},\Lambda ,L;k=\mathrm{1,2},\Lambda ,N_c)\\ f_{\mathrm{sc}}^r={\left(\frac{1}{P+R}\underset{i=1}{\overset{P+R}{\Sigma }}{\left(\frac{{\sigma }_{\mathrm{vm}}^{s_i,r}}{g^w\left({\rho }_{s_i}\right){\sigma }_Y^{s_i}}\right)}^{q_2}\right)}^{1/q_2}\\ {\bar{\rho }}_{n_v}^1\le {\rho }_{n_v}\le {\bar{\rho }}_{n_v}=1,v=\mathrm{1,2},\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\approx \underset{v=1}{\overset{Q}{\Sigma }}{\rho }_{n_v}{V}_0^{n_v}+\underset{p=1}{\overset{P}{\Sigma }}{\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；${\rho }_a^m=f_k^{-1}(\underset{j=1}{\overset{m^h}{\Sigma }}{f}_k\left({\rho }_{m_j}^{\left(l\right)}\right)/m^h)---\left(15\right)$$N_a=\left\{m_j\right|{\rho }_{m_j}^{\left(l\right)}\ge \max ({\rho }_a^m,{\bar{\rho }}_m^1)\mathrm{and}({\rho }_{m_j}^{\left(l\right)}-{\rho }_{m_j}^{(l-1)})\ge 0.0\}---\left(16\right)$式中m^h为当前迭代步初始有效结构上人工材料单元个数，m_j为当前迭代步初始有效结构上人工材料单元编号，为第五步获得的当前迭代步的第m_j单元拓扑变量，是函数的逆函数；设计空间的减缩和扩展操作公式表示为式(17)：${\rho }_{h_i}^{*}=\left\{\begin{array}{cc}{\rho }_{h_i}^{\left(l\right)}& h_i\notin N_d\\ 0& h_i\in N_d\end{array}\right.,$${\rho }_{m_j}^{*}=\left\{\begin{array}{cc}{\rho }_{m_j}^{\left(l\right)}& m_j\in N_a\\ 0& m_j\notin 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)确定单元集合${\bar{N}}_d^T=\left\{h_i\right|{\rho }_{h_i}^{\left(l\right)}<{\rho }_d^{\mathrm{Tp}}\}---\left(19\right)$式中是一预先指定的阀值，选取该阶段设计空间的减缩和扩展操作公式表示为：${\rho }_{h_i}^{*}=\left\{\begin{array}{cc}1& h_i\notin {\bar{N}}_d^T\\ 0& h_i\in {\bar{N}}_d^T\end{array}\right.---\left(20\right)$根据，执行方程(20)的操作，并将当前迭代步有效结构所有人工材料单元的拓扑变量赋零，获得的优化结构作为第二个优化求解阶段的初始结构；将设计空间减缩和扩展得到的结构拓扑变量矢量仍记为ρ^(l)，置迭代步数l＝l+1，转向第二步继续优化迭代。6.3)第二个优化求解阶段的结构设计空间减缩和扩展根据第五步获得的拓扑变量解，按式(21a)和(21b)计算两个参数，而后按式(22)确定单元集合和，形成的单元集合含有较大拓扑变量值的单元；${\rho }_d^{\left(l\right)}=\frac{1}{2}\{f_k^{-1}(\underset{i=1}{\overset{N}{\Sigma }}{f}_k\left({\rho }_i^{\left(l\right)}\right)/N)+\underset{i=1}{\overset{N}{\Sigma }}{\rho }_i^{\left(l\right)}/N\}---\left(21a\right)$${\rho }_a^{\left(l\right)}=\frac{1}{2}\{f_k^{-1}(\underset{j=1}{\overset{P+R}{\Sigma }}{f}_k\left({\rho }_{s_j}^{\left(l\right)}\right)/(P+R))+\underset{i=1}{\overset{N}{\Sigma }}{\rho }_i^{\left(l\right)}/N\}---\left(21b\right)$${\bar{N}}_d=\left\{h_i\right|{\rho }_{h_i}^{\left(l\right)}\le {\rho }_d^{\left(l\right)}\},$${\bar{N}}_a=\left\{m_j\right|{\rho }_{m_j}^{\left(l\right)}\ge {\rho }_a^{\left(l\right)}\}---\left(22\right)$式中s_j(j＝1，2，Λ，(P+R))是当前迭代步初始有效结构的所有保留单元编号，P+R为当前迭代步初始有效结构所有保留单元个数；为第五步获得的当前迭代步的第i单元拓扑变量；是函数的逆函数；该阶段设计空间的减缩和扩展操作公式表示为式(23)：${\rho }_{h_i}^{*}=\left\{\begin{array}{cc}1& h_i\notin {\bar{N}}_d\\ 0& h_i\in {\bar{N}}_d\end{array}\right.,$${\rho }_{m_j}^{*}=\left\{\begin{array}{cc}1& m_j\in {\bar{N}}_a\\ 0& m_j\notin {\bar{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},\Lambda ,L}{\max }\{1-{\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，转向第二步继续优化迭代。