多孔介质模型镁合金流变铸轧LBM模拟

张 莹，马 强，陈杨华，倪泽联

（南昌大学机电工程学院，南昌330031）

多孔介质模型镁合金流变铸轧LBM模拟

张 莹，马 强，陈杨华，倪泽联

（南昌大学机电工程学院，南昌330031）

根据REV尺度多孔介质格子－Boltzmann方法基本理论，建立求解高固相率半固态浆料渗流过程的热流耦合数值模型．应用模型对变孔隙率Poiseuille流动及纯扩散凝固进行模拟．LBM模拟结果与有限差分模拟吻合．分析了压差作用下浆料在平板间渗流冷却与热流耦合过程．模拟结果表明浆料渗流速度随温度降低而减慢，相同时刻下凝固潜热较大的浆料渗流Re数高．利用提出的模型初步探索了AZ91D镁合金半固态流变铸轧过程的热流场分布．

格子－Boltzmann；多孔介质；半固态；热流耦合

作为半固态金属凝固过程的最后阶段，高固相率浆料内的传输行为很大程度上影响了材料凝固的最终质量．目前，有部分学者通过多孔介质模型去研究了高固相率半固态浆料的挤压成形过程中的应力和两相分离行为并取得成功［1－2］，但是基于多孔介质模型对半固态对流凝固过程的热流场分析尚没有相关的文献报道．从孔隙尺度上说，研究液态合金在固相晶粒组织间的流动传热传质过程，分析晶粒生长环境，是获得晶粒组织形态，控制凝固质量的基础［3］．从表征体元（REV）尺度上说，研究浆料的固相率和渗透率、物性参数对热流场的影响，揭示固相率的半固态浆料热流场分布规律，对发展半固态铸造成型工艺有着十分重要的意义．

在多孔介质数值模拟的方法中，Lattice⁃Boltzmann方法（LBM）无需直接求解复杂的偏微分方程，处理边界条件容易，计算并行好等特点在孔隙尺度［4］和REV尺度［5］都得到了很好的应用，并且也成功的用于处理多孔介质中的相变问题［6］．本文以变孔隙率多孔介质流动传热模型，利用LBM为计算工具，探索高固相率半固态浆料的传热和流动耦合作用关系．

1 数学模型

1．1 变孔隙率多孔介质流动传热宏观方程组

在半固态浆料传热过程中，随着浆料固相率的增加，孔隙率也相应的降低，而固相率是温度的函数，因此可以建立孔隙率与温度的函数关系．根据不可压缩多孔介质的广义N－S方程，在忽略粘性耗散热并且考虑相变时得到

u，T，p为液相的流速温度和压强等宏观物理量，ρ，υ为液相的密度、运动黏度和比热容等物性参数在本文中均视为常数．Q为相变潜热源项，ε，f为介质孔隙率和固相率，满足ε＋f＝1．σ，α为有效比热容和有效热扩散率，可以由固相和液相两者热物性参数加权得到

cp，kp，cs，ks分别为液相和固相的比热容和导热系数．

REV尺度研究多孔介质内的流动问题时，固相的骨架与占据骨架间空隙的流体相互作用是通过流体所受的阻力来表现的．多孔介质的外力项包括Darcy阻力，Brinkman黏性力和非线性的Forchheimaer阻力，根据文献［5］可以表示为

G为流体所受的体积力，νs为流体黏性系数，Fε为多孔介质的形状函数，K为渗透率，其表达式根据Ergun关联式［7］得到，具体表达形式如公式（5）所示．

其中dp为颗粒物直径．固相率与温度的关系受到合金物性和凝固条件的影响十分复杂加上热力学数据不完善，在平衡假设下常采用Scheil［8］方程求解：

TM为纯金属的熔点，TL为合金液相线温度．k溶质分配系数，当k＝0时方程变为线性方程．通过求导代换把潜热项合并入有效比热容项中：

1．2 格子Boltzmann速度分布函数演化方程

根据郭照立等人提出的GLBE［9］模型，广义N－S动量方程所对应的演化方程为

