CN114943054A - 一种流场耦合相场的确定多相烧结的形貌和物理场的方法 - Google Patents

Abstract

本发明提出了一种流场耦合相场的确定多相烧结的形貌和物理场的方法，主要包括烧结相场模型的建立以及算法的提出，具体步骤如下：构建描述多相烧结过程的多项式型自由能；构建Navier‑Stokes方程、Cahn‑Hillard方程与Allen‑Cahn方程耦合的多相烧结相场模型；对Navier‑Stokes方程、Cahn‑Hillard方程与Allen‑Cahn方程构造空间差分格式和时间差分格式进行有限差分的求解；利用SIMPLE算法编写程序得到耦合模型的数值解；输出保守型相场、速度场变量、压强场变量模拟结果。本发明模型将相场理论与流场Navier‑Stokes方程耦合起来，能够定量分析烧结过程中形貌演化和速度场变化方式，为模拟颗粒重排、熔体渗透等烧结过程奠定了技术基础。

Description

一种流场耦合相场的确定多相烧结的形貌和物理场的方法

技术领域

本发明涉及烧结过程形貌模拟技术领域，具体涉及一种流场耦合相场的确定多相烧结的形貌和物理场的方法。

背景技术

烧结过程是通过加热使粉末颗粒结合成块状整体材料的过程。烧结过程可以分为三类，晶体材料的烧结，非晶材料的烧结(即粘性烧结)以及液相烧结。表面张力作用下引起的颗粒变形和凝聚是晶体材料固相烧结的主要驱动力，经过多年的发展固相烧结的相关理论模型已经比较成熟。然而大部分的烧结过程均涉及流体相，烧结过程中液相能够润湿固体颗粒并提供毛细力，从而加快烧结颗粒的演化，改善烧结组织和物理性能。作为实验研究的有力补充，理论模型和计算机模拟可以动态地再现液相烧结过程的烧结颈增长，自由流动，气孔消失的机理，并预测液相烧结速度。然而液相烧结过程涉及烧结扩散理论与Navier-Stokes流体方程的耦合问题难以解决；模拟固液气不同相和固相内不同晶粒的形貌演化较为困难，导致目前仍然缺少对液相烧结完整过程的形貌演变行为进行研究和预测的模型。

近十年来，液相烧结过程的模拟受到了越来越多的关注和研究。Aldazabal等人于2004年发表论文“Simulation of liquid phase sintering using the Monte Carlomethod”采取蒙特卡洛的方法，主要使用基体中溶质的浓度，溶质的扩散和表面能作为变量对液相烧结过程进行模拟。Seong Jin Park等人于2006年发表论文“Finite ElementSimulation of Liquid Phase Sintering with Tungsten Heavy Alloys”采用了有限元方法进行模拟。Abdul Malik Tahir等人于2015年发表论文“Modeling of the primaryrearrangement stage of liquid phase sintering”，采用了相场模型模拟液相烧结最初的晶粒的重排过程。

相较于以上所使用的方法，对于模拟固液气多相烧结过程，相场方法物理意义明确，模拟结果与实验结果吻合较好。W.Villanueva等人于2009年发表论文“Multicomponentand multiohase simulation of liquid-phase sintering”，在这篇文章中，他们采用了Allen-Cahn方程去区分不同相的烧结过程，但是没有对Allen-Cahn方程区分固相中不同晶粒的现象进行模拟。

发明内容

有鉴于此，本发明目的在于提出一种流场耦合相场的确定多相烧结的形貌和物理场的方法，该方法能够展现多相烧结晶粒形貌和流体颗粒内部物理场的演变过程，能够分析多颗粒中气孔的演化规律。

为了达到上述目的，本发明采用如下技术方案：

一种流场耦合相场的确定多相烧结的形貌和物理场的方法，包括以下步骤：

S1：通过实验方法获取需要的材料物性参数，定义保守型相场变量、相场序参量、速度场和压强场；考虑固相晶粒、液相颗粒和气孔三相的共存相互作用，采用相场变量进行多相耦合，构建多相烧结的多项式型体系自由能；

S2：根据S1构建的体系自由能，建立保守型相场变量演化Cahn-Hilliard方程，用于描述固液气不同相，考虑流体速度场的作用，加入对流项；取向场演化Allen-Cahn方程，用于描述固相晶粒与液相颗粒以及气孔的相互作用；速度场、压强场演化的Navier-Stokes方程和外加应力项设置为界面处表面张力作用，得到演化方程组形成多相烧结相场模型；

S3：根据S2构建的多相烧结相场模型，先进行无量纲化，再使用半交错网格进行有限差分处理，半交错网格中不同网格使用中心格式进行投影；对于方程组进行时间空间的差分处理：时间差分上，选用二阶龙格库塔法；空间差分上，对流项使用迎风差分格式，其余项使用中心差分格式进行求解；

S4：根据步骤S2构建的多相烧结相场模型与S3中构建的差分格式，设定初始条件和边界条件；基于SIMPLE算法耦合求解形貌演化和物理场演化，其中速度和压强计算，需先求解预测步的速度和压强，再进行修正，进行迭代直到满足质量守恒方程的限制条件，输出相场变量、速度场变量、压强场变量的精确值；

S5：根据步骤S4的计算方法，输出的相场变量、速度场变量、压强场变量的精确值，对输出结果进行可视化；根据可视化的微观组织结构，观察烧结过程中颗粒的形貌演化，液相颗粒内部速度场，探究多相烧结过程演化规律。

步骤S1中，所述构建多相烧结过程的多项式型体系自由能的方法为：

通过第一性原理与实验的方法获取需要的材料物性参数，这些材料物性参数包括材料晶界能、表面能、扩散率、密度、动力粘性系数和表面张力大小；

