CN116230142B - 一种铝合金凝固动力学过程的介观尺度预测方法 - Google Patents

Abstract

本发明公开了一种铝合金凝固动力学过程的介观尺度预测方法，包括以下步骤：S1、初始化网络；S2、基于连续形核模型生成形核点；S3、计算流场；S4、计算溶质场；S5、计算温度场；S6、基于修正的偏心正方形算法捕获规则，将液相单元转变为界面单元；S7、界面单元生长；S8、更新元胞状态；S9、计算糊状区的渗透率；S10、判断是否满足结束条件，若不满足，进行步骤S3，若满足，输出模拟结果。本发明采用上述一种铝合金凝固动力学过程的介观尺度预测方法，能够在节省计算成本的情况下更加准确有效地预测铝合金的凝固过程，同时提高计算效率和数值稳定性。

Description

一种铝合金凝固动力学过程的介观尺度预测方法

技术领域

本发明涉及铝合金铸造技术领域，尤其是涉及一种铝合金凝固动力学过程的介观尺度预测方法。

背景技术

铸造铝合金因其具有高强度重量比、良好的延展性和耐腐蚀性、低成本以及优异的铸造性能而成为各行各业的重要合金材料。金属合金的宏观性能主要由其微观组织结构决定，提高材料性能的各种方法本质上都是从控制凝固组织结构入手，凝固组织的形貌、尺寸和分布影响铝合金的最终性能，研究和控制金属材料的凝固组织已成为材料科学领域的热点问题之一。而凝固作为一个跨尺度现象，宏观上涉及热量、质量传输，介观上表现为枝晶与晶粒长大，微观上表现为形核生长，枝晶的形核与生长决定了微观组织形貌。随着人们对铸件的可靠性等要求越来越高，同时对合金综合性能和特种性能的要求也不断提高，所以铝合金铸造以及铝合金的结晶行为和凝固组织的预测显得尤为重要。

由于铝合金凝固过程的复杂性，实验参数不易精确控制，实验条件易受干扰，一些学者通过数值模拟的方式再现凝固过程。但研究大多采用简化的热质耦合或流质耦合，由于铝合金的凝固过程是一个热质流多物理场交互作用的复杂演变过程，一些学者通过元胞自动机耦合有限元法或者格子玻尔兹曼方法研究热质流耦合，现有的模型不能反映出晶粒的真实形貌，二次、三次枝晶臂不突出；另外一些学者通过相场耦合格子玻尔兹曼的方式研究热质流耦合，这种模型计算量大，计算效率低，只能模拟较小区域内的枝晶生长。

发明内容

本发明的目的是提供一种铝合金凝固动力学过程的介观尺度预测方法，解决上述背景技术中提出的问题。

为实现上述目的，本发明提供了一种铝合金凝固动力学过程的介观尺度预测方法，包括以下步骤：

S1、初始化网络；

S2、基于连续形核模型生成形核点；

S3、计算流场；

S4、计算溶质场；

S5、计算温度场；

S6、基于修正的偏心正方形算法捕获规则，将液相单元转变为界面单元；

S7、界面单元生长；

S8、更新元胞状态；

S9、计算糊状区的渗透率；

S10、判断是否满足结束条件，若不满足，进行步骤S3，若满足，输出模拟结果，凝固结束。

优选的，步骤S1中，将固定的二维计算区域均匀划分为正方形网格，每个正方形网格代表一个单元胞，初始化单元胞，为每个单元胞的状态参数、物性参数、固相率、液相溶质浓度、固相溶质浓度、速度、分布函数赋予初始值；其中单元胞的状态包括三个，即液相胞、界面胞和固相胞，分别对应三个状态参数值State(i,j)＝0、State(i,j)＝1和State(i,j)＝2。

优选的，步骤S2中，过冷度ΔT通过晶粒形核密度n(ΔT)对应一个一定的晶粒形核概率，通过基于高斯分布的连续形核模型生成形核点；

形核发生在一系列的形核位置上，形核质点为连续而非离散的分布函数，晶粒形核的模型实质为过冷度ΔT关于晶粒形核密度n(ΔT)的函数，过冷度ΔT与晶粒形核密度n(ΔT)一一对应，晶粒形核密度n(ΔT)代表着晶粒形核的概率大小是一定的，所以过冷度ΔT通过晶粒形核密度n(ΔT)对应一个一定的晶粒形核的概率；新生的晶核密度可表示为：

其中晶粒形核密度随过冷度的关系式为：

式中：n_max为最大晶粒密度，ΔT为元胞的过冷度，ΔT_max为最大形核过冷度，ΔT_σ为过冷度标准差，这两个式子的组合即为晶粒的形核模型。

优选的，步骤S3中，基于格子玻尔兹曼模型将铝合金液体流动过程分解成碰撞和迁移两个部分，基于D2Q9格子玻尔兹曼算法计算流场的控制方程；根据不同的对流条件对流场的控制方程在枝晶边界处应用不同的边界条件，得到枝晶边界处节点的动量分布函数值；对节点的动量分布函数值进行迁移，并在迁移之后施加边界条件，得到各节点的流体速度；

利用D2Q9格子玻尔兹曼算法计算流场的控制方程为铝合金液体流动过程中的节点碰撞方程：

式中，r是位置矢量，t表示时间，dt是时间步长，α＝0,1,...,8表示9个离散速度方向，其中e_α为离散速度分布式：

f_α ^eq表示速度的平衡态分布函数，如下式所示：

式中ρ(x,t)代表t时刻LBM模型中格子的密度，dx为空间步长，u为节点的流体速度，w_α是权系数，如下式所示：

无量纲松弛时间τ_f与流体动力黏度之间的关系为：

F_i为计算流场产生的浮力所添加的源项，当对流条件为强迫对流时，F_i＝0；当对流条件为自然对流时，F_i的计算式如下：

式中，F为自然对流引起的浮力，由温度差和浓度差决定，计算式如下：

F＝gρ₀β_T(T-T₀)+gρ₀β_C(C-C₀)

