CN110245449A - 一种镁合金铸造件成分不均匀性数值预测方法 - Google Patents

Abstract

一种镁合金铸造件成分不均匀性数值预测方法，本发明涉及镁合金铸造件成分不均匀性数值预测方法。本发明的目的是为了解决现有方法无法准确预测镁合金铸件宏观偏析形成的问题。过程为：一、采用元胞自动机法模拟具有不同生长取向的α‑Mg枝晶生长，获得枝晶比表面积随固相分数变化曲线；二、对铸造系统进行网格剖分；三、针对所有下角标char＝0的网格，计算能量、成分和动量守恒方程，获得铸件内温度、平均成分和速度的分布；四、针对所有下角标char非0的网格，计算能量守恒方程，获得温度分布；五、重复二、三和四，直到凝固结束，输出铸件内平均成分场。本发明用于镁合金铸造件成分预测领域。

Description

一种镁合金铸造件成分不均匀性数值预测方法

技术领域

本发明涉及镁合金铸造件成分不均匀性数值预测方法。

背景技术

镁合金具有比重轻、比强度和比刚度高及阻尼减震性能好等一系列优点，是轻量化的重要材料。镁合金铸件在国防、航空航天和汽车领域中具有广阔应用前景，铸造法是其成形的主要工艺，因此提高镁合金铸件质量对我国工业发展具有重要意义。成分不均匀是铸件凝固过程中易形成的一种铸造缺陷，也称偏析。宏观尺度上的偏析为宏观偏析，凝固结束后合金成分高于初始合金成分则称为正偏析，低于初始合金成分则称为负偏析。宏观偏析的形成一方面会增加铸件热裂倾向，另一方面导致铸件力学性能不均匀，同时宏观偏析的形成会引起铸件局部凝固滞后，夹杂物会富集该区。后续热处理无法消除宏观偏析，补焊工艺只适用于消除形成在铸件表层的宏观偏析。宏观偏析的形成是铸件废品率增加的一个重要原因。因此控制宏观偏析缺陷，对于提高镁合金铸件成品率具有重要意义。

三种常用的研究宏观偏析的手段是理论计算、实验和数值模拟。常用的宏观偏析理论计算模型为夏尔模型或杠杆模型，这个两个模型可以准确描述固相以及固液界面前液相中溶质分布特点，理论计算更适合一维凝固情况，由于需要对模型进行解析求解，不适合描述铸件二维或三维凝固过程中的宏观偏析形成。采用实验手段揭示不同凝固条件下宏观偏析形成特点需要依赖大量实验数据。由于合金凝固过程不可视，基于实验结果反推具体过程会具有盲目性和不确定性。开展大量实验对人力、物力和能源是一种浪费，同时也延长了研发周期。借助数值模拟方法预测宏观偏析形成在铸造工艺开发或工艺改进过程中具有明显优势。基于经典理论模型和数值方法，数值模拟可以实时再现凝固过程中的温度、成分和动量传输，可以获得宏观偏析形成的起始时间以及不同工艺参数对偏析分布的影响规律。基于不同物理量随时间变化特点，对不同物理量在偏析形成过程中所起的作用进行评价。采用数值模拟手段可以加速工艺开发或工艺改进进程，同时节约人力、物力和能源。

凝固过程中降温速度缓慢、液体流动是宏观偏析形成的主要原因。降温速度慢一方面会存进长时间的液体流动，另一方面会增加枝晶二次臂间距(即增加了枝晶网络的渗透率)，从而促进了溶质的长程迁移，最终引发宏观偏析。目前对宏观偏析的数值模拟通常采用简化方法处理枝晶二次臂间距信息，即基于实验金相测量，获得二次枝晶臂间距特征值，将该值作为输入参数计算渗透率。该方法存在的问题是：实际凝固过程中，多个枝晶沿三维方向生长且生长取向任意，而金相测量法针对试样被切割后的二维截面，对试样进行切割时很难找到某一枝晶的对称面，二维截面上无法准确判断两个先析出相之间的关系(是同属于一个晶粒，还是属于不同晶粒)，因此所获得的二次枝晶臂间距特征值存在很大误差。可以采用微观元胞自动机技术对枝晶形貌进行模拟计算，从而获得二次枝晶臂间距特征值。但是镁为密排六方结构，枝晶生长过程中呈现六重对称特性，模拟计算过程中受计算网格各向异性的影响很大。已有研究中，采用六边形网格可以较好再现六重对称枝晶形貌，但是不能计算枝晶以任意角度生长情况；采用正方形网格可以较好再现枝晶形貌且能计算枝晶以任意角度生长，但是引入的随机捕获因子缺少物理含义。

发明内容

本发明的目的是为了解决现有方法无法准确预测镁合金铸件宏观偏析形成的问题，而提出一种镁合金铸造件成分不均匀性数值预测方法。

一种镁合金铸造件成分不均匀性数值预测方法具体过程为：

步骤一、采用元胞自动机法模拟具有不同生长取向的α-Mg枝晶生长，获得枝晶比表面积随固相分数变化曲线；

步骤二、对铸造系统进行网格剖分，对铸造系统直角坐标系X轴方向和Y轴方向采用相同的网格剖分步长Δx米＝Δy米；

Δx＝X_max-X_min，Δy＝Y_max-Y_min；

计算标号为(js，ks)_char的网格，下角标char＝2表示铸型网格，下角标char＝0表示铸件网格；

