留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

基于PDE灵敏度滤波器的算法研究

孟换利,  张岐良,  王杰

文章导航 > 应用数学和力学  > 2023 >  44(1): 80-92
孟换利，张岐良，王杰. 基于PDE灵敏度滤波器的算法研究 [J]. 应用数学和力学，2023，44（1）：80-92 doi: 10.21656/1000-0887.430064
引用本文: 孟换利，张岐良，王杰. 基于PDE灵敏度滤波器的算法研究 [J]. 应用数学和力学，2023，44（1）：80-92 doi: 10.21656/1000-0887.430064
shu
MENG Huanli, ZHANG Qiliang, WANG Jie. Algorithm Research Based on the PDE Sensitivity Filter[J]. Applied Mathematics and Mechanics, 2023, 44(1): 80-92. doi: 10.21656/1000-0887.430064
Citation: MENG Huanli, ZHANG Qiliang, WANG Jie. Algorithm Research Based on the PDE Sensitivity Filter[J]. Applied Mathematics and Mechanics, 2023, 44(1): 80-92. doi: 10.21656/1000-0887.430064
shu

基于PDE灵敏度滤波器的算法研究

doi: 10.21656/1000-0887.430064
  • 1.

    中山大学 系统科学与工程学院，广州 510006

  • 2.

    天津航天长征火箭制造有限公司，天津 300462

基金项目: 中央高校基本科研业务费（76210-31610007）
详细信息
    作者简介:

    孟换利（1995—），女，硕士生（E-mail：menghli@mail2.sysu.edu.cn）

    张岐良（1985—），男，助理教授，博士，硕士生导师（通讯作者. E-mail：zhangqliang@mail.sysu.edu.cn）

  • 中图分类号: O241.82

Algorithm Research Based on the PDE Sensitivity Filter

  • 1.

    School of Systems Science and Engineering, SUN Yat-Sen University, Guangzhou 510006, P.R.China

  • 2.

    Tianjin Aerospace Long March Rocket Manufacturing Co. Ltd. , Tianjin 300462, P.R.China

