CN112100882A - 一种具有光滑边界的连续体结构密度演化拓扑优化方法 - Google Patents

Abstract

本发明公开的一种具有光滑边界的连续体结构密度演化拓扑优化方法，包括以下步骤：建立连续体结构的几何初始设计域；基于无惩罚函数的材料插值模型建立几何初始设计域的有限元模型；建立连续体结构拓扑优化问题数学模型；计算优化目标和约束条件对设计变量的灵敏度，并过滤；采用优化算法根据灵敏度更新设计变量，并对密度进行分级过滤，达到收敛和约束条件，得到优化结构的有限元模型；构建节点应变能水平集，对边界单元进行处理，得到基于节点密度或节点应变能水平集显示的光滑边界的结构拓扑优化模型；本发明提出无惩罚函数的材料插值模型和密度分级过滤方法，对结构边界单元进行处理，得到了带有光滑边界的优化结果，高效的解决拓扑优化问题。

Description

一种具有光滑边界的连续体结构密度演化拓扑优化方法

技术领域

本发明涉及结构优化的研究领域，特别涉及一种具有光滑边界的连续体结构密度演化拓扑优化方法。

背景技术

结构拓扑优化是想办法把有限材料合理地分布到设计域上使结构达到最优的性能。结构拓扑优化设计方法是新型的数字化结构设计方式，通过建立包含目标函数和约束方程的数学模型，对结构进行有限元数值分析，在设计域中按照优化准则或者数学规划方法迭代求解出满足目标要求的材料分布。

结构拓扑优化设计方法最开始只是局限于机械结构设计问题，现在已经逐步扩展到物理学，流体，声学，电磁学，光学和热传导等结构设计问题。

和Kikuchi早在1988年提出了拓扑优化的均匀化方法，随后几十年，在均匀化方法的基础上发展了很多不同的方法，包括密度法(SIMP)，水平集法(LSM)和渐进演化法(ESO)，形状导数法和相场法等等。以上这些方法都是基于有限元理论，把设计域离散为若干单元，然后根据优化程序来确定在结构性能最优的情况下，哪些单元是实体材料单元(单元密度为1)，哪些单元是空心单元(单元密度为0)或中间单元(单元密度介于0和1之间)。

迄今为止，拓扑优化的理论已经相当成熟。在结构拓扑优化设计领域，密度法和演化法以单元密度或者节点密度作为设计变量，水平集方法则以结构形状倒数来确定结构的边界或者形状。

目前，常用的拓扑优化方法是SIMP方法(固体各项同性材料惩罚函数法)，为了保证算法的稳定性，惩罚函数通常为大于等于2的整数，SIMP法优化过程中材料相对密度介于0～1之间。变密度法相对于其他方法的最大优势是：原理简单，且对任意形状的可行设计域具有普适性。此外，渐进演化方法具有算法简单，与有限元分析程序连接容易等优点，最终得到的优化结果的有限元模型只有0/1的单元；水平集方法可以得到的优化结果则具有光滑的几何边界，因为惩罚函数没有明确的物理含义，在优化问题与刚度和质量同时相联系的时候，会严重影响优化算法的稳定。

发明内容

本发明的目的在于克服现有技术的缺点与不足，提供一种具有光滑边界的连续体结构密度演化拓扑优化方法，提出了无惩罚函数的材料插值模型和密度分级过滤方法，通过该方法可以得到只有0/1单元的最优拓扑结构的有限元模型。然后采用节点应变能零水平集的思想，对结构边界单元进行处理，最后得到了带有光滑边界的优化结果；其目的高效的解决结构优化问题，该方法既延续了变密度法的优点，又综合具有渐进演化方法和水平集方法的优点，可以实现三种方法取长补短的完美结合。

本发明的目的通过以下的技术方案实现：

一种具有光滑边界的连续体结构密度演化拓扑优化方法，其特征在于，包括以下步骤：

S1、通过画图建模模块建立连续体结构的几何初始设计域；

S2、基于无惩罚函数的材料插值模型建立几何初始设计域的有限元模型；把结构划分为若干个有限元单元网格，通过单元密度的有无确定优化的结果；