定义保守型相场变量来区分不同的相：定义保守型相场变量c，在液相中为1，在固相和气相中为0；定义变量φ，在固相中为1，在气相和液相中为0；定义一系列相场序参量η_i，i＝1,2,3…m，代表m个不同取向的晶粒，在第i个固体晶粒中，η_i为1，在其他晶粒中η_i为0；因此有固相第i个晶粒中φ＝1,c＝0,η_i＝1,η_j＝0,j≠i；液相中φ＝0,c＝1,η_i＝0,i＝1,2,3…m；气相中φ＝0,c＝0,η_i＝0,i＝1,2,3…m，界面处φ,c,η连续变化；

基于Landau相变理论构建的描述固液气多相烧结过程的多项式型局部自由能密度，其表达式如下：

式中，f(c,φ,η₁,η₂,…η_m)表示局部自由能密度；其中，A₁,A₂,A₃,A₄,A₅,A₆,A₇为唯象参数，与材料的物理参数有关；

利用构建的局部自由能密度f(c,φ,η₁,η₂,…η_m)，基于热力学理论构建多相烧结过程的多项式型体系自由能，构建体系自由能如下：

式中F，为体系自由能；其中，积分号内部第一项表示多相体系的局部自由能密度函数；第二项至第四项表示界面梯度能密度；其中κ_c、κ_φ、κ_η表示梯度项系数。

步骤S2包括以下步骤：

S2.1建立区分不同相的保守型相场变量演化Cahn-Hillard方程，在保守型相场变量演化Cahn-Hillard方程中加入对流项

计算公式如下：

其中M_c表示液相迁移系数；v表示流体微团的速度；t表示时间；

对于保守型相场变量φ，由于固体不存在流体对流引起的相迁移，直接利用保守型相场变量演化Cahn-Hilliard方程，计算公式如下：

其中M_φ表示固相迁移系数；

S2.2建立区分固相内不同晶粒的向场演化Allen-Cahn方程，对于相场序参量η，直接利用向场演化Allen-Cahn方程，计算公式如下：

其中L是自由界面的迁移率；

S2.3对于空间场流体内速度和压强的演化，根据Navier-Stokes方程动量守恒与质量守恒方程式得出，将多相烧结过程中流体当作不可压缩粘性牛顿流体，因此Navier-Stokes方程计算公式如下：

其中，ρ(c,φ)表示流体微团的密度；p表示压强；τ表示流体微团受到粘滞应力，其表达式为

其中μ表示动力粘性系数，在流体相内部使用实验测试的结果，在相界面处使用保守型相场变量的插值求得，界面的密度和动力粘性系数的大小表达式如下所示：

当φ≤0.8时，

当φ＞0.8时，

其中ρ_gas、ρ_liquid、ρ_solid分别表示气相、液相、固相的密度大小；μ_gas、μ_liquid、μ_solid分别表示气相、液相、固相的动力粘性系数大小；气液两相为流体，存在动力粘性系数，使用倒数的插值杜绝模拟过程中由于计算误差引起的某空间点保守型相场变量不在0到1之间，导致该点物理参数超过单相内的物理参数的范围引起计算的弥散；假设固相动力粘性系数的大小为无穷大，粘滞力过大，固体颗粒内部流场没有速度的演化；

Navier-Stokes方程中，S表示流体微团受到的外加应力，多相烧结过程中流体微团受到的外加应力主要来自相界面的表面张力大小，其表达式如下：

若c≠0且φ≠0，

若φ＝0,

若c＝0,

其中若c＝0，为固气界面，

表示气液界面受到张力大小，固液界面附加应力的大小；若φ＝0，为气液界面，

表示气液界面附加应力的大小；若c≠0且φ≠0，为固液界面，两个保守型相场变量共同作用，

表示固液界面附加应力的大小。

步骤S3包括以下步骤：

S3.1先对不同参数进行无量纲化过程，利用选定基础参数对每个物理量除以选定参数得到同量级、模拟截断误差小的无量纲方程组；

S3.2进行时间差分处理，使用二阶精度的龙格库塔法，设每个计算时间步长为Δt，令n为演化的第n个时间步长；将Cahn-Hilliard方程、Navier-Stokes方程、Allen-Cahn方程进行移项，有

分别对三者使用二阶精度的龙格库塔法表达式如下：

v⁽¹⁾＝vⁿ+Δtf(vⁿ)

vⁿ⁺¹＝1/2vⁿ+1/2(v⁽¹⁾+Δtf(v⁽¹⁾))

c⁽¹⁾＝cⁿ+Δtg(cⁿ)

cⁿ⁺¹＝1/2cⁿ+1/2(c⁽¹⁾+Δtg(c⁽¹⁾))

η⁽¹⁾＝ηⁿ+Δth(ηⁿ)

ηⁿ⁺¹＝1/2ηⁿ+1/2(η⁽¹⁾+Δth(η⁽¹⁾))

其中vⁿ、cⁿ、ηⁿ分别表示n时刻速度场的值、保守型相场变量的值、相场序参量的值；v⁽¹⁾、c⁽¹⁾、η⁽¹⁾分别表示时间半推进步速度场的值、保守型相场变量的值、相场序参量的值；vⁿ⁺¹、cⁿ⁺¹、ηⁿ⁺¹分别表示n+1时刻速度场的值、保守型相场变量的值、相场序参量的值；

S3.3空间差分，为防止棋盘式的速度分布，使用半交错网格进行空间网格划分，标量使用坐标网格(I,J)；u表示x方向速度分量，使用坐标网格(i,J)；w表示y方向速度分量，使用坐标网格(I,j)；坐标网格(i,J)是坐标网格(I,J)向水平速度方向错开半个网格、坐标网格(I,j)是坐标网格(I,J)向竖直速度方向错开半个网格；

