CN113239584A - 一种优化增材制造方法及系统 - Google Patents

一种优化增材制造方法及系统 Download PDF

Info

Publication number
CN113239584A
CN113239584A CN202110451954.4A CN202110451954A CN113239584A CN 113239584 A CN113239584 A CN 113239584A CN 202110451954 A CN202110451954 A CN 202110451954A CN 113239584 A CN113239584 A CN 113239584A
Authority
CN
China
Prior art keywords
unit
stress data
mechanical model
volume fraction
sensitivity
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202110451954.4A
Other languages
English (en)
Other versions
CN113239584B (zh
Inventor
王宪杰
赵阳
周关根
张帆
王兆毅
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Yunnan University YNU
Original Assignee
Yunnan University YNU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Yunnan University YNU filed Critical Yunnan University YNU
Priority to CN202110451954.4A priority Critical patent/CN113239584B/zh
Publication of CN113239584A publication Critical patent/CN113239584A/zh
Application granted granted Critical
Publication of CN113239584B publication Critical patent/CN113239584B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • G06F30/23Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/20Finite element generation, e.g. wire-frame surface description, tesselation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00Details relating to the application field
    • G06F2113/10Additive manufacturing, e.g. 3D printing
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Graphics (AREA)
  • Software Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明涉及一种优化增材制造方法及系统,该方法包括:构建待优化结构的力学模型,将力学模型通过网格划分为多个单元;将边界条件、荷载和预设的应力数据添加到力学模型中,获得输出力学模型;将各单元区分为实体单元和软单元;对当前体积分数进行缩减;对滤波后的各单元的应力数据进行排序,根据缩减后的体积分数确定应力数据阈值;根据应力数据阈值将输出力学模型中各单元重新区分为实体单元和软单元;判断体积分数是否不等于目标体积分数且未达到收敛条件;若是,则返回步骤“获得当前输出力学模型中各单元的应力数据”;若否,则将输出力学模型中软单元进行删除,获得最终的力学模型。本发明提高了结构刚度。

Description

一种优化增材制造方法及系统
技术领域
本发明涉及结构优化技术领域,特别是涉及一种优化增材制造方法及系统。
背景技术
1992年谢亿民和Steven等人提出了一种基于启发式规则的拓扑优化方法,即渐进结构优化法(evolutionary structural optimization,ESO),这个方法通过逐渐把结构中的低效材料删除,从而实现结构逐渐达到最优状态。其后,为了完善ESO方法只能删除单元的缺陷,又提出双向渐进优化方法(Bi-directional Evolutionary StructureOptimization,BESO),使得在优化过程中能够自由地删除或添加单元。BESO作为目前最常用的连续体拓扑优化方法之一,能够与现成的商用有限元分析软件相连接。
在增材制造(3D打印)技术出现以前,由于结构在拓扑优化后造型复杂,存在传统的加工技术(车、铣、刨、磨)无法加工的地方,因此拓扑优化仅用于结构找型。相对于传统制造的减法工艺而言,增材制造革命般的使用了加法制造技术,通过软件将模型切片,并以层层叠加的方式进行制造。目前,业界使用拓扑优化求解器的软件有Ansys、Nastran、美国澳汰尔公司的Optistruct以及德国Fe-design公司的Tosca等,虽然大部分软件有良好的界面可供使用者操作,但在使用需求上依然不够开放,此类软件在拓扑过程中并没有加入基于制造工艺的约束,因此使用3D成型技术时也有颇多制造难点。
发明内容
本发明的目的是提供一种优化增材制造方法及系统,在结构优化过程中具有很好的收敛性和有效性,提高了结构刚度。
为实现上述目的,本发明提供了如下方案:
一种优化增材制造方法,包括:
提取待优化结构工况中的边界条件和荷载;
构建待优化结构的力学模型,将所述力学模型通过网格划分为多个单元;将所述边界条件和所述荷载添加到所述力学模型中,并将预设的应力数据添加到所述力学模型中,获得输出力学模型;
基于材料差值方法,将各所述单元区分为实体单元和软单元,各所述实体单元构成实体单元集,各所述软单元构成软单元集;获得体积分数,所述体积分数为实体单元的单元数占单元总数的比值;
获得当前输出力学模型中各单元的应力数据;
对各单元的应力数据进行滤波,获得滤波后的各单元的应力数据;
对当前体积分数进行缩减;
对滤波后的各单元的应力数据进行排序,根据缩减后的体积分数确定应力数据阈值;
根据所述应力数据阈值将输出力学模型中各单元重新区分为实体单元和软单元;
判断体积分数是否不等于目标体积分数且未达到收敛条件;
若是,则返回步骤“获得当前输出力学模型中各单元的应力数据”;
若否,则将输出力学模型中软单元进行删除,获得最终的力学模型。
可选地,所述对各单元的应力数据进行滤波,获得滤波后的各单元的应力数据,具体包括:
计算各单元的灵敏度;
获得各单元的应力数据中最大应力值;
依次计算以单元i为中心的设定半径范围内所有单元的灵敏度加权均值,用该灵敏度加权均值更新单元i的灵敏度,实现对所有单元进行灵敏度的修正,其中i表示第i个单元;
将单元i的灵敏度与所述最大应力值的平均值作为单元i的应力数据,获得滤波后的各单元的应力数据;
灵敏度的计算公式为:
Figure BDA0003039047410000021
其中,a表示灵敏度,E1表示实体单元的弹性模量,E2表示软单元的弹性模量,U表示位移,K表示刚度值,T表示转置操作,Xmin表示软单元乘子,P表示惩罚指数。
可选地,所述收敛条件包括收敛公式,所述收敛公式表示为:
Figure BDA0003039047410000031
其中,e为相对误差,K为迭代次数,N为正常数,i′∈[1,N],CK-i′+1表示第K次迭代时各单元的应力数据中最大应力值,当e小于1%时,达到收敛条件。
可选地,所述对当前体积分数进行缩减,具体包括:
根据公式V0=V0′*(1-er)对当前体积分数进行缩减,其中,V0表示缩减后的体积分数,V0′*(1-er)表示当前体积分数,er表示体积缩减率。
可选地,所述力学模型为通过ABAQUS有限元软件构建的有限元模型。
本发明还公开了一种优化增材制造系统,包括:
约束条件提取模块,用于提取待优化结构工况中的边界条件和荷载;
输出力学模型构建模块,用于构建待优化结构的力学模型,将所述力学模型通过网格划分为多个单元;将所述边界条件和所述荷载添加到所述力学模型中,并将预设的应力数据添加到所述力学模型中,获得输出力学模型;
单元区分模块,用于基于材料差值方法,将各所述单元区分为实体单元和软单元,各所述实体单元构成实体单元集,各所述软单元构成软单元集;获得体积分数,所述体积分数为实体单元的单元数占单元总数的比值;
应力数据获得模块,用于获得当前输出力学模型中各单元的应力数据;
应力数据滤波模块,用于对各单元的应力数据进行滤波,获得滤波后的各单元的应力数据;
体积分数缩减模块,用于对当前体积分数进行缩减;
应力数据阈值获得模块,用于对滤波后的各单元的应力数据进行排序,根据缩减后的体积分数确定应力数据阈值;
单元重新区分模块,用于根据所述应力数据阈值将输出力学模型中各单元重新区分为实体单元和软单元;
迭代判断模块,用于判断体积分数是否不等于目标体积分数且未达到收敛条件;
返回模块,用于当迭代判断模块判断为是时,返回应力数据获得模块;
最终的力学模型获得模块,用于当迭代判断模块判断为否时,将输出力学模型中软单元进行删除,获得最终的力学模型。
可选地,所述应力数据滤波模块,具体包括:
灵敏度计算单元,用于计算各单元的灵敏度,其中,灵敏度的计算公式为:
Figure BDA0003039047410000041
其中,a表示灵敏度,E1表示实体单元的弹性模量,E2表示软单元的弹性模量,U表示位移,K表示刚度值,T表示转置操作,Xmin表示软单元乘子,P表示惩罚指数;
最大应力值获得单元,用于获得各单元的应力数据中最大应力值;
灵敏度修正单元,用于依次计算以单元i为中心的设定半径范围内所有单元的灵敏度加权均值,用该灵敏度加权均值更新单元i的灵敏度,实现对所有单元进行灵敏度的修正,其中i表示第i个单元;
应力数据修正单元,用于将单元i的灵敏度与所述最大应力值的平均值作为单元i的应力数据,获得滤波后的各单元的应力数据。;
可选地,所述收敛条件包括收敛公式,所述收敛公式表示为:
Figure BDA0003039047410000042
其中,e为相对误差,K为迭代次数,N为正常数,i′∈[1,N],CK-i′+1表示第K次迭代时各单元的应力数据中最大应力值,当e小于1%时,达到收敛条件。
可选地,所述体积分数缩减模块,具体包括:
体积分数缩减单元,用于根据公式V0=V0′*(1-er)对当前体积分数进行缩减,其中,V0表示缩减后的体积分数,V0′*(1-er)表示当前体积分数,er表示体积缩减率。
可选地,所述力学模型为通过ABAQUS有限元软件构建的有限元模型。
根据本发明提供的具体实施例,本发明公开了以下技术效果:
本发明根据设定的目标体积分数,对初始体积分数进行缩减,每次缩减对输出力学模块中单元中的实体单元和软单元进行重新划分,直到达到收敛条件或目标体积分数,在满足体积约束的情况下提升结构整体的刚度。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案,下面将对实施例中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图仅仅是本发明的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动性的前提下,还可以根据这些附图获得其他的附图。
图1为本发明一种优化增材制造方法流程示意图;
图2为本发明一种优化增材制造系统结构示意图;
图3为结构模型节点受荷分析示意图;
图4为本发明rhino建模示意图;
图5为本发明ABAQUS网格划分示意图;
图6为本发明增材制造优化结果示意图;
图7为ABAQUS的Tosca模块优化结果示意图;
图8为本发明优化结果有限元提取示意图;
图9为本发明设计域边界光滑处理示意图;
图10为本发明打印软件操作界面示意图;
图11为本发明获得的最终力学模型的应力性能分析图一;
图12为本发明获得的最终力学模型的应力性能分析图二;
图13为本发明获得的最终力学模型的应力性能分析图三。
具体实施方式
下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例仅仅是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
本发明的目的是提供一种优化增材制造方法及系统,在结构优化过程中具有很好的收敛性和有效性,提高了结构刚度。
为使本发明的上述目的、特征和优点能够更加明显易懂,下面结合附图和具体实施方式对本发明作进一步详细的说明。
图1为本发明一种优化增材制造方法流程示意图,如图1所示,一种优化增材制造方法,包括:
步骤101:提取待优化结构工况中的边界条件和荷载。
整体工况采集(以网壳结构为例):通过ABAQUS有限元软件直接建立整体网壳结构的有限元模型,基于实际工况(按照规范设置荷载)在ABAQUS中进行有限元分析,计算后采集与需要拓扑优化的结构部件相连位置处的边界条件和荷载工况(轴力、剪力、弯矩等),这一步在ABAQUS内部进行单元节点值查询即可,提取结果值应用于步骤102中的需要拓扑优化的部件的荷载添加。
步骤102:构建待优化结构的力学模型,将所述力学模型通过网格划分为多个单元;将所述边界条件和所述荷载添加到所述力学模型中,并将预设的应力数据添加到所述力学模型中,获得输出力学模型。所述力学模型为通过ABAQUS有限元软件构建的有限元模型。
步骤102建立需要优化部件的力学模型具体为:在三维软件中(这里使用rhino)进行优化部件模型的建模,将得到的模型导入ABAQUS中进行网格划分,(网格划分后单元序号自动生成),获得各单元,基于步骤101中采集的荷载以及边界条件添加到ABAQUS力学模型中,并添加好需要的场输出,场输出包括应力、应变和能量等,另外基于材料差值方法预先划分好set1(集合1)和set2(集合2),set1为实体单元集,set2为软单元集,将远离受荷载区域的几个单元划入set2,其余划入set1,然后导出模型的.inp文件,inp文件将用于应力值计算有限元分析。
步骤103:基于材料差值方法,将各所述单元区分为实体单元和软单元,各所述实体单元构成实体单元集,各所述软单元构成软单元集;获得体积分数,所述体积分数为实体单元的单元数占单元总数的比值。
步骤104:获得当前输出力学模型中各单元的应力数据。通过编程软件(MATLAB、Python等)调用ABAQUS分析模块对.inp有限元文件进行后台计算,计算后得到存有计算结果值的.ODB文件,编程提取.ODB文件中优化模型每个单元格的计算结果值信息(应力、应变、位移和应变能等)。
步骤105:对各单元的应力数据进行滤波,获得滤波后的各单元的应力数据。通过编程将优化部分所有单元的计算结果值进行灵敏度敏度分析,对分析结果进行滤波,通过加入一个灵敏度算子rmin进行应力均匀处理,均匀处理后的应力值将应用于材料划分,将使得模型结果更为光滑有效。
所述对各单元的应力数据进行滤波,获得滤波后的各单元的应力数据,具体包括:
计算各单元的灵敏度;
获得各单元的应力数据中最大应力值;
依次计算以单元i为中心的设定半径范围内所有单元的灵敏度加权均值,用该灵敏度加权均值更新单元i的灵敏度,实现对所有单元进行灵敏度的修正,其中i表示第i个单元;
将单元i的灵敏度与所述最大应力值的平均值作为单元i的应力数据,获得滤波后的各单元的应力数据;
灵敏度的计算公式为:
Figure BDA0003039047410000071
其中,a表示灵敏度,E1表示实体单元的弹性模量,E2表示软单元的弹性模量,U表示位移,K表示刚度值,T表示转置操作,Xmin表示软单元乘子,P表示惩罚指数。
步骤106:对当前体积分数进行缩减。
所述对当前体积分数进行缩减,具体包括:
根据公式V0=V0′*(1-er)对当前体积分数进行缩减,其中,V0表示缩减后的体积分数,V0′*(1-er)表示当前体积分数,er表示体积缩减率。
步骤107:对滤波后的各单元的应力数据进行排序,根据缩减后的体积分数确定应力数据阈值。
步骤108:根据所述应力数据阈值将输出力学模型中各单元重新区分为实体单元和软单元。
步骤109:判断体积分数是否不等于目标体积分数且未达到收敛条件。
所述收敛条件包括收敛公式,所述收敛公式表示为:
Figure BDA0003039047410000081
其中,e为相对误差,K为迭代次数,N为正常数,i′∈[1,N],CK-i′+1表示第K次迭代时各单元的应力数据中最大应力值,当e小于1%时,达到收敛条件。
若是,则返回步骤104。
若否,则执行步骤110。
步骤110:将输出力学模型中软单元进行删除,获得最终的力学模型。具体为:输出力学模型的单元划分只是在.inp文件内部进行了形式上的单元划分,而实际需要删除的单元信息还存在,因此通过编程在.inp文件中将需要删除的单元及其信息进行实际性的删除,然后进入ABAQUS将三维OBJ模型提取出来,而后进行光滑处理及打印。
光滑处理及打印:将得到的OBJ模型放入三维建模软件中进行光滑处理,进而导出STL文件,放入三维打印软件中进行打印,即可完成增材制造。
本发明使得算法在结构优化过程中具有很好的收敛性和有效性,并通过衔接ABAQUS求解器以及MATLAB进行二次开发,其简单灵活的方式为后续的算法改进提供了广阔的平台。
下面以具体实施例说明本发明一种优化增材制造方法。
Step1:整体工况采集
对需要进行拓扑优化的结构进行工况分析,这里所用的案例为云南大学图书馆穹顶,如图3所示,提取工况中的边界条件和受荷方式。将模型在ABAQUS中进行建模,基于实际工况进行有限元分析,然后提取设计域附近连接杆件的荷载施加情况及其边界条件,这里提取优化节点附近杆件的轴力情况,所提取据将应用于Step2中节点处(即设计域)的边界条件和荷载添加。
Step2:建立需要优化部件的力学模型
在Rhino中进行建模,建模采用40*40*40的正六面体,如图4所示,由于BESO算法中单元能够自由删除或添加,这样建模能够提供模型较大的优化空间。在ABAQUS的设计域part中设置两个以材料属性为划分原则的set集,分别为软单元集(set2)和实体单元集(set1),用来存放后面划分好的单元。图5为将建好的模型导入Abaqus进行网格划分后的示意图。将Step1中所得的荷载情况和边界条件进行添加。设计域采用了计算性能较好的六面体网格。在研究材料的最优分布时,常用以结构最小柔度为设计目标,柔度即应变能,应变能越小则变形越小,刚度越大。在分析中将场输出设置为应变能密度ESEDEN。最后输出设置好的有限元.inp文件。
表1优化模型属性参数表
Figure BDA0003039047410000091
Step3:设置目标函数
设置初始体积比V0,初始体积比V0根据Step2中set1中单元数占总单元数的比值设置;目标体积比(目标体积分数)V=50%,目标体积按照需求自行设置即可;进化率er=0.01,此参数为体积缩减率。
当前体积收缩按照V0=V0*(1-er)进行缩减。
取当前收缩后的体积比V0和目标体积比V两者的最大值作为当前收缩体积比。
Step4:ABAQUS的MATLAB调用
在MATLAB中通过dos命令调用ABAQUS对Step2中的.inp文件进行分析,得到分析后的ODB数据文件,基于Python语言中有相应的ODB数据库,通过MATLAB调用Python程序,对ODB数据文件进行应力值的提取,提取后数据存入TXT文本中,这里的提取的是应变能密度ESEDEN,接下来进入Step5分析。
Step5:单元应力滤波
编程读取Step4中TXT文件中的应力数据后,应力数据与ABAQUS中模型的单元序号排列位置对应放入矩阵中。
为了使得到的模型更加的光滑,接下来进行灵敏度分析。算法中,单相材料优化模型中软单元将被视作第二材料,分别对实体单元添加弹性模量E1,软单元添加弹性模量E2,并加入软单元乘子Xmin(一个趋向于0的极小值)和惩罚指数P,使得模型在有限元计算中不引起奇异现象,增加模型的光滑程度。
其灵敏度模型定义如下:
Figure BDA0003039047410000101
其中,U表示位移,K表示刚度值,T表示转置操作,UTKU表示单元应变能。
随着惩罚指数P的不断增加,上述灵敏度模型可等价为:
Figure BDA0003039047410000102
上述公式子中,a即为根据Step4中初始应力状态灵敏度进行灵敏度分析后所得的单元灵敏度,使应力大的有效单元的应力值将变得更大,应力小的有效单元的应力值将变得更小,进而使得所有单元更加清晰,在一定程度上有效解决棋盘格现象。随后进行有效单元的灵敏度修正。
在本实施例中,为了识别参与计算修正节点灵敏度的单元是否为有效单元,引入一个半径为rmin的圆,将圆内所有单元(L个)包含进去进行单元ijk的灵敏度计算。
通过for循环对包含所有经过公式(2)惩罚后的单元的应力值的三维矩阵ijk中进行半径为rmin范围内单元应力值的敏度计算,将圆内所有单元(L个)包含进去进行单元的灵敏度计算,将ijk位置处的单元的灵敏度数修改,i表示x轴的坐标,j表示y的坐标,k表示z轴的坐标。如下式:
Figure BDA0003039047410000111
式中,权重因子ω(rn)的计算公式为ω(rn)=rmin-rn,rn为矩阵内第n个单元与ijk位置处矩阵的距离,an表示第n个单元的灵敏度,aijk单元ijk修改后的灵敏度。
最终获得单元应力滤波后的单元应力值矩阵stress_all,单元应力值矩阵stress_all中个元素为各单元修改后的灵敏度。
记录下滤波前的单元应力矩阵最大应力值ss_max。
Figure BDA0003039047410000112
将上一步循环得到的平均后的应力数值stress_all记为old_stress,与当前循环所得的应力值相加后平均确定修改后的应力值,降低应力突变。Sj′表示当次循环单元的应力值,Sj′-1表示上一步循环单元的应力值。
Step6:去留单元划分准则
单元的删添可通过判断删添单元数量是否满足目标体积V来进行控制,取平均后的单元灵敏度的最大值αmax和最小值αmin的一半作为添加(实体单元)和删除(软单元)的灵敏度阈值th,当体积多于或少于目标体积分数时,通过在最大最小值范围内做差值逼近,调整阈值(灵敏度阈值)上下浮动,直到找到满足目标体积的阈值。
Figure BDA0003039047410000121
这一步将得到需要保留或删除的单元的单元序号,αi表示第i个单元的灵敏度值,VAadd表示保留(实体单元)的体积分数,VAdel表示删除(软单元)的体积分数,V*表示当前需要缩减到的体积分数,VAadd>V*时则阈值增大,VAdel>V*时阈值减小。
将划分好的单元的单元序号重新写入.inp文件中进行下一次循环计算,写入时分别写入.inp文件中的两个单元集set1和set2中,这两个单元集根据材料差值进行单元集属性划分,划分为近似不存在的软单元和高效受力的实体单元。随后进行下一步迭代。
Step7:收敛判断
未达到Step3中目标体积分数时,将Step6划分单元后的.inp文件重新代入Step4,同时重复Step4-Step7。若达到Step3中目标体积分数时,则对比上一循环应力值和当前循环应力值的误差,根据目标函数ss_max的历史变化量来判断优化是否收敛,收敛准则定义为:
Figure BDA0003039047410000122
其中e为相对误差,K为迭代次数,N取正常数。
当误差值小于1%时(接近于0),则收敛完成,即得到正确结果,否则重复Step4-Step7。正确结果将在Step8中提取出来。
Step8:几何模型的.OBJ文件提取
计算完成后,基于算法最后一步迭代的Step6划分好的单元序号进行矩阵交集计算,将需要保留的单元提取出来,放入.inp文件中,在ABAQUS中将模型的OBJ几何模型提取出来。
这里的几何模型提取分为两步骤:1、修改单元,通过将两种材料的单元集做交集,最终将不需要的单元序号反选掉,剩余的写入文件中;2、修改节点,节点中需要注意的是两种材料的共用节点不能删除。模型修改方法定义如下:
Figure BDA0003039047410000131
式中[M]为第一种材料的单元序号与所有单元序号交集的矩阵,[M2]第一种材料的单元序号矩阵,[Me]为所有单元序号的矩阵,[N]为第一种材料单元的节点序号与第二种材料单元的节点序号交集的矩阵,[N’]为与第一种材料无关联的所有节点序号矩阵,[N1]为第一种材料(set1)实体单元单元的节点序号,[N2]为第二种材料(set2)软单元的节点序号。
对于拓扑优化来说,删去低效单元后,最终的结果应该是体积减少,同时总体应变能减少,这样刚度才会变大。如图6、图7所示,此方法得出结果与ABAQUS的Tosca模块结果基本一致。随着迭代进行,最大应变能减小,最小应变能增加,结构的单元应变能更加趋向于分布均匀,从应力云图为模型优化后的结果进行分析,分析结果如图11-图13所示,可以看到应力分布均匀,整体属于高效受力状态,该优化算法属于良性优化,结果是可靠的。
最终得到提取出来的ABAQUS有限元模型的inp文件。
Step9:优化模型光滑及打印
图8为模型的优化结果提取示意图,将Step8提取出来的inp文件模型导出OBJ格式,在rhino中读取OBJ文件,通过调整平滑系数和平滑次数进行边界处理,控制光滑后模型的光滑度及其体积,可以看到处理后模型的锯齿形状有所钝化,如图9所示,得到光滑后的三维模型的.stl文件格式。
然后在打印软件中进行切片打印。图10为设计域的优化结果在打印软件中的操作示意图。将其导入FlashPrint中,调整模型的摆放位置,进行打印工艺参数设定。可以设置的工艺参数为:温度、打印速度、打印高度、填充率等。
本方法通过使用三维建模软件、有限元分析软件以及编程工具等众多平台软件,较好地完成了模型设计-制造的全过程,本发明基于材料插值方法,通过引入惩罚指数P,对传统的BESO方法的灵敏度算法进行了合理的改进,结构优化过程中具有很好的收敛性和有效性,节省了分析时间,极大地提高了效率,另外本发明能够较为开放的对优化的需求进行人为干预,方便后续各种工艺约束的加入,同时程序的优化结果与ABAQUS自带的Tosca模块优化的结果基本相似,由应力云图以及性能指标分析可知,目前算法基本能够达到所需的优化目标,在满足体积约束的情况下提升了结构整体刚度,程序本身具有不错的鲁棒性。
图2为本发明一种优化增材制造系统结构示意图,如图2所示,本发明还公开了一种优化增材制造系统,包括:
约束条件提取模块201,用于提取待优化结构工况中的边界条件和荷载;
输出力学模型构建模块202,用于构建待优化结构的力学模型,将所述力学模型通过网格划分为多个单元;将所述边界条件和所述荷载添加到所述力学模型中,并将预设的应力数据添加到所述力学模型中,获得输出力学模型;
单元区分模块203,用于基于材料差值方法,将各所述单元区分为实体单元和软单元,各所述实体单元构成实体单元集,各所述软单元构成软单元集;获得体积分数,所述体积分数为实体单元的单元数占单元总数的比值;
应力数据获得模块204,用于获得当前输出力学模型中各单元的应力数据;
应力数据滤波模块205,用于对各单元的应力数据进行滤波,获得滤波后的各单元的应力数据;
体积分数缩减模块206,用于对当前体积分数进行缩减;
应力数据阈值获得模块207,用于对滤波后的各单元的应力数据进行排序,根据缩减后的体积分数确定应力数据阈值;
单元重新区分模块208,用于根据所述应力数据阈值将输出力学模型中各单元重新区分为实体单元和软单元;
迭代判断模块209,用于判断体积分数是否不等于目标体积分数且未达到收敛条件;
返回模块210,用于当迭代判断模块209判断为是时,返回应力数据获得模块204;
最终的力学模型获得模块211,用于当迭代判断模块209判断为否时,将输出力学模型中软单元进行删除,获得最终的力学模型。
所述应力数据滤波模块205,具体包括:
灵敏度计算单元,用于计算各单元的灵敏度,其中,灵敏度的计算公式为:
Figure BDA0003039047410000151
其中,a表示灵敏度,E1表示实体单元的弹性模量,E2表示软单元的弹性模量,U表示位移,K表示刚度值,T表示转置操作,Xmin表示软单元乘子,P表示惩罚指数。
最大应力值获得单元,用于获得各单元的应力数据中最大应力值。
灵敏度修正单元,用于依次计算以单元i为中心的设定半径范围内所有单元的灵敏度加权均值,用该灵敏度加权均值更新单元i的灵敏度,实现对所有单元进行灵敏度的修正,其中i表示第i个单元。
应力数据修正单元,用于将单元i的灵敏度与所述最大应力值的平均值作为单元i的应力数据,获得滤波后的各单元的应力数据。
所述收敛条件包括收敛公式,所述收敛公式表示为:
Figure BDA0003039047410000152
其中,e为相对误差,K为迭代次数,N为正常数,i′∈[1,N],CK-i′+1表示第K次迭代时各单元的应力数据中最大应力值,当e小于1%时,达到收敛条件。
所述体积分数缩减模块206,具体包括:
体积分数缩减单元,用于根据公式V0=V0′*(1-er)对当前体积分数进行缩减,其中,V0表示缩减后的体积分数,V0′*(1-er)表示当前体积分数,er表示体积缩减率。
所述力学模型为通过ABAQUS有限元软件构建的有限元模型。
本说明书中各个实施例采用递进的方式描述,每个实施例重点说明的都是与其他实施例的不同之处,各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言,由于其与实施例公开的方法相对应,所以描述的比较简单,相关之处参见方法部分说明即可。
本文中应用了具体个例对本发明的原理及实施方式进行了阐述,以上实施例的说明只是用于帮助理解本发明的方法及其核心思想;同时,对于本领域的一般技术人员,依据本发明的思想,在具体实施方式及应用范围上均会有改变之处。综上所述,本说明书内容不应理解为对本发明的限制。

Claims (10)

1.一种优化增材制造方法,其特征在于,包括:
提取待优化结构工况中的边界条件和荷载;
构建待优化结构的力学模型,将所述力学模型通过网格划分为多个单元;将所述边界条件和所述荷载添加到所述力学模型中,并将预设的应力数据添加到所述力学模型中,获得输出力学模型;
基于材料差值方法,将各所述单元区分为实体单元和软单元,各所述实体单元构成实体单元集,各所述软单元构成软单元集;获得体积分数,所述体积分数为实体单元的单元数占单元总数的比值;
获得当前输出力学模型中各单元的应力数据;
对各单元的应力数据进行滤波,获得滤波后的各单元的应力数据;
对当前体积分数进行缩减;
对滤波后的各单元的应力数据进行排序,根据缩减后的体积分数确定应力数据阈值;
根据所述应力数据阈值将输出力学模型中各单元重新区分为实体单元和软单元;
判断体积分数是否不等于目标体积分数且未达到收敛条件;
若是,则返回步骤“获得当前输出力学模型中各单元的应力数据”;
若否,则将输出力学模型中软单元进行删除,获得最终的力学模型。
2.根据权利要求1所述的优化增材制造方法,其特征在于,所述对各单元的应力数据进行滤波,获得滤波后的各单元的应力数据,具体包括:
计算各单元的灵敏度;
获得各单元的应力数据中最大应力值;
依次计算以单元i为中心的设定半径范围内所有单元的灵敏度加权均值,用该灵敏度加权均值更新单元i的灵敏度,实现对所有单元进行灵敏度的修正,其中i表示第i个单元;
将单元i的灵敏度与所述最大应力值的平均值作为单元i的应力数据,获得滤波后的各单元的应力数据;
灵敏度的计算公式为:
Figure FDA0003039047400000021
其中,a表示灵敏度,E1表示实体单元的弹性模量,E2表示软单元的弹性模量,U表示位移,K表示刚度值,T表示转置操作,Xmin表示软单元乘子,P表示惩罚指数。
3.根据权利要求1所述的优化增材制造方法,其特征在于,所述收敛条件包括收敛公式,所述收敛公式表示为:
Figure FDA0003039047400000022
其中,e为相对误差,K为迭代次数,N为正常数,i′∈[1,N],CK-i′+1表示第K次迭代时各单元的应力数据中最大应力值,当e小于1%时,达到收敛条件。
4.根据权利要求1所述的优化增材制造方法,其特征在于,所述对当前体积分数进行缩减,具体包括:
根据公式V0=V0′*(1-er)对当前体积分数进行缩减,其中,V0表示缩减后的体积分数,V0′*(1-er)表示当前体积分数,er表示体积缩减率。
5.根据权利要求1所述的优化增材制造方法,其特征在于,所述力学模型为通过ABAQUS有限元软件构建的有限元模型。
6.一种优化增材制造系统,其特征在于,包括:
约束条件提取模块,用于提取待优化结构工况中的边界条件和荷载;
输出力学模型构建模块,用于构建待优化结构的力学模型,将所述力学模型通过网格划分为多个单元;将所述边界条件和所述荷载添加到所述力学模型中,并将预设的应力数据添加到所述力学模型中,获得输出力学模型;
单元区分模块,用于基于材料差值方法,将各所述单元区分为实体单元和软单元,各所述实体单元构成实体单元集,各所述软单元构成软单元集;获得体积分数,所述体积分数为实体单元的单元数占单元总数的比值;
应力数据获得模块,用于获得当前输出力学模型中各单元的应力数据;
应力数据滤波模块,用于对各单元的应力数据进行滤波,获得滤波后的各单元的应力数据;
体积分数缩减模块,用于对当前体积分数进行缩减;
应力数据阈值获得模块,用于对滤波后的各单元的应力数据进行排序,根据缩减后的体积分数确定应力数据阈值;
单元重新区分模块,用于根据所述应力数据阈值将输出力学模型中各单元重新区分为实体单元和软单元;
迭代判断模块,用于判断体积分数是否不等于目标体积分数且未达到收敛条件;
返回模块,用于当迭代判断模块判断为是时,返回应力数据获得模块;
最终的力学模型获得模块,用于当迭代判断模块判断为否时,将输出力学模型中软单元进行删除,获得最终的力学模型。
7.根据权利要求6所述的优化增材制造系统,其特征在于,所述应力数据滤波模块,具体包括:
灵敏度计算单元,用于计算各单元的灵敏度,其中,灵敏度的计算公式为:
Figure FDA0003039047400000031
其中,a表示灵敏度,E1表示实体单元的弹性模量,E2表示软单元的弹性模量,U表示位移,K表示刚度值,T表示转置操作,Xmin表示软单元乘子,P表示惩罚指数;
最大应力值获得单元,用于获得各单元的应力数据中最大应力值;
灵敏度修正单元,用于依次计算以单元i为中心的设定半径范围内所有单元的灵敏度加权均值,用该灵敏度加权均值更新单元i的灵敏度,实现对所有单元进行灵敏度的修正,其中i表示第i个单元;
应力数据修正单元,用于将单元i的灵敏度与所述最大应力值的平均值作为单元i的应力数据,获得滤波后的各单元的应力数据。
8.根据权利要求6所述的优化增材制造系统,其特征在于,所述收敛条件包括收敛公式,所述收敛公式表示为:
Figure FDA0003039047400000041
其中,e为相对误差,K为迭代次数,N为正常数,i′∈[1,N],CK-i′+1表示第K次迭代时各单元的应力数据中最大应力值,当e小于1%时,达到收敛条件。
9.根据权利要求6所述的优化增材制造系统,其特征在于,所述体积分数缩减模块,具体包括:
体积分数缩减单元,用于根据公式V0=V0′*(1-er)对当前体积分数进行缩减,其中,V0表示缩减后的体积分数,V0′*(1-er)表示当前体积分数,er表示体积缩减率。
10.根据权利要求6所述的优化增材制造系统,其特征在于,所述力学模型为通过ABAQUS有限元软件构建的有限元模型。
CN202110451954.4A 2021-04-26 2021-04-26 一种优化增材制造方法及系统 Active CN113239584B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110451954.4A CN113239584B (zh) 2021-04-26 2021-04-26 一种优化增材制造方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110451954.4A CN113239584B (zh) 2021-04-26 2021-04-26 一种优化增材制造方法及系统

Publications (2)

Publication Number Publication Date
CN113239584A true CN113239584A (zh) 2021-08-10
CN113239584B CN113239584B (zh) 2022-03-11

Family

ID=77129204

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110451954.4A Active CN113239584B (zh) 2021-04-26 2021-04-26 一种优化增材制造方法及系统

Country Status (1)

Country Link
CN (1) CN113239584B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113704888A (zh) * 2021-08-23 2021-11-26 中国飞机强度研究所 一种单元应力筛选方法
CN115630412A (zh) * 2022-10-25 2023-01-20 浙江大学 基于多轴3d打印的自支撑结构优化设计及制造方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110069800A (zh) * 2018-11-17 2019-07-30 华中科技大学 具有光滑边界表达的三维结构拓扑优化设计方法及设备
CN110414165A (zh) * 2019-08-01 2019-11-05 华东交通大学 一种基于全局应力约束的多相材料柔顺机构拓扑优化方法
CN111737839A (zh) * 2020-05-19 2020-10-02 广州大学 基于动态进化率和自适应网格的beso拓扑优化方法及其应用
WO2020215533A1 (zh) * 2019-04-26 2020-10-29 大连理工大学 一种基于材料场缩减级数展开的结构拓扑优化方法
CN112100882A (zh) * 2020-08-27 2020-12-18 华南理工大学 一种具有光滑边界的连续体结构密度演化拓扑优化方法
CN112100774A (zh) * 2020-09-16 2020-12-18 哈尔滨理工大学 一种基于变密度法的应力和应变能双约束的拓扑优化方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110069800A (zh) * 2018-11-17 2019-07-30 华中科技大学 具有光滑边界表达的三维结构拓扑优化设计方法及设备
WO2020215533A1 (zh) * 2019-04-26 2020-10-29 大连理工大学 一种基于材料场缩减级数展开的结构拓扑优化方法
CN110414165A (zh) * 2019-08-01 2019-11-05 华东交通大学 一种基于全局应力约束的多相材料柔顺机构拓扑优化方法
CN111737839A (zh) * 2020-05-19 2020-10-02 广州大学 基于动态进化率和自适应网格的beso拓扑优化方法及其应用
CN112100882A (zh) * 2020-08-27 2020-12-18 华南理工大学 一种具有光滑边界的连续体结构密度演化拓扑优化方法
CN112100774A (zh) * 2020-09-16 2020-12-18 哈尔滨理工大学 一种基于变密度法的应力和应变能双约束的拓扑优化方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
X. HUANG ET AL.: "Bi-directional evolutionary topology optimization of continuum structures with one or multiple materials", 《COMPUT MECH》 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113704888A (zh) * 2021-08-23 2021-11-26 中国飞机强度研究所 一种单元应力筛选方法
CN113704888B (zh) * 2021-08-23 2024-02-23 中国飞机强度研究所 一种单元应力筛选方法
CN115630412A (zh) * 2022-10-25 2023-01-20 浙江大学 基于多轴3d打印的自支撑结构优化设计及制造方法
CN115630412B (zh) * 2022-10-25 2023-04-28 浙江大学 基于多轴3d打印的自支撑结构优化设计及制造方法

Also Published As

Publication number Publication date
CN113239584B (zh) 2022-03-11

Similar Documents

Publication Publication Date Title
CN110110413B (zh) 一种基于材料场缩减级数展开的结构拓扑优化方法
CN113239584B (zh) 一种优化增材制造方法及系统
CN111737839B (zh) 基于动态进化率和自适应网格的beso拓扑优化方法及其应用
CN111027110B (zh) 一种连续体结构拓扑与形状尺寸综合优化方法
CN110069800B (zh) 具有光滑边界表达的三维结构拓扑优化设计方法及设备
CN108647370B (zh) 基于双环迭代的无人直升机气动外形优化设计方法
CN111859733B (zh) 一种基于蚁群算法的汽车排气系统可靠性优化方法
CN112836411A (zh) 加筋板壳结构的优化方法、装置、计算机设备和存储介质
CN108763658A (zh) 基于等几何方法的组合薄壁结构固有频率设计方法
CN110414127B (zh) 一种面向增材制造的支撑体积约束拓扑优化方法
CN112257897B (zh) 基于改进多目标粒子群的电动汽车充电优化方法和系统
CN113536623B (zh) 一种材料不确定性结构稳健性拓扑优化设计方法
CN112287447B (zh) 钢结构框架结构智能优化系统及方法
CN113591191A (zh) 基于整体下滑方向的二维边坡稳定矢量和分析方法及系统
CN113051796B (zh) 一种应用于增材制造的结构拓扑优化设计方法
US11217013B2 (en) Method for preserving shapes in solid model when distributing material during topological optimization
CN111539138A (zh) 基于阶跃函数的结构动力学峰值时域响应灵敏度求解方法
CN116562166A (zh) 一种基于ihba的配电网分布式电源选址定容方法
CN111274624A (zh) 一种基于rbf代理模型的多工况异形节点拓扑优化设计方法
CN114239363A (zh) 一种基于ABAQUS二次开发Python语言的变密度拓扑优化的方法
CN113239457A (zh) 一种基于灰聚类算法模型的多工况车架拓扑优化方法
CN112766609A (zh) 一种基于云计算的用电量预测方法
CN113515824B (zh) 一种筋条布局与基板形状协同的拓扑优化设计方法
CN115660488A (zh) 基于混沌粒子群遗传算法的区域水资源承载力评价方法
CN117150851A (zh) 一种基于智能交互的结构设计分析优化一体化方法及系统

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant