基于元胞自動機-格子玻爾茲曼模型的枝晶碰撞行為模擬*

張士杰 王穎明 王琦 李晨宇 李日

(河北工業大學材料科學與工程學院,天津 300401)

合金凝固過程中,游離枝晶在熔體中的運動行為是研究合金凝固組織形成過程的關鍵問題之一.元胞自動機-格子玻爾茲曼耦合模型是近年來進行凝固微觀組織數值模擬的主要數值模型.本文改進了模擬枝晶生長的元胞自動機和格子玻爾茲曼模型,使其能夠模擬過冷熔體中等軸晶的運動行為.改進的模型采用伽利略不變的動量交換法計算流體力,通過求解質心運動方程計算枝晶的運動位移,并通過動網格技術實現枝晶的運動,運用硬球模型處理枝晶的碰撞.采用該模型模擬了Al-4.7%Cu 合金過冷熔體中單枝晶的沉降、牛頓流體中兩圓形粒子的沉降和兩枝晶的彈性碰撞.模擬結果表明,本模型在計算枝晶生長運動過程時可以很好地維持枝晶的形貌.采用本模型計算等軸枝晶的碰撞過程表明,枝晶的運動會擾動其周圍的熔體,造成周圍熔體濃度顯著變化,進而影響枝晶的生長,加劇枝晶生長的不對稱性.

1 引言

合金凝固過程中,鑄型壁面的表層細晶粒可能會脫落,已經凝固的枝晶臂(特別是二次枝晶臂)也可能發生熔斷現象,再加上在熔體中隨機出現的數量可觀的自由等軸晶,所以在液相區會出現大量游離的等軸晶,這些游離的等軸晶在重力和對流的推動下在合金熔體中運動(平移和旋轉),且常常會發生碰撞現象[1].熔體中游離枝晶的運動對鑄件的組織形貌和成分分布會有很大的影響[2].

當前,對合金凝固過程中枝晶生長的數值模擬研究已經比較成熟,可以有效地計算出純擴散和對流作用下枝晶的生長過程[3-5],并已發展到了三維(3 dimension,3D)[6-8].而關于枝晶運動的模擬,由于在計算枝晶運動的同時還要同時計算枝晶生長,所以有一定難度.目前該方面的研究比較少,已發表的枝晶運動的文獻中大多采用相場(phase field,PF)方法[9-13].采用PF 計算出的界面是連續擴散界面,對運動界面的處理有天然的優勢,但是計算量非常大,計算效率較低,計算規模小.而具有計算簡單快速、天然并行等優點的元胞自動機(cellular automata,CA)模型在模擬微觀組織方面的應用,大多數都只模擬了靜止枝晶的生長[14-17],對于枝晶運動的模擬研究很少,且不能很好地維持枝晶的形貌[18,19].本文旨在改進CA 模型,使其在模擬枝晶運動時能夠較好地維持枝晶的形貌和成分,在提升效率的同時提高模擬精度;同時,在模型中還加入了碰撞處理機制,使該模型能夠處理枝晶運動引起的碰撞現象.

2 模擬方法

采用格子玻爾茲曼(lattice Boltzmann,LB)模型計算流場、溫度場和濃度場,采用CA 方法模擬微觀枝晶生長,通過求解運動方程來計算枝晶的運動位移(平動和轉動),采用硬球模型計算枝晶碰撞(假設枝晶碰撞為完全彈性碰撞),耦合上述方法建立起枝晶運動生長的數值模型以研究枝晶在熔體中的運動行為.

2.1 LB 模型

采用單松弛時間的D2Q9 雙分布函數模型[20],流場、溫度場和濃度場的玻爾茲曼方程為

其中ωi是權重函數;c是格子速度;u是流體宏觀速度,溫度場和濃度場的格式和流場一樣.流體的流場、溫度場和濃度場的松弛時間τf,τt,τc分別和動力學黏度ν、熱擴散系數α、溶質擴散系數D有關:

流場、溫度場和濃度場的源項為

其中 Δfs是固相率增量;L是潛熱;Cp是比熱容;Cl是液相濃度;k是平衡分配系數;F是粒子受到的合力,由Boussinesq 近似[21]給出:

其中g為重力加速度;ρ0,T0,C0分別為初始密度、溫度和濃度;βT和βC分別為溫度膨脹系數和溶質膨脹系數.

流體宏觀密度ρ,速度u,溫度T和濃度C分別為

2.2 CA 模型

溶質擴散模型最早由Rappaz和Thévoz[22]于1987 年提出的,之后經過多次改進.2007 年,Zhu和Stefanescu[23]提出了用于定量化計算合金凝固過程中樹枝晶生長的虛擬前沿跟蹤模型.該模型認為合金枝晶生長是由局部液體平衡成分和局部液體實際成分的差異所驅動的.根據局部平衡的熱力學概念,界面平衡溫度T*定義為

fs為元胞的固相分數.根據Gibbs-Thomson公式,界面能各向異性函數可以表示為

其中δ為各向異性系數.因此,界面平衡濃度可以計算為

如果當前時刻界面元胞液相濃度小于界面平衡濃度,則該元胞固相增加.固相率增量為

2.3 枝晶運動方程

將運動枝晶看作剛體處理,枝晶在熔體中運動所受到的流體力采用伽利略不變的動量交換法[24]計算,流體對α枝晶的邊界點x的作用力為

其中Ω是α枝晶的邊界區域;ρl和ρs分別為流體和固體的密度,固體質量

枝晶的速度和角速度為

其中Iα為慣性張量.

2.4 碰撞模型

游離的枝晶在熔體中碰撞的情況很復雜,初步假定枝晶間的碰撞為完全彈性碰撞,因此采用硬球模型來處理碰撞步驟.將枝晶看做剛性粒子處理.對枝晶應用動量定理和動量矩定理:

其中m為枝晶質量;v0和v是碰撞前、后的速度;ω0和ω是碰撞前、后的角速度;S為內力沖量.

枝晶碰撞的瞬間,由于兩枝晶碰撞所受內力沖量相等,根據沖量定理和沖量矩定理,可以求得

其中ω10,ω20和ω1,ω2分別為碰撞前、后枝晶1和枝晶2 的角速度;v10,v20和v1,v2分別為碰撞前、后枝晶1和枝晶2 的速度;I1,I2為枝晶1,2 的轉動慣量;r1和r2分別為碰撞點與枝晶1和枝晶2 質心的位置矢量.

再引入恢復系數e,并與(26)式聯立寫成矩陣AX=b形式:

2.5 動網格方案

雖然CA 方法計算效率高,但是在處理移動邊界問題上存在困難.由于CA 方法的界面是尖銳界面,所以其界面在空間上的移動是不連續的,然而濃度在空間上的傳輸必須保證是連續的,這導致了在處理枝晶運動時溶質濃度在枝晶周圍積聚,人為地改變了枝晶的環境,影響了枝晶的生長.為了解決這個問題,本文改進了CA-LB 模型,采用局部動網格技術處理枝晶的運動.

每個枝晶都擁有一個局部的可移動網格,動網格會隨著枝晶的長大而增大,動網格的角速度和平動速度與內部枝晶一樣(如圖1所示).在背景網格中求解流場和溫度場,計算枝晶受到的力和轉矩.在每個動網格和背景網格中同時求解濃度場,將動網格中的溶質濃度覆蓋背景網格,動網格中流體節點的速度是相對速度,

圖1 動網格方案示意圖Fig.1.Schematic diagram of dynamic grid scheme.

其中,M(θ)是旋轉矩陣,v1是流體節點在動網格中的速度,vg是流體節點在背景網格中的速度,vd是動網格在該點的合速度,由方程(19)求出.旋轉矩陣被定義為

其中,θ是枝晶的轉動角度.動網格的邊界從全局網格插值得出.枝晶的生長和捕獲在動網格中計算,將枝晶的固相分數在每一步插值到背景網格中.