S3.4对于步骤S2.1中Cahn-Hillard方程中加入的对流项

由于半交错网格，速度v、速度u向保守型相场变量c网格点上进行插值计算得出计算表达形式为：

使用迎风差分格式，离散格式如下所示；

若

若

若

若

式中，u_(i+2,J)表示(i+2,J)坐标点处x方向速度、u_(i+1,J)表示(i+1,J)坐标点处x方向速度、u_(i,J)表示(i,J)坐标点处x方向速度、u_(i-1,J)表示(i-1,j)坐标点处x方向速度；w_(I,j+2)表示(I,j+2)坐标点处y方向速度、w_(I,j+1)表示(I,j+1)坐标点处y方向速度、w_(I,j)表示(I,j)坐标点出处y方向速度、w_(I,j-1)表示(I,j-1)坐标点处y方向速度；c_(I,J+1)表示(I,J+1)坐标点处保守型相场变量c的值、c_(I,J)表示(I,J)坐标点处保守型相场变量c的值、c_(I,J-1)表示(I,J-1)坐标点处保守型相场变量c的值。

S3.5对于步骤S2.4中Navier-Stokes方程中对流项

根据运动速度方向的不同选择相应的差分表达，将差分表达水平竖直方向进行分解，x方向表示为(uu_x+wu_y)，y方向表示为(uw_x+ww_y)，式中u_x表示x方向速度u对x方向的偏导数、u_y表示x方向速度u对y方向的偏导数；w_x表示y方向速度w对x方向的偏导数、w_y表示y方向速度w对y方向的偏导数；由于半交错网格，速度w向速度u网格点上进行插值计算，得出计算表达形式为

同理

式中u_(I,j)表示x方向速度u向速度w网格点(I,j)的投影、w_(i,J)表示y方向速度w向速度u网格点(i,J)的投影；u_(i-1,J)、u_(i-1,J+1)、u_(i,J)、u_(i,J+1)分别表示x方向速度u在网格点(i-1,J)、(i-1,J+1)、(i,J)、(i,J+1)；w_(I,j-1)、w_(I+1,j-1)、w_(I,j)、w_(I+1,j)分别表示y方向速度w在网格点(I,j-1)、(I+1,j-1)、(I,j)、(I+1,j)；即x，y方向计算形式如下：

其中导数项差分格式如下：

S3.6对于建立的多相烧结相场模型，其余项均采用二阶中心差分格式，拉普拉斯算子采用五点差分法。

步骤S4包括以下步骤：

S4.1设定初始条件为：随机分布固液气三相，其中初始速度场均为0，初始压强场全部设定为大气压强，在第一个时间步求解得到满足Navier-Stokes方程限制条件的压强场；设定边界条件为周期性边界条件；

S4.2在时间步内，先根据上一时间步的值，计算保守型相场变量和相场序参量在空间的分布，边界条件使用虚网格法下的周期性边界条件；

S4.3代入步骤S4.2得出的保守型相场变量在空间的分布，先忽略压强场p的作用求解Navier-Stokes方程，得出压强修正速度场的预测步值，边界条件使用虚网格法下的周期性边界条件；

S4.4代入步骤S4.3得出的压强修正速度场的预测步值，利用SIMPLE算法得出，压强场的值和压强修正速度场的推进步值表达如下所示：

其中v^*n+1表示压强修正速度场的推进步值；v^*表示压强修正速度场的预测步值；ρⁿ表示n时间步内密度场的值；pⁿ⁺¹表示第n+1时间步内压强场的值，得出第n+1时间步压强场的值pⁿ⁺¹和压强修正推进步值v^*n+1；

S4.5根据步骤S4.4求解出的压强修正速度场的推进步值v^*n+1，代入到Navier-Stokes方程的质量守恒方程

计算推进步速度场的散度值是否满足小于设定的限制参数θ＝0.001，若满足输出速度场和压强场的空间分布，若不满足重复S4.3-S4.5直到满足限制条件；

S4.6根据步骤S4.2求解出的保守型相场变量的空间分布；步骤S4.5求解出的速度场和压强场的空间分布；重复步骤S4.2-步骤4.5，带入龙格库塔法的公式中得出下一时间步n+1时间步的相场变量、速度场、压强场的分布；

S4.7根据步骤S4.6多相烧结过程模拟结果，得出时间步推进后的变量的解，带入S4.2进行下一个时间步的求解，时间推进直到平衡状态。

步骤S5包括以下步骤：

S5.1输出保守型相场变量，相场序参量变量，压强场变量，x、y方向的速度场变量的txt文件；

S5.2将步骤S5.1中得到的txt文件，画图得出空间场内的保守型相场变量，相场序参量变量，压强场变量，x、y方向的速度场变量分布；

S5.3根据可视化的微观组织结构，观察烧结过程中颗粒的形貌演化，液相颗粒内部速度场和压强场的演化规律，探究多相多颗粒烧结致密化过程的演化规律。

和现有技术相比较，本发明具备如下优点：

(1)本发明提供的一种流场耦合相场的确定多相烧结的形貌和物理场的方法，将相场方法与流体Navier-Stokes方程耦合起来；可以灵活的描述任意组成的混合状态，便于统一处理不同相在网格内部所占比例，便于模拟不同相之间的转变和形貌演化。传统的相场方法模拟固态晶粒的形貌以及晶界的演化，而与Navier-Stokes方程的耦合使得研究气泡与液相颗粒内部的速度场和压强场的变化更加有效，获得了准确的、具有预测能力的多相烧结过程的介观尺度模型。