摘要

    摘要
  • 摘要:

    采用PDE灵敏度滤波器可以消除连续体结构拓扑优化结果存在的棋盘格现象、数值不稳定等问题，且PDE灵敏度滤波器的实质是具有Neumann边界条件的Helmholtz偏微分方程。针对大规模PDE灵敏度滤波器的求解问题，有限元分析得到其代数方程，分别采用共轭梯度算法、多重网格算法和多重网格预处理共轭梯度算法对代数方程进行求解，并且研究精度、过滤半径以及网格数量对拓扑优化效率的影响。结果表明：与共轭梯度算法和多重网格算法相比，多重网格预处理共轭梯度算法迭代次数最少，运行时间最短，极大地提高了拓扑优化效率。

    Abstract:

    The PDE sensitivity filter can eliminate the checkerboard patterns and numerical instability existing in the topology optimization results of continuum structures, and the essence of the PDE sensitivity filter is the Helmholtz partial differential equation with Neumann boundary conditions. To solve the large-scale PDE sensitivity filter problem, the conjugate gradient algorithm, the multigrid algorithm and the multigrid preconditioned conjugate gradient algorithm were used to solve the algebraic equations obtained by finite element analysis, and the effects of accuracy, filter radius and grid numbers on the efficiency of topology optimization were studied. The results show that, compared with the conjugate gradient algorithm and the multi-grid algorithm, the multi-grid preconditioned conjugate gradient algorithm has the least number of iterations and the shortest running time, which greatly improves the efficiency of topology optimization.

    • HTML全文

  • 引　　言

    为应对越来越高的力学性能要求和越来越严酷的服役环境，在结构设计中应用拓扑优化技术已成为一种必然选择。在结构的概念设计阶段通过结构拓扑优化确定材料的最优布局形式，可以帮助设计者发现结构最佳的应力传导路径，从而以最少的材料实现最优的结构性能。Deaton和Grandhi[1]调查了2000年至2012年连续体结构拓扑优化的研究现状与应用。Zhu等[2]综述了面向增材制造的拓扑优化技术的典型应用场景、研究热点和挑战。卫志军等[3]对双隔板进行了拓扑优化，可以有效地抑制流体晃荡。

    然而，连续体结构的拓扑优化结果往往存在棋盘格现象、数值不稳定等问题，Sigmund[4]以及Bruns、Tortorelli[5]使用的灵敏度滤波器和密度滤波器可以有效地解决上述问题。研究者们在此基础上发展了多种过滤方法，如灰度单元等效转换方法[6]、动态参数调整方法[7]和灰度单元分层双重惩罚方法[8]等。上述方法在一定程度上可以有效地抑制灰度单元、消除棋盘格现象，防止边界灰度扩散，但没有关注过滤效率。与通过获取邻域的单元信息来修正中心单元的灵敏度滤波器不同，PDE滤波器的实质是将灵敏度过滤值表示成Neumann边界条件下的Helmholtz方程的解，其求解只需考虑有限元离散的网格信息，避免了占用大量内存，极大地提高了过滤效率。然而随着网格的不断细分，该滤波器的代数方程维数越来越高，采用以Gauss消元法为基础的直接法求解，会导致耗时过长和内存需求量增加等问题，从而影响该滤波器的效率。此外，随着计算机性能的提高和算法技术的发展，精细的数值模拟成为可能，使得待求解的线性方程组的维数越来越高，研发高效的线性求解器加速算法[9-11]已变得至关重要。

    多重网格是一种高效的求解算法，近期在偏微分方程参数反演问题方面有许多应用[12-13]，并广泛地应用到图像处理[14]、流体模拟[15]和结构拓扑优化[16-17]等领域。为提高大规模结构拓扑优化灵敏度过滤的效率，本文基于体积约束下柔度值最小的固体各向同性惩罚微结构(SIMP)插值模型，针对大规模PDE灵敏度滤波器的求解问题，对其进行有限元分析得到代数方程，再采用共轭梯度算法、多重网格算法以及多重网格预处理共轭梯度算法进行求解，得到了过滤域的节点灵敏度值，通过节点到单元的映射得到了单元灵敏度过滤值。

    1.   优 化 模 型

    在变密度法中，SIMP模型是目前工程领域应用最广泛的材料插值模型，为了使模型更具有代表性，基于SIMP插值模型，建立在体积约束下柔度值最小的连续体结构拓扑优化模型：

    $$\left\{ \begin{aligned} & \mathrm{min}\;{\boldsymbol{C}}={{\boldsymbol{F}}}^{\rm{T}}{\boldsymbol{U}}={{\boldsymbol{U}}}^{\rm{T}}{\boldsymbol{KU}}={\displaystyle \sum _{i=1}^{N}{{\boldsymbol{u}}}_{i}^{\rm{T}}{{\boldsymbol{k}}}_{i}}{{\boldsymbol{u}}}_{i}={\displaystyle \sum _{i=1}^{N}({{\boldsymbol{E}}}_{\mathrm{min}} + {x}_{i}^{{\rm{p}}}({{\boldsymbol{E}}}_{0}-{{\boldsymbol{E}}}_{\mathrm{min}})){{\boldsymbol{u}}}_{i}^{\rm{T}}}{{\boldsymbol{k}}}_{0}{{\boldsymbol{u}}}_{i}, \\&\qquad \text{s}\text{.t}\text{.}\quad V=f{V}_{0}={\displaystyle \sum _{i=1}^{N}{x}_{i}}{v}_{i}\le {V}^{*},\\ &\quad\qquad\quad {\boldsymbol{F}}={\boldsymbol{KU}},\\ &\qquad\quad\quad 0 < {x}_{\mathrm{min}}\le {x}_{i}\le {x}_{\mathrm{max}}\le 1, \end{aligned} \right.$$ (1)

    其中，目标函数C为结构的总体柔度，F为力向量，$ {\boldsymbol{U}} $为位移，K为结构的总刚度矩阵，$ {x_i} $和$ {{\boldsymbol{u}}_i} $是第$ i $个单元的相对密度和节点位移，$ {{\boldsymbol{k}}_i} $和$ {{\boldsymbol{k}}_0} $是单元刚度矩阵，$ V $是优化后的结构体积，$ {V^*} $为结构体积约束，$ {V_0} $为整个设计域的初始体积，$ f $为优化体积比，$ {v_i} $为第$ i $个单元的体积，$ {x_{\min }} $和$ {x_{\max }} $为单元相对密度的最小和最大极限值，$ N $为结构离散单元总数。

    SIMP法又称幂律法，其材料插值模型以密度为设计变量，单元材料弹性模量和密度之间的关系为

    $$ {{\boldsymbol{E}}_i} = {{\boldsymbol{E}}_{\min }} + x_i^{\text{p}}({{\boldsymbol{E}}_0} - {{\boldsymbol{E}}_{\min }}) , $$ (2)

    其中$ {{\boldsymbol{E}}_0} $和$ {{\boldsymbol{E}}_{\min }} $是固体和空洞部分的材料弹性模量，$ {{\boldsymbol{E}}_i} $为第$ i $个单元插值后的弹性模量。

    2.   拓扑优化算法

    2.1   总体框架

    连续体结构拓扑优化实现流程如图1所示。首先，初始化优化参数：定义问题的设计区域、约束和载荷。其次，建立有限元模型：划分有限元网格，将结构离散化，建立式(1)所示的连续体结构拓扑优化数学模型。另外，有限元分析：计算刚度矩阵、节点位移、目标柔度值和灵敏度。之后，PDE灵敏度过滤：解决数值不稳定问题和消除棋盘格现象。然后，采用二分法和OC准则更新Lagrange乘子和设计变量，判断是否满足体积约束条件：满足则输出设计变量；若不满足，则回到上一步再次更新Lagrange乘子和设计变量。最后，检查是否满足收敛准则：若满足，则输出拓扑优化结果；若不满足，则再次进行有限元分析，进行下一轮迭代。

    图  1  基于SIMP插值模型连续体结构拓扑优化实现流程图
    Figure  1.  The flow chart of continuum structure topology optimization based on the SIMP interpolation model
    下载: 全尺寸图片 幻灯片

    为解决传统PDE灵敏度滤波器过滤慢的问题，本文采用共轭梯度算法求解该滤波器的对称正定代数方程，为进一步解决网格数量增加导致系数矩阵维数增加和趋近最优解时，该算法收敛速度减缓的问题，结合多重网格可以均匀衰减误差，加快算法收敛速度的优点，使用多重网格算法和预处理共轭梯度算法求解该方程，提升过滤效率。此外，在预处理共轭梯度算法中，采用多重网格预处理降低系数矩阵的条件数。

    2.2   PDE灵敏度滤波器

    本文采用PDE灵敏度滤波器[18-19]，其核心思想是将过滤域的节点灵敏度值隐式表示为具有Neumann边界条件的Helmholtz方程的隐式解，然后通过节点到单元的映射得到单元灵敏度过滤结果：

    $$ \left\{ \begin{aligned} & { - R_{\min }^2{\nabla ^2}\tilde \varphi + \tilde \varphi = \varphi ,} \\ & \frac{{\partial \tilde \varphi }}{{\partial {\boldsymbol{n}}}} = {\boldsymbol{0}}, \end{aligned}\right. $$ (3)

    其中$ \varphi = x\dfrac{{\partial c}}{{\partial x}},\tilde \varphi = x\dfrac{{\partial \tilde c}}{{\partial x}} $，x是单元相对密度，$ \dfrac{{\partial c}}{{\partial x}} $和$ \dfrac{{\partial \tilde c}}{{\partial x}} $是过滤前后的灵敏度，$ {\boldsymbol{n}} $是边界的法向量，$ {R_{\min }} $是长度参数，用来确定过滤区域，它与标准的过滤半径$ {r_{\min }} $之间满足的关系为

    $$ {R_{\min }} = \frac{{{r_{\min }}}}{{2\sqrt 3 }}{\text{.}} $$ (4)

    2.3   PDE灵敏度滤波器的求解

    PDE灵敏度滤波器的实质是以Neumann边界为条件的Helmholtz偏微分方程，其求解步骤如下：

    第1步 有限元离散得到Helmholtz方程(3)的代数方程：

    $$ {{\boldsymbol{K}}_{\text{F}}}{\tilde x_{\text{N}}} = {{\boldsymbol{T}}_{\text{F}}}\left(x\frac{{\partial c}}{{\partial x}}\right), $$ (5)

    其中$ {{\boldsymbol{K}}_{\text{F}}} $是标量问题的标准有限元刚度矩阵，$ {{\boldsymbol{T}}_{\text{F}}} $是将单元相对密度映射到具有节点值的向量的矩阵。针对代数方程(5)的系数矩阵是对称正定稀疏病态的特点，本文拟采用共轭梯度算法、多重网格算法和多重网格预处理共轭梯度算法求解该代数方程，从而得到过滤区域内的节点灵敏度值$ {\tilde x_{\text{N}}} $。

    第2步 通过对过滤域的节点灵敏度值映射得到单元灵敏度过滤值：

    $$ \tilde {\boldsymbol{x}} = {\boldsymbol{T}}_{\text{F}}^{\rm{T}}{\tilde x_{\text{N}}}{\text{.}}$$ (6)
    2.3.1   多重网格算法

    在固定网格中采用常规的迭代方法求解线性方程组，初始时收敛速度很快，趋近最优解时，收敛速度慢，计算时间长。采用Fourier级数表示误差，得到两种不同的误差：波长较短、波形频率高的高频分量误差和波长较长、波形频率低的低频分量误差。由于波形频率的高低是相对的，把细网格中的低频波形放到粗网格中就有可能变成高频波形。因此采用多重网格来消除误差分量，使各种频率的误差得到比较均匀的衰减，从而加快迭代收敛速度。假设有L重网格，k=1为最细的网格层，$ {\boldsymbol{u}}_k^{m + 1} $是$ {{\boldsymbol{A}}_k}{\boldsymbol{u}}_k^m = {{\boldsymbol{f}}_k} $在k层上的多重网格迭代解，多重网格V循环具体过程如下（伪代码见算法1）：

    步骤1 前光滑。采用常规的迭代方法消除误差高频分量，具体步骤为：在细网格上，采取2 ~ 3次的阻尼Jacobi迭代方法得到线性方程组$ {{\boldsymbol{A}}_k}{\boldsymbol{u}}_k^m = {{\boldsymbol{f}}_k} $的近似解。迭代公式的通式为

    $$ {\boldsymbol{u}}_k^{m + 1} = {\boldsymbol{u}}_k^m + {\boldsymbol{S}}_k^{ - 1}({{\boldsymbol{f}}_k} - {{\boldsymbol{A}}_k}{\boldsymbol{u}}_k^m), $$ (7)

    其中$ {{\boldsymbol{S}}_k} $为光滑子，$ {{\boldsymbol{S}}_k} = \omega {{\boldsymbol{D}}_k} $，$ \omega $为阻尼因子，$ {{\boldsymbol{D}}_k} $为$ {{\boldsymbol{A}}_k} $的主对角元素矩阵。

    步骤2 粗网格校正。计算在粗网格上误差的近似解，具体步骤为：首先计算细网格上的残差${\boldsymbol{r}}_k^m = {{\boldsymbol{f}}_k} - {{\boldsymbol{A}}_k}{\boldsymbol{u}}_k^m$；然后使用限制算子$ {\boldsymbol{R}}_k^{k + 1} $把细网格上的残差限制到粗网格上，得到粗网格上的残差$ {\boldsymbol{r}}_{k + 1}^m = {\boldsymbol{R}}_k^{k + 1}{\boldsymbol{r}}_k^m $；最后根据网格的层数求解线性方程组$ {{\boldsymbol{A}}_{k + 1}}\tilde {\boldsymbol{e}}_{k + 1}^m = {\boldsymbol{R}}_k^{k + 1}{\boldsymbol{r}}_k^m $得到粗网格上误差的近似解，若是最粗的网格层，则利用直接法求解$ \tilde {\boldsymbol{e}}_{k + 1}^m = {\boldsymbol{A}}_{k + 1}^{ - 1}{\boldsymbol{r}}_{k + 1}^m $，否则利用多重网格V循环求解$ \tilde {\boldsymbol{e}}_{k + 1}^m{\text{ = }}V \cdot {\text{cycle(}}{{\boldsymbol{A}}_{k + 1}},{\boldsymbol{r}}_{k + 1}^m,0,k + 1) $，其中粗网格算子$ {{\boldsymbol{A}}_{k + 1}} $取Galerkin算子$ {{\boldsymbol{A}}_{k + 1}} = {\boldsymbol{R}}_k^{k + 1}{{\boldsymbol{A}}_k}{\boldsymbol{P}}_{k + 1}^k $。

    步骤3 插值。计算细网格上的解，具体步骤为：先使用插值算子$ {\boldsymbol{P}}_{k + 1}^k $把粗网格上的近似误差插值到细网格上$ {\boldsymbol{P}}_{k + 1}^k\tilde {\boldsymbol{e}}_{k + 1}^m = \tilde {\boldsymbol{e}}_k^m $，然后使用公式$ {\boldsymbol{u}}_k^{m + 1} = {\boldsymbol{u}}_k^{m + 1} + \tilde {\boldsymbol{e}}_k^m $得到在细网格的解。本文采用的插值算子$ {\boldsymbol{P}}_{k + 1}^k $为

    $$ \left[ {\begin{array}{*{20}{c}} {{v_1}} \\ {{v_2}} \\ {{v_3}} \\ {{v_4}} \\ {{v_5}} \\ {{v_6}} \\ {{v_7}} \\ {{v_8}} \\ {{v_9}} \end{array}} \right] = {\left[ {\begin{array}{*{20}{c}} 1&{\dfrac{1}{2}}&0&{\dfrac{1}{2}}&{\dfrac{1}{4}}&0&0&0&0 \\ 0&{\dfrac{1}{2}}&1&0&{\dfrac{1}{4}}&{\dfrac{1}{2}}&0&0&0 \\ 0&0&0&{\dfrac{1}{2}}&{\dfrac{1}{4}}&0&1&{\dfrac{1}{2}}&0 \\ 0&0&0&0&{\dfrac{1}{4}}&{\dfrac{1}{2}}&0&{\dfrac{1}{2}}&1 \end{array}} \right]^{\text{T}}}\left[ {\begin{array}{*{20}{c}} {{u_1}} \\ {{u_2}} \\ {{u_3}} \\ {{u_4}} \end{array}} \right]{\text{ = }}{\boldsymbol{P}}_{k + 1}^k\left[ {\begin{array}{*{20}{c}} {{u_1}} \\ {{u_2}} \\ {{u_3}} \\ {{u_4}} \end{array}} \right]{\text{.}} $$ (8)

    限制算子为$ {\boldsymbol{P}}_{k + 1}^k = {({\boldsymbol{P}}_{k + 1}^k)^{\rm{T}}} $，其中$ {u}_{1}，{u}_{2}，{u}_{3}，{u}_{4} $是粗网格上的节点，$ {v_i}\left( {i = 1,2, \cdots ,9} \right) $是细网格上的节点。

    步骤4 后光滑。为保证整个迭代对称，在细网格上，采取具有相同参数的阻尼Jacobi迭代法求解线性方程组$ {{\boldsymbol{A}}_k}{\boldsymbol{u}}_k^{m + 1} = {{\boldsymbol{f}}_k} $，即$ {\boldsymbol{u}}_k^{m + 1} = {\boldsymbol{u}}_k^{m + 1} + {\boldsymbol{S}}_k^{ - 1}({{\boldsymbol{f}}_k} - {{\boldsymbol{A}}_k}{\boldsymbol{u}}_k^{m + 1}) $。

    　　算法1
    　　1 input:
    　　2　level k;
    　　3　initial guess $ {\boldsymbol{u}}_k^m $;
    　　4　right hand side vector $ {{\boldsymbol{f}}_k} $;
    　　5　coefficient matrix under level $ k $ grid $ {{\boldsymbol{A}}_k} $;
    　　6 output:
    　　7　updated solution $ {\boldsymbol{u}}_k^{m + 1} $
    　　8 pre-smooth: $ {\boldsymbol{u}}_k^{m + 1} = {{\boldsymbol{S}}^{{\mu _1}}}({{\boldsymbol{A}}_k},{{\boldsymbol{f}}_k},{\boldsymbol{u}}_k^m) $
    　　9 coarse grid correction:
    　　10　compute residual: $ {\boldsymbol{r}}_k^m = {{\boldsymbol{f}}_k} - {{\boldsymbol{A}}_k}{\boldsymbol{u}}_k^m $
    　　11　restrict the residual to coarse-grid: $ {\boldsymbol{r}}_{k + 1}^m = {\boldsymbol{R}}_k^{k + 1}{\boldsymbol{r}}_k^m $
    　　12　if k + 1=L then
    　　13　　solve approximate error solution $ \tilde {\boldsymbol{e}}_{k + 1}^m = {\boldsymbol{A}}_{k + 1}^{ - 1}{\boldsymbol{r}}_{k + 1}^m $
    　　14　else
    　　15　　$ \tilde {\boldsymbol{e}}_{k + 1}^m{\text{ = }}V \cdot {\text{cycle(}}{{\boldsymbol{A}}_{k + 1}},{\boldsymbol{r}}_{k + 1}^m,0,k + 1) $
    　　16　end
    　　17 interpolate:
    　　18　$ \tilde {\boldsymbol{e}}_k^m = P_{k + 1}^k\tilde {\boldsymbol{e}}_{k + 1}^m $
    　　19　$ {\boldsymbol{u}}_k^{m + 1} = {\boldsymbol{u}}_k^{m + 1} + \tilde {\boldsymbol{e}}_k^m $
    　　20 post-smooth:
    　　21　$ {\boldsymbol{u}}_k^{m + 1} = {{\boldsymbol{S}}^{{\mu _2}}}({{\boldsymbol{A}}_k},{{\boldsymbol{f}}_k},{\boldsymbol{u}}_k^{m + 1}) $。

    2.3.2   多重网格预处理共轭梯度算法

    考虑到代数方程的系数矩阵是大型、稀疏、对称正定的，有时甚至为病态的，则重点研究了多重网格预处理共轭梯度算法（伪代码见算法2），采用多重网格作为预处理方法降低系数矩阵的条件数。

    　　算法2
    　　input:
    　　　symmetric positive definite coefficient matrix $ {\boldsymbol{A}} $;
    　　　right hand side vector $ {\boldsymbol{b}} $;
    　　　initiate a guess for $ {{\boldsymbol{x}}_0} \in {\mathbb{R}^n} $
    　　output:
    　　　solution $ {\boldsymbol{x}} $;
    　　1 compute the residual on the finest grid $ {{\boldsymbol{r}}_0} \mathop {\rm{ = }}\limits^\Delta {\boldsymbol{b}} - {\boldsymbol{A}}{{\boldsymbol{x}}_0} $;
    　　2 for $ i = 1,2,3, \cdots $ do
    　　3 solve $ {\boldsymbol{A}}{{\boldsymbol{z}}_0} = {{\boldsymbol{r}}_0} $ by the multigrid V-cycle algorithm
    　　4　initialize $ {{\boldsymbol{p}}_0} = {{\boldsymbol{z}}_0} $
    　　5　compute a step length the $ {\alpha _i} = ({\boldsymbol{r}}_{i - 1}^{\text{T}}{{\boldsymbol{z}}_{i - 1}})/({\boldsymbol{p}}_{i - 1}^{\text{T}}{\boldsymbol{A}}{{\boldsymbol{p}}_{i - 1}}) $
    　　 　 (where $ {{\boldsymbol{z}}_{i - 1}} $ can be solved by the multigrid V-cycle algorithm)
    　　6　update the approximate solution $ {{\boldsymbol{x}}_i} = {{\boldsymbol{x}}_{i - 1}} + {\alpha _i}{{\boldsymbol{p}}_{i - 1}} $
    　　7　update the residual $ {{\boldsymbol{r}}_i} = {{\boldsymbol{r}}_{i - 1}} - {\alpha _i}{\boldsymbol{A}}{{\boldsymbol{p}}_{i - 1}} $
    　　8　solve $ {\boldsymbol{A}}{{\boldsymbol{z}}_i} = {{\boldsymbol{r}}_i} $ by the multigrid V-cycle algorithm
    　　9　compute a gradient correction factor $ {\beta _i} = ({\boldsymbol{r}}_i^{\rm{T}}{{\boldsymbol{z}}_i})/({\boldsymbol{r}}_{i - 1}^{\rm{T}}{{\boldsymbol{z}}_{i - 1}}) $
    　　10　set the new search direction $ {{\boldsymbol{p}}_i} = {{\boldsymbol{z}}_i} + {\beta _i}{{\boldsymbol{p}}_{i - 1}} $
    　　11 end
    　　12 until convergence。

    2.4   设计变量的更新

    优化准则法（OC准则）具有收敛速度快、易于实现等特点。本文采取OC准则求解优化模型(1)，基本思路为，先引入Lagrange乘子把约束问题(1)转化为无约束问题，再根据极值点存在的条件即K-T条件以及不动点迭代得到设计变量的迭代公式为

    $$ {({x_i})^{k + 1}} = \left\{ \begin{aligned} & \min \{ {{({x_i})}^k+m},1\} ,\qquad {{({x_i})}^k}{{(D_i^{(k)})}^\eta }\geqslant\min \{ {{({x_i})}^k}+m,1\} , \\ & {{({x_i})}^k}{{(D_i^{(k)})}^\eta },\qquad\max \{ {{({x_i})}^k}-m,{x_{\min }}\} \leqslant {{({x_i})}^k}{{(D_i^{(k)})}^\eta } < \min \{ {{({x_i})}^k}+m,1\} , \\ & \max \{ {{({x_i})}^k}-m,{x_{\min }}\} ,\qquad{{({x_i})}^k}{{(D_i^{(k)})}^\eta } < \max \{ {{({x_i})}^k}-m,{x_{\min }}\} , \end{aligned} \right. $$ (9)

    其中${(D_i^{(k)})^\eta } = \dfrac{{p{{({x_i})}^{p - 1}}{\boldsymbol{u}}_i^{\text{T}}{{\boldsymbol{k}}_0}{{\boldsymbol{u}}_i}}}{{\mu {v_i}}}$，$ \eta $为阻尼因子，$ m $为移动极限。

    2.5   收敛准则

    包括外循环和内循环，外循环包括有限元分析和更新单元灵敏度，迭代收敛条件如下：

    $$ \left| {x_{\max }^{(k + 1)} - x_{\max }^{(k)}} \right| \leqslant \varepsilon,\qquad \varepsilon = 0.01, $$ (10)

    其中$ \varepsilon $为收敛控制因子，$ x_{\max }^{(k)} $和$ x_{\max }^{(k + 1)} $分别为第$ k $次和第$ k + 1 $次设计变量最大值。

    内循环包括更新设计变量、Lagrange乘子和过滤域的节点值。采用二分法和OC准则更新Lagrange乘子和设计变量，以满足模型(1)中的体积约束和设计变量的边界约束，迭代收敛条件如下：

    $$ \frac{{{l_2} - {l_1}}}{{{l_2} + {l_1}}} \leqslant \eta, \qquad \eta = 0.001, $$ (11)

    其中$ \eta $为收敛控制因子，$ {l_2} $和$ {l_1} $为较大和较小的Lagrange乘子。采用共轭梯度算法、多重网格算法以及多重网格预处理共轭梯度算法求解代数方程，构造迭代收敛准则为

    $$ \frac{{{{\left\| {{{\boldsymbol{r}}_k}} \right\|}_2}}}{{{{\left\| {{{\boldsymbol{r}}_0}} \right\|}_2}}} \leqslant \tau , $$ (12)

    其中$ \tau $为收敛控制因子，$ {{\boldsymbol{r}}_k} $和$ {{\boldsymbol{r}}_0} $为第$ k $次和初始时代数方程的残量。

    3.   数 值 算 例

    为了验证本文拓扑优化方法的可行性、可靠性和高效性，针对长960 mm、高320 mm的悬臂梁和长960 mm、高480 mm Michell结构的棋盘格现象、数值不稳定、灵敏度过滤慢等问题开展了数值实验，如图2和3所示（算例优化参数为体积约束比为0.5、弹性模量为1 GPa、Poisson比为0.3、惩罚因子为3、移动极限为0.2、外载荷F为1 kN）。拓扑优化程序基于文献[16,20]的MATLAB代码进行编写，在PDE滤波器的求解中引入了共轭梯度（CG）算法、多重网格（MG）算法以及多重网格预处理共轭梯度（MGCG）算法，并比较了三种算法对结构拓扑优化效率的影响，相关数值结果如下。

    图  2  悬臂梁
    Figure  2.  The cantilever beam
    下载: 全尺寸图片 幻灯片
    图  3  Michell结构
    Figure  3.  The Michell structure
    下载: 全尺寸图片 幻灯片

    3.1   精度的影响

    依次取4种不同的误差精度P(P=1E–2,1E–3,1E–4,1E–5)，采用CG算法、网格层数为5的MG算法和MGCG算法求解经有限元分析Helmholtz方程得到的代数方程，并对整个结构分别划分为$ 480 \times 160 $和$ 480 \times 240 $个单元的悬臂梁和Michell结构进行半径为6的PDE灵敏度过滤拓扑优化实验。其中，不同精度下悬臂梁和Michell结构的优化结构、目标柔度值C、整体迭代次数I以及整体优化时间T的分布如图4、图5、表1、图6、图7所示。

    图  4  不同精度下CG、MG和MGCG算法的悬臂梁的优化结构：(a) P=1E−2 (CG)；(b) P=1E−2 (MG)；(c) P=1E−2 (MGCG)；(d) P=1E−3 (CG)；(e) P=1E−3 (MG)；(f) P=1E−3 (MGCG)；(g) P=1E−4 (CG)；(h) P=1E−4 (MG)；(i) P=1E−4 (MGCG)；(j) P=1E−5 (CG)；(k) P=1E−5 (MG)；(l) P=1E−5 (MGCG)
    Figure  4.  The optimized cantilever beams by CG, MG and MGCG algorithms with different degrees of precision: (a) P=1E−2 (CG); (b) P=1E−2 (MG); (c) P=1E−2 (MGCG); (d) P=1E−3 (CG); (e) P=1E−3 (MG); (f) P=1E−3 (MGCG); (g) P=1E−4 (CG); (h) P=1E−4 (MG); (i) P=1E−4 (MGCG); (j) P=1E−5 (CG); (k) P=1E−5 (MG); (l) P=1E−5 (MGCG)
    下载: 全尺寸图片 幻灯片
    图  5  不同精度下CG、MG和MGCG算法的Michell结构的优化结构：(a) P=1E−2 (CG)；(b) P=1E−2 (MG)；(c) P=1E−2 (MGCG)；(d) P=1E−3 (CG)；(e) P=1E−3 (MG)；(f) P=1E−3 (MGCG)；(g) P=1E−4 (CG)；(h) P=1E−4 (MG)；(i) P=1E−4 (MGCG)；(j) P=1E−5 (CG)；(k) P=1E−5 (MG)； (l) P=1E−5 (MGCG)
    Figure  5.  The optimized Michell structures by CG, MG and MGCG algorithms with different degrees of precision: (a) P=1E−2 (CG); (b) P=1E−2 (MG); (c) P=1E−2 (MGCG); (d) P=1E−3 (CG); (e) P=1E−3 (MG); (f) P=1E−3 (MGCG); (g) P=1E−4 (CG); (h) P=1E−4 (MG); (i) P=1E−4 (MGCG); (j) P=1E−5 (CG); (k) P=1E−5 (MG); (l) P=1E−5 (MGCG)
    下载: 全尺寸图片 幻灯片
    表  1  两种结构的拓扑优化迭代次数和柔度值
    Table  1.  Total numbers of iterations and compliance for topology optimization of 2 structures
    structureprecision PCG algorithmMG algorithmMGCG algorithm
    compliance Ctime T/siterations Icompliance Ctime T/siterations Icompliance Ctime T/siterations I
    cantilever beam1E−2186.885 287170187.306 996167184.738 364114
    1E−3186.883 482155186.939 3100160184.414 754100
    1E−4186.884 284155186.939 2105160184.014 075133
    1E−5186.883 984155186.939 1111160185.215 564113
    Michell structure1E−217.125 233138716.909 5707816.923 86878
    1E−317.123 423127117.119 926027716.923 96778
    1E−417.123 123027117.119 926927716.952 68093
    1E−517.123 123127117.074 726726716.943 66271
    下载: 导出CSV 
    | 显示表格
    图  6  悬臂梁结构优化时间分布图
    Figure  6.  The optimized time distribution of the cantilever beam
    下载: 全尺寸图片 幻灯片
    图  7  Michell结构优化时间分布图
    Figure  7.  The optimized time distribution of the Michell structure
    下载: 全尺寸图片 幻灯片

    从图4、图5和表1可以看出，在过滤半径为6的PDE灵敏度滤波器中使用不同精度的CG算法、MG算法和MGCG算法得到的悬臂梁和Michell结构的优化结构图基本相同，并且目标函数值非常接近，至多相差1%，说明PDE滤波器求解算法的精度对拓扑优化的结果影响不大。此外，在PDE滤波器中使用MGCG算法进行结构拓扑优化的迭代次数整体低于其他两种算法。

    从图6和图7可以看出，同一精度下，在PDE滤波器中使用MGCG算法进行结构拓扑优化的时间整体低于CG算法和MGCG算法，且随着精度的增加，在PDE滤波器中使用CG算法进行结构拓扑优化的时间逐渐降低并趋于稳定，在PDE滤波器中使用MG算法进行结构拓扑优化的时间逐渐增加并趋于稳定。此外，在PDE滤波器中采用MGCG算法对悬臂梁和Michell结构拓扑优化的时间分别比CG算法少了11% ~ 34%和65% ~ 80%，比MG算法少了29% ~ 46%，3%和77%。因此，从结构的整体优化效率看，MGCG算法高于CG算法和MG算法。

    3.2   过滤半径的影响

    依次取4种不同的过滤半径R(R=3，6，9，12), 采用CG算法、网格层数为5的MG算法以及MGCG算法求解有限元分析Helmholtz方程得到的代数方程，并对整个结构分别划分为$ 480 \times 160 $和$ 480 \times 240 $个单元的悬臂梁和Michell结构进行控制精度为1E–4的PDE灵敏度过滤的拓扑优化实验。其中，不同过滤半径下悬臂梁和Michell结构的优化结构、目标柔度值C、整体迭代次数I以及整体优化时间T的分布分别如图8、图9、表2、图10、图11所示。

    图  8  不同过滤半径下CG、MG和MGCG算法的悬臂梁的优化结构：(a) R=3 (CG)；(b) R=3 (MG)；(c) R=3 (MGCG)；(d) R=6 (CG)；(e) R=6 (MG)；(f) R=6 (MGCG)；(g) R=9 (CG)；(h) R=9 (MG)；(i) R=9 (MGCG)；(j) R=12 (CG)；(k) R=12 (MG)；(l) R=12 (MGCG)
    Figure  8.  The optimized cantilever beams by CG, MG and MGCG algorithms with different filter radius: (a) R=3 (CG); (b) R=3 (MG); (c) R=3 (MGCG); (d) R=6 (CG); (e) R=6 (MG); (f) R=6 (MGCG); (g) R=9 (CG); (h) R=9 (MG); (i) R=9 (MGCG); (j) R=12 (CG); (k) R=12 (MG); (l) R=12 (MGCG)
    下载: 全尺寸图片 幻灯片
    图  9  不同过滤半径下CG、MG和MGCG算法的Michell结构的优化结构：(a) R=3 (CG)；(b) R=3 (MG)；(c) R=3 (MGCG)；(d) R=6 (CG)；(e) R=6 (MG)；(f) R=6 (MGCG)；(g) R=9 (CG)；(h) R=9 (MG)；(i) R=9 (MGCG)；(j) R=12 (CG)；(k) R=12 (MG)；(l) R=12 (MGCG)
    Figure  9.  The optimized Michell structures by CG, MG and MGCG algorithms with different filter radius: (a) R=3 (CG); (b) R=3 (MG); (c) R=3 (MGCG); (d) R=6 (CG); (e) R=6 (MG); (f) R=6 (MGCG); (g) R=9 (CG); (h) R=9 (MG); (i) R=9 (MGCG); (j) R=12 (CG); (k) R=12 (MG); (l) R=12 (MGCG)
    下载: 全尺寸图片 幻灯片

    从图8、图9和表2可以看出，在PDE灵敏度滤波器中使用精度为1E–4的CG算法、MG算法和MGCG算法得到的悬臂梁和Michell结构的优化结构图、目标函数值的特征随着过滤半径的变化呈现类似的变化规律。随着过滤半径的增大，拓扑优化结构的清晰度越来越低，细支结构越来越少，说明该两种算法对结构承力路径的识别精度与过滤半径间体现出负相关性，且柔度值越来越大。且从表2可以看出，随着过滤半径的改变，在PDE滤波器中使用MG算法和MGCG算法进行结构拓扑优化的迭代次数整体低于CG算法。

    表  2  两种结构的拓扑优化迭代次数和柔度值(不同过滤半径)
    Table  2.  Total numbers of iterations and compliance for topology optimization of 2 structures (different filter radius)
    structurefilter radius RCG algorithmMG algorithmMGCG algorithm
    compliance Ctime T/siterations Icompliance Ctime T/siterations Icompliance Ctime T/siterations I
    cantilever beam3180.759 1505917182.218 3147225182.352 058103
    6186.884 282155186.939 2101160184.014 073133
    9192.522 998183193.053 764100193.096 75397
    12197.418 773136197.669 567105197.641 962115
    Michell structure316.738 584092616.864 516516816.894 18090
    617.123 123927117.119 928027716.952 68493
    917.382 911713317.322 611511717.168 65967
    1217.696 49911417.696 510711417.686 4100114
    下载: 导出CSV 
    | 显示表格
    图  10  悬臂梁结构优化时间分布图(不同过滤半径)
    Figure  10.  The optimized time distribution of the cantilever beam (different filter radius)
    下载: 全尺寸图片 幻灯片
    图  11  Michell结构优化时间分布图(不同过滤半径)
    Figure  11.  The optimized time distribution of the Michell structure (different filter radius)
    下载: 全尺寸图片 幻灯片

    从图10和图11可以看出，同一过滤半径下，在PDE滤波器中使用MG算法和MGCG算法进行结构拓扑优化的时间整体低于CG算法，且随着过滤半径增大，MGCG算法的优势逐渐降低。此外，在PDE滤波器中采用MGCG算法对悬臂梁和Michell结构拓扑优化的时间分别比CG算法少了11%~89%和50%~91%，比MG算法少了7%~61%，7%和70%。因此，从结构的整体优化效率看，MGCG算法高于CG算法和MG算法。

    3.3   网格数量的影响

    分别对网格划分为$ 960 \times 320 $，过滤半径为12，层数为6；网格划分为$ 480 \times 160 $，过滤半径为6，层数为5；网格划分为$ 240 \times 80 $，过滤半径为3，层数为4的悬臂梁和网格划分为$ 960 \times 480 $，过滤半径为12，层数为6；网格划分为$ 480 \times 240 $，过滤半径为6，层数为5；网格划分为$ 240 \times 120 $，过滤半径为3，层数为4的Michell结构进行控制精度为1E–4的PDE灵敏度过滤的拓扑优化实验。其中，不同网格数量下悬臂梁和Michell结构的优化结构、目标柔度值C、整体迭代次数I以及整体优化时间T的分布分别如图12、图13、表3、图14、图15所示。

    图  12  不同网格数量下CG、MG和MGCG算法的悬臂梁的优化结构：(a) 网格为240 × 80 (CG)；(b) 网格为240 × 80 (MG)；(c) 网格为240 × 80 (MGCG)；(d) 网格为480 × 160 (CG)；(e) 网格为480 × 160 (MG)；(f) 网格为480 × 160 (MGCG)；(g) 网格为960 × 320 (CG)；(h) 网格为960 × 320 (MG)；(i) 网格为960 × 320 (MGCG)
    Figure  12.  The optimized cantilever beams by CG, MG and MGCG algorithms with different grids: (a) grid 240 × 80 (CG); (b) grid 240 × 80 (MG); (c) grid 240 × 80 (MGCG); (d) grid 480 × 160 (CG); (e) grid 480 × 160 (MG); (f) grid 480 × 160 (MGCG); (g) grid 960 × 320 (CG); (h) grid 960 × 320 (MG); (i) grid 960 × 320 (MGCG)
    下载: 全尺寸图片 幻灯片
    图  13  不同网格数量下CG、MG和MGCG算法的Michell结构的优化结构：(a) 网格为240 × 120 (CG)；(b) 网格为240 × 120 (MG)；(c) 网格为240 × 120 (MGCG)；(d) 网格为480 × 240 (CG)；(e) 网格为480 × 240 (MG)；(f) 网格为480 × 240 (MGCG)；(g) 网格为960 × 480 (CG)；(h) 网格为960 × 480 (MG)；(i) 网格为960 × 480 (MGCG)
    Figure  13.  The optimized Michell structures by CG, MG and MGCG algorithms with different grids: (a) grid 240 × 120 (CG); (b) grid 240 × 120 (MG); (c) grid 240 × 120 (MGCG); (d) grid 480 × 240 (CG); (e) grid 480 × 240 (MG); (f) grid 480 × 240 (MGCG); (g) grid 960 × 480 (CG); (h) grid 960 × 480 (MG); (i) grid 960 × 480 (MGCG)
    下载: 全尺寸图片 幻灯片
    表  3  两种结构的拓扑优化迭代次数和柔度值(不同网格数量)
    Table  3.  Total numbers of iterations and compliance for topology optimization of 2 structures (different grids)
    structuregrid sizeCG algorithmMG algorithmMGCG algorithm
    compliance Ctime T/siterations Icompliance Ctime T/siterations Icompliance Ctime T/siterations I
    cantilever beam240 × 80186.484 016133186.315 319133183.148 41296
    480 × 160186.884 282155186.939 2101160184.014 073133
    960 × 320188.099 5430167188.513 5495173185.695 5419155
    Michell structure240 × 12015.995 44223315.867 0209115.816 01583
    480 × 24017.123 123527117.119 927427716.952 68093
    960 × 48018.293 9131929218.288 6148 130618.095 936191
    下载: 导出CSV 
    | 显示表格
    图  14  悬臂梁结构优化时间分布图(不同网格数量)
    Figure  14.  The optimized time distribution of the cantilever beam (different grids)
    下载: 全尺寸图片 幻灯片
    图  15  Michell结构优化时间分布图(不同网格数量)
    Figure  15.  The optimized time distribution of the Michell structure (different grids)
    下载: 全尺寸图片 幻灯片

    从图12、图13和表3可以看出，在PDE灵敏度滤波器中使用精度为1E–4的CG算法、MG算法和MGCG算法得到不同网格数量的悬臂梁和Michell结构的优化结构基本相同，并且随着网格数量增加，柔度值越来越大。可见PDE滤波器求解算法对不同网格数量的优化结构影响不大，且可有效地消除网格依赖性和棋盘格现象。此外，在PDE滤波器中使用MGCG算法进行结构拓扑优化的迭代次数整体低于其他两种算法。

    从图14和图15可以看出，同一网格数量下，在PDE滤波器中使用MGCG算法进行结构拓扑优化的时间整体低于CG算法和MGCG算法。特别地，在PDE滤波器中采用MGCG算法对悬臂梁和Michell结构拓扑优化的时间分别比CG算法少了3% ~ 25%和64% ~ 73%，比MG算法少了15% ~ 37%，25%和76%。因此，从结构的整体优化效率看，MGCG算法高于CG算法和MG算法。

    4.   结　　论

    本文在PDE灵敏度滤波器中引入了CG算法、MG算法和MGCG算法，对悬臂梁和Michell结构进行拓扑优化，并比较了算法精度、过滤半径和网格数量对优化结果和优化效率的影响。实验结果表明：

    1） 不同精度的三种算法得到的优化结构、目标函数值基本相同。随着精度的提升，MG算法的结构优化时间不断增加，CG算法总体趋于稳定，MGCG算法有微小的波动，且MGCG算法结构优化时间整体少于其他两种算法。

    2） 随着过滤半径的增大，CG算法和MG算法的结构优化时间大幅减少，整体呈下降趋势，MGCG算法有微小的波动，但优化时间整体少于其他两种算法，且过滤半径越小，优势越明显。

    3） 不同网格数量的三种算法得到的优化结构、目标函数值基本相同。随着网格数量的增加，三种算法的结构优化时间和迭代次数整体呈上升趋势，且随着网格数量的增加，MGCG算法的优势越明显。

    因此，从结构的整体优化效率看，MGCG算法优于CG算法和MG算法。

  • 图  1  基于SIMP插值模型连续体结构拓扑优化实现流程图

    Figure  1.  The flow chart of continuum structure topology optimization based on the SIMP interpolation model

    下载: 全尺寸图片 幻灯片

    图  2  悬臂梁

    Figure  2.  The cantilever beam

    下载: 全尺寸图片 幻灯片

    图  3  Michell结构

    Figure  3.  The Michell structure

    下载: 全尺寸图片 幻灯片

    图  4  不同精度下CG、MG和MGCG算法的悬臂梁的优化结构：(a) P=1E−2 (CG)；(b) P=1E−2 (MG)；(c) P=1E−2 (MGCG)；(d) P=1E−3 (CG)；(e) P=1E−3 (MG)；(f) P=1E−3 (MGCG)；(g) P=1E−4 (CG)；(h) P=1E−4 (MG)；(i) P=1E−4 (MGCG)；(j) P=1E−5 (CG)；(k) P=1E−5 (MG)；(l) P=1E−5 (MGCG)

    Figure  4.  The optimized cantilever beams by CG, MG and MGCG algorithms with different degrees of precision: (a) P=1E−2 (CG); (b) P=1E−2 (MG); (c) P=1E−2 (MGCG); (d) P=1E−3 (CG); (e) P=1E−3 (MG); (f) P=1E−3 (MGCG); (g) P=1E−4 (CG); (h) P=1E−4 (MG); (i) P=1E−4 (MGCG); (j) P=1E−5 (CG); (k) P=1E−5 (MG); (l) P=1E−5 (MGCG)

    下载: 全尺寸图片 幻灯片

    图  5  不同精度下CG、MG和MGCG算法的Michell结构的优化结构：(a) P=1E−2 (CG)；(b) P=1E−2 (MG)；(c) P=1E−2 (MGCG)；(d) P=1E−3 (CG)；(e) P=1E−3 (MG)；(f) P=1E−3 (MGCG)；(g) P=1E−4 (CG)；(h) P=1E−4 (MG)；(i) P=1E−4 (MGCG)；(j) P=1E−5 (CG)；(k) P=1E−5 (MG)； (l) P=1E−5 (MGCG)

    Figure  5.  The optimized Michell structures by CG, MG and MGCG algorithms with different degrees of precision: (a) P=1E−2 (CG); (b) P=1E−2 (MG); (c) P=1E−2 (MGCG); (d) P=1E−3 (CG); (e) P=1E−3 (MG); (f) P=1E−3 (MGCG); (g) P=1E−4 (CG); (h) P=1E−4 (MG); (i) P=1E−4 (MGCG); (j) P=1E−5 (CG); (k) P=1E−5 (MG); (l) P=1E−5 (MGCG)

    下载: 全尺寸图片 幻灯片

    图  6  悬臂梁结构优化时间分布图

    Figure  6.  The optimized time distribution of the cantilever beam

    下载: 全尺寸图片 幻灯片

    图  7  Michell结构优化时间分布图

    Figure  7.  The optimized time distribution of the Michell structure

    下载: 全尺寸图片 幻灯片

    图  8  不同过滤半径下CG、MG和MGCG算法的悬臂梁的优化结构：(a) R=3 (CG)；(b) R=3 (MG)；(c) R=3 (MGCG)；(d) R=6 (CG)；(e) R=6 (MG)；(f) R=6 (MGCG)；(g) R=9 (CG)；(h) R=9 (MG)；(i) R=9 (MGCG)；(j) R=12 (CG)；(k) R=12 (MG)；(l) R=12 (MGCG)

    Figure  8.  The optimized cantilever beams by CG, MG and MGCG algorithms with different filter radius: (a) R=3 (CG); (b) R=3 (MG); (c) R=3 (MGCG); (d) R=6 (CG); (e) R=6 (MG); (f) R=6 (MGCG); (g) R=9 (CG); (h) R=9 (MG); (i) R=9 (MGCG); (j) R=12 (CG); (k) R=12 (MG); (l) R=12 (MGCG)

    下载: 全尺寸图片 幻灯片

    图  9  不同过滤半径下CG、MG和MGCG算法的Michell结构的优化结构：(a) R=3 (CG)；(b) R=3 (MG)；(c) R=3 (MGCG)；(d) R=6 (CG)；(e) R=6 (MG)；(f) R=6 (MGCG)；(g) R=9 (CG)；(h) R=9 (MG)；(i) R=9 (MGCG)；(j) R=12 (CG)；(k) R=12 (MG)；(l) R=12 (MGCG)

    Figure  9.  The optimized Michell structures by CG, MG and MGCG algorithms with different filter radius: (a) R=3 (CG); (b) R=3 (MG); (c) R=3 (MGCG); (d) R=6 (CG); (e) R=6 (MG); (f) R=6 (MGCG); (g) R=9 (CG); (h) R=9 (MG); (i) R=9 (MGCG); (j) R=12 (CG); (k) R=12 (MG); (l) R=12 (MGCG)

    下载: 全尺寸图片 幻灯片

    图  10  悬臂梁结构优化时间分布图(不同过滤半径)

    Figure  10.  The optimized time distribution of the cantilever beam (different filter radius)

    下载: 全尺寸图片 幻灯片

    图  11  Michell结构优化时间分布图(不同过滤半径)

    Figure  11.  The optimized time distribution of the Michell structure (different filter radius)

    下载: 全尺寸图片 幻灯片

    图  12  不同网格数量下CG、MG和MGCG算法的悬臂梁的优化结构：(a) 网格为240 × 80 (CG)；(b) 网格为240 × 80 (MG)；(c) 网格为240 × 80 (MGCG)；(d) 网格为480 × 160 (CG)；(e) 网格为480 × 160 (MG)；(f) 网格为480 × 160 (MGCG)；(g) 网格为960 × 320 (CG)；(h) 网格为960 × 320 (MG)；(i) 网格为960 × 320 (MGCG)

    Figure  12.  The optimized cantilever beams by CG, MG and MGCG algorithms with different grids: (a) grid 240 × 80 (CG); (b) grid 240 × 80 (MG); (c) grid 240 × 80 (MGCG); (d) grid 480 × 160 (CG); (e) grid 480 × 160 (MG); (f) grid 480 × 160 (MGCG); (g) grid 960 × 320 (CG); (h) grid 960 × 320 (MG); (i) grid 960 × 320 (MGCG)

    下载: 全尺寸图片 幻灯片

    图  13  不同网格数量下CG、MG和MGCG算法的Michell结构的优化结构：(a) 网格为240 × 120 (CG)；(b) 网格为240 × 120 (MG)；(c) 网格为240 × 120 (MGCG)；(d) 网格为480 × 240 (CG)；(e) 网格为480 × 240 (MG)；(f) 网格为480 × 240 (MGCG)；(g) 网格为960 × 480 (CG)；(h) 网格为960 × 480 (MG)；(i) 网格为960 × 480 (MGCG)

    Figure  13.  The optimized Michell structures by CG, MG and MGCG algorithms with different grids: (a) grid 240 × 120 (CG); (b) grid 240 × 120 (MG); (c) grid 240 × 120 (MGCG); (d) grid 480 × 240 (CG); (e) grid 480 × 240 (MG); (f) grid 480 × 240 (MGCG); (g) grid 960 × 480 (CG); (h) grid 960 × 480 (MG); (i) grid 960 × 480 (MGCG)

    下载: 全尺寸图片 幻灯片

    图  14  悬臂梁结构优化时间分布图(不同网格数量)

    Figure  14.  The optimized time distribution of the cantilever beam (different grids)

    下载: 全尺寸图片 幻灯片

    图  15  Michell结构优化时间分布图(不同网格数量)

    Figure  15.  The optimized time distribution of the Michell structure (different grids)

    下载: 全尺寸图片 幻灯片
    　　算法1
    　　1 input:
    　　2　level k;
    　　3　initial guess $ {\boldsymbol{u}}_k^m $;
    　　4　right hand side vector $ {{\boldsymbol{f}}_k} $;
    　　5　coefficient matrix under level $ k $ grid $ {{\boldsymbol{A}}_k} $;
    　　6 output:
    　　7　updated solution $ {\boldsymbol{u}}_k^{m + 1} $
    　　8 pre-smooth: $ {\boldsymbol{u}}_k^{m + 1} = {{\boldsymbol{S}}^{{\mu _1}}}({{\boldsymbol{A}}_k},{{\boldsymbol{f}}_k},{\boldsymbol{u}}_k^m) $
    　　9 coarse grid correction:
    　　10　compute residual: $ {\boldsymbol{r}}_k^m = {{\boldsymbol{f}}_k} - {{\boldsymbol{A}}_k}{\boldsymbol{u}}_k^m $
    　　11　restrict the residual to coarse-grid: $ {\boldsymbol{r}}_{k + 1}^m = {\boldsymbol{R}}_k^{k + 1}{\boldsymbol{r}}_k^m $
    　　12　if k + 1=L then
    　　13　　solve approximate error solution $ \tilde {\boldsymbol{e}}_{k + 1}^m = {\boldsymbol{A}}_{k + 1}^{ - 1}{\boldsymbol{r}}_{k + 1}^m $
    　　14　else
    　　15　　$ \tilde {\boldsymbol{e}}_{k + 1}^m{\text{ = }}V \cdot {\text{cycle(}}{{\boldsymbol{A}}_{k + 1}},{\boldsymbol{r}}_{k + 1}^m,0,k + 1) $
    　　16　end
    　　17 interpolate:
    　　18　$ \tilde {\boldsymbol{e}}_k^m = P_{k + 1}^k\tilde {\boldsymbol{e}}_{k + 1}^m $
    　　19　$ {\boldsymbol{u}}_k^{m + 1} = {\boldsymbol{u}}_k^{m + 1} + \tilde {\boldsymbol{e}}_k^m $
    　　20 post-smooth:
    　　21　$ {\boldsymbol{u}}_k^{m + 1} = {{\boldsymbol{S}}^{{\mu _2}}}({{\boldsymbol{A}}_k},{{\boldsymbol{f}}_k},{\boldsymbol{u}}_k^{m + 1}) $。
    下载: 导出CSV
    　　算法2
    　　input:
    　　　symmetric positive definite coefficient matrix $ {\boldsymbol{A}} $;
    　　　right hand side vector $ {\boldsymbol{b}} $;
    　　　initiate a guess for $ {{\boldsymbol{x}}_0} \in {\mathbb{R}^n} $
    　　output:
    　　　solution $ {\boldsymbol{x}} $;
    　　1 compute the residual on the finest grid $ {{\boldsymbol{r}}_0} \mathop {\rm{ = }}\limits^\Delta {\boldsymbol{b}} - {\boldsymbol{A}}{{\boldsymbol{x}}_0} $;
    　　2 for $ i = 1,2,3, \cdots $ do
    　　3 solve $ {\boldsymbol{A}}{{\boldsymbol{z}}_0} = {{\boldsymbol{r}}_0} $ by the multigrid V-cycle algorithm
    　　4　initialize $ {{\boldsymbol{p}}_0} = {{\boldsymbol{z}}_0} $
    　　5　compute a step length the $ {\alpha _i} = ({\boldsymbol{r}}_{i - 1}^{\text{T}}{{\boldsymbol{z}}_{i - 1}})/({\boldsymbol{p}}_{i - 1}^{\text{T}}{\boldsymbol{A}}{{\boldsymbol{p}}_{i - 1}}) $
    　　 　 (where $ {{\boldsymbol{z}}_{i - 1}} $ can be solved by the multigrid V-cycle algorithm)
    　　6　update the approximate solution $ {{\boldsymbol{x}}_i} = {{\boldsymbol{x}}_{i - 1}} + {\alpha _i}{{\boldsymbol{p}}_{i - 1}} $
    　　7　update the residual $ {{\boldsymbol{r}}_i} = {{\boldsymbol{r}}_{i - 1}} - {\alpha _i}{\boldsymbol{A}}{{\boldsymbol{p}}_{i - 1}} $
    　　8　solve $ {\boldsymbol{A}}{{\boldsymbol{z}}_i} = {{\boldsymbol{r}}_i} $ by the multigrid V-cycle algorithm
    　　9　compute a gradient correction factor $ {\beta _i} = ({\boldsymbol{r}}_i^{\rm{T}}{{\boldsymbol{z}}_i})/({\boldsymbol{r}}_{i - 1}^{\rm{T}}{{\boldsymbol{z}}_{i - 1}}) $
    　　10　set the new search direction $ {{\boldsymbol{p}}_i} = {{\boldsymbol{z}}_i} + {\beta _i}{{\boldsymbol{p}}_{i - 1}} $
    　　11 end
    　　12 until convergence。
    下载: 导出CSV

    表  1  两种结构的拓扑优化迭代次数和柔度值

    Table  1.   Total numbers of iterations and compliance for topology optimization of 2 structures

    structureprecision PCG algorithmMG algorithmMGCG algorithm
    compliance Ctime T/siterations Icompliance Ctime T/siterations Icompliance Ctime T/siterations I
    cantilever beam1E−2186.885 287170187.306 996167184.738 364114
    1E−3186.883 482155186.939 3100160184.414 754100
    1E−4186.884 284155186.939 2105160184.014 075133
    1E−5186.883 984155186.939 1111160185.215 564113
    Michell structure1E−217.125 233138716.909 5707816.923 86878
    1E−317.123 423127117.119 926027716.923 96778
    1E−417.123 123027117.119 926927716.952 68093
    1E−517.123 123127117.074 726726716.943 66271
    下载: 导出CSV

    表  2  两种结构的拓扑优化迭代次数和柔度值(不同过滤半径)

    Table  2.   Total numbers of iterations and compliance for topology optimization of 2 structures (different filter radius)

    structurefilter radius RCG algorithmMG algorithmMGCG algorithm
    compliance Ctime T/siterations Icompliance Ctime T/siterations Icompliance Ctime T/siterations I
    cantilever beam3180.759 1505917182.218 3147225182.352 058103
    6186.884 282155186.939 2101160184.014 073133
    9192.522 998183193.053 764100193.096 75397
    12197.418 773136197.669 567105197.641 962115
    Michell structure316.738 584092616.864 516516816.894 18090
    617.123 123927117.119 928027716.952 68493
    917.382 911713317.322 611511717.168 65967
    1217.696 49911417.696 510711417.686 4100114
    下载: 导出CSV

    表  3  两种结构的拓扑优化迭代次数和柔度值(不同网格数量)

    Table  3.   Total numbers of iterations and compliance for topology optimization of 2 structures (different grids)

    structuregrid sizeCG algorithmMG algorithmMGCG algorithm
    compliance Ctime T/siterations Icompliance Ctime T/siterations Icompliance Ctime T/siterations I
    cantilever beam240 × 80186.484 016133186.315 319133183.148 41296
    480 × 160186.884 282155186.939 2101160184.014 073133
    960 × 320188.099 5430167188.513 5495173185.695 5419155
    Michell structure240 × 12015.995 44223315.867 0209115.816 01583
    480 × 24017.123 123527117.119 927427716.952 68093
    960 × 48018.293 9131929218.288 6148 130618.095 936191
    下载: 导出CSV
    • 参考文献(20)
  • [1] DEATON J D, GRANDHI R V. A survey of structural and multidisciplinary continuum topology optimization: post 2000[J]. Structural and Multidisciplinary Optimization, 2014, 49: 1-38. doi: 10.1007/s00158-013-0956-z
    [2] ZHU J H, ZHOU H, WANG C, et al. A review of topology optimization for additive manufacturing: status and challenges[J]. Chinese Journal of Aeronautics, 2021, 34(1): 91-110. doi: 10.1016/j.cja.2020.09.020
    [3] 卫志军, 申利敏, 关晖, 等. 拓扑优化技术在抑制流体晃荡中的数值模拟研究[J]. 应用数学和力学, 2021, 42(1): 49-57

    WEI Zhijun, SHEN Limin, GUAN Hui, et al. Numerical simulation of topology optimization technique fortank sloshing suppression[J]. Applied Mathematics and Mechanics, 2021, 42(1): 49-57.(in Chinese)
    [4] SIGMUND O. Design of material structures using topology optimization[D]. PhD Thesis. Copenhagen : Technical University of Denmark, 1994.
    [5] BRUNS T E, TORTORELLI D A. Topology optimization of non-linear elastic structures and compliant mechanisms[J]. Computer Methods in Applied Mechanics & Engineering, 2001, 190(26/27): 3443-3459.
    [6] 吴一帆, 郑百林, 何旅洋, 等. 结构拓扑优化变密度法的灰度单元等效转换方法[J]. 计算机辅助设计与图形学学报, 2017, 29(4): 759-767

    WU Yifan, ZHENG Bailin, HE Lüyang, et al. Equivalent conversion method of gray-scale elements for SIMP in structures[J]. Journal of Computer-Aided Design & Computer Graphics, 2017, 29(4): 759-767.(in Chinese)
    [7] XING J, QIE L. A simple way to achieve black-and-white designs in topology optimization[J]. Journal of Physics Conference Series, 2021, 1798(1): 012043. doi: 10.1088/1742-6596/1798/1/012043
    [8] 廉睿超, 敬石开, 何志军, 等. 拓扑优化变密度法的灰度单元分层双重惩罚方法[J]. 计算机辅助设计与图形学学报, 2020, 32(8): 1349-2583

    LIAN Ruichao, JING Shikai, HE Zhijun, et al. A hierarchical double penalty method of gray-scale elements for SIMP in Topology optimization[J]. Journal of Computer-Aided Design & Computer Graphics, 2020, 32(8): 1349-2583.(in Chinese)
    [9] 雷阳, 封建湖. 基于参数化水平集法的材料非线性子结构拓扑优化[J]. 应用数学和力学, 2021, 42(11): 1150-1160

    LEI Yang, FENG Jianhu. Topology optimization of nonlinear material structures based on parameterized level set and substructure methods[J]. Applied Mathematics and Mechanics, 2021, 42(11): 1150-1160.(in Chinese)
    [10] LIU B, WU Y, GUO J, et al. Finite difference Jacobian based Newton-Krylov coupling method for solving multi-physics nonlinear system of nuclear reactor[J]. Annals of Nuclear Energy, 2020, 148: 107670. doi: 10.1016/j.anucene.2020.107670
    [11] 肖文可, 陈星玎. 求解PageRank问题的重启GMRES修正的多分裂迭代法[J]. 应用数学和力学, 2022, 43(3): 330-340

    XIAO Wenke, CHEN Xingding. A modified multi-splitting iterative method with the restarted GMRES to solve the PageRank problem[J]. Applied Mathematics and Mechanics, 2022, 43(3): 330-340.(in Chinese)
    [12] LIU T. A nonlinear multigrid method for inverse problem in the multiphase porous media flow[J]. Applied Mathematics and Computation, 2018, 320(C): 271-281.
    [13] LIU T. A multigrid-homotopy method for nonlinear inverse problems[J]. Computers & Mathematics With Applications, 2020, 79(6): 1706-1717.
    [14] DENG L J, HUANG T Z, ZHAO X L, et al. Signal restoration combining Tikhonov regularization and multilevel method with thresholding strategy[J]. Journal of the Optical Society of America A: Optics, Image Science, and Vision, 2013, 30(5): 948-955. doi: 10.1364/JOSAA.30.000948
    [15] WAN F, YIN Y, ZHANG Q, et al. Analysis of parallel multigrid methods in real-time fluid simulation[J]. International Journal of Modeling, Simulation, and Scientific Computing, 2017, 8(4): 121-128.
    [16] AMIR O, AAGE N, LAZAROV B S. On multigrid-CG for efficient topology optimization[J]. Structural and Multidisciplinary Optimization, 2014, 49(5): 815-829. doi: 10.1007/s00158-013-1015-5
    [17] LIAO Z, ZHANG Y, WANG Y, et al. A triple acceleration method for topology optimization[J]. Structural and Multidisciplinary Optimization, 2019, 60(2): 727-744. doi: 10.1007/s00158-019-02234-6
    [18] LAZAROV B S, SIGMUND O. Filters in topology optimization based on Helmholtz-type differential equations[J]. International Journal for Numerical Methods in Engineering, 2011, 86(6): 765-781. doi: 10.1002/nme.3072
    [19] LAZAROV B S, SIGMUND O. Sensitivity filters in topology optimisation as a solution to Helmholtz type differential equation[C]// Proceedings of the Eighth World Congress on Structural and Multidisciplinary Optimization. Lisbon, Portugal, 2009.
    [20] ANDREASSEN E, CLAUSEN A, SCHEVENELS M, et al. Efficient topology optimization in MATLAB using 88 lines of code[J]. Structural and Multidisciplinary Optimization, 2011, 43(1): 1-16. doi: 10.1007/s00158-010-0594-7
    • 相关文章
    • 施引文献
  • 资源附件(0)
  • 加载中
图(15) / 表(5)
计量
  • 文章访问数:  780
  • HTML全文浏览量:  763
  • PDF下载量:  146
  • 被引次数: 0
出版历程
  • 收稿日期:  2022-03-01
  • 修回日期:  2022-12-15
  • 网络出版日期:  2023-01-06
  • 刊出日期:  2023-01-01

目录

摘要
HTML全文

  • 引　　言

  • 1.   优 化 模 型

  • 2.   拓扑优化算法

  • 2.1   总体框架

  • 2.2   PDE灵敏度滤波器

  • 2.3   PDE灵敏度滤波器的求解

  • 2.4   设计变量的更新

  • 2.5   收敛准则

  • 3.   数 值 算 例

  • 3.1   精度的影响

  • 3.2   过滤半径的影响

  • 3.3   网格数量的影响

  • 4.   结　　论

参考文献(20)
相关文章
施引文献
资源附件(0)

/

返回文章
返回
联系我们

版权所有 © 《应用数学和力学》编辑部

地址：重庆市南岸区学府大道66号 重庆交通大学90号信箱邮编：400074

电话/传真: 023-62652450　

电子邮箱: amm1980@vip.163.com； amm1980@163.com； applmathmech@cqjtu.edu.cn

期刊在线
邮件订阅 RSS

工信部备案号：渝ICP备11007697号-3

公安备案号：渝公网安备50010802005915号

本系统由 北京仁和汇智信息技术有限公司 开发