基于PF-LBM模型的自然对流影响枝晶生长模拟

2022-09-05朱昶胜曹雅星马芳兰

兰州理工大学学报 2022年4期

朱昶胜, 曹雅星, 雷鹏, 冯力, 马芳兰

(1. 兰州理工大学计算机与通信学院, 甘肃兰州 730050; 2. 兰州理工大学省部共建有色金属先进加工与再利用国家重点实验室, 甘肃兰州 730050; 3. 兰州理工大学网络与信息中心, 甘肃兰州 730050)

确定性方法通过确定过冷度函数来预测枝晶大小及生长速度[1].由于其本身的特性无法对枝晶生长过程中的随机部分进行模拟，应用受到较大限制.随后Brown等[2-3]建立的元胞自动机模型，考虑了枝晶局部界面曲率的影响，增大了模拟尺度，定性地观察了枝晶生长过程，但此方法无法模拟复杂枝晶的准确形貌，且对于体积过大、计算量大的铸件无法达到计算要求.国内外学者利用相场法(PF)对微观组织进行模拟研究，定量模拟了宏观作用下枝晶生长，较好地捕捉了微观组织动力学演变过程的同时，研究各个参数对枝晶生长的影响，简化了计算步骤，提高了计算效率.Wheeler等[4]在建立的WBM模型的基础上，引入相场与溶质场的修正项，建立新的模型，模拟了金属凝固现象.朱昶胜等[5]研究了Al-Cu合金的倾斜枝晶生长动力学和定向凝固过程中的形态转变，特别分析了冷却速率和主晶间距对凝固组织的影响，研究表明，冷却速率和主晶间距会在一定程度上影响倾斜枝晶的生长角度.

在液态金属凝固过程中，对流对枝晶形貌及周围温度场、溶质场等都有明显影响，对于传统流体力学计算方法，现已存在大量研究[6-7]，其中，基于Navier-Stoke方程[8]的流场数值计算方法是描述粘性不可压缩流体动量守恒的运动方程，可进行离散求解，但仅适用于单向流问题，且不能满足连续介质的假设条件，这些缺点将会造成流场计算不易收敛、计算困难等问题.目前大多数研究是基于外力场作用导致的强迫对流，忽略了由于温度和溶质分布变化引起的自然对流，从而导致计算结果出现偏差，动量守恒方程存在计算量大和计算效率低等问题.格子玻尔兹曼方法(LBM)集合大量虚拟粒子，在离散格子中进行碰撞、迁移，形成流动现象，该方法自身具有并行计算的特点，程序易于实施，具备流体结构互相作用描述简单的优势，可以处理复杂的边界和外力项问题.Miller等[9]提出了一种在对流情况下保持各向同性不变的动力学格子Bolzmann方法，耦合相场法，研究了不同瑞利数和普朗特数下纯金属Ga的生长过程以及固、液相变与对流之间的相互作用.

本文将计算枝晶生长的相场法与计算流场、浓度场和温度场的格子玻尔兹曼方法进行耦合，建立了二维PF-LBM模型，模拟了 Al-4.5% Cu枝晶在自然对流作用下的凝固过程，研究了枝晶形貌、流场、温度场和溶质场的变化以及各向异性强度对枝晶生长的影响.

1 数学模型

1.1 枝晶生长相场模型

相场模型是一种建立在热力学基础上，描述系统演化动力学的模型，其核心思想就是引入一个或多个变量，用弥散界面代替尖锐界面来描述界面[10].在固、液两相系统中引入一个序参量，也就是相场变量φ(r,t)来标识物质状态，液相取值为0或-1，固相取值为1.Ginzburg-Landau方程的关键在于具体体系给出的适当自由能泛函，本文采用的自由能泛函表达式如下：

(1)

根据热力学第二定律，随着凝固进行，体系自由势能逐渐减小，将高温液态金属看作受力后极易变形，且切应力与变形速率成正比的牛顿流体，其数值计算要考虑浓度、温度分布影响，需要加入连续方程和动量方程[11]，利用先行不可逆动力学方程推导出自由能密度函数为