(2)本发明提供的一种流场耦合相场的确定多相烧结的形貌和物理场的方法，能够计算不同初始颗粒分布下液相烧结的形貌演变；能够计算包含多相烧结的任意过程，无须区分不同相变机理、区分不同的烧结方式而使用不同的模拟方式：能够计算颗粒曲率和压强的演化，方便计算颗粒内部应力的大小；能够动态再现晶间气泡消失的过程。

(3)本发明提供的一种计算Navier-Stokes方程、Cahn-Hillard方程与Allen-Cahn方程耦合模型的算法，提供对流项半交错网格下构造的二阶精度的迎风差分格式；构造其余项的中心差分格式，利用虚网格法进行边界条件的处理；采用SIMPLE算法快速、准确的求解压强场的分布；采用迭代计算，对于Navier-Stokes方程、Cahn-Hillard方程与Allen-Cahn方程进行耦合求解。

附图说明

图1是本发明模拟流程图。

图2是空间场内相场变量分布。

图3是半交错网格分布。

图4a是改变迁移系数烧结过程中两个颗粒颈球比随时间变化曲线。

图4b是改变迁移系数烧结过程中两个颗粒收缩率随时间变化曲线。

图5是液相两颗粒轴线上速度分布和相场分布图。

图6是接触角180°三颗粒接触的局部形貌图，其中，图(a)、(b)、(c)、(d)、(e)、(f)分别代表t＝0、t＝5×10⁴、t＝1×10⁶、t＝8×10⁶、t＝2×10⁷、t＝3.5×10⁷时间步长的保守型相场变量c的演化结果。

图7是接触角88°三颗粒接触的局部形貌图，其中图(a)、(b)、(c)、(d)、(e)、(f)分别代表t＝0、t＝2×10⁴、t＝2×10⁵、t＝2×10⁶、t＝1×10⁷、t＝3×10⁷时间步长的保守型相场变量c的演化结果。

图8是接触角60°三颗粒接触的局部形貌图，其中，图(a)、(b)、(c)、(d)、(e)、(f)分别代表t＝0、t＝6×10⁴、t＝2×10⁵、t＝1×10⁶、t＝6×10⁶、t＝2×10⁷时间步长的保守型相场变量c的演化结果。

图9是不同三颗粒接触角，两两颗粒之间颈球比随时间变化关系

图10是随机分布多颗粒气孔演化局部形貌图。图(a)、(b)、(c)、(d)、(e)、(f)、(g)、(h)、(i)分别代表t＝0、t＝2×10³、t＝2×10⁴、t＝6×10⁴、t＝1×10⁵、t＝1×10⁶、t＝5×10⁶、t＝1×10⁷、t＝2×10⁷时间步长的保守型相场变量c的演化结果。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚地描述；以下实例用于说明本发明，但不能用来限制本发明的使用范围；

本实施方式提出了一种流场耦合相场的确定多相烧结的形貌和物理场的方法。该方法主要有：基于Landau相变理论构建的描述固液气多相烧结过程的多项式型自由能；构建Navier-Stokes方程、Cahn-Hillard方程与Allen-Cahn方程耦合的固液气多相烧结相场模型；对Navier-Stokes方程、Cahn-Hillard方程与Allen-Cahn方程构造空间差分格式和时间差分格式进行有限差分的求解；利用SIMPLE算法编写程序得到耦合模型的数值解；输出相场变量、速度场变量、压强场变量计算结果并对结果进行分析处理，如图1所示，具体如下：

(1)构建描述固液气多相烧结过程多项式型体系自由能

定义保守型相场变量来区分不同的相：定义保守型相场变量c，在液相中为1，在固相和气相中为0；定义变量φ，在固相中为1，在气相和液相中为0；定义一系列相场序参量η_i，i＝1,2,3…m，代表m个不同取向的晶粒，在第i个固体晶粒中，η_i为1，在其他晶粒中η_i为0；因此有固相第i个晶粒中φ＝1,c＝0,η_i＝1,η_j＝0(j≠i)；液相中φ＝0,c＝1,η_i＝0(i＝1,2,3…m)；气相中φ＝0,c＝0,η_i＝0(i＝1,2,3…m)，界面处φ,c,η连续变化，空间场每一点均存在对应的φ，c，η_i，如附图2所示；基于Landau相变理论构建的描述固液气多相烧结过程的局部自由能密度，其表达式如下：

式中，f(c,φ,η₁,η₂,…η_m)表示局部自由能密度；其中，A₁,A₂,A₃,A₄,A₅,A₆,A₇为唯象参数，与材料的物理参数有关；

利用构建的局部自由能密度f(c,φ,η₁,η₂,…η_m)，基于热力学理论构建多相烧结过程的多项式型体系自由能，构建体系自由能如下：

式中，F为多项式型体系自由能；其中积分号内部第一项表示多相体系的局部自由能密度函数；第二项至第四项表示界面梯度能；其中κ_c、κ_φ、κ_η表示梯度项系数。

(2)构建耦合的固液气多相烧结相场模型：

建立区分不同相的保守型相场变量演化Cahn-Hillard方程，在保守型相场变量演化Cahn-Hillard方程中加入对流项

对于保守型相场变量φ，代表固相内部保守场变量变化，不存在由于对流保守型相场变量的变化，直接利用保守型相场变量演化Cahn-Hilliard方程，计算公式如下：

其中M_c、M_φ分别表示表示液相迁移系数、固相迁移系数；v表示流体微团的速度；t表示时间；

对于相场演化Allen-Cahn方程，利用相场序参量区分固相内不同晶粒，对于相场序参量η_i，直接利用相场演化Allen-Cahn方程，计算公式如下：

其中L是自由界面的迁移率；

对于空间场流体内速度和压强的演化，根据Navier-Stokes方程动量守恒与质量守恒方程式得出，将多相烧结过程中流体当作不可压缩粘性牛顿流体，因此Navier-Stokes方程计算公式如下：因此质量守恒方程式，密度不随时间变化；且动量守恒方程中存在粘性项；其计算公式如下：