其中js和ks均为整数，js为在铸造系统直角坐标系X轴上的坐标，ks为在铸造系统直角坐标系Y轴上的坐标；js的取值范围是1～M，ks的取值范围是1～N；

X_max为铸造系统直角坐标系X轴方向最大值，X_min为铸造系统直角坐标系X轴方向最小值，Y_max为铸造系统直角坐标系Y轴方向最大值，Y_min为铸造系统直角坐标系Y轴方向最小值，单位均为米；

步骤三、针对所有下角标char＝0的网格(js，ks)_char＝0，计算能量、成分和动量守恒方程，获得铸件内温度、平均成分和速度的分布；

步骤四、针对所有下角标char非0的网格(js，ks)_char≠0，计算能量守恒方程，获得温度T_c分布；

T_c为非铸件网格单元的温度，即铸型网格的温度；

步骤五、重复步骤二、步骤三和步骤四，直到凝固结束，即下角标char＝0的计算网格的f_sr(js,ks)_char＝0＝1.0，输出铸件内平均成分场[C]；f_sr为固相分数。

本发明的有益效果为：

微观尺度计算能够提供微观偏析和枝晶形貌信息，因为没有考虑液相流动，因此不能预测宏观偏析。基于微观尺度计算，开展宏观尺度传热、传质和液体流动计算，可以对宏观偏析(一种铸造缺陷)形成进行预测且宏观偏析预测更符合实际。

相比于二次枝晶臂间距，本发明比表面积可以更为准确的表征枝晶网络对流体的阻碍影响。本发明对镁合金铸件凝固过程中宏观偏析形成进行数值模拟，基于正方形网格计算镁枝晶形貌，需引入随机因子，从而获得枝晶比表面积特征值，将该信息引入渗透率计算模型，更为准确的计算宏观偏析形成。

本发明设计了镁合金铸件中宏观偏析数值模拟的方法，基于更加真实的枝晶形貌，计算枝晶比表面积随固相分数变化曲线，从而获得枝晶结构渗透率，更为准确地预测了铸件中宏观偏析形成，解决了目前宏观偏析预测中，基于实验数据获得枝晶二次臂间距存在较大误差的问题，为铸造件成形过程中宏观偏析缺陷控制供了理论指导。

由附图3可见，尽管枝晶生长取向不同，但每个枝晶形貌均呈六重对称状，即一次臂之间的夹角基本接近60°，符合密排六方晶体结构特点。

本发明适用于镁合金各类尺寸铸件的砂型和金属型铸造过程中宏观偏析的预测。利用本发明可以更为准确的预测偏析形成，从多个方面为工艺开发和改进提供理论支持，市场应用潜力巨大，一旦被广泛采用，将有几十亿元以上的产值。

附图说明

图1为本发明正方形网格附近邻居单元定义示意图；

图2为本发明邻居单元之间的配对关系示意图；

图3为本发明模拟所得三个生长取向分别为0°、30°和45°的α-Mg枝晶生长过程示意图；

图4为本发明Mg-4wt％Y合金铸件铸造系统示意图；

图5为本发明实验所得Mg-4wt％Y合金铸件中间截面金相组织照片示意图；

图6为本发明铸件中所取试样金相组织实验照片示意图；

图7为本发明铸件所取试样位置附近的冷却曲线示意图；

图8为本发明模拟所得二维截面晶粒组织示意图；

图9为本发明模拟计算所得平均比表面积随固相分数(f_s)变化曲线示意图；

图10a为Mg-4wt％Y合金铸件宏观偏析分布实验结果示意图；

图10b为本发明基于平均比表面积随固相分数变化曲线，模拟所得Mg-4wt％Y合金铸件宏观偏析分布模拟结果图；

图10c为本发明基于给定枝晶臂间距600μm，模拟所得Mg-4wt％Y合金铸件宏观偏析分布结果图。

具体实施方式

具体实施方式一：本实施方式一种镁合金铸造件成分不均匀性数值预测方法具体过程为：

步骤一、采用元胞自动机法模拟具有不同生长取向的α-Mg枝晶生长，获得枝晶比表面积随固相分数变化曲线；

步骤二、对铸造系统进行宏观尺度网格剖分，对铸造系统直角坐标系X轴方向和Y轴方向采用相同的网格剖分步长Δx米＝Δy米；

Δx＝X_max-X_min，Δy＝Y_max-Y_min；

计算标号为(js，ks)_char的网格，下角标char＝2表示铸型网格，下角标char＝0表示铸件网格，下角标char＝4、5、6、7和8分别表示内冷铁网格、外冷铁网格、冒口套网格、保温材料网格和绝热材料网格；

其中js和ks均为整数，js为在铸造系统直角坐标系X轴上的坐标，ks为在铸造系统直角坐标系Y轴上的坐标；js的取值范围是1～M，ks的取值范围是1～N；

X_max为铸造系统直角坐标系X轴方向最大值，X_min为铸造系统直角坐标系X轴方向最小值，Y_max为铸造系统直角坐标系Y轴方向最大值，Y_min为铸造系统直角坐标系Y轴方向最小值，单位均为米；

步骤三、针对所有下角标char＝0的网格(js，ks)_char＝0，计算能量、成分和动量守恒方程，获得铸件内温度、平均成分和速度的分布；

步骤四、针对所有下角标char非0的网格(js，ks)_char≠0，计算能量守恒方程，获得温度T_c分布；

T_c为非铸件网格单元的温度，即铸型网格的温度；

