EI、Scopus 收录
中文核心期刊

基于自适应泡泡法的薄壳结构拓扑优化设计

张华麟, 杨东, 史之君, 蔡守宇

张华麟, 杨东, 史之君, 蔡守宇. 基于自适应泡泡法的薄壳结构拓扑优化设计. 力学学报, 2023, 55(5): 1165-1173. DOI: 10.6052/0459-1879-22-562
引用本文: 张华麟, 杨东, 史之君, 蔡守宇. 基于自适应泡泡法的薄壳结构拓扑优化设计. 力学学报, 2023, 55(5): 1165-1173. DOI: 10.6052/0459-1879-22-562
Zhang Hualin, Yang Dong, Shi Zhijun, Cai Shouyu. Topology optimization of thin shell structures based on adaptive bubble method. Chinese Journal of Theoretical and Applied Mechanics, 2023, 55(5): 1165-1173. DOI: 10.6052/0459-1879-22-562
Citation: Zhang Hualin, Yang Dong, Shi Zhijun, Cai Shouyu. Topology optimization of thin shell structures based on adaptive bubble method. Chinese Journal of Theoretical and Applied Mechanics, 2023, 55(5): 1165-1173. DOI: 10.6052/0459-1879-22-562
张华麟, 杨东, 史之君, 蔡守宇. 基于自适应泡泡法的薄壳结构拓扑优化设计. 力学学报, 2023, 55(5): 1165-1173. CSTR: 32045.14.0459-1879-22-562
引用本文: 张华麟, 杨东, 史之君, 蔡守宇. 基于自适应泡泡法的薄壳结构拓扑优化设计. 力学学报, 2023, 55(5): 1165-1173. CSTR: 32045.14.0459-1879-22-562
Zhang Hualin, Yang Dong, Shi Zhijun, Cai Shouyu. Topology optimization of thin shell structures based on adaptive bubble method. Chinese Journal of Theoretical and Applied Mechanics, 2023, 55(5): 1165-1173. CSTR: 32045.14.0459-1879-22-562
Citation: Zhang Hualin, Yang Dong, Shi Zhijun, Cai Shouyu. Topology optimization of thin shell structures based on adaptive bubble method. Chinese Journal of Theoretical and Applied Mechanics, 2023, 55(5): 1165-1173. CSTR: 32045.14.0459-1879-22-562

基于自适应泡泡法的薄壳结构拓扑优化设计

基金项目: 国家自然科学基金(11702254), 河南省科技攻关(212102210068)和国家重点研发计划(2017YFB1102800)资助项目
详细信息
    通讯作者:

    蔡守宇, 副教授, 主要研究方向为结构拓扑优化设计. E-mail: caishouyu@zzu.edu.cn

  • 中图分类号: O343.2

