CN109871574A - 一种基于代理模型的多尺度拓扑优化方法 - Google Patents

Abstract

本发明属于结构优化设计技术领域，并公开了一种基于代理模型的多尺度拓扑优化方法。包括下列步骤：(1)采用基于参数化水平集的拓扑优化方法优化得到多种原型微结构的最优拓扑构型；(2)对优化后的原型微结构的水平集函数进行插值，以此获得非均匀微观结构的等效属性，随后将这些非均匀的微观结构当作样本点构建预测模型以预测宏观结构域内所有宏观单元的等效属性；(3)采用变厚度法优化宏观结构域内的材料分布，使宏观结构的性能最优。通过本发明，实现宏观结构几何、逐点不同的宏观单元构型及其在宏观结构域内位置的联合优化，以较少的计算成本，最大限度地发挥材料潜力，节约了耗材，降低了成本。

Description

一种基于代理模型的多尺度拓扑优化方法

技术领域

本发明属于结构优化设计技术领域，更具体地，涉及一种基于代理模型的多尺度拓扑优化方法。

背景技术

多尺度结构是一种典型的多孔结构，具有高比刚度、高比强度、吸能、热绝缘等优异的力学性能，广泛应用于航空航天、汽车工业、通信电子等相关结构设计领域。多尺度拓扑优化设计是新型的数字化结构设计方式，能充分发掘宏微观两个尺度上的设计潜能，以最少的材料用量或者最低的成本实现结构的最佳性能。

目前，本领域相关技术人员已做出了一些研究，如文献1：“Wang Y,Wang MY,ChenF(2016)Structure-material integrated design by level sets Structural andMultidisciplinary Optimization 54:1145-1156.”公开了一种基于水平集方法的多尺度结构设计方法。该方法假设宏观结构由一种微结构组成。又如文献2：“Li H,Luo Z,Gao L,Qin Q(2018a)Topology optimization for concurrent design of structures withmulti-patch microstructures by level sets Computer Methods in AppliedMechanics and Engineering 331:536-561.”公开了一种基于水平集方法的多尺度结构设计方法。该方法假设宏观结构由多种不同的微结构组成，相较于前一种假设，该方法进一步扩大了多尺度设计空间。然而，为充分发掘多尺度设计空间，应根据宏观结构域内逐点不同的材料属性需求设计相应的微观结构，即宏观结构域内的微结构是逐点不同的，从这一点上说，上述两种方法均限制了多尺度的设计空间。如文献3：“Xia L,Breitkopf P(2014)Concurrent topology optimization design of material and structure withinFE2nonlinear multiscale analysis framework Computer Methods in AppliedMechanics and Engineering278:524-542”公开了一种微结构逐点设计的多尺度设计方法。但该方法需要对宏观设计域内每一个微结构进行单独优化，充分发掘了宏微观尺度的设计空间，计算成本高昂。

因此针对上述问题，以较低的计算成本，设计宏观结构域内微结构属性逐点变化的多尺度结构，以充分发掘多尺度的设计空间，最大限度发挥材料潜力，是当前亟待解决的研究热点问题。

发明内容

针对现有技术的以上缺陷或改进需求，本发明提供了一种基于代理模型的多尺度拓扑优化方法，通过建立原型微结构并优化其构型，然后利用原型微结构差值后获取样本点并利用样本点建立预测模型，以此获得宏观单元的等效属性，最后建立优化模型获得宏观结构的构型使得宏观结构的性能最优，实现宏观结构几何、逐点不同的宏观单元构型及其在宏观结构域内位置的联合优化，充分发掘宏微观多尺度的设计空间，以较少的计算成本，最大限度地发挥材料潜力，提升结构性能节约了耗材，降低了成本。

为实现上述目的，按照本发明，提供了一种基于代理模型的多尺度拓扑优化方法，其特征在于，该方法包括下列步骤：

(a)对待处理对象的宏观结构进行网格划分，使得该宏观结构被划分为多个宏观单元，对每个所述宏观单元赋予初始等效密度值，设定多种原型微结构的初始构型和初始等效密度，根据每个所述宏观单元的初始等效密度值进行聚类，以此确定每个所述宏观单元对应的所述原型微观结构，采用基于参数化水平集的拓扑优化方法优化每个所述原型微结构的初始构型，以此获得每个所述原型微结构的最优拓扑结构构型；

