基于LBM-CA模型的自然对流下液氢中固空枝晶生长模拟

戴闻骁 雷 刚 郑晓红 梁文清 疏志勇 钱 华

(1 东南大学能源与环境学院 南京 210096)

(3 东南大学航天低温推进剂国家重点实验室合作基地 南京 210096)

(3 航天低温推进剂国家重点实验室 北京 100028)

1 引 言

液氢的储存形式具有综合密度高、方便使用、技术成熟的优点,尤其在航天用的火箭发动机上具有良好的应用前景。由于在液氢的加注,转运等过程中不可避免渗入空气,研究表明,液氢中积累的固空是富氧的[1],当固空中的氧浓度达到一定比例时,液氢系统易发生爆炸[2-4]。因此研究液氢系统中固空凝固过程对提高液氢系统安全性能有重要应用价值。

当液氢流动停止时,凝固释放的潜热和溶质导致局部密度不均将引发自然对流。相关研究表明[5],自然对流将对凝固过程中枝晶的生长产生影响。由于液氢的低温特性,且固空晶体在系统中随机位置形核生长,现有实验设备难以直接清晰观察到固空颗粒的微观形成过程,更难对其成分进行定量测量分析。刘海生[2-3]通过理论计算和掺杂-复温的实验方法测得在液氢系统复温过程中氮氧的浓度变化,得到固空晶体中心贫氧外表富氧的宏观定性结论,并观察到雪花状固空晶体,但未能定量给出氮氧浓度分布情况。

鉴于难以进行微观的实验观察和定量测量,本研究采用数值模拟的方法研究自然对流作用下固空的微观结构。传统的CFD 方法难以模拟凝固过程中晶体具体的微观组织,为将微观组织形成过程与流场、温度场、浓度场相耦合,选用格子玻尔兹曼(LBM)模型耦合元胞自动机(CA)的模拟方法。孙东科[6]耦合LBM 与CA 方法模拟了在复杂流场下枝晶生长的特性,该方法具备了计算速度快、天然并行、能够处理复杂边界等优点。但由于四边形网格的局限性,LBMCA模型多用于模拟铁、铜、铝[7-8]等具有面心立方和体心立方晶体结构的四重对称微观枝晶生长。周伟煜[9]采用LBM-CA 方法模拟了固空四重对称枝晶生长的过程,虽然得到氧浓度在固空颗粒中的分布情况,但与刘海生观察到的雪花状结构不符,即仍未解决网格局限性问题。在自然对流条件下流场受枝晶结构影响较大,因此LBM-CA 方法用于模拟自然对流作用下固空晶体六重对称生长时,网格各向异性将影响计算的准确性。本文模拟中引入减弱因子[10]并优化差分方法[11]减小正方形网格各向异性及离散各向异性的影响。

本研究将空气简化为N2-21%O2的氮氧混合物,以过冷度作为生长驱动力,采用改进的LBM-CA 方法研究在自然对流条件下,溶质析出与潜热扩散对枝晶形貌的影响,探究固空颗粒中氮氧浓度的分布情况。

2 计算模型

本研究中使用LBM-CA模型将固空枝晶生长与流场的演化结合到一起,CA 方法模拟枝晶生长过程,LBM模型模拟流场的演化,两种模型相互耦合模拟在复杂流场下固空的六重对称枝晶生长过程。

2.1 CA模型

网格点的凝固程度用固相率分数fs表示,根据网格点的固相率分数fs,将网格分为液相网格(L),固相网格(S)及界面网格(S/L)。

对仍未完全凝固的网格,即界面网格,以过冷度作为枝晶生长的动力,界面网格的固相率不断增加最终成为新的固相网格,并按照Moore 元胞类型,捕获周围液相网格使其成为新的界面网格(图1)。

图1 Moore 元胞捕获原则Fig.1 Moore cell capture principle

界面处固相率的增长Δfs与枝晶生长速度vn有关:

式中:Δt与Δx分别为CA模型下的时间步长和空间步长。bred为形状减弱因子,早期的CA 方法中,枝晶生长的各向异性依赖于网格,枝晶的生长方向还是沿着网格线的方向。早期这种人为的各向异性重现了枝晶生长存在的各向异性,但是当枝晶的择优生长方向与网格不重合时,CA 方法存在的网格各向异性问题就凸显出来。本研究引入了削弱因子(reduction factor)以减弱四边形网格的影响[10]。