式中：fi（x，t）为粒子速度分布函数；（x，t）平衡态分布函数，ei为离散速度，i为离散速度方向．松弛时间平衡态分布方程和外力项的表达式为

式中：cs为格子声速，ωi为权系数．相应的宏观变量的表达式可以定义为

速度量的具体求解方法可参考文献［10］计算步骤获得．

1．3 格子Boltzmann温度分布函数演化方程

对于能量方程，可以改写成标准无源项对流扩散方程形式：

针对本文的宏观方程，式中

根据文献［11］可以得到变化

将对流扩散方程改写为

由文献［12］得到对应的演化方程为

式中C为平衡态的二阶矩：

类似动量方程可以写出松弛时间和温度宏观表达式：

在本文中采用固相率与温度为线性关系简化，则α′是与温度无关的常数，τg为常数．本文计算中的速度边界条件采用非平衡态反弹格式［13］，热边界条件则仿照非平衡态外推格式［14］写出．

2 数值模拟与验证

2．1 数值验证

文献11对上述GLBE模型模拟平板间均质多孔介质Poiseuille流动的数值验证，结合本文研究的具体对象，将变孔隙率条件下GLBE模型对速度场模拟与有限差分解做了对比验证．Poi⁃seuille流动的控制方程为

数值验证中设置孔隙率ε沿y方向以三种不同孔隙分布，设y∗＝y/Hy，三种条件下ε与y∗的函数关系式分别为

模拟条件为Re＝100，由图1可知，在三种情况下GLBE的模拟结果与有限差分解的模拟结果基本一致，通过（a）和（b）条件的模拟表明随着多孔介质孔隙率的线性增大速度非线性的增加，这体现了强制对流非线性阻力项的作用．

模拟区域中没有对流作用时，其凝固传热变为纯扩散问题．本文将一维非稳态纯扩散方程的有限差分方法（FDM）和热LBM对半固态凝固过程模拟做了对比分析．模拟设定冷边界TS＝0．88，热边界温度TH＝1．0，半固态合计的凝固温度区间为0．90～0．95，初始状态T0＝1．0，表示凝固潜热与显热之比的斯蒂芬数以及表示无量纲时间的傅里叶数为：Ste＝cpL/ΔT，Fo＝ατ/L2，ΔT为固液相的温差，L为特征长度．定义固相率为0．5的界面为合金的凝固前沿位置．图2可以看出LBM对半固态凝固过程模拟结果与有限差分解的吻合良好．Ste数越小潜热的作用越大，凝固的过程需要释放潜热大，凝固过程较慢，因此在同一Fo数时温度较高．相应的Ste数小的合金的凝固过程也会相对慢一些（如图3所示），图3中还可以得到在凝固过程初始时场内平均温度等于热边界温度，与冷边界温差较大，因此在定温条件下冷边界热流大，凝固过程较快．随着凝固过程进行，凝固前沿向热端推进，模拟区域内的平均温度降低，冷边界热流量不断减少，热边界的热流量增大，致使凝固速度逐渐减慢．

图1 不同孔隙率结构速度分布LBM模拟结果与差分解对比

图2 凝固过程温度分布LBM模拟与差分解对比

图3 凝固前沿位置LBM模拟与差分解对比

2．2 压差作用下平板间半固态浆料渗流过程模拟

对有些不伴随浆料较大塑性变形的铸造工艺，如液芯轻压下、半固态连铸、机械振动、充型铸造过程，可以近似认为固相骨架没有遭到破坏，半固态的液相在重力或压力作用下发生轻微的液固分离现象，此时可近似认为是半固态浆料在未凝固多孔固相骨架内的渗流过程．

本节通过LBM模拟研究了半固态浆料在外部压差作用下在两平板间流动传热过程，探索多孔介质模型在此类热流场模拟中的可行性．平板间距离为H、板长L、L/H＝8，两板有恒定温度TS＝0．92，并保持静止．半固态浆料中的固相假设与平板无相对位移，即速度为0．液相在两端有一定压差作用下流动，液固相间产生一定的速度差．定义义表征板长两端压差关系欧拉数为Eu＝Δp/（ρu02），u0流动初始时刻平均速度．液相入口段温度恒定并等于初始温度T0＝0．97，半固态合金的凝固区间为0．9～1．0，定义无量纲温度θ＝（TTS）/（T0－TS）．