S3、建立连续体结构拓扑优化问题数学模型，以结构应变能最小为目标，所用材料体积为约束建立优化的数学模型；

S4、计算优化目标和约束条件对设计变量的灵敏度，并对灵敏度进行过滤；

S5、采用优化算法根据灵敏度更新设计变量，并对密度进行分级过滤，直至达到收敛和约束条件，得到优化结构的有限元模型；

S6、构建节点应变能水平集，对边界单元进行处理，得到基于节点谜底或节点应变能水平集显示的光滑边界的结构拓扑优化模型，并保存。

进一步地，所述步骤S1具体为：画图建模模块根据连续体结构的结构形状画出连续体结构的初始几何实体图形，进而建立几何初始设计域，即采用画图建模软件，画出结构的初始几何实体图形。

进一步地，所述步骤S2具体为：将几何初始设计域划分为N个有限元单元，具体的单元类型根据优化的问题来定义，指定每个对应有限元单元的相对密度为x_e，按照有限元的方法生成单元刚度矩阵，再组装成全局刚度矩阵，形成初始设计域的有限元模型。有限元模型中用到的材料的弹性模量的插值模型为：

E(x_e)＝E_min+x_e(E₀-E_min)，

其中，E₀为材料的实际弹性模量；E_min是为了防止单元刚度矩阵奇异所给定的一个无穷小量，物理含义为当单元密度为0时材料的弹性模量；x_e为第e个有限元单元密度。

进一步地，所述步骤S3具体为：设置静力问题以结构应变能c最小作为目标函数，体积比f为约束条件，同时需要满足静力平衡方程，其优化问题数学模型如下：

其中，c为结构应变能，U为全局位移，F为负荷向量，K为全局刚度矩阵，u_e为单元位移向量，k₀为单元刚度矩阵，x_e为设计变量，即有限元单元密度；x为有限元单元密度向量，N为几何初始设计域采用有限元离散的单元总数，v(x)为优化结构可用材料体积，v₀为几何初始设计域的总体积，f为体积分数。

进一步地，所述步骤S4具体如下：

先采用伴随矩阵法求解目标函数对有限元单元密度x_e的灵敏度：

其中，c为结构应变能，E₀为材料的实际弹性模量；E_min为无穷小量，物理含义为当单元密度为0时材料的弹性模量；x_e为第e个有限元单元密度，也为设计变量；u_e为单元位移向量，k₀为单元刚度矩阵；

再计算约束条件对有限元单元密度x_e的灵敏度：

其中，v为结构总体积，x_e为第e个有限元单元的密度；

为了避免棋盘格现象，对敏感度进行过滤：

其中，n_e为在过滤半径r_min内与有限元单元e中心点的直线距离的所有有限元单元集合，γ为单元虚密度值的下限值；H_ei为权重系数，表示在过滤半径r_min内各个有限元单元所占的权重；

权重系数定义如下：

H_ei＝max(0,r_min-Δ(e,i))，

其中，Δ(e,i)为有限元单元e中心点到有限元单元i中心点的直线距离。这样处理的目的一方面是为了避免棋盘格现象。另一方面，通过调整过滤半径，也可以得到细节尺寸不一样的拓扑结构，进而满足工艺生产需求

进一步地，所述步骤S5具体为：采用密度分级过滤方法进行更新设计变量，具体的优化算法可以采用OC算法，MMA算法和GCMMA算法来更新设计变量，然后对密度进行过滤。

传统的密度法中，为了更快的达到收敛条件，一般采用惩罚函数来实现。即使这样，最后还会存在很多的中间密度单元,而且需要迭代更多的步数才能收敛。本实施例提出一种密度分级过滤的方法，可以完美解决灰度单元的问题，并且快速收敛，这样得到的密度单元不是0就是1。计算如下：

其中，a₀为中间参数，G_f为密度过滤分级因子，

为更新的设计变量，loop为当前循环步数，a₁为当前循环步数中设计变量小于a0(loop)对应的单元序号；