對于全域而言,可能會出現多個動網格部分重疊的情況.對于重疊的部分:如果都是流體節點,方程(3)中的源項為0,方程退化為簡單的對流擴散方程,因為動網格中的速度是相對速度,每個動網格中的流動相對于背景網格都是一樣的,所以每個動網格中的濃度值都是相等的;如果有界面節點和固體節點,則選定該值為背景網格的值.本方法的優點是將不連續的枝晶移動轉化為連續的熔體流動,通過局部網格中流體的相對流動體現枝晶的運動,避免了直接處理枝晶時溶質的積聚現象.

2.6 算法

模型的耦合算法流程如圖2所示,算法如下:

圖2 算法流程圖Fig.2.Algorithm flowchart.

1) 初始化計算域,給計算域中每個元胞的密度、溫度和濃度屬性賦初值,設置晶核在計算域中的位置和擇優生長方向,計算枝晶的質量、質心和轉動慣量;

2) 求解方程(1)計算流場,根據方程(18)—方程(21)計算枝晶受到的力和扭矩;

3) 計算全局溫度場和溶質場,方程(2)和方程(3);

4) 采用雙線性插值方法將背景網格中的濃度、溫度和速度插值到局部動態網格中;

5) 計算局部網格中的溫度場和濃度場,在局部網格中求解方程(13)—方程(17)計算枝晶的生長和捕獲,更新枝晶的質量、質心和轉動慣量;

6) 根據方程(22)和方程(23)計算枝晶的位移和轉動角,據此移動局部網格;

7) 將局部網格中的固相分數fs、溫度和濃度插值到背景網格;

8)枝晶碰撞處理;

9) 重復(2)—(8)步,直至程序結束.

3 模擬結果與討論

選取Al-4.7%Cu(質量含量)合金為模擬材料,其物性參數如表1所列.

表1 Al-4.7%Cu(質量含量)合金的熱物性參數[26]Table 1.Physical properties of Al-4.7%Cu (weight percent) alloy[26].

3.1 自然對流下單枝晶沉降

將計算區域劃分為300× 600 個網格,網格間距為0.5 μm,時間步長為 2.5×10-6s,初始過冷度為7 K,流場邊界為非平衡外推邊界條件,溫度場和濃度場采用無擴散邊界條件,固液邊界為無滑移邊界條件.初始時刻,在計算域(150,450)處設置1 個擇優生長角為45°的晶核.圖3為自然對流下單枝晶沉降過程中溶質濃度和流體速度演化.在枝晶運動早期,由于枝晶運動速度較小,運動及熔體對流對枝晶生長的影響較小,枝晶生長比較均衡,四個一次枝晶臂形貌比較接近,如圖3(a)所示.隨著枝晶的下落,枝晶迎流端生長加快,二次枝晶臂相比于順流端發達,枝晶的下游熔體的溶質濃度高于上游的溶質濃度,如圖3(b)和圖3(c)所示.

將運動對枝晶生長的影響進行定量分析.圖4為圖3(b)中O點(150,250)豎直方向上的溶質濃度變化.圖中橫坐標r為距O點的距離,O點水平線下邊為枝晶的迎流端上邊為順流端.從圖4可以看出,迎流端的峰值比順流端的高,這是因為枝晶呈45°生長時,隨著枝晶的下落溶質會積聚在枝晶凹形的枝晶臂間;在溶質變化平穩后,順流端的溶質濃度略高于迎流端,即如圖3(b)所示,在枝晶后會形成一段淡淡的拖尾.圖5為枝晶尖端生長速度隨時間的變化.從圖5中可以看出,在枝晶凝固初期,無論是上游尖端還是下游尖端的生長速率基本一致,都處于快速下降,這是由于凝固初期枝晶生長所需的過冷以熱過冷為主.隨著生長進入穩定階段,可以看出上游尖端的生長速度明顯高于下游尖端,上游尖端的穩態生長速度平均約為2.49 mm/s,下游尖端的穩態生長速度平均約為1.47 mm/s,上游尖端的生長速度比下游尖端的生長速度快了約1.7 倍.從圖3可以看出,相比于之前的研究[18,19],采用本文提出的動網格方案計算枝晶移動可以很好地保持枝晶的形貌.對于單個枝晶而言,小的計算域便可以展示枝晶的運動情況,因此設置的計算域較小,應用本模型在Intel Reon E5-2650 v42.20 GHz CPU,內存64 G 的服務器上單線程計算時間為1085 s,而采用openMP加速后,12 線程計算共用時246 s.

