一种基于Abaqus二次开发散热器主片的仿真方法与流程

一种基于abaqus二次开发散热器主片的仿真方法
技术领域
1.本发明涉及汽车发动机水箱散热器主片成型制造技术领域，具体涉及一种基于abaqus二次开发散热器主片的仿真方法。


背景技术：

2.汽车散热器是汽车的重要的零部件之一，而散热器主片是连接散热器芯子和水室的零件，其成型质量直接决定了散热器的质量。主片如果过度减薄则可能受力开裂导致散热器泄露，同时主片如果过度起皱则可能导致装配应力。为缩短开发周期，降低成本，近年来冲压成型分析软件autoform、pamstamp等逐渐兴盛，能够实现冲压件从生产产到应用全过程的有限元仿真分析。冲压成型仿真是模拟冲压件从落料到成型的一系列过程，可以对零件减薄率、起皱风险、成型性能等进行分析，从而指导生产。
3.公开号为cn104881531b，名称为一种autoform冲压成型信息到碰撞仿真模型的映射方法，其公开了采用autoform进行冲压成型过程仿真，但是汽车散热器用主片采用autoform软件分析冲压成型，由于其设置简单、计算速度快和自适应网格技术等，然而，autoform内置的barlat89屈服准则并不能很好的预测铝合金起皱现象，导致了仿真结果无法指导实际生产，同时由于autoform软件封闭性强，无法对内置准则进行二次开发。
4.综上所述，如何实现复杂本构模型和仿真求解器的结合，从而实现冲压成型过程中对成型质量的高精度预测，是目前本领域需要解决的重要技术问题。


技术实现要素：