(2)

式(1、2)中：φ为相场变量；F为体系自由能；f为双稳态事函数；H为体系量纲的焓；c为浓度；ε为相场项系数；δ为溶质场梯度系数;D(φ)是溶质扩散率，fc为自由能一阶偏导数，fcc为自由能二阶偏导数.

此外，各向异性是需要考虑的重要因素之一，表达式如下:

ε(θ)=ε(1+γcoskθ)

(3)

加入各向异性得到的相场方程为

(4)

其中:方程(3)中ε为方向矢量;θ为法向与生长主轴夹角;k为各向异性模数;γ是各向异性强度，数值越大强度越大，方程(4)中ε′、ε″表示对θ的一阶、二阶导数，φxy、φyy表示对x、y的二阶偏导数.

1.2 LBM流场模型

LBM是一种基于介观模拟尺度的计算流体力学的方法，直接从离散模型出发，应用质量、动量和能量守恒的规律，从不同角度建立起宏观与微观、连续与离散之间的联系.格子方程模型中DnQm[12-13]系列最具代表性，其中n为空间维数，m为离散速度数.采用基于单松弛(BGK)模型的D2Q9模型，流体碰撞过程中，对应的演化方程表达式为

fi(x+ciδt,t+δt)-fi(x,t)=

(5)

二维D2Q9模型和三维D3Q19常用来计算流动问题，如图1所示，D2Q9模型是根据不同网格确定的平衡分布函数，其中一个虚拟粒子有9个离散速度，由于边界条件不参与演化，因此要根据已知边界条件来确定粒子分布函数之后才能进行下一步计算，对应的离散速度表达式为

(6)

式中：ωi为权函数，i取0～8，cs为格子声速.

图1 晶格结构Fig.1 Lattice structure

LBM用于计算受对流和扩散影响的浓度场和温度场，其分布函数演化方程与流体分布函数类似，如下式所示[14-15]：

gi(x+eiΔt,t+Δt)-gi(x,t)=

(7)

hi(x+eiΔt,t+Δt)-hi(x,t)=

(8)

式中：gi(x,t)为浓度场分布函数，hi(x,t)为温度场分布函数，Gi、Hi分别为枝晶生长过程中浓度梯度、温度梯度造成的溶质源和温度源项.

1.3 数值模拟计算

由于铜铝合金属于面心立方晶系，具有高耐热性和高韧性，因此选取Al-4.5%Cu二元合金，表1为合金物性参数.初始状态下，取溶质浓度为c0，温度为T0，模拟区域为800×800的正方形网格，具有初始浓度和初始温度，模拟中心为400×400，四壁设为温度恒定的固定边界，采用非平衡外推格式.

在计算区域内放置一个或多个半径为R的晶核，并选择最优方向生长.通过式(4)计算截面平衡成分，式(5)计算固液相流动，凝固场中浓度分布发生改变通过式(2)表示，式(1)完成相场的迁移，计算枝晶生长过程.当完成一个时间步长的传输计算后，由于自然对流作用下枝晶生长还需考虑排出溶质和释放潜热所造成的温度和浓度分布不均匀，导致固相分数、液相浓度增加，若界面网格未完全凝固，则利用式(7)将释放出的溶质增量Δc加到网格中剩余液相溶质分布函数，利用式(8)将释放的温度增量ΔT加到同网格的温度分布函数，若完全凝固，则将Δc和ΔT加到相邻网格液相溶质分布函数中，获得每个时间步长的流场、温度场以及浓度场分布，其中流场的固液界面使用反弹格式(5、6)实现枝晶对速度分布影响.此模型包含了耦合液相流动、溶质和热量传输与枝晶生长的物理机制，通过重复上述步骤直到模拟结束就可实现流场计算和相场模拟的耦合.

2 模拟结果与分析

2.1 LBM模拟自然对流验证

