一种稳定成孔的改进水平集拓扑优化方法与流程

本发明属于结构拓扑优化的技术领域，更具体地，涉及一种稳定成孔的改进水平集拓扑优化方法。

背景技术

结构拓扑优化是一种新型的数字化结构设计方法，通过建立包含设计变量、目标函数与约束条件的数学模型，对结构进行有限元分析，在设计域中按照优化准则或者数学规划方法迭代出满足约束条件并达到设计目标的材料布局。结构拓扑优化在汽车、航空航天、高端装备等行业有重大应用价值。

水平集拓扑优化方法和beso(双向渐进结构优化)方法是结构拓扑优化中常用的两种方法。水平集拓扑优化方法通过高一维空间的水平集函数的零水平集来表示结构边界，设计变量为结构边界，建立边界与目标函数的关系，通过优化算法演化边界来得到最优结构，beso方法是由eso(渐进结构优化)方法发展而来，通过逐步地去除低应力单元或添加已删除单元实现结构最优化，该方法是一种基于梯度的数学优化方法，具有基于历史迭代信息的敏度过滤和稳定性。

在传统水平集拓扑优化方法中，由于汉密尔顿-雅克比方程的性质，结构内部无法生成孔洞。因此，水平集拓扑优化中结构的初始设计需要有许多孔洞，这导致该方法依赖于初始设计，尤其对于二维问题，为了克服对初始设计的依赖性，水平集方法已经做了很多改进，比如结合拓扑导数或利用径向基函数，但这些改进措施在初始设计中需要在设计域定义很多孔洞，其中孔洞的位置、数量和大小均会对最终的优化结果产生很大的影响；另外，这些改进措施概念复杂，不易进行数值分析，难以应用在实际工程中，目前尚未有针对水平集拓扑优化稳定成孔并在工程应用中可行的有效方法。

技术实现要素：

针对现有技术的以上缺陷或改进需求，本发明提供了一种稳定成孔的改进水平集拓扑优化方法，通过利用两次优化多次迭代的方式，从结构内部去除无效材料的特点稳定地生成孔洞，得到水平集拓扑优化中每一次迭代的初始设计，进而迭代出最优结构，尤其是其中的第一次优化，采用了通过设定敏度数并通过敏度数删除初始赋值中设定的实体区域，以此形成孔洞，即双向渐进结构优化beso的方法，无需对初始孔洞进行复杂的设计，有效地克服了水平集方法对初始设计的依赖性，概念简单，便于实施。

为实现上述目的，按照本发明的一个方面，提供了一种稳定成孔的改进水平集拓扑优化方法，其特征在于，该方法包括下列步骤：

(a)对待处理对象的设计域进行有限元网格划分获得多个网格单元，建立与每个网格单元相一一对应的水平集函数，给每个网格单元的水平集函数赋予初始值，使得所述设计域划分为实体区域和孔洞区域；

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

(c)利用所述设计域的位移向量、弹性矩阵分别计算所述实体和孔洞区域相应的速度场，利用所述设计域的位移向量和所述实体区域的刚度矩阵计算所述实体区域中每个网格单元的敏度数，并根据该敏度数设定敏度数的可接受阈值；

(d)将所述实体区域的每个网格单元对应的敏度数与所述可接受阈值比较，删除所述敏度数小于所述可接受阈值对应的网格单元，即进行第一次优化，由此使得所述设计域的实体区域和网格区域重新划分，同时根据重新划分后获得的实体区域和孔洞区域，更新所述每个网格单元的水平集函数；

(e)利用步骤(c)中获得的所述实体和孔洞区域各自的速度场，分别更新所述第一次优化获得的实体和孔洞区域中每个网格单元的水平集函数，即进行第二次优化，所述更新后的水平集函数对应新的实体区域和孔洞区域，即第二次优化的实体和孔洞区域；

(f)计算所述第二次优化的实体区域的柔度和体积误差，并分别与预设误差阈值进行比较，若所述柔度和体积误差均小于所述误差阈值，则所述第二次优化的实体和孔洞区域为最终所需的结果，否则返回步骤(b)。

进一步优选地，在步骤(b)中，所述计算弹性矩阵和刚度矩阵优选采用下列表达式进行，

其中，d是弹性矩阵，e是弹性模量，ν是泊松比，k^e是刚度矩阵，b为是每个网格单元的应变矩阵，ae是每个网格单元的面积，da是面积微元。

进一步优选地，在步骤(b)中，所述计算设计域的位移向量优选采用下列表达式，

ku＝f

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

进一步优选地，在步骤(c)中，所述计算速度场优选采用下列表达式，

vnormal＝dε(u)·ε(u)-λ

其中，ε(u)＝bu，vnormal是速度场，λ为拉格朗日乘子。

进一步优选地，在步骤(c)中，所述计算敏度数优选采用下列表达式，

