CN111523270B - 一种改进的连续体结构拓扑优化后处理方法 - Google Patents

Abstract

本发明公开了一种改进的连续体结构拓扑优化后处理方法，其特征在于：所述一种改进的连续体结构拓扑优化后处理方法采用较大的网格尺寸，即有限元分析单元尺寸大于默认值1，对目标模型进行拓扑优化时可有效避免因网格数目过多引起的细小柱状结构；对优化模型进行二值化处理并分析优化模型且提取锯齿形边界线，通过数值计算得到目标离散角点集，将其作为样本点进行曲线拟合，对得到的新边界线进一步插值处理，生成边界光滑的拓扑结构，并对其进行实体化处理。本发明提出的方法有效消除了优化结构的锯齿边界，保证优化结构的原始设计要求，降低了计算难度，减少了计算时长，增强了模型可制造性，并且有效降低了拓扑优化后处理工序的难度。

Description

一种改进的连续体结构拓扑优化后处理方法

技术领域

本发明涉及结构优化相关技术领域，特别涉及一种改进的连续体结构拓扑优化后处理方法。

背景技术

拓扑优化方法是最具前景的结构优化方法之一，该方法能够有效的利用材料，弥补尺寸、形状优化无法改变结构拓扑形式的局限性。因其在优化过程中产生和消除小空腔的灵活性，在近些年中得到特别关注和迅速发展。

连续体结构的拓扑优化是优化设计的一个重要问题，目前最常用的连续体结构拓扑优化方法有均匀化方法、变密度法、渐进结构优化方法(ESO)、水平集法(Levelset)、独立连续映射方法(ICM)等。

拓扑优化结构由于依托网格进行有限元分析，优化结果采用单元体的有无来表征连续体结构中的材料特征，导致结果呈锯齿形。若使用更多的单元数目，理论上得到边界描述详细且结果理想的有限元模型，同时避免不同的优化结果的产生。但实际上，由于网络依赖性及数值不稳定等问题，更多的单元数目将产生无穷多的设计变量，使优化结果中柱状结构的尺寸变小、数量增大，增加了几何复杂性，同时增加优化时长。与此同时，锯齿状的结构边界，使得拓扑优化的结果变得极难加工，且不便于导入其他CAD软件中进行优化后处理。为得到光滑边界，需进一步进行形状优化，若采用映射关系将栅格模型转化为光滑的曲线或曲面模型，映射关系的处理难度通常较大，导致实现结果费时费力，不甚理想。

发明内容

本发明的目的在于克服现有技术的缺点与不足，提出一种改进的连续体结构拓扑优化后处理方法，该方法采用较大的网格尺寸，即有限元分析单元尺寸大于默认值1，对目标模型进行拓扑优化时可有效避免因网格数目过多引起的细小柱状结构；对优化模型进行二值化处理并分析优化模型且提取锯齿形边界线，通过数值计算得到目标离散角点集，将其作为样本点进行曲线拟合，有效避免后处理模型与原始拓扑优化模型产生较大的偏差，影响结构性能；对得到的新边界线进一步插值处理，生成边界光滑的拓扑结构，并对其进行实体化处理。该方法有效保证优化结构原始设计要求，并避免网格依赖性、棋盘格现象。

为达成上述目的，本发明采用如下技术方案，一种改进的连续体结构拓扑优化后处理方法，包括如下步骤：

S1，对目标模型进行拓扑优化处理得到拓扑优化结构。

S2，对得到的拓扑优化结果进行分析提取边界线，并提取锯齿形边界线的角点得到点集P。

S3，对步骤S2得到的离散角点进行光顺连接。

S4，对得到的曲线封闭区域进行裁剪处理，得到光顺的曲线封闭区域，进而生成边界光滑的拓扑结构。

S5，检测生成的拓扑结构是否满足要求，满足执行步骤6，否则重设参数重复执行步骤S3至步骤S4。

S6，得到最终的拓扑结构，对其进行实体化处理。

优选的，所述步骤S1中连续体结构拓扑优化方法采用如变密度法、渐进结构优化方法(ESO)、独立连续映射方法(ICM)等。

优选的，所述步骤S2中对拓扑优化结果进行二值化处理及均值滤波处理，并将二值化结果转化为0-1的数值矩阵。将得到的边界提取结果记为Z(x,y)，对锯齿形边界各角点进行识别计算。设定坐标原点，以预先划分的离散单元大小为坐标单元尺寸，进而确定各角点坐标定义为P_j(x_j,y_j)，构成点集P＝{P₁,P₂,···,P_n}。