图4为Ste＝0．5．Eu＝3×106时不同时刻x/L＝0．5时平板内H方向无量纲温度分布，可以看到随着时间推移，平板间的流体在流动过程中不断被冷却，温度逐步降低．随着相应的固相率增加，液相的合金不断的凝固构成多孔介质骨架的一部分，多孔介质的孔隙率降低．图5反映了不同Fo数时x/L＝0．5时平板内H方向Re数的变化，随着半固态区间的温度降低，两平板间的液相流速也不断降低，这是由于凝固，孔隙率的减少液相流动阻力增大的结果．并且Re数的降低速度比温度降低程度更为迅速，分析其原因为在对流条件下孔隙率的降低对流动阻力有非线性作用，即前文所提到的Forchheimaer力作用，在Fo＝0．5时，流场中心的Re数基本和边界区域Re数差距已经很小．

图4 不同Fo数平板间的温度分布

通过计算流场内不同时刻流速和温度的平均值就可以得到半固态浆料渗流过程平均Re数和无量纲温度随时间变化的趋势，图6为Eu＝3× 106时不同Ste数下无量纲温度θ随Fo数的变化规律，由图可见，初始温度相同半固态浆料由于Ste数不同形成不同的冷却曲线．同一时刻Ste数高的合金温度较低，这与数值验证部分所说明的潜热大的半固态合金温度下降较慢是一致的．图7为Eu＝3×106时Re数随Fo数变化趋势，图7表明随着浆料Ste的增大同计算区域内的同一时刻平均Re数也随之降低，这是由于合金的凝固潜热影响温度场，也造成固相率增加速度的差异，Ste数小较难凝固的半固态合金下相同时间时孔隙率也较高，从而在相同压差下流动阻力小，流场内的液相流动速度也相对高一些．从两图对比也可以看出因为非线性力的作用温度曲线的曲率较速度曲线的曲率要小一些．

图5 不同Fo时平板间的Re数分布

图6 不同Ste下无量纲温度随Fo数变化

图7 不同Ste下Re数随Fo数变化

3 半固态合金双辊连续铸轧成形过程LBM模拟

在薄带连续铸轧技术是一种近些年发展迅速的短流程板带成形技术，双辊薄带连续铸轧技术是众多该类成形技术中极具前途的一种［15］．半固态合金铸轧工艺如图8所示，制备好的半固态浆料通过铸嘴送入一对反向等角速度旋转的轧辊中，由于轧辊中有强制水冷，使得浆料在轧辊缝隙内迅速凝固形成薄带．整个过程在极短的时间内连续完成，温度场的分布和热流的传递直接决定了镁合金的凝固质量．本节以镁合金AZ91D半固态铸轧成形过程为对象，探索利用本文中提出的LBM方法和多孔介质模型来描述在半固态浆料在双棍间的热流耦合凝固过程．

图8 双辊薄带连续铸轧原理示意图

3．1 模拟参数的选取与无量纲化

根据文献［16］列出的实际工艺参数本节采用的模拟几何参数和速度边界为：轧辊直径：d＝210 mm，辊缝宽度d0＝4 mm，轧辊速度um＝0．6 m/min，进口温度825～855 K．

为了完成LBM模拟需要对实际的物理参数进行无量纲转，根据文献［16］，模拟涉及的AZ91D物性参数如表1所示．定义Re数为Re＝um·d0/ν，Pr数为Pr＝ρcpν/λ，重新定义无量纲温度θ＝（TTL）/（TH－TL），将表1中数据带入计算得到，Re＝20，Pr＝0．052，Ste＝350．87，Eu＝0，入口处无量纲温度θ0＝0．656～0．896．

表1 AZ91D热物理性能

3．2 液固相间相互作用力确定