依次步骤4至步骤5循环，直至达到收敛和约束条件，得到只包含0/1单元的优化结构的有限元模型，这一步得到的优化结果类似于传统的SIMP法和ESO方法得到的优化结果。

进一步地，所述步骤S6具体为：

构建节点应变能水平集，对边界单元进行处理，最终得到基于节点密度或节点应变能零水平集显示的光滑边界的结构拓扑优化模型，并保存为stl格式的3维模型文件。具体步骤分为三步：1.单元细分；2.构造节点应变能的水平集；3.优化结果零水平集显示和存储。

S601、单元细分；单元细分存在以下两种情况：需要得到高分辨率的优化结果，则需要过滤半径小于1.5；对于三维结构，其优化过程计算量大；两种情况都需要进行单元细分；

细分单元不是必须的，在以下两种情况下需要做单元细分，一是如果需要得到高分辨率的优化结果，此时过滤半径很小，必须基于0/1单元优化结果的有限元模型进行一次细分再进行边界单元处理。二是对于三维结构，优化过程的计算量非常大，可以先划分粗网格得到初步的优化结果的有限元模型，然后再进行单元细分，得到更加精细的优化结果。对于平面四节点矩形单元而言，细分前后单元密度转换如下：

其中x_i表示第i个单元细分前的密度，细分为4个单元后，对应的密度分别为x_i1,x_i2,x_i3,x_i4；原单元节点编号分别为1,2,3,4，细分后节点编号为1,2,...,9，如图2所示；

细分后节点的位移可以重构有限元计算(只需计算一次)，为了节省计算时间，也可以基于细分前的粗网格节点的位移通过线性插值直接得到，计算公式如下：

其中u_i(i＝1,2,3,4)表示细分前单元节点的位移，

表示细分后单元节点的位移。对于三维问题，可采用类似的方法得到细分单元的密度和节点位移。

S602、构造节点应变能的水平集；单元的应变能为：

其中，x_e为第e个有限元单元的密度，u_i为第e个单元对应节点的位移，k₀为第e个单元的单元刚度矩阵；

节点应变能通过插值得到，对于二维平面问题，单元上任一点的应变能水平集函数表示为：

其中，Φ(o)为任一点o的水平集的值，

为i个单元的第j个节点上的应变能，N_j(o)为插值函数；

如果对应平面矩形4节点单元的问题，其表达式为：

N_j(o)＝(1+ξ_oξ_j)(1+η_oη_j)，

上式中(ξ_o,η_o)对应o点的节点坐标，(ξ_j,η_j)为j节点的坐标，具体如图3所示。对于三维问题，可用类似的方法构造8节点单元的基于节点应变能的水平集函数。

基于应变能水平集函数，则可根据不同的水平集值把设计区域分为实体区域Ω、空心区域D\Ω、边界区域Γ，区分如下：

在迭代过程中，为了得到光滑的边界，优化后的有限元模型分为实心单元、空心单元、边界单元；对单元密度的处理如下：

其中，

为i个单元的第j个节点上的应变能的值，S为水平集值；

对应边界单元，0＜x_e＜1，把边界单元细分成20×20的精细网格，每个单元的应变能可以根据四个角节点的应变能采用插值计算得到，统计单元应变能数值大于水平集值S的个数记为N_lb，单元总数记为N_mb，则边界单元的密度为N_lb/N_mb；

采用二分法确定水平集值S的值，水平集值的上界S_ub和下界S_lb的初始值分别对应节点应变能的最大值和最小值，则水平集值为：S＝(S_ub+S_lb)/2；

当单元密度更新后，为了满足体积约束，通过对比优化设计模型的体积V和材料约束体积V₀，水平集值的上界和下届也随着更新，其更新策略为：

水平集值的上界S_ub和下界S_lb在迭代的过程中不断更新，直到S_ub和S_lb的值S_ub-S_lb＜10^-9，最后结构优化设计模型的体积V就会达到约束值V₀；

