CN110069800B - 具有光滑边界表达的三维结构拓扑优化设计方法及设备 - Google Patents

Abstract

本发明公开了一种具有光滑边界表达的三维结构拓扑优化设计方法及设备，属于结构拓扑优化领域。该方法基于双向渐进拓扑优化建立优化问题的数学模型，结合目标结构的有限元模型进行参数初始化及有限元分析，得到单元位移矩阵及单元目标函数，从而求取单元敏度数；并对单元的节点敏度数进行过滤与平均，然后进行节点进化，生成新的设计；对新的设计判断是否同时满足体积比及收敛条件；均满足则获得最优设计方案，否则，返回有限元分析步骤，开始循环迭代优化。本发明可在划分较小网格尺寸的条件下，使所设计的三维结构始终具有光滑的结构边界，无需曲线/面的拟合等后处理操作，有效规避现有方法严重依赖网格尺寸的弊端。

Description

具有光滑边界表达的三维结构拓扑优化设计方法及设备

技术领域

本发明属于结构拓扑优化技术领域，更具体地，涉及一种具有光滑边界表达的三维结构拓扑优化设计方法。

背景技术

结构是指由梁、板、壳、三维实体等组成的能够承受和传递外荷载的系统。

结构拓扑是指构成材料在设计空间中的分布方式以及构件之间的连接方式。结构拓扑优化旨在满足材料体积、承载能力、刚度指标等约束条件的前提下，寻求具有某种最优性能的结构拓扑形式。因为拓扑优化不依赖于结构初始构型及工程师经验，可获得完全意想不到的优质构型，己成为结构、材料创新设计的重要工具。增材制造技术(智能制造领域)的发展进一步解放了设计人员的工艺限制，使得拓扑优化技术(智能设计领域)的发展更具广阔前景，受到学者与工程人员的广泛关注。

拓扑优化方法根据设计变量的不同可以分为两类：一种是基于几何信息描述的方法，代表性的有水平集法(LSM)，可移动变形组件法(MMC)以及可移动变形孔洞法；另一种是基于单元状态描述的方法，代表性的有变密度法(SIMP、RAMP)，渐进结构优化法(ESO)以及双向渐进结构优化法(BESO)；两类方法各有利弊。基于几何信息描述的方法程序实现繁琐，敏度信息推导复杂，计算效率低；但所设计的结构边界清晰，光滑。而基于单元状态描述的方法设计的结构边界不清晰，不光滑，且具有棋盘格，网格依赖性等弊端；此外，数值计算精度也不高。但该类方法的程序便于实现，敏度信息便于推导，计算效率高，尤其适用于三维结构的优化设计。

作为基于单元状态描述的重要方法之一，BESO凭借其概念简单，程序实现简便、计算高效等优点已被广泛应用于先进结构与材料的设计中。虽然BESO方法获得了长足的发展与良好的应用，但依然具有以下不足：设计结果依赖于网格尺寸，结构边界不光滑，带有锯齿状，数值计算结果不精确等。如何能在渐进结构优化理论的框架下，使得结构设计边界在不依赖于网格尺寸的条件下，总能得到光滑的结构设计边界以及良好地数值计算精度，一直是研究者亟待解决的问题。尤其是工程中的三维结构优化与设计，对于该需求尤为迫切。

发明内容

针对现有技术存在的以上缺陷或改进需求，本发明提供了一种具有光滑边界表达的三维结构拓扑优化设计方法，其目的在于，基于BESO方法的渐进优化框架，在划分较小网格尺寸的条件下，使所设计的三维结构始终具有光滑的结构边界，无需曲线/面的拟合等后处理操作，便可与现如今流行的工程制图软件直接连接，使得基于拓扑优化的智能设计技术与基于增材制造的智能制造技术形成闭环。

为实现上述目的，按照本发明的一个方面，提供了一种具有光滑边界表达的三维结构拓扑优化设计方法，包括以下步骤：

步骤1：基于双向渐近拓扑优化法，构建具有指定的最优性能的结构拓扑优化模型；

步骤2：给定目标三维结构的初始设计域、施加载荷与边界条件，划分有限元网格，建立有限元模型；