优选的，所述步骤S3中由下式判断各直线段长度：

其中，为直线段P_iP_i+1的长度，k为距离系数，Δt为预设最小离散单元的大小。

对步骤S3所得到的离散角点集合P采用最小二乘法的二次及三次样条曲线拟合形成光顺的曲线。对经光顺处理的曲线进行光顺程度判别，其标准采用能量判别法，能量值越小，对应点集的光顺性越好。并观察曲率图，保证曲率的单调性变化尽可能的小，曲率变化尽可能趋于线性，变化均匀。对得到的所有线进一步采用二次样条插值方法，形成曲线封闭区域的初步结果。

优选的，所述步骤S6中将经程序生成的光滑曲线边界线的数据信息写入特定TXT文件，将其导入三维绘图软件转换为实体引用，以便后续进一步操作。

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

1、本发明所用数值计算方法简单，减少计算时长。

2、本发明有效保证拓扑优化结构各项特性基本稳定，进而保证优化结构原始设计要求在许用范围内。

3、本发明有效消除优化后结构锯齿边界，使之平滑过渡，模型可制造性大大增强，能够直接用于3D打印(增材制造技术)、数控加工等，并且有效降低了拓扑优化后处理工序的难度。

附图说明

图1是本发明一种改进的连续体结构拓扑优化后处理方法的流程图。

图2是本发明的实施例结构示意图。

图3是本发明的实施例拓扑优化结果图。

图4是本发明的实施例二值化处理优化结果图。

图5是本发明的实施例边界提取结果图。

图6是本发明的实施例锯齿边界角点提取图，其中，“*”为角点位置。

图7是本发明的实施例光滑边界的拓扑优化结果。

图8是本发明的实施例实体化模型。

具体实施方式

为了更好的阐述本发明的技术方案，下面结合附图对本发明的具体实施方式进一步详细叙述。

本发明的目的在于克服现有技术的缺点与不足，提出一种改进的连续体结构拓扑优化后处理方法，该方法采用较少的网格划分数目，即有限元分析单元尺寸大于默认值1，对目标模型进行拓扑优化时可有效避免因网格数目过多引起的细小柱状结构；对优化模型进行二值化处理并分析优化模型且提取锯齿形边界线，通过数值计算得到目标离散角点集，将其作为样本点进行曲线拟合，有效避免后处理模型与原始拓扑优化模型产生较大的偏差，影响结构性能；对得到的新边界线进行采取进一步插值处理，进而生成边界光滑的拓扑结构，并对其进行实体化处理。该方法有效保证优化结构原始设计要求，并避免网格依赖性、棋盘格现象。

如附图1所示为本发明的流程图，从图可知，本发明提出一种改进的连续体结构拓扑优化后处理方法，包括以下步骤：

本实施例中以附图2所示二维平面应力结构来进行说明，该结构设计区域为150mm×40mm，对设计区域左右两端采用固定约束，在附图2所示位置受到竖直方向载荷作用。将结构离散成为6000个四结点双线性正四边形单元，材料的弹性模量取E＝10⁸，泊松比取ν＝0.3，体积比设定为0.5。因实施例结构为左右对称结构，为较少计算量，在程序计算时以二分之一结构做分析计算。

S1，对目标模型进行拓扑优化处理得到拓扑优化结构。

连续体结构拓扑优化方法采用如变密度法、渐进结构优化方法(ESO)、独立连续映射方法(ICM)等。本发明专利采用变密度法，以最小体积权值为优化目标，以结构刚度及位移为约束条件，对实施例进行拓扑优化处理，其优化结果如附图3所示。

优化问题可表述为：

subject to：

其中C(ρ)为给定拓扑的柔度，U是全局位移矢量，F是全局负载向量，K是全局刚度矩阵，是单位杨氏模量的元素刚度矩阵，V(ρ)是材料体积，V₀是设计域的体积，f是预先设定的体积分数，ρ_min是一个包含最低允许相对密度的向量。

S2，对得到的拓扑优化结果进行分析并提取边界线，并提取锯齿形边界线的角点得到点集P。

S2-1，对拓扑优化结果进行二值化处理及均值滤波处理，得到边界分明的优化结果，如附图4所示。进一步对优化结果经过腐蚀膨胀边缘处理并进行边界提取，调用graythresh(x)函数自动获取阈值，对提取的边界结果再一次二值化处理，结果如附图5所示，并将二值化结果转化为0-1的数值矩阵。

S2-2，将得到的边界提取结果记为Z(x,y)，对锯齿形边界各角点进行识别计算。

