CN109840348B

CN109840348B - 一种三重加速的拓扑优化方法

Info

Publication number: CN109840348B
Application number: CN201811537582.1A
Authority: CN
Inventors: 王英俊; 廖中源
Original assignee: South China University of Technology SCUT
Current assignee: South China University of Technology SCUT
Priority date: 2018-12-15
Filing date: 2018-12-15
Publication date: 2023-06-20
Anticipated expiration: 2038-12-15
Also published as: CN109840348A

Abstract

本发明公开了一种三重加速的拓扑优化方法。为了加速拓扑优化的进程，该方法从多重网格密度映射、基于预处理与初值的共轭梯度法、局部更新三个方面进行改进，多重网格密度映射将设计域划分为几层粗细不同的网格，先在粗网格层上进行拓扑求解，将其得到的优化结果映射到下一层细网格作为初始值，使得细网格层上的拓扑优化进程跳过冗长的迭代，加快拓扑优化的进程；基于预处理和初值的共轭梯度法是在拓扑优化方程的迭代求解中加入预处理调节器以及初值，加快了求解速度；局部更新是选择更新单元，减少迭代时单元更新的数量来达到加速效果。本发明基于上述方法，实现了拓扑优化的三重加速，显著加速拓扑优化进程，获得高精度低计算成本的优化结果。

Description

一种三重加速的拓扑优化方法

技术领域

本发明属于拓扑优化技术领域，特别涉及一种三重加速的拓扑优化方法，

技术背景

拓扑优化是一种重要的结构设计方法，它是通过给定的约束条件，在设计域内寻找最优的结构拓扑形状。它是一个新兴的交叉领域，涉及数学、力学和计算机科学等学科，对制造行业和建筑行业，特别是航空航天，汽车以及能源等行业都有着重要的影响。随着计算能力的飞速提升，拓扑优化也得到了蓬勃发展。

随着工程技术的不断发展，对拓扑优化的要求也越来越高。其中一个难题就是大规模的拓扑优化问题。对于大规模问题，传统的优化方法效率低下，迭代次数和计算精度往往不能满足现实的优化要求。为了快速获得精度较高的优化结果，必须对拓扑优化的过程进行加速。

因此，三重加速的拓扑优化方法可以大大加快拓扑优化的进程，减少迭代次数以及计算时间，同时获得精度较高的优化结果。

发明内容

本发明的目的在于提供一种三重加速的拓扑优化方法，以解决传统拓扑优化方法在大规模问题上的低效率问题。其包含了三重加速：分别是多重网格密度映射，基于预处理和初值的共轭梯度法，局部更新方法。

本发明通过如下技术方案实现。

一种三重加速的拓扑优化方法，其包括多重网格密度映射、基于预处理与初值的共轭梯度法、局部更新三个方面；多重网格密度映射将原本的设计域划分为几层粗细不同的网格，先在粗网格层上进行拓扑求解，将其得到的优化结果映射到下一层细网格作为初始值，从而使得细网格层上的拓扑优化进程跳过冗长的迭代，加快拓扑优化的进程；基于预处理和初值的共轭梯度法是在拓扑优化方程的迭代求解中加入预处理调节器以及初始值，从而加快求解速度；局部更新则是选择性地更新单元，减少迭代时单元更新的数量达到加速效果。具体实现过程为以下步骤：

(1)定义拓扑优化的基本参数。

(2)定义初始网格和层数i，循环迭代i次。

(3)在每层网格中决定应用局部更新方法的时机。

(4)求解节点和单元的映射关系。

(5)将单元密度转化为节点密度。

(6)利用节点密度确定更新范围。

(7)将小密度单元的单元密度转化为0。

(8)在位移矩阵的迭代求解中应用基于预处理调节器和初值的共轭梯度法进行求解。

(9)在更新范围内更新单元密度。

(10)判断是否收敛。若不收敛，跳到步骤(5)；若收敛则跳出。

(11)判断是否达到多重网格的层数。若不达到，跳到步骤(12)。若达到，输出数据。

(12)在多重网格的两层之间，需要细分单元，由粗变细。

(13)在多重网格的两层之间，粗单元的优化结果将被下一层细网格继承，并转至步骤(3)。

进一步的，步骤(1)的网格在所有步骤的网格中最粗，网格单元数量也最少

进一步的，步骤(3)所述的局部更新时机与两次迭代之间单元密度的相对差异有关；根据两次迭代之间的单元密度差异收敛性和设定的启动参数，当两次迭代之间的单元密度差异收敛到某个数值时，启动局部更新方法。

进一步的，步骤(4)所述的节点与单元的映射关系指的是每个单元包含的节点以及每个节点周围的单元。