圖3 自然對流下單枝晶沉降過程中溶質濃度和流體速度的時間演化 (a) t=0.005 s;(b) t=0.015 s;(c) t=0.025 sFig.3.Evolution of solute concentration and fluid velocity during precipitation of single dendrite under natural convection:(a) t=0.005 s;(b) t=0.015 s;(c) t=0.025 s.

圖4 t=0.015 s 時,模擬域O 點處豎直方向上的溶質濃度變化Fig.4.When t=0.015 s,the solute concentration change in the vertical direction at point O in the simulation domain.

圖5 枝晶尖端生長速度隨時間的變化Fig.5.The change of dendrite tip growth rate with time.

3.2 牛頓流體中兩圓形粒子的沉降

為了驗證碰撞程序的可行性,模擬了牛頓流體中兩圓形粒子的沉降.模擬域設置為寬W=2 cm,高H=8 cm,兩密度ρp=1.01 g·cm-3,直徑D=0.2 cm 粒子位于模擬域水平方向的中間,豎直方向兩粒子分別在7.2和6.8 cm 處.域中充滿密度ρf=1 g·cm-3,運動黏度νf=0.01 cm2/s 的水.在LB 模型中,將模擬域剖分為200× 800 的網格,松弛時間τ=0.65,上下邊界采用非平衡外推格式,左右邊界采用無滑移邊界條件.這里使用的條件與Feng和Michaelides[27]使用的條件一樣.初始時刻,兩粒子和流體都處于靜止狀態,粒子受到重力和浮力的作用開始運動.

圖6(a)—(d)分別是t=1.0 s,t=1.5 s,t=2.5 s,t=4.0 s 時兩粒子沉降過程的瞬時渦量輪廓.圖7(a)和圖7(b)分別為粒子質心橫、縱坐標隨時間變化的圖像.多粒子沉降時會出現漂移-接觸-翻滾(drafting-kissing-tumblin,DKT)現象,由圖6中可以看出,運動早期粒子1 被粒子2 沉降產生的尾渦捕獲使其沉降速度加快,到1.5 s 時與粒子1 接觸,直至2.5 s 時兩粒子發生翻滾.運用本文的模型成功模擬再現了Feng 的結果.在圖7中比較了Feng和Michaelides[27]的結果,可以看出粒子下落的軌跡是一致的,存在的一些誤差是本文對流體力和碰撞的處理與Feng和Michaelides 的不同,本文計算流體力采用的動量交換法滿足伽利略不變性,排除了邊界速度的影響,Feng和Michaelides[27]采用的動量交換法不滿足伽利略不變性.由于粒子邊界速度的影響造成計算流體力時存在誤差.

圖6 兩個球形粒子沉降過程中,四個時間的瞬時渦量輪廓 (a) 1.0 s;(b) 1.5 s;(c) 2.5 s;(d) 4.0 sFig.6.The instantaneous vorticity profiles of two spherical particles in the process of settling at four times:(a) 1.0 s;(b) 1.5 s;(c) 2.5 s;(d) 4.0 s.

圖7 粒子質心坐標 (a) 橫坐標;(b) 縱坐標Fig.7.Coordinates of the particles centroid:(a) Transverse coordinates;(b) longitudinal coordinates.

3.3 雙枝晶彈性碰撞