式中，g是重力加速度，ρ₀是流体的密度，β_T是温度膨胀系数，β_C是浓度膨胀系数，T代表流体的温度，T₀代表流体的初始温度，C代表流体的浓度，C₀代表流体的初始浓度；

当对流条件为强迫对流时，左右边界采用充分发展边界，并给左边界一个恒定的入流速度，为了保证计算区域的质量和动量守恒，网格的上下边界采用周期性边界：

f_α(x,0,t+dt)＝f'_α(x,N_y,t)

f_α(x,N_y,t+dt)＝f'_α(x,0,t)

式中，N_y表示y方向的网格数；

强迫对流下对应宏观的密度和速度可由下式确定：

当对流条件为自然对流时，网格的边界条件采用非平衡外推格式：

f_α(r+e_α,t+dt)＝f_α ^(eq)(r+e_α,t+dt)-f_α ^(ne)(r+e_α,t+dt)

＝f_α ^(eq)(r+e_α,t+dt)+[f_α(r,t)-f_α ^eq(r,t)]

其中，f_α ^(eq)是平衡态部分，f_α ^(ne)是非平衡态部分；

对应宏观的密度和速度可由下式确定：

优选的，步骤S4中，把流场f_α替换成溶质场g_α，将流场的源项F_i替换为溶质场的源项w_αC_l(1-k)Δf_s，就能得到带源项的计算溶质场的标准格子玻尔兹曼方程，从而解决溶质输运的问题，在枝晶边界处采用非平衡外推的边界条件，得到枝晶边界处节点的浓度分布函数值；对节点的浓度分布函数值进行迁移，并在迁移之后施加边界条件，得到各节点的浓度分布；

利用D2Q9格子玻尔兹曼算法建立溶质场模型，溶质场演化方程为：

式中，表示流体浓度的平衡态分布函数，如下式所示：

式中，C(x,t)代表t时刻LBM模型中格子的浓度；

无量纲松弛时间τ_g与溶质扩散系数D之间的关系为：

网格边界条件采用非平衡外推格式，对应宏观的浓度可由下式确定：

源项中，Δf_s是界面元胞在一个时间步长内的固相增量，与步骤S7有关，k为平衡分配系数，C_l为液相溶质浓度，与步骤S7有关，步骤S3得到的流体密度和速度也影响步骤S4中的浓度。

优选的，步骤S5中，把流场f_α替换成温度场h_α，将流场的源项F_i替换为温度场的源项w_αΔf_sL/c_p，就能得到计算温度场的标准格子玻尔兹曼方程，从而解决热量输运的问题，在枝晶边界处采用非平衡外推的边界条件，得到枝晶边界处节点的温度分布函数值；对节点的温度分布函数值进行迁移，并在迁移之后施加边界条件，得到各节点的流体温度；

利用D2Q9格子玻尔兹曼算法建立温度场模型，温度场演化方程为：

式中，表示速度的平衡态分布函数，如下式所示：

式中，T(x,t)代表t时刻LBM模型中格子的温度；

无量纲松弛时间τ_h与热扩散系数λ之间的关系为：

网格边界条件采用非平衡外推格式，对应宏观的温度可由下式确定：

源项中，L为潜热，c_p为热容，步骤S3得到的流体密度和速度也影响步骤S5中的温度。

优选的，步骤S6中，基于修正的偏心正方形算法捕获规则，采用一系列对角线方向与枝晶生长方向相同的正方形，辅助元胞进行在生长方向上的界面胞捕获和状态传递，将相邻的状态值为State(i,j)＝0液相单元转变为状态值为State(i,j)＝1的界面单元，利用此算法可获得任意生长角度的枝晶；

当计算域中某元胞形核后，其状态参数值由State(i,j)＝0变为状态参数值State(i,j)＝1，可指定一个也可随机赋予一个优先生长角度θ₀，形核后的(i,j)单元受到枝晶生长动力学的驱动，固相分数不断增加，在形核元胞的中心根据元胞的固相分数布置一个有倾斜角度θ₀的正方形，其对角线长度随生长发生偏移，该正方形的半对角线长可由下式计算：

L₁(t+dt)＝L₁(t)+V_n·dt

式中，V_n为固液界面法向生长速度，n为界面法向，

n＝n_xi+n_yj

式中，n_x,n_y为界面法向n的分量，计算式如下：

式中，和/>分别表示固相率在x轴和y轴上的一阶偏导，/>和/>分别为固相率在x轴和y轴上的二阶偏导，其中f_s可通过步骤S7计算；

假设在界面胞中固液界面为平界面，V_n的计算式如下：

V_n＝V_xn_x+V_yn_y

式中，V_x,V_y为固液界面在x,y方向的速度分量，计算式如下：

式中，k为平衡分配系数，D_s为固相扩散系数，D_l为液相扩散系数，为界面(i,j)处液相浓度，/>为界面(i,j)处固相浓度，C_s(i-1,j)为节点(i-1,j)处固相浓度，f_s(i-1,j)为节点(i-1,j)处固相分数，C_s(i+1,j)为节点(i+1,j)处固相浓度，f_s(i+1,j)为节点(i+1,j)处固相分数，C_l(i-1,j)为节点(i-1,j)处液相浓度，C_l(i+1,j)为节点(i+1,j)处液相浓度，C_s(i,j-1)为节点(i,j-1)处固相浓度，f_s(i,j-1)为节点(i,j-1)处固相分数，C_s(i,j+1)为节点(i,j+1)处固相浓度，f_s(i,j+1)为节点(i,j+1)处固相分数，C_l(i,j-1)为节点(i,j-1)处液相浓度，C_l(i,j+1)为节点(i,j+1)处液相浓度；

随着(i,j)单元固相分数的增加，偏移的正方形的四个顶点触碰到周围的液相单元，把相邻的液相单元(i,j+1)，(i,j-1)，(i-1,j)，(i+1,j)的状态值变为界面单元，其状态参数值变为State(i,j)＝1，并把优先生长角度θ₀赋予给这四个被捕获的单元，保持生长方向与初始生长方向的一致性，被捕获的单元又以同样的方式捕获周围单元中的液相元胞，其顶点形成一个新的次级正方形，并继续生长，为避免元胞的八邻居元胞出现未被捕获的情况，当某元胞的固相率为1时，将八邻居元胞中的液相胞全部捕获为界面胞；通过对正方形顶点的元胞进行标记，不断计算其位置，正方形的分布受界面生长动力学的显著影响，利用此修正的偏心正方形算法可良好地反映出等轴晶的形貌，且枝干密集。