步骤3：定义步骤1的优化设计模型的初始化设计参数，包括：进化率er、过滤半径r_min、惩罚因子p、体分比V^*以及单元内插值点数ngrid；

步骤4：基于步骤2的有限元模型以及步骤3定义的初始化参数进行有限元分析，得到单元位移矩阵以及单元应变能；

步骤5：将单元应变能进行组装，合成结构整体应变能；并基于步骤4得到的单元位移矩阵进行灵敏度分析，求得单元敏度数；

步骤6：将单元敏度数转化为节点敏度数，并进行节点敏度数的过滤与平均操作；

步骤7：对经步骤6过滤及平均后的各个节点敏度数按由大到小的规则进行排序，随后确定设计域内插值点材料的添加或去除阈值、每个单元八个节点的水平集函数值以及每个单元内各个插值点材料的水平集函数值；依据设计域内各个插值点材料的添加或删除，实现结构的演化与更新；

步骤8：判断结构当前体积V_k是否满足预设的目标体积V^*，若满足，则V_k＝V^*，转步骤9进行收敛条件的判断；若不满足，则令

转至步骤4继续循环迭代；其中，V_k与

分别表示结构当前迭代的体积与上一步迭代的体积，er为结构进化率；

步骤9：判断是否满足预设的收敛条件，若满足，则优化循环结束，输出当前优化结果作为最优设计结果；否则，转至步骤4继续循环迭代。

进一步地，步骤1中的指定性能为结构柔度、全局或局部应力、全局或局部频率中的任意一种，对于结构柔度最小化问题，其拓扑优化模型如下：

其中，C表示结构柔度，U与K分别表示结构的位移向量与总刚度矩阵；V是结构当前的体分比，N为设计域内离散的单元总数；V_i表示第i个单元的体积，V_i ^f表示优化过程中第i个单元所占的体分比，V^*是预设的结构目标体分比；s是二进制变量，s＝1表示设计域内任一插值点材料的实状态，s＝0表示设计域内任一插值点材料的虚状态，D表示结构设计域。

进一步地，步骤1中的V_i ^f有三种情况：当一个单元内的所有插值点状态为s＝1时，该单元为实单元，V_i ^f＝1；当一个单元内的所有插值点状态为s＝0时，该单元为虚单元，V_i ^f＝0；当一个单元内的部分插值点状态为s＝1，而另一部分的插值点状态为s＝0时，该单元为边界单元，V_i ^f∈(0,1)。

进一步地，步骤5中求取单元敏度数的方法如下：

其中，k_i为第i个单元的刚度矩阵，

表示第i个单元内因插值点材料的添加或删除而不断变化的单元状态，x_i＝1对应于实单元，

对应于虚单元，x_i∈(x_min,1)对应于边界单元，x_min为预设的远小于1的参数，以避免有限元分析中总刚度矩阵K的奇异；p为惩罚因子；k₀为24*24的矩阵常量；

其中，α_i为单元i的敏度数，V_i ^f为第i个单元的体积分数，u_i为第i个单元的位移矩阵。

进一步地，当材料的弹性模量为E₀，泊松比为ν时，k₀的计算公式如下：

其中，元素k₁、k₂、k₃、k₄、k₅、k₆、

与

的计算公式如下：

在k₁、k₂、k₃、k₄、k₅、k₆、

与

中的各矩阵元素分别为：

e₁＝-(6ν-4)/9，e₂＝1/12，e₃＝-1/9，

e₄＝-(4ν-1)/12，e₅＝(4ν-1)/12，e₆＝1/18，

e₇＝1/24，e₈＝-1/12，e₉＝(6ν-5)/36，

e₁₀＝-(4ν-1)/24，e₁₁＝-1/24，e₁₂＝(4ν-1)/24，

e₁₃＝(3ν-1)/18，e₁₄＝(3ν-2)/18。

进一步地，步骤6中，进行节点敏度数的过滤与平均操作的方法如下：

ω_ij＝max(0,r_min-r(i,j)) (8)

其中，

表示第j个节点过滤后的敏度数，ω_ij是表示距离j节点远近的线性权重因子，α_i表示第i个节点的敏度数；

表示第j个单元第k次迭代过滤后的敏度数，