进一步的，步骤(5)所述的单元密度转化为节点密度是将节点周围的单元密度加权平均得到，节点周围的单元映射通过步骤(4)得到。

进一步的，步骤(6)所述的确定更新范围通过如下步骤实现：(1)定义节点密度的临界值范围，并把临界值范围与所有节点密度比较；当节点密度属于该临界值范围时，将节点周围的单元纳入更新范围；(2)定义扩散圈数；更新单元可在原范围的基础上，扩散圈数的扩大相应圈数；实现如下：首先将求得旧更新范围内所有节点的集合，再将属于旧更新范围节点周围的单元纳入更新范围，这就实现了扩散一圈；重复上述过程，实现扩散圈数。

进一步的，步骤(7)所述的小密度单元的转化是通过定义一个阈值，当单元密度小于阈值时，该单元密度将直接转化为0；当单元密度大于该阈值时，单元密度转化为1；值得注意的是，该阈值不能与体积约束等约束条件冲突。

进一步的，步骤(8)使用的初值是上一次迭代求解得到的方程未知量的解，在初值附近寻找最优解。

进一步的，步骤(12)所述的细分网格是将原始网格加密即增加网格数。

进一步的，步骤(13)单元所述的优化结果包含但不限于单元密度，节点密度，单元位移和节点位移等；其中单元密度可直接继承上一层粗网格；节点位移等节点量为节点的优化结果，可通过形函数插值得到单元内部的节点量。

本发明对比传统拓扑优化方法主要的优点如下：

(1)求解时间短。

(2)对于大规模网格拓扑优化的加速效果非常好。网格规模越大，加速效果越明显。

(3)对于三维模型的拓扑优化加速效果更加显著。

附图说明

图1为三重加速的拓扑优化方法的流程图。

图2是多重网格方法的一种网格细化示意图。

图3a至图3c是不同网格规模下的传统拓扑优化方法结果对比图。

图3d至图3f是不同网格规模下的三重加速方法的最终优化结果对比图。

图4a是三维情况下三重加速的拓扑优化方法效果图。

图4b是三维情况下传统的拓扑优化方法效果图。

具体实施方式

下面结合附图和实施例对本发明做进一步的说明，但本发明要求保护的范围并不局限于实施例表述的范围。此外，本发明中的各个实施方式中涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

本发明可用于多种基于网格单元的拓扑优化方法，下面是以基于变密度法(SIMP)的三重加速拓扑优化方法为例。变密度法是将每个单元的弹性模量与其相对密度有着确定的函数关系，以单元相对密度作为设计变量。基于变密度法的最小柔度优化问题的基本模型如下：

Findx＝(x₁,x₂…x_N)^T

x_e是设计变量(单元相对密度)，其取值范围0～1；

x_min是最小相对密度(避免出现奇异)；

x是设计变量的向量；

N是设计变量的个数；

c是总体柔度；

E_e是单元弹性模量(胞元等效弹性模量)；

U是总体位移矩阵；

F是总体受力矩阵；

K是总体刚度矩阵；

k_e是单元刚度矩阵；

u_e是单元位移矩阵；

V(x),V₀是分别是实体体积和设计域总体积；

f是体积分数；

T表示转置矩阵，N表示单元ID编号。

基于经典的变密度法，本发明的三重加速的拓扑优化方法主要包括以下步骤：

(1)定义拓扑优化的基本参数；定义设计域，载荷，约束等边界条件，设定体积约束和目标函数如上面的模型所示

是目标函数，使得柔度最小。

是约束条件。其中实施例的体积约束是0.5。

(2)定义初始网格和层数i，循环迭代i次(i＝1,2,……)。

(3)在每层网格中决定应用局部更新方法的时机。设定启动参数ST,两次迭代之间单元密度的相对差异值为change，当change＜ST，此时拓扑优化已得到近似的结构，启动局部更新。启动参数ST可以使得局部更新具有一个合适的初始结构，加快求解速度，提高优化精度。

(4)求解节点和单元的映射关系。映射关系是指每个单元包含多少节点，以及每个节点周围有多少单元，可通过单元与节点的编号规律求得。

(5)将单元密度转化为节点密度，节点密度由节点周围的单元密度加权平均得到，节点周围的单元映射通过步骤(4)得到。

(6)利用节点密度确定更新范围。设定一个临界值ρ_c，当节点密度大于ρ_c时，将该节点周围的单元纳入更新范围。再定义扩散层数D_num。首先将更新域中所有节点密度赋1，然后判断节点密度是否等于1，将等于1的节点周围的单元纳入更新范围，这样就实现了更新域向外扩散一层。重复D_num次，扩散D_num层。