优选的，步骤S7中，对于状态值为State(i,j)＝1的界面单元，采用元胞自动机模型建立晶粒的生长模型，计算界面曲率、界面各向异性、溶质浓度场、固相分数的增加量，当许多元胞状态发生改变宏观上就变现为晶粒的增长，溶质的扩散过程是伴随晶粒生长产生的；

采用平均计数法求解界面曲率：

式中，dx为空间步长，N为当前元胞的邻胞数，对于此二维模型计算域，N＝9，为第k₁个邻胞的固相分数；

铝合金凝固过程中的枝晶为具有四重对称轴的立方晶体，界面各向异性函数的计算如下：

f(θ,θ₀)＝1-15εcos[4(θ-θ₀)]

式中，ε为界面能各向异性强度，θ为固液界面法向向量与水平方向的夹角，θ₀为枝晶最优生长角；

采用平均浓度C表征枝晶生长过程中的溶质分布：

C＝(1-f_s)C_l+f_sC_s

式中，C_s与C_l分别为固、液相溶质浓度，界面胞的C为固、液相浓度按固相率的插值，固液界面的溶质浓度为：

式中，与/>分别为界面固相溶质浓度和界面液相溶质浓度，k为溶质平衡分配系数；固液界面单元界面局部平衡认为界面液相溶质浓度与界面实际溶质成分的溶质差是枝晶生长的驱动力，界面液相溶质浓度/>的计算式如下：

式中，T^*为界面温度，T_l为平衡液相线温度，m_l为液相线的斜率，Γ为Gibbs-Thomson系数；

固相分数的计算式如下，固相分数增加的过程即为晶粒的生长过程：

式中，L₁(t+dt)通过步骤S6计算，θ为界面法向与水平方向的夹角，计算式如下：

优选的，步骤S8中，当状态值为State(i,j)＝1的界面元胞的固相分数f_s＝1时，将其状态值修改为State(i,j)＝2，并将界面元胞的液相浓度修改为0；

步骤S9中，结合Darcy定律，计算铝合金液相在糊状区的渗透率；

Darcy定律认为渗透率与一次枝晶臂间距有关，公式如下：

式中，λ₁为一次枝晶臂间距；

结合Darcy定律，推导得出格子单位下多孔介质渗透率的计算式为：

式中，表示流体动力粘度，Δρ表示入口和出口的密度差，u_x表示x方向的速度，N_x,N_y表示x、y方向的网格数。

因此，本发明采用上述一种铝合金凝固动力学过程的介观尺度预测方法，有益效果在于：本发明的模型能够在节省计算成本的情况下更真实有效地再现铝合金凝固过程中等轴晶、柱状晶形核长大过程中的流场、浓度场、温度场分布，并能在不同流场存在的情况下对等轴晶、柱状晶的具体尺寸、二次枝晶臂间距、枝晶臂间凹槽处的凝固生长进行定量分析，且能计算出铝合金液相在糊状区的渗透率，能够更加准确有效地预测铝合金的凝固过程，同时提高了计算效率和数值稳定性。

下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。

附图说明

图1为本发明的流程示意图；

图2为本发明修正的偏心正方形算法的示意图；

图3为本发明算法的类视图；

图4为本发明无流场时枝晶最优生长角θ₀为45度时枝晶的浓度场分布图；

图5为本发明流场条件为强迫对流时枝晶最优生长角θ₀为45度时枝晶的浓度场分布图；

图6为本发明流场条件为自然对流时枝晶最优生长角θ₀为45度时枝晶的浓度场分布图；

图7为本发明不同取向枝晶在强迫对流条件下形貌的示意图；

图8为本发明择优取向角θ₀为45度时在强迫对流条件下不同入流速度下的形貌的示意图，；

图9为本发明择优取向角θ₀为45度时在强迫对流条件下不同冷却速率下的形貌的示意图；

图10为本发明择优取向角θ₀为45度时在强迫对流条件下不同过冷度下的形貌的示意图；

图11为本发明择优取向角θ₀为0度时在自然对流条件下形貌的示意图；

图12为本发明随机多个枝晶在强迫对流条件下不同时刻下的形貌的示意图；

图13为本发明随机多个枝晶在自然对流条件下不同时刻下的形貌的示意图；

图14为本发明随机多个柱状晶在强迫对流条件下不同时刻下的形貌的示意图；

图15为本发明随机多个柱状晶在自然对流条件下不同时刻下的形貌的示意图。

具体实施方式

以下通过附图和实施例对本发明的技术方案作进一步说明。

除非另外定义，本发明使用的技术术语或者科学术语应当为本发明所属领域内具有一般技能的人士所理解的通常意义。本发明中使用的“第一”、“第二”以及类似的词语并不表示任何顺序、数量或者重要性，而只是用来区分不同的组成部分。“包括”或者“包含”等类似的词语意指出现该词前面的元件或者物件涵盖出现在该词后面列举的元件或者物件及其等同，而不排除其他元件或者物件。术语“设置”、“安装”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。“上”、“下”、“左”、“右”等仅用于表示相对位置关系，当被描述对象的绝对位置改变后，则该相对位置关系也可能相应地改变。

实施例

如图1所示，本发明所述的一种铝合金凝固动力学过程的介观尺度预测方法，包括以下步骤：

步骤S1：初始化网络。

将固定的二维计算区域均匀划分为正方形网格，每个正方形网格代表一个单元胞，初始化单元胞，为每个单元胞的状态参数、物性参数、固相率、液相溶质浓度、固相溶质浓度、速度、分布函数赋予初始值。其中单元胞的状态包括三个，即液相胞、界面胞和固相胞，分别对应三个状态参数值State(i,j)＝0、State(i,j)＝1和State(i,j)＝2。