步骤五、重复步骤二、步骤三和步骤四，直到凝固结束，即下角标char＝0的计算网格的f_sr(js，ks)_char＝0＝1.0，输出铸件内平均成分场[C]；f_sr为固相分数。

具体实施方式二：本实施方式与具体实施方式一不同的是，所述步骤一中采用元胞自动机法模拟具有不同生长取向的α-Mg枝晶生长，获得枝晶比表面积随固相分数变化曲线；具体过程为：

步骤一(1)、将枝晶生长计算域(已经确定的)进行微观尺度网格剖分，采用边长尺寸为Δstep的正方形网格，每一个正方形网格用(j，k)标识，j的取值范围为[1，n]，k的取值范围为[1，m]，j和k均为整数，n、m均为整数，j为在直角坐标系X轴上的坐标，k为在直角坐标系Y轴上的坐标；

步骤一(2)、为计算域中每一个正方形网格赋予邻居对象，通过对邻居对象的捕捉完成枝晶生长；

标号为(j∈[1，2]，k∈[1，m])、(j∈[n-1，n]，k∈[1，m])、(j∈[3，n-2]，k∈[1，2])、(j∈[3，n-2)]，k∈[m-1，m])的正方形网格定义为边界网格，不参与枝晶生长，因此不需赋予邻居对象；

其余每一个标号为(j，k)的正方形网格都有12个邻居：

(j-1，k)_left，(j-2，k)_leftleft，(j+1，k)_right，(j+2，k)_rightright，(j-1，k+1)_upleft，(j-1，k+2)_upleftup，(j+1，k+1)_upright，(j+1，k+2)_uprightup，(j-1，k-1)_downleft，(j-1，k-2)_downleftdown，(j+1，k-1)_downright，(j+1，k-2)_downrighdown；

步骤一(3)、确定凝固过程中每个正方形网格的状态；

步骤一(4)、给定冷却速度单位为℃/s，计算域中每个网格都具有相同的温度，当凝固时间为t_i时，单位为s，网格(j，k)对应的温度为T(j，k，t_i)；

步骤一(5)、等轴枝晶形核核心数为Neg个，由实验测量提供。等轴枝晶生长取向角θ随机选择、等轴枝晶形核位置随机选择，等轴枝晶形核过冷度为0.5℃；

某一时刻t_i，当等轴枝晶形核核心所处网格的温度T(j，k，t_i)＜(T_L-0.5)，等轴枝晶形核瞬间发生；等轴枝晶形核核心所处网格对应的物理量为：f_s＝1，C_s＝C_ok_par，C_l＝0，C_aver＝C_sf_s+C_l(1-f_s)，state＝2，icolor＝Random[1，Neg]，growthθ＝Random[0，90°]；

其中，f_s为固相分数，最小值为0表示液相，最大值为1表示固相；C_s为固相成分；C_o为合金初始成分；k_par为平衡分配系数；C_l为液相成分；C_aver为平均成分，成分的单位为wt％；T_L为合金液相线温度，单位为℃；

每一个等轴枝晶形核核心赋予一个标号icolor，用来区分不同晶粒，Random[1，Neg]表示在1～Neg之间选取一个随机数；

state为状态；growthθ为等轴枝晶生长取向角度，Random[0，90°]表示在0°～90°之间随机选取一个值；

当等轴枝晶形核核心所处网格的温度T(j，k，t_i)≥(T_L-0.5)，无形核现象发生；

步骤一(6)、通过等轴枝晶形核而发生凝固的网格称为母核心网格(j，k)_moth，母核心网格(j，k)_moth将捕获周围12个邻居网格，如果某一邻居网格对应的state不等于0，则不被捕获；如果某一邻居网格对应的state等于0，则被捕获；

被捕获的网格(j，k)_son将由液态变为生长状态，即state(j，k)_son由0变为1，被捕获网格(j，k)_son具有与母核心网格(j，k)_moth相同的icolor和growthθ值，表明属于同一枝晶；

在12个邻居网格中，如果(j-1，k)_left和(j-2，k)_leftleft均同时被捕获，则这两个网格配对，赋予配对标号1，否则这两个网格之间没有配对关系；如图2；

如果(j+1，k)_right和(j+2，k)_rightright均同时被捕获，则这两个网格配对，赋予配对标号2，否则这两个网格之间没有配对关系；

如果(j-1，k+1)_upleft和(j-1，k+2)_upleftup均同时被捕获，则这两个网格配对，赋予配对标号3，否则这两个网格之间没有配对关系；

如果(j+1，k+1)_upright和(j+1，k+2)_uprightup均同时被捕获，则这两个网格配对，赋予配对标号4，否则这两个网格之间没有配对关系；

如果(j-1，k-1)_downleft和(j-1，k-2)_downleftdown均同时被捕获，则这两个网格配对，赋予配对标号5，否则这两个网格之间没有配对关系；

如果(j+1，k-1)_downright和(j+1，k-2)_downrighdown均同时被捕获，则这两个网格配对，赋予配对标号6，否则这两个网格之间没有配对关系；

当配对的两个网格固相分数都等于1时，这两个网格的状态均变为2且在后续计算中不再改变，一直都是2(代表已经凝固)；

如果配对无法完成，则配对标号为-1，该网格当f_s＝1时即可变为固态state＝2；

作为等轴枝晶形核核心的网格不具有配对网格；

步骤一(7)、某时刻t_i下，如果state(j，k)＝1且f_s(j，k)≠1，执行步骤一(8)；

如果网格state(j，k)＝1且f_s(j，k)＝1，不执行任何计算，在计算域中继续寻找state(j，k)＝1且f_s(j，k)≠1的网格；

