CN110321611A - 一种多材料结构拓扑优化方法 - Google Patents

Abstract

本发明属于结构优化设计领域，并具体公开了一种多材料结构拓扑优化方法。该方法包括：对待处理对象的设计域进行有限元网格划分获得预设数量的网格单元并建立与实体材料一一对应的水平集函数；计算弹性矩阵、刚度矩阵和位移向量；计算速度场和敏度数，并根据敏度数对网格单元中的实体材料进行转换或删除，获得当前迭代下第一次优化后的水平集函数；利用速度场，更新获得当前迭代下第二次优化后的水平集函数；根据柔度和体积误差判别第二次优化的结果是否收敛。本发明克服了经典水平集拓扑优化方法对初始设计的依赖性，以及双向渐进优化方法无法获得光滑边界的缺陷，并在一定程度上提高了结构刚度，提供的方法稳定并且有效。

Description

一种多材料结构拓扑优化方法

技术领域

本发明属于结构优化设计领域，更具体地，涉及一种多材料结构拓扑优化方法。

背景技术

多材料结构的拓扑优化是指在规定的空间域内将具有不同力学特性的多种材料合理地配置在一起，以获得优异性能的结构。近年来，随着增材制造技术的出现，多材料结构的制造成本相对降低，因此这样的结构在工程应用中受到越来越多的关注。

多材料结构拓扑优化的设计方法中，水平集方法利用多个水平集函数表示被不同材料占据的区域以及它们之间的交界，再进行交界的演化，优化得到结构的拓扑与形状；双向渐近优化(BESO)方法通过逐渐移除和添加多种材料来实现结构中材料的最优配置。

利用经典水平集方法解决多材料结构优化问题，优化过程中的拓扑变化包括孔洞的消失、生成、合并，区域的分离，以及不同材料之间的转换，这样的材料转换被称为广义成孔，即孔洞不再是空区域，而成为一种由假想实体材料占据的区域。然而，在经典水平集拓扑优化方法中，结构内部无法生成孔洞，这导致二维问题的最优解依赖于初始设计，三维拓扑优化的灵活性会受到不利影响。BESO方法的概念简单，易于数值实现，便于工程人员使用，但其缺陷在于优化得到的结果边界不光滑。以上两种方法在多材料结构拓扑优化中均存在不足之处，无法有效地解决实际的工程问题。

发明内容

针对现有技术的上述缺点和/或改进需求，本发明提供了一种多材料结构拓扑优化方法，其中将BESO方法的材料转换格式集成到基于水平集的多材料结构拓扑优化中，充分利用了BESO方法便于广义成孔的特点，有效地克服了经典水平集方法对初始设计的依赖性，保证了最优结构的边界光滑。

为实现上述目的，本发明提出了一种多材料结构拓扑优化方法，该方法包括如下步骤：

(a)对待处理对象的设计域进行有限元网格划分获得预设数量的网格单元，利用N种实体材料和一种弱材料占据每个所述网格单元，其中被所述实体材料占据的区域定义为实体区域，被所述弱材料占据的区域定义为孔洞区域，然后根据所述设计域中实体区域和孔洞区域的分布，建立与所述实体材料一一对应的水平集函数；

(b)根据待处理对象的材料属性，分别设定所述实体材料和弱材料的弹性模量和泊松比，以此计算所述实体材料和弱材料对应的弹性矩阵，并获得所述实体区域中各网格单元对应的刚度矩阵，然后计算所述设计域的位移向量；

(c)利用所述实体材料的弹性矩阵和设计域的位移向量计算速度场，并利用所述网格单元的刚度矩阵和位移向量计算所述实体区域内各网格单元的敏度数，根据所述实体区域中各网格单元的敏度数对该网格单元中的实体材料进行转换或删除，从而更新所述实体区域和孔洞区域的分布，由此获得当前迭代下第一次优化后的水平集函数；

(d)利用步骤(c)中获得的所述速度场，更新当前迭代下所述第一次优化后的水平集函数，获得当前迭代下第二次优化后的水平集函数；

(e)计算当前迭代下所述第二次优化后的实体区域的柔度和体积误差，并分别与其预设误差阈值进行比较，判断所述柔度和体积误差是否均小于其所述预设误差阈值，若是，则将当前迭代下所述第二次优化后的水平集函数作为最优结构输出，若否，则返回步骤(b)。

