1.本发明属于金属材料铸件微观组织及力学性能仿真技术领域，具体涉及一种基于实验数据的金属铸件不均匀微观组织及力学性能云图计算方法，即基于冷速不均匀性计算金属铸件不均匀微观组织及力学性能分布云图的方法。
背景技术：2.铸造是一种低成本、高效率的金属复杂结构件的成形技术，广泛应用于金属结构件制造领域；绝大多数金属铸件均存在壁厚不均匀问题。金属液注入型腔后，这种壁厚不均匀性必然导致铸件各部位熔体的冷却速率存在差异，而金属熔体冷却速率的差异则会导致铸件微观组织及力学性能的不均匀。因此，因熔体冷却速率差异导致的微观组织及力学性能不均匀性普遍存在于金属铸件(尤其是金属型铸造件)中。
3.零部件微观组织及力学性能的不均匀性对铸件的使役性能、受力仿真、结构设计等至关重要。例如，铸件微观组织的不均匀分布会严重影响铸件的抗疲劳性能；另一方面，基于有限元方法的零部件受力仿真及结构设计，大多基于零部件力学性能均匀的假设，忽略铸件力学性能的不均匀性必然会大幅降低受力仿真及结构设计的可靠性及精度。因此，如何定量计算并预测铸件不均匀微观组织及力学性能分布，并绘制关键微观组织参量及力学性能参量在铸件不同位置分布的云图，有利于对铸件的使役性能进行预测、结构设计进行优化。
4.对于不均匀微观组织的模拟，目前的方法有元胞自动机。例如，公告日为2021年3月12日、公告号为cn110245449b的中国专利文件公开了一种镁合金铸造件成分不均匀性数值预测方法，其中采用元胞自动机法模拟具有不同生长取向的α-mg枝晶生长，获得枝晶比表面积随固相分数变化曲线。目前的研究多集中与一维和二维元胞自动机，对于三维实体铸件还不能整体预测。二维元胞自动机通常以规则的空间单元(如四方格网、等边三角形等)划分，按照von
·
neumann型、moore型和margolus型等方法进行规则处理，这直接导致元胞状态转换规则一般不能应用于更远的单元。
5.对于预测力学性能的方法，目前可以通过组织预测。例如，公告日为2015年9月9日、公告号为cn104899412a的中国专利公开了一种铝合金铸件力学性能预测方法，其中采用孔隙率与二次枝晶臂间距来预测力学性能，该方法中孔隙率、二次枝晶臂间距与力学性能的关系是通过实验得到的，但孔隙率、二次枝晶臂间距是通过建立凝固过程的数学模型得到的，并未与实际实验结合。
6.目前的微观组织预测和力学性能预测都没有与实验结合，且理论参数难以获取，无法精确三维预测。在零部件的实际铸造过程中无法观察到冷速，也无法判断零部件的微观组织和力学性能。
7.因此，亟需探索一种新的金属铸件不均匀微观组织及力学性能云图计算方法。
技术实现要素：8.本发明针对金属铸件不均匀组织及力学性能分布难以精确、定量预测的难题，提供了一种基于实验数据的金属铸件不均匀微观组织及力学性能云图计算方法，具体是基于铸件熔体冷却速率与微观组织参量及力学性能参量间定量关系，结合铸件充型模拟，计算金属铸件不均匀微观组织参量及力学性能(强度、塑性)参量的方法。
9.本发明的发明构思：
10.首先通过实验和/或理论计算的方法，建立熔体冷却速率与微观组织参量及力学性能参量间的定量关系；然后，通过模流分析软件(easycast、procast等)模拟充型过程中铸件不同位置的冷却速率，将包含铸件位置坐标信息及冷却速率的原始文件导出；进一步通过编程，基于所获得的熔体冷却速率与微观组织参量及力学性能参量间的定量关系、铸件位置坐标信息及冷却速率，计算得到铸件不同位置对应的微观组织参量及力学性能；最后，通过matlab编程，计算获得铸件微观组织及力学性能的定量分布云图。
11.为实现上述目的，本发明所提供的技术解决方案是：
12.一种基于实验数据的金属铸件不均匀微观组织及力学性能云图计算方法，其特殊之处在于：
13.s1、建立样品微观组织参量以及力学性能参量与熔体冷却速率间的定量关系：
14.s1.1制备不同冷速的实验样品：
15.通过设计不同壁厚的浇铸模具，使其冷速分布在0～120k/s之间，获得7组以上不同冷却速率的实验样品；
16.s1.2对s1.1制备出的7组以上不同冷却速率的实验样品分别进行微观组织测量，获得各微观组织参量与熔体冷却速率间的对应数据，进一步通过数据分析(在origin软件里进行数据分析)，建立微观组织参量与熔体冷却速率间的定量关系；
17.对s1.1制备出的7组以上不同冷却速率的实验样品分别进行标准力学性能测试，获得力学性能参量与熔体冷却速率间的对应数据，进一步通过数据分析(在origin软件里进行数据分析)，建立力学性能参量与熔体冷却速率间的定量关系；
18.s2、在模流分析软件里获得铸件位置坐标与冷却速率的对应关系(铸件与实验样品材料相同)：
19.通过模流分析软件，分析不同铸造工艺条件下，熔体充满铸件型腔后，铸件不同位置的温度场变化(不同的铸造工艺是指重力铸造、低压铸造等这种铸造条件不同的工艺，在这些模流分析软件中通过设置铸造参数可以得到不同铸造工艺的温度场)，选择铸件凝固发生前熔体的冷却速率作为各位置的冷却速率，获得铸件位置坐标与冷却速率的对应关系(具体在模流软件里获得铸件模型与冷却速率的云图，具体的数据需要利用matlab软件提取)，获得零部件铸件三维冷速分布；
20.s3、利用s1建立的微观组织参量与熔体冷却速率间的定量关系、力学性能参量与熔体冷却速率间的定量关系，以及s2得到铸件位置坐标与冷却速率的对应关系，通过matlab程序，获得铸件位置坐标与微观组织参量和力学性能参量的对应关系(即得到铸件不同位置对应的微观组织参量及力学性能)；之后，利用scatter3的函数，获得金属铸件不均匀微观组织及力学性能云图。
21.具体是：利用matlab本身自带的函数进行编程提取铸件位置坐标和冷却速率信
息，并代入s1建立的关系，进而得到铸件位置坐标与对应的微观组织参量和力学性能参量，最后再利用scatter3的函数，将铸件位置坐标与微观组织参量和力学性能参量的对应关系在三维空间内作图，可获得各参量在铸件整体空间内的分布，至此获得金属铸件不均匀微观组织及力学性能云图。
22.进一步地，s1中，所述微观组织参量是指二次枝晶臂间距；
23.所述力学性能参量是指屈服强度或塑性。
24.进一步地，s2中，所述模流软件为easycast、procast。
25.进一步地，所述s1和s2顺序可以调换。
26.本发明的优点是：
27.1.本发明提出了一种铸件不均匀微观组织及力学性能在铸件空间分布的计算和预测方法。与当前商用软件中的微观组织预测方法相比，基于实验测量建立的组织及力学参量与冷却速率的定量关系具有可靠性高、精度高，可避免理论模型计算过程中计算参数难获取，参数值与真实值偏离大的问题，可更为有效的实现组织及力学参数的计算及预测。
28.2.本发明应用过程中，使用者还可根据实际工况，针对不同材料自主建立数据库，通过数据分析，修正和建立组织及力学参量与冷却速率的定量关系，在实际应用过程中不断优化此定量关系，进一步提高计算精度及可靠性。
29.3.基于本发明计算的力学性能云图数据还可通过数据编程导入abaqus等商用有限元仿真软件，进行基于力学性能不均匀分布的铸件的受力仿真及结构优化设计。相比基于假设铸件力学性能处处相同的受力分析及结构优化设计，这一优化处理可显著提升力学仿真的精度，从而更为有效的指导铸件的结构设计。
30.4.利用本发明方法，可定量计算出金属铸件微观组织参量及力学性能参量在铸件不同位置的分布，可为铸件工艺设计、使役性能预测及铸件结构优化设计提供有效支撑。
31.5.本发明是在模流分析软件中获得零部件的冷速，再从实验中获得微观组织参量与冷却速率、力学性能参量与冷却速率的定量关系，最终获得零部件的微观组织与力学性能云图。本发明提供方法的精确度在于在模流分析软件中对零部件网格划分的大小，网格越细则越精确，尝试过最小的网格大小为0.2mm。
附图说明
32.图1为本发明一种基于实验数据的金属铸件不均匀微观组织及力学性能云图计算方法的实施流程图。
33.图2为实施例1和实施例2中实验测量a356铸造铝合金二次枝晶臂间距与熔体冷却速率的对应关系及通过数据拟合获得的定量关系模型。(参考文献为附录1)
34.图3为实施例1和实施例2中实验测量a356铸造铝合金屈服强度与熔体冷却速率的对应关系及通过数据拟合获得的定量关系模型。
35.图4为实施例1中基于easycast模流分析软件获得铸件冷却速率分布后，计算得到的某典型a356铝合金铸件的二次枝晶臂间距及屈服强度分布云图。
36.图5为实施例2中基于procast模流分析软件获得铸件冷却速率分布后，计算得到的某典型a356铝合金铸件的二次枝晶臂间距及屈服强度分布云图。
具体实施方式
37.以下结合附图和具体实施例对本发明的内容作进一步的详细描述：实施例1和实施例2都是在重力铸造下获得的a356铸造铝合金微观组织参量与冷速、力学性能参量与冷速的定量关系下获得云图。通过本发明的方法，可改变铸造方法和合金来获得其他铸造合金零部件的微观组织和力学性能云图。
38.本发明一种基于实验数据的金属铸件不均匀微观组织及力学性能云图计算方法，包括以下步骤：
39.s1、建立样品微观组织参量以及力学性能参量与熔体冷却速率间的定量关系
40.s1.1制备不同冷速的实验样品：
41.通过设计不同壁厚的浇铸模具，使其冷速分布在0～120k/s之间，获得7组以上不同冷却速率的实验样品；
42.s1.2对s1.1制备出的7组以上不同冷却速率的实验样品分别进行微观组织测量，获得各微观组织参量与熔体冷却速率间的对应数据，进一步通过数据分析(在origin软件里进行数据分析)，建立微观组织参量与熔体冷却速率间的定量关系；所述微观组织参量是指二次枝晶臂间距；
43.对s1.1制备出的7组以上不同冷却速率的实验样品分别进行标准力学性能测试，获得力学性能参量与熔体冷却速率间的对应数据，进一步通过数据分析(在origin软件里进行数据分析)，建立力学性能参量与熔体冷却速率间的定量关系；所述力学性能参量是指屈服强度或塑性；
44.s2、在模流分析软件(easycast、procast)里获得零部件铸件(与实验样品材料相同)位置坐标与冷却速率的对应关系：
45.通过模流分析软件，分析不同铸造工艺条件下，熔体充满铸件型腔后，铸件不同位置的温度场变化(不同的铸造工艺是指重力铸造、低压铸造等这种铸造条件不同的工艺，在这些模流分析软件中通过设置铸造参数可以得到不同铸造工艺的温度场)，选择铸件凝固发生前熔体的冷却速率作为各位置的冷却速率，获得铸件位置坐标与冷却速率的对应关系(具体在模流软件里获得铸件模型与冷却速率的云图，具体的数据需要利用matlab软件提取)，获得零部件铸件三维冷速分布；
46.s3、利用s1建立的微观组织参量与熔体冷却速率间的定量关系、力学性能参量与熔体冷却速率间的定量关系，以及s2得到铸件位置坐标与冷却速率的对应关系，通过matlab程序，获得铸件位置坐标与微观组织参量和力学性能参量的对应关系(即得到铸件不同位置对应的微观组织参量及力学性能)；之后，利用scatter3的函数，获得金属铸件不均匀微观组织及力学性能云图。
47.具体是：利用matlab本身自带的函数进行编程提取铸件位置坐标和冷却速率信息，并代入s1建立的关系，进而得到铸件位置坐标与对应的微观组织参量和力学性能参量，最后再利用scatter3的函数，将铸件位置坐标与微观组织参量和力学性能参量的对应关系在三维空间内作图，可获得各参量在铸件整体空间内的分布，至此获得金属铸件不均匀微观组织及力学性能云图。上述s1和s2顺序可以调换。
48.实施例1
49.基于easycast模流分析软件，计算某典型a356铝合金铸件的二次枝晶臂间距及屈
服强度分布云图，具体步骤包括：
50.步骤1、从easycast完成凝固模拟后，获得冷速文件.vel文件。
51.步骤2、分析.vel文件的数据：
52.步骤2-1、从.vel文件的第二行获得网格大小，实施例1的网格大小为0.5mm。
53.步骤2-2、easycast是有限差分法，网格是六面体，从.vel文件的第三行开始提供的数据可以确定中心点的坐标和冷速，例如“5 46 5 9.677527”代表的是x方向上第5个0.5mm、y方向上第46个0.5mm、z方向上第5个0.5mm的位置，9.677527代表的是冷速，因此用中心点的冷速来代表整个网格的冷速。
54.步骤3、删掉.vel文件的前两行和最后四行，只留下能确定中心点的数据和冷速的四列数据，获得基础信息文件lengsu.txt。
55.步骤4、打开easycast_exchang.m的程序文件(附录2)，这个程序提取中心点的信息和冷速，将冷速代入二次枝晶臂间距与熔体冷却速率的关系式以及屈服强度与熔体冷却速率的关系式，如图2和图3所示。修改基础信息文件lengsu.txt、二次枝晶臂间距文件sdas_easycast.txt和屈服强度文件qdt_easycast.txt的路径、文件名和文件格式，点击运行后就可以得到带有二次枝晶臂间距和屈服强度信息的模型文件。
56.步骤5、打开easycast_yuntu.m的程序文件(附录3)，将数据文件(二次枝晶臂间距sdas_easycast.txt、屈服强度qdt_easycast.txt)放在同一文件夹下，修改load后的文件名和文件格式，在x、y、z后修改网格大小，设置刻度轴范围、视图视角以及标题，最后运行得到二次枝晶臂间距及屈服强度分布云图，如图4所示。
57.实施例2
58.基于procast模流分析软件，计算某典型a356铝合金铸件的二次枝晶臂间距及屈服强度分布云图，具体步骤包括：
59.步骤1、在procast获得节点(node)的编号和坐标。铸造模拟完成后，在viewer版块打开g.unf文件，然后在左上角点击mesh，在file下拉菜单点击export，输出g.inp文件。
60.步骤2、在procast获得冷速文件。在procast的viewer版块，在results的下拉菜单点击metallurgical tools，选择r,g,l，根据nyiama criterion计算准则，a填1，b填0，c填1，d填-0.5，l upper temperature填615(tliquidus+2)，l lower temperature填548(tsolidus)，r,gtemperature填554.5(tsolidus+0.1*(tliquidus-tsolidus))，然后计算得到procast的铸造模拟冷速分布，在selection菜单栏里选择node，框选住模型(只选择铸件)，在file下拉菜单点击export as，然后选择patran,i-deas,stl,g3d
…
这一项，输出带有节点编号的冷速文件.ntl。
61.步骤3、打开冷速文件.ntl，删掉前四行，获得只有节点编号和冷速两列信息的冷速文件lengsu_pro.txt。
62.步骤4、获得二次枝晶臂间距文件和屈服强度文件。打开procast_exchange的程序文件(附录4)，这个程序提取节点(node)的坐标和冷速，将冷速代入二次枝晶臂间距与熔体冷却速率的关系式以及屈服强度与熔体冷却速率的关系式，如图2和图3所示。修改g.inp文件、冷速文件lengsu_pro.txt、二次枝晶臂间距文件sdas_procast.txt和屈服强度文件qdt_procast.txt的路径、文件名和文件格式，点击运行后就可以得到带有二次枝晶臂间距和屈服强度信息的模型文件。
63.步骤5、打开procast_yuntu.m的程序文件(附录5)，将数据文件(二次枝晶臂间距sdas_procast.txt、屈服强度qdt_procast.txt)放在同一文件夹下，修改load后的文件名和文件格式，设置刻度轴范围、视图视角以及标题，最后运行得到二次枝晶臂间距及屈服强度分布云图，如图5所示。
64.实施例1和实施例2在获得二次枝晶臂间距及屈服强度分布云图中需要用到的微观组织参量以及力学性能参量与熔体冷却速率间的定量关系，均按照前述方法获得。
65.综上，利用本发明所提出的一种铸件不均匀微观组织及力学性能在铸件空间分布的计算和预测方法，将金属铸件不均匀微观组织参量及力学性能(强度、塑性)参量可视化。
66.图2的参考文献
67.[1]kai z,yucheng s,suolin w,et al.effects of cooling rate on high-temperature mechanical properties of a356 alloys[j].rare metal materials and engineering,2011,40(s3):63-68.
[0068]
[2]ni h,sun b,jiang h,et al.effect of jdn-i flux on das of a356 alloy at different cooling rate[j].materials ence&engineering a,2003,348(1-2):1-5.
[0069]
[3]刘颖卓,党波,刘峰.冷却速率对a356铝合金显微组织和微观硬度的影响[j].西安工业大学学报,2013(02):128-133.
[0070]
[4]s,boontein,n,et al.reduction in secondary dendrite arm spacing in cast aluminium alloy a356 by sb addition[j].international journal of cast metals research,2011.
[0071]
[5]党波.a356铝合金非平衡凝固与固态相变一体化研究[d].西北工业大学,2017.
[0072]
附录2
[0073]
easycast_exchang.m程序代码
[0074]
easycast_exchang.m
[0075]
clear all；
[0076]
clc；
[0077]
line_number＝112828；％冷速文件lengsu.txt总行数
[0078]
m＝1；
[0079]
n＝0；
[0080]
a＝{}；
[0081]
fw1＝fopen('lengsu.txt','r')；％打开冷速文件lengsu.txt为只读
[0082]
fw6＝fopen('sdas_easycast.txt','a+')；％新建储存二次枝晶臂信息和模型信息的文件sdas_easycast.txt
[0083]
fw7＝fopen('qdt_easycast.txt','a+')；％新建储存屈服强度信息和模型信息的文件
[0084]
qdt_easycast.txt
[0085]
while m《＝line_number
[0086]
n＝n+length(a)；
[0087]
fseek(fw1,n,-1)；％指针
[0088]
a＝fgets(fw1)；％取一行数字
[0089]
b＝regexp(a,”,'split')；％通过空格分隔数据得到一个cell
[0090]
c＝str2num(char(b))；％通过str2num和char函数转换格式
[0091]
sdas(m)＝3*10^(-5)*(c(4))^(-0.302)；％sdas与冷速的公式，sdas单位是m，冷速单位是k/s
[0092]
sdas(m)＝sdas(m)*10^6；％修改单位为μm
[0093]
qft(m)＝0.0011*(c(4))^2+0.3998*c(4)+171.12；％屈服强度单位是mpa
[0094]
fprintf(fw6,'％d％d％d％6.3f\n',c(1),c(2),c(3),sdas(m))；
[0095]
fprintf(fw7,'％d％d％d％6.3f\n',c(1),c(2),c(3),qft(m))；
[0096]
m＝m+1；
[0097]
end
[0098]
fclose('all')；
[0099]
附录3：
[0100]
easycast_yuntu.m程序代码
[0101]
easycast_yuntu.m
[0102]
clear all；
[0103]
a＝load('qdt_easycast.txt')；％载入储存屈服强度信息和模型信息的信息文件
[0104]
x＝a(:,1)*0.5；％中心点x轴的坐标
[0105]
y＝a(:,2)*0.5；％中心点y轴的坐标
[0106]
z＝a(:,3)*0.5；％中心点z轴的坐标
[0107]
v＝a(:,4)；％屈服强度
[0108]
scatter3(x,y,z,20,v,'filled','s')；％20代表散点的大小，’s’代表散点的形状用立方体表示
[0109]
ax＝gca；％返回当前图形的轴或表
[0110]
axis equal；％保持刻度轴大小相同
[0111]
axis([0 60 0 60 0 20])；％设置刻度轴范围，分别是x轴最小值、x轴最大值、y轴最小值、y轴最大值、z轴最小值、z轴最大值
[0112]
view(30,60)；％调整视图的方位角和仰角
[0113]
xlabel('x')；％为x轴添加标签
[0114]
ylabel('y')；％为y轴添加标签
[0115]
zlabel('z')；％为z轴添加标签
[0116]
colorbar；％添加色标
[0117]
title('ys(mpa)')；％设置图的标题
[0118]
附录4：procast_exchang.m程序代码
[0119]
procast_exchang.m
[0120]
clear all；
[0121]
clc；
[0122]
fw1＝fopen('g.inp','r')；％打开g.inp文件
[0123]
line_number＝73926；％g.inp文件的行数
[0124]
z＝{}；
[0125]
w＝{}；
[0126]
p＝0；
[0127]
o＝1；
[0128]
while o《＝line_number％此循环是为了寻找g.inp文件里*node的位置
[0129]
if isempty(z)
[0130]
p＝p+length(w)；
[0131]
fseek(fw1,p,-1)；％指针
[0132]
w＝fgets(fw1)；％取一行数字
[0133]
z＝strfind(w,'*node')；％寻找*node
[0134]
end
[0135]
o＝o+1；
[0136]
end
[0137]
p＝p+length(w)；
[0138]
nodenum＝13643；％node的数量
[0139]
fw5＝fopen('sdas_procast.txt','a+')；％新建储存二次枝晶臂信息和模型信息的文件sdas_procast.txt
[0140]
fw6＝fopen('qdt_procast.txt','a+')；％新建储存屈服强度信息和模型信息的文件
[0141]
qdt_procast.txt
[0142]
y＝1；
[0143]
cr＝load('lengsu_pro.txt')；％载入procast的冷速文件
[0144]
while y《＝nodenum％此循环是把冷速提取出来代入屈服强度与熔体冷却速率的关系式得到屈服强度
[0145]
fseek(fw1,p,-1)；％指针
[0146]
c＝fgets(fw1)；％取一行数字
[0147]
p＝p+length(c)；
[0148]
jg＝regexp(c,',','split')；％通过空格分隔数据得到一个cell
[0149]
d＝str2num(char(jg))；％通过str2num和char函数转换格式
[0150]
ncr＝cr(y,2)；
[0151]
sdas(y)＝3*10^(-5)*(ncr)^(-0.302)；％％sdas与冷速的公式，sdas单位是m，冷速单位是k/s
[0152]
sdas(y)＝sdas(y)*10^6；％修改单位为μm
[0153]
qft(y)＝0.0011*(ncr)^2+0.3998*ncr+171.12；％屈服强度单位是mpa
[0154]
fprintf(fw5,'％9.6f％9.6f％9.6f％8.6f\n',d(2),d(3),d(4),sdas(y))；
[0155]
fprintf(fw6,'％9.6f％9.6f％9.6f％8.6f\n',d(2),d(3),d(4),qft(y))；
[0156]
y＝y+1；
[0157]
end
[0158]
fclose('all')；
[0159]
附录5：procast_yuntu.m程序代码
[0160]
procast_yuntu.m
[0161]
clear all；
[0162]
a＝load('qdt_procast.txt')；％载入储存屈服强度信息和模型信息的信息文件
[0163]
x＝a(:,1)；
[0164]
y＝a(:,2)；
[0165]
z＝a(:,3)；
[0166]
v＝a(:,4)；
[0167]
scatter3(x,y,z,80,v,'filled')；％20代表散点的大小
[0168]
ax＝gca；％返回当前图形的轴或表
[0169]
axis equal；％保持刻度轴大小相同
[0170]
axis([-40 40-40 40-10 30])；％设置刻度轴范围，分别是x轴最小值、x轴最大值、y轴最小值、y轴最大值、z轴最小值、z轴最大值
[0171]
view(30,70)；％调整视图的方位角和仰角
[0172]
xlabel('x')；％为x轴添加标签
[0173]
ylabel('y')；％为y轴添加标签
[0174]
zlabel('z')；％为z轴添加标签
[0175]
colorbar；％添加色标
[0176]
title('ys(mpa)')；％设置图的标题
[0177]
上述说明出示并描述了发明的若干优选实施例，但如前所示，须理解本发明并非局限于本文所披露的形式，不应看作是对其他实施例的排除，而可用于各种其他组合、修改和环境，并能够在本文所述发明构想范围内，通过上述教导或相关领域的技术或知识进行改动。而本领域人员所进行的改动和变化不脱离发明的精神和范围，则都应在发明所附权利要求的保护范围内。
