横向磁场下侧壁加热方腔熔化的数值模拟研究1)

孙思睿 张 杰 倪明玖

* (中国科学院大学工程科学学院,北京 100049)

† (西安交通大学航天学院,机械结构强度与振动国家重点实验室,西安 710049)

引言

磁场由于其非接触性以及对流动的调控能力,被广泛应用于包含导电材料的固-液相变过程中,如电磁冶金[1-2]、晶体生长[3-4]、增材制造[5-6]等等.而方腔熔化模型作为研究磁场作用下的固-液相变问题的基础问题,通常作为研究这类问题的一种重要方式.研究磁场作用下的方腔熔化问题,对相关工程实际应用具有重要的指导意义.

首先,热源作用下的方腔熔化问题作为一种典型工况,是探究相变材料熔化现象的一种有效手段,已经有学者开展了众多实验和数值模拟方面的研究[7-11].其中,Zhu 等[12]通过数值模拟计算了侧壁加热的方腔熔化问题,发现了在三维和二维情形下固-液界面存在显著区别.Hu 等[13]通过尺度分析和数值模拟验证,发现只有Pr≪1 且长宽比远小于1 的情况下,三维方腔熔化问题才可以简化为二维问题.并且文中给出的熔化体积分数标度关系对二维和三维结果给出了良好的预测.

近年来,由于电磁冶金行业的发展,各种类型的电磁场在冶金工艺中的应用日益广泛.磁场对液态金属的凝固/熔化的影响得到了越来越多的关注.对于外加磁场下的方腔熔化问题,Dulikravich 等[14]考虑了三维不可压缩黏性流体在磁场和不同重力条件下的凝固和熔化过程,发现磁场的强度和方向会显著影响流动的速度和涡量.Veilleux 等[15]和Zhang等[16]通过数值模拟和实验研究了低重力环境下的侧壁加热方腔熔化问题.他们发现磁场可以用来模拟在实际的低重力环境中发现的关键熔化特征.文献[17]通过数值模拟研究了在低温侧壁面作用下三维方腔内的自然对流和凝固过程.发现不同方向磁场对流动和传热的影响存在显著区别,尤其是平行温度梯度方向的磁场对流动的制稳作用最强.Feng 等[18]采用LBM 方法计算了沿平面内不同方向磁场作用下底面加热二维方腔内的固态镓(Ga)的熔化.研究发现小用于描述洛伦兹力和黏性力的比值,其中B0表示磁场强度,H表示方腔尺寸,σe表示金属电导率,μ 表示流体黏性),竖直磁场和水平磁场对熔化过程的阻碍作用相似;而在大Ha下,在熔化前期竖直磁场的阻碍作用更强而在熔化后期水平磁场的阻碍作用更强.Ghalambaz 等[19]从实验和数值角度详细讨论了洛伦兹力对倾斜方腔内熔化流体流动的抑制作用.文献[20]考虑了方腔底面中心方形截面条状热源加热熔化的模型,讨论了大磁场强度下洛伦兹力对流动及传热的抑制作用.以上的讨论主要关注均匀磁场对固-液相变的影响,对于其他类型磁场(如移动磁场[21-22]、四极磁场[23]、非均匀磁场[24]等)也有相应的文献研究.

此外,对于不同方向磁场的影响,在MHD(magnetohydrodynamics)自然对流中也有许多文献进行过研究分析.文献[25]通过实验测量了不同加热速率和不同方向、大小磁场作用下熔融金属镓的自然对流传热效率.揭示了垂直于环流方向的磁场对流动的阻滞作用明显弱于另外两个方向.文献[26-27]分别进行了不同磁场方向下的三维自然对流数值模拟并发现了施加垂直环流磁场能够提高传热效率.同时发现此方向磁场对速度的阻碍作用呈Gr/Ha用于描述自然对流浮力和黏性力之比,其中g为重力加速度,β 为热膨胀系数,ΔT为冷热源温差,κ 为热扩散系数,ν 为运动黏性系数.)的量级而非其他情形下的Gr/Ha2的量级.Chen 等[28]在此基础上,运用文献[29]提出的准二维模型与三维模型计算结果进行对比,发现垂直环流方向磁场作用下,在小Ha时二次流动因受洛伦兹力的抑制而减小,从而使得主流区的速度增强;在大Ha下三维模型与准二维模型的计算结果接近一致,并且通过准二维模型计算验证了磁场对流动和传热的阻碍作用的标度关系.