步骤S2：基于连续形核模型生成形核点。

过冷度ΔT通过晶粒形核密度n(ΔT)对应一个一定的晶粒形核概率，通过基于高斯分布的连续形核模型生成形核点。

形核发生在一系列的形核位置上，形核质点为连续而非离散的分布函数，晶粒形核的模型实质为过冷度ΔT关于晶粒形核密度n(ΔT)的函数，过冷度ΔT与晶粒形核密度n(ΔT)一一对应，晶粒形核密度n(ΔT)代表着晶粒形核的概率大小是一定的，所以过冷度ΔT通过晶粒形核密度n(ΔT)对应一个一定的晶粒形核的概率。新生的晶核密度可表示为：

其中晶粒形核密度随过冷度的关系式为：

式中：n_max为最大晶粒密度，ΔT为元胞的过冷度，ΔT_max为最大形核过冷度，ΔT_σ为过冷度标准差，这两个式子的组合即为晶粒的形核模型。

步骤S3：计算流场。

基于格子玻尔兹曼模型将铝合金液体流动过程分解成碰撞和迁移两个部分，基于D2Q9格子玻尔兹曼算法计算流场的控制方程；根据不同的对流条件对流场的控制方程在枝晶边界处应用不同的边界条件，得到枝晶边界处节点的动量分布函数值；对节点的动量分布函数值进行迁移，并在迁移之后施加边界条件，得到各节点的流体速度。

利用D2Q9格子玻尔兹曼算法计算流场的控制方程为铝合金液体流动过程中的节点碰撞方程：

式中，r是位置矢量，t表示时间，dt是时间步长，α＝0,1,...,8表示9个离散速度方向，其中e_α为离散速度分布式：

f_α ^eq表示速度的平衡态分布函数，如下式所示：

式中ρ(x,t)代表t时刻LBM模型中格子的密度，dx为空间步长，u为节点的流体速度，w_α是权系数，如下式所示：

无量纲松弛时间τ_f与流体动力黏度之间的关系为：

F_i为计算流场产生的浮力所添加的源项，当对流条件为强迫对流时，F_i＝0；当对流条件为自然对流时，F_i的计算式如下：

式中，F为自然对流引起的浮力，由温度差和浓度差决定，计算式如下：

F＝gρ₀β_T(T-T₀)+gρ₀β_C(C-C₀)

式中，g是重力加速度，ρ₀是流体的密度，β_T是温度膨胀系数，β_C是浓度膨胀系数，T代表流体的温度，T₀代表流体的初始温度，C代表流体的浓度，C₀代表流体的初始浓度。

当对流条件为强迫对流时，左右边界采用充分发展边界，并给左边界一个恒定的入流速度，为了保证计算区域的质量和动量守恒，网格的上下边界采用周期性边界：

f_α(x,0,t+dt)＝f'_α(x,N_y,t)

f_α(x,N_y,t+dt)＝f'_α(x,0,t)

式中，N_y表示y方向的网格数。

强迫对流下对应宏观的密度和速度可由下式确定：

当对流条件为自然对流时，网格的边界条件采用非平衡外推格式：

f_α(r+e_α,t+dt)＝f_α ^(eq)(r+e_α,t+dt)-f_α ^(ne)(r+e_α,t+dt)

＝f_α ^(eq)(r+e_α,t+dt)+[f_α(r,t)-f_α ^eq(r,t)]

其中，f_α ^(eq)是平衡态部分，f_α ^(ne)是非平衡态部分。

对应宏观的密度和速度可由下式确定：

步骤S4：计算溶质场。

把流场f_α替换成溶质场g_α，将流场的源项F_i替换为溶质场的源项w_αC_l(1-k)Δf_s，就能得到带源项的计算溶质场的标准格子玻尔兹曼方程，从而解决溶质输运的问题，在枝晶边界处采用非平衡外推的边界条件，得到枝晶边界处节点的浓度分布函数值。对节点的浓度分布函数值进行迁移，并在迁移之后施加边界条件，得到各节点的浓度分布。

利用D2Q9格子玻尔兹曼算法建立溶质场模型，溶质场演化方程为：

式中，表示流体浓度的平衡态分布函数，如下式所示：

式中，C(x,t)代表t时刻LBM模型中格子的浓度。

无量纲松弛时间τ_g与溶质扩散系数D之间的关系为：

网格边界条件采用非平衡外推格式，对应宏观的浓度可由下式确定：

源项中，Δf_s是界面元胞在一个时间步长内的固相增量，与步骤S7有关，k为平衡分配系数，C_l为液相溶质浓度，与步骤S7有关，步骤S3得到的流体密度和速度也影响步骤S4中的浓度。

步骤S5：计算温度场。

把流场f_α替换成温度场h_α，将流场的源项F_i替换为温度场的源项w_αΔf_sL/c_p，就能得到计算温度场的标准格子玻尔兹曼方程，从而解决热量输运的问题，在枝晶边界处采用非平衡外推的边界条件，得到枝晶边界处节点的温度分布函数值。对节点的温度分布函数值进行迁移，并在迁移之后施加边界条件，得到各节点的流体温度。

利用D2Q9格子玻尔兹曼算法建立温度场模型，温度场演化方程为：

式中，表示速度的平衡态分布函数，如下式所示：

式中，T(x,t)代表t时刻LBM模型中格子的温度；

无量纲松弛时间τ_h与热扩散系数λ之间的关系为：

网格边界条件采用非平衡外推格式，对应宏观的温度可由下式确定：

源项中，L为潜热，c_p为热容，步骤S3得到的流体密度和速度也影响步骤S5中的温度。

步骤S6：基于修正的偏心正方形算法捕获规则，将液相单元转变为界面单元。

基于修正的偏心正方形算法捕获规则，采用一系列对角线方向与枝晶生长方向相同的正方形，辅助元胞进行在生长方向上的界面胞捕获和状态传递，将相邻的状态值为State(i,j)＝0液相单元转变为状态值为State(i,j)＝1的界面单元，利用此算法可获得任意生长角度的枝晶。