(b)对于每个原型微结构，根据该原型微结构的初始等效密度设定其中非均匀微结构的数量，按照该数量对每个所述原型微结构的最优拓扑结构构型进行插值，以此获得每个非均匀微结构的构型，根据该构型采用均匀化方法计算每个所述非均匀微结构的等效属性，将每个所述非均匀微结构作为一个样本点，由所有的所述非均匀微结构组成的样本集构建关于等效密度和等效属性的预测模型，根据每个所述宏观单元的初始等效密度，利用所述预测模型预测每个所述宏观单元的等效属性；

(c)根据(b)预测获得的每个所述宏观单元的等效属性，建立关于所述宏观结构的等效密度和特征频率的优化模型，以此获得在所述宏观结构的特征频率最大的情况下对应的每个所述宏观单元的等效密度和所述宏观结构的最大特征频率；

(d)判断所述宏观结构的最大特征频率是否收敛，

若收敛，按照每个宏观单元的等效密度确定该宏观单元对应的原型微结构，并以该原型微结构的最优拓扑结构构型进行插值，由此获得每个所述宏观单元的构型，进而获得所述宏观结构的构型，即该宏观结构的多尺度拓扑优化结果；

若不收敛，将当前每个所述宏观单元的等效密度作为初始等效密度，当前每个所述原型微结构的最优拓扑结构构型作为原型微结构的初始构型，返回步骤(a)，直至所述宏观结构的最大特征频率收敛。

进一步优选地，在步骤(a)中，所述采用基于参数化水平集的拓扑优化方法优化每个所述原型微结构的初始构型，以此获得每个所述原型微结构的最优拓扑结构构型，优选按照下列步骤进行：

(a1)按照下列表达式计算所述原型微结构的等效属性

其中，是第m个原型微结构的等效属性，H是均匀化，D_pqrs是微观设计域内任意一点的弹性系数矩阵，，ij,kl与pq,rs均是施加单元测试测试应变的方向，是指在pq方向下的单元测试应变场，是由pq方向下的单元测试应变场引起的未知应变场，表示在ij方向下的微观位移场，是指在rs方向下的单元测试应变场，是由rs方向下的单元测试应变场引起的未知应变场，表示在kl方向下的微观位移场，H(·)是Heaviside函数，是第m个原型微结构的水平集函数，是第m个原型微结构的设计域，MI表示微观尺度下的物理量。

(a2)按照下列方程求解所述宏观结构的第k阶特征频率ω_k对应的特征向量u_k，

其中，u_k是第k阶特征频率ω_k的特征向量，K是宏观结构的整体刚度矩阵，M是宏观结构的整体质量矩阵；

(a3)以宏观结构基频最大为优化目标，各个所述原型微结构的等效密度为给定的材料用量，采用基于参数化水平集拓扑优化方法构建所述原型微结构的优化模型，该优化模型如下：

其中，是微观设计变量，MI表示α_m,n是微观尺度下的物理量，R是原型微结构的种类，N是一个原型微结构的节点数量，是第m个原型微结构在其设计域内的水平集函数，B是应变-位移矩阵，α^MI是微观设计变量的向量表示形式，是宏观结构域内被第m个原型微结构占据的区域，G_m是第m个原型微结构的体积约束，是第m个原型微结构的等效密度，是宏观固体单元的体积，是正则化的设计变量的上边界，是正则化的设计变量,，是正则化的设计变量,的下边界，N_s是形状函数矩阵，K是宏观结构的整体刚度矩阵，M是宏观结构的整体质量矩阵；

(a4)按照下列表达式计算目标函数及约束条件对微观尺度设计变量的敏度，具体如下：

其中，δ(·)是Dirac函数，是紧支撑径向基插值函数。

(a5)采用基于梯度的优化算法更新所述微观尺度设计变量并输出优化后的各原型微结构构型。

进一步优选地，在步骤(a)中，所述聚类优选按照下列关系式进行：

其中，是密度聚类后的第Ne个宏观单元密度值，表示密度聚类前的第Ne个宏观单元密度值，和分别表示第m-1个和m个原型微结构的等效密度。