本工作研究垂直主环流方向的横向磁场影响,并通过加热左侧壁面来驱动方腔熔化产生的自然对流.本文仅考虑Rayleigh 数Ra≤106的情况表示热传导和热对流之间的比值),此时方腔内存在流动不稳定性但仍未出现湍流,故壁面湍流的复杂情形[30]不在本文的讨论范围内.另外熔化过程中可能存在杂质颗粒在流体中的流动,类似于椭球的低速流动[31],由于其对流动整体影响较小故本文亦不做考虑.由于磁场的存在,熔融金属流体的环形流动将产生感应电流和洛伦兹力,使得主流区的速度受到抑制并沿磁场方向变得更加均匀.然而,当磁场足够强时,垂直磁场方向的两侧壁面附近的剪切层具有更大的速度梯度,一般将此剪切层称为Hartmann 层.在大Ha下熔融液态金属的自然对流三维数值模拟中,若要精确解析Hartmann层内的电流分布,需要非常密集的网格和大量的计算资源,此时采用准二维模型进行等效模拟可以大大提高计算效率.本文对于固-液相变的模拟采用焓方法处理[32]而对于感应电势采用相容守恒方法处理[33].通过三维数值模拟发现小Ha下增大磁场会使得自然对流增强和熔化速率加快,同时流动及界面形状趋向于准二维.此外,本文采用准二维模型分析了大Ha下,磁场大小和自然对流强度对方腔内的流动形态、熔化速率以及固-液界面形状的影响.并且通过尺度分析得到了大Ha下,熔化速率和垂直最大流速随不同Ra和Ha的变化关系.

1 数值模型

1.1 控制方程

本文所研究的横向磁场作用下侧壁加热熔化的三维方腔计算模型如图1 所示.立方体方腔边长为H,各壁面均为绝缘壁面.左侧高温壁面温度为Th,方腔内固相初始温度和右侧壁面温度均为Ti略低于熔点温度Tm.重力方向竖直向下,均匀磁场B0垂直于主环流且方向沿z轴正向.

图1 物理模型示意图Fig.1 Sketch of the physical model

本文采用液态金属镓(Ga)作为工质,材料属性参数如表1 所示.该液态金属流体可视为单一组分的牛顿流体,其流动状态为非定常不可压缩层流流动.由于金属镓的固-液两相在熔点温度附近的密度差很小(见表1,约为3%),因此除浮力项采用Boussinesq 假设来表征密度差引起的对流外,方程中其他位置中的密度均视为常数 ρl,数值模拟过程中不考虑相变过程的体积膨胀.注意此简化在前人的数值模拟中已广泛采用[15,18].由于液态金属流动的低磁雷诺数属性(Rem=μmσeLu≪1 ,μm表示导磁率),导电流体流动所产生的感应磁场亦无需考虑.和壁面热源相比,由焦耳耗散所产生的热量可以忽略不计.综上,本文中采用的无量纲形式的控制方程组如下

表1 物性参数表Table 1 Material physical properties

其中,u,j,Φ,θ 分别表示无量纲形式的速度,电流密度,电势和温度.方程中各变量对应的特征尺度分别为长度H,质量 ρH3,速度 κ/H,时间H2/κ,磁场强度B0,电势 κB0,电流密度 σeκB0/H,相应地,各无量纲数的定义分别为: Fourier 数用于描述非稳态导热过程的无量纲时间,Prandtl 数用于描述动量和热扩散系数之比,Stefan 数用于描述显热变化量与相变潜热的比值,Hartmann数Ha 如前面定义.