枝晶的生长速度vn为局部过冷度的函数:

其中:μk为各向异性系数:

式中:μk0为平均向异性系数,εμ为各向异性因子,θx-θ0表示界面法向量与择优生长方向的夹角,(°);m表示不同的多重对称性,在6 重对称枝晶生长中m取6。

总过冷度ΔT由热过冷Tt(K),成分过冷Tc(K)和曲率过冷TR(K)组成,即:

式中:Clx为当地液相浓度,%;C0为初始浓度,%;当液相线斜率ml为负时局部浓度越高越不利于凝固的进行。Γ为Gibbs-Thomson 系数,K为界面曲率,其可根据物理意义表示为:

2.2 LBM模型

格子玻尔兹曼方法(LBM)根据离散的Boltzmann方程进行演化,在演化过程中分为迁移和碰撞步,因其天然并然性大大加快了计算速度。

本文采用常用基于单松弛时间的D2Q9模型即LBGK[12]方法,其演化方程为:

式中:fi(x,t)为粒子分布函数为粒子平衡态分布函数,δt为格子时间步长,ei为离散速度,τf为松弛时间Fi为外力源项,可表示为:

式中:cs为格子声速,等于。c,c为格子速度,c=δx/δt。δx、δt分别为格子步长和格子时间步长,F为自然对流的浮力,根据Boussinesp 近似,该项可表示:

式中:ρ0为流体的初始密度,kg/m3;βT,βC分别为温度膨胀系数和浓度膨胀系数;g为重力加速度,m/s2;C0为初始浓度,%;T0为初始温度,K。

离散速度ei,wi为权重系数,当i=0 时,wi=4/9,当i=1—4 时,wi=1/9,当i=5—8 时,wi=1/36。

类似地,溶质场及温度场的演化也可用对应的格子Boltzmann 方程计算:

式中:ΔCl为凝固过程中由于潜热释放和溶质析出对周围浓度,%;ΔTH为凝固过程中由于潜热释放和溶质析出对周围浓度温度,K。

式中:k为溶质分配系数;Lh为单位质量溶液所含凝固潜热,kJ/kg;cp为等压热容,kJ/(kg·K)。

速度场,温度场,浓度场的平衡态分布函数分别为:

流场、浓度场、温度场的松弛时间为运动粘度ν、溶质扩散系数DL和热扩散系数α的函数:

根据LBGK模型将粒子分布函数恢复到密度ρ,(kg/m3);流速u,(m/s);浓度C,%;温度T,K;等宏观物理量。

2.3 CA 与LBM 耦合

初始状态下,在计算区域内放置一个晶核作为枝晶生长起点,根据上述捕获规则,固相网格捕获周围液相网格,使其成为固液界面网格。通过LBM模型得到流场的速度、温度、浓度分布,决定该时间步长内固液界面的浓度和温度。根据界面的温度及浓度,利用CA 法计算该时间步长内的枝晶生长。凝固释放出的潜热和溶质作为源项添加到LBM 演化方程中,在自然对流过程中,因浓度和温度的不均衡,在重力作用下引发的自然对流也需考虑。模拟程序流程图如图2 所示。

图2 LBM-CA模拟流程图Fig.2 Flow chart for simulating dendritic growth

2.4 边界条件

在计算区域充满了一定过冷度和一定浓度的氮氧混合液,热边界均采用非平衡外推,凝固潜热最终释放到外界坏境中保证凝固过程一直进行,浓度场采用无扩散边界条件。计算区域四周及因枝晶生长在计算区域内部形成的复杂固液界面均设为无滑移边界,并在LBM 计算中采用标准反弹格式处理。由于相对于溶质在液相中的扩散系数,固相中扩散系数小3 个量级,所以忽略溶质在固相中的扩散,因而溶质无法以扩散的形式穿过固液界面。

3 模拟结果与讨论

将计算区域划分为300 ×300 的网格,网格尺寸Δx=0.5 μm。LBM模型中的松弛因子取τf=1,因此格子粘度为1/6。t=0 时刻,计算区域中心放置一半径为5 格子的圆形晶核,初始过冷度为0.6 K,冷却速率为5 K/s。模拟中所需的参数如表1 所示[9,13]。