在枝晶生长过程中，流动促进溶质和热量的传输，对浓度、温度分布产生影响，流场直接影响耦合模型和算法的正确性，因此需要对计算LBM模型进行验证.封闭方腔是一个水平放置的二维腔体，如图2所示，方腔左壁温度高于右壁温度并保持恒定，上下两壁设置为绝热，高为H，方腔左壁被加热至Th，右壁被冷却至Tc并保持恒定，Th>Tc.以方腔左下角作为坐标原点，u的方向为x轴方向，重力的反方向为y轴方向.其中两个最基本无量纲参数为瑞利数Ra和普朗特数Pr，Ra描述方腔中自然对流和换热强度参数，表示自然对流的强度，Ra越大，自然对流强度越大，普朗特数记为Pr=v/α.平均Nu数和Ra之间存在着幂律关系：Nu=a(Ra)b，式中a=0.142和b=0.299.

图2 封闭方腔示意图Fig.2 Schematic diagram of closed square cavity

在开始阶段，不论是否存在对流，枝晶生长速度都很大，如图3所示，随着枝晶不断向界面前沿释放溶质，导致溶质富集，各个尖端生长速度下降，扩散达到平衡时，速度趋于稳定，稳态的速度与LGK的解析值非常接近，因此，应用LBM模型和算法计算自然对流是可行的.

图3 PF-LBM模型与LGK模型理论值对比

2.2 自然对流对单晶粒生长的影响

将择优取向与水平夹角为0度的晶核放置在计算区域中心，设定半径为R=10 mm，初始过冷度ΔT0=10 K、溶质初始c0=0.45 wt.%，并以枝晶中心为原点建立平面直角坐标系，将两条数轴分别置于水平位置与垂直位置，将枝晶分为4个象限，如图4a所示.枝晶生长是一种动态过程，需要考虑自然对流对枝晶的影响.在相变过程中，因潜热释放和溶质再分配会使液相中溶质浓度和温度分布不均匀，导致凝固场中产生温度梯度、浓度梯度，使流体密度不均匀，从而在重力作用下引起流动，产生漩涡.枝晶生长初期，自然对流强度较弱，液体流动方向由方腔底部绕枝晶至顶部，由于枝晶的生长规律以及边界条件的限制，流场出现左右对称的两个涡流，如图4a1和图4b1所示，此时枝晶尚不发达，枝晶整体较为对称.随着潜热释放以及溶质排出，界面温度和浓度升高，流动达到一定强度时，产生四个漩涡，四个垂直方向生长的主枝晶臂均逐渐变粗，枝晶对称性遭到破坏，垂直向上的一次枝晶臂受到抑制，垂直向下的一次枝晶臂生长受到促进，水平方向一次枝晶臂与纯扩散条件下的枝晶生长趋势近似，如图4a2和图4b2所示.

图4 自然对流条件下,不同时刻单晶粒相场、溶质场演变 Fig.4 The evolution of single dendrite phase field and concentration field at different times under the condition of natural convection

枝晶臂尖端速度和浓度随时间变化曲线如图5所示.开始阶段，晶粒处于一个低浓度、低温度的情况下，容易获得较大的生长动力，一次枝晶臂的尖端以较大速度开始生长.同时，随着主枝晶臂向固液界面处延伸，枝晶排出的溶质进行再分配，潜热的释放富集在固/液界面，导致溶质成分增加，枝晶前沿液相浓度也随之升高，并使固相溶质浓度低于液相初始浓度，过冷度随着液相浓度增加而降低，枝晶生长驱动力与开始时相比有明显降低，生长速度也明显下降.一段时间后，固相初始温度与液相初始温度基本达到平衡状态，当t=0.52×10-3s时，生长驱动力基本保持在相对稳定区间，枝晶尖端生长速度基本达到稳定值.由图5a可知，垂直向下的一次枝晶臂尖端生长速度最高，达到0.021×10-3m/s，垂直向上的一次枝晶臂尖端生长速度最低，为0.19×10-3m/s，水平方向一次枝晶臂生长速度相较于纯扩散条件下的情况略高一些.在此过程中，液相流动将溶质从方腔底部冲刷至顶部，垂直向下一次枝晶臂生长速度快，排出溶质多，富集主要集中在枝晶前沿，即使液相流动会将溶质从底部传输到顶部，垂直向下枝晶臂生长快，前沿富集的溶质还未传送出去，由图5b可知，垂直向下一次枝晶臂尖端稳态浓度依旧明显高于垂直向上一次枝晶臂.

