張 瑩，馬 強，陳楊華，倪澤聯(lián)

（南昌大學機電工程學院，南昌330031）

多孔介質(zhì)模型鎂合金流變鑄軋LBM模擬

張 瑩，馬 強，陳楊華，倪澤聯(lián)

（南昌大學機電工程學院，南昌330031）

根據(jù)REV尺度多孔介質(zhì)格子－Boltzmann方法基本理論，建立求解高固相率半固態(tài)漿料滲流過程的熱流耦合數(shù)值模型．應用模型對變孔隙率Poiseuille流動及純擴散凝固進行模擬．LBM模擬結(jié)果與有限差分模擬吻合．分析了壓差作用下漿料在平板間滲流冷卻與熱流耦合過程．模擬結(jié)果表明漿料滲流速度隨溫度降低而減慢，相同時刻下凝固潛熱較大的漿料滲流Re數(shù)高．利用提出的模型初步探索了AZ91D鎂合金半固態(tài)流變鑄軋過程的熱流場分布．

格子－Boltzmann；多孔介質(zhì)；半固態(tài)；熱流耦合

作為半固態(tài)金屬凝固過程的最后階段，高固相率漿料內(nèi)的傳輸行為很大程度上影響了材料凝固的最終質(zhì)量．目前，有部分學者通過多孔介質(zhì)模型去研究了高固相率半固態(tài)漿料的擠壓成形過程中的應力和兩相分離行為并取得成功［1－2］，但是基于多孔介質(zhì)模型對半固態(tài)對流凝固過程的熱流場分析尚沒有相關(guān)的文獻報道．從孔隙尺度上說，研究液態(tài)合金在固相晶粒組織間的流動傳熱傳質(zhì)過程，分析晶粒生長環(huán)境，是獲得晶粒組織形態(tài)，控制凝固質(zhì)量的基礎(chǔ)［3］．從表征體元（REV）尺度上說，研究漿料的固相率和滲透率、物性參數(shù)對熱流場的影響，揭示固相率的半固態(tài)漿料熱流場分布規(guī)律，對發(fā)展半固態(tài)鑄造成型工藝有著十分重要的意義．

在多孔介質(zhì)數(shù)值模擬的方法中，Lattice?Boltzmann方法（LBM）無需直接求解復雜的偏微分方程，處理邊界條件容易，計算并行好等特點在孔隙尺度［4］和REV尺度［5］都得到了很好的應用，并且也成功的用于處理多孔介質(zhì)中的相變問題［6］．本文以變孔隙率多孔介質(zhì)流動傳熱模型，利用LBM為計算工具，探索高固相率半固態(tài)漿料的傳熱和流動耦合作用關(guān)系．

1 數(shù)學模型

1．1 變孔隙率多孔介質(zhì)流動傳熱宏觀方程組

在半固態(tài)漿料傳熱過程中，隨著漿料固相率的增加，孔隙率也相應的降低，而固相率是溫度的函數(shù)，因此可以建立孔隙率與溫度的函數(shù)關(guān)系．根據(jù)不可壓縮多孔介質(zhì)的廣義N－S方程，在忽略粘性耗散熱并且考慮相變時得到

u，T，p為液相的流速溫度和壓強等宏觀物理量，ρ，υ為液相的密度、運動黏度和比熱容等物性參數(shù)在本文中均視為常數(shù)．Q為相變潛熱源項，ε，f為介質(zhì)孔隙率和固相率，滿足ε＋f＝1．σ，α為有效比熱容和有效熱擴散率，可以由固相和液相兩者熱物性參數(shù)加權(quán)得到

cp，kp，cs，ks分別為液相和固相的比熱容和導熱系數(shù)．

REV尺度研究多孔介質(zhì)內(nèi)的流動問題時，固相的骨架與占據(jù)骨架間空隙的流體相互作用是通過流體所受的阻力來表現(xiàn)的．多孔介質(zhì)的外力項包括Darcy阻力，Brinkman黏性力和非線性的Forchheimaer阻力，根據(jù)文獻［5］可以表示為