S603、最终优化的结果通过Φ(o)-S＝0的节点应变能零水平集可直接显示，其物理含义是当边界节点上的应变能都等于S时，对应的零水平集的结构满足体积约束条件，此时结构应变能最小；优化结果也可通过把边界单元经过处理后的有限元模型的单元密度转化为节点密度，采用节点密度ψ(o)-0.5＝0的水平集显示；得到的结果具有显式的光滑边界，存储为stl格式的文件。

本发明与现有技术相比，具有如下优点和有益效果：

本发明基于材料插值模型，去掉了物理意义不明确的惩罚函数，使基于密度优化的优化算法更加稳定；提出了一种密度分级过滤方法，得到只有0/1单元的拓扑结构，使用该过滤方法可以显著减少迭代步数，提高优化算法的计算效率；采用节点应变能零水平集，把边界单元经过处理后，可以得到具有光滑边界的拓扑优化结构模型，可以直接3D打印。

附图说明

图1为本发明所述一种具有光滑边界的连续体结构密度演化拓扑优化方法的流程图；

图2为本发明所述实施例1中悬臂梁的初始的几何设计域、约束条件和荷载情况示意图；

图3为本发明所述实施例1中优化目标和约束函数的迭代收敛曲线和边界密度演化过程示意图；

图4为本发明所述实施例1中只包含0/1单元的优化结果的有限元模型示意图；

图5为本发明所述实施例1中进行边界处理时的单元细分示意图；

图6为本发明所述实施例1中进行边界处理时基于节点应变能构造水平集，自然坐标系下任意点和单元节点坐标示意图；

图7为本发明所述实施例1中边界单元经过处理后的优化结果的有限元模型示意图；

图8为本发明所述实施例1中边界单元处理过程中节点应变能的水平集示意图；

图9为本发明所述实施例1中边界单元经过处理后的有限元模型对应的节点应变能零水平集示意图；

图10为本发明所述实施例1中边界单元经过处理后有限元模型对应的节点密度在0-1之间变化的水平集示意图；

图11为本发明所述实施例1中边界单元经过处理后的有限元模型对应的节点密度ψ(o)-0.5＝0的水平集示意图；

图12为本发明所述实施例2中三维结构桥梁的初始几何设计域，约束条件和荷载情况示意图；

图13为本发明所述实施例2中桥梁优化过程中目标函数和约束条件的迭代收敛曲线示意图；

图14为本发明所述实施例2中桥梁优化结果的有限元模型示意图；

图15为本发明所述实施例2中桥梁结构的优化最终显示结果示意图。

具体实施方式

下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。

实施例1：

一种具有光滑边界的连续体结构密度演化拓扑优化方法，如图1所示，包括以下步骤：

S1、通过画图建模模块建立连续体结构的几何初始设计域；

S2、基于无惩罚函数的材料插值模型建立几何初始设计域的有限元模型，把结构划分为若干个有限元单元网格，通过单元密度的有无确定优化的结果；

S3、建立连续体结构拓扑优化问题数学模型，以结构应变能最小为目标，所用材料体积为约束建立优化的数学模型；

S4、计算优化目标和约束条件对设计变量的灵敏度，并对灵敏度进行过滤；

S5、采用优化算法根据灵敏度更新设计变量，并对密度进行分级过滤，直至达到收敛和约束条件，得到优化结构的有限元模型；

S6、构建节点应变能水平集，对边界单元进行处理，得到基于节点谜底或节点应变能水平集显示的光滑边界的结构拓扑优化模型，并保存。

具体如下：

步骤1：如图2所示悬臂梁，建立连续体结构拓扑优化几何模型。将几何尺寸定为设计域，右端中心荷载为F，优化目标为使其应变能最小，体积约束为0.5。但本发明的设计域不限于此。

步骤2：建立连续体结构拓扑优化问题有限元模型，建立不带惩罚函数的材料插值模型。每个有限单元的相对密度预先设为1。

步骤3：基于以上模型建立连续体结构拓扑优化问题数学模型，最小化结构应变能c作为目标函数，体积比f为约束条件，其优化数学模型如下：

subject to:v(x)/v₀＝f