5.本发明的目的是针对现有技术的不足，提供一种基于abaqus二次开发散热器主片的仿真方法，解决现有仿真软件不能很好的预测铝合金起皱现象，仿真结果无法指导实际生产，封闭性强，无法对内置准则进行二次开发的问题。
6.本发明为实现上述目的所采用的技术方案是：
7.一种基于abaqus二次开发散热器主片的仿真方法，其包括以下步骤：
8.(1)在ug中构建散热器主片的几何模型及模具模型；
9.(2)将其导入hypermesh软件进行前处理网格划分，得到inp文件导入abaqus求解器中；
10.(3)通过拉伸实验得到铝合金3003材料的性能参数和应力应变曲线；
11.(4)利用fortran语言编写umat子程序，将barlat91屈服准则嵌入至abaqus中；
12.(5)在abaqus求解器中求解并通过试模实验进行本构模型准确性的验证。
13.所述步骤(1)具体包括：在ug软件设计散热器主片的几何模型，根据经验调整凹凸模圆角，模具间隙，压边间隙等几何参数，将在ug中设计好的凹模、凸模、压料片和片体分别导出igs格式并保存。
14.所述步骤(2)具体包括：
15.(2.1)依次将凹模、凸模、压料片和片体igs三维格式导入hypermesh软件，由于igs
的转换误差，先进行模型修复，保证模型网格划分质量；
16.(2.2)进行整体种子布局的设置，然后在圆角处和加强筋处进行局部加密网格，进行前处理网格划分，最后得到inp文件导入abaqus求解器中。
17.所述步骤(3)具体包括：
18.(3.1)根据国标设计拉伸试样，通过激光切割得到1.5mm厚的材料样板，再将样板按照0
°
、45
°
、90
°
三个方向进行切割，得到拉伸样条；
19.(3.2)根据单轴拉伸试验获取各组分材料的工程应力应变数据，并经公式转换为真实应力应变数据；根据单轴拉伸试验获取各组分材料在0
°
、45
°
和90
°
三个方向的塑性应变比；
20.(3.3)根据上述材料参数，利用遗传算法求解barlat91的各向异性系数。
21.当试验是为了获得材料应力应变关系时，试验方法应为将试样拉伸至断裂，拉伸速率为5mm/min；当试验是为了获得材料塑性应变比时，试验方法应为控制试样的拉伸变形达到20％，拉伸速率为1mm/min。
22.所述步骤(4)具体包括：
23.(4.1)将abaqus、visual studio和intel visual fortran三个软件进行关联，为子程序编写和调用提供环境，其中abaqus为分析软件，visual studio为代码运行环境，intel visual fortran为编译器；
24.(4.2)读取当前增量步对应的材料参数、开始应力张量、应变增量和状态变量。其中材料参数包括弹性模量e、泊松比u、硬化参数(a、σs和n)和屈服准则各向异性系数，状态变量为等效塑性应变
25.(4.3)判断当前增量步是否进入弹塑性阶段，根据实际更新阶段求解弹塑性刚度矩阵和更新状态变量，采用改进欧拉算法进行应力更新；
26.(4.4)完成后得到一个用户子程序文件，用于植入后续仿真计算。
27.所述步骤(5)具体包括：
28.(5.1)在abaqus求解器中设置分析步、边界条件、载荷，调用umat子程序文件提交计算；
29.(5.2)通过试模实验测量厚度值，从而进行本构模型准确性的验证。
30.本发明的有益效果是：本发明的仿真方法，通过abaqus二次开发，将barlat91屈服函数嵌入本构模型，从而实现对铝合金成型质量的高精度预测，能够提高散热器主片起皱的仿真预测精度，能够准确地表征散热器主片工艺参数和成型质量的非线性映射，便于指导实际生产，减少废品率和提高开发效率，减少后期售后服务成本，能广泛应用于铝合金散热器主片冲压成型制造技术领域。
附图说明
31.图1为本发明实施例的具体流程图；
32.图2为本发明实施例切割样板示意图；
33.图3为本发明实施例拉伸样条尺寸图；
34.图4为本发明实施例样条三个方向应力应变曲线拟合；
35.图5为本发明实施例三维示意图；
36.图6为本发明实施例abaqus后厚度结果云图；
37.图7为主片试模产品图；
38.图8为主片试模切割后验证位置点示意图；
39.图9为本构模型仿真与实际试模的相对误差图。
具体实施方式
40.实施例：参见图1至图9，本实施例提供一种基于abaqus二次开发散热器主片的仿真方法，其包括以下步骤：
41.(1)在ug中构建散热器主片的几何模型及模具模型；
42.(2)将其导入hypermesh软件进行前处理网格划分，得到inp文件导入abaqus求解器中；
43.(3)通过拉伸实验得到铝合金3003材料的性能参数和应力应变曲线；
44.(4)利用fortran语言编写umat子程序，将barlat91屈服准则嵌入至abaqus中；
45.(5)在abaqus求解器中求解并通过试模实验进行本构模型准确性的验证。
46.所述步骤(1)具体包括：在ug软件设计散热器主片的几何模型，可以看到在散热器主片的中间部位存在许多压包，成型结构复杂，有较多加强筋及转角区域。根据经验调整凹凸模圆角，模具间隙，压边间隙等几何参数，将在ug中设计好的凹模、凸模、压料片和片体分别导出igs格式并保存。
47.所述步骤(2)具体包括：
48.(2.1)依次将凹模、凸模、压料片和片体igs三维格式导入hypermesh软件，由于igs的转换误差，先进行模型修复，保证模型网格划分质量；
49.(2.2)进行整体种子布局的设置，然后在圆角处和加强筋处进行局部加密网格，进行前处理网格划分，最后得到inp文件导入abaqus求解器中。
50.所述步骤(3)具体包括：
51.(3.1)根据国标设计拉伸试样，如图3所示，通过激光切割得到1.5mm厚的材料样板，再将样板按照0
°
、45
°
、90
°
三个方向进行切割，如图2所示，得到拉伸样条；
52.(3.2)所采用的拉伸设备为岛津万能材料试验机ag-x，应力应变拉伸数据按0.01s的时间间隔进行采集。根据单轴拉伸试验获取各组分材料的工程应力应变数据，并经公式转换为真实应力应变数据；如图4所示，根据单轴拉伸试验获取各组分材料在0
°
、45
°
和90
°
三个方向的塑性应变比，如表1所示；
53.表1 3003铝合金塑性应变比
54.方向0
°
45
°
90
°
塑性应变比r0.810.700.78
55.(3.3)根据表2材料参数，利用遗传算法求解barlat91的各向异性系数。
56.表2材料参数
57.角度屈服强度/mpa抗拉强度/mpa均匀伸长率0
°
42.15111.780.26545
°
41.89105.580.27190
°
40.6298.990.283
58.所述barlat91屈服函数公式为：
[0059][0060]
式中：s1、s2和s3是应力张量矩阵s的特征值；当材料为体心立方材料取6，当材料为面心立方材料取8；
[0061][0062][0063][0064]sxy
＝hσ
xy
[0065]szx
＝gσ
zx
[0066]syz
＝fσ
yz
[0067]
式中：σ
xx
、σ
yy
和σ
zz
为正应力；σ
xy
、σ
zx
和σ
yz
为切应力；a、b、c、f、g和h为材料的各向异性系数；其中f和g由于试验条件无法测量简化为按照各向同性处理(即赋值为1)，拉伸试验沿轧制方向0
°
，45
°
，90
°
的屈服应力σ0，σ
45
，σ
90
以及抗拉强度σb，其应力张量可以表示为：
[0068]
σ0＝[σ0,0,0,0,0,0]
[0069]
σ
45
＝[σ
45
/2,σ
45
/2,0,σ
45
/2,0,0]
[0070]
σ
90
＝[0,σ
90
,0,0,0,0]
[0071]
σb＝[σb,σb,0,0,0,0]
[0072]
将其代入barlat91屈服函数得到四个非线性方程：
[0073][0074][0075][0076][0077]
利用遗传算法进行迭代求解，其数学模型如下：
[0078][0079]
最终求解到的6个各项异性系数如表3所示：
[0080]
表3 barlat91屈服准则各向异性系数
[0081]
参数abcfgh数值0.91120430.61247231.3012608110.5496971
[0082]
当试验是为了获得材料应力应变关系时，试验方法应为将试样拉伸至断裂，拉伸速率为5mm/min；当试验是为了获得材料塑性应变比时，试验方法应为控制试样的拉伸变形达到20％，拉伸速率为1mm/min。
[0083]
所述步骤(4)具体包括：
[0084]
(4.1)将abaqus、visual studio和intel visual fortran三个软件进行关联，为子程序编写和调用提供环境，其中abaqus为分析软件，visual studio为代码运行环境，intel visual fortran为编译器；
[0085]
(4.2)读取当前增量步对应的材料参数、开始应力张量、应变增量和状态变量。其中材料参数包括弹性模量e、泊松比u、硬化参数a、屈服应力σs、硬化指数n和barlat91屈服准则对应的六个各向异性系数，状态变量为等效塑性应变
[0086]
(4.3)判断当前增量步是否进入弹塑性阶段，根据实际更新阶段求解弹塑性刚度矩阵和更新状态变量，采用改进欧拉算法进行应力更新；
[0087]
屈服函数的建立如下：
[0088]
屈服准则和硬化准则分别为barlat91屈服准则和ludwik模型，屈服函数如下所示：
[0089][0090]
式中，σ
ij
为应力张量；为等效应力；为微元等效塑性应变；σs为初始屈服应力。
[0091]
应力更新方式
[0092]
dε＝dεe+dε
p
[0093]
式中，dεe为微元弹性应变增量；dε
p
为微元塑性应变增量；dε为微元应变增量。
[0094]
应力更新如下：
[0095]
σ
m+1
＝σm+dσ
[0096]
dσ＝de(dε-dε
p
)
[0097]
式中，σm为初始时的应力张量；σ
m+1
为结束时的应力张量；dσ为微元应力增量；
[0098]de
为弹性刚度矩阵，其一般由弹性模量e和泊松比u求得，如下所示：
[0099]
[0100]
该微元应变增量dε对应的微元塑性应变增量dε
p
[0101][0102]
屈服函数对应力的偏导公式
[0103][0104]
式中，为微元塑性应变增量分量；dλ一般等于微元等效塑性应变增量
[0105]
判断当前增量步是否进入弹塑性阶段，若为完全弹性阶段则子程序调用结束，若为完全弹塑性阶段，采用改进欧拉算法更新应力，同时进行状态变量和弹塑性刚度矩阵的更新。
[0106]
第一步，获取abaqus主程序给出的变量数据，主要包括当前增量步开始时的应力张量σn、当前增量步对应的应变增量δεn、等效塑性应变和相对精度δ。
[0107]
第二步，将应变增量δεn分解成m个微元应变增量dε，求解m：
[0108]
σa＝σn+de(
△
εn‑△
ε
pn
)
[0109][0110][0111]
式中，δε
pn
为将应力张量σn代入式(3-41)求得的塑性应变增量；同理δε
pa
为应力张量σa对应的塑性应变增量。
[0112]
第三步：求微元应变增量dε：
[0113][0114]
第四步：对每个微元应变增量dε进行应力更新：
[0115]
σf＝σi+de(dε
i-dε
pi
)
[0116][0117]
式中，dεi代表应变增量δεn中的第i个微元应变增量；σi代表第i个微元应变增量开始时对应的应力张量；σ
i+1
代表第i个微元应变增量结束时对应的应力张量；dε
pi
为将应力张量σi代入式(3-41)求得的微元塑性应变增量；同理dε
pf
为应力张量σf对应的微元塑性应变增量。
[0118]
第五步：根据式(3-56)将m个微元应变增量都进行计算得到当前增量步结束时对应的应力张量σ
n+1
：
[0119][0120]
否则利用二分法求解比例因子r，
[0121]
由比例因子r可得弹性阶段对应的应变增量δε1：
[0122]
δε1＝rδε
[0123]
将应力更新到弹塑性阶段与弹性阶段的边界处：
[0124]
σ
div
＝σn+derδε
[0125]
式中，σ
div
为弹性阶段与弹塑性阶段分界处的屈服应力张量。
[0126]
得关于比例因子r的等式：
[0127][0128]
对应的零点就是比例因子r的解。采用传统的二分法求解，最后输出r。
[0129]
(4.4)完成后得到一个用户子程序文件，用于植入后续仿真计算。
[0130]
所述步骤(5)具体包括：
[0131]
(5.1)在求解器中设置分析步为静态隐式分析、设置位移载荷约束，下模固定，压料面行程25mm，凸模行程50mm，设置常态库伦摩擦模型，最后调用umat子程序文件提交计算，图6为abaqus仿真结果厚度云图；
[0132]
(5.2)通过表4中工艺参数进行试模实验测量厚度值，如图7所示，从而选取5个位置点，如图8所示进行本构模型准确性的验证，表5为试模和仿真验证结果。
[0133]
表4汽车散热器主片拉深成形工艺参数
[0134][0135]
表4中的模具圆角为凸模上的结构圆角.表5试模和仿真验证结果
[0136][0137]
以上所述，仅为本发明的较佳可行实施例，并非用以局限本发明的专利范围，故凡运用本发明说明书内容所作的方法步骤变化，均包含在本发明的保护范围。