G為流體所受的體積力，νs為流體黏性系數(shù)，F(xiàn)ε為多孔介質(zhì)的形狀函數(shù)，K為滲透率，其表達式根據(jù)Ergun關(guān)聯(lián)式［7］得到，具體表達形式如公式（5）所示．

其中dp為顆粒物直徑．固相率與溫度的關(guān)系受到合金物性和凝固條件的影響十分復雜加上熱力學數(shù)據(jù)不完善，在平衡假設下常采用Scheil［8］方程求解：

TM為純金屬的熔點，TL為合金液相線溫度．k溶質(zhì)分配系數(shù)，當k＝0時方程變?yōu)榫€性方程．通過求導代換把潛熱項合并入有效比熱容項中：

1．2 格子Boltzmann速度分布函數(shù)演化方程

根據(jù)郭照立等人提出的GLBE［9］模型，廣義N－S動量方程所對應的演化方程為

式中：fi（x，t）為粒子速度分布函數(shù)；（x，t）平衡態(tài)分布函數(shù)，ei為離散速度，i為離散速度方向．松弛時間平衡態(tài)分布方程和外力項的表達式為

式中：cs為格子聲速，ωi為權(quán)系數(shù)．相應的宏觀變量的表達式可以定義為

速度量的具體求解方法可參考文獻［10］計算步驟獲得．

1．3 格子Boltzmann溫度分布函數(shù)演化方程

對于能量方程，可以改寫成標準無源項對流擴散方程形式：

針對本文的宏觀方程，式中

根據(jù)文獻［11］可以得到變化

將對流擴散方程改寫為

由文獻［12］得到對應的演化方程為

式中C為平衡態(tài)的二階矩：

類似動量方程可以寫出松弛時間和溫度宏觀表達式：

在本文中采用固相率與溫度為線性關(guān)系簡化，則α′是與溫度無關(guān)的常數(shù)，τg為常數(shù)．本文計算中的速度邊界條件采用非平衡態(tài)反彈格式［13］，熱邊界條件則仿照非平衡態(tài)外推格式［14］寫出．

2 數(shù)值模擬與驗證

2．1 數(shù)值驗證

文獻11對上述GLBE模型模擬平板間均質(zhì)多孔介質(zhì)Poiseuille流動的數(shù)值驗證，結(jié)合本文研究的具體對象，將變孔隙率條件下GLBE模型對速度場模擬與有限差分解做了對比驗證．Poi?seuille流動的控制方程為

數(shù)值驗證中設置孔隙率ε沿y方向以三種不同孔隙分布，設y?＝y(tǒng)/Hy，三種條件下ε與y?的函數(shù)關(guān)系式分別為

模擬條件為Re＝100，由圖1可知，在三種情況下GLBE的模擬結(jié)果與有限差分解的模擬結(jié)果基本一致，通過（a）和（b）條件的模擬表明隨著多孔介質(zhì)孔隙率的線性增大速度非線性的增加，這體現(xiàn)了強制對流非線性阻力項的作用．