表示第j个单元第k-1次迭代过滤后的敏度数，

表示第j个单元第k次迭代平均后的敏度数。

进一步地，步骤7中，按照公式(10)得出设计域内插值点材料的添加或去除阈值，按照公式(11)得出每个单元八节点的水平集函数值，按照公式(12)计算出每个单元内各个插值点材料的水平集函数值：

th表示插值点的添加与删除阈值，dcnd_max表示单元节点敏度数的最大值，dcnd_min表示单元节点敏度数的最小值；

Φ＝dcnd-th (11)

Φ表示单元节点的水平集函数值，dcnd表示单元节点的敏度数；

Φ(ξ)表示插值点ξ的水平集函数值，Φ(N_j)表示插值点ξ所在单元第j个节点的水平集函数值，

表示插值点ξ的形状函数，j＝1～8；

x_ξ、y_ξ、z_ξ分别表示插值点ξ在x、y、z方向的坐标；

当插值点的水平集函数值Φ＞0时，为实点；当插值点的水平集函数值Φ＝0时，为边界点；当插值点的水平集函数值Φ＜0时，为虚点。

进一步地，步骤(9)中的收敛条件定义为：

其中，var表示目标函数的变化率，C_k-i+1与C_k-M-i+1分别表示第k-i+1次、第k-M-i+1次迭代计算所得的结构柔度，k为当前迭代次数，τ是允许的目标函数收敛误差，M为预设的正整数。

为了实现上述目的，另一方面，本发明还提供了一种计算机可读存储介质，该计算机可读存储介质上存储有计算机程序，该计算机程序被处理器执行时实现如前所述的任意一种三维结构拓扑优化设计方法。

为了实现上述目的，另一方面，本发明还提供了一种具有光滑边界表达的三维结构拓扑优化设计设备，包括如前所述的计算机可读存储介质以及处理器，处理器用于调用和处理计算机可读存储介质中存储的计算机程序。

总体而言，本发明所构思的以上技术方案与现有技术相比，能够取得下列有益效果：

(1)BETO方法继承了传统BESO方法的启发式优化思想，概念简单，程序容易实现，具有良好地数值计算精度等突出优点，可广泛应用于二维与三维结构的优化设计；

(2)众所周知，每步迭代中的有限元分析操作对拓扑优化的计算效率有较大影响，而BETO方法通过在单元内均布点材料的方式，使得插值点的敏度数计算独立于有限元分析，因而可在较稀疏网格尺寸条件下，有效规避BESO方法设计结果严重依赖网格尺寸的弊端，并得到光滑的结构边界表达以及良好地数值计算精度，拓扑结果可与现今流行的工程制图软件(如Pro/E，Solidworks等CAD软件)直接进行连接，而无需曲线/面拟合等额外的后处理操作。

附图说明

图1是本发明较佳实施提供的一种具有光滑边界表达的三维拓扑优化设计方法流程图；

图2采用图2中的一种设计周期性多材料结构的拓扑优化方法进行优化的悬臂梁结构的设计域与边界条件；

图3是任一立方体单元内的节点编号示意图；

图4是图2中的悬臂梁结构优化后得到的最优结构拓扑图，其中，(a)为BESO方法计算得到的结果，(b)为BETO方法计算得到的结果。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

图1为本发明实施例提供的一种具有光滑边界表达的三维结构拓扑优化设计方法流程图。

(1)基于双向渐近拓扑优化法(简称：BETO法)构建三维结构的拓扑优化模型；

其中，C代表结构柔度，U与K分别表示结构的位移向量与总刚度矩阵。V是结构当前的体分比，N意味着设计域内离散的单元总数。一般情况下V_i＝1，表示第i个单元的体积，V_i ^f表示优化过程中结构在第i个单元所占的体分比，V^*是设计人员预设的结构目标体分比；s是二进制变量，表示设计域内任一插值点的实状态(s＝1)或虚状态(s＝0)，D是设计域。

由于刚度与柔度是一对相反的物理量，该模型能够最小化结构柔度或最大化结构刚度。

此外，由于本发明的方法具有通用性，还适用于最小化结构全局应力、最大化结构基频等其他指定性能的优化设计，只需要替换函数(1)中的目标函数(min)及约束条件(subject to)为相应性能的计算公式即可。