与2．2节的分析不同，半固态浆料的连续铸轧过程中，合金已经凝固的部分和未完全凝固的固相骨架在轧辊的带动下向出口方向运动．为了描述这种固相运动过程中的多孔介质传热流动过程，假设多孔介质的固相部分满足无滑移条件，即固相部分与轧辊边具有相同的线速度um，定义液相与固相的相对速度为ur＝u－um．相应的公式（4）修改为

由公式可以看出ur＞0时，固相骨架对液相的流动表现出阻力作用，而ur＜0时固相对液相的流动具有推动作用，并设定轧辊壁面具有一定的过冷度θs＝－0．2．

3．3 模拟结果与分析

作为LBM在材料凝固成形过程中的初步探索，本文仅不同的入口温度下AZ91D合金连铸铸轧热流场进行模拟，以验证研究方法的可行性．图9为不同浇铸入口温度下达到稳态是的半固态浆料固相率分布，从图中可看到随着浇铸温度的提高，合金在出口处的固相率逐渐减低．无量纲温度θ0＝0．656时，浆料还未到辊缝出口处变完全凝固为固相，此时固相板带在轧辊的压力作用下容易开裂，造成不合格产品，θ0＝0．816时浆料还未完全凝固就在轧辊推动下从辊缝送出，导致铸轧失败．而在本文所采取的冷却边界条件下，无量纲θ0＝0．736时浆料在辊缝处恰好完全凝固，是最理想的工艺条件．图10为θ0＝0．736是半固态镁合金在铸轧区域内速度绝对值分布情况．

图9 不同入口温度时的半固态浆料温度场分布（a－θ0＝0．656，b－θ0＝0．736，c－θ0＝0．816）

图10 θ0＝0．736半固态浆料的速度场分布

可以看到在出入口压力相等的情况下，半固态浆料一边被轧辊凝固一边在轧辊的作用下向出口方向运动．在最贴近轧辊的薄层内已经完全凝固的浆料在无滑移条件下与轧辊保持相同速度．而入口处因为固相率低固相骨架对液相的推动作用较弱，流动区域最宽，所以速度最小．随着凝固的进行固相率不断升高，并且双辊宽度逐渐缩小，半固态浆料的速度不断升高，直至进入辊缝位置是速度和轧辊速度保持一致，在完全凝固成板带后送出铸轧区域．

4 结 论：

本文通过将孔隙率与温度耦合，得到模拟半固态相变现象的变孔隙率LBM模型．将相变源项并入非稳态项避免了求解有源项的对流扩散方程，得到求解半固态相变传热问题的热LBM方程．对平板间外部压差作用下高固相率浆料的渗流冷却过程和板带铸轧成型过程进行热流耦合模拟，具体结论如下：

1）所提出提出热流耦合LBM模型的相关模拟结果与有限差分解一致，验证了模型的正确性．

2）半固态浆料平板间渗流过程中，半固态浆料液相流动速度随冷却时间增加而降低．浆料的相变潜热越大，传热降温过程越慢，相同压差作用下流场Re数越高．Re数随温度的变化趋势表现出非线性．

3）利用多孔介质模型和LBM方法在通过无量纲处理后，对AZ91D镁合金半固态连续铸轧过程的温度场分布和流场分布进行了初步的模拟初步探索．得到了不同浇铸温度时铸轧区域内的固相率分布变化，和一定冷却条件下的最佳浇铸温度，分析了铸轧区域内的速度分布．

［1］KHOSRAVANIA，AASHURI H，DAVAMI P．Liquid segregation behavior of semi⁃solid AZ91 alloy during back extrusion test［J］．Journal of Alloys and Com⁃pounds，2009，477：822－827．

［2］SUN Z，LOG R E，BERNACKIM．3D finite element model of semi⁃solid permeability in an equiaxed granu⁃lar structure［J］．Computational Materials Science，2010，49：158－170．

［3］ZHU M F，LEE S Y，HONG C P．Modified cellular automatonmodel for the prediction of dendritic growth with melt convection［J］．Phys Rev E，2004，69（6）：061610．

［4］RAULM．Numerical simulations of surface reaction in porousmedia with lattice Boltzmann［J］．Chemical En⁃gineering Science，2012，69：628－643．