步骤一(8)、针对state(j，k)＝1且f_s(j，k)≠1的网格计算固液界面曲率k_curve、固液界面法向角生长过冷度ΔT、生长速度V_tip、几何因子GG、固相分数f_s、液相成分C_l和固相成分C_s；

步骤一(9)、针对state(j，k)＝0和state(j，k)＝1且f_s(j，k)≠1的网格计算溶质扩散方程；

步骤一(10)、计算每个枝晶的比表面积S_v(i)，针对第i个晶粒，计算域中state(j，k)＝1且icolor(j，k)＝i的网格的总数为n-stateone；计算域中state(j，k)＝2且icolor(j，k)＝i的网格的总数为n-statetwo；

建立平均比表面积随固相分数变化曲线，为后续宏观偏析模拟所用。

其它步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是，所述步骤一(3)中确定凝固过程中每个正方形网格的状态；具体过程为：

步骤一(2)中所定义的边界网格只能为液态，即state(j，k)＝0，state为网格的状态；其余网格具有三种状态：

当f_s(j，k)＝0时，网格(j，k)为液态，即state(j，k)＝0；

当0＜f_s(j，k)＜1时，网格(j，k)为生长态，即state(j，k)＝1；

当f_s(j，k)＝1时且配对网格(配对就是后面的1-6的配对方式)的f_s也等于1，则网格(j，k)为固态，即state(j，k)＝2；

f_s(j，k)为网格(j，k)的固相分数(将枝晶生长计算域(已经确定的)进行微观尺度网格剖分的网格(j，k)的固相分数)。

其它步骤及参数与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是，所述步骤一(4)中

其中，T_L为合金液相线温度，单位为℃；冷却速度由实验测量提供。

其它步骤及参数与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是，所述步骤一(8)中针对state(j，k)＝1且f_s(j，k)≠1的网格计算固液界面曲率k_curve、固液界面法向角生长过冷度ΔT、生长速度V_tip、几何因子GG、固相分数f_s、液相成分C_l和固相成分C_s；具体过程为：

捕获state(j，k)＝1且f_s(j，k)≠1的网格的母核心网格为(j_moth，k_moth)；当前state(j，k)＝1的网格的4近邻网格分别为：(j-1，k)、(j+1，k)、(j，k-1)和(j，k+1)；

其中，T_M为纯镁熔点，单位为℃，m_l为液相线斜率，单位为℃/wt％，为平均生长动力学系数，单位为m/s/℃，δ_k为动力学各向异性系数，Γ为吉布斯系数，单位为℃m，为在直角坐标系X轴方向的导数，为在直角坐标系Y轴方向的导数，T为网格温度，由步骤一(4)中给出；

其中，为上一时刻t-Δt的固相分数，为当前时刻t的固相分数，Δt为时间步长，t-Δt为上一时刻；Δstep为边长尺寸；V_tip为生长速度，单位为m/s；

当时，液相成分C_l和固相成分C_s计算如下：

当时，液相成分C_l和固相成分C_s计算如下：

其中，为时刻t的液相成分，为上一时刻t-Δt的液相成分，为时刻t的固相成分，为时刻t-Δt的固相成分；

在t＝0s时，C_s＝0，f_s＝0，C_l＝C_o；

如果state(j，k)＝1且f_s(j，k)≠1的网格所对应的配对标号不等于-1(等于1-6中的某一个数)，则表明state(j，k)＝1且f_s(j，k)≠1的网格和与state(j，k)＝1且f_s(j，k)≠1的网格配对的网格要同时凝固，即状态同时变为2；

时刻t下，配对标号相同的两个网格，当f_s均等于1时，这两个网格均由生长态变为固态；否则，两个配对网格中的一个网格即使f_s＝1也仍将处于生长状态(state＝1)，不能捕捉邻居网格，直至另一个网格的f_s曾加至1，这两个网格均由生长态变为固态。

其它步骤及参数与具体实施方式一至四之一相同。

具体实施方式六：本实施方式与具体实施方式一至五之一不同的是，所述步骤一(9)中针对state(j，k)＝0和state(j，k)＝1且f_s(j，k)≠1的网格计算溶质扩散方程：

其中D_l为溶质元素在液相中的扩散系数，单位为m²/s。

其它步骤及参数与具体实施方式一至五之一相同。

具体实施方式七：本实施方式与具体实施方式一至六之一不同的是，所述步骤一(10)中计算每个枝晶的比表面积S_v(i)，针对第i个晶粒，计算域中state(j，k)＝1且icolor(j，k)＝i的网格的总数为n-stateone；计算域中state(j，k)＝2且icolor(j，k)＝i的网格的总数为n-statetwo；建立平均比表面积随固相分数变化曲线，为后续宏观偏析模拟所用；具体过程为：

计算每个枝晶的比表面积S_v(i)，单位为m²/m³，i表示第i个晶粒；针对第i个晶粒，计算域中state(j，k)＝1且icolor(j，k)＝i的网格的总数为n-stateone；计算域中state(j，k)＝2且icolor(j，k)＝i的网格的总数为n-statetwo；

建立平均比表面积随固相分数变化曲线，为后续宏观偏析模拟所用。

其它步骤及参数与具体实施方式一至六之一相同。

具体实施方式八：本实施方式与具体实施方式一至七之一不同的是，所述步骤二中Δx和Δy的取值范围为：1×10^-3米～4×10^-3米。