其中，αe是实体区域中每个网格单元的敏度数，ue是实体区域中每个网格单元的位移向量，是实体区域的刚度矩阵。

进一步优选地，在步骤(c)中，所述设定敏度数的可接受阈值优选采用下列方法：

(c1)将所述实体区域所有网格单元的敏度数按照从小到大的顺序排列，预设实体区域的删除率，按照下列表达式计算第一个试探临界值

其中，p是删除率，ns是实体区域中网格单元的总数量；

(c2)计算实体区域中与孔洞区域交界处的网格单元的敏度数的平均值预设小于1的比例参数μ，按照下列表达式计算另一个试探临界值

(c3)取所述试探临界值和中的较小值作为敏度数的可接受阈值。

进一步优选地，在步骤(e)中，所述更新所述第一次优化获得的实体和孔洞区域中每个网格单元的水平集函数，优选采用下列表达式进行更新，

其中，φnew(x)是第一次优化获得的水平集函数，t是时间。

进一步优选地，在步骤(f)中，所述柔度和体积误差优选按照下列表达式计算，

其中，cerr、verr分别表示柔度和体积的误差，c^k-i+1、c^k-5-i+1分别表示第(k-i+1)次和第(k-5-i+1)次迭代的柔度，v^k为第k次迭代的实体区域的体积，为最终优化后实体区域的总体积，i是设定的计算次数，取1～5之间的任意正整数，k是迭代次数。

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

1、本发明中通过采用双向渐进结构化方法(beso)，通过利用初始赋值中设定的实体单元的每个网格单元对应的敏度数，设定敏度数的可接受阈值，然后利用该可接受阈值删除初始赋值中设定的实体区域，以此形成新的孔洞和实体区域，现有技术中，最终的优化结果与初始设计中孔洞和实体区域的分布息息相关，但是本发明中对初始中孔洞和实体区域可任意设置，且对最终优化结果没有影响，即有效地克服了水平集方法对初始设计的依赖性；

2、本发明提供了一种稳定成孔的改进水平集拓扑优化方法，其利用beso方法在结构内部稳定地生成孔洞，采用两次优化的方式获得最终所需的被划分的结构，使得最终的优化结构，可靠性高，实用性强；

3、本发明提供的方法简单，可操作性强，具有工程应用的可行性，便于理解和实施，适用于所有受到载荷的力学结构，适用范围广，对工程中后期的结构设计具有指导性的作用。

附图说明

图1是按照本发明的优选实施例所构建的方法流程图；

图2是按照本发明的优选实施例所构建的待优化的平面悬臂梁结构示意图；

图3是按照本发明的优选实施例所构建的平面悬臂梁水平集拓扑优化中的初始设计；

图4是按照本发明的优选实施例所构建的平面悬臂梁水平集拓扑优化结果；

图5是按照本发明的优选实施例所构建的平面悬臂梁水平集拓扑优化目标函数和体积比的收敛曲线；

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

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

具体实施方式

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

图1是按照本发明的优选实施例所构建的方法流程图，如图1所示，一种稳定成孔的改进水平集拓扑优化方法，该方法包括以下步骤：

图2是按照本发明的优选实施例所构建的待优化的平面悬臂梁结构示意图，如图2所示，以带有集中载荷的悬臂梁柔度最小化的优化问题为例来进一步说明本发明。给定2m×1m矩形设计域d，区域左端固定，区域右端中点处施加竖直向下的集中载荷1n。对所述悬臂梁的结构进行优化设计，使其柔度最小化，即刚度最大化。

s1选用四节点正方形单元对矩形设计域d划分有限元网格，网格数量为160×80，网格的边长为6.25mm，节点数量为161×81，每个网格的单元中心点坐标为x，定义水平集函数φ(x)如下，按照每个网格的中心点的坐标，赋值给每个网格单元对应的水平集函数φ(x)，以此将设计域划分为实体区域和孔洞区域，本实施例中φ(x)的取值范围是[-3,3]，

s2根据悬臂梁所用材料的属性，设定实体区域的弹性模量es为1pa和泊松比νs为0.3，孔洞区域的弹性模量应趋于0以及泊松比应与实体区域一致，则设定孔洞区域的弹性模量ev为1×10^-3pa和泊松比νv为0.3，由此建立各个网格单元的刚度矩阵k^e，其计算式如下：

式中，b为每个网格单元的应变矩阵，通过每个网格单元的4个节点的坐标和该单元的边长，可获得该网格单元对应的形函数和位移函数，利用每个网格单元的形函数和位移函数计算获得每个网格单元的应变矩阵b，应变矩阵b的计算方法为现有技术，此处不再累述，d为弹性矩阵，da为面积微元，ae为网格单元e的面积，弹性矩阵d计算如下：