模擬區(qū)域中沒有對流作用時，其凝固傳熱變?yōu)榧償U散問題．本文將一維非穩(wěn)態(tài)純擴散方程的有限差分方法（FDM）和熱LBM對半固態(tài)凝固過程模擬做了對比分析．模擬設定冷邊界TS＝0．88，熱邊界溫度TH＝1．0，半固態(tài)合計的凝固溫度區(qū)間為0．90～0．95，初始狀態(tài)T0＝1．0，表示凝固潛熱與顯熱之比的斯蒂芬數(shù)以及表示無量綱時間的傅里葉數(shù)為：Ste＝cpL/ΔT，F(xiàn)o＝ατ/L2，ΔT為固液相的溫差，L為特征長度．定義固相率為0．5的界面為合金的凝固前沿位置．圖2可以看出LBM對半固態(tài)凝固過程模擬結(jié)果與有限差分解的吻合良好．Ste數(shù)越小潛熱的作用越大，凝固的過程需要釋放潛熱大，凝固過程較慢，因此在同一Fo數(shù)時溫度較高．相應的Ste數(shù)小的合金的凝固過程也會相對慢一些（如圖3所示），圖3中還可以得到在凝固過程初始時場內(nèi)平均溫度等于熱邊界溫度，與冷邊界溫差較大，因此在定溫條件下冷邊界熱流大，凝固過程較快．隨著凝固過程進行，凝固前沿向熱端推進，模擬區(qū)域內(nèi)的平均溫度降低，冷邊界熱流量不斷減少，熱邊界的熱流量增大，致使凝固速度逐漸減慢．

圖1 不同孔隙率結(jié)構(gòu)速度分布LBM模擬結(jié)果與差分解對比

圖2 凝固過程溫度分布LBM模擬與差分解對比

圖3 凝固前沿位置LBM模擬與差分解對比

2．2 壓差作用下平板間半固態(tài)漿料滲流過程模擬

對有些不伴隨漿料較大塑性變形的鑄造工藝，如液芯輕壓下、半固態(tài)連鑄、機械振動、充型鑄造過程，可以近似認為固相骨架沒有遭到破壞，半固態(tài)的液相在重力或壓力作用下發(fā)生輕微的液固分離現(xiàn)象，此時可近似認為是半固態(tài)漿料在未凝固多孔固相骨架內(nèi)的滲流過程．

本節(jié)通過LBM模擬研究了半固態(tài)漿料在外部壓差作用下在兩平板間流動傳熱過程，探索多孔介質(zhì)模型在此類熱流場模擬中的可行性．平板間距離為H、板長L、L/H＝8，兩板有恒定溫度TS＝0．92，并保持靜止．半固態(tài)漿料中的固相假設與平板無相對位移，即速度為0．液相在兩端有一定壓差作用下流動，液固相間產(chǎn)生一定的速度差．定義義表征板長兩端壓差關(guān)系歐拉數(shù)為Eu＝Δp/（ρu02），u0流動初始時刻平均速度．液相入口段溫度恒定并等于初始溫度T0＝0．97，半固態(tài)合金的凝固區(qū)間為0．9～1．0，定義無量綱溫度θ＝（TTS）/（T0－TS）．

圖4為Ste＝0．5．Eu＝3×106時不同時刻x/L＝0．5時平板內(nèi)H方向無量綱溫度分布，可以看到隨著時間推移，平板間的流體在流動過程中不斷被冷卻，溫度逐步降低．隨著相應的固相率增加，液相的合金不斷的凝固構(gòu)成多孔介質(zhì)骨架的一部分，多孔介質(zhì)的孔隙率降低．圖5反映了不同F(xiàn)o數(shù)時x/L＝0．5時平板內(nèi)H方向Re數(shù)的變化，隨著半固態(tài)區(qū)間的溫度降低，兩平板間的液相流速也不斷降低，這是由于凝固，孔隙率的減少液相流動阻力增大的結(jié)果．并且Re數(shù)的降低速度比溫度降低程度更為迅速，分析其原因為在對流條件下孔隙率的降低對流動阻力有非線性作用，即前文所提到的Forchheimaer力作用，在Fo＝0．5時，流場中心的Re數(shù)基本和邊界區(qū)域Re數(shù)差距已經(jīng)很小．

圖4 不同F(xiàn)o數(shù)平板間的溫度分布