作为进一步优选地，在步骤(a)和步骤(d)中，根据网格单元节点坐标到区域边界的距离以及符号确定所述水平集函数，其中所述符号的确定方式为：

式中，Φ_i(x)为水平集函数，为设计域中任意一点的坐标x，为第i种材料占据的实体区域i，为实体区域i的区域边界，为设计域D内除实体区域i及其区域边界的其他区域。

作为进一步优选地，在步骤(b)中，采用下式计算所述实体材料和弱材料对应的弹性矩阵：

式中，i为材料的编号，其中i取1～N时为实体材料，i取N+1时为弱材料，D_i为第i种材料的弹性矩阵，E_i为第i种材料的弹性模量，ν_i为第i种材料的泊松比；

同时，采用下式计算每个网格单元对应的刚度矩阵：

式中，K^e为网格单元的刚度矩阵，A_e为网格单元的面积，B为应变矩阵，dA为面积的微分；

此外，采用下式计算所述设计域的位移向量：

Ku＝f

式中，K是设计域的刚度矩阵，f是设计域受到的外力向量，u是设计域的位移向量。

作为进一步优选地，在步骤(c)中，采用下式计算速度场：

式中，V_n(x)为速度场，ε(u)为应变张量，λ_i、λ_j为拉格朗日乘子，为实体区域i与孔洞区域的交界，Γ_ij为实体区域i与实体区域j的交界；

同时，采用下式计算所述实体区域内各单元的敏度数：

式中，α_e为实体区域中网格单元的敏度数，u_e为实体区域中网格单元的位移向量，为实体区域i中网格单元的刚度矩阵。

作为进一步优选地，在步骤(c)中，根据所述实体区域中各网格单元的敏度数对其实体材料进行转换或删除的具体过程为：

将实体区域i中所有网格单元的敏度数按照从小到大的顺序排列，预设该实体区域i的进化体积比c_i和第i种实体材料的体积比上限并计算所述实体区域i的总体积V_i，当满足时，将实体区域i中敏度数在第c_i×N_Si位之前的网格单元的实体材料由第i种实体材料转换为第i+1种实体材料，N_Si为实体区域i的单元总数，当i＝N时，将该实体材料删除转换为所述弱材料。

作为进一步优选地，在步骤(d)中，采用下式对当前迭代下所述第一次优化后的水平集函数进行更新：

式中，为当前迭代下第一次优化后的水平集函数，t为时间，V_n为速度场。

作为进一步优选地，在步骤(e)中，采用下式对所述柔度的误差进行计算：

式中，C_err为柔度的误差，n为预设的计算次数，k为当前迭代次数，C^k-n+1为第(k-n+1)次迭代的柔度，C^k-5-n+1为第(k-5-n+1)次迭代的柔度；

采用下式对所述体积的误差进行计算：

式中，V_ierr为体积的误差，V_i ^k为第k次迭代的实体区域i的体积，为第i种实体材料的体积比上限。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，主要具备以下的技术优点：

1.本发明将BESO方法的材料转换格式集成到基于水平集的多材料结构拓扑优化方法中，充分利用了BESO方法便于广义成孔的特点，有效地克服了经典水平集方法对初始设计的依赖性，同时利用水平集方法的特点可得到具有光滑边界的最优结构，便于设计和制造；

2.尤其是，本发明通过对实体材料的转换和删除方法进行优化，能够获得较佳的优化效果，提高优化效率；

3.本发明提供的方法有工程应用的可行性，适用于所有受到载荷的力学结构，适用范围广，对多材料结构设计具有指导性作用。

附图说明

图1是本发明实施例提供的多材料结构拓扑优化方法的流程图；

图2是按照本发明优选实施例构建的待优化的桥型结构示意图；

图3是按照本发明优选实施例构建的桥型结构的初始设计示意图；

图4是按照本发明优选实施例构建的桥型结构的拓扑优化结果；

图5是按照本发明优选实施例构建的桥型结构的优化收敛曲线；

图6是现有技术中经典水平集拓扑优化的初始设计；

图7是现有技术中经典水平集拓扑优化的优化结果。

具体实施方式

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

如图1所示，本发明实施例提供了一种多材料结构拓扑优化方法，该方法包括如下步骤：