需要说明的是,对于相变过程,本文将采用焓方法进行模拟,将固-液两相视为单一相流体共同计算温度场和速度场.对于温度场,利用工质内总比焓与显热和相变(熔化)潜热的关系式h=cpΔT+ΔH,其中 ΔH=f L,f为液相所占的体积分数,L为相变潜热.并将该关系式代入有量纲形式的能量方程,可将方程改写为

其中,对流项中关于相变潜热的分项只在界面附近有效,且界面附近速度值很小固可将该项忽略.而对于速度场,则采用Carman-Kozeny 关系式对固相以及混合网格施加一个等效黏性,使得速度趋近于0[34].其中C为Carman-Kozeny 常数,其值取为 106,ε 为防止分母除0 的计算常数,取为 10-3.

1.2 准二维模型

根据文献[29]得到的准二维模型,可以得到考虑固液相变的控制方程组的无量纲形式如下

上述方程中,下标 ⊥ 表示该方程的求解域为主流区中任一垂直于磁场方向的截面.通常来说,当Ha≫1时,Hartmann 层的厚度H/Ha显著小于黏性边界层厚度,层内的黏性耗散也由于速度梯度的增大而增强并可将其等效为对主流区流动的线性阻力(如式(8)的最后一项所示).同时,主流区的流场沿磁场方向趋于均匀,Ha≫1 流动由三维状态转变为准二维状态.

1.3 边界条件

对于三维和准二维模型,各壁面速度均采用无滑移无穿透边界条件.左右两侧壁面温度采用Dirichlet边界条件,其余壁面温度采用绝热边界条件.各壁面的电势采用绝缘边界条件.横向均匀磁场产生的洛伦兹力项作为体积力作用在全场.对于三维模型,洛伦兹力的计算需要通过先行求解关于 Φ 的电势泊松方程来得到.而对于准二维模型,此时洛伦兹力项作为线性项可直接代入动量方程计算.

2 数值算法及验证

本文中的数值模拟研究基于开源自适应计算平台Basilisk[35],N-S 方程的计算采用两步近似投影法[36],电势泊松方程的计算采用相容守恒格式.Hartmann层厚度为H/Ha且在本文计算中至少包含了6 个网格,以精确解析电流边界层.通过网格无关性验证,可以发现为使液相体积分数趋于收敛,主流区网格尺寸应满足 Δr≤H/128,固-液界面附近网格尺寸应满足 Δr≤H/256.

2.1 二维/三维无磁场侧壁加热方腔熔化验证

本文首先验证了固-液相变过程模拟的精确性.图2 比较了三维矩形腔中,左侧壁面加热熔化作用下,不同时刻的固-液界面形状以及位置.该工况中矩形腔的长高比为D/H=1.4,宽高比为W/H=0.6,各无量纲参数为Pr=0.024 4,Ra=6.6×105,以及S te=0.04.本文模拟的三维结果与文献[7]的实验测量结果和文献[12]的三维模拟结果符合良好.同时还针对二维模拟结果与文献[12]的结果进行了比较,发现由三维和二维模型计算得到的界面形状存在显著差异,这种差异主要是由三维情形下沿第三个方向产生的二次流动造成的.

图2 固-液界面形状及位置比较.空心点为文献[7]的实验结果,虚线为文献[12]的数值模拟结果,实线为本文的计算结果Fig.2 Profile comparison of the melting fronts.The hollow dots are measured by Ref.[7] experimentally,the dashed lines are the numerical results by Ref.[12],while the solid lines are the present numerical solutions

2.2 三维MHD 方腔自然对流验证

随后,本文验证了磁场作用下三维方腔内的自然流动问题.图3 展示了Ra=1×105时,不同Ha下的平均竖直速度在靠近加热壁面中线(x=0.05,y=0.5)的分布情况,各Ha下的结果均与文献[28]的模拟结果符合.