当计算域中某元胞形核后，其状态参数值由State(i,j)＝0变为状态参数值State(i,j)＝1，可指定一个也可随机赋予一个优先生长角度θ₀，形核后的(i,j)单元受到枝晶生长动力学的驱动，固相分数不断增加，在形核元胞的中心根据元胞的固相分数布置一个有倾斜角度θ₀的正方形，其对角线长度随生长发生偏移，该正方形的半对角线长可由下式计算：

L₁(t+dt)＝L₁(t)+V_n·dt

式中，V_n为固液界面法向生长速度，n为界面法向，

n＝n_xi+n_yj

式中，n_x,n_y为界面法向n的分量，计算式如下：

式中，和/>分别表示固相率在x轴和y轴上的一阶偏导，/>和/>分别为固相率在x轴和y轴上的二阶偏导，其中f_s可通过步骤S7计算。

假设在界面胞中固液界面为平界面，V_n的计算式如下：

V_n＝V_xn_x+V_yn_y

式中，V_x,V_y为固液界面在x,y方向的速度分量，计算式如下：

式中，k为平衡分配系数，D_s为固相扩散系数，D_l为液相扩散系数，为界面(i,j)处液相浓度，/>为界面(i,j)处固相浓度，C_s(i-1,j)为节点(i-1,j)处固相浓度，f_s(i-1,j)为节点(i-1,j)处固相分数，C_s(i+1,j)为节点(i+1,j)处固相浓度，f_s(i+1,j)为节点(i+1,j)处固相分数，C_l(i-1,j)为节点(i-1,j)处液相浓度，C_l(i+1,j)为节点(i+1,j)处液相浓度，C_s(i,j-1)为节点(i,j-1)处固相浓度，f_s(i,j-1)为节点(i,j-1)处固相分数，C_s(i,j+1)为节点(i,j+1)处固相浓度，f_s(i,j+1)为节点(i,j+1)处固相分数，C_l(i,j-1)为节点(i,j-1)处液相浓度，C_l(i,j+1)为节点(i,j+1)处液相浓度。

随着(i,j)单元固相分数的增加，偏移的正方形的四个顶点触碰到周围的液相单元，把相邻的液相单元(i,j+1)，(i,j-1)，(i-1,j)，(i+1,j)的状态值变为界面单元，其状态参数值变为State(i,j)＝1，并把优先生长角度θ₀赋予给这四个被捕获的单元，保持生长方向与初始生长方向的一致性，被捕获的单元又以同样的方式捕获周围单元中的液相元胞，其顶点形成一个新的次级正方形，并继续生长，为避免元胞的八邻居元胞出现未被捕获的情况，当某元胞的固相率为1时，将八邻居元胞中的液相胞全部捕获为界面胞。通过对正方形顶点的元胞进行标记，不断计算其位置，正方形的分布受界面生长动力学的显著影响，利用此修正的偏心正方形算法可良好地反映出等轴晶的形貌，且枝干密集。

步骤S7：界面单元生长。

对于状态值为State(i,j)＝1的界面单元，采用元胞自动机模型建立晶粒的生长模型，计算界面曲率、界面各向异性、溶质浓度场、固相分数的增加量，当许多元胞状态发生改变宏观上就变现为晶粒的增长，溶质的扩散过程是伴随晶粒生长产生的。

采用平均计数法求解界面曲率：

式中，dx为空间步长，N为当前元胞的邻胞数，对于此二维模型计算域，N＝9，为第k₁个邻胞的固相分数；

铝合金凝固过程中的枝晶为具有四重对称轴的立方晶体，界面各向异性函数的计算如下：

f(θ,θ₀)＝1-15εcos[4(θ-θ₀)]

式中，ε为界面能各向异性强度，θ为固液界面法向向量与水平方向的夹角，θ₀为枝晶最优生长角。

采用平均浓度C表征枝晶生长过程中的溶质分布：

C＝(1-f_s)C_l+f_sC_s

式中，C_s与C_l分别为固、液相溶质浓度，界面胞的C为固、液相浓度按固相率的插值，固液界面的溶质浓度为：

式中，与/>分别为界面固相溶质浓度和界面液相溶质浓度，k为溶质平衡分配系数；固液界面单元界面局部平衡认为界面液相溶质浓度与界面实际溶质成分的溶质差是枝晶生长的驱动力，界面液相溶质浓度/>的计算式如下：/>

式中，T^*为界面温度，T_l为平衡液相线温度，m_l为液相线的斜率，Γ为Gibbs-Thomson系数。

固相分数的计算式如下，固相分数增加的过程即为晶粒的生长过程：

式中，L₁(t+dt)通过步骤S6计算，θ为界面法向与水平方向的夹角，计算式如下：

步骤S8：更新元胞状态。

当状态值为State(i,j)＝1的界面元胞的固相分数f_s＝1时，将其状态值修改为State(i,j)＝2，并将界面元胞的液相浓度修改为0。

S9、计算糊状区的渗透率。

结合Darcy定律，计算铝合金液相在糊状区的渗透率。

Darcy定律认为渗透率与一次枝晶臂间距有关，公式如下：

结合Darcy定律，推导得出格子单位下多孔介质渗透率的计算式为：

式中，表示流体动力粘度，Δρ表示入口和出口的密度差，u_x表示x方向的速度，N_x,N_y表示x、y方向的网格数。

S10、不断进行迭代计算，判断是否满足结束条件。若不满足，返回进行步骤S3，若满足，输出模拟结果，凝固结束。

本实例提供的方法可提高捕获模型的准确性，大大改善程序的可行性。本实例提供的方法可研究等轴晶、柱状晶的生长过程及其周围的溶质偏析情况以及温度场的分布情况。本实例提供的方法可模拟不同取向下流场对等轴晶、柱状晶的形貌变化。本实例提供的方法可研究不同过冷度、流场强度、冷却速率、对等轴晶生长的影响。

图2为本发明修正的偏心正方形算法的示意图。图3为本发明算法的类视图，包含了编程计算的所有方法。

图4为本发明无流场时枝晶最优生长角θ₀为45度时枝晶的浓度场分布，可以清晰地看到枝晶的二、三次枝晶臂。