(a)对待处理对象的设计域进行有限元网格划分获得预设数量的网格单元，利用N种实体材料和一种弱材料占据每个网格单元，其中被实体材料占据的区域定义为实体区域，被弱材料占据的区域定义为孔洞区域，然后根据设计域中实体区域和孔洞区域的分布，建立与实体材料一一对应的水平集函数；

(b)根据待处理对象的材料属性，分别设定实体材料和弱材料的弹性模量和泊松比，以此计算实体材料和弱材料对应的弹性矩阵，并获得实体区域中各网格单元对应的刚度矩阵，然后计算设计域的位移向量；

其中，采用式(1)计算实体材料和弱材料对应的弹性矩阵：

式中，i为材料的编号，其中i取1～N时为实体材料，i取N+1时为弱材料，D_i为第i种材料的弹性矩阵，E_i为第i种材料的弹性模量，ν_i为第i种材料的泊松比；

采用式(2)计算每个网格单元对应的刚度矩阵：

式中，K^e为网格单元的刚度矩阵，A_e为网格单元的面积，B为应变矩阵，dA为面积的微分；

采用式(3)计算设计域的位移向量：

Ku＝f (3)

式中，K是设计域的刚度矩阵，f是设计域受到的外力向量，u是设计域的位移向量；

(c)利用实体材料的弹性矩阵和设计域的位移向量计算速度场，并利用网格单元的刚度矩阵和位移向量计算实体区域内各网格单元的敏度数，根据实体区域中各网格单元的敏度数对该网格单元中的实体材料进行转换或删除，从而更新实体区域和孔洞区域的分布，由此获得当前迭代下第一次优化后的水平集函数；

其中，采用式(4)计算速度场：

式中，V_n(x)为速度场，ε(u)为应变张量，λ_i、λ_j为拉格朗日乘子，为实体区域i与孔洞区域的交界，Γ_ij为实体区域i与实体区域j的交界；

采用式(5)计算实体区域内各单元的敏度数：

式中，α_e为实体区域中网格单元的敏度数，u_e为实体区域中网格单元的位移向量，为实体区域i中网格单元的刚度矩阵；

根据实体区域中各网格单元的敏度数对其实体材料进行转换或删除的具体过程为：

将实体区域i中所有网格单元的敏度数按照从小到大的顺序排列，预设该实体区域i的进化体积比c_i和第i种实体材料的体积比上限并计算实体区域i的总体积V_i，当满足时，将实体区域i中敏度数在第c_i×N_Si位之前的网格单元的实体材料由第i种实体材料转换为第i+1种实体材料，N_Si为实体区域i的单元总数，当i＝N时，将该实体材料删除转换为弱材料；

(d)利用步骤(c)中获得的速度场，采用式(6)对当前迭代下第一次优化后的水平集函数进行更新，获得当前迭代下第二次优化后的水平集函数；

式中，为当前迭代下第一次优化后的水平集函数，t为时间，V_n为速度场；

(e)计算当前迭代下第二次优化后的实体区域的柔度和体积误差，并分别与其预设误差阈值进行比较，判断柔度和体积误差是否均小于其预设误差阈值，若是，则将当前迭代下第二次优化后的水平集函数作为最优结构输出，若否，则返回步骤(b)，根据当前迭代下第二次优化后的水平集函数重复步骤(b)～(e)；

更具体地，采用式(7)对柔度的误差进行计算：

式中，C_err为柔度的误差，n为预设的计算次数，k为当前迭代次数，C^k-n+1为第(k-n+1)次迭代的柔度，C^k-5-n+1为第(k-5-n+1)次迭代的柔度；

采用式(8)对体积的误差进行计算：

式中，V_ierr为体积的误差，V_i ^k为第k次迭代的实体区域i的体积，为第i种实体材料的体积比上限。

进一步，在步骤(a)和步骤(d)中，根据网格单元节点坐标到区域边界的距离以及符号确定水平集函数，其中符号的确定方式为：

式中，Φ_i(x)为水平集函数，为设计域中任意一点的坐标x，为第i种材料占据的实体区域i，为实体区域i的区域边界，为设计域D内除实体区域i及其区域边界的其他区域。