技术特征：
1.一种基于abaqus二次开发散热器主片的仿真方法，其特征在于，其包括以下步骤：(1)在ug中构建散热器主片的几何模型及模具模型；(2)将其导入hypermesh软件进行前处理网格划分，得到inp文件导入abaqus求解器中；(3)通过拉伸实验得到铝合金3003材料的性能参数和应力应变曲线；(4)利用fortran语言编写umat子程序，将barlat91屈服准则嵌入至abaqus中；(5)在abaqus求解器中求解并通过试模实验进行本构模型准确性的验证。2.根据权利要求1所述的基于abaqus二次开发散热器主片的仿真方法，其特征在于，所述步骤(1)具体包括：在ug软件设计散热器主片的几何模型，根据经验调整凹凸模圆角，模具间隙，压边间隙等几何参数，将在ug中设计好的凹模、凸模、压料片和片体分别导出igs格式并保存。3.根据权利要求1所述的基于abaqus二次开发散热器主片的仿真方法，其特征在于，所述步骤(2)具体包括：(2.1)依次将凹模、凸模、压料片和片体igs三维格式导入hypermesh软件，并进行模型修复，保证模型网格划分质量；(2.2)进行整体种子布局的设置，然后在圆角处和加强筋处进行局部加密网格，进行前处理网格划分，最后得到inp文件导入abaqus求解器中。4.根据权利要求1所述的基于abaqus二次开发散热器主片的仿真方法，其特征在于，所述步骤(3)具体包括：(3.1)根据国标设计拉伸试样，通过激光切割得到1.5mm厚的材料样板，再将样板按照0
°
、45
°
、90
°
三个方向进行切割，得到拉伸样条；(3.2)根据单轴拉伸试验获取各组分材料的工程应力应变数据，并经公式转换为真实应力应变数据；根据单轴拉伸试验获取各组分材料在0
°
、45
°
和90
°
三个方向的塑性应变比；(3.3)根据上述材料参数，利用遗传算法求解barlat91的各向异性系数。5.根据权利要求1所述的基于abaqus二次开发散热器主片的仿真方法，其特征在于，所述步骤(4)具体包括：(4.1)将abaqus、visual studio和intel visual fortran三个软件进行关联，为子程序编写和调用提供环境，其中abaqus为分析软件，visual studio为代码运行环境，intel visual fortran为编译器；(4.2)读取当前增量步对应的材料参数、开始应力张量、应变增量和状态变量。其中材料参数包括弹性模量e、泊松比u、硬化参数(a、σ
s
和n)和屈服准则各向异性系数，状态变量为等效塑性应变(4.3)判断当前增量步是否进入弹塑性阶段，根据实际更新阶段求解弹塑性刚度矩阵和更新状态变量，采用改进欧拉算法进行应力更新；(4.4)完成后得到一个用户子程序文件，用于植入后续仿真计算。6.根据权利要求1所述的基于abaqus二次开发散热器主片的仿真方法，其特征在于，所述步骤(5)具体包括：(5.1)在abaqus求解器中设置分析步、边界条件、载荷，调用umat子程序文件提交计算；(5.2)通过试模实验测量厚度值，从而进行本构模型准确性的验证。

技术总结
本发明公开了一种基于Abaqus二次开发散热器主片的仿真方法，其包括以下步骤：(1)在UG中构建散热器主片的几何模型及模具模型；(2)将其导入Hypermesh软件进行前处理网格划分，得到inp文件导入Abaqus求解器中；(3)通过拉伸实验得到铝合金3003材料的性能参数和应力应变曲线；(4)利用Fortran语言编写UMAT子程序，将Barlat91屈服准则嵌入至Abaqus中；(5)在Abaqus求解器中求解并通过试模实验进行本构模型准确性的验证。通过本发明的仿真方法，能够提高散热器主片起皱的仿真预测精度，能够准确地表征散热器主片工艺参数和成型质量的非线性映射，便于指导实际生产，减少废品率和提高开发效率，能广泛应用于铝合金散热器主片冲压成型制造技术领域。压成型制造技术领域。压成型制造技术领域。


技术研发人员：袁熙 余晖跃 李耀聪
受保护的技术使用者：瀚蕾智造（广东）科技有限公司
技术研发日：2023.04.03
技术公布日：2023/8/14