图5为本发明流场条件为强迫对流时枝晶最优生长角θ₀为45度时枝晶的浓度场分布，可以清晰地看到与无流场时枝晶形貌的区别。由于流动对铝合金晶粒的影响，沿流场方向枝晶的生长得到促进，而逆流端由于溶质聚集枝晶生长受到抑制，且可以清晰地看到枝晶的二、三次枝晶臂。

图6为本发明流场条件为自然对流时枝晶最优生长角θ₀为45度时枝晶的浓度场分布，可以清晰地看到与无流场时枝晶形貌的区别。由于温度差和浓度差引起的自然对流对铝合金晶粒的影响，上侧枝晶(迎流侧)的生长受到促进，下侧枝晶(背流侧)的生长普遍受到抑制，且可以清晰地看到枝晶的二、三次枝晶臂，还可以看到因为温度差和浓度差产生自然对流而导致的漩涡。

图7为本发明不同取向枝晶在强迫对流条件下形貌的示意图。其中，图7(a)为强迫对流情况下不同取向枝晶生长过程中的溶质场分布图；图7(b)为强迫对流情况下不同取向枝晶生长过程中的流场分布图；图7(c)为强迫对流情况下不同取向枝晶生长过程中的温度场分布图，可以看到很明显的二次枝晶臂，并可测量二次枝晶臂的间距。

图8为本发明择优取向角θ₀为45度时在强迫对流条件下不同入流速度下的形貌的示意图。其中，图8(a)为强迫对流情况下不同入流速度下枝晶生长过程中的溶质场分布图；图8(b)为强迫对流情况下不同入流速度下枝晶生长过程中的流场分布图；图8(c)为强迫对流情况下不同入流速度下枝晶生长过程中的温度场分布图，入流速度越大，枝晶的非对称性越强，可以看到很明显的二次枝晶臂，并可测量二次枝晶臂的间距。

图9为本发明择优取向角θ₀为45度时在强迫对流条件下不同冷却速率下的形貌的示意图。其中，图9(a)为强迫对流情况下不同冷却速率下枝晶生长过程中的溶质场分布图；图9(b)为强迫对流情况下不同冷却速率下枝晶生长过程中的流场分布图；图9(c)为强迫对流情况下不同冷却速率下枝晶生长过程中的温度场分布图，冷却速率越大，一次枝晶臂越粗，二次枝晶臂越多。

图10为本发明择优取向角θ₀为45度时在强迫对流条件下不同过冷度下的形貌的示意图。其中，图10(a)为强迫对流情况下不同过冷度下枝晶生长过程中的溶质场分布图；图10(b)为强迫对流情况下不同过冷度下枝晶生长过程中的流场分布图；图10(c)为强迫对流情况下不同过冷度下枝晶生长过程中的温度场分布图，过冷度越大，枝晶生长加快，枝晶臂不断粗化，二次枝晶生长更加发达，并造成枝晶尖端固液界面处溶质浓度偏高，加重了溶质偏析。

图11为本发明择优取向角θ₀为0度时在自然对流条件下形貌的示意图。

其中，图11(a)为自然对流情况下不同瑞利数时枝晶生长过程中的溶质场分布图；图11(b)为自然对流情况下不同瑞利数时枝晶生长过程中的流场分布图；图11(c)为自然对流情况下不同瑞利数时枝晶生长过程中的温度场分布图，自然对流强度越大，枝晶的非对称性越强，可以看到很明显的二次枝晶臂，并可测量二次枝晶臂的间距，并且可以看到清晰的溶质羽流与因为温度差和浓度差产生自然对流而造成的涡。

图12为本发明随机多个枝晶在强迫对流条件下不同时刻下的形貌的示意图。其中，图12(a)为强迫对流情况下迭代步数为500时枝晶生长过程中的流场、溶质场、温度场分布图；图12(b)为强迫对流情况下迭代步数为3000时枝晶生长过程中的流场、溶质场、温度场分布图。

图13为本发明随机多个枝晶在自然对流条件下不同时刻下的形貌的示意图。其中，图13(a)为自然对流情况下迭代步数为500时枝晶生长过程中的流场、溶质场、温度场分布图；图13(b)为自然对流情况下迭代步数为3000时枝晶生长过程中的流场、溶质场、温度场分布图。

图14为本发明随机多个柱状晶在强迫对流条件下不同时刻下的形貌的示意图。其中，图14(a)为随机多个柱状晶在强迫对流条件下迭代步数为5000时枝晶生长过程中的溶质场、流场、温度场分布图；图14(b)为随机多个柱状晶在强迫对流条件下迭代步数为9000时枝晶生长过程中的溶质场、流场、温度场分布图。

图15为本发明随机多个柱状晶在自然对流条件下不同时刻下的形貌的示意图。其中，图15(a)为随机多个柱状晶在自然对流条件下迭代步数为5000时枝晶生长过程中的溶质场、流场、温度场分布图；图15(b)为随机多个柱状晶在自然对流条件下迭代步数为9000时枝晶生长过程中的溶质场、流场、温度场分布图。

因此，本发明采用上述一种铝合金凝固动力学过程的介观尺度预测方法，能够在节省计算成本的情况下更真实有效地再现铝合金凝固过程中等轴晶、柱状晶形核长大过程的流场、浓度场、温度场分布，并能在不同流场存在的情况下对等轴晶、柱状晶的具体尺寸、二次枝晶臂间距、枝晶臂间凹槽处的凝固生长进行定量分析，且能计算出铝合金液相在糊状区的渗透率，能够更加准确有效地预测铝合金在凝固过程中的特征，同时提高了计算效率和数值稳定性。

最后应说明的是：以上实施例仅用以说明本发明的技术方案而非对其进行限制，尽管参照较佳实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对本发明的技术方案进行修改或者等同替换，而这些修改或者等同替换亦不能使修改后的技术方案脱离本发明技术方案的精神和范围。

Claims (6)

1.一种铝合金凝固动力学过程的介观尺度预测方法，其特征在于：包括以下步骤：

S1、初始化网络；

S2、基于连续形核模型生成形核点；