其它步骤及参数与具体实施方式一至七之一相同。

具体实施方式九：本实施方式与具体实施方式一至八之一不同的是，所述步骤三中针对所有下角标char＝0的网格(js，ks)_char＝0，计算能量、成分和动量守恒方程，获得铸件内温度、平均成分和速度的分布；具体过程为：

步骤三(1)、计算能量守恒方程，求得铸件内温度T_r分布：

h_s＝c_pLT_r

h_l＝c_pLT_r+ΔH

[H]＝f_srh_s+(1-f_sr)h_l

其中，h_s和h_l分别为固相和液相热焓，c_pL为液相比热，[H]为混合热焓，ρ_L为液体密度，λ_L为液体导热系数，ΔH为结晶潜热，T_r为铸件内温度，f_sr为固相分数(宏观计算得到的固相分数)，为液体流动速度在直角坐标系X和直角坐标系Y方向上的矢量和，在0s时的值为0ms^-1；

步骤三(2)、计算固相分数f_sr：

由夏尔公式：以及步骤三(1)中[H]＝c_pLT_r+(1-f_sr)ΔH的等式关系，得到：

该方程中只有一个未知数T_r，求解该方程可获得铸件内温度T_r，从而得到f_sr；

步骤三(3)、计算成分守恒方程，求得铸件内平均成分[C]分布：

其中

[C]＝f_sr·C_lr·k_par+(1-f_sr)·C_lr

其中，[C]为平均成分，C_lr为液相成分，在0s时的值为0ms^-1；

步骤三(4)、计算动量守恒方程，求得速度：

其中，U_x为直角坐标系X轴方向液体流动速度，单位为m/s；U_y为直角坐标系Y轴方向液体流动速度，单位为m/s；P为压力，单位为Pa；μ_l为液体粘度，单位为Pa.s；g为重力加速度，单位为m/s²，角标x代表X轴，下角标y代表Y轴，β_T和β_C分别为温度膨胀系数(1/℃)和溶质膨胀系数(1/wt％)；K_x和K_y分别为凝固组织在直角坐标系X和直角坐标系Y方向上的渗透率，单位为m²；

平均比表面积随固相分数f_sr变化关系曲线由步骤一得到，采用插值算法获得步骤三(2)计算所得固相分数对应的比表面积；固相分数f_sr、f_s不一样，是为了区分一个是微观计算得到的固相分数，一个是宏观计算得到的固相分数，但是比表面积随固相分数变化的关系曲线是相同的。

其它步骤及参数与具体实施方式一至八之一相同。

具体实施方式十：本实施方式与具体实施方式一至九之一不同的是，所述步骤四中针对所有下角标char非0的计算网格(js，ks)_char≠0，计算能量守恒方程，获得温度T_c分布；T_c为非铸件网格单元的温度，即铸型网格的温度；具体过程为：

[H]＝c_pT_c

其中c_p为比热，ρ_c为密度，λ_c为导热系数。

其它步骤及参数与具体实施方式一至九之一相同。

采用以下实施例验证本发明的有益效果：

实施例一：

选取Mg-4wt％Y合金为研究对象，采用元胞自动机模型模拟α-Mg枝晶生长，模拟计算域尺寸3435μm×3435μm，由500×500个正方形网格组成。数值模拟树脂砂型中放置石墨冷铁铸造过程中Mg-4wt％Y合金宏观偏析形成，二维计算域尺寸254cm×76cm，由127×38个正方形网格组成。计算所需输入参数以及数值模拟砂型铸造过程中Mg-4wt％Y合金宏观偏析形成的相关参数列于表1。

表1 Mg-4wt％Y合金数值模拟热物性参数和计算参数

附图说明：

(1)由附图1可见，正方形网格(j，k)的邻居网格共有12个：(j-1，k)，(j-2，k)，(j+1，k)，(j+2，k)，(j-1，k+1)，(j-1，k+2)，(j+1，k+1)，(j+1，k+2)，(j-1，k-1)，(j-1，k-2)，(j+1，k-1)，(j+1，k-2)。当正方形网格(j，k)的状态变为固态后(state＝2)，如果其有配对网格且配对网格的状态也为固态，就可捕捉周围的邻居网格，从而实现枝晶的生长。

(2)由附图3可见，尽管枝晶生长取向不同，但每个枝晶形貌均呈六重对称状，即一次臂之间的夹角基本接近60°，符合密排六方晶体结构特点。

(3)由附图4可见，铸件在树脂砂型中凝固，铸件端部放置一块石墨冷铁，“TC4”代表测温热电耦安放位置。

(4)由附图5可见，凝固组织由柱状晶区和等轴晶区组成。柱状晶区形成在铸件上部及靠近冷铁处；其余位置为等轴晶区。在等轴晶区切取一块试样(由线框标识)，统计晶粒个数。该试样选取位置靠近热电耦“TC4”的位置。

(5)针对附图5中所选取的位置，采用线切割方法切得试样，对试样表面进行抛光和腐蚀处理，在显微镜下进行金相组织观察，结果示于附图6。在显微镜下可以观察到晶界，用白色线勾勒晶界，从而获得晶粒个数(14个)作为模拟枝晶形貌的输入参数。

(6)实验所获得的温度随时间变化曲线(冷却曲线)示于附图7，温度对时间求一阶导数可得冷却速度，取液相线和固相线之间冷却速度的平均值(约为0.17℃/s)，作为模拟枝晶形貌的输入参数。