其他物性参数保持不变，加入起伏结构和能量起伏观察枝晶生长过程.观察发现，起伏结构和能量起伏为晶核形成提供了能量，凝固界面在起伏的作用下失去稳定性，形成二次或更高次枝晶臂.无起伏时，一次枝晶臂表面光滑，未出现侧向分枝，加入起伏后，一次枝晶臂溶质含量最低，与过冷液体界面之间稳定性降低，一次枝晶臂上出现凸起，形成侧向分枝，随着凝固进行,二次枝晶臂变粗，间距增大.

图5 单晶粒一次枝晶臂尖端速度，浓度对比图Fig.5 Comparison diagram of single dendrite primary dendrite arm tip velocity and concentration

整个模拟过程中枝晶外部区域受自然对流影响显著，涡流带动溶质扩散，导致上游区域枝晶臂尖端固/液界面溶质浓度降低，释放出的溶质增多，浓度略高于下游区域枝晶臂尖端浓度，如图6b和图6d所示.由于整个枝晶曲率基本相同，可知上游区域枝晶臂总过冷度大于下游区域，过冷度较小处界面前沿排出的溶质和热量传输至远离界面区域时间较长，上游区域枝晶生长驱动力更大，时间越久，枝晶的不对称性越明显，由图6a和图6c可知，3、4象限内枝晶臂尖端生长速度大于1、2象限内枝晶臂，其中垂直向下一次枝晶臂与3、4象限内二次枝晶臂生长速度近似，垂直向上的一次枝晶臂与1、2象限内二次枝晶臂生长速度近似.

图6 单晶粒二次枝晶臂尖端生长速度、浓度对比图Fig.6 Comparison diagram of growth speed and concentration at the tip of single grain secondary dendrite arm

2.3 自然对流情况下各向异性强度对单晶粒生长的影响

各向异性强度是指界面表面张力、界面厚度以及界面动力学的各向异性程度，记为γ.图7为自然对流情况下t=0.1 s时，各向异性强度分别为0.3、0.05、0.07时枝晶生长形貌.在自然对流存在且各向异性强度较大时，枝晶沿着〈100〉方向生长，随着界面前沿愈加不稳定，枝晶周围温度梯度、热扩散速度相继减小，表面张力以及液膜厚度等相关参数使枝晶周围扰动变明显，枝晶尺寸变大，枝晶干变细长，内部凹陷和尖端尖锐程度变明显，二次枝晶臂间距增大，幅值减小，根部出现“颈缩”现象.自然对流增强了各向异性强度对枝晶形貌的影响，各向异性强度为0.05时出现明显“颈缩”现象.

图7 t=0.1 s各向异性强度不同时单晶粒相场、溶质场示意图

由于下游区域溶质来不及扩散，大量溶质富集在枝晶尖端，下游区域枝晶臂尖端凝固，析出溶质浓度变高，远离枝晶臂的液相区域溶质浓度几乎不受流动影响，溶质浓度近似，直到溶质扩散和对流迁移达到动态平衡时，枝晶尖端溶质以及浓度也趋于稳定，上游枝晶生长受到促进，生长速度较为平缓，下游枝晶受到抑制，生长速度有明显下降趋势，如图8所示.各向异性强度在0.01～0.03时枝晶速率斜率最大，生长速度变化较大，取值在0.03～0.07时斜率趋于一致，生长速度基本相同，取值在0.07之后，枝晶形貌出现失真现象，与理论值的变化不符，如图9所示.

图8 各向异性强度不同时枝晶生长速度、浓度对比图

图9 t=0.1 s时各向异性强度与生长速度对比图Fig.9 Comparison diagram of anisotropic strength and growth rate at t = 0.1 s

2.4 自然对流对多晶粒生长的影响