其中ρ表示流体微团的密度；p表示压强；ρ表示密度、μ表示动力粘性系数。由于模拟作用在整个空间场内，因此应对模拟中多相弥散界面处和单相内部的物理参数值需分别给出。在流体相内部使用实验测试的结果，在相界面处使用保守型相场变量的插值求得，界面的密度和动力粘性系数的插值表达式如下所示：

当φ≤0.8时，

当φ＞0.8时，

其中ρ_gas、ρ_liquid、ρ_solid分别表示气相、液相、固相的密度大小；μ_gas、μ_liquid、μ_solid分别表示气相、液相、固相的动力粘性系数大小；气液两相为流体，存在动力粘性系数，使用倒数的插值杜绝模拟过程中由于计算误差引起的某空间点保守型相场变量不在0到1之间，导致该点物理参数超过单相内的物理参数的范围引起计算的弥散；假设固相动力粘性系数的大小为无穷大，粘滞力过大，固体颗粒内部“流场”没有速度的演化。

(3)构造空间差分格式和时间差分格式

首先进行无量纲化过程，使得计算数量级相近，减小计算机造成的截断误差。显式欧拉格式精度较低，在不同相界面处存在突变会引起非物理振荡的出现；隐式格式由于非线性项的存在，较难求解，因此使用二阶龙格库塔法进行时间差分，既方便求解，又可以得到足够精确的值，计算时先得到预测步内不同空间场的值，再将预测步空间场的值输入得到最终下一个时间步的值。

空间网格选择，若在一套均匀直角网格上进行差分，计算过程中会出现速度的棋盘式分布，且在之后计算的每个时间步会持续棋盘式分布无法继续演化。因此应使用半交错网格进行空间网格划分，标量使用坐标网格(I,J)；u表示x方向速度分量，使用坐标网格(i,J)；w表示y方向速度分量，使用网格(I,j)；坐标网格(i,J)是坐标网格(I,J)向水平速度方向推半个网格、坐标网格(I,j)是坐标网格(I,J)向竖直速度方向推半个网格，见附图3，运用半交错网格，如果出现棋盘式的速度分布，在下一个时间步长计算后，棋盘式速度分布会被抹平。构建好半交错网，对于不同坐标下的转化利用中心格式计算。

空间差分格式的构造，需要根据方程的数学特性构造精确度高、不会产生非物理解的差分格式。由于Navier-Stokes方程与加入对流项的Cahn-Hillard方程差分化后对于方程求解稳定性影响最大的为对流项，对于等式右边其余项一般的作用可以使得方程的振荡减少。因此对Navier-Stokes方程与加入对流项的Cahn-Hillard方程的对流项需要使用迎风差分格式进行计算；对于Navier-Stokes方程、Cahn-Hillard方程与Allen-Cahn方程的其余项均使用二阶中心差分格式；拉普拉斯算子使用五点差分法进行计算。迎风格式的构造差分格式时，需要判断该点速度方向，根据不同的速度方向得到不同的差分格式。应先在半交错网格下，构造出速度值在另一网格上的值，再判断其是否大于0后，得到相应的迎风差分格式。

(4)利用SIMPLE算法编写程序得到耦合模型的数值解

首先设定初始场中随机分布一些固液颗粒，得到初始场内相场变量的分布；初始速度场均为0；初始压强场每点压强均为大气压强，第一个时间步内求解得出合理的压强场分布；边界条件采用周期性边界条件。

先根据上一时间步的值，计算保守型相场变量和相场序参量在空间的空间分布，边界条件使用虚网格法下的周期性边界条件；将相场变量在空间的分布预测值带入到忽略压强场p的作用的Navier-Stokes方程进行求解，得出速度场变量的预测步的值，边界条件使用虚网格法下的周期性边界条件。

利用SIMPLE算法得出压强场和速度场的修正表达如下所示：

其中V^*n+1表示压强修正速度场的推进步值；V^*表示压强修正速度场的预测步值；ρⁿ表示n时间步内密度场的值；pⁿ⁺¹表示第n+1时间步内压强场的值，得出第n+1时间步压强场的值pⁿ⁺¹和压强修正推进步值v^*n+1。

将得出的速度场预测步值，带入到上修正表达式中求解出压强修正速度场的推进步值V^*n+1，代入到Navier-Stokes方程的质量守恒方程中▽·(v)＝0，计算推进步速度场的散度值是否满足小于设定的限制参数θ＝0.001，若满足输出速度场和压强场的空间分布，若不满足重复将得出的新的速度场和压强场重新带入到忽略压强场p的作用的Navier-Stokes方程进行求解，重复上述步骤，直到满足限制条件。

根据求解出的保守型相场变量、速度场和压强场的空间分布；重复上述过程，带入龙格库塔法的公式中得出下一时间步n+1时间步的相场变量、速度场、压强场的分布。一个时间步推进完成，重新带入算法第一步进行下一个时间步的推进求解，设定输出每隔一定时间步推进的值。

(5)输出各变量计算结果并对结果进行分析处理

根据步骤(4)输出保守型相场变量，相场序参量变量，压强场变量，x、y方向的速度场变量的txt文件；将得到的txt文件，作图得出空间场内的保守型相场变量，相场序参量变量，压强场变量，x、y方向的速度场变量分布；各变量的空间部分得到的可视化微观组织结构，观察烧结过程中颗粒的形貌演化，液相颗粒内部速度场和压强场的演化规律，探究多相多颗粒烧结致密化过程的演化规律。

实施示例一

本实例按照具体实施方式中的流程，使用尼龙12参数，对两个颗粒的烧结过程的形貌演变和多个物理场进行分析。