表1 模拟时所用物性参数Table 1 Parameters used in present work

3.1 单个固空枝晶在自然对流下的生长

本部分研究了单个固空枝晶在自然对流下的生长情况。取温度瑞利数Ra=5 000,浓度瑞利数Ra=5 000,择优生长角度θ0与重力方向夹角为0°及90°时,自然对流与纯扩散条件下枝晶生长情况对比如图3、图4 所示(枝晶形貌由黑色轮廓线表示)。在纯扩散条件下,凝固析出的溶质及释放的潜热仅依靠扩散向四周传递,各枝晶界面浓度、温度几乎一致,呈对称分布。由于界面枝晶生长各向异性,晶核沿着特定的方向生长,固-液界面由圆形发展为正六边形。与纯扩散条件相比,在t=0.1 s 时,自然对流条件下枝晶生长表现出的不对称性已相当明显,当择优生长方向不同时,固空枝晶表现出的生长不对称性也不尽相同,但均呈现出下侧枝晶的生长受到促进,上侧枝晶的生长普遍受到抑制。

图3 择优生长角θ0=90°,t=0.1 s 时枝晶形貌及氧浓度分布Fig.3 Dendritic morphology and oxygen concentration at t=0.1 s with prefered growth orientation θ0=90°

图4 择优生长角θ0=0°,t=0.1 s 时枝晶形貌及氧浓度分布Fig.4 Dendritic morphologies and oxygen concentration at distribution t=0.1 s with prefered growth orientation θ0=0°

由于枝晶生长过程中在界面处析出溶质和释放潜热,导致其附近流体密度与较远处流体密度有较大区别,在重力的作用下引起自然对流,左右两侧形成对称分布的两个顺时针方向的旋涡。局部形成了小的旋涡,在小漩涡处出现了贫氧区。本研究主要考虑浓度场的变化引起的成分过冷度不同,从图3、图4所示的自然对流条件下浓度分布可以看出,自然对流促使溶质再分配出现局部差异。由式(9)可知成分过冷与界面浓度有关,各生长方向界面浓度不一致,因而导致枝晶生长速度上的差异。在自然对流的作用下,下侧枝晶处于迎流位置,析出溶质及凝固潜热易被流场带到枝晶上侧,从而使流场上游的枝晶尖端界面处浓度更低,过冷度更大,促进了该部分枝晶的生长。相反地,流场下游即枝晶上侧尖端,溶质富集,扩散层较厚,溶质不易向四周扩散,生长受到抑制。该结果与对流条件下的枝晶生长结论[14]一致。

进一步分析自然对流对浓度场的影响,图3、图4可以看出自然对流更利于溶质的扩散,使溶质分布更均匀。图5 为择优生长方向为0°,t=0.1 s 时,模拟区域中线(y=75 μm)处氧浓度的变化曲线。从图5可以看出,固空颗粒呈现出外表富氧,中心贫氧的特性,与刘海生等[3]宏观试验结果一致。与刘海生等人的试验相比,本文的模拟给出了微观的固空枝晶结构及更为精确的浓度分布结果。在t=0.1 s 时,纯扩散条件下的氧浓度峰值为23.71%,自然对流条件下的浓度峰值为23.14%且高氧浓度区域更少,LTTCHFIELD 等[15]通过枪击试验所得的氧浓度安全阈值为24%,可见自然对流将降低液氢系统中固空带来的安全风险。

图5 t=0.1 s 时y=75 μm 处不同条件下氧浓度分布Fig.5 Oxygen concentration distribution under different conditions with y=75 μm and at t=0.1 s

3.2 自然对流强度对固空晶体的影响

衡量对流强度的两个无量纲参数是浓度瑞利数和温度瑞利数,由于在枝晶生长过程中浓度瑞利数的影响较小,因此本文仅分析温度瑞利数对固空枝晶的影响。图6 给出了不同对流强度下枝晶生长情况。当温度瑞利数Ra从3 000 增加至7 000,自然对流对下侧枝晶生长的促进越来越明显。t=0.1 s 时,下侧枝晶的长度分别为:36.5 μm,41 μm,48 μm。同时,上侧枝晶臂随着流场强度的增加逐渐消失。