KU＝F

0≤x≤1

式中c是结构应变能，U和F分别是全局位移和荷载向量，K为全局刚度矩阵，u_e是单元位移向量，k₀是单元刚度矩阵，x是设计变量向量(即单元密度向量)，N是设计域采用有限元离散的单元总数，v(x)和v₀分别是可用的材料体积和设计域的总体积，f是体积分数。

步骤4：优化问题敏度分析计算，采用伴随矩阵法求解目标函数和约束函数的灵敏度，并对灵敏度进行过滤。其目标函数敏度分析为：

约束函数敏度分析为：

步骤5：此实例采用优化准侧(OC)算法求解，更新设计变量x_i，然后对密度进行分级过滤，依次步骤4至步骤5循环，直至达到收敛条件，得到只包含0/1单元的优化结构的有限元模型。体积约束为0.5，过滤半径取2，敏度过滤分级因子取50，最后应变能为31.30。应变能，体积迭代收敛曲线和边界密度演化过程如图3所示，图中(a),(b),(c),(d)分别为第10步，第20步，第30步，第40步密度演化优化的结果。图4为密度演化收敛后的有限元模型，只包含0/1单元。

步骤6：基于优化结构的有限元模型，构建节点应变能水平集，对边界单元进行自动化处理，得到基于节点密度或节点应变能零水平集显示的光滑边界的结构拓扑优化模型。最后显示的优化结果具有显式的光滑边界，对于三维模型，可存储为stl格式的文件直接用于3D打印。

具体步骤分为三步：1.单元细分；2.构造节点应变能的水平集；3.优化结果零水平集显示和存储。

S601、单元细分；单元细分存在以下两种情况：需要得到高分辨率的优化结果，则需要过滤半径小于1.5；对于三维结构，其优化过程计算量大；两种情况都需要进行单元细分；

细分单元不是必须的，在以下两种情况下需要做单元细分，一是如果需要得到高分辨率的优化结果，此时过滤半径很小，必须基于0/1单元优化结果的有限元模型进行一次细分再进行边界单元处理。二是对于三维结构，优化过程的计算量非常大，可以先划分粗网格得到初步的优化结果的有限元模型，然后再进行单元细分，得到更加精细的优化结果。对于平面四节点矩形单元而言，细分前后单元密度转换如下：

其中x_i表示第i个单元细分前的密度，细分为4个单元后，对应的密度分别为x_i1,x_i2,x_i3,x_i4；原单元节点编号分别为1,2,3,4，细分后节点编号为1,2,...,9，如图5所示；

细分后节点的位移可以重构有限元计算(只需计算一次)，为了节省计算时间，也可以基于细分前的粗网格节点的位移通过线性插值直接得到，计算公式如下：

其中u_i(i＝1,2,3,4)表示细分前单元节点的位移，

表示细分后单元节点的位移。对于三维问题，可采用类似的方法得到细分单元的密度和节点位移。

S602、构造节点应变能的水平集；单元的应变能为：

其中，x_e为第e个有限元单元的密度，u_i为第e个单元对应节点的位移，k₀为第e个单元的单元刚度矩阵；

节点应变能通过插值得到，对于二维平面问题，单元上任一点的应变能水平集函数表示为：

其中，Φ(o)为任一点o的水平集的值，

为i个单元的第j个节点上的应变能，N_j(o)为插值函数；

如果对应平面矩形4节点单元的问题，其表达式为：

N_j(o)＝(1+ξ_oξ_j)(1+η_oη_j)，

上式中(ξ_o,η_o)对应o点的节点坐标，(ξ_j,η_j)为j节点的坐标，具体如图6所示。对于三维问题，可用类似的方法构造8节点单元的基于节点应变能的水平集函数。

基于应变能水平集函数，则可根据不同的水平集值把设计区域分为实体区域Ω、空心区域D\Ω、边界区域Γ，区分如下：

在迭代过程中，为了得到光滑的边界，优化后的有限元模型分为实心单元、空心单元、边界单元；对单元密度的处理如下：

其中，