进一步优选地，在步骤(b)中，所述预测模型优选采用Kriging模型。

进一步优选地，在步骤(c)中，所述优化模型优选按照下列表达式：

其中，是宏观单元相对密度，NE是宏观单元总数，u_k是第k阶特征频率ω_k的特征向量，K是宏观结构的整体刚度矩阵，M是宏观结构的整体质量矩阵，G是给定的材料用量约束条件，是宏观固体单元的体积，Ω^MA是整个宏观结构设计域，是限定的材料使用量，ρ_min是设计变量的最小值，f_pm是通过插值原型微结构得到的非均匀微结构的最小密度，ρ_max是设计变量的最大值。

进一步优选地，在步骤(d)中，所述判断所述宏观结构的最大特征频率是否收敛是判断前后两次计算获得的最大特征频率的差值是否小于预设阈值，若小于预设阈值则收敛。

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

1、本发明提供的方法采用了形状插值技术，与现有的技术相比，保证宏观结构域内所有宏观单元较好的连接性，采用Kriging模型预测宏观结构域内的所有宏观单元等效属性，计算成本低廉；

2、本发明提供的方法优化后的宏观结构由逐点变化的宏观单元组成，极大增加了优化设计的可行空间，同时由于了采用多种原型微结构，更进一步的扩大了该优势；

3、本发明提供的方法，实现了宏观结构几何、逐点不同的宏观单元构型及其在宏观结构域内位置的联合优化，本发明相较于传统的多尺度设计，不仅大大减少了计算成本，也极大的扩展了多尺度设计空间，大幅提升了结构性能。

附图说明

图1是本发明实施例中提供的一种基于代理模型的多尺度拓扑优化方法的流程图；

图2是采用图1中于代理模型的宏观结构基频最大化的多尺度拓扑优化方法来设计的宏观结构设计域、载荷及边界条件示意图；

图3是图2中宏观结构材料布局优化后的示意图；

图4是图3中宏观结构的第一阶模态形状示意图；

图5是图2中宏观结构的最优多尺度设计结构示意图；

图6是四种原型微结构的最优构型示意图图；

图7是图6中四种原型微结构插值得到的非均匀微结构等效属性相对于等效密度的变化示意图；

图8是采用本发明所提供的方法进行多尺度优化时宏观目标函数与约束条件的迭代过程示意图。

具体实施方式

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

一种基于代理模型的多尺度拓扑优化方法，其流程图如图1所示，主要包括如下步骤：

步骤一，有限元网格划分，获得多个宏观单元，设定多种原型微结构的初始构型，采用均匀化理论计算多种原型微结构的等效属性；基于参数化水平集的拓扑优化方法优化得到多种原型微结构的最优拓扑构型，具体包括以下子步骤：

(1)初始化设计参数，包括材料属性，原型微结构设计域的长宽、水平集函数及其有限元模型，以及优化设计参数；

(2)采用均匀化方法计算多种原型微结构的等效属性D_ijH_kl，所述均匀化方法计算公式如下：

其中,m表示是第m个原型微结构的等效属性，H表示均匀化，D_pqrs表示微观设计域内任意一点的弹性系数矩阵，i,j,k,l与p,q,r,s均指施加单元测试测试应变的方向，ij＝11,22,12，分别表示水平方向、竖直方向、剪切方向，kl,pq,rs与ij类似，是指在pq方向下的单元测试应变场，是由pq方向下的单元测试应变场引起的未知应变场，表示在ij方向下的微观位移场，是指在rs方向下的单元测试应变场，是由rs方向下的单元测试应变场引起的未知应变场，表示在kl方向下的微观位移场，H(·)表示Heaviside函数，用于表述结构形式的特征函数

(3)基于宏观单元密度分布进行密度聚类，即用具有等效密度的原型微结构代替密度区间的所有宏观单元，所述聚类机制如下：

其中表示密度聚类后的第Ne个宏观单元密度值，表示密度聚类前的第Ne个宏观单元密度值，和分别表示第m-1个和m个原型微结构的等效密度。