图6 t=0.1 s 时不同对流强度下枝晶形貌及氧浓度分布Fig.6 Dendritic morphology and oxygen concentration distribution under different convection intensities at t=0.1 s

虽然上侧枝晶臂在流场的作用下生长受到抑制,从图7 不同瑞利数下固相率随时间变化曲线可以看出,自然对流对枝晶生长的固相率净增长仍起到促进作用。随着自然对流强度的增加,固相率增长更快。在初始阶段,由于流场较弱,枝晶生长与纯扩散的条件差异较小,不同条件下固相率增速无较大区别。随着溶质的不断析出与潜热释放,区域的浓度差及温度差差异增大,自然对流强度增强,流场速度更快,自然对流对枝晶生长的影响越来越明显。自然对流使模拟区域的溶质分布更均匀,局部氧高浓度区减小,枝晶界面处的浓度扩散层更薄,促进了枝晶的生长。

图7 不同对流强度下固相率随时间变化趋势Fig.7 Solid-phase fraction changes with time under different convection intensities

从图8 不同条件下氧浓度峰值随时间变化也可以看出,纯扩散条件下的氧浓度峰值在0.02 s 后均高于存在自然对流,进一步说明自然对流能促进溶质的扩散,促进固相率的增长。注意到Ra=3 000 时,氧浓度峰值在t=0.01 s 时出现极大值达到22.58%,由于在凝固初期自然对流强度较小,溶质未充分向远处低浓度区扩散,导致溶质富集在枝晶上侧,随着自然对流强度增强,该峰值浓度下降。Ra=5 000 时凝固后期溶质峰值高于Ra=3 000,这是因为自然对流提高了凝固速率,更多的溶质析出,然而其自然对流强度不及Ra=7 000 时,不能及时向周围传递,导致浓度逐渐超过Ra=3 000 时的峰值浓度。

图8 不同对流强度下浓度峰值随时间变化Fig.8 Peak concentration changes with time under different convection intensities

3.3 多晶粒的影响

在中心晶粒四周放置4 个晶核,研究多晶粒生长时各晶粒之间的相互影响,模拟结果如图9 所示。凝固初始阶段,晶核均对称生长,如图9a。随着凝固的进行,自然对流强度增强,在枝晶附近形成了大小不等的旋涡,溶质扩散更充分,氧浓度在整个计算区域内分布更均匀,各枝晶的生长逐渐发展。图9b、图9c可以看出,随着枝晶生长,不断有溶质析出,中心处枝晶析出的溶质由于被周围晶粒阻挡,难以向周围扩散,同时又由于周围枝晶排出的溶质也向中心聚集,在中心区域形成了溶质富集区,抑制了中心晶粒的生长,向柱状晶转变。外部的枝晶生长受自然对流的影响较为显著,下侧枝晶的生长速度明显快于上侧枝晶,即自然对流将下游析出溶质带到上游,促使了下游枝晶的生长。当t=0.15 s 时,枝晶的增长使得流动区域更为狭窄,形成了更多的小漩涡,加速了溶质的扩散,浓度分布更为均匀。

图9 多个固空颗粒在自然对流作用下枝晶形貌随时间变化Fig.9 Dendrite morphology of several solid particles changes with time under natural convection

4 结 论

(1)建立改进的LBM-CA模型,减小网格各向异性影响,模拟自然对流条件下液氢中固空六重对称枝晶生长过程。结果进一步验证了固空颗粒外表富氧中心贫氧的特点,自然对流下固空枝晶生长不再呈现出对称性,上游枝晶生长得到促进,下游受溶质聚集的影响枝晶生长受到抑制。

(2)自然对流有利于析出溶质在区域内分布更均匀,氧浓度峰值较纯扩散条件有所下降,这有利于降低液氢系统的安全风险。自然对流条件下枝晶界面浓度扩散层更薄,促进了凝固过程,且自然对流强度越强,对固空的生长更有利。

(3)对多晶粒在自然对流条件下的生长模拟结果表明,受周围枝晶的影响,中心颗粒的枝晶生长受到明显抑制。