计算Z(x,y)在x、y两方向的梯度I_x、I_y，如下所示：

计算Z(x,y)在x、y两方向梯度的乘积，如下所示：

I_xy＝I_x·I_y

通过使用高斯函数对I_x、I_y及I_xy高斯加权处理，其中σ取1，从而确定出矩阵Q的元素X、Y及XY，如下所示：

XY＝(I_xy)

计算每个节点的响应值R，并对小于阈值t的响应值R置为零，响应值R可定义为：

R＝{R:detQ-α(traceQ)²＜t}

在邻域内对其进行非极大值抑制，其中局部最大值所在点即为所计算角点。

以附图5中边界线左下角位置作为坐标原点，以预先划分的离散单元大小为坐标单元尺寸，进而确定各角点坐标定义为P_j(x_j,y_j)，构成点集P＝{P₁,P₂,···,P_n}。

S3，对步骤S2得到的离散角点进行光顺连接。

S3-1，由下式判断各直线段长度：

其中，为直线段P_iP_i+1的长度，k为距离系数，Δt为预设最小离散单元的大小。本实施例中取k＝4，Δt＝1。

对满足条件的直线段储存于线段集中。

S3-2，步骤2中所得到的离散角点集合P＝{C(x_m,y_m),m＝1,2,···,n}，对其采用最小二乘法的二次及三次样条曲线拟合形成光顺的曲线。设y＝f(x)使得误差平方和最小，即：

其中，φ＝span{φ₁(x),φ₂(x),···,φ_m(x)}

p(x)＝a₁φ₁(x)+a₂φ₂(x)+···+a_mφ_m(x)(m＜n)

对经光顺处理的曲线进行光顺程度判别，其标准采用能量判别法，能量值越小，对应点集的光顺性越好。并观察曲率图，保证曲率的单调性变化尽可能的小，曲率变化尽可能趋于线性，变化均匀。能量值E可由下式计算：

其中，e_i为弦P_i-1P_i方向上的单位向量，k_i为P_i-1P_i的弦长，|e_i+1-e_i|表示弦向量P_iP_i+1与P_i-1P_i之间方向的变化量，k_i+1+k_i为两段弦长度之和。

将上述光顺样条曲线储存于曲线集/>中。

S3-3，将得到的线集L₁、L₂进行交叉合并得到总线集：

并进一步对总线集L采用插值处理，形成曲线封闭区域的初步结果。

S4，对得到的曲线封闭区域进行裁剪处理，得到光顺的曲线封闭区域，进而生成边界光滑的拓扑结构。

S4-1，遍历线集L＝{l₁,l₂,···,l_s}，计算线集L中任意两条线l_i与l_i+1的交点，依次记为Q_i(x_i,y_i)，并将其存入点集Q＝{Q₁,Q₂,···,Q_n}。

S4-2，从交点集Q中任取一点Q_i作为起始点，以线集L中各线作为路径，初始以顺时针方向进行搜索，当搜索至下一交点Q_i+1，进行判断，若为单一路径，则继续沿此路径搜索，若为多路径，分别搜索各路径，仅当搜索至下一交点Q_i+2的路径为有效路径，舍去其余无效路径，重复上述操作继续搜索，直至搜索至初始点Q_i，将循环搜索的有效路径数据记录在数据矩K中。当对所以中点及路径搜索完成后结束操作，得到最终数据集K。

S4-3，将数据集K重新生成曲线，可得到光顺的曲线封闭区域，进而生成边界光滑的拓扑结构，如附图7所示。

S5，检测生成的拓扑结构是否满足要求，满足执行步骤6，否则重设参数重复执行步骤S3至步骤S4。

S6，得到最终的拓扑结构，对其进行实体化处理。

将经程序生成的光滑边界曲线的数据信息写入特定TXT文件，在三维软件(以SolidWorks为例)使用曲线功能导入TXT数据，生成实体曲线；将生成的曲线转换为实体引用，以便后续如倒角、倒圆、增删改的操作，从而可直接进行3D打印或数控加工等。需说明在SolidWorks中，此导入方式生成的曲线不能直接编辑，需新建草图，选择转换实体引用功能，将生成的曲线转化成可编辑的实体，如附图8所示。

上述实施例的描述仅用于进一步详细解释本发明的内容，但并不构成对本发明保护范围的限定。在不脱离本发明的前提下，通过合乎逻辑的分析、推理或对本发明实施例或其中一部分技术特征的修改、等同替换或改进等，均应当视为属于本发明的保护范围之内。

Claims (1)