(7)基于实验测量所提供的晶粒个数和冷却速度，针对形核位置随机以及生长取向随机的14个晶粒进行数值模拟，模拟所得枝晶结构示于附图8，不同颜色代表不同晶粒，即具有不同生长取向和形核位置的晶粒。多个枝晶凝固过程中，由于溶质场叠加及机械碰撞，枝晶不能得到充分的发展，同时还会发生枝晶臂弯曲现象，这与枝晶生长凝固理论相符合。

(8)枝晶组织模拟过程中，记录不同时刻下固相分数(f_s)和其对应的平均比表面积值，获得随f_s变化曲线，示于附图9。随着固相分数增加，枝晶平均比表面积降低，符合凝固理论中所描述枝晶生长特点。

(9)针对附图4中的铸造系统，进行宏观偏析数值模拟。采用附图9所提供的曲线计算枝晶网络渗透率。随着凝固的进行(固相分数增加)，枝晶网络逐渐发达对液体流动的阻碍也增强，因此减小了渗透率。

附图10b模拟所得铸件中心截面宏观偏析分布与图10a实验测量进行对比，可见两者较好吻合，均呈现如下趋势：正偏析([C]＞C_o)区主要集中在铸件中部和下部，负偏析([C]＜C_o)区主要分布在铸件上部。这是因为凝固过程中析出的Y元素密度高于先析出相α-Mg，液体在靠近左、右铸型处向下流动，在中心部位向上流动，因此Y元素富集在铸件底部。尽管图10c中宏观偏析趋势与图10a相似，但是从偏析程度上看，图10c中正偏析区域明显增大，与实验存在较大差别。

本发明还可有其它多种实施例，在不背离本发明精神及其实质的情况下，本领域技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。

Claims (10)

1.一种镁合金铸造件成分不均匀性数值预测方法，其特征在于：所述方法具体过程为：

步骤一、采用元胞自动机法模拟具有不同生长取向的α-Mg枝晶生长，获得枝晶比表面积随固相分数变化曲线；

步骤二、对铸造系统进行网格剖分，对铸造系统直角坐标系X轴方向和Y轴方向采用相同的网格剖分步长Δx米＝Δy米；

Δx＝X_max-X_min，Δy＝Y_max-Y_min；

计算标号为(js，ks)_char的网格，下角标char＝2表示铸型网格，下角标char＝0表示铸件网格；

其中js和ks均为整数，js为在铸造系统直角坐标系X轴上的坐标，ks为在铸造系统直角坐标系Y轴上的坐标；js的取值范围是1～M，ks的取值范围是1～N；

X_max为铸造系统直角坐标系X轴方向最大值，X_min为铸造系统直角坐标系X轴方向最小值，Y_max为铸造系统直角坐标系Y轴方向最大值，Y_min为铸造系统直角坐标系Y轴方向最小值，单位均为米；

步骤三、针对所有下角标char＝0的网格(js，ks)_char＝0，计算能量、成分和动量守恒方程，获得铸件内温度、平均成分和速度的分布；

步骤四、针对所有下角标char非0的网格(js，ks)_char≠0，计算能量守恒方程，获得温度T_c分布；

T_c为非铸件网格单元的温度，即铸型网格的温度；

步骤五、重复步骤二、步骤三和步骤四，直到凝固结束，即下角标char＝0的计算网格的f_sr(js，ks)_char＝0＝1.0，输出铸件内平均成分场[C]；f_sr为固相分数。

2.根据权利要求1所述一种镁合金铸造件成分不均匀性数值预测方法，其特征在于：所述步骤一中采用元胞自动机法模拟具有不同生长取向的α-Mg枝晶生长，获得枝晶比表面积随固相分数变化曲线；具体过程为：

步骤一(1)、将枝晶生长计算域进行网格剖分，采用边长尺寸为Δstep的正方形网格，每一个正方形网格用(j，k)标识，j的取值范围为[1，n]，k的取值范围为[1，m]，j和k均为整数，n、m均为整数，j为在直角坐标系X轴上的坐标，k为在直角坐标系Y轴上的坐标；

步骤一(2)、为计算域中每一个正方形网格赋予邻居对象，通过对邻居对象的捕捉完成枝晶生长；

标号为(j∈[1，2]，k∈[1，m])、(j∈[n-1，n]，k∈[1，m])、(j∈[3，n-2]，k∈[1，2])、(j∈[3，n-2)]，k∈[m-1，m])的正方形网格定义为边界网格，不参与枝晶生长，因此不需赋予邻居对象；

其余每一个标号为(j，k)的正方形网格都有12个邻居：

(j-1，k)_left，(j-2，k)_leftleft，(j+1，k)_right，(j+2，k)_rightright，(j-1，k+1)_upleft，(j-1，k+2)_upleftup，(j+1，k+1)_upright，(j+1，k+2)_uprightup，(j-1，k-1)_downleft，(j-1，k-2)_downleftdown，(j+1，k-1)_downright，(j+1，k-2)_downrighdown；

步骤一(3)、确定凝固过程中每个正方形网格的状态；

步骤一(4)、给定冷却速度单位为℃/s，计算域中每个网格都具有相同的温度，当凝固时间为t_i时，单位为s，网格(j，k)对应的温度为T(j，k，t_i)；

步骤一(5)、等轴枝晶形核核心数为Neg个，等轴枝晶生长取向角θ随机选择、等轴枝晶形核位置随机选择，等轴枝晶形核过冷度为0.5℃；