通過計算流場內(nèi)不同時刻流速和溫度的平均值就可以得到半固態(tài)漿料滲流過程平均Re數(shù)和無量綱溫度隨時間變化的趨勢，圖6為Eu＝3× 106時不同Ste數(shù)下無量綱溫度θ隨Fo數(shù)的變化規(guī)律，由圖可見，初始溫度相同半固態(tài)漿料由于Ste數(shù)不同形成不同的冷卻曲線．同一時刻Ste數(shù)高的合金溫度較低，這與數(shù)值驗證部分所說明的潛熱大的半固態(tài)合金溫度下降較慢是一致的．圖7為Eu＝3×106時Re數(shù)隨Fo數(shù)變化趨勢，圖7表明隨著漿料Ste的增大同計算區(qū)域內(nèi)的同一時刻平均Re數(shù)也隨之降低，這是由于合金的凝固潛熱影響溫度場，也造成固相率增加速度的差異，Ste數(shù)小較難凝固的半固態(tài)合金下相同時間時孔隙率也較高，從而在相同壓差下流動阻力小，流場內(nèi)的液相流動速度也相對高一些．從兩圖對比也可以看出因為非線性力的作用溫度曲線的曲率較速度曲線的曲率要小一些．

圖5 不同F(xiàn)o時平板間的Re數(shù)分布

圖6 不同Ste下無量綱溫度隨Fo數(shù)變化

圖7 不同Ste下Re數(shù)隨Fo數(shù)變化

3 半固態(tài)合金雙輥連續(xù)鑄軋成形過程LBM模擬

在薄帶連續(xù)鑄軋技術(shù)是一種近些年發(fā)展迅速的短流程板帶成形技術(shù)，雙輥薄帶連續(xù)鑄軋技術(shù)是眾多該類成形技術(shù)中極具前途的一種［15］．半固態(tài)合金鑄軋工藝如圖8所示，制備好的半固態(tài)漿料通過鑄嘴送入一對反向等角速度旋轉(zhuǎn)的軋輥中，由于軋輥中有強制水冷，使得漿料在軋輥縫隙內(nèi)迅速凝固形成薄帶．整個過程在極短的時間內(nèi)連續(xù)完成，溫度場的分布和熱流的傳遞直接決定了鎂合金的凝固質(zhì)量．本節(jié)以鎂合金AZ91D半固態(tài)鑄軋成形過程為對象，探索利用本文中提出的LBM方法和多孔介質(zhì)模型來描述在半固態(tài)漿料在雙棍間的熱流耦合凝固過程．

圖8 雙輥薄帶連續(xù)鑄軋原理示意圖

3．1 模擬參數(shù)的選取與無量綱化

根據(jù)文獻［16］列出的實際工藝參數(shù)本節(jié)采用的模擬幾何參數(shù)和速度邊界為：軋輥直徑：d＝210 mm，輥縫寬度d0＝4 mm，軋輥速度um＝0．6 m/min，進口溫度825～855 K．

為了完成LBM模擬需要對實際的物理參數(shù)進行無量綱轉(zhuǎn)，根據(jù)文獻［16］，模擬涉及的AZ91D物性參數(shù)如表1所示．定義Re數(shù)為Re＝um·d0/ν，Pr數(shù)為Pr＝ρcpν/λ，重新定義無量綱溫度θ＝（TTL）/（TH－TL），將表1中數(shù)據(jù)帶入計算得到，Re＝20，Pr＝0．052，Ste＝350．87，Eu＝0，入口處無量綱溫度θ0＝0．656～0．896．

表1 AZ91D熱物理性能

3．2 液固相間相互作用力確定

與2．2節(jié)的分析不同，半固態(tài)漿料的連續(xù)鑄軋過程中，合金已經(jīng)凝固的部分和未完全凝固的固相骨架在軋輥的帶動下向出口方向運動．為了描述這種固相運動過程中的多孔介質(zhì)傳熱流動過程，假設多孔介質(zhì)的固相部分滿足無滑移條件，即固相部分與軋輥邊具有相同的線速度um，定義液相與固相的相對速度為ur＝u－um．相應的公式（4）修改為

由公式可以看出ur＞0時，固相骨架對液相的流動表現(xiàn)出阻力作用，而ur＜0時固相對液相的流動具有推動作用，并設定軋輥壁面具有一定的過冷度θs＝－0．2．