按照本发明的一个优选实施例，利用本发明提供的方法优化桥型结构，图2是待优化的桥型结构，以带有集中载荷的结构柔度最小化的优化问题为例对本发明作进一步说明，给定2m×1m矩形设计域D，区域左下角固定，右下角只被限制竖直方向的自由度，下边中点处施加竖直向下的集中载荷1N，对所述桥型结构进行优化设计，使其柔度最小化，即刚度最大化。

(a)由于结构的对称性，将矩形设计域D的左半部分(或右半部分)划分为有限个网格单元，用于有限元分析，单元类型选为四节点双线性正方形单元，单元边长定为12.5mm，则单元数量为80×80；同时，在矩形设计域D的左半部分(或右半部分)设置结构化网格，用于水平集计算，网格节点数量为81×81，有限元的节点与水平集网格的节点重合，计算水平集网格的节点坐标x，根据网格节点坐标到区域边界的距离以及符号判定，定义多个水平集函数Φ_i(x)，其符号确定如下：

利用两种实体材料和一种弱材料占据每个网格单元，故i＝1,2；

(b)根据桥型结构所用材料的属性，设定第1种实体材料的弹性模量E₁为5Pa和泊松比ν₁为0.3，第2种实体材料的弹性模量E₂为1Pa和泊松比ν₂为0.3，假设用于占据孔洞区域的弱材料的弹性模量E₃为1×10^-3Pa和泊松比ν₃为0.3，被实体材料占据的区域被统称为实体区域，被第1种实体材料占据的区域被统称为实体区域一，被第2种实体材料占据的区域被统称为实体区域二，被弱材料占据的区域被统称为孔洞区域，由此建立网格单元的刚度矩阵K^e，其计算式如下：

式中，B为网格单元的应变矩阵，根据每个单元的节点坐标和该单元的边长，可获得该单元对应的形函数和位移函数，利用每个单元的形函数和位移函数计算获得每个单元的应变矩阵B，应变矩阵B的计算方法为现有技术，此处不再累述，D_i为第i种材料的弹性矩阵，dA为面积的微分，A_e为网格单元的面积，弹性矩阵D_i计算如下：

将E₁和ν₁代入上式可得第1种实体材料的弹性矩阵D₁，利用刚度矩阵的计算方法计算获得实体区域一中每个网格单元对应的刚度矩阵同样地，将E₂和ν₂代入上式可得第二种实体材料的弹性矩阵D₂，进而得到将E₃和ν₃代入上式可得人为假想的弱材料的弹性矩阵D₃，进而得到当单元被不同区域分割时，计算单元被第1种实体材料、第二种实体材料及弱材料占据的区域面积与单元总面积的比值分别为Ar₁、Ar₂和Ar₃，则被两种实体材料占据的网格单元对应的刚度矩阵被第1种实体材料和弱材料占据的网格单元对应的刚度矩阵而被第二种实体材料和弱材料占据的网格单元对应的刚度矩阵同时被三种材料占据的网格单元对应的刚度矩阵将全部的单元刚度矩阵 组装后形成得到设计域整体刚度矩阵K，该组装方法为现有技术，在此不再赘述，根据Ku＝f计算得到设计域的整体位移向量u，其中，f为设计域受到的外力向量，本实施例的外力向量即为设计域下边中点处施加的竖直向下的1N的集中载荷；

(c)根据设计域的整体位移向量u每两个元素代表一个单元的位移向量，获得设计域内各网格单元的位移向量u_e；利用不同材料的弹性矩阵D₁、D₂、D₃和整体位移向量u计算设计域内的速度场V_n；利用实体区域中每个单元的位移向量u_e和实体区域的刚度矩阵计算实体区域中每个单元的敏度数α_e，根据该敏度数设定材料转换的敏度数阈值；

通过形状导数分析计算速度场V_n，其计算式如下：

式中，ε(u)为应变张量，其计算式为ε(u)＝Bu，λ₁、λ₂均为拉格朗日乘子(初始值均设定为3×10^-3)，表示被第i种材料占据区域与孔洞区域的交界，Γ₁₂表示实体区域一与实体区域二的交界；

利用实体区域中每个单元的位移向量u_e，计算各实体单元的敏度数α_e，其计算公式如下：