(7)将小密度单元的单元密度转化为0。设定单元密度阈值，当单元密度小于该阈值时，将该单元密度转化为0。阈值不能大于体积约束。小密度单元置零可以减少灰度单元，提高收敛速度。

(8)对于有限元方程KU＝F，K是刚度矩阵，U是位移矩阵，F是外部载荷，在求解位移矩阵U中应用基于预处理调节器和初值的共轭梯度法进行求解。在每次的位移矩阵求解中，上一次求得的位移矩阵U_pre会作为初始值应用到本次求解中，在初始值附近寻找最优解，减少迭代次数减少迭代次数以及计算时间。

基于预处理调节器和初值的共轭梯度法的迭代求解过程如下所示：

输入：K,M,F,U_pre

步骤一：初始化

U₀＝U_pre％％U_pre表示前一次优化求解结果

r₀＝F-KU₀,z₀＝M^-1r₀

步骤二：循环计算

While r_i>r_t

U_i＝U_i-1+α_ip_i；r_i＝r_i-1-α_iw

i＝i+1.

End while

输出：U_i；

其中U₀为给定初值(在该算法中是前一次优化的结果)；M为预处理矩阵；i为迭代次数；p_i为第i次迭代方向向量，α_i，β_i，ρ_i，z_i，w是中间变量，r_i是残差。当r_i小于r_t时，可以认为优化结果U_i已达到需要的精度，迭代停止。在算法中，下标0,1,2…i表示迭代次数，粗体加斜体的符号表示矩阵或向量，M^-1表示M的逆矩阵，

中的T表示转置矩阵。

(9)在更新范围内更新单元密度。

(10)判断是否收敛。若不收敛，跳到步骤(5)；若收敛则跳出。收敛准则可根据实际需要定义，最常用的是迭代前后的最大单元密度变化量，当密度变化小于阈值，则收敛。

(11)判断是否达到多重网格的最细网格层。根据步骤(2)预先设定的具体层数进行判断，循环次数若等于i，输出数据。否则跳到步骤(12)。所输出的数据包括：输出目标函数值，体积分数，密度分布以及最终的拓扑结构等等。

(12)在多重网格的两层之间，需要细分单元，由粗变细。如图2所示，对于二维网格，是1个网格细化为4个；对于三维网格，是1个网格细化为8个。

(13)在多重网格的两层之间，粗单元的优化单元密度以及节点位移将被下一层细网格继承。由于细网格是由粗网格细分出来，故其单元密度直接等于原始单元，拓扑优化结果映射到下一层细网格作为初始值。关于节点的优化结果是通过形函数插值得到单元内部的节点值。而节点位移可由如下公式插值得到：

表示细化后的网格的位移矩阵，其中n是该网格的所属的层数，i表示下层节点ID,k为上层网格节点数。N_j(X_i)表示上层网格j号节点密度的形函数在下层网格i节点坐标处的数值X_i表示下层节点i的坐标。j表示上层节点ID，/>

表示上层网格的位移矩阵。形函数可以通过多种方法构造，最常用的是拉格朗日插值法。对于一个粗网格单元，包含4个插值节点(2D)或8个插值节点(3D)，单元内每一点的数值可通过粗网格插值节点计算得到，上层粗网格j号插值节点密度的形函数可由如下公式得到：x_i，y_i，z_i分别表示除了节点j其他插值节点在x，y，z方向的坐标值，x_j，y_j，z_j表示节点j在x，y，z方向的坐标值。x,y,z是单元内部任意一点的坐标值。在2D单元中，z方向的坐标值全部为0。

图3a至图3c分别为40×20、80×40、160×80网格规模下的传统拓扑优化方法，图3d至图3f分别为40×20、80×40、160×80网格规模下的三重加速方法的最终优化结果，图4a和图4b分别是三维情况下传统优化方法和三重加速的优化结果。可以看到，三重加速有着黑白分明的结构，灰度单元更少。而且三重加速所耗时间更少，能在更短时间内获得精度较高的拓扑优化结构。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所做的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims

1.一种三重加速的拓扑优化方法，其特征在于，包括多重网格密度映射、基于预处理与初值的共轭梯度法、局部更新三个方面；多重网格密度映射将原本的设计域划分为几层粗细不同的网格，先在粗网格层上进行拓扑求解，将其得到的优化结果映射到下一层细网格作为初始值，从而使得细网格层上的拓扑优化进程跳过冗长的迭代，加快拓扑优化的进程；基于预处理和初值的共轭梯度法是在拓扑优化方程的迭代求解中加入预处理调节器以及初始值，从而加快求解速度；局部更新则是选择性地更新单元，减少迭代时单元更新的数量达到加速效果；具体实现过程为以下步骤：