本实施例的理论模型为：

所采用的物理参数如下表所示：

表 尼龙12的物理参数

计算两等大颗粒烧结颈和颈球比随时间变化曲线如图4a和图4b所示，颈球比表达为X/D₀，收缩率表达为(4r₀-L)/2r₀，其中r₀为初始颗粒的半径，L为双颗粒演化过程中的长度。可以看到二者同时变化，变化趋势均为开始变化较快随着时间步长的推进变慢。由于开始时不平衡，驱动力较大，之后变化速度明显变慢。M增大，整体演化速率加快，整体演化趋势不变，且最终的稳定的值并不发生改变。理论上计算满足质量守恒定律，双等大颗粒最终形成一个圆形颗粒，该颗粒的半径应为初始半径的

倍，即颈球比最终应等于

收缩率最终应等于

模拟结果中，颈球比收敛于1.411接近

收缩率收敛于0.588接近

因此本文提出的模型满足质量守恒以及相场变量守恒

分析提出模型和算法对于速度场分布的结果，液相颗粒轴线上速度分布和相场分布如图5所示。在液体颗粒接触界面上存在一个滞止点(stagnation point)，该点速度为0且在该点附近水平速度的大小与距滞止点的距离成线性关系，线性变化率与轴向应变率有关，该区域称为速度的接触应变区。远离接触界面的区域称为刚体运动区域，其变化相较于接触应变区变化较慢。再向两边为可以看到为气相区域，此时速度很小且迅速减为0。速度模拟结果与的理论分析结果一致。

实施示例二

本实例按照具体实施方式中的流程，使用尼龙12参数，对三个颗粒的烧结过程进行研究。

本实施例的理论模型为：

改变初始条件，对不同三颗粒初始颗粒接触角的颗粒形貌分布演化结果如图6(a)、(b)、(c)、(d)、(e)、(f)、、图7(a)、(b)、(c)、(d)、(e)、(f)、图8(a)、(b)、(c)、(d)、(e)、(f)所示。

不同接触角两两颗粒之间的颈球比如图9所示，初始阶段两两颗粒之间演化独立，不受颗粒初始排列方式的影响。对于180°接触角的情况，两烧结颈演化各自独立因此其变化趋势与双颗粒极为接近。对于88°接触角的情况，演化速度大于双颗粒演化速度，主要由于三颗粒之间存在扭转，颗粒向中心扭转，二者之间的凹状孔隙造成曲率梯度迁移增大，引起烧结颈形成速度的加快。对于60°接触角的情况，演化速度小于双颗粒演化速度，主要由于颗粒之间互相接触，其中气孔对物质迁移存在阻碍作用，引起烧结颈形成速度减慢。

实施示例三

本实例按照具体实施方式中的流程，使用尼龙12参数，对多个颗粒的烧结过程进行研究。

本实施例的理论模型为：

随机分布多颗粒的演化过程如图10(a)、(b)、(c)、(d)、(e)、(f)、(g)、(h)、(i)所示，演化过程中，烧结初始阶段，烧结颈形成，接触点扩大。中间阶段，孔隙先进行球化，之后被液相填充。中间阶段的颈部生长较快，通过相邻颈部的融合，使得孔隙进行变化，较小孔隙迅速被液相填充，较大孔隙则发生球化现象。中间阶段结束时，孔隙闭合圆状，总孔隙率下降，此时部分大孔隙并未形成完美的圆形。最终阶段，孔隙闭合致密度增加，导致大孔隙被挤压成封闭的圆形，此时孔隙数量减少平均尺寸增加，最终演化为表面能最小的圆形大颗粒，与理论结果吻合较好。

Claims (6)

1.一种流场耦合相场的确定多相烧结的形貌和物理场的方法，其特征在于：包括以下步骤：

S1：通过实验方法获取需要的材料物性参数，定义保守型相场变量、相场序参量、速度场和压强场；考虑固相晶粒、液相颗粒和气孔三相的共存相互作用，采用相场变量进行多相耦合，构建多相烧结的多项式型体系自由能；

S2：根据S1构建的体系自由能，建立保守型相场变量演化Cahn-Hilliard方程，用于描述固液气不同相，考虑流体速度场的作用，加入对流项；取向场演化Allen-Cahn方程，用于描述固相晶粒与液相颗粒以及气孔的相互作用；速度场、压强场演化的Navier-Stokes方程和外加应力项设置为界面处表面张力作用，得到演化方程组形成多相烧结相场模型；

S3：根据S2构建的多相烧结相场模型，先进行无量纲化，再使用半交错网格进行有限差分处理，半交错网格中不同网格使用中心格式进行投影；对于方程组进行时间空间的差分处理：时间差分上，选用二阶龙格库塔法；空间差分上，对流项使用迎风差分格式，其余项使用中心差分格式进行求解；

S4：根据步骤S2构建的多相烧结相场模型与S3中构建的差分格式，设定初始条件和边界条件；基于SIMPLE算法耦合求解形貌演化和物理场演化，其中速度和压强计算，需先求解预测步的速度和压强，再进行修正，进行迭代直到满足质量守恒方程的限制条件，输出相场变量、速度场变量、压强场变量的精确值；

S5：根据步骤S4的计算方法，输出的相场变量、速度场变量、压强场变量的精确值，对输出结果进行可视化；根据可视化的微观组织结构，观察烧结过程中颗粒的形貌演化，液相颗粒内部速度场，探究多相烧结过程演化规律。

2.根据权利要求1所述的一种流场耦合相场的确定多相烧结的形貌和物理场的方法，其特征在于：步骤S1中，所述构建多相烧结过程的多项式型体系自由能的方法为：