TOPOLOGY OPTIMIZATION OF THIN SHELL STRUCTURES BASED ON ADAPTIVE BUBBLE METHOD

  • 摘要: 为有效解决薄壳结构拓扑优化设计难题, 并满足其对分析模型精度和优化结果质量的高要求, 结合等几何壳体分析方法提出一种基于自适应泡泡法的新型拓扑优化设计框架. 等几何分析技术在薄壳分析方面具有天然的优势: 一方面可为薄壳结构建立起精确的NURBS分析模型, 避免了模型转换操作及误差; 另一方面还可保证待分析物理场的高阶连续性, 无需设置转角自由度等. 为了在给定壳面上实现结构的拓扑演化, 借助NURBS曲面(即等几何分析中的薄壳中面)的映射关系, 仅需在规则的二维参数区域内改变结构拓扑即可. 鉴于此, 采用自适应泡泡法在壳面参数区域内开展拓扑优化, 该方法包含孔洞建模、孔洞引入和固定网格分析3个模块, 其在当前工作中分别基于闭合B样条、拓扑导数理论和有限胞元法实现. 其中, 闭合B样条兼具参数和隐式两种表达形式, 参数形式便于在CAD系统中直接生成精确的结构模型; 隐式形式不仅便于开展孔洞的融合/分离操作, 还能与有限胞元法有机结合以替代繁琐的修剪曲面分析方法. 理论分析和数值算例表明, 所提优化设计框架将复杂的薄壳结构拓扑优化问题转化为简单的二维结构拓扑优化问题, 在保证足够分析精度的基础上使用相对很少的设计变量就可得到具有清晰光滑边界且便于导入到CAD系统的优化结果.
    Abstract: Combined with the isogeometric shell analysis (IGA) method, a new optimization design framework based on the adaptive bubble method (ABM) is proposed in this work, in order to effectively solve the topology optimization design problem of thin shell structures, and meet the high standards for the accuracy of analysis models as well as the quality of optimization results. The IGA technique has its natural advantages in thin shell analysis: on one hand, precise analysis models for thin shell structures are established with IGA, and model transformation operations and the resulting errors could therefore be avoided; on the other hand, the high-order continuity of physical fields to be solved can be guaranteed without setting the rotational degrees of freedom. Thanks to the mapping relationship related to the NURBS surface (i.e., the middle surface of thin shell), the structural topology evolution of a given shell surface can be achieved easily in the 2D regular parametric domain. In view of this, the ABM is adopted to carry out topology optimization in the parametric domain, and it contains three modules: the modeling of holes with closed B-splines (CBS), the insertion of holes via the topological derivative theory, and the fixed-grid analysis based on the finite cell method (FCM). It should be noted that holes are expressed in both parametric and implicit forms with CBS. The parametric form makes it convenient to import the structural model into the CAD system exactly. The implicit form not only facilitates the merging and separating operations of holes, but also can be well combined with the FCM which is far more convenient than trimming surface analysis (TSA). Theoretical analysis and numerical examples indicate that the proposed design framework could convert the complicated thin shell structural topology optimization problem into the simple one in the 2D domain, and optimized results with clear and smooth boundaries can be obtained with relatively few design variables.
  • 薄壳结构在自然界中随处可见, 如起到保护作用的贝壳、鸡蛋壳和花生壳等. 由于其具有比强度/刚度高、空间承载性能好等优点, 现已广泛应用于航空、航天、航海、机械、交通和建筑等工程实际装备和产品中. 在航空航天工业快速发展、装备制造业同质化竞争加剧以及大力推进资源节约集约利用的今天, 自重小、跨度大、用料少和形式多的薄壳结构成为当前优化技术的重要设计对象[1-5].

    作为材料节省效果最为显著的致力于充分挖掘材料性能的优化技术, 拓扑优化自20世纪80年代Bendsøe等[6]提出基于均匀化理论的实现方法以来, 进入了快速发展时期并形成了各具特色的方法[7-10]. 多类拓扑优化方法已在壳结构轻量化设计等方面得到应用[11-22], 其中变密度法由于物理概念清晰简洁且已集成到各大商用CAE软件, 而在壳结构拓扑优化设计中被普遍采用[4, 11-15]. 然而, 相关优化工作大都使用传统有限元分析方法对壳结构进行力学响应计算, 结构分析模型(CAE模型)与几何模型(CAD模型)由于采用了不同的描述方式而存在不一致性: 前者由有限个离散的单元构成且常使用Lagrange, Serendipity或Hermite等插值函数; 后者一般由连续解析的样条语言(如NURBS) 描述. 鉴于壳结构(尤其是薄壳结构)的分析结果对模型误差异常敏感, 在基于有限元分析方法的优化设计中不仅要采用十分精细的单元划分, 还要进行繁琐的模型转换操作.

    等几何分析(isogeometric analysis, IGA)方法实现了CAD与CAE模型的统一[23], 然而其常用的张量积形式的NURBS模型具有极大的拓扑局限性, 难以描述拓扑优化中必然出现的多孔结构. 为克服IGA的不足并实现等几何拓扑优化设计, Seo等[24]结合裁剪曲面分析(trimmed surface analysis, TSA)技术对被引入多个孔洞的拓扑复杂壳结构进行高精度的力学响应分析. Kang等[25]进一步开展了基于TSA的壳结构拓扑优化设计. Zhang等[26-27]将移动可变形孔洞(moving morphable void, MMV)方法与TSA技术相结合, 解决了孔洞边界变形异常(如Zigzag边界和自相交现象)等多个问题, 形成了具有较好数值稳定性的壳结构拓扑优化方法.

    前述基于TSA的等几何拓扑优化方法在壳结构设计方面的技术优势除了设计变量少、分析精度高以及优化结果质量好之外, 还包括: 在简单二维参数区域内即可实现对复杂壳体的拓扑优化设计; 无需设置转角自由度就能保证待分析物理场的高阶连续性; 使用单个NURBS面片即可构建拓扑复杂的结构模型; 避免了模型转换操作及由之带来的模型误差. 但是, 此类方法存在以下问题.

    (1) TSA的两个核心步骤较为繁琐, 不利于动辄需要迭代几百甚至上千步的拓扑优化设计的开展. 首先是基于传统B样条修剪曲线[24-25]或凸多边形B样条曲线[26-27]的裁剪单元识别步骤, 远难于基于隐式曲线的判别过程; 其次是关于裁剪单元的积分运算步骤, 需要对含曲线边界的三角形细分单元使用复杂的NURBS增强积分策略.

    (2) 以传统或凸多边形B样条曲线为边界的孔洞在融合与分离时, 需要为新的孔洞构建B样条曲线边界, 这种复杂操作在拓扑优化中会反复进行.

    为解决上述两个问题, 并保留基于TSA的等几何拓扑优化方法的众多技术优势, 本工作在利用IGA进行薄壳结构分析的基础上, 采用自适应泡泡法(adaptive bubble method, ABM)[28]在壳面的二维参数区域内开展拓扑优化. ABM的孔洞建模模块使用兼具参数和隐式两种表达形式的闭合B样条(closed B-splines, CBS)[29]. 参数形式使优化结果能够方便精确地转换成CAD模型; 隐式形式不仅便于裁剪单元的识别以及孔洞的融合与分离操作, 还能配合ABM的固定网格分析模块——有限胞元法(finite cell method, FCM)[30]对裁剪单元进行高效率、高精度的自适应数值积分. 此外, 得益于ABM的孔洞引入模块, 当前设计框架不存在初始布局依赖性问题.

    IGA通常采用NURBS曲面对薄壳结构进行精确建模[23]. NURBS是一种非常优秀的建模方式, 国际标准化组织在1991年正式颁布的工业产品数据交换标准中, 将之定义为自由型曲线曲面的唯一数学描述方法. NURBS曲面的表达式为

    $$ S(\xi ,\eta ) = \sum\limits_{i = 1}^n {\sum\limits_{j = 1}^m {{R_{i,j}}(\xi ,\eta ){{\boldsymbol{P}}_{i,j}}} } $$ (1)

    式中, Pi,j为曲面控制点, Ri,j(ξ, η)为定义如下的二维NURBS基函数

    $$ {R_{i,j}}(\xi ,\eta ) = \frac{{{N_{i,p}}(\xi ){N_{j,q}}(\eta ){w_{i,j}}}}{{\displaystyle\sum\limits_{i = 1}^n {\displaystyle\sum\limits_{j = 1}^m {{N_{i,p}}(\xi ){N_{j,q}}(\eta )} } {w_{i,j}}}}{\text{ , }}0 \leqslant \xi ,\eta \leqslant 1 $$ (2)

    式中, wi, j为权因子, (ξ, η)表示参数域坐标, 参数域由两个一维节点向量Ξ = {ξ1, ξ2, ξ3,···, ξn + p + 1}和Η = {η1, η2, η3,···, ηm+q+1}组成, Ni, p(ξ)和Nj, q(η)分别为定义在ΞΗ上的p次和q次B样条基函数, 它们的个数分别为nm, 且可由Cox-deBoor递推公式计算得到[31].

    由于三维物理域内的NURBS曲面是由其二维矩形参数域自然映射而成, 则对于带孔薄壳结构, 一种建模思想如图1所示. 首先构建出表达薄壳中面的NURBS模型, 进而利用NURBS曲面的映射关系, 在二维参数域中开孔即可实现三维物理域中的薄壳结构开孔. 在此基础上, 通过调节参数域中的孔洞形状即可改变带孔薄壳曲面上的孔洞边界, 而且孔洞边界在变化过程中始终位于薄壳曲面上. 如果进一步改变参数域中的孔洞数量, 则实现了薄壳结构的拓扑变化. 基于该思想, 本文将复杂的三维薄壳结构拓扑优化问题转化为简单的二维参数域内拓扑优化问题.

    图  1  带孔薄壳结构建模示意图
    Figure  1.  Schematic diagram of modeling an opening thin shell

    本文采用Kirchhoff-Love壳的IGA理论[32]对薄壳结构进行力学响应分析, 薄壳变形由中面的变形描述. 中面上的点位移可定义为

    $$ {\boldsymbol{u}}({\theta ^1},{\theta ^2}) = {\boldsymbol{x}}({\theta ^1},{\theta ^2}) - \bar {\boldsymbol{x}}({\theta ^1},{\theta ^2}) $$ (3)

    式中, $\bar {\boldsymbol{x}}$x分别是中面上任意点在参考构型和实际构型的位置向量, θ1θ2是曲线坐标. 由于忽略了横向剪切应变, 薄壳的格林−拉格朗日应变张量可表示为

    $$ {E_{\alpha \beta }} = {\varepsilon _{\alpha \beta }} + {\theta ^3}{\kappa _{\alpha \beta }} $$ (4)

    式中, αβ从1和2中取值, θ3为厚度方向上的坐标, 常数部分εαβ表示薄膜应变, 线性部分καβ表示弯曲应变[33].

    对薄壳结构的总势能求变分可得

    $$ \delta \varPi = \delta {W_{{{\rm{int}}} }} + \delta {W_{{\rm{ext}}}} = 0 $$ (5)
    $$\begin{split} & \delta {W_{{{\rm{int}}} }} = \mathop \int \nolimits_\varOmega \left( {{\boldsymbol{n}}:\delta {\boldsymbol{\varepsilon}} + {\boldsymbol{m}}:\delta {\boldsymbol{\kappa}} } \right){\rm{d}}\varOmega = \\ &\qquad \mathop \int \nolimits_\varOmega \left( {t\delta {{\boldsymbol{\varepsilon}} ^{\rm{T}}}{\boldsymbol{C}}{\boldsymbol{\varepsilon}} + \frac{{{t^3}}}{{12}}\delta {{\boldsymbol{\kappa}} ^{\rm{T}}}{\boldsymbol{C}}{\boldsymbol{\kappa}} } \right){\rm{d}}\varOmega \end{split} $$ (6)
    $$ \delta {W_{{\rm{ext}}}} = \int_\varOmega {({{\boldsymbol{f}}_q} \cdot \delta {\boldsymbol{u}}){\rm{d}}\varOmega } + \int_\varGamma {({{\boldsymbol{f}}_N} \cdot \delta {\boldsymbol{u}}){\rm{d}}\varGamma } $$ (7)

    式中, WintWext分别为内力功和外力功, nm分别为薄膜内力张量和弯曲内力张量, C是材料的本构张量, t是薄壳厚度, fq为单位面积的均布载荷, fN为单位长度的轴向力, u是式(3)定义的薄壳中面上的点位移, εκ是与式(4)中εαβκαβ分别相对应的薄膜应变张量和弯曲应变张量.

    利用式(2)定义的NURBS基函数对薄壳结构的位移场进行离散, 并简化排序规则, 将Ri,j记为RI, 则任意点的近似位移可表示为

    $$ {{\boldsymbol{u}}^h} = \left( {\begin{array}{*{20}{c}} u \\ v \\ w \end{array}} \right) = \sum\limits_{I = 1}^{np} {{R_I}(\xi ,\eta )} {{\boldsymbol{u}}_I} $$ (8)

    式中, uI (I = 1, 2, ···, np)是控制点的位移, np是控制点的个数. 该点的薄膜应变和弯曲应变分别为

    $$ {{\boldsymbol{\varepsilon}} ^h} = \sum\limits_{I = 1}^{np} {{{\boldsymbol{M}}_I}\left( {\xi ,\eta } \right)} {{\boldsymbol{u}}_I},\quad {{\boldsymbol{\kappa}} ^h} = \sum\limits_{I = 1}^{np} {{{\boldsymbol{B}}_I}\left( {\xi ,\eta } \right){{\boldsymbol{u}}_I}} $$ (9)

    式中, MIBI分别为与第I个控制点相对应的薄膜应变矩阵M和弯曲应变矩阵B的分块子矩阵. 将式(8)和式(9)代入式(5) ~ 式(7)得到平衡方程

    $$ {\boldsymbol{KU}}{\text{ = }}{\boldsymbol{F}} $$ (10)

    其中, 薄壳结构刚度矩阵和载荷向量为

    $$ {\boldsymbol{K}}{\text{ = }}\int_\varOmega {\left( {{{\boldsymbol{M}}^{\rm{T}}}{\boldsymbol{DM}}t + \frac{{{t^3}}}{{12}}{{\boldsymbol{B}}^{\rm{T}}}{\boldsymbol{DB}}} \right)} {\rm{d}}\varOmega $$ (11)
    $$ {\boldsymbol{F}}{\text{ = }}\int_\varOmega {{{\boldsymbol{R}}^{\text{T}}}{{\boldsymbol{f}}_q}} {\rm{d}}\varOmega + \int_\varGamma {{{\boldsymbol{R}}^{\text{T}}}{{\boldsymbol{f}}_N}} {\rm{d}}\varGamma $$ (12)

    式中, D是材料的弹性矩阵, R是由NURBS基函数RI组建而成的形函数矩阵.

    在薄壳曲面参数域内开展拓扑优化的自适应泡泡法包含3个模块: 孔洞建模、孔洞引入和固定网格分析, 下面将逐一进行介绍.

    作为光滑变形隐式曲线(smoothly deformable implicit curves, SDIC)[28]的简化形式, CBS在参数域中的隐式表达式为

    $$ {\xi ^2} + {\eta ^2} = R{(\theta )^2} $$ (13)

    其参数形式为

    $$ \left.\begin{split} & \xi = R(\theta )\cos \theta \\ & \eta = R(\theta )\sin \theta \end{split} \right\} $$ (14)

    式中, R(θ)为半轴长函数, θ为定义如下的极角

    $$ \theta = \arctan \frac{\eta }{\xi } + H( - \xi )\text{π} $$ (15)

    其中, H(·)为Heaviside函数, 可见θ的值域为[−π/2, 3π/2). 为使R(θ)能够随θ的改变而自由光滑地变化, 借助B样条基函数Nk,r(ζ)将其定义为

    $$ R\left( \theta \right) = \sum\limits_{k = 1}^{nb} {{N_{k,r}}(\zeta )} {p_k} = \sum\limits_{k = 1}^{nb} {{N_{k,r}}\left( {\frac{{\theta + \text{π} /2}}{{2\text{π} }}} \right){p_k}} $$ (16)

    式中, nbr分别是B样条基函数的个数和次数, pk是控制参数, ζ的取值范围为[0, 1). 本工作使用节点均匀分布的unclamped类型B样条曲线表示R(θ), 见图2(a), 这样每个控制参数pk的作用可看作是等价的. 为形成图2(b)所示的CBS, 需使前r个和后r个控制点相重合, 也即pλ = pnb−r+λ (λ = 1, 2, ···, r). 采用该方式构建的CBS的光滑性处处相同.

    图  2  基于CBS的孔洞表示方式
    Figure  2.  CBS-based hole representation

    在拓扑优化中孔洞是被逐步引入的, 根据式(13), 可将第i个引入的孔洞用隐函数表示如下

    $$ {\varphi _i} = \sqrt {{{\left( {\xi - {\xi _i}} \right)}^2} + {{\left( {\eta - {\eta _i}} \right)}^2}} - {R_i}\left( {{\theta _i}} \right) $$ (17)

    式中, (ξi, ηi)T为该孔洞的插入点坐标. 参见图1, 若将参数域中Dp的隐函数记为ΦDp, 则Ωp的隐函数为

    $$ {\varPhi _{\varOmega p}} = \min \left\{ {{\varPhi _{Dp}},\left\{ {{\varphi _i}} \right\}_{i = 1}^{num}} \right\} $$ (18)

    式中, num表示孔的数量, 其值在优化迭代中会动态增加. ΦΩp满足

    $$ \left.\begin{split} & {\varPhi _{\varOmega p}}\left( {\boldsymbol{\xi}} \right) = 0,\forall {\boldsymbol{\xi}} \in \partial {\varOmega _p} \\ & {\varPhi _{\varOmega p}}\left( {\boldsymbol{\xi}} \right) > 0,\forall {\boldsymbol{\xi}} \in {\varOmega _p}{{\backslash }}\partial {\varOmega _p} \\ & {\varPhi _{\varOmega p}}\left( {\boldsymbol{\xi}} \right) < 0,\forall {\boldsymbol{\xi}} \in {D_p}{{\backslash }}{\varOmega _p} \end{split} \right\} $$ (19)

    式中, ξ表示参数域中的任意点, 根据隐函数ΦΩp可直接判断出点ξ是否位于区域Ωp内.

    拓扑导数反映无限小单一区域的扰动(引入孔洞、夹杂和裂纹等)对给定函数的影响. 当前优化工作的目标函数为薄壳结构柔顺度

    $$\begin{split} & C = \mathop \int \nolimits_\varOmega \left( {{\boldsymbol{n}}\left( {\boldsymbol{u}} \right):{\boldsymbol{\varepsilon}} \left( {\boldsymbol{u}} \right) + {\boldsymbol{m}}\left( {\boldsymbol{u}} \right):{\boldsymbol{\kappa}} \left( {\boldsymbol{u}} \right)} \right){\rm{d}}\varOmega= \\ &\qquad \mathop \int \nolimits_{{\varOmega _p}} \left( {{\boldsymbol{n}}\left( {\boldsymbol{u}} \right):{\boldsymbol{\varepsilon}} \left( {\boldsymbol{u}} \right) + {\boldsymbol{m}}\left( {\boldsymbol{u}} \right):{\boldsymbol{\kappa}} \left( {\boldsymbol{u}} \right)} \right)\left| {\boldsymbol{J}} \right|{\rm{d}}\xi {\rm{d}}\eta \end{split} $$ (20)

    式中, J为反映薄壳中面参数域与物理域之间映射关系的雅可比矩阵, 其余参数定义同前. 参考文献[34-35], 若不考虑体积力等设计相关性载荷, 则上式表达的柔顺度的拓扑导数为

    $$ \begin{split} &{D}_{{\rm{T}}}\left({{\boldsymbol{\xi}} }_{\varOmega p}\right) = \underset{{r}_{t}\to 0}{\mathrm{lim}}\frac{C\left({\varOmega }_{p}\backslash B({{\boldsymbol{\xi}} }_{\varOmega p},{r}_{t})\right)-C({\varOmega }_{p})}{\text{π} {r}_{t}{}^{2}}\text{ = }\\ &\qquad \frac{4\left|{\boldsymbol{J}}\right|}{1 + v}\left[{\boldsymbol{n}}\left({\boldsymbol{u}}\left({{\boldsymbol{\xi}} }_{\varOmega p}\right)\right):{\boldsymbol{\varepsilon}} \left({\boldsymbol{u}}\left({{\boldsymbol{\xi}} }_{\varOmega p}\right)\right)\right] + \\ &\qquad \frac{4\left|{\boldsymbol{J}}\right|}{1 + v}\left[{\boldsymbol{m}}\left({\boldsymbol{u}}\left({{\boldsymbol{\xi}} }_{\varOmega p}\right)\right):{\boldsymbol{\kappa}} \left({\boldsymbol{u}}\left({{\boldsymbol{\xi}} }_{\varOmega p}\right)\right)\right]-\\ &\qquad \frac{1-3v}{1-{v}^{2}}\left|{\boldsymbol{J}}\right|{\rm{tr}}\left[{\boldsymbol{n}}\left({\boldsymbol{u}}\left({{\boldsymbol{\xi}} }_{\varOmega p}\right)\right)\right]{\rm{tr}}\left[{\boldsymbol{\varepsilon}} \left({\boldsymbol{u}}\left({{\boldsymbol{\xi }}}_{\varOmega p}\right)\right)\right]-\\ &\qquad \frac{1-3v}{1-{v}^{2}}\left|{\boldsymbol{J}}\right|{\rm{tr}}\left[{\boldsymbol{m}}\left({\boldsymbol{u}}\left({\xi }_{\varOmega p}\right)\right)\right]{\rm{tr}}\left[{\boldsymbol{\kappa}} \left({\boldsymbol{u}}\left({{\boldsymbol{\xi}} }_{\varOmega p}\right)\right)\right]\end{split}$$ (21)

    式中, B(ξΩp, rt)表示区域Ωp内以任意点ξΩp为圆心且以rt为半径的圆孔, v为材料泊松比. 由于孔洞的引入会削弱结构刚度, 也即导致柔顺度的增加, 因此拓扑导数DT是非负的. 若某点处DT值很小或接近于零, 则表示在该点引入孔洞对结构柔顺度的影响可忽略不计, 因而可优先在此处插入孔洞.

    参照文献[28]的早期研究, 为保证拓扑优化过程的稳健性, 设置拓扑导数阈值${{\overline {D_{\rm{T}}}}}$以限制每步优化迭代的开孔数量; 为使设计变量尽量被充分利用, 对已插入的孔洞设置影响区域Ii以避免在其附近引入新孔. 图3展示了单步优化迭代的孔洞引入过程, 首先对Ωp内采样点的DT值从小到大进行排序, 并将DT小于${{\overline {D_{\rm{T}}}}}$的采样点作为潜在插入点(图中标序号的点), 进而结合影响区域逐步引入孔洞.

    图  3  孔洞自适应引入机制示意图
    Figure  3.  Adaptive introduction mechanism of holes

    FCM本质上一种是采用高阶形函数插值逼近待求物理场的虚拟区域法[30], 其结合式(18)定义的隐函数ΦΩp开展力学响应分析的方式见图4. 参数域中的Ωp由规则设计区域Dp开孔而得, 则根据FCM的思想, 可将Dp作为计算域并用矩形胞元离散, 进而结合ΦΩp将胞元快速区分为虚拟胞元、物理胞元和边界胞元, 以分类进行计算. 其中, 对于被∂Ωp切割的边界胞元需采用四叉树自适应积分策略, 如图4(c)所示, 值得一提的是基于ΦΩp的四叉树细化操作十分简便. 由式(11)和式(12)得到FCM的刚度矩阵和载荷向量如下

    图  4  有限胞元方法示意图
    Figure  4.  Schematic diagram of the finite cell method
    $$ {\boldsymbol{K}}{\text{ = }}\int_{{D_p}} {\left( {{{\boldsymbol{M}}^{\rm{T}}}{\boldsymbol{DM}}t + \frac{{{t^3}}}{{12}}{{\boldsymbol{B}}^{\rm{T}}}{\boldsymbol{DB}}} \right)H\left( {{\varPhi _{\varOmega p}}} \right)} \left| {\boldsymbol{J}} \right|{\rm{d}}\xi {\rm{d}}\eta $$ (22)
    $$ {\boldsymbol{F}}{\text{ = }}\int_{{D_p}} {{{\boldsymbol{R}}^{\text{T}}}{{\boldsymbol{f}}_q}} H\left( {{\varPhi _{\varOmega p}}} \right)\left| {\boldsymbol{J}} \right|{\rm{d}}\xi {\rm{d}}\eta + \int_{{\varGamma _{Np}}} {{{\boldsymbol{R}}^{\text{T}}}{{\boldsymbol{f}}_N}\left| {{{\boldsymbol{J}}_\varGamma }} \right|} {\rm{d}}{\varGamma _p} $$ (23)

    其中, JΓ为反映参数域中Ωp边界∂Ωp和物理域中薄壳结构边界∂Ω之间映射关系的雅可比矩阵. 由式(19)可知, 当点ξ位于Ωp内部时H(ФΩp(ξ))为1, 反之则H(ФΩp(ξ))为0.

    本文旨在利用有限量材料实现薄壳结构的刚度最大化(柔顺度最小化)设计, 这种典型拓扑优化问题的数学模型为

    $$ \begin{split} &\text{find }{\left\{{\xi }_{i},{\eta }_{i},{\left\{{p}_{i,k}\right\}}_{k = 1}^{nb\text-r}\right\}}_{i = 1}^{num}\\ &\text{min }C\text{ = }{{\boldsymbol{F}}}^{\text{T}}{\boldsymbol{U}}\\ &\text{s}\text{.t}\text{. }{\boldsymbol{KU}}\text{ } = \text{ }{\boldsymbol{F}}\\ &V = \text{ }{\displaystyle {\int }_{{D}_{p}}H\left({\varPhi }_{\varOmega p}\right)\left|{\boldsymbol{J}}\right|{\rm{d}}\xi {\rm{d}}\eta \leqslant {V}_{\mathrm{lim}}},\\ &\qquad{p}_{i,k}\text{ }>\text{0, }i = 1,2,\cdots ,num,\text{ }k\text{ = 1, 2,}\cdots ,nb-r\end{split} $$ (24)

    式中, V为薄壳结构体积, Vlim为体积上限, 其余参数定义同前, 每个孔洞有nbr + 2个设计变量, 现将变量统一记为

    $$ {\boldsymbol{v}} = \left\{ {{v_s}} \right\}_{s = 1}^{num \times (nb - r + 2)} = \left\{ {{\xi _i},{\eta _i},\left\{ {{p_{i,k}}} \right\}_{k = 1}^{nb - r}} \right\}_{i = 1}^{num} $$ (25)

    薄壳结构柔顺度C对设计变量vs的灵敏度为

    $$ \frac{{\partial C}}{{\partial {v_s}}} = {\left( {\frac{{\partial {\boldsymbol{F}}}}{{\partial {v_s}}}} \right)^{\rm{T}}}{\boldsymbol{U}} + {{\boldsymbol{F}}^{\rm{T}}}\frac{{\partial {\boldsymbol{U}}}}{{\partial {v_s}}} $$ (26)

    F是设计无关载荷, 即上式右端第一项为0, 则结合平衡方程KU = F可得

    $$ \frac{{\partial C}}{{\partial {v_s}}} = {{\boldsymbol{F}}^{\rm{T}}}\frac{{\partial {\boldsymbol{U}}}}{{\partial {v_s}}}{\text{ = }}{{\boldsymbol{F}}^{\rm{T}}}\left( { - {{\boldsymbol{K}}^{ - 1}}\frac{{\partial {\boldsymbol{K}}}}{{\partial {v_s}}}{\boldsymbol{U}}} \right) = - {{\boldsymbol{U}}^{\rm{T}}}\frac{{\partial {\boldsymbol{K}}}}{{\partial {v_s}}}{\boldsymbol{U}} $$ (27)

    为使推导过程更加通用, 将结构刚度和体积函数统一表达为

    $$ G = \int_{{D_p}} {gH({\varPhi _{\varOmega p}})\left| {\boldsymbol{J}} \right|{\rm{d}}\xi {\rm{d}}\eta } $$ (28)

    g = BTDBt3/12 + MTDMtG表示刚度矩阵; 当g = 1时G表示体积. G的灵敏度推导如下

    $$ \begin{split} & \frac{{\partial G}}{{\partial {v_s}}} = \int_{{D_p}} {g\frac{{\partial H\left( {{\varPhi _{\varOmega p}}} \right)}}{{\partial {v_s}}}} \left| {\boldsymbol{J}} \right|{\rm{d}}\xi {\rm{d}}\eta = \\ &\qquad \int_{{D_p}} {g\frac{{{\rm{d}}H\left( {{\varPhi _{\varOmega p}}} \right)}}{{{\rm{d}}{\varPhi _{\varOmega p}}}}\frac{{\partial {\varPhi _{\varOmega p}}}}{{\partial {v_s}}}} \left| {\boldsymbol{J}} \right|{\rm{d}}\xi {\rm{d}}\eta = \\ &\qquad \int_{{D_p}} {g\frac{{\partial {\varPhi _{\varOmega p}}}}{{\partial {v_s}}}\frac{{\left| {\boldsymbol{J}} \right|}}{{\left\| {\nabla {\varPhi _{\varOmega p}}} \right\|}}\left( {\frac{{{\rm{d}}H\left( {{\varPhi _{\varOmega p}}} \right)}}{{{\rm{d}}{\varPhi _{\varOmega p}}}}\left\| {\nabla {\varPhi _{\varOmega p}}} \right\|} \right){\rm{d}}\xi {\rm{d}}\eta } = \\ &\qquad \int_{{D_p}} {g\frac{{\partial {\varPhi _{\varOmega p}}}}{{\partial {v_s}}}\frac{{\left| {\boldsymbol{J}} \right|}}{{\left\| {\nabla {\varPhi _{\varOmega p}}} \right\|}}\hat \delta {\rm{d}}\xi {\rm{d}}\eta }= \\ &\qquad \int_{\partial {\varOmega _p}} {g\frac{{\partial {\varPhi _{\varOmega p}}}}{{\partial {v_s}}}\frac{{\left| {\boldsymbol{J}} \right|}}{{\left\| {\nabla {\varPhi _{\varOmega p}}} \right\|}}{\rm{d}}} {\varGamma _p} \\[-12pt]\end{split} $$ (29)

    式中, $ \hat \delta $为Heaviside函数的法向导数, 称为Dirac delta函数, 它能将区域积分转换为边界积分[36], 关于∂ΦΩp/∂vs$\nabla $ΦΩp的计算方式见文献[28].

    现使用新型设计框架对典型薄壳结构(壳中面为球面的一部分)进行拓扑优化设计, 并研究壳面曲率对优化结果拓扑构型的影响. 材料的弹性模量和泊松比分别取2.0 × 1011 Pa和0.3, 每个CBS孔洞使用22个设计变量, 体积上限Vlim设为初始体积的40%, 优化算法采用移动渐近线优化算法(method of moving asymptotes, MMA)[37].

    薄壳结构几何模型和边界条件如图5所示, 壳面顶点受40 kN的集中载荷, 4个角点被固定, 顶点相对角点的高度为1.2 m, 壳厚为0.05 m. 对于此对称问题, 本文仅将薄壳结构的1/4作为优化对象.

    图  5  受集中垂直载荷的薄壳结构
    Figure  5.  Thin shell subjected to a vertical load

    对壳面的参数域(即图4中的设计区域或计算域Dp)采用32 × 32的有限胞元网格划分, 网格中点即为计算拓扑导数的采样点, 拓扑导数阈值${{\overline {D_{\rm{T}}}}}$和影响区域Ii分别设为3%和0.25. 图6(a)展示了ABM在参数域的实施过程, 图6(b)为由参数域映射到物理域的壳面拓扑演化过程, 优化过程中共引入10个孔洞, 开孔位置合理且孔洞融合/分离自然, 设计变量共220个. 图7给出了设计结果的位移云图. 目标函数和约束函数的迭代收敛曲线如图8所示, 优化结果的体积严格满足约束. 图9为优化结果的CAD模型, 其可借助CBS参数形式(14)和NURBS的映射操作在CAD系统中生成.

    图  6  薄壳结构拓扑优化过程
    Figure  6.  Topological evolution of the thin shell
    图  7  设计结果的位移云图
    Figure  7.  Displacement contours of the design result
    图  8  柔顺度及体积的收敛曲线
    Figure  8.  The convergent curves for compliance and volume
    图  9  优化结果的CAD模型
    Figure  9.  CAD model of the optimized result

    由于薄壳结构的刚度对壳面曲率的变化非常敏感, 因而不同曲率下的最优拓扑构型可能不同. 为验证这一猜想, 并反映薄壳曲面曲率对拓扑优化结果的影响, 下面将图5壳面的控制点高度坐标值分别缩小2倍和扩大2倍, 并采用本文所提设计框架开展拓扑优化设计. 最终优化结果如图10所示, 薄壳结构的曲面曲率自上而下依次增大, 设计结果的拓扑构型有着较大的差异.

    图  10  不同壳面曲率下的拓扑优化结果对比图(由上至下曲率逐步增大)
    Figure  10.  Comparison of topology optimization design results for different shell curvatures (the curvature gradually increases from top to bottom)

    本文提出了一种新型的薄壳结构拓扑优化设计框架, 该框架采用IGA保证了薄壳结构分析模型的精确性, 并能将三维空间中的薄壳结构优化问题转化到简单二维参数域中进行求解; 采用ABM不仅便于在被孔洞修剪的参数域上开展孔洞融合/分离操作和高精度的固定网格分析计算, 还可消除优化结果的初始布局依赖性. 数值算例表明, 所提设计框架使用少量变量就能得到具有清晰光滑边界且便于导入到CAD系统的优化结果. 进一步的研究将在图10测试结果的基础上, 通过把壳面控制点增添为设计变量, 开展薄壳结构的曲面形状和拓扑构型的协同优化设计, 从而尽可能提升其力学性能.

  • 图  1   带孔薄壳结构建模示意图

    Figure  1.   Schematic diagram of modeling an opening thin shell

    图  2   基于CBS的孔洞表示方式

    Figure  2.   CBS-based hole representation

    图  3   孔洞自适应引入机制示意图

    Figure  3.   Adaptive introduction mechanism of holes

    图  4   有限胞元方法示意图

    Figure  4.   Schematic diagram of the finite cell method

    图  5   受集中垂直载荷的薄壳结构

    Figure  5.   Thin shell subjected to a vertical load

    图  6   薄壳结构拓扑优化过程

    Figure  6.   Topological evolution of the thin shell

    图  7   设计结果的位移云图

    Figure  7.   Displacement contours of the design result

    图  8   柔顺度及体积的收敛曲线

    Figure  8.   The convergent curves for compliance and volume

    图  9   优化结果的CAD模型

    Figure  9.   CAD model of the optimized result

    图  10   不同壳面曲率下的拓扑优化结果对比图(由上至下曲率逐步增大)

    Figure  10.   Comparison of topology optimization design results for different shell curvatures (the curvature gradually increases from top to bottom)

  • [1] 张卫红, 周涵, 李韶英等. 航天高性能薄壁构件的材料−结构一体化设计. 航空学报, 2022, 出版中

    Zhang Weihong, Zhou Han, Li Shaoying, et al. Material-structure integrated design for high-performance aerospace thin-walled component. Acta Aeronautica et Astronautica Sinica, 2022, in press (in Chinese))

    [2] 杨利鑫, 何东泽, 陈强等. 高温环境下大尺度薄壁结构的电性能优化设计. 宇航学报, 2021, 42(9): 1099-1107 (Yang Lixin, He Dongze, Chen Qiang, et al. Electrical performance optimization design of large-scale thin-wall structures in thermal environment. Journal of Astronautics, 2021, 42(9): 1099-1107 (in Chinese)
    [3] 王博, 周子童, 周演等. 薄壁结构多层级并发加筋拓扑优化研究. 计算力学学报, 2021, 38(4): 487-497 (Wang Bo, Zhou Zitong, Zhou Yan, et al. Concurrent topology optimization hierarchical stiffened thin-walled structures. Chinese Journal of Computational Mechanics, 2021, 38(4): 487-497 (in Chinese)
    [4]

    Träff EA, Sigmund O, Aage N. Topology optimization of ultra high resolution shell structures. Thin-Walled Structures, 2021, 160: 107349 doi: 10.1016/j.tws.2020.107349

    [5]

    Jiang BS, Zhang JY, Ohsaki M. Shape optimization of free-form shell structures combining static and dynamic behaviors. Structures, 2021, 29: 1791-1807 doi: 10.1016/j.istruc.2020.12.045

    [6]

    Bendsøe MP, Kikuchi N. Generating optimal topologies in structural design using a homogenization method. Computer Methods in Applied Mechanics and Engineering, 1988, 71(2): 197-224 doi: 10.1016/0045-7825(88)90086-2

    [7]

    Bendsøe MP, Sigmund O. Topology Optimization: Theory, Methods and Applications. Berlin: Springer Verlag, 2003: 9-24

    [8]

    Xie YM, Steven GP. A simple evolutionary procedure for structural optimization. Computers and Structures, 1993, 49(5): 885-896 doi: 10.1016/0045-7949(93)90035-C

    [9] 隋允康, 叶红玲. 连续体结构拓扑优化的 ICM 方法. 北京: 科学出版社, 2013

    Sui Yunkang, Ye Hongling. ICM Method for Topology Optimization of Continuum Structures. Beijing: Science Press, 2013 (in Chinese))

    [10]

    Wang MY, Wang XM, Guo DM. A level set method for structural topology optimization. Computer Methods in Applied Mechanics and Engineering, 2003, 192(1-2): 227-246 doi: 10.1016/S0045-7825(02)00559-5

    [11]

    Maute K, Ramm E. Adaptive topology optimization of shell structures. American Institute of Aeronautics and Astronautics, 1997, 35: 1767-1773 doi: 10.2514/2.25

    [12]

    Ansola R, Canales J, Tarrago JA, et al. On simultaneous shape and material layout optimization of shell structures. Structural and Multidisciplinary Optimization, 2002, 24(3): 175-184 doi: 10.1007/s00158-002-0227-x

    [13] 粟华, 陈伟俊, 龚春林等. 环向肋增稳的薄壁圆筒结构定向拓扑优化方法. 宇航学报, 2022, 43(3): 374-382 (Su Hua, Chen Weijun, Gong Chunlin, et al. Oriented topology optimization of thin-walled structure with circumferential rib enhancement of stability. Journal of Astronautics, 2022, 43(3): 374-382 (in Chinese)
    [14]

    Hassani B, Tavakkoli SM, Ghasemnejad H. Simultaneous shape and topology optimization of shell structures. Structural and Multidisciplinary Optimization, 2013, 48: 221-233

    [15] 牟佳信, 邢彬, 郭梅等. 考虑热、弹、流耦合的航空发动机齿轮箱壳体拓扑优化分析. 机械传动, 2022, 46(2): 127-134 (Mu Jiaxin, Xing Bin, Guo Mei, et al. Topology optimization analysis of aero-engine gearbox housing considering the coupling of thermal elastohydrodynamic lubrication. Journal of Mechanical Transmission, 2022, 46(2): 127-134 (in Chinese)
    [16] 朱润, 隋允康. 基于ICM方法的多工况位移约束下板壳结构拓扑优化设计. 固体力学学报, 2012, 33(1): 81-90 (Zhu Run, Sui Yunkang. Topological optimization design of plate-shell structure under multi-condition displacement constraints based on ICM method. Chinese Journal of Solid Mechanics, 2012, 33(1): 81-90 (in Chinese)
    [17]

    Huo WD, Liu C, Du ZL, et al. Topology optimization on complex surfaces based on the moving morphable component method and computational conformal mapping. Journal of Applied Mechanics, 2022, 89(5): 051008 doi: 10.1115/1.4053727

    [18]

    Xu XQ, Gu XD, Chen SK. Shape and topology optimization of conformal thermal control structures on free-form surfaces: A dimension reduction level set method (DR-LSM). Computer Methods in Applied Mechanics and Engineering, 2022, 398: 115183

    [19]

    Ho-Nguyen-Tan T, Kim HG. An efficient method for shape and topology optimization of shell structures. Structural and Multidisciplinary Optimization, 2022, 65: 119 doi: 10.1007/s00158-022-03213-0

    [20]

    Ho-Nguyen-Tan T, Kim HG. Level set-based topology optimization for compliance and stress minimization of shell structures using trimmed quadrilateral shell meshes. Computers and Structures, 2022, 259: 106695 doi: 10.1016/j.compstruc.2021.106695

    [21]

    Meng XC, Xiong YL, Xie YM, et al. Shape-thickness-topology coupled optimization of free-form shells. Automation in Construction, 2022, 142: 104476 doi: 10.1016/j.autcon.2022.104476

    [22]

    Jiang XD, Zhang WS, Liu C, et al. An explicit approach for simultaneous shape and topology optimization of shell structures. Applied Mathematical Modelling, 2023, 113: 613-639 doi: 10.1016/j.apm.2022.09.028

    [23]

    Hughes TJR, Cottrell JA, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 2005, 194(39-41): 4135-4195 doi: 10.1016/j.cma.2004.10.008

    [24]

    Seo YD, Kim HJ, Youn SK. Isogeometric topology optimization using trimmed spline surfaces. Computer Methods in Applied Mechanics and Engineering, 2010, 199(49-52): 3270-3296 doi: 10.1016/j.cma.2010.06.033

    [25]

    Kang P, Youn SK. Isogeometric topology optimization of shell structures using trimmed NURBS surfaces. Finite Elements in Analysis and Design, 2016, 120: 18-40 doi: 10.1016/j.finel.2016.06.003

    [26]

    Zhang WS, Li DD, Kang P, et al. Explicit topology optimization using IGA-based moving morphable void (MMV) approach. Computer Methods in Applied Mechanics and Engineering, 2019, 360: 112685

    [27]

    Zhang WS, Jiang S, Liu C, et al. Stress-related topology optimization of shell structures using IGA/TSA-based moving morphable void (MMV) approach. Computer Methods in Applied Mechanics and Engineering, 2020, 366: 113036 doi: 10.1016/j.cma.2020.113036

    [28] 蔡守宇, 张卫红, 高彤等. 基于固定网格和拓扑导数的结构拓扑优化自适应泡泡法. 力学学报, 2019, 51(4): 1235-1244 (Cai Shouyu, Zhang Weihong, Gao Tong, et al. Adaptive bubble method using fixed mesh and topological derivative for structural topology optimization. Chinese Journal of Theoretical and Applied Mechanics, 2019, 51(4): 1235-1244 (in Chinese)
    [29]

    Zhang WH, Zhao LY, Gao T, et al. Topology optimization with closed B-splines and Boolean operations. Computer Methods in Applied Mechanics and Engineering, 2017, 315: 652-670 doi: 10.1016/j.cma.2016.11.015

    [30]

    Parvizian J, Düster A, Rank E. Finite cell method. Computational Mechanics, 2007, 41(1): 121-133 doi: 10.1007/s00466-007-0173-y

    [31]

    Piegl L, Tiller W. The NURBS Book. 2nd Edition. New York: Springer Verleg, 1997: 117-139

    [32]

    Kiendl J, Bletzinger KU, Linhard J, et al. Isogeometric shell analysis with Kirchhoff-love elements. Computer Methods in Applied Mechanics and Engineering, 2009, 198(49-52): 3902-3914 doi: 10.1016/j.cma.2009.08.013

    [33]

    Cirak F, Ortiz M, Schröder P. Subdivision surfaces: a new paradigm for thin-shell finite-element analysis. International Journal for Numerical Methods in Engineering, 2000, 47(12): 2039-2072 doi: 10.1002/(SICI)1097-0207(20000430)47:12<2039::AID-NME872>3.0.CO;2-1

    [34]

    Nvotny AA, Sokolowski J. Topological Derivatives in Shape Optimization. Heidelberg: Springer-Verlag, 2013

    [35]

    Sokolowski J, Zochowski A. On the topological derivative in shape optimization. SIAM Journal on Control and Optimization, 1999, 37(4): 1251-1272 doi: 10.1137/S0363012997323230

    [36]

    Cai SY, Zhang WH, Zhu JH, et al. Stress constrained shape and topology optimization with fixed mesh: A B-spline finite cell method combined with level set function. Computer Methods in Applied Mechanics and Engineering, 2014, 278(15): 361-387

    [37]

    Svanberg K. The method of moving asymptotes—a new method for structural optimization. International Journal for Numerical Methods in Engineering, 1987, 24(2): 359-373 doi: 10.1002/nme.1620240207

  • 期刊类型引用(2)

    1. 杨龙,李龙龙. 基于工程限值约束的结构拓扑优化设计分析. 中国建筑金属结构. 2024(02): 39-41 . 百度学术
    2. 于淼,赵文武,单兴业,毕道明. 一种基于应力控制的金属薄壳结构变形修复方法研究. 测控技术. 2024(04): 89-94 . 百度学术

    其他类型引用(1)

图(10)
计量
  • 文章访问数:  835
  • HTML全文浏览量:  298
  • PDF下载量:  169
  • 被引次数: 3
出版历程
  • 收稿日期:  2022-11-27
  • 录用日期:  2023-04-02
  • 网络出版日期:  2023-04-03
  • 刊出日期:  2023-05-17

目录

/

返回文章
返回