［5］RONG Fumei，GUO Zhaoli，CHAI Zhenhua，et al．A Lattice Boltzmann Model for axsymmetriy thermal flows through Porousmedia［J］．International Jounal of Heat and Mass Transfer，2010，53：5519－5527．

［6］杲东彦，陈振乾．格子Boltzmann方法模拟高Darcy数多孔介质内融化传热过程［J］．化工学报，2011，62（2）：321－328．GAO Dongyan，CHEN Zhenqian．Lattice Boltzmann method for heat transfer of melting in porous media with high Darcy number［J］．Journal of Chemical In⁃dustry and Engineering（China），2011，62（2）：321－328．

［7］ERGUN S．Flow through packed columns［J］．Chem．Eng．Prog．．1952，48：89－94．

［8］FLEMINGSM C．Behavior ofmetal alloys in the semi⁃solid state［J］．Metal．Tras，1991，22B：269－293．

［9］GUO Z L，ZHAO T S．Lattice Boltzmannmodel for in⁃compressible flows through porous media［J］．Phys Rev．E，2002，66，036304．

［10］GUO Z L，ZHAO T S．A latticce Boltzmann model for convection heat transfer in porousmedia［J］．Numer Heat Trans Par t B，2005，47：157－177．

［11］XIANG Xiugiao，WANG Zhenghua，Shi Baochang．Modified lattice Boltzmann scheme for nonlinear con⁃vection diffusion equations［J］．Communications in Nonlinear Science and Numerical Simulation，2012，17（6）：2415－2425．

［12］SHIB C，GUO Z L．Lattice Boltzmannmodel for non⁃linear convection diffusion equations［J］．P hys Rev E，2009，79：016701．

［13］ZOU Q S，HE X Y．On pressure and velocity boundary conditions for the lattice Boltzmann BGK model［J］．Physics of Fluids，1995，7（12）：2928－2930．

［14］TANG G H，TAO W Q，HE Y．LThermal boundary condition for the thermal lattice Boltzmann equation［J］．Phys Rev．E，2005/72（1）：016703－1－6．

［15］HENDRIEKSC．Strip Casting Technology⁃revolution in the Steel in dustly［C］．MPT intemational，1995，18（3）：42－46．

［16］谢水生，杨浩强，黄国杰，等．半固态镁合金连续铸轧过程的数值模拟［J］．塑性过程学报，2007，14（1）：80－84．XIE Shuisheng，YANG Haoqiang，HUANG Guojie，et al．Numerical simulation of semi⁃solid Magnesium alloy in continuous casting process［J］．Journal of Plasticity Engineering，2007，14（1）：80－84．

（编辑 张积宾）

LBM simulation ofmagnesium alloy rheological cast⁃rolling w ith porousmedium model

ZHANG Ying，MA Qiang，CHEN Yanghua，NIZelian

（School of Mechanical and Electronic Engineering，Nanchang University，Nanchang 330031，China）

A thermal⁃fluid coupling numericalmodel was built to analyze the process of high solid fraction semi⁃solid slurry seepage according to the lattice⁃Boltzmann method of REV scale porousmedium．Themodel was applied to simulate the Poiseuille flow with various porosity and the pure diffusion solidification．The simulation results of LBM is consistent with the simulation of finite difference，and it′s correctness has been verified．The seepage cooling and the thermal⁃fluid coupling process of slurry between the flat was analyzed under differential pressure．The simulation results show that the slurry flow velocity decreases as temperature decreases．The RE number of slurry seepage gets higher as the solidification latent heat is raising at the same time．The proposed model has been utilized to tentatively explore the heat⁃flow field distribution in the process of AZ91D magnesium alloy semi⁃solid casting and rolling．

lattice⁃boltzmann；porousmedia；semi⁃solid；thermal⁃fluid coupling

O359

A

1005－0299（2014）02－0117－06

2012－08－12．

国家自然科学基金（51066005）；江苏省研究生创新专项基金资助（YC2013－S028）．

张 莹（1970－），女，博士，副教授．

陈杨华，E⁃mail：chyhjx＠126．com．