多晶粒生长过程与单晶粒相比，生长情况更复杂，不仅会受到自身生长动力学和自然对流的影响，而且还受到相邻枝晶的制约和限制.为了更好地观察自然对流对多晶粒的影响，将5个枝晶置于计算区域内，其中1个放置在区域中心，其余4个对称分布在中心枝晶周围，各向异性强度γ=0.07，设置择优取向与水平夹角为0°，其他条件与单晶粒生长过程模拟条件相同.初始阶段，自然对流不明显，每个枝晶生长不受相邻枝晶的影响.当凝固时间为0.03 s时，5个晶粒生长速度大致相同，分别在固/液界面处，形成溶质富集边界层，色卡显示最高体积分数到达约为5.0%，此时5个晶粒生长不受影响，溶质富集边界没有接触，如图10所示.当t=0.07 s时，枝晶生长至尖端相互接触，各自排出的溶质富集融合成更大的边界层，中心枝晶周围溶质大量富集，不能有效地传输出去，生长受到抑制，外围4个枝晶外侧区域受自然对流影响明显，不断排出的溶质凝固潜热富集在枝晶尖端，体积分数约为4.0%，固相排出的溶质被涡流带到液相区域，枝晶生长速度较快.随着凝固时间推移，枝晶臂生长成熟，中间空隙部分溶质体积分数升高至5.2%.另外，与单晶粒在自然对流下生长过程相同，计算区域上游位置的枝晶明显比下游位置的枝晶生长速度快，不同的是，单晶粒在自然对流的情况下产生左右对称的四个涡流，而多晶粒生长时枝晶外侧凹陷处出现大小不等的多个涡流，加速溶质向远处液相扩散，溶质边界层富集的溶质在对流的冲刷下，导致上游枝晶的溶质边界层较薄，部分溶质被冲走促进生长.

图10 自然对流条件下多晶粒相场和溶质场示意图

将枝晶择优取向与水平夹角设置为30°，其他条件与无角度多晶粒生长过程模拟条件相同.由于重力垂直向下，自然对流使上游区域以及枝晶右侧枝晶臂生长较快，较为发达，如图11a和图11b所示.随着凝固进行，由于枝晶不同的生长取向及相对位置，枝晶互相制约对方生长，晶粒在生长方向与水平夹角为30°和120°上生长速度逐渐增大，溶质在中心枝晶周围无法扩散，体积分数高达5.2%，枝晶间内部区域的枝晶臂主要受到相邻枝晶的影响，枝晶外部区域受自然对流影响较为明显，上游枝晶的主枝晶臂右侧二次枝晶臂发达于左侧，甚至出现了三次枝晶臂，如图11c、d所示，自然对流流向由于枝晶生长方向影响，周围漩涡发生改变，不再呈现左右对称状态.

图11 择优取向与水平夹角为30°时多晶粒生长规律示意图 Fig.11 Schematic diagram of multi grain growth law with the angle of 30 degrees between preferred orientation and horizontal

2.5 自然对流情况下各向异性强度对多晶粒生长的影响

当γ=0.03时，由于自然对流，主枝晶臂a1生长得到强化，速度加快，枝晶臂变粗壮，二次枝晶臂a3受到a2以及其他相邻枝晶臂挤压，溶质富集在枝晶尖端，导致a3生长驱动力变小，抑制其生长，如图12a所示.随着各向异性强度增大，枝晶在不同方向上的生长速度明显不同，树枝状逐渐明显，枝晶臂间距越来越小，枝晶臂越来越粗壮，具体表现为La1>Lb1>Lc1，La2>Lb2>Lc2，La3>Lb3>Lc3，如图12b和图12c所示.由于枝晶间相互抑制、挤压，外围枝晶排出的溶质和潜热富集，中心溶液浓度及温度均升高，导致在各向异性强度为0.03时，出现明显“颈缩”现象.

图12 各向异性强度不同时择优取向与水平夹角为30°时多晶粒溶质场示意图 Fig.12 Schematic diagram of multi grain concentration field with the angle of 30 degrees between preferred orientation and horizontal at different anisotropic strength