(4)根据聚类后的密度分布，将多种原型微结构的等效属性代入宏观平衡方程求解宏观结构的第k阶特征频率ω_k对应的特征向量u_k，宏观平衡方程如下：

其中u_k是第k阶特征频率ω_k的特征向量，K是宏观结构的整体刚度矩阵，M是宏观结构的整体质量矩阵。

(5)以宏观结构基频最大为优化目标，各原型微结构的等效密度为给定的材料用量，即约束条件，基于参数化水平集拓扑优化方法构建原型微结构的优化模型，优化模型如下：

其中，是微观设计变量，MI表示α_m,n是微观尺度下的物理量，R是原型微结构的种类，N是一个原型微结构的节点数量，是第m个原型微结构在其设计域内的水平集函数，可通过CSRBFs插值表示为:

其中，是第m个原型微结构的CSRBF插值函数，是CSRBF插值函数的扩展系数，即第m个原型微结构的设计变，u_k是第k阶特征频率ω_k的特征向量，通过将R个原型微结构的等效属性代入宏观平衡方程求解，T表示矩阵的转置，B是应变-位移矩阵，N是形状函数矩阵，表示宏观结构域内被第m个原型微结构占据的区域，MA表示Ω_m是宏观尺度下的物理量，表示宏观固体单元的体积，G_m表示第m个原型微结构的体积约束，满足给定的体积是第m个原型微结构的等效密度，和是的上下边界，正则化的设计变量，便于后续的优化算法求解。

(6)计算目标函数及约束条件对微观尺度设计变量的敏度信息，所述敏度计算公式如下：

其中

δ(·)为Dirac函数，用于计算Heaviside函数的偏导数。

(7)采用基于梯度的优化算法更新微观尺度设计变量并输出优化后的各原型微结构构型。

步骤二，采用形状插值技术对优化后的原型微结构的水平集函数进行插值，以得到一系列非均匀的微观结构，以这些非均匀的微观结构为样本点构建Kriging模型，随后基于构建的Kriging模型预测宏观结构域内所有宏观单元的等效属性，具体包括以下子步骤：

(1)对于每一种原型微结构，设定非均匀微结构数量，根据该数量插值到原型微结构中通过插值各原型微结构，以此获得每一个非均匀微结构的构型；

(2)采用均匀化方法计算插值得到的非均匀微结构的等效属性，并将每个非均匀微结构作为一个样本点，以此作为构建Kriging模型的样本点；

(3)构建每个原型微结构的Kriging模型，并用构建的Kriging模型预测宏观结构域内所有宏观单元的等效属性，其中，每一个原型微结构的Kriging模型具有指定的预测密度区间，所有原型微结构的Kriging模型的密度预测区间之和覆盖区间[f_pm,1]，f_pm＝0.05是最小的可预测密度值。。

步骤三，基于Kriging模型预测的微观结构属性，在给定的约束条件下，采用变厚度法优化宏观结构域内的材料分布，以使宏观结构基频最大，具体包括以下子步骤：

(1)采用预测的所有宏观单元的等效属性，实施宏观有限元分析，以得到宏观结构的第k阶特征频率ω_k对应的特征向量u_k；

(2)以宏观结构基频最大为优化目标，给定的材料用量为约束，基于变厚度法构建宏观材料分布优化模型，所述优化模型如下：

其中为宏观尺度设计变量，即宏观单元相对密度，NE表示宏观单元总数，u_k是第k阶特征频率ω_k的特征向量，K是宏观结构的整体刚度矩阵，M是宏观结构的整体质量矩阵，G表示给定的材料用量约束条件，表示宏观固体单元的体积，Ω^MA表示整个宏观结构设计域，为限定的材料使用量，ρ_min为设计变量的最小值，一般取0.001，目的是避免刚度矩阵，f_pm表示可通过插值原型微结构得到的非均匀微结构的最小密度，一般取0.05，ρ_max是设计变量的最大值，一般取1。

(3)计算目标函数及约束条件对宏观尺度设计变量的敏度信息，所述计算敏度信息公式如下：

其中

(4)采用基于梯度的优化算法更新宏观尺度设计变量并输出优化后的宏观材料布局形式。