(2)如图2所示，给定结构初始设计域，施加载荷与边界条件，划分有限元网格；

(3)定义BETO方法初始化设计参数，主要包括：进化率er、过滤半径rmin、惩罚因子p、体分比V^*以及单元内插值点数ngrid；

(4)根据边界条件与载荷施加的工况，对步骤(2)的有限元模型进行有限元分析，得到单元位移矩阵；

F_i＝k_iu_i (2)

其中，F_i为第i个单元的载荷，k_i为第i个单元的刚度矩阵，u_i为第i个单元的位移矩阵；

(5)基于步骤(4)得到的单元位移矩阵进行灵敏度分析，求解单元敏度数；

其中，k_i为第i个单元的刚度矩阵，

代表第i个单元内因插值点的添加或删除而变化的单元状态，进一步地，x_i＝1对应实单元，

对应虚单元，x_i∈(x_min,1)对应边界单元。x_min为预设的远小于1的参数，以避免有限元分析中总刚度矩阵K的奇异，一般情况下，可以取x_min＝0.001。p为惩罚因子3，α_i为单元i的敏度数，V_i ^f为第i个单元的体积分数。对某种指定的各向同性材料来说，k₀为24*24的矩阵常量。

当材料的弹性模量为E₀，泊松比为ν时，k₀的计算公式如下：

其中，元素k₁、k₂、k₃、k₄、k₅、k₆、

与

的计算公式如下：

在k₁、k₂、k₃、k₄、k₅、k₆、

与

中的各矩阵元素分别为：

e₁＝-(6ν-4)/9，e₂＝1/12，e₃＝-1/9，

e₄＝-(4ν-1)/12，e₅＝(4ν-1)/12，e₆＝1/18，

e₇＝1/24，e₈＝-1/12，e₉＝(6ν-5)/36，

e₁₀＝-(4ν-1)/24，e₁₁＝-1/24，e₁₂＝(4ν-1)/24，

e₁₃＝(3ν-1)/18，e₁₄＝(3ν-2)/18。

(6)将单元敏度数转化为节点敏度数，并分别根据以下公式进行节点敏度数的过滤与平均操作；

ω_ij＝max(0,r_min-r(i,j)) (8)

其中，

代表第j个节点过滤后的敏度数，ω_ij是一个表示距离j节点远近的线性权重因子，α_i表示第i个节点的敏度数。

表示第j个单元第k次迭代过滤后的敏度数，

表示第j个单元第k-1次迭代过滤后的敏度数，

表示第j个单元第k次迭代平均后的敏度数。

(7)对节点敏度数按由大到小的规则进行排序，按照以下公式得出设计域内插值点的添加或去除阈值。

th代表插值点的添加与删除阈值，dcnd_max表示单元节点敏度数的最大值，dcnd_min表示单元节点敏度数的最小值。利用下述公式得出每个单元8节点的水平集函数值，

Φ＝dcnd-th (11)

Φ代表单元节点的水平集函数值，dcnd表示单元节点的敏度数。利用以下公式计算出每个单元内各个插值点的水平集函数值。

Φ(ξ)代表插值点ξ的水平集函数值，

表示插值点ξ的形状函数，Φ(N₁),Φ(N₂)...Φ(N₈)表示插值点ξ所在单元8个节点的水平集函数值，j＝1～8，节点编号如图3所示。

x_ξ、y_ξ、z_ξ分别表示插值点ξ在x、y、z方向的坐标；

当节点的水平集函数值Φ(ξ)＞0时，为实点；当节点的水平集函数值Φ(ξ)＝0时，为边界点；当节点的水平集函数值Φ(ξ)＜0时，为虚点。依据设计域内各个节点的添加或删除，从而实现结构的演化与更新。值得说明的是假定每个单元内沿着x、y、z方向分别均布10个点，则一个立方体单元内就均布了1000个节点。

(8)判断结构当前体积V_k是否满足预设的目标体积V^*，若满足，则V_k＝V^*，进入下一步收敛条件的判断；若不满足，则令

转至步骤4继续循环迭代。其中，V_k与

分别表示结构当前迭代的体积与上一步迭代的体积，er为结构进化率。