图3 Ra=1×105,不同 Ha 下位于 x=0.05,y=0.5 线段处的平均竖直速度分布图.点状线为文献[28]的计算结果,实线为本文的计算结果Fig.3 The distribution of mean vertical velocity along x=0.05,y=0.5for Ra=1×105 with different Ha.Dotted lines are the results simulated by Ref.[28],solid lines are the present numerical results

此外,本文还比较了Ra=1×105时,三维和准二维情形下不同Ha下的时空平均Nu,如表2 所示发现各个工况均与文献[28]的计算结果吻合良好.由此充分证明了本文计算模型和算法的准确性.

表2 不同 Ha 下的 与文献[28]的计算结果比较Table 2 Comparison of against the numerical results provided by Ref.[28] at different Ha

表2 不同 Ha 下的 与文献[28]的计算结果比较Table 2 Comparison of against the numerical results provided by Ref.[28] at different Ha

3 结果分析与讨论

3.1 小 Ha 下的三维侧壁加热方腔熔化

文献[26-27]分别进行了不同磁场方向下的三维自然对流数值模拟并发现了施加垂直环流磁场能够提高传热效率.这是由于小Ha下洛伦兹力的整流作用造成的.在此Ha区间,平行磁场方向的二次流动受到了抑制使得流动被迫向主环流发展,同时Hartmann层中的黏性耗散没有明显增强,流动仍然由浮升力主导.本节中三维情形下的计算均以Ra=105,S te=0.05时的工况进行分析讨论.

图4 展示了Fo=1.6,在靠近加热壁面中线(x=0.05,y=0.5)的竖直速度和竖直洛伦兹力分布.当Ha=0时,流动结构受到三维效应的影响变得混乱,主环流的速度存在较大波动.当Ha=50 时,靠近壁面位置的速度由于洛伦兹力的作用有所增大,主流区的洛伦兹力接近于0,此时流动受到的洛伦兹力在Hartmann 层内的驱动作用大于在主流区中的焦耳耗散,使得速度整体增大并且趋于均匀.

图4 Fo=1.6,x=0.05 ,y=0.5 处的瞬时竖直速度 uy (黑线)和洛伦兹力分量 Fly (蓝线)的分布Fig.4 The distribution of the instantaneous vertical velocity (black lines) and y-component Lorentz force (blue line) at Fo=1.6,x=0.05 and y=0.5

图5 对比展示了Fo=1.6时,Ha=0 (z≤0.5) 和Ha=50(z≥0.5) 两工况中,靠近壁面热源的平面(x=0.05)内的瞬时竖直速度云图和平面速度矢量.从平面速度矢量的变化可以看出Ha=50 时横向的二次流动明显受到了抑制,从瞬时竖直速度云图可以发现Ha=50 时主环流中的最大速度明显增大且分布更均匀.这是由于洛伦兹力对Hartmann 层流动的驱动使得液相的流动更多地被导向主环流方向从而使得横向二次流动减少,因此整体上表现出对流场的整流效果.为了确定整个流场中速度的变化情况,本文进一步统计了全场总动能随时间的变化情况,如图6 所示.可以看到各工况下总动能随时间均存在一定波动,而随着液相区域的增大,Ha=50 工况下的总动能明显大于Ha=0 工况下的总动能,这说明小Ha下洛伦兹力的作用以Hartmann 层内的驱动力为主导,对流体做正功,使得流动速度整体增大.

图5 Fo=1.6,x=0.05 处的瞬时竖直速度云图,分别对应Ha=0(z ≤0.5)和 Ha=50 (z ≥0.5).图中矢量表示 z-y 平面内的速度Fig.5 The contour maps of the instantaneous vertical velocity (uy) for Ha=0(z ≤0.5) and Ha=50 (z ≥0.5) at x=0.05.The vectors present the velocity in the z-y plane

图6 Ha=0 和 Ha=50 时,全场总动能随时间的变化情况Fig.6 Evolution of total kinetic energy in the whole domain for Ha=0andHa=50