通过第一性原理与实验的方法获取需要的材料物性参数，这些材料物性参数包括材料晶界能、表面能、扩散率、密度、动力粘性系数和表面张力大小；

定义保守型相场变量来区分不同的相：定义保守型相场变量c，在液相中为1，在固相和气相中为0；定义变量φ，在固相中为1，在气相和液相中为0；定义一系列相场序参量η_i，i＝1,2,3…m，代表m个不同取向的晶粒，在第i个固体晶粒中，η_i为1，在其他晶粒中η_i为0；因此有固相第i个晶粒中φ＝1,c＝0,η_i＝1,η_j＝0,j≠i；液相中φ＝0,c＝1,η_i＝0,i＝1,2,3…m；气相中φ＝0,c＝0,η_i＝0,i＝1,2,3…m，界面处φ,c,η连续变化；

基于Landau相变理论构建的描述固液气多相烧结过程的多项式型局部自由能密度，其表达式如下：

式中，f(c,φ,η₁,η₂,…η_m)表示局部自由能密度；其中，A₁,A₂,A₃,A₄,A₅,A₆,A₇为唯象参数，与材料的物理参数有关；

利用构建的局部自由能密度f(c,φ,η₁,η₂,…η_m)，基于热力学理论构建多相烧结过程的多项式型体系自由能，构建体系自由能如下：

式中F，为体系自由能；其中，积分号内部第一项表示多相体系的局部自由能密度函数；第二项至第四项表示界面梯度能密度；其中κ_c、κ_φ、κ_η表示梯度项系数。

3.根据权利要求1所述的一种流场耦合相场的确定多相烧结的形貌和物理场的方法，其特征在于：步骤S2包括以下步骤：

S2.1建立区分不同相的保守型相场变量演化Cahn-Hillard方程，在保守型相场变量演化Cahn-Hillard方程中加入对流项▽·(vc)，计算公式如下：

其中M_c表示液相迁移系数；v表示流体微团的速度；t表示时间；

对于保守型相场变量φ，由于固体不存在流体对流引起的相迁移，直接利用保守型相场变量演化Cahn-Hilliard方程，计算公式如下：

其中M_φ表示固相迁移系数；

S2.2建立区分固相内不同晶粒的向场演化Allen-Cahn方程，对于相场序参量η，直接利用向场演化Allen-Cahn方程，计算公式如下：

其中L是自由界面的迁移率；

S2.3对于空间场流体内速度和压强的演化，根据Navier-Stokes方程动量守恒与质量守恒方程式得出，将多相烧结过程中流体当作不可压缩粘性牛顿流体，因此Navier-Stokes方程计算公式如下：

▽·(v)＝0

其中，ρ(c,φ)表示流体微团的密度；p表示压强；τ表示流体微团受到粘滞应力，其表达式为▽·τ＝▽·(μ(c,φ)(▽v+▽v^T))；其中μ表示动力粘性系数，在流体相内部使用实验测试的结果，在相界面处使用保守型相场变量的插值求得，界面的密度和动力粘性系数的大小表达式如下所示：

当φ≤0.8时，

当φ＞0.8时，

其中ρ_gas、ρ_liquid、ρ_solid分别表示气相、液相、固相的密度大小；μ_gas、μ_liquid、μ_solid分别表示气相、液相、固相的动力粘性系数大小；气液两相为流体，存在动力粘性系数，使用倒数的插值杜绝模拟过程中由于计算误差引起的某空间点保守型相场变量不在0到1之间，导致该点物理参数超过单相内的物理参数的范围引起计算的弥散；假设固相动力粘性系数的大小为无穷大，粘滞力过大，固体颗粒内部流场没有速度的演化；

Navier-Stokes方程中，S表示流体微团受到的外加应力，多相烧结过程中流体微团受到的外加应力主要来自相界面的表面张力大小，其表达式如下：

若c≠0且φ≠0，

若φ＝0,

若c＝0,

其中若c＝0，为固气界面，

表示气液界面受到张力大小，固液界面附加应力的大小；若φ＝0，为气液界面，

表示气液界面附加应力的大小；若c≠0且φ≠0，为固液界面，两个保守型相场变量共同作用，

表示固液界面附加应力的大小。

4.根据权利要求3所述的一种流场耦合相场的确定多相烧结的形貌和物理场的方法，其特征在于：步骤S3包括以下步骤：

S3.1先对不同参数进行无量纲化过程，利用选定基础参数对每个物理量除以选定参数得到同量级、模拟截断误差小的无量纲方程组；

S3.2进行时间差分处理，使用二阶精度的龙格库塔法，设每个计算时间步长为Δt，令n为演化的第n个时间步长；将Cahn-Hilliard方程、Navier-Stokes方程、Allen-Cahn方程进行移项，有

分别对三者使用二阶精度的龙格库塔法表达式如下：

v⁽¹⁾＝vⁿ+Δtf(vⁿ)

vⁿ⁺¹＝1/2vⁿ+1/2(v⁽¹⁾+Δtf(v⁽¹⁾))

c⁽¹⁾＝cⁿ+Δtg(cⁿ)

cⁿ⁺¹＝1/2cⁿ+1/2(c⁽¹⁾+Δtg(c⁽¹⁾))

η⁽¹⁾＝ηⁿ+Δth(ηⁿ)

ηⁿ⁺¹＝1/2ηⁿ+1/2(η⁽¹⁾+Δth(η⁽¹⁾))

其中vⁿ、cⁿ、ηⁿ分别表示n时刻速度场的值、保守型相场变量的值、相场序参量的值；v⁽¹⁾、c⁽¹⁾、η⁽¹⁾分别表示时间半推进步速度场的值、保守型相场变量的值、相场序参量的值；vⁿ⁺¹、cⁿ⁺¹、ηⁿ⁺¹分别表示n+1时刻速度场的值、保守型相场变量的值、相场序参量的值；