1.一种改进的连续体结构拓扑优化后处理方法，其特征在于：包括如下步骤：

S1，对目标模型进行拓扑优化处理得到拓扑优化结构；

采用变密度法进行拓扑优化处理，以最小体积权值为优化目标，以结构刚度及位移为约束条件，

优化问题可表述为：

subject to：

其中，C(ρ)为给定拓扑的柔度，U是全局位移矢量，F是全局负载向量，K是全局刚度矩阵，是单位杨氏模量的元素刚度矩阵，V(ρ)是材料体积，V₀是设计域的体积，f是预先设定的体积分数，ρ_min是一个包含最低允许相对密度的向量；

S2，对得到的拓扑优化结果进行分析提取边界线，并提取锯齿形边界线的角点得到点集P；

S2-1，对拓扑优化结果进行二值化处理及均值滤波处理，得到边界分明的优化结果，进一步对优化结果经过腐蚀膨胀边缘处理并进行边界提取，调用graythresh(x)函数自动获取阈值，对提取的边界结果再一次二值化处理，并将二值化结果转化为0-1的数值矩阵；

S2-2，将得到的边界提取结果记为Z(x,y)，对锯齿形边界各角点进行识别计算；

计算Z(x,y)在x、y两方向的梯度I_x、I_y，

计算Z(x,y)在x、y两方向梯度的乘积，I_xy＝I_x·I_y，

通过使用高斯函数对I_x、I_y及I_xy高斯加权处理，确定出矩阵Q的元素X、Y及XY，XY＝(I_xy)；

计算每个节点的响应值R，并对小于阈值t的响应值R置为零，响应值R可定义为：R＝{R:detQ-α(traceQ)²<t}；

在邻域内对其进行非极大值抑制，其中局部最大值所在点即为所计算角点；

以边界线左下角位置作为坐标原点，以预先划分的离散单元大小为坐标单元尺寸，进而确定各角点坐标定义为P_j(x_j,y_j)，构成点集P＝{P₁,P₂,…,P_n}；

S3，对步骤S2得到的离散角点进行光顺连接；

S3-1，判断各直线段长度： 为直线段P_iP_i+1的长度，k为距离系数，Δt为预设最小离散单元的大小；对满足条件的直线段储存于线段集中；

S3-2，点集P＝{C(x_m,y_m),m＝1,2,…,n}，采用最小二乘法的二次及三次样条曲线拟合形成光顺的曲线；

设y＝f(x)使得误差平方和最小，即：

其中，φ＝span{φ₁(x),φ₂(x),···,φ_m(x)}，

p(x)＝a₁φ₁(x)+a₂φ₂(x)+···+a_mφ_m(x)(m＜n)，

对经光顺处理的曲线进行光顺程度判别，并观察曲率图，

能量值E可由下式计算：

其中，e_i为弦P_i-1P_i方向上的单位向量，k_i为P_i-1P_i的弦长，|e_i+1-e_i|表示弦向量P_iP_i+1与P_i-1P_i之间方向的变化量，k_i+1+k_i为两段弦长度之和，

将光顺样条曲线储存于曲线集/>中；

S3-3，将得到的线集L₁、L₂进行交叉合并得到总线集

并进一步对总线集L采用插值处理，形成曲线封闭区域的初步结果；

S4，对得到的曲线封闭区域进行裁剪处理，得到光顺的曲线封闭区域，进而生成边界光滑的拓扑结构；

S4-1，遍历线集L＝{l₁,l₂,…,l_s}，计算线集L中任意两条线l_i与l_i+1的交点，依次记为Q_i(x_i,y_i)，并将其存入点集Q＝{Q₁,Q₂,…,Q_n}；

S4-2，从交点集Q中任取一点Q_i作为起始点，以线集L中各线作为路径，初始以顺时针方向进行搜索，当搜索至下一交点Q_i+1，进行判断，若为单一路径，则继续沿此路径搜索，若为多路径，分别搜索各路径，仅当搜索至下一交点Q_i+2的路径为有效路径，舍去其余无效路径，重复上述操作继续搜索，直至搜索至初始点Q_i，将循环搜索的有效路径数据记录在数据矩K中；当对所以中点及路径搜索完成后结束操作，得到最终数据集K；

S4-3，将数据集K重新生成曲线，可得到光顺的曲线封闭区域，进而生成边界光滑的拓扑结构；

S5，检测生成的拓扑结构是否满足要求，满足执行步骤6，否则重设参数重复执行步骤S3至步骤S4；

S6，得到最终的拓扑结构，对其进行实体化处理。