将es和νs代入上式可得实体区域的弹性矩阵ds，利用刚度矩阵的计算方法计算获得实体区域中每个网格单元对应的刚度矩阵，实体区域所有网格单元的刚度矩阵组合形成实体区域的刚度矩阵同样地，将ev和νv代入上式可得孔洞区域的弹性矩阵dv，进而得到将实体区域和孔洞区域的刚度矩阵和组装后形成得到设计域整体刚度矩阵k，该组装方法为现有技术，在此不再赘述，根据ku＝f计算得到设计域整体位移向量u，其中，f为设计域受到的外力向量，本实施例的外力向量即为设计域右端中点处施加的竖直向下的1n的集中载荷；

s3根据设计域的位移向量获得实体区域、孔洞区域以及所有网格单元中每个网格单元各自对应的位移向量；利用实体和孔洞区域各自的弹性矩阵和位移向量分别计算相应的速度场vnormal；利用实体区域中每个网格单元的位移向量和实体区域的刚度矩阵计算实体区域中每个网格单元的敏度数αe，根据该敏度数设定敏度数的可接受阈值；

整体位移向量u中每个元素代表一个网格单元的位移向量，从步骤s2获得的整体位移向量u获知每个网格单元的位移向量，实体区域包括的网格单元的位移向量形成实体区域的位移向量，同理获得孔洞区域的位移向量，

通过形状导数分析计算速度场vnormal，vnormal的方向沿实体区域与孔洞区域的交界线的向孔洞方向的法向方向，其计算式如下：

vnormal＝dε(u)·ε(u)-λ

式中，ε(u)＝bu为应变张量，λ为拉格朗日乘子(初始值设定为10^-3)；

利用实体区域中的每个网格单元的位移向量ue，计算全部实体区域中每个网格单元的敏度数αe，其计算公式如下：

按照从小到大的升序方式给全部实体区域的敏度数排序，可得ns为当前结构的实体区域数，设定实体区域的删除率为0.04，排序中第0.04ns位的敏度数记为将第0.04ns位之前的实体区域删除，计算αe的一个试探临界值

φ(x)＝0时对应的实体区域与孔洞区域交界处的实体区域，计算该交界处的实体区域敏度数的平均值则可得αe的另一个试探临界值μ为设定的参数，需满足0≤μ≤1，这里设置μ＝0.5；比较两个试探敏度值取两者间的小值作为αe的敏度数可接受阈值αth；

s4比较全部实体区域的敏度数αe和其阈值αth，当满足αe＜αth时，去除该实体区域，即插入孔洞，获得新的孔洞和实体区域的设计域，新生成的孔洞大小为删除的网格单元的大小，根据设计域中新的孔洞和实体区域的分布，重新赋值给水平集函数，由此得到新的水平集函数φnew(x)；

s5根据实体和孔洞区域各自的速度场，求解汉密尔顿-雅克比方程即再次更新水平集函数φnew(x)，也即将其重新初始化使得实体区域与孔洞区域的交界处平滑过渡，避免水平集函数的梯度过大或过小，具体过程与传统水平集拓扑优化方法的相关步骤一致，在此不赘述，根据再次更新获得的水平集函数，设计域也被再次更新划分为新的实体和孔洞区域；

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

c＝f^tu

其中，cerr、verr分别表示实体区域柔度和体积的误差，c^k-i+1、c^k-5-i+1分别表示第(k-i+1)次和第(k-5-i+1)次迭代的柔度(如何计算)，柔度c＝f^tu，f为设计域受到的外力向量，u是整体位移向量，从步骤s2到s5记为一次迭代，v^k为第k次迭代的实体区域的体积，为优化后实体区域结构的体积上限，本实施例中取为结构初始体积的1/2，若收敛则优化过程结束，若不收敛则重复s2到s5直到收敛，i是设定的次数，取1～5之间的任意正整数，k是迭代次数。。

图3是按照本发明的优选实施例所构建的平面悬臂梁水平集拓扑优化中的初始设计，如图3所示，图中将设计域的上下端设置为孔洞区域，即图中的白色部分，设计域的中间部分为实体区域，见图中的灰色部分，图4是按照本发明的优选实施例所构建的平面悬臂梁水平集拓扑优化结果，如图4所示，获得图中的最优实体和孔洞区域后，可知在现有的载荷下，该悬臂梁在结构设计中哪些区域适合作为实体区域，哪些区域适合设计为孔洞区域。

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

图6是现有技术中水平集拓扑优化中的初始设计；图7是现有技术水平集拓扑优化中的优化结果，如图6和7所示，可以看出，本发明可以克服传统水平集拓扑优化方法对初始设计的依赖，获得基本相同的优化结果。

本发明提供的一种在水平集拓扑优化中利用双向渐进结构优化稳定成孔的方法，其利用beso方法从结构内部去除无效材料的特点稳定地生成孔洞，得到水平集拓扑优化中每一次迭代的初始设计，进而迭代出最优结构，克服了传统水平集拓扑优化方法对初始设计的依赖，具有工程应用的可行性，便于理解和实施。

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