(9)判断收敛条件是否满足，若满足，则优化循环结束；否则，转至步骤4进行循环。收敛条件定义为：

其中，var表示目标函数的变化率，C_k-i+1与C_k-M-i+1分别代表第k-i+1次、

第k-M-i+1次迭代计算所得的结构柔度，k为当前迭代次数，τ是允许的目标函数收敛误差，本实施例取τ＝0.0001，M＝5。

如图4所示的BESO方法计算结果与BETO方法计算结果，对比可知，传统BESO方法的计算结果具有较为明显的锯齿状边界，而采用本发明的BETO方法，其通过单元内插值点材料的方式，使得三维结构计算结果具有光滑的边界。此外，本发明提供的BETO方法还具备良好地数值计算精度。

通过对比可知，本发明所提出的BETO优化法，沿袭BESO方法的渐进结构优化概念简单，程序便于实现，计算高效等优点的基础上，基于较小的网格尺寸以及较低的计算成本，设计出具有光滑边界表达的三维结构，并提高其数值计算精度。设计出的结构无需后处理操作，可直接与现今流行的CAD软件无缝连接，真正实现优质结构智能设计与智能制造的闭环。

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

Claims (10)

1.一种具有光滑边界表达的三维结构拓扑优化设计方法，其特征在于，包括以下步骤：

步骤1：基于双向渐近拓扑优化法，构建具有指定的最优性能的结构拓扑优化模型；

步骤2：给定目标三维结构的初始设计域、施加载荷与边界条件，划分有限元网格，建立有限元模型；

步骤3：定义步骤1的优化设计模型的初始化设计参数，包括：进化率er、过滤半径r_min、惩罚因子p、体分比V^*以及单元内插值点数ngrid；

步骤4：基于步骤2的有限元模型以及步骤3定义的初始化参数进行有限元分析，得到单元位移矩阵以及单元应变能；

步骤5：将单元应变能进行组装，合成结构整体应变能；并基于步骤4得到的单元位移矩阵进行灵敏度分析，求得单元敏度数；

步骤6：将单元敏度数转化为节点敏度数，并进行节点敏度数的过滤与平均操作；

步骤7：对经步骤6过滤及平均后的各个节点敏度数按由大到小的规则进行排序，随后确定设计域内插值点材料的添加或去除阈值、每个单元八个节点的水平集函数值以及每个单元内各个插值点材料的水平集函数值；依据设计域内各个插值点材料的添加或删除，实现结构的演化与更新；

步骤8：判断结构当前体积V_k是否满足预设的目标体积V^*，若满足，则V_k＝V^*，转步骤9进行收敛条件的判断；若不满足，则令

转至步骤4继续循环迭代；其中，V_k与

分别表示结构当前迭代的体积与上一步迭代的体积，er为结构进化率；

步骤9：判断是否满足预设的收敛条件，若满足，则优化循环结束，输出当前优化结果作为最优设计结果；否则，转至步骤4继续循环迭代。

2.如权利要求1所述的一种具有光滑边界表达的三维结构拓扑优化设计方法，其特征在于，步骤1中的指定性能为结构柔度、全局或局部应力、全局或局部频率中的任意一种，对于结构柔度最小化问题，其拓扑优化模型如下：

其中，C表示结构柔度，U与K分别表示结构的位移向量与总刚度矩阵；V是结构当前的体分比，N为设计域内离散的单元总数；V_i表示第i个单元的体积，V_i ^f表示优化过程中第i个单元所占的体分比，V^*是预设的结构目标体分比；s是二进制变量，s＝1表示设计域内任一插值点材料的实状态，s＝0表示设计域内任一插值点材料的虚状态，D表示结构设计域。

3.如权利要求2所述的一种具有光滑边界表达的三维结构拓扑优化设计方法，其特征在于，步骤1中的V_i ^f有三种情况：当一个单元内的所有插值点状态为s＝1时，该单元为实单元，V_i ^f＝1；当一个单元内的所有插值点状态为s＝0时，该单元为虚单元，V_i ^f＝0；当一个单元内的部分插值点状态为s＝1，而另一部分的插值点状态为s＝0时，该单元为边界单元，V_i ^f∈(0,1)。