S3、计算流场；

S4、计算溶质场；

S5、计算温度场；

S6、基于修正的偏心正方形算法捕获规则，将液相单元转变为界面单元；

S7、界面单元生长；

S8、更新元胞状态；

S9、计算糊状区的渗透率；

S10、判断是否满足结束条件，若不满足，进行步骤S3，若满足，输出模拟结果，凝固结束；

步骤S3中，基于格子玻尔兹曼模型将铝合金液体流动过程分解成碰撞和迁移两个部分，基于D2Q9格子玻尔兹曼算法计算流场的控制方程；根据不同的对流条件对流场的控制方程在枝晶边界处应用不同的边界条件，得到枝晶边界处节点的动量分布函数值；对节点的动量分布函数值进行迁移，并在迁移之后施加边界条件，得到各节点的流体速度；

利用D2Q9格子玻尔兹曼算法计算流场的控制方程为铝合金液体流动过程中的节点碰撞方程：

式中，r是位置矢量，t表示时间，dt是时间步长，α＝0,1,...,8表示9个离散速度方向，其中e_α为离散速度分布式：

表示速度的平衡态分布函数，如下式所示：

式中ρ(x,t)代表t时刻LBM模型中格子的密度，dx为空间步长，u为节点的流体速度，w_α是权系数，如下式所示：

无量纲松弛时间τ_f与流体动力黏度之间的关系为：

F_i为计算流场产生的浮力所添加的源项，当对流条件为强迫对流时，F_i＝0；当对流条件为自然对流时，F_i的计算式如下：

式中，F为自然对流引起的浮力，由温度差和浓度差决定，计算式如下：

F＝gρ₀β_T(T-T₀)+gρ₀β_C(C-C₀)

式中，g是重力加速度，ρ₀是流体的密度，β_T是温度膨胀系数，β_C是浓度膨胀系数，T代表流体的温度，T₀代表流体的初始温度，C代表流体的浓度，C₀代表流体的初始浓度；

当对流条件为强迫对流时，左右边界采用充分发展边界，并给左边界一个恒定的入流速度，为了保证计算区域的质量和动量守恒，网格的上下边界采用周期性边界：

f_α(x,0,t+dt)＝f'_α(x,N_y,t)

f_α(x,N_y,t+dt)＝f'_α(x,0,t)

式中，N_y表示y方向的网格数；

强迫对流下对应宏观的密度和速度可由下式确定：

当对流条件为自然对流时，网格的边界条件采用非平衡外推格式：

其中，是平衡态部分，/>是非平衡态部分；

对应宏观的密度和速度可由下式确定：

步骤S4中，把流场f_α替换成溶质场g_α，将流场的源项F_i替换为溶质场的源项w_αC_l(1-k)Δf_s，就能得到带源项的计算溶质场的标准格子玻尔兹曼方程，从而解决溶质输运的问题，在枝晶边界处采用非平衡外推的边界条件，得到枝晶边界处节点的浓度分布函数值；对节点的浓度分布函数值进行迁移，并在迁移之后施加边界条件，得到各节点的浓度分布；

利用D2Q9格子玻尔兹曼算法建立溶质场模型，溶质场演化方程为：

式中，表示流体浓度的平衡态分布函数，如下式所示：

式中，C(x,t)代表t时刻LBM模型中格子的浓度；

无量纲松弛时间τ_g与溶质扩散系数D之间的关系为：

网格边界条件采用非平衡外推格式，对应宏观的浓度可由下式确定：

源项中，Δf_s是界面元胞在一个时间步长内的固相增量，与步骤S7有关，k为平衡分配系数，C_l为液相溶质浓度，与步骤S7有关，步骤S3得到的流体密度和速度也影响步骤S4中的浓度；

步骤S5中，把流场f_α替换成温度场h_α，将流场的源项F_i替换为温度场的源项w_αΔf_sL/c_p，就能得到计算温度场的标准格子玻尔兹曼方程，从而解决热量输运的问题，在枝晶边界处采用非平衡外推的边界条件，得到枝晶边界处节点的温度分布函数值；对节点的温度分布函数值进行迁移，并在迁移之后施加边界条件，得到各节点的流体温度；

利用D2Q9格子玻尔兹曼算法建立温度场模型，温度场演化方程为：

式中，表示速度的平衡态分布函数，如下式所示：

式中，T(x,t)代表t时刻LBM模型中格子的温度；

无量纲松弛时间τ_h与热扩散系数λ之间的关系为：

网格边界条件采用非平衡外推格式，对应宏观的温度可由下式确定：

源项中，L为潜热，c_p为热容，步骤S3得到的流体密度和速度也影响步骤S5中的温度。

2.根据权利要求1所述的铝合金凝固动力学过程的介观尺度预测方法，其特征在于：步骤S1中，将固定的二维计算区域均匀划分为正方形网格，每个正方形网格代表一个单元胞，初始化单元胞，为每个单元胞的状态参数、物性参数、固相率、液相溶质浓度、固相溶质浓度、速度、分布函数赋予初始值；其中单元胞的状态包括三个，即液相胞、界面胞和固相胞，分别对应三个状态参数值State(i,j)＝0、State(i,j)＝1和State(i,j)＝2。

3.根据权利要求2所述的铝合金凝固动力学过程的介观尺度预测方法，其特征在于：步骤S2中，过冷度ΔT通过晶粒形核密度n(ΔT)对应一个一定的晶粒形核概率，通过基于高斯分布的连续形核模型生成形核点；

形核发生在一系列的形核位置上，形核质点为连续而非离散的分布函数，晶粒形核的模型实质为过冷度ΔT关于晶粒形核密度n(ΔT)的函数，过冷度ΔT与晶粒形核密度n(ΔT)一一对应，晶粒形核密度n(ΔT)代表着晶粒形核的概率大小是一定的，所以过冷度ΔT通过晶粒形核密度n(ΔT)对应一个一定的晶粒形核的概率；新生的晶核密度可表示为：

其中晶粒形核密度随过冷度的关系式为：