(1)基于变密度法的最小柔度优化问题的基本模型，定义拓扑优化的基本参数，基于变密度法的最小柔度优化问题的基本模型如下：

Find x＝(x₁,x₂…x_N)^T

x_e是设计变量，其取值范围0～1；

x_min是最小相对密度；

x是设计变量的向量；

N是设计变量的个数；

c是总体柔度；

E_e是单元弹性模量；

U是总体位移矩阵；

F是总体受力矩阵；

K是总体刚度矩阵；

k_e是单元刚度矩阵；

u_e是单元位移矩阵；

V(x),V₀是分别是实体体积和设计域总体积；

f是体积分数；

T表示转置矩阵，x₁,x₂…x_N表示单元ID编号为1～N的设计变量；

定义设计域，载荷，约束边界条件，设定体积约束和目标函数：

是目标函数，使得柔度最小；

是约束条件；

(2)定义初始网格和层数i，循环迭代i次；

(3)在每层网格中决定应用局部更新时机：设定启动参数ST,两次迭代之间单元密度的相对差异值为change，当change＜ST，此时拓扑优化已得到近似的结构，启动局部更新；

(4)求解节点和单元的映射关系；

(5)将单元密度转化为节点密度；

(6)利用节点密度确定更新范围：设定一个临界值ρ_c，当节点密度大于ρ_c时，将该节点周围的单元纳入更新范围；再定义扩散层数D_num；首先将更新域中所有节点密度赋1，然后判断节点密度是否等于1，将等于1的节点周围的单元纳入更新范围，这样就实现了更新域向外扩散一层；重复D_num次，扩散D_num层；

(7)将小密度单元的单元密度转化为0：当单元密度小于设置的阈值时，将该单元密度转化为0；

(8)在位移矩阵的迭代求解中应用基于预处理调节器和初值的共轭梯度法进行求解：对于有限元方程KU＝F，K是刚度矩阵，U是位移矩阵，F是外部载荷，在求解位移矩阵U中应用基于预处理调节器和初值的共轭梯度法进行求解；

(9)在更新范围内更新单元密度；

(10)判断是否收敛；若不收敛，跳到步骤(5)；若收敛则跳出；

(11)判断是否达到多重网格的层数；若没达到，跳到步骤(12)；若达到，输出数据；

(12)在多重网格的两层之间，需要细分单元，由粗变细；

2.根据权利要求1所述的一种三重加速的拓扑优化方法，其特征在于，步骤(1)的网格在所有步骤的网格中最粗，网格单元数量也最少。

3.根据权利要求1所述的一种三重加速的拓扑优化方法，其特征在于，步骤(3)所述的局部更新时机与两次迭代之间单元密度的相对差异有关；根据两次迭代之间的单元密度差异收敛性和设定的启动参数，启动局部更新方法。

4.根据权利要求1所述的一种三重加速的拓扑优化方法，其特征在于，步骤(4)所述的节点与单元的映射关系指的是每个单元包含的节点以及每个节点周围的单元。

5.根据权利要求1所述的一种三重加速的拓扑优化方法，其特征在于，步骤(5)所述的节点密度是将节点周围的单元密度加权平均得到，节点周围的单元映射通过步骤(4)得到。

6.根据权利要求1所述的一种三重加速的拓扑优化方法，其特征在于，步骤(6)所述的确定更新范围通过如下步骤实现：1)定义节点密度的临界值范围，并把临界值范围与所有节点密度比较；当节点密度属于该临界值范围时，将节点周围的单元纳入更新范围；2)定义扩散圈数；更新单元在原范围的基础上，扩大相应圈数；扩散圈数的实现如下：首先将求得旧更新范围内所有节点的集合，再将属于旧更新范围节点周围的单元纳入更新范围，即扩散了一圈；重复上述过程，实现扩散圈数。

7.根据权利要求1一种三重加速的拓扑优化方法，其特征在于，步骤(7)所述的小密度单元的单元密度转化是通过定义一个阈值，当单元密度小于阈值时，该单元密度将直接转化为0；当单元密度大于该阈值时，单元密度转化为1；该阈值不能与约束条件冲突。

8.根据权利要求1一种三重加速的拓扑优化方法，其特征在于，步骤(8)的初值是上一次迭代求解得到的方程未知量的解，在初值附近寻找最优解。

9.根据权利要求1一种三重加速的拓扑优化方法，其特征在于，步骤(12)所述的细分网格是将原始网格加密即增加网格数。

10.根据权利要求1一种三重加速的拓扑优化方法，其特征在于，步骤(13)所述的粗单元优化结果包含单元密度，节点密度，单元位移和节点位移；其中单元密度直接继承上一层粗网格；节点位移为节点的优化结果，通过插值算法得到单元内部的节点位移。