网格单元敏度数α_e需进行过滤处理和平均处理得到现有的BESO方法中已详述这一过程，在此不再赘述，然后按照从小到大的升序方式给实体区域一的单元敏度数排序，可得N_S1为当前实体区域一的单元总数，设定实体区域一的进化体积比c₁为0.04，排序中第0.04N_S1位的敏度数记为α_N，计算实体区域一的单元敏度数临界值同时设定第1种实体材料的体积比为0.5，计算实体区域一的单元总体积V₁，当满足时，将敏度数小于的单元材料转换为第2种材料；

同样地，按照从小到大的升序方式给实体区域二的单元敏度数排序，N_S2为当前实体区域二的单元总数，设定实体区域二的进化体积比c₂为0.04，将敏度数在第0.04N_S2位之前的单元删除，即插入孔洞，获得包含新的实体区域和孔洞区域的设计域，新生成的孔洞半径为水平集网格的边长，其值与有限单元边长相等；

根据设计域中新的实体区域和孔洞区域的分布，重新计算网格节点坐标到区域边界的距离，其符号的判定同步骤(a)，获得当前迭代下第一次优化的

(d)根据步骤(c)中获得的速度场V_n，求解汉密尔顿-雅克比方程即再次更新获得当前迭代下第二次优化的水平集函数具体过程与经典水平集拓扑优化方法的相关步骤一致，在此不赘述，根据再次更新获得的水平集函数，设计域也被再次更新划分为新的实体区域一、实体区域二和孔洞区域；

(e)判断优化结果是否收敛，收敛准则表示如下：

其中，C_err、V_1err、V_2err分别表示结构柔度误差、第1种材料的体积误差和第2种材料的体积误差，C^k-n+1、C^k-5-n+1分别表示第(k-n+1)次和第(k-5-n+1)次迭代的柔度，柔度C＝f^Tu，f为设计域受到的外力向量，u是整体位移向量，从步骤(b)到(e)记为一次迭代，V₁ ^k为第k次迭代的实体区域一的体积，为人为设定的实体区域一的体积上限，为第k次迭代的实体区域一的体积，为人为设定的实体区域二的体积上限，若收敛则优化过程结束，若不收敛则重复(b)到(e)直到收敛，n是设定的次数，取1～5之间的任意正整数，k是当前迭代次数。

图3是按照本发明的优选实施例所构建的桥型结构的初始设计，如图3所示，设计域由第1种材料(灰色)全部占据，图4是按照本发明的优选实施例所构建的桥型结构的拓扑优化结果，如图4所示，获得图中的最优实体和孔洞区域后，可知在现有的载荷下，该桥型结构在结构设计中哪些区域适合被第1种材料(灰色)占据，哪些区域适合被第2种材料(黑色)占据，哪些区域适合设计为孔洞区域(白色)。

图5是按照本发明的优选实施例所构建的桥型结构在拓扑优化过程中的目标函数(柔度值)和体积比的收敛曲线，如图5所示，最终目标函数和体积比趋于稳定，即均收敛，即采用本发明提供的方法获得的最终的优化结果是最优划分结构。

图6是现有技术中经典水平集拓扑优化中的初始设计；图7是现有技术经典水平集拓扑优化中的优化结果，如图6和7所示，可以看出，本发明可以克服经典水平集拓扑优化方法对初始设计的依赖，获得更优的优化结果，最终的结构柔度值减小了7.4％，即结构刚度更大。

本发明提供的一种多材料结构拓扑优化方法，将BESO方法的材料转换格式集成到基于水平集的多材料结构拓扑优化中，充分利用了BESO方法便于广义成孔的特点，有效地克服了经典水平集方法对初始设计的依赖性，保证了最优结构的边界光滑，并在一定程度上提高了结构刚度。

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

Claims (7)

1.一种多材料结构拓扑优化方法，其特征在于，该方法包括如下步骤：

(a)对待处理对象的设计域进行有限元网格划分获得预设数量的网格单元，利用N种实体材料和一种弱材料占据每个所述网格单元，其中被所述实体材料占据的区域定义为实体区域，被所述弱材料占据的区域定义为孔洞区域，然后根据所述设计域中实体区域和孔洞区域的分布，建立与所述实体材料一一对应的水平集函数；

(b)根据待处理对象的材料属性，分别设定所述实体材料和弱材料的弹性模量和泊松比，以此计算所述实体材料和弱材料对应的弹性矩阵，并获得所述实体区域中各网格单元对应的刚度矩阵，然后计算所述设计域的位移向量；