为了确认小Ha下磁场对传热的增强作用,比较了不同Ha下液相体积分数和Nu随时间的变化.图7展示了Ha=0,Ha=25 和Ha=50 时液相体积分数随时间的变化,进一步验证了小Ha下洛伦兹力对流场结构的整流作用.主流区速度受的洛伦兹力的抑制作用较小,同时由于剪切边界层内洛伦兹力做正功加快了附近流体的速度从而使得主流区增大,而因此出现传热和熔化增强的现象.

图7 不同 Ha 下,(a)液相体积分数随时间的变化和(b)Nu随时间的变化Fig.7 (a) Time evolution of the liquid volume fraction and (b) time evolution of Nu at differentHa

图8 通过流线图对比Ha=0 和Ha=50 时流线和固-液界面形状,发现随着磁场增大,流动逐渐趋于稳定和二维化,同时,从界面上云图的x轴坐标值可以看到界面形状也趋于二维化.因此采用准二维模型模拟大横向磁场下自然对流熔化问题是有效且可靠的.

图8 Fo=1.6,界面形状及流线图比较,界面上云图的颜色表示x轴坐标值Fig.8 Comparison of the solid-liquid interface and the streamlines at Fo=1.6The contours on the interface are the x coordinate

3.2 大 Ha 下的准二维侧壁加热方腔熔化

由上文可知,随着Ha的增大,方腔内熔融流体的流动形态逐渐从混乱(不稳定)的三维流动转变为稳定的准二维环流.同时,固-液界面也从中间凸起的三维界面转变为准二维界面.运用提出的准二维模型,本节分别对S te=0.05 时,1 04≤Ra≤106和100 ≤Ha≤6400范围内的方腔熔化自然对流过程中流动、传热和界面形状进行了分析.

图9 展示了Ha=400 时,不同时刻和Ra下界面的形状、流场和温度场.可以看到Fo=0.2 时,熔融流体区域所占宽度只有方腔尺寸的15%,此时流场对应的瞬时Ra很小,传热和熔化受热传导主导.而其中Ra=1×106时已经出现4 个涡流结构使得对流换热的作用开始增强.而对于Fo=1 时,Ra较小的两个工况中均出现单一主环流,而Ra=1 ×106中的流动由于受到浮力的驱动作用更强仍然表现出主环流和多个小环流共存的情况.当Fo=3 时,对流作用对界面形状的影响更为明显.

图9 Ha= 400 时,不同 Ra 下的温度云图及流线图比较(从左至右依次为 Ra=1×104,Ra=1×105, Ra=1×106)Fig.9 Comparison of the temperature contour and the streamlines at different instants (Ra=1×104,Ra=1×105, Ra=1×106 from left to right) at Ha= 400

图9 Ha=400 时,不同Ra 下的温度云图及流线图比较(从左至右依次为 Ra=1×104,Ra=1×105, Ra=1×106)(续)Fig.9 Comparison of the temperature contour and the streamlines at different instants (Ra=1×104,Ra=1×105, Ra=1×106 from left to right) at Ha= 400 (continued)

图10 展示了Ra=1×106时,固-液界面随Ha增大的变化趋势,所取时刻均为Fo=2.可以看到当Ha逐渐增大时,由于自然对流受到的抑制增强从而使得固-液界面形态逐渐由大曲率界面向平直界面转化.

图10 Ra=1×106,Fo=2 固-液界面随 Ha 增大的变化趋势Fig.10 Comparison of the solid-liquid interfaces at Ra=1×106 and Fo=2when Ha is varied