某一时刻t_i，当等轴枝晶形核核心所处网格的温度T(j，k，t_i)＜(T_L-0.5)，等轴枝晶形核瞬间发生；等轴枝晶形核核心所处网格对应的物理量为：f_s＝1，C_s＝C_ok_par，C_l＝0，C_aver＝C_sf_s+C₁(1-f_s)，state＝2，icolor＝Random[1，Neg]，growthθ＝Random[0，90°]；

其中，f_s为固相分数，最小值为0表示液相，最大值为1表示固相；C_s为固相成分；C_o为合金初始成分；k_par为平衡分配系数；C_l为液相成分；C_aver为平均成分，成分的单位为wt％；T_L为合金液相线温度，单位为℃；

每一个等轴枝晶形核核心赋予一个标号icolor，用来区分不同晶粒，Random[1，Neg]表示在1～Neg之间选取一个随机数；

state为状态；growthθ为等轴枝晶生长取向角度，Random[0，90°]表示在0°～90°之间随机选取一个值；

当等轴枝晶形核核心所处网格的温度T(j，k，t_i)≥(T_L-0.5)，无形核现象发生；

步骤一(6)、通过等轴枝晶形核而发生凝固的网格称为母核心网格(j，k)_moth，母核心网格(j，k)_moth将捕获周围12个邻居网格，如果某一邻居网格对应的state不等于0，则不被捕获；如果某一邻居网格对应的state等于0，则被捕获；

被捕获的网格(j，k)_son将由液态变为生长状态，即state(j，k)_son由0变为1，被捕获网格(j，k)_son具有与母核心网格(j，k)_moth相同的icolor和growthθ值，表明属于同一枝晶；

在12个邻居网格中，如果(j-1，k)_left和(j-2，k)_leftleft均同时被捕获，则这两个网格配对，赋予配对标号1，否则这两个网格之间没有配对关系；

如果(j+1，k)_right和(j+2，k)_rightright均同时被捕获，则这两个网格配对，赋予配对标号2，否则这两个网格之间没有配对关系；

如果(j-1，k+1)_upleft和(j-1，k+2)_upleftup均同时被捕获，则这两个网格配对，赋予配对标号3，否则这两个网格之间没有配对关系；

如果(j+1，k+1)_upright和(j+1，k+2)_uprightup均同时被捕获，则这两个网格配对，赋予配对标号4，否则这两个网格之间没有配对关系；

如果(j-1，k-1)_downleft和(j-1，k-2)_downleftdown均同时被捕获，则这两个网格配对，赋予配对标号5，否则这两个网格之间没有配对关系；

如果(j+1，k-1)_downright和(j+1，k-2)_downrighdown均同时被捕获，则这两个网格配对，赋予配对标号6，否则这两个网格之间没有配对关系；

当配对的两个网格固相分数都等于1时，这两个网格的状态均变为2且在后续计算中不再改变；

如果配对无法完成，则配对标号为-1，该网格当f_s＝1时即可变为固态state＝2；

作为等轴枝晶形核核心的网格不具有配对网格；

步骤一(7)、某时刻t_i下，如果state(j，k)＝1且f_s(j，k)≠1，执行步骤一(8)；

如果网格state(j，k)＝1且f_s(j，k)＝1，不执行任何计算，在计算域中继续寻找state(j，k)＝1且f_s(j，k)≠1的网格；

步骤一(8)、针对state(j，k)＝1且f_s(j，k)≠1的网格计算固液界面曲率k_curve、固液界面法向角生长过冷度ΔT、生长速度V_tip、几何因子GG、固相分数f_s、液相成分C_l和固相成分C_s；

步骤一(9)、针对state(j，k)＝0和state(j，k)＝1且f_s(j，k)≠1的网格计算溶质扩散方程；

步骤一(10)、计算每个枝晶的比表面积S_v(i)，针对第i个晶粒，计算域中state(j，k)＝1且icolor(j，k)＝i的网格的总数为n-stateone；计算域中state(j，k)＝2且icolor(j，k)＝i的网格的总数为n-statetwo；

建立平均比表面积随固相分数变化曲线。

3.根据权利要求2所述一种镁合金铸造件成分不均匀性数值预测方法，其特征在于：所述步骤一(3)中确定凝固过程中每个正方形网格的状态；具体过程为：

步骤一(2)中所定义的边界网格只能为液态，即state(j，k)＝0，state为网格的状态；其余网格具有三种状态：

当f_s(j，k)＝0时，网格(j，k)为液态，即state(j，k)＝0；

当0＜f_s(j，k)＜1时，网格(j，k)为生长态，即state(j，k)＝1；

当f_s(j，k)＝1时且配对网格的f_s也等于1，则网格(j，k)为固态，即state(j，k)＝2；

f_s(j，k)为网格(j，k)的固相分数。

4.根据权利要求3所述一种镁合金铸造件成分不均匀性数值预测方法，其特征在于：所述步骤一(4)中

其中，T_L为合金液相线温度，单位为℃。

5.根据权利要求4所述一种镁合金铸造件成分不均匀性数值预测方法，其特征在于：所述步骤一(8)中针对state(j，k)＝1且f_s(j，k)≠1的网格计算固液界面曲率k_curve、固液界面法向角生长过冷度ΔT、生长速度V_tip、几何因子GG、固相分数f_s、液相成分C_l和固相成分C_s；具体过程为：

捕获state(j，k)＝1且f_s(j，k)≠1的网格的母核心网格为(j_moth，k_moth)；当前state(j，k)＝1的网格的4近邻网格分别为：(j-1，k)、(j+1，k)、(j，k-1)和(j，k+1)；