(c)利用所述实体材料的弹性矩阵和设计域的位移向量计算速度场，并利用所述网格单元的刚度矩阵和位移向量计算所述实体区域内各网格单元的敏度数，根据所述实体区域中各网格单元的敏度数对该网格单元中的实体材料进行转换或删除，从而更新所述实体区域和孔洞区域的分布，由此获得当前迭代下第一次优化后的水平集函数；

(d)利用步骤(c)中获得的所述速度场，更新当前迭代下所述第一次优化后的水平集函数，获得当前迭代下第二次优化后的水平集函数；

(e)计算当前迭代下所述第二次优化后的实体区域的柔度和体积误差，并分别与其预设误差阈值进行比较，判断所述柔度和体积误差是否均小于其所述预设误差阈值，若是，则将当前迭代下所述第二次优化后的水平集函数作为最优结构输出，若否，则返回步骤(b)。

2.如权利要求1所述的多材料结构拓扑优化方法，其特征在于，在步骤(a)和步骤(d)中，根据网格单元节点坐标到区域边界的距离以及符号确定所述水平集函数，其中所述符号的确定方式为：

式中，Φ_i(x)为水平集函数，为设计域中任意一点的坐标x，为第i种材料占据的实体区域i，为实体区域i的区域边界，为设计域D内除实体区域i及其区域边界的其他区域。

3.如权利要求1或2所述的多材料结构拓扑优化方法，其特征在于，在步骤(b)中，采用下式计算所述实体材料和弱材料对应的弹性矩阵：

式中，i为材料的编号，其中i取1～N时为实体材料，i取N+1时为弱材料，D_i为第i种材料的弹性矩阵，E_i为第i种材料的弹性模量，ν_i为第i种材料的泊松比；

同时，采用下式计算每个网格单元对应的刚度矩阵：

式中，K^e为网格单元的刚度矩阵，A_e为网格单元的面积，B为应变矩阵，dA为面积的微分；

此外，采用下式计算所述设计域的位移向量：

Ku＝f

式中，K是设计域的刚度矩阵，f是设计域受到的外力向量，u是设计域的位移向量。

4.如权利要求1所述的多材料结构拓扑优化方法，其特征在于，在步骤(c)中，采用下式计算速度场：

式中，V_n(x)为速度场，ε(u)为应变张量，λ_i、λ_j为拉格朗日乘子，为实体区域i与孔洞区域的交界，Γ_ij为实体区域i与实体区域j的交界；

同时，采用下式计算所述实体区域内各单元的敏度数：

式中，α_e为实体区域中网格单元的敏度数，u_e为实体区域中网格单元的位移向量，为实体区域i中网格单元的刚度矩阵。

5.如权利要求1所述的多材料结构拓扑优化方法，其特征在于，在步骤(c)中，根据所述实体区域中各网格单元的敏度数对其实体材料进行转换或删除的具体过程为：

将实体区域i中所有网格单元的敏度数按照从小到大的顺序排列，预设该实体区域i的进化体积比c_i和第i种实体材料的体积比上限并计算所述实体区域i的总体积V_i，当满足时，将实体区域i中敏度数在第c_i×N_Si位之前的网格单元的实体材料由第i种实体材料转换为第i+1种实体材料，N_Si为实体区域i的单元总数，当i＝N时，将该实体材料删除转换为所述弱材料。

6.如权利要求1所述的多材料结构拓扑优化方法，其特征在于，在步骤(d)中，采用下式对当前迭代下所述第一次优化后的水平集函数进行更新：

式中，为当前迭代下第一次优化后的水平集函数，t为时间，V_n为速度场。

7.如权利要求1～6任一项所述的多材料结构拓扑优化方法，其特征在于，在步骤(e)中，采用下式对所述柔度的误差进行计算：

式中，C_err为柔度的误差，n为预设的计算次数，k为当前迭代次数，C^k-n+1为第(k-n+1)次迭代的柔度，C^k-5-n+1为第(k-5-n+1)次迭代的柔度；

采用下式对所述体积的误差进行计算：

式中，V_ierr为体积的误差，V_i ^k为第k次迭代的实体区域i的体积，为第i种实体材料的体积比上限。