(5)判断目标函数是否满足收敛条件，若满足收敛条件，则输出当前宏观结构的最优多尺度设计结果，否则转至步骤一。

请参阅图2至图8，以下以两端夹持的宏观结构的设计来进一步说明本发明。如图2所示，两端夹持的宏观结构的设计域尺寸为长1.4m，高0.2m，设计域左右两端所有自由度被固定约束，设计域中间位置施加一集中质量点，所用材料的材料属性为弹性模量E₀＝201GPa，泊松比μ＝0.3，密度ρ₀＝7.8×10³kg/m³，原型微结构的设计域采用50×50＝2500个四边形单元进行网格划分，宏观设计域采用210×30＝6300个四边形单元进行网格划分，优化目标为宏观结构的特征频率最大，限定的材料用量为40％。

如图3所示为通过变厚度法优化后的宏观材料分布图，图4为优化后的宏观结构的第一阶模态形状图，图5是宏观结构的最优多尺度设计，具有最大化的特征频率136.5924HZ，由图5可知，宏观结构的最优多尺度设计包含逐点不同的宏观单元，高密度的宏观单元主要分布在宏观结构的外部区域，而低密度的宏观单元主要分布在宏观结构的内部区域，该种宏观单元分布模式能有效抵抗如图4中所示的第一阶振动模态的变形。图6是四种原型微结构的最优构型，分别具有等效密度0.2、0.4、0.6和0.8。图7是图6中四种原型微结构插值得到的非均匀微结构等效属性相对于等效密度的变化图。由图6和图7可知，通过插值原型微结构得到的所有非均匀微结构均具有显著的正交各向异性属性，能灵活地根据宏观结构域内各点处的材料属性需求提供不同的等效属性，最大程度地利用材料潜能，使宏观结构性能最优。图8是采用本发明所提供的方法进行多尺度优化时宏观目标函数与约束条件的迭代过程示意图，在迭代115步后，目标函数与约束条件均稳定收敛，显示了该设计方法的高效性。本发明提供的一种基于代理模型的多尺度拓扑优化方法，实现了宏观结构几何、逐点不同的宏观单元构型及其在宏观结构域内位置的联合优化，本发明相较于传统的多尺度设计，不仅大大减少了计算成本，也极大的扩展了多尺度设计空间，大幅提升了结构性能，同时也可保证宏观结构域内逐点不同宏观单元的连接性。

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

Claims (6)

1.一种基于代理模型的多尺度拓扑优化方法，其特征在于，该方法包括下列步骤：

(a)对待处理对象的宏观结构进行网格划分，使得该宏观结构被划分为多个宏观单元，对每个所述宏观单元赋予初始等效密度值，设定多种原型微结构的初始构型和初始等效密度，根据每个所述宏观单元的初始等效密度值进行聚类，以此确定每个所述宏观单元对应的所述原型微观结构，采用基于参数化水平集的拓扑优化方法优化每个所述原型微结构的初始构型，以此获得每个所述原型微结构的最优拓扑结构构型；

(b)对于每个原型微结构，根据该原型微结构的初始等效密度设定其中非均匀微结构的数量，按照该数量对每个所述原型微结构的最优拓扑结构构型进行插值，以此获得每个非均匀微结构的构型，根据该构型采用均匀化方法计算每个所述非均匀微结构的等效属性，将每个所述非均匀微结构作为一个样本点，由所有的所述非均匀微结构组成的样本集构建关于等效密度和等效属性的预测模型，根据每个所述宏观单元的初始等效密度，利用所述预测模型预测每个所述宏观单元的等效属性；

(c)根据(b)预测获得的每个所述宏观单元的等效属性，建立关于所述宏观结构的等效密度和特征频率的优化模型，以此获得在所述宏观结构的特征频率最大的情况下对应的每个所述宏观单元的等效密度和所述宏观结构的最大特征频率；

(d)判断所述宏观结构的最大特征频率是否收敛，

若收敛，按照每个宏观单元的等效密度确定该宏观单元对应的原型微结构，并以该原型微结构的最优拓扑结构构型进行插值，由此获得每个所述宏观单元的构型，进而获得所述宏观结构的构型，即该宏观结构的多尺度拓扑优化结果；