为i个单元的第j个节点上的应变能的值，S为水平集值；

对应边界单元，0＜x_e＜1，把边界单元细分成20×20的精细网格，每个单元的应变能可以根据四个角节点的应变能采用插值计算得到，统计单元应变能数值大于水平集值S的个数记为N_lb，单元总数记为N_mb，则边界单元的密度为N_lb/N_mb；

采用二分法确定水平集值S的值，水平集值的上界S_ub和下界S_lb的初始值分别对应节点应变能的最大值和最小值，则水平集值为：S＝(S_ub+S_lb)/2；

当单元密度更新后，为了满足体积约束，通过对比优化设计模型的体积V和材料约束体积V₀，水平集值的上界和下届也随着更新，其更新策略为：

水平集值的上界S_ub和下界S_lb在迭代的过程中不断更新，直到S_ub和S_lb的值S_ub-S_lb＜10^-9，最后结构优化设计模型的体积V就会达到约束值V₀；

S603、最终优化的结果通过Φ(o)-S＝0的节点应变能零水平集可直接显示，其物理含义是当边界节点上的应变能都等于S时，对应的零水平集的结构满足体积约束条件，此时结构应变能最小；优化结果也可通过把边界单元经过处理后的有限元模型的单元密度转化为节点密度，采用节点密度ψ(o)-0.5＝0的水平集显示；得到的结果具有显式的光滑边界，存储为stl格式的文件。

图7是本发明方法实施例中的边界单元经过处理后的优化结果的有限元模型；图8是本发明方法实施例中边界单元处理过程中节点应变能的水平集；图9是本发明方法实施例中的边界单元经过处理后的有限元模型对应的节点应变能零水平集；图10是本发明方法实施例中的边界单元经过处理后有限元模型对应的节点密度在0-1之间变化的水平集；图11是本发明方法实施例中的边界单元经过处理后的有限元模型对应的节点密度为0.5对应的水平集。图9和图11的结果非常接近，但图11更加边界更加光滑，是因为节点密度等值线显示相当于把单元密度做了一次线性插值。实际设计中，可以根据不同的需求选择图9和图11的显示结果。

此实例中，为了验证本文提出的边界密度演化方法(BDE)方法的稳定性，分别考虑了过滤半径，密度过滤分级因子和有限元单元网格尺寸对优化结果的影响。对比结果由表1，表2和表3给出：

表1.过滤半径对优化算法稳定性的影响

表1可以看出，在体积分数和密度过滤分级因子不变的情况下，过滤半径对优化目标值和迭代步数影响不大，主要是局部细节尺寸的改变，也就是当过滤半径大小不同时，可以得到目标值相近但形态不同的优化结果。此外当3D打印有局部尺寸工艺限制时，可以通过改变过滤半径来实现。

表2.密度过滤分级因子对优化算法稳定性的影响。

在体积分数和过滤半径不变的情况下，密度过滤分级因子对优化目标值影响不大。但是密度过滤分级因子越小，迭代步数越少；密度过滤分级因子越大，迭代步数越多，同时计算耗时越长。与表1类似，当改变密度过滤分级因子时，也可以得到目标值相近但形态不同的优化结果。

表3.有限元单元网格尺寸对优化算法稳定性的影响。

在体积分数和过滤半径不变的情况下，有限元网格划分越多，也就是单元尺寸越小的情况下，密度过滤分级因子和有限元单元网格尺寸对优化目标值影响不大。但是如果需要得到高分辨率的优化结果，需要把有限元单元网格尽量划小，同时加大密度过滤分级因子。

此实例中，为了进一步验证本文提出的边界密度演化拓扑优化方法的优点，表2分别列出了SIMP,BESO,LSM和BDE四种方法在不同单元网格和过滤半径下的结果。对比结果由表4给出：

表4.SIMP,BESO,LSM和BDE四种方法的结果对比