S3.3空间差分，为防止棋盘式的速度分布，使用半交错网格进行空间网格划分，标量使用坐标网格(I,J)；u表示x方向速度分量，使用坐标网格(i,J)；w表示y方向速度分量，使用坐标网格(I,j)；坐标网格(i,J)是坐标网格(I,J)向水平速度方向错开半个网格、坐标网格(I,j)是坐标网格(I,J)向竖直速度方向错开半个网格；

S3.4对于步骤S2.1中Cahn-Hillard方程中加入的对流项▽·(vc)，由于半交错网格，速度v、速度u向保守型相场变量c网格点上进行插值计算得出计算表达形式为：

使用迎风差分格式，离散格式如下所示；

若

若

若

若

式中，u_(i+2,J)表示(i+2,J)坐标点处x方向速度、u_(i+1,J)表示(i+1,J)坐标点处x方向速度、u_(i,J)表示(i,J)坐标点处x方向速度、u_(i-1,J)表示(i-1,j)坐标点处x方向速度；w_(I,j+2)表示(I,j+2)坐标点处y方向速度、w_(I,j+1)表示(I,j+1)坐标点处y方向速度、w_(I,j)表示(I,j)坐标点出处y方向速度、w_(I,j-1)表示(I,j-1)坐标点处y方向速度；c_(I,J+1)表示(I,J+1)坐标点处保守型相场变量c的值、c_(I,J)表示(I,J)坐标点处保守型相场变量c的值、c_(I,J-1)表示(I,J-1)坐标点处保守型相场变量c的值。

S3.5对于步骤S2.4中Navier-Stokes方程中对流项ρv▽·v，根据运动速度方向的不同选择相应的差分表达，将差分表达水平竖直方向进行分解，x方向表示为(uu_x+wu_y)，y方向表示为(uw_x+ww_y)，式中u_x表示x方向速度u对x方向的偏导数、u_y表示x方向速度u对y方向的偏导数；w_x表示y方向速度w对x方向的偏导数、w_y表示y方向速度w对y方向的偏导数；由于半交错网格，速度w向速度u网格点上进行插值计算，得出计算表达形式为

同理

式中u_(I,j)表示x方向速度u向速度w网格点(I,j)的投影、w_(i,J)表示y方向速度w向速度u网格点(i,J)的投影；u_(i-1,J)、u_(i-1,J+1)、u_(i,J)、u_(i,J+1)分别表示x方向速度u在网格点(i-1,J)、(i-1,J+1)、(i,J)、(i,J+1)；w_(I,j-1)、w_(I+1,j-1)、w_(I,j)、w_(I+1,j)分别表示y方向速度w在网格点(I,j-1)、(I+1,j-1)、(I,j)、(I+1,j)；即x，y方向计算形式如下：

其中导数项差分格式如下：

S3.6对于建立的多相烧结相场模型，其余项均采用二阶中心差分格式，拉普拉斯算子采用五点差分法。

5.根据权利要求1所述的一种流场耦合相场的确定多相烧结的形貌和物理场的方法，其特征在于：步骤S4包括以下步骤：

S4.1设定初始条件为：随机分布固液气三相，其中初始速度场均为0，初始压强场全部设定为大气压强，在第一个时间步求解得到满足Navier-Stokes方程限制条件的压强场；设定边界条件为周期性边界条件；

S4.2在时间步内，先根据上一时间步的值，计算保守型相场变量和相场序参量在空间的分布，边界条件使用虚网格法下的周期性边界条件；

S4.3代入步骤S4.2得出的保守型相场变量在空间的分布，先忽略压强场p的作用求解Navier-Stokes方程，得出压强修正速度场的预测步值，边界条件使用虚网格法下的周期性边界条件；

S4.4代入步骤S4.3得出的压强修正速度场的预测步值，利用SIMPLE算法得出，压强场的值和压强修正速度场的推进步值表达如下所示：

其中v^*n+1表示压强修正速度场的推进步值；v^*表示压强修正速度场的预测步值；ρⁿ表示n时间步内密度场的值；pⁿ⁺¹表示第n+1时间步内压强场的值，得出第n+1时间步压强场的值pⁿ⁺¹和压强修正推进步值v^*n+1；

S4.5根据步骤S4.4求解出的压强修正速度场的推进步值v^*n+1，代入到Navier-Stokes方程的质量守恒方程▽·(V)＝0，计算推进步速度场的散度值是否满足小于设定的限制参数θ＝0.001，若满足输出速度场和压强场的空间分布，若不满足重复S4.3-S4.5直到满足限制条件；

S4.6根据步骤S4.2求解出的保守型相场变量的空间分布；步骤S4.5求解出的速度场和压强场的空间分布；重复步骤S4.2-步骤4.5，带入龙格库塔法的公式中得出下一时间步n+1时间步的相场变量、速度场、压强场的分布；

S4.7根据步骤S4.6多相烧结过程模拟结果，得出时间步推进后的变量的解，带入S4.2进行下一个时间步的求解，时间推进直到平衡状态。

6.根据权利要求1所述的一种流场耦合相场的确定多相烧结的形貌和物理场的方法，其特征在于：步骤S5包括以下步骤：

S5.1输出保守型相场变量，相场序参量变量，压强场变量，x、y方向的速度场变量的txt文件；

S5.2将步骤S5.1中得到的txt文件，画图得出空间场内的保守型相场变量，相场序参量变量，压强场变量，x、y方向的速度场变量分布；

S5.3根据可视化的微观组织结构，观察烧结过程中颗粒的形貌演化，液相颗粒内部速度场和压强场的演化规律，探究多相多颗粒烧结致密化过程的演化规律。

Images