若不收敛，将当前每个所述宏观单元的等效密度作为初始等效密度，当前每个所述原型微结构的最优拓扑结构构型作为原型微结构的初始构型，返回步骤(a)，直至所述宏观结构的最大特征频率收敛。

2.如权利要求1所述的一种基于代理模型的多尺度拓扑优化方法，其特征在于，在步骤(a)中，所述采用基于参数化水平集的拓扑优化方法优化每个所述原型微结构的初始构型，以此获得每个所述原型微结构的最优拓扑结构构型，优选按照下列步骤进行：

(a1)按照下列表达式计算所述原型微结构的等效属性

其中，是第m个原型微结构的等效属性，H是均匀化，D_pqrs是微观设计域内任意一点的弹性系数矩阵，，ij,kl与pq,rs均是施加单元测试测试应变的方向，是指在pq方向下的单元测试应变场，是由pq方向下的单元测试应变场引起的未知应变场，表示在ij方向下的微观位移场，是指在rs方向下的单元测试应变场，是由rs方向下的单元测试应变场引起的未知应变场，表示在kl方向下的微观位移场，H(·)是Heaviside函数，是第m个原型微结构的水平集函数，是第m个原型微结构的设计域，MI表示微观尺度下的物理量。

(a2)按照下列方程求解所述宏观结构的第k阶特征频率ω_k对应的特征向量u_k，

其中，u_k是第k阶特征频率ω_k的特征向量，K是宏观结构的整体刚度矩阵，M是宏观结构的整体质量矩阵；

(a3)以宏观结构基频最大为优化目标，各个所述原型微结构的等效密度为给定的材料用量，采用基于参数化水平集拓扑优化方法构建所述原型微结构的优化模型，该优化模型如下：

其中，是微观设计变量，MI表示α_m,n是微观尺度下的物理量，R是原型微结构的种类，N是一个原型微结构的节点数量，是第m个原型微结构在其设计域内的水平集函数，B是应变-位移矩阵，α^MI是微观设计变量的向量表示形式，是宏观结构域内被第m个原型微结构占据的区域，G_m是第m个原型微结构的体积约束，是第m个原型微结构的等效密度，是宏观固体单元的体积，是正则化的设计变量的上边界，是正则化的设计变量,，是正则化的设计变量，的下边界，N_s是形状函数矩阵，K是宏观结构的整体刚度矩阵，M是宏观结构的整体质量矩阵；

(a4)按照下列表达式计算目标函数及约束条件对微观尺度设计变量的敏度，具体如下：

其中，δ(·)是Dirac函数，是紧支撑径向基插值函数。

(a5)采用基于梯度的优化算法更新所述微观尺度设计变量并输出优化后的各原型微结构构型。

3.如权利要求1或2所述的一种基于代理模型的多尺度拓扑优化方法，其特征在于，在步骤(a)中，所述聚类优选按照下列关系式进行：

其中，是密度聚类后的第Ne个宏观单元密度值，表示密度聚类前的第Ne个宏观单元密度值，和分别表示第m-1个和m个原型微结构的等效密度。

4.如权利要求1-3任一项所述的一种基于代理模型的多尺度拓扑优化方法，其特征在于，在步骤(b)中，所述预测模型优选采用Kriging模型。

5.如权利要求1-4任一项所述的一种基于代理模型的多尺度拓扑优化方法，其特征在于，在步骤(c)中，所述优化模型优选按照下列表达式：

其中，是宏观单元相对密度，NE是宏观单元总数，u_k是第k阶特征频率ω_k的特征向量，K是宏观结构的整体刚度矩阵，M是宏观结构的整体质量矩阵，G是给定的材料用量约束条件，是宏观固体单元的体积，Ω^MA是整个宏观结构设计域，是限定的材料使用量，ρ_min是设计变量的最小值，f_pm是通过插值原型微结构得到的非均匀微结构的最小密度，ρ_max是设计变量的最大值。

6.如权利要求1-5任一项所述的一种基于代理模型的多尺度拓扑优化方法，其特征在于，在步骤(d)中，所述判断所述宏观结构的最大特征频率是否收敛是判断前后两次计算获得的最大特征频率的差值是否小于预设阈值，若小于预设阈值则收敛。