从表4看出，在网格划分不一样的时候，四种方法的得到的优化目标值非常接近，因为四种方法均为数值算法，得到的都是局部最优解。从本质上来讲，BDE方法就是在SIMP法的基础上去掉了惩罚函数并加上了密度分级过滤和边界光滑处理。相比SIMP法而言，收敛更快。而且第一阶段的优化结果的有限元模型与BESO法类似，没有灰度单元。同时由于加上了边界光滑处理，最后可以得到与LSM法完全类似的具有光滑边界的优化结果。总的来说BDE方法的优点是一可以得到具有光滑几何边界的优化结果，二是网格划分比较小的情况下可以得到高分辨率的结果。从优化图形上来看，BDE方法得到的优化结果也更接近Michell桁架的结果。

实施例2：

BDE优化方法可以很容易扩展到三维的情况，在另一个三维结构实例中采用BSE方法辅助设计一座桥梁。其初始设计域如图12所示，桥面为非设计域，桥面受到的节点荷载为10000KN，底面四点铰支，所用材料弹性模量为20Gpa，泊松比为0.2，体积约束0.2。将初始结构按对称取1/2模型进行计算，此时设计域离散为52800个8节点单元。

在家用IntelI5-CPU3.3GHz的电脑上计算，最后迭代步数为60步，计算耗时2375秒，最后优化结果应变能为36.91，目标和约束收敛曲线如图13所示。优化后的有限元模型如图14所示。为了边界更加光滑，图14计算的有限元细分为421400个单元，经过边界单元处理最终的桥梁结构优化模型如图15所示，图15的优化结果模型可直接采用3D打印。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (7)

1.一种具有光滑边界的连续体结构密度演化拓扑优化方法，其特征在于，包括以下步骤：

S1、通过画图建模模块建立连续体结构的几何初始设计域；

S2、基于无惩罚函数的材料插值模型建立几何初始设计域的有限元模型；把结构划分为若干个有限元单元网格，通过单元密度的有无确定优化的结果；

S3、建立连续体结构拓扑优化问题数学模型，以结构应变能最小为目标，所用材料体积为约束建立优化的数学模型；

S4、计算优化目标和约束条件对设计变量的灵敏度，并对灵敏度进行过滤；

S5、采用优化算法根据灵敏度更新设计变量，并对密度进行分级过滤，直至达到收敛和约束条件，得到优化结构的有限元模型；

S6、构建节点应变能水平集，对边界单元进行处理，得到基于节点谜底或节点应变能水平集显示的光滑边界的结构拓扑优化模型，并保存。

2.根据权利要求1所述的一种具有光滑边界的连续体结构密度演化拓扑优化方法，其特征在于，所述步骤S1具体为：画图建模模块根据连续体结构的结构形状画出连续体结构的初始几何实体图形，进而建立几何初始设计域。

3.根据权利要求1所述的一种具有光滑边界的连续体结构密度演化拓扑优化方法，其特征在于，所述步骤S2具体为：将几何初始设计域划分为N个有限元单元，指定每个对应有限元单元的相对密度为x_e，按照有限元的方法生成单元刚度矩阵，再组装成全局刚度矩阵，形成初始设计域的有限元模型，有限元模型中用到的材料的弹性模量的插值模型为：

E(x_e)＝E_min+x_e(E₀-E_min)，

其中，E₀为材料的实际弹性模量；E_min为无穷小量，物理含义为当单元密度为0时材料的弹性模量；x_e为第e个有限元单元密度。

4.根据权利要求1所述的一种具有光滑边界的连续体结构密度演化拓扑优化方法，其特征在于，所述步骤S3具体为：设置静力问题以结构应变能c最小作为目标函数，体积比f为约束条件，同时需要满足静力平衡方程，其优化问题数学模型如下：

其中，c为结构应变能，U为全局位移，F为负荷向量，K为全局刚度矩阵，u_e为单元位移向量，k₀为单元刚度矩阵，x_e为设计变量，即有限元单元密度；x为有限元单元密度向量，N为几何初始设计域采用有限元离散的单元总数，v(x)为优化结构可用材料体积，v₀为几何初始设计域的总体积，f为体积分数。