图11 展示了在不同Fo,Ha和Ra下,熔融金属液相体积分数f和传热Nu的演化情况.如图11(a)所示,磁场强度对熔化速度的影响与Ra有关.当Ra=1×104时,液相在熔化过程中受到的浮升力过小,熔化界面随时间的变化与一维热传导的解析解(图中空心点)的结果基本一致,说明此时熔化完全由热传导主导.Hu 等[13]通过尺度分析得到热传导占主导时的标度关系为fl～(S teFo)12,对流换热占主导时的标度关系fl～S teFo.因此当S te为常数时,fl与Fo之间仍然应当满足上述标度关系.当Ra=1×105和Ra=1×106时,由于液相中对流作用的增强使得在Ha较小时仍能在对流阶段保持液相体积分数与Fo的一次方关系,而随着Ha逐渐增大,凝固速率由于洛伦兹力对主环流的阻滞作用逐渐下降,并分别在Ha=3200 和Ha=12 800 时回复到由热传导主导的液相体积分数与Fo的1/2 次方关系.

图11(b)展示了不同Ra下,磁场强度对Nu的影响.当Ra=1×104时,不同Ha下Nu基本没有变化,当Ra=1×105和Ra=1×106时,由于自然对流随Ha而减弱,使得Nu逐渐减小.值得注意的是当Ra=1×106时,尽管流动已趋于二维,在小Ha下平面内的强烈且不稳定的自然对流仍然会引起壁面处Nu的震荡.

图11 不同 Ha 下,(a)液相体积分数和(b)左侧壁面 Nu 随时间的变化,从左至右依次为 Ra=1×104,Ra=1×105, Ra=1×106Fig.11 Time evolution of (a) the liquid volume fraction and (b) the Nusselt number at different Hartmann number,while the sub-panels correspond to Ra=1×104,Ra=1×105 and Ra=1×106 from left to right

图11 不同 Ha 下,(a)液相体积分数和(b)左侧壁面 Nu 随时间的变化,从左至右依次为 Ra=1×104,Ra=1×105, Ra=1×106 (续)Fig.11 Time evolution of (a) the liquid volume fraction and (b) the Nusselt number at different Hartmann number,while the sub-panels correspond to Ra=1×104,Ra=1×105 and Ra=1×106 from left to right (continued)

如Tagawa 等[27]的文中所分析的,对于磁场方向垂直主环流方向的MHD 问题,侧层中的电流是用来表征洛伦兹力的阻滞作用的一个重要方式.并且由电流与速度的关系可以得知方腔侧壁剪切层内的最大速度应在Ha/Gr,即HaPr/Ra的量级.图12展示了熔化过程中熔融液体的最大竖向速度uy,max与无量纲参数组合HaPr/Ra的变化关系.可以看到在Ha足够大的情况下,本文计算中所考虑的不同Ha及Ra下各工况中的uy,max均满足同一标度关系.

图12 熔化过程中的最大速度(uy,max)与无量纲参数组合HaPr/Ra的关系Fig.12 The relation between the maximum velocity (uy,max) during the melting process and the combination of dimensionless numbers ofHaPr/Ra

4 结论

本文通过数值模拟的方法研究了垂直环流方向横向磁场对侧壁加热方腔内熔化和自然对流的影响,并分别对小Ha和大Ha两个区间的工况采用三维模型和准二维模型进行计算和分析,得到了以下主要结论.

本文采用三维数值模拟分析了小Ha下的情形.在该区间下的磁场强度较小,洛伦兹力对流动的影响有限,其对流动的影响主要表现为对流动结构的整流作用,而熔化和流动过程仍然受浮升力主导.随着磁场增大,二次流动受到抑制,且主环流速度有所增大,从而出现熔化和传热速率增大的情况,同时熔化过程中固-液界面的存在也使得流动形态进一步趋向于二维化,而二维的流动形态也反过来决定了固-液界面也将趋于二维化.

对于大Ha下的情形,随着磁场增大,方腔内的流动形态和界面形状均趋于二维化.经验证在此区间准二维模型的计算结果与三维模型非常接近,采用准二维模型计算能在保证准确的前提下大大缩短计算时间.可以发现不同Ra下的对流强度都受到洛伦兹力的抑制而显著减小,固-液界面形态也趋于平直.液相的熔化速率和壁面处的传热效率都呈单调下降趋势.而熔化过程中的最大速度仍然满足Tagawa 等[27]提出的无量纲参数组合HaPr/Ra之间的标度关系.