作为最轻的金属结构材料，镁合金因其优异的物理特性、较低的成本、出色的力学性能及吸能效果，在汽车和航空航天等领域有着广阔的应用前景，更是享有“21 世纪的绿色结构材料”美誉[1-2]。当前，镁合金零部件主要通过铸造工艺生产，例如座椅框架、方向盘轴等[3-4]。镁合金铸造成形工艺主要有砂型铸造、金属型铸造、挤压铸造、压铸等，无论采用哪种成形工艺，气孔缺陷都难以避免[5-6]。气孔来源于未从金属液中逸出的气泡，如果凝固速度大于气泡的外逸速度，气泡就会残留在金属中形成气孔[7]。按照形成原因，气孔可分为卷入性气孔、侵入性气孔、反应性气孔和析出性气孔等[8]。气孔会对铸件的性能产生负面影响，如减少材料的有效承载面积、产生应力集中、恶化力学性能等，从而影响铸件的使用寿命。
在实际工况中，气泡的动力学行为极其复杂，不仅受浮力、表面张力等的作用，其形貌还会受到熔体中杂质颗粒、固体组织等障碍物的影响，如凝固过程中的枝晶组织会与气泡发生复杂的相互作用[9-10]。因此，研究铸造过程熔体中的气泡动力学行为对减少和消除气孔缺陷至关重要[11]。然而，在镁合金零部件的铸造过程中，熔体温度高、流速大，再加上铸造模具不透明、合金易氧化等原因，采用实验手段对铸造过程中的气泡动力学行为进行直接的观察和研究存在较大的困难[12]。
通过数值模拟，将气泡在镁合金熔体中的动力学行为抽象为气泡在含障碍物微通道内的动力学行为，将凝固过程中的枝晶组织抽象为不同尺寸和形状的障碍物，可以对气泡和障碍物间的相互作用进行定量化研究。本文在对镁合金气孔缺陷实验表征的基础上，借助数值模拟技术，将熔体中的气液两相流行为抽象化，通过在计算域中设置可变宽度的障碍物，研究障碍物尺寸对气泡动力学的影响，为研究镁合金熔体中的气泡动力学行为及气孔缺陷的控制提供指导。研究成果对提高镁合金铸件性能和扩大镁合金在汽车等行业中的应用具有重要意义。
1 实验和模型方法
1.1 实验方法
本文选用铸态Mg-30%Al(质量分数)合金作为研究对象。实验中所用的原材料为高纯Mg(99.99%)和纯铝(99.9%)。合金在RPM-20-9 坩埚熔化电阻炉中进行熔炼，炉膛预热温度为500 ℃，熔炼时通入CO2∶SF6 为99∶1 的混合保护气体，熔炼结束后采用圆柱形铸造模具进行浇注，模具内壁均匀涂抹氮化硼-无水乙醇悬浮液，模具预热温度为200 ℃。
将制备的圆柱铸锭沿中轴线进行取样，样品尺寸为10 mm×10 mm×5 mm。使用碳化硅砂纸对试样进行水磨，研磨完成后，先后用蒸馏水和无水乙醇冲洗试样，再用风干机吹干备用。采用金相腐蚀液(2 mL 冰乙酸+10 mL 无水乙醇+过饱和苦味酸)对试样表面进行腐蚀，并用无水乙醇冲洗合金试样表面的腐蚀产物，吹干后进行显微组织观察。使用JEOL JSM-7800F 场发射扫描电子显微镜(SEM)对铸态合金的显微组织和气孔的分布、形貌等进行观察。合金实际成分采用电感耦合等离子体发射光谱仪(ICP-OES)进行测定。
1.2 相场-格子玻尔兹曼模型
1.2.1 保守相场模型
引入非守恒序参量，即相场变量φ，刻画体系状态。当φ=0 时，表示液相；当φ=1 时，表示气相；当φ在0～1 之间取值时，表示气液两相弥散界面。通过获得每个时刻φ 的分布来确定体系的演变。φ 的控制方程为[13]：
式中，M 是迁移率；δ 是界面厚度参数；ν 是流速矢量。
在平衡状态下，沿界面法线方向，φ 的变化用双曲正切来描述，即：
式中，ε 是沿界面法线方向的坐标。
1.2.2 格子玻尔兹曼模型
当气泡接触障碍物时，会出现较大的变形，如挤压、分裂等，影响数值求解的稳定性[14]。考虑到相界面本质上是介观的，而格子玻尔兹曼模型也是介观的，因此采用基于介观动力学的格子玻尔兹曼方法计算流体速度。相比基于连续介质的Navier-Stokes方程，格子玻尔兹曼模型假设真实的流体由一系列粒子组成，稳定性高，易于处理复杂边界。基于Bhatnagar-Gross-Krook(BGK)碰撞模型[15]，分布函数的演化描述如下：
式中，fi(r,t)是分布函数；τ 是弛豫时间，它与运动黏度υkin 有关：
(r,t)是平衡分布函数，表示为：
式中，ρ 是粒子密度；c=δx/δt=1 是晶格速度；δx 和δt分别是晶格间距和格子玻尔兹曼模型中的时间步长；wi 是权重系数。图1 为格子玻尔兹曼模型中的三维十九速度（D3Q19）离散速度模型，其中3 是指空间维数，19 是指离散速度模型中格子分速度方向的数目。因D3Q19 模型具有较高的精确度和数值稳定性，所以采用D3Q19 模型。对于D3Q19 模型，w0=1/3，w1-6=1/18，w7-18=1/36，沿各自方向的离散速度ci 见表1。
表1 D3Q19模型中ci的3个速度分量ci=(u, v, w)
Tab.1 Three velocity components of ci=(u, v, w) in the D3Q19 model
i 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 u 0 1 -1 0 0 0 0 1 -1 -1 1 1 -1 -1 1 0 0 0 0 ν 0 0 0 1 -1 0 0 1 1 -1 -1 0 0 0 0 1 -1 -1 1 w 0 0 0 0 0 1 -1 0 0 0 0 1 1 -1 -1 1 1 -1 -1
图1 D3Q19 离散速度模型
Fig.1 D3Q19 discrete velocity model
Fi(r,t)是离散外力：
式中，Fg 是驱动气泡运动的浮力(或惯性力)密度；g是沿y 轴负方向的重力加速度；Fs 是表面张力密度；μ 是化学势，它取决于φ 的梯度[16-17]，即：
式中，σ 是表面张力系数。
这2 种力揭示了气泡上升动量和张力效应间的竞争关系。较大的浮力(或惯性力)意味着更快的上升速度和更大的变形，而较大的表面张力对应于接近球形的形状，即变形较小。
考虑这2 个力的影响，流速的计算公式如下：
将相场法和格子玻尔兹曼模型相结合，建立了相场-格子玻尔兹曼模型，相场法用于刻画气液界面，格子玻尔兹曼模型用于刻画气液两相流，该方法结合了相场法的弥散界面特征和格子玻尔兹曼模型的稳定性和边界易处理等特点，提高了气泡动力学行为模拟的稳定性。
1.2.3 离散化
为了减小离散误差，式(1)在空间上采用显式有限差分法离散，在时间上用四阶Runge-Kutta 法离散。在四阶Runge-Kutta 法中，相场值的增量由4 个估计斜率的加权平均值确定，其权重分别为1/6，2/6，2/6，1/6，即：
式中，φn 是φ(n)在n 时刻的Runge-Kutta 近似值，
当4 个估计斜率均等于k1 时，四阶Runge-Kutta法还原成一阶欧拉法。
为了保持高稳定性，避免可能导致数值不稳定的梯度急剧变化，式(16～17)中的散度和拉普拉斯算子采用D3Q19 晶格上的18 个最近邻点进行离散[18]。
式中，i，j 和k 分别是沿x，y 和z 轴的单位向量。
2 结果与讨论
2.1 实验结果分析
通过ICP-OES测得样品实际成分为Mg-30.37%Al(质量分数)，SEM 组织如图2 所示，气孔形貌相对比较规则，周边平滑，呈现圆形或近圆形。从分布上来看，气孔集中于α-Mg 和共晶组织的交界处，呈现单个或几个气孔同时存在的特征。此外，气孔还会与枝晶发生相互作用。气孔影响枝晶的生长，阻碍枝晶发展，枝晶会影响气体的运动，改变气孔形貌。
图2 铸态Mg-30%Al(质量分数)合金的SEM 组织图片
Fig.2 SEM images of the as-cast Mg-30%Al(mass fraction)alloys
通过对不同位置试样的组织形貌进行对比，发现气孔沿铸锭高度方向的分布具有一定规律性，在铸锭底部区域气孔数量较少，随着取样位置的升高，气孔数量增加。这一方面是由于位置相对靠上的熔体温度较低，气体的溶解度较小，另一方面是由于气泡形成后会上浮，容易在铸锭上部积聚。
2.2 数值模拟分析
2.2.1 初始条件和边界设置
图3 为三维数值模拟中障碍物的设置和气泡初始位置。球形气泡初始直径d0 为32 dx，中心位于(0.5X,0.5Y,0.2Z)，其中X=160 dx，Y=160 dx，Z=384 dx，表示矩形微通道的大小，dx 为网格尺寸。障碍物为长方体形状，中心位于(0.5X,0.5Y,0.5Z)，沿X-Z 中央纵截面截取障碍物，其二维截面中障碍物纵向尺寸为ver/d0=1，横向尺寸设置6 种宽度，hor/d0 分别为1/32、1/16、1/8、1/4、1/2 和1。无量纲参数设置如下：液相密度ρl=1.0，气相密度ρg=1.0×10-3，迁移率M=1.0×10-2，液体运动黏度υkin=1.500×10-1，表面张力σ=5.5×10-2，界面厚度δ=3.2，dx=0.8。对应的Reynolds、Eötvös 和Morton 3 个无量纲数的大小分别为34.12、14.88 和2.43×10-3，位于Bhaga 和Weber[19]提出的气泡形状图谱中的球形区域。
图3 初始化示意图
Fig.3 Initialization diagram
在计算域所有壁面和障碍物表面，对相场变量采用零诺伊曼边界条件，对速度场变量采用无滑移边界条件。流动发生在液相和气相中，气泡碰到障碍物表面和壁面后发生挤压变形。
2.2.2 不同障碍物设置下的气泡动力学
图4 为不同障碍物设置下的气泡演变过程模拟结果。气泡轮廓根据φ=0.5 提取。在惯性力、黏性力和表面张力的综合作用下，气泡沿中心对称轴自底向上运动，在hor/d0≤1/4 时，提取9 个时刻；当hor/d0＞1/4 时，因障碍物宽度较大，气泡到达障碍物附近后无法继续向上运动，提取4 个时刻。障碍物通过影响流场分布，改变了气泡动力学。当矩形障碍物宽度hor/d0≤1/4 时，气泡接触障碍物后，受障碍物切割作用发生分裂，形成两个子气泡继续上升，因左右两侧流场分布存在差异，气泡具有不同的上升速度，如图4(a～d)所示；当hor/d0＞1/4 时，气泡到达障碍物底部后，无法继续上升，发生挤压变形，如图4(e～f)所示。
图4 不同障碍物设置下的气泡演变过程
Fig.4 Bubble evolution under different obstacle settings
为了表征气泡形状的演变，引入三维形状因子FF3 和二维形状因子FF2 两个形状参数量化气泡形状与流场条件、障碍物几何尺寸间的关系。
三维形状因子FF3 定义为气泡实测体积与具有和气泡相同表面积的球的体积之比，表示气泡形状和球之间的相似性，取值范围(0,1]，三维形状因子FF3 越接近1，气泡越接近球形，其表达式为：
式中，V 为气泡的实测体积；S 为气泡的实测表面积。
二维形状因子FF2 定义为气泡在y=0.5Y 面上的横截面积与具有和气泡横截面相同周长的圆的面积之比，表示气泡横截面和圆之间的相似性，取值范围(0,1]，二维形状因子FF2 越接近1，气泡横截面越接近圆形，其表达式为：
式中，A 为气泡在y=0.5Y 面上的横截面积；P 为气泡在y=0.5Y 面上的横截周长。
图5 为在计算域内设置不同类型的障碍物时，气泡的三维形状因子FF3 和二维形状因子FF2 随时间的变化，其中t*=1 000 dt。在气泡接触到障碍物之前，形状因子-时间曲线重合；接触到障碍物后，曲线发生较大波动。当矩形障碍物宽度hor/d0＞1/4时，气泡碰到障碍物后发生分裂，分裂后的气泡表面积、横截周长等发生较大变化，在t/t*=1.1 时，FF3和FF2 出现明显下降。分裂后新形成的两个子气泡绕过障碍物后，继续上升，形貌变化较小，FF3 和FF2 出现微小波动。当hor/d0>1/4 时，气泡受到障碍物阻挡无法继续上浮，在障碍物底部发生挤压变形，FF3 和FF2 在1.1≤t/t*≤2.2 时出现明显波动；在后续过程中，气泡形状维持近球形，FF3 和FF2 趋于稳定；由于气泡没有发生分裂，气泡的形状因子比hor/d0≤1/4 时大。
图5 不同障碍物设置下气泡的形状参数与时间的关系，t*=1 000 dt
Fig.5 Shape parameters vs time under different obstacle settings,t*=1 000 dt
图6 为形状参数FF3 和FF2 与矩形障碍物宽度hor 的关系。当矩形障碍物宽度hor 较小时，气泡受到障碍物的切割作用，形状发生较大变化。hor 越小，气泡碰到障碍物后发生分裂的倾向越明显，FF3和FF2 越小。当矩形障碍物宽度hor/d0＞1/4 时，气泡受到障碍物阻挡。hor 越大，阻挡作用越强，但受限于气泡直径和障碍物间的相对大小，气泡变形度并不是越来越大。以图4(e)为例，气泡虽被障碍物阻挡，但碰到障碍物后，会将障碍物包裹，产生比图4(f)更大的变形，FF3 和FF2 随hor 的增大而增大。
图6 形状参数与障碍物宽度的关系
Fig.6 Relationship between shape parameters and obstacle horizontal width
气泡在计算域顶部的变形也被考虑，通过和在障碍物附近的变形进行对比，FF3 和FF2 的值均变小，且在障碍物宽度变化时，其值基本没有变化，平均值分别为0.68 和0.48，表明气泡碰到计算域顶部后虽会发生进一步的变形，但该变形受障碍物宽度的影响较小。
微通道中的气泡动力学与障碍物的几何尺寸有很大关系。随着障碍物宽度的增大，障碍物对气泡运动的影响由切割作用变为阻挡作用。结合图4 中的气泡轮廓图及图5～6 中形状参数的变化，障碍物的切割作用会导致更大的气泡变形。这为研究凝固过程中气泡与凝固组织间的相互作用提供了指导。以凝固过程为例，枝晶尖端半径对气泡的上升过程会产生明显影响，通过调控凝固条件，改变枝晶尺寸，将改变气泡的分裂和上升行为，从而影响气泡从熔体中的逸出过程。
3 结论
(1)铸态Mg-Al 合金中的气孔缺陷集中于α-Mg和共晶组织的交界处，其表面较为光滑，形貌一般呈圆形或近圆形。气孔会与枝晶相互作用，形貌发生变化，影响枝晶生长。
(2)气泡在微通道内运动时，受障碍物影响，形状发生变化，障碍物切割作用导致的气泡形状改变明显大于其阻挡作用。
(3)当矩形障碍物宽度较小时，受障碍物切割作用，气泡会发生分裂，障碍物宽度越小，发生分裂的倾向越明显。当矩形障碍物宽度较大时，障碍物的影响会变为阻挡作用，宽度越大，对气泡运动的阻挡作用越大。
[1]潘复生，吴国华.新型合金材料：镁合金[M].北京：中国铁道出版社，2017.
[2]SONG J F, CHEN J, XIONG X M, et al.Research advances of magnesium and magnesium alloys worldwide in 2021[J].Journal of Magnesium and Alloys,2022,10(4):863-898.
[3]肖旅，侯正全，吴国华，等.高强韧稀土镁合金大型复杂铸件制造技术研究现状及展望[J].特种铸造及有色合金，2021,41(7):793-801.
[4]薛斌，许忠斌，张小岩，等.轻量化精密铸造成型技术在航空航天关键部件中的应用[J].铸造技术，2022,43(4):290-294.
[5]张占领，张艳琴，刘真，等.镁合金压铸件常见缺陷及改进措施[J].铸造技术，2019,40(7):718-721.
[6]蒋斌，刘文君，肖旅，等.航空航天用镁合金的研究进展[J].上海航天，2019,36(2):22-30.
[7]张昂，郭志鹏，蒋斌，等.合金凝固组织和气孔演变相场模拟研究进展[J].中国有色金属学报，2021,31(11):2976-3009.
[8]陈荣石，周波，李吉林，等.铸造高强耐热Mg-Y-Nd(-Gd)-Zr 和Mg-Gd-Y-Zr 系镁合金组织性能和铸造缺陷对比[J].铸造，2021,70(1):15-20.
[9]ZHANG A, GUO Z P, JIANG B, et al.Multiphase and multiphysics modeling of dendrite growth and gas porosity evolution during solidification[J].Acta Materialia,2021,214:117005.
[10]ZHANG A, DU J L, ZHANG X P, et al.Phase-field modeling of microstructure evolution in the presence of bubble during solidification[J].Metallurgical and Materials Transactions A,2020,51:1023-1037.
[11]喻兵，贾征，李又佳，等.镁合金熔体净化技术研究进展[J].铸造技术，2021,42(7):635-643,650.
[12]苏东泊.铸造镁铝合金气孔缺陷及气泡动力学模拟研究[D].重庆：重庆大学，2022.
[13]ZHANG A, DU J L, GUO Z P, et al.Conservative phase-field method with a parallel and adaptive-mesh-refinement technique for interface tracking[J].Physical Review E,2019,100(2):023305.
[14]ZHANG A,SU D B,LI C M,et al.Investigation of bubble dynamics in a micro-channel with obstacles using a conservative phase-field lattice Boltzmann method[J].Physics of Fluids,2022,34(4):043312.
[15]BHATNAGAR P L, GROSS E P, KROOK M, et al.A Model for collision processes in gases.I.small amplitude processes in charged and neutral one-component systems[J].Physical Review,1954,94(3):511-525.
[16]FAKHARI A, GEIER M, LEE T, et al.A mass-conserving lattice Boltzmann method with dynamic grid refinement for immiscible two-phase flows [J].Journal of Computational Physics,2016,315:434-457.
[17]ZHANG A, GUO Z P, WANG Q G, et al.Three-dimensional numerical simulation of bubble rising in viscous liquids:A conservative phase-field lattice-Boltzmann study[J].Physicsof Fluids,2019,31(6):063106.
[18]CHENG M, HUA J S, LOU J.Simulation of bubble-bubble interaction using a lattice Boltzmann method[J].Computers & Fluids,2010,39(2):260-270.
[19]BHAGA D, WEBER M E.Bubbles in viscous liquids: shapes,wakes and velocities[J].Journal of Fluid Mechanics, 1981, 105:61-85.