式中：n_max为最大晶粒密度，ΔT为元胞的过冷度，ΔT_max为最大形核过冷度，ΔT_σ为过冷度标准差，这两个式子的组合即为晶粒的形核模型。

4.根据权利要求1所述的铝合金凝固动力学过程的介观尺度预测方法，其特征在于：步骤S6中，基于修正的偏心正方形算法捕获规则，采用一系列对角线方向与枝晶生长方向相同的正方形，辅助元胞进行在生长方向上的界面胞捕获和状态传递，将相邻的状态值为State(i,j)＝0液相单元转变为状态值为State(i,j)＝1的界面单元，利用此算法可获得任意生长角度的枝晶；

当计算域中某元胞形核后，其状态参数值由State(i,j)＝0变为状态参数值State(i,j)＝1，可指定一个也可随机赋予一个优先生长角度θ₀，形核后的(i,j)单元受到枝晶生长动力学的驱动，固相分数不断增加，在形核元胞的中心根据元胞的固相分数布置一个有倾斜角度θ₀的正方形，其对角线长度随生长发生偏移，该正方形的半对角线长可由下式计算：

L₁(t+dt)＝L₁(t)+V_n·dt

式中，V_n为固液界面法向生长速度，n为界面法向，

n＝n_xi+n_yj

式中，n_x,n_y为界面法向n的分量，计算式如下：

式中，和/>分别表示固相率在x轴和y轴上的一阶偏导，/>和/>分别为固相率在x轴和y轴上的二阶偏导，其中f_s可通过步骤S7计算；

假设在界面胞中固液界面为平界面，V_n的计算式如下：

V_n＝V_xn_x+V_yn_y

式中，V_x,V_y为固液界面在x,y方向的速度分量，计算式如下：

式中，k为平衡分配系数，D_s为固相扩散系数，D_l为液相扩散系数，为界面(i,j)处液相浓度，/>为界面(i,j)处固相浓度，C_s(i-1,j)为节点(i-1,j)处固相浓度，f_s(i-1,j)为节点(i-1,j)处固相分数，C_s(i+1,j)为节点(i+1,j)处固相浓度，f_s(i+1,j)为节点(i+1,j)处固相分数，C_l(i-1,j)为节点(i-1,j)处液相浓度，C_l(i+1,j)为节点(i+1,j)处液相浓度，C_s(i,j-1)为节点(i,j-1)处固相浓度，f_s(i,j-1)为节点(i,j-1)处固相分数，C_s(i,j+1)为节点(i,j+1)处固相浓度，f_s(i,j+1)为节点(i,j+1)处固相分数，C_l(i,j-1)为节点(i,j-1)处液相浓度，C_l(i,j+1)为节点(i,j+1)处液相浓度；

随着(i,j)单元固相分数的增加，偏移的正方形的四个顶点触碰到周围的液相单元，把相邻的液相单元(i,j+1)，(i,j-1)，(i-1,j)，(i+1,j)的状态值变为界面单元，其状态参数值变为State(i,j)＝1，并把优先生长角度θ₀赋予给这四个被捕获的单元，保持生长方向与初始生长方向的一致性，被捕获的单元又以同样的方式捕获周围单元中的液相元胞，其顶点形成一个新的次级正方形，并继续生长，为避免元胞的八邻居元胞出现未被捕获的情况，当某元胞的固相率为1时，将八邻居元胞中的液相胞全部捕获为界面胞；通过对正方形顶点的元胞进行标记，不断计算其位置，正方形的分布受界面生长动力学的显著影响，利用此修正的偏心正方形算法可良好地反映出等轴晶的形貌，且枝干密集。

5.根据权利要求4所述的铝合金凝固动力学过程的介观尺度预测方法，其特征在于：步骤S7中，对于状态值为State(i,j)＝1的界面单元，采用元胞自动机模型建立晶粒的生长模型，计算界面曲率、界面各向异性、溶质浓度场、固相分数的增加量，当许多元胞状态发生改变宏观上就变现为晶粒的增长，溶质的扩散过程是伴随晶粒生长产生的；

采用平均计数法求解界面曲率：

式中，dx为空间步长，N为当前元胞的邻胞数，对于此二维模型计算域，N＝9，为第k₁个邻胞的固相分数；

铝合金凝固过程中的枝晶为具有四重对称轴的立方晶体，界面各向异性函数的计算如下：

f(θ,θ₀)＝1-15εcos[4(θ-θ₀)]

式中，ε为界面能各向异性强度，θ为固液界面法向向量与水平方向的夹角，θ₀为枝晶最优生长角；

采用平均浓度C表征枝晶生长过程中的溶质分布：

C＝(1-f_s)C_l+f_sC_s

式中，C_s与C_l分别为固、液相溶质浓度，界面胞的C为固、液相浓度按固相率的插值，固液界面的溶质浓度为：

式中，与/>分别为界面固相溶质浓度和界面液相溶质浓度，k为溶质平衡分配系数；固液界面单元界面局部平衡认为界面液相溶质浓度与界面实际溶质成分的溶质差是枝晶生长的驱动力，界面液相溶质浓度/>的计算式如下：

式中，T^*为界面温度，T_l为平衡液相线温度，m_l为液相线的斜率，Γ为Gibbs-Thomson系数；

固相分数的计算式如下，固相分数增加的过程即为晶粒的生长过程：

式中，L₁(t+dt)通过步骤S6计算，θ为界面法向与水平方向的夹角，计算式如下：

6.根据权利要求5所述的铝合金凝固动力学过程的介观尺度预测方法，其特征在于：步骤S8中，当状态值为State(i,j)＝1的界面元胞的固相分数f_s＝1时，将其状态值修改为State(i,j)＝2，并将界面元胞的液相浓度修改为0；

步骤S9中，结合Darcy定律，计算铝合金液相在糊状区的渗透率；

Darcy定律认为渗透率与一次枝晶臂间距有关，公式如下：

式中，λ₁为一次枝晶臂间距；

结合Darcy定律，推导得出格子单位下多孔介质渗透率的计算式为：

式中，表示流体动力粘度，Δρ表示入口和出口的密度差，u_x表示x方向的速度，N_x,N_y表示x、y方向的网格数。