3．3 模擬結(jié)果與分析

作為LBM在材料凝固成形過程中的初步探索，本文僅不同的入口溫度下AZ91D合金連鑄鑄軋熱流場進行模擬，以驗證研究方法的可行性．圖9為不同澆鑄入口溫度下達到穩(wěn)態(tài)是的半固態(tài)漿料固相率分布，從圖中可看到隨著澆鑄溫度的提高，合金在出口處的固相率逐漸減低．無量綱溫度θ0＝0．656時，漿料還未到輥縫出口處變完全凝固為固相，此時固相板帶在軋輥的壓力作用下容易開裂，造成不合格產(chǎn)品，θ0＝0．816時漿料還未完全凝固就在軋輥推動下從輥縫送出，導致鑄軋失?。诒疚乃扇〉睦鋮s邊界條件下，無量綱θ0＝0．736時漿料在輥縫處恰好完全凝固，是最理想的工藝條件．圖10為θ0＝0．736是半固態(tài)鎂合金在鑄軋區(qū)域內(nèi)速度絕對值分布情況．

圖9 不同入口溫度時的半固態(tài)漿料溫度場分布（a－θ0＝0．656，b－θ0＝0．736，c－θ0＝0．816）

圖10 θ0＝0．736半固態(tài)漿料的速度場分布

可以看到在出入口壓力相等的情況下，半固態(tài)漿料一邊被軋輥凝固一邊在軋輥的作用下向出口方向運動．在最貼近軋輥的薄層內(nèi)已經(jīng)完全凝固的漿料在無滑移條件下與軋輥保持相同速度．而入口處因為固相率低固相骨架對液相的推動作用較弱，流動區(qū)域最寬，所以速度最?。S著凝固的進行固相率不斷升高，并且雙輥寬度逐漸縮小，半固態(tài)漿料的速度不斷升高，直至進入輥縫位置是速度和軋輥速度保持一致，在完全凝固成板帶后送出鑄軋區(qū)域．

4 結(jié) 論：

本文通過將孔隙率與溫度耦合，得到模擬半固態(tài)相變現(xiàn)象的變孔隙率LBM模型．將相變源項并入非穩(wěn)態(tài)項避免了求解有源項的對流擴散方程，得到求解半固態(tài)相變傳熱問題的熱LBM方程．對平板間外部壓差作用下高固相率漿料的滲流冷卻過程和板帶鑄軋成型過程進行熱流耦合模擬，具體結(jié)論如下：

1）所提出提出熱流耦合LBM模型的相關(guān)模擬結(jié)果與有限差分解一致，驗證了模型的正確性．

2）半固態(tài)漿料平板間滲流過程中，半固態(tài)漿料液相流動速度隨冷卻時間增加而降低．漿料的相變潛熱越大，傳熱降溫過程越慢，相同壓差作用下流場Re數(shù)越高．Re數(shù)隨溫度的變化趨勢表現(xiàn)出非線性．

3）利用多孔介質(zhì)模型和LBM方法在通過無量綱處理后，對AZ91D鎂合金半固態(tài)連續(xù)鑄軋過程的溫度場分布和流場分布進行了初步的模擬初步探索．得到了不同澆鑄溫度時鑄軋區(qū)域內(nèi)的固相率分布變化，和一定冷卻條件下的最佳澆鑄溫度，分析了鑄軋區(qū)域內(nèi)的速度分布．

［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數(shù)多孔介質(zhì)內(nèi)融化傳熱過程［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］謝水生，楊浩強，黃國杰，等．半固態(tài)鎂合金連續(xù)鑄軋過程的數(shù)值模擬［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）；江蘇省研究生創(chuàng)新專項基金資助（YC2013－S028）．

張 瑩（1970－），女，博士，副教授．

陳楊華，E?mail：chyhjx＠126．com．