其中，T_M为纯镁熔点，单位为℃，m₁为液相线斜率，单位为℃/wt％，为平均生长动力学系数，单位为m/s/℃，δ_k为动力学各向异性系数，Γ为吉布斯系数，单位为℃m，为在直角坐标系X轴方向的导数，为在直角坐标系Y轴方向的导数，T为网格温度；

其中，为上一时刻t-Δt的固相分数，为当前时刻t的固相分数，Δt为时间步长，t-Δt为上一时刻；Δstep为边长尺寸；V_tip为生长速度，单位为m/s；

当时，液相成分C_l和固相成分C_s计算如下：

当时，液相成分C_l和固相成分C_s计算如下：

其中，为时刻t的液相成分，为上一时刻t-Δt的液相成分，为时刻t的固相成分，为时刻t-Δt的固相成分；

在t＝0s时，C_s＝0，f_s＝0，C_l＝C_o；

如果state(j，k)＝1且f_s(j，k)≠1的网格所对应的配对标号不等于-1，则表明state(j，k)＝1且f_s(j，k)≠1的网格和与state(j，k)＝1且f_s(j，k)≠1的网格配对的网格要同时凝固，即状态同时变为2；

时刻t下，配对标号相同的两个网格，当f_s均等于1时，这两个网格均由生长态变为固态；否则，两个配对网格中的一个网格即使f_s＝1也仍将处于生长状态state＝1，不能捕捉邻居网格，直至另一个网格的f_s曾加至1，这两个网格均由生长态变为固态。

6.根据权利要求5所述一种镁合金铸造件成分不均匀性数值预测方法，其特征在于：所述步骤一(9)中针对state(j，k)＝0和state(j，k)＝1且f_s(j，k)≠1的网格计算溶质扩散方程：

其中D_l为溶质元素在液相中的扩散系数，单位为m²/s。

7.根据权利要求6所述一种镁合金铸造件成分不均匀性数值预测方法，其特征在于：所述步骤一(10)中计算每个枝晶的比表面积S_v(i)，针对第i个晶粒，计算域中state(j，k)＝1且icolor(j，k)＝i的网格的总数为n-stateone；计算域中state(j，k)＝2且icolor(j，k)＝i的网格的总数为n-statetwo；建立平均比表面积随固相分数变化曲线；具体过程为：

计算每个枝晶的比表面积S_v(i)，单位为m²/m³，i表示第i个晶粒；

针对第i个晶粒，计算域中state(j，k)＝1且icolor(j，k)＝i的网格的总数为n-stateone；计算域中state(j，k)＝2且icolor(j，k)＝i的网格的总数为n-statetwo；

建立平均比表面积随固相分数变化曲线。

8.根据权利要求7所述一种镁合金铸造件成分不均匀性数值预测方法，其特征在于：所述步骤二中Δx和Δy的取值范围为：1×10^-3米～4×10^-3米。

9.根据权利要求8所述一种镁合金铸造件成分不均匀性数值预测方法，其特征在于：所述步骤三中针对所有下角标char＝0的网格(js，ks)_char＝0，计算能量、成分和动量守恒方程，获得铸件内温度、平均成分和速度的分布；具体过程为：

步骤三(1)、计算能量守恒方程，求得铸件内温度T_r分布：

h_s＝c_pLT_r

h_l＝c_pLT_r+ΔH

[H]＝f_srh_s+(1-f_sr)h_l

其中，h_s和h_l分别为固相和液相热焓，c_pL为液相比热，[H]为混合热焓，ρ_L为液体密度，λ_L为液体导热系数，ΔH为结晶潜热，T_r为铸件内温度，f_sr为固相分数，为液体流动速度在直角坐标系X和直角坐标系Y方向上的矢量和，在0s时的值为0m s^-1；

步骤三(2)、计算固相分数f_sr：

由夏尔公式：以及步骤三(1)中[H]＝c_pLT_r+(1-f_sr)ΔH的等式关系，得到：

求解该方程获得铸件内温度T_r，从而得到f_sr；

步骤三(3)、计算成分守恒方程，求得铸件内平均成分[C]分布：

其中

[C]＝f_sr·C_lr·k_par+(1-f_sr)·C_lr

其中，[C]为平均成分，C_lr为液相成分，在0s时的值为0m s^-1；

步骤三(4)、计算动量守恒方程，求得速度：

其中，U_x为直角坐标系X轴方向液体流动速度，单位为m/s；U_y为直角坐标系Y轴方向液体流动速度，单位为m/s；P为压力，单位为Pa；μ_l为液体粘度，单位为Pa.s；g为重力加速度，单位为m/s²，角标x代表X轴，下角标y代表Y轴，β_T和β_C分别为温度膨胀系数和溶质膨胀系数；K_x和K_y分别为凝固组织在直角坐标系X和直角坐标系Y方向上的渗透率，单位为m²；

平均比表面积随固相分数f_sr变化关系曲线由步骤一得到，采用插值算法获得步骤三(2)计算所得固相分数对应的比表面积。

10.根据权利要求9所述一种镁合金铸造件成分不均匀性数值预测方法，其特征在于：所述步骤四中针对所有下角标char非0的计算网格(js，ks)_char≠0，计算能量守恒方程，获得温度T_c分布；T_c为非铸件网格单元的温度，即铸型网格的温度；具体过程为：

[H]＝c_pT_c

其中c_p为比热，ρ_c为密度，λ_c为导热系数。