5.根据权利要求4所述的一种具有光滑边界的连续体结构密度演化拓扑优化方法，其特征在于，所述步骤S4具体如下：

采用伴随矩阵法求解目标函数对有限元单元密度x_e的灵敏度：

其中，c为结构应变能，E₀为材料的实际弹性模量；E_min为无穷小量，物理含义为当单元密度为0时材料的弹性模量；x_e为第e个有限元单元密度，也为设计变量；u_e为单元位移向量，k₀为单元刚度矩阵；

再计算约束条件对有限元单元密度x_e的灵敏度：

其中，v为结构总体积，x_e为第e个有限元单元的密度；

为了避免棋盘格现象，对敏感度进行过滤：

其中，n_e为在过滤半径r_min内与有限元单元e中心点的直线距离的所有有限元单元集合，γ为单元虚密度值的下限值；H_ei为权重系数，表示在过滤半径内各个有限元单元所占的权重；

权重系数定义如下：

H_ei＝max(0,r_min-Δ(e,i))，

其中，Δ(e,i)为有限元单元e中心点到有限元单元i中心点的直线距离。

6.根据权利要求1所述的一种具有光滑边界的连续体结构密度演化拓扑优化方法，其特征在于，所述步骤S5具体为：采用密度分级过滤方法进行更新设计变量，计算如下：

其中，a₀为中间参数，G_f为密度过滤分级因子，

为更新的设计变量，loop为当前循环步数，a₁为当前循环步数中设计变量小于a0(loop)对应的单元序号；

直至达到收敛和约束条件，得到只包含0/1单元的优化结构的有限元模型。

7.根据权利要求1所述的一种具有光滑边界的连续体结构密度演化拓扑优化方法，其特征在于，所述步骤S6具体为：

S601、单元细分；单元细分存在以下两种情况：需要得到高分辨率的优化结果，则需要过滤半径小于1.5；对于三维结构，其优化过程计算量大；两种情况都需要进行单元细分；

S602、构造节点应变能的水平集；单元的应变能为：

其中，x_e为第e个有限元单元的密度，u_i为第e个单元对应节点的位移，k₀为第e个单元的单元刚度矩阵；

节点应变能通过插值得到，对于二维平面问题，单元上任一点的应变能水平集函数表示为：

其中，Φ(o)为任一点o的水平集的值，

为i个单元的第j个节点上的应变能，N_j(o)为插值函数；具体表达式为：

N_j(o)＝(1+ξ_oξ_j)(1+η_oη_j)，

其中(ξ_o,η_o)为节点o在坐标系中对应的坐标；

基于应变能水平集函数，则可根据不同的水平集值把设计区域分为实体区域Ω、空心区域D\Ω、边界区域Γ，区分如下：

在迭代过程中，优化后的有限元模型分为实心单元、空心单元、边界单元；对单元密度的处理如下：

其中，

为i个单元的第j个节点上的应变能的值，S为水平集值；

统计单元应变能数值大于水平集值S的个数记为N_lb，单元总数记为N_mb，则边界单元的密度为N_lb/N_mb；

采用二分法确定水平集值S的值，水平集值的上界S_ub和下界S_lb的初始值分别对应节点应变能的最大值和最小值，则水平集值为：S＝(S_ub+S_lb)/2；

当单元密度更新后，为了满足体积约束，通过对比优化设计模型的体积V和材料约束体积V₀，水平集值的上界和下届也随着更新，其更新策略为：

水平集值的上界S_ub和下界S_lb在迭代的过程中不断更新，直到S_ub和S_lb的值S_ub-S_lb＜10^-9，最后结构优化设计模型的体积V就会达到约束值V₀；

S603、最终优化的结果通过Φ(o)-S＝0的节点应变能零水平集可直接显示，其物理含义是当边界节点上的应变能都等于S时，对应的零水平集的结构满足体积约束条件，此时结构应变能最小；优化结果也可通过把边界单元经过处理后的有限元模型的单元密度转化为节点密度，采用节点密度ψ(o)-0.5＝0的水平集显示；得到的结果具有显式的光滑边界，存储为stl格式的文件。

Images