4.如权利要求1所述的一种具有光滑边界表达的三维结构拓扑优化设计方法，其特征在于，步骤5中求取单元敏度数的方法如下：

其中，k_i为第i个单元的刚度矩阵，

表示第i个单元内因插值点材料的添加或删除而不断变化的单元状态，x_i＝1对应于实单元，

对应于虚单元，x_i∈(x_min,1)对应于边界单元，x_min为预设的远小于1的参数，以避免有限元分析中总刚度矩阵K的奇异；p为惩罚因子；k₀为24*24的矩阵常量；

其中，α_i为单元i的敏度数，V_i ^f为第i个单元的体积分数，u_i为第i个单元的位移矩阵。

5.如权利要求4所述的一种具有光滑边界表达的三维结构拓扑优化设计方法，其特征在于，当材料的弹性模量为E₀，泊松比为ν时，k₀的计算公式如下：

其中，元素k₁、k₂、k₃、k₄、k₅、k₆的计算公式如下：

在k₁、k₂、k₃、k₄、k₅、k₆中的各矩阵元素分别为：

e₁＝-(6ν-4)/9， e₂＝1/12， e₃＝-1/9，

e₄＝-(4ν-1)/12， e₅＝(4ν-1)/12， e₆＝1/18，

e₇＝1/24， e₈＝-1/12， e₉＝(6ν-5)/36，

e₁₀＝-(4ν-1)/24， e₁₁＝-1/24， e₁₂＝(4ν-1)/24，

e₁₃＝(3ν-1)/18， e₁₄＝(3ν-2)/18。

6.如权利要求1～5任意一项所述的一种具有光滑边界表达的三维结构拓扑优化设计方法，其特征在于，步骤6中，进行节点敏度数的过滤与平均操作的方法如下：

ω_ij＝max(0,r_min-r(i,j)) (8)

其中，

表示第j个节点过滤后的敏度数，ω_ij是表示距离j节点远近的线性权重因子，α_i表示第i个节点的敏度数；

表示第j个单元第k次迭代过滤后的敏度数，

表示第j个单元第k-1次迭代过滤后的敏度数，

表示第j个单元第k次迭代平均后的敏度数。

7.如权利要求1～5任意一项所述的一种具有光滑边界表达的三维结构拓扑优化设计方法，其特征在于，步骤7中，按照公式(10)得出设计域内插值点材料的添加或去除阈值，按照公式(11)得出每个单元八节点的水平集函数值，按照公式(12)计算出每个单元内各个插值点材料的水平集函数值：

th表示插值点的添加与删除阈值，dcnd_max表示单元节点敏度数的最大值，dcnd_min表示单元节点敏度数的最小值；

Φ＝dcnd-th (11)

Φ表示单元节点的水平集函数值，dcnd表示单元节点的敏度数；

Φ(ξ)表示插值点ξ的水平集函数值，Φ(N_j)表示插值点ξ所在单元第j个节点的水平集函数值，

表示插值点ξ的形状函数，j＝1～8；

x_ξ、y_ξ、z_ξ分别表示插值点ξ在x、y、z方向的坐标；

当插值点的水平集函数值Φ＞0时，为实点；当插值点的水平集函数值Φ＝0时，为边界点；当插值点的水平集函数值Φ＜0时，为虚点。

8.如权利要求1～5任意一项所述的一种具有光滑边界表达的三维结构拓扑优化设计方法，其特征在于，步骤(9)中的收敛条件定义为：

其中，var表示目标函数的变化率，C_k-i+1与C_k-M-i+1分别表示第k-i+1次、第k-M-i+1次迭代计算所得的结构柔度，k为当前迭代次数，τ是允许的目标函数收敛误差，M为预设的正整数。

9.一种计算机可读存储介质，其特征在于，该计算机可读存储介质上存储有计算机程序，该计算机程序被处理器执行时实现如权利要求1～8任一项所述的三维结构拓扑优化设计方法。

10.一种具有光滑边界表达的三维结构拓扑优化设计设备，其特征在于，包括如权利要求9所述的计算机可读存储介质以及处理器，处理器用于调用和处理计算机可读存储介质中存储的计算机程序。

Images