計算區域為600× 400 個網格,初始過冷度為7 K,溫度場邊界為絕熱邊界條件,溶質場為無擴散邊界條件.為了便于研究枝晶的碰撞行為,將初始流場設置為靜止狀態,在計算域(200,200)和(400,230)的位置放置兩個晶核,相對于水平方向的擇優生長角分別為0°和45°,組成為Cs=kCl.當枝晶生長一定大小時,枝晶不在生長并直接賦給兩個枝晶0.1和—0.1 (lattice unit)的水平速度,忽略重力的影響.

兩枝晶運動前的位置和形貌如圖8(a),左邊是枝晶1 右邊是枝晶2,兩個枝晶受到熔體的阻力作減速運動直至發生碰撞如圖8(b),碰撞前枝晶1,2 的水平速度分別為0.08和—0.082,碰撞后速度分別為—0.029和0.024,碰后角速度為0和0.0036 rad.圖8(c)和圖8(d)顯示兩枝晶碰撞后的運動情況.碰撞時,由于碰撞點和枝晶1 的質心在同一高度,所以碰撞產生的沖量對枝晶1 的轉矩為0,枝晶1不發生轉動.

圖8 枝晶偏心碰撞形貌圖 (a) 0.008 s;(b) 0.0089 s;(c) 0.009 s;(d) 0.00915 s.箭頭表示速度v 的大小和方向,灰度表示流體溶質濃度CFig.8.The morphology of dendrite eccentric collision:(a) 0.008 s;(b) 0.0089 s;(c) 0.009 s;(d) 0.00915 s.The arrow indicates the size and direction of velocity v,and the gray scale indicates the concentration of solute C.

圖9為0.008和0.009 s 時,枝晶2 周圍熔體的平均溶質濃度Cave,x-軸表示距離枝晶2 質心的元胞個數.Cave的計算規則如圖10所示,圖10中白色格子表示枝晶的質心,x表示距離質心的元胞數.當x=1 時,Cave等于淺灰色格子中所有液相格點的平均濃度;當x=2 時,Cave等于淺黑色格子中所有液相格點的平均濃度;其他與x值對應的Cave同理計算.從圖9中可以看出,在0.008 s時枝晶處于靜止生長狀態,枝晶2 周圍的溶質濃度在界面處最大隨后緩慢的減少直至穩定.這是因為枝晶生長時在固液界面處向外排出溶質,排出的溶質向周圍均勻的擴散,所以形成了這樣的濃度梯度.在0.009 s 時,溶質濃度的變化變得不均衡,在枝晶尖端附近濃度最高.因為枝晶2 在碰撞后發生了轉動,對周圍熔體進行了擾動,使周圍的熔體濃度變化混亂,將溶質甩在了枝晶的尖端附近.通過對雙枝晶彈性碰撞的計算可以證明本文所用碰撞模型可以表達不規則等軸枝晶的碰撞過程.

圖9 枝晶2 附近熔體平均溶質濃度Fig.9.Average solute concentration of melt near dendrite 2.

圖10 熔體平均溶質濃度計算規則Fig.10.Calculation rules of melt average solute concentration.

4 結論

本文建立了一個改進的CA-LB 模型,該模型采用動網格技術實現枝晶的運動,每個枝晶擁有一個可移動的局部網格,在局部網格內通過流體的相對流動來體現枝晶的移動;模型中還加入了處理枝晶的碰撞的機制.應用該模型模擬了自然對流下單枝晶的沉降過程,模擬結果表明,枝晶運動的上游端濃度變化梯度大,生長速度加快,二次枝晶臂發達;同時表明該模型可以很好地模擬枝晶的運動過程,并維持枝晶的形貌.對模型的碰撞機制進行了測試,模擬了牛頓流體中兩圓形粒子的沉降,成功再現了DKT 現象,證明了碰撞算法的可行性.計算了兩枝晶的彈性碰撞,模擬結果表明該模型可以處理不規則枝晶的偏心碰撞.