朱昶勝, 曹雅星, 雷 鵬, 馮 力, 馬芳蘭

(1. 蘭州理工大學 計算機與通信學院, 甘肅 蘭州 730050; 2. 蘭州理工大學 省部共建有色金屬先進加工與再利用國家重點實驗室, 甘肅 蘭州 730050; 3. 蘭州理工大學 網(wǎng)絡(luò)與信息中心, 甘肅 蘭州 730050)

確定性方法通過確定過冷度函數(shù)來預(yù)測枝晶大小及生長速度[1].由于其本身的特性無法對枝晶生長過程中的隨機部分進行模擬，應(yīng)用受到較大限制.隨后Brown等[2-3]建立的元胞自動機模型，考慮了枝晶局部界面曲率的影響，增大了模擬尺度，定性地觀察了枝晶生長過程，但此方法無法模擬復(fù)雜枝晶的準確形貌，且對于體積過大、計算量大的鑄件無法達到計算要求.國內(nèi)外學者利用相場法(PF)對微觀組織進行模擬研究，定量模擬了宏觀作用下枝晶生長，較好地捕捉了微觀組織動力學演變過程的同時，研究各個參數(shù)對枝晶生長的影響，簡化了計算步驟，提高了計算效率.Wheeler等[4]在建立的WBM模型的基礎(chǔ)上，引入相場與溶質(zhì)場的修正項，建立新的模型，模擬了金屬凝固現(xiàn)象.朱昶勝等[5]研究了Al-Cu合金的傾斜枝晶生長動力學和定向凝固過程中的形態(tài)轉(zhuǎn)變，特別分析了冷卻速率和主晶間距對凝固組織的影響，研究表明，冷卻速率和主晶間距會在一定程度上影響傾斜枝晶的生長角度.

在液態(tài)金屬凝固過程中，對流對枝晶形貌及周圍溫度場、溶質(zhì)場等都有明顯影響，對于傳統(tǒng)流體力學計算方法，現(xiàn)已存在大量研究[6-7]，其中，基于Navier-Stoke方程[8]的流場數(shù)值計算方法是描述粘性不可壓縮流體動量守恒的運動方程，可進行離散求解，但僅適用于單向流問題，且不能滿足連續(xù)介質(zhì)的假設(shè)條件，這些缺點將會造成流場計算不易收斂、計算困難等問題.目前大多數(shù)研究是基于外力場作用導(dǎo)致的強迫對流，忽略了由于溫度和溶質(zhì)分布變化引起的自然對流，從而導(dǎo)致計算結(jié)果出現(xiàn)偏差，動量守恒方程存在計算量大和計算效率低等問題.格子玻爾茲曼方法(LBM)集合大量虛擬粒子，在離散格子中進行碰撞、遷移，形成流動現(xiàn)象，該方法自身具有并行計算的特點，程序易于實施，具備流體結(jié)構(gòu)互相作用描述簡單的優(yōu)勢，可以處理復(fù)雜的邊界和外力項問題.Miller等[9]提出了一種在對流情況下保持各向同性不變的動力學格子Bolzmann方法，耦合相場法，研究了不同瑞利數(shù)和普朗特數(shù)下純金屬Ga的生長過程以及固、液相變與對流之間的相互作用.

本文將計算枝晶生長的相場法與計算流場、濃度場和溫度場的格子玻爾茲曼方法進行耦合，建立了二維PF-LBM模型，模擬了 Al-4.5% Cu枝晶在自然對流作用下的凝固過程，研究了枝晶形貌、流場、溫度場和溶質(zhì)場的變化以及各向異性強度對枝晶生長的影響.

1 數(shù)學模型

1.1 枝晶生長相場模型

相場模型是一種建立在熱力學基礎(chǔ)上，描述系統(tǒng)演化動力學的模型，其核心思想就是引入一個或多個變量，用彌散界面代替尖銳界面來描述界面[10].在固、液兩相系統(tǒng)中引入一個序參量，也就是相場變量φ(r,t)來標識物質(zhì)狀態(tài)，液相取值為0或-1，固相取值為1.Ginzburg-Landau方程的關(guān)鍵在于具體體系給出的適當自由能泛函，本文采用的自由能泛函表達式如下：

(1)

根據(jù)熱力學第二定律，隨著凝固進行，體系自由勢能逐漸減小，將高溫液態(tài)金屬看作受力后極易變形，且切應(yīng)力與變形速率成正比的牛頓流體，其數(shù)值計算要考慮濃度、溫度分布影響，需要加入連續(xù)方程和動量方程[11]，利用先行不可逆動力學方程推導(dǎo)出自由能密度函數(shù)為

(2)

式(1、2)中：φ為相場變量；F為體系自由能；f為雙穩(wěn)態(tài)事函數(shù)；H為體系量綱的焓；c為濃度；ε為相場項系數(shù)；δ為溶質(zhì)場梯度系數(shù);D(φ)是溶質(zhì)擴散率，fc為自由能一階偏導(dǎo)數(shù)，fcc為自由能二階偏導(dǎo)數(shù).

此外，各向異性是需要考慮的重要因素之一，表達式如下:

ε(θ)=ε(1+γcoskθ)

(3)

加入各向異性得到的相場方程為

(4)

其中:方程(3)中ε為方向矢量;θ為法向與生長主軸夾角;k為各向異性模數(shù);γ是各向異性強度，數(shù)值越大強度越大，方程(4)中ε′、ε″表示對θ的一階、二階導(dǎo)數(shù)，φxy、φyy表示對x、y的二階偏導(dǎo)數(shù).

1.2 LBM流場模型

LBM是一種基于介觀模擬尺度的計算流體力學的方法，直接從離散模型出發(fā)，應(yīng)用質(zhì)量、動量和能量守恒的規(guī)律，從不同角度建立起宏觀與微觀、連續(xù)與離散之間的聯(lián)系.格子方程模型中DnQm[12-13]系列最具代表性，其中n為空間維數(shù)，m為離散速度數(shù).采用基于單松弛(BGK)模型的D2Q9模型，流體碰撞過程中，對應(yīng)的演化方程表達式為

fi(x+ciδt,t+δt)-fi(x,t)=

(5)

二維D2Q9模型和三維D3Q19常用來計算流動問題，如圖1所示，D2Q9模型是根據(jù)不同網(wǎng)格確定的平衡分布函數(shù)，其中一個虛擬粒子有9個離散速度，由于邊界條件不參與演化，因此要根據(jù)已知邊界條件來確定粒子分布函數(shù)之后才能進行下一步計算，對應(yīng)的離散速度表達式為

(6)

式中：ωi為權(quán)函數(shù)，i取0～8，cs為格子聲速.

圖1 晶格結(jié)構(gòu)Fig.1 Lattice structure

LBM用于計算受對流和擴散影響的濃度場和溫度場，其分布函數(shù)演化方程與流體分布函數(shù)類似，如下式所示[14-15]：

gi(x+eiΔt,t+Δt)-gi(x,t)=

(7)

hi(x+eiΔt,t+Δt)-hi(x,t)=

(8)

式中：gi(x,t)為濃度場分布函數(shù)，hi(x,t)為溫度場分布函數(shù)，Gi、Hi分別為枝晶生長過程中濃度梯度、溫度梯度造成的溶質(zhì)源和溫度源項.

1.3 數(shù)值模擬計算

由于銅鋁合金屬于面心立方晶系，具有高耐熱性和高韌性，因此選取Al-4.5%Cu二元合金，表1為合金物性參數(shù).初始狀態(tài)下，取溶質(zhì)濃度為c0，溫度為T0，模擬區(qū)域為800×800的正方形網(wǎng)格，具有初始濃度和初始溫度，模擬中心為400×400，四壁設(shè)為溫度恒定的固定邊界，采用非平衡外推格式.

在計算區(qū)域內(nèi)放置一個或多個半徑為R的晶核，并選擇最優(yōu)方向生長.通過式(4)計算截面平衡成分，式(5)計算固液相流動，凝固場中濃度分布發(fā)生改變通過式(2)表示，式(1)完成相場的遷移，計算枝晶生長過程.當完成一個時間步長的傳輸計算后，由于自然對流作用下枝晶生長還需考慮排出溶質(zhì)和釋放潛熱所造成的溫度和濃度分布不均勻，導(dǎo)致固相分數(shù)、液相濃度增加，若界面網(wǎng)格未完全凝固，則利用式(7)將釋放出的溶質(zhì)增量Δc加到網(wǎng)格中剩余液相溶質(zhì)分布函數(shù)，利用式(8)將釋放的溫度增量ΔT加到同網(wǎng)格的溫度分布函數(shù)，若完全凝固，則將Δc和ΔT加到相鄰網(wǎng)格液相溶質(zhì)分布函數(shù)中，獲得每個時間步長的流場、溫度場以及濃度場分布，其中流場的固液界面使用反彈格式(5、6)實現(xiàn)枝晶對速度分布影響.此模型包含了耦合液相流動、溶質(zhì)和熱量傳輸與枝晶生長的物理機制，通過重復(fù)上述步驟直到模擬結(jié)束就可實現(xiàn)流場計算和相場模擬的耦合.

2 模擬結(jié)果與分析

2.1 LBM模擬自然對流驗證

在枝晶生長過程中，流動促進溶質(zhì)和熱量的傳輸，對濃度、溫度分布產(chǎn)生影響，流場直接影響耦合模型和算法的正確性，因此需要對計算LBM模型進行驗證.封閉方腔是一個水平放置的二維腔體，如圖2所示，方腔左壁溫度高于右壁溫度并保持恒定，上下兩壁設(shè)置為絕熱，高為H，方腔左壁被加熱至Th，右壁被冷卻至Tc并保持恒定，Th>Tc.以方腔左下角作為坐標原點，u的方向為x軸方向，重力的反方向為y軸方向.其中兩個最基本無量綱參數(shù)為瑞利數(shù)Ra和普朗特數(shù)Pr，Ra描述方腔中自然對流和換熱強度參數(shù)，表示自然對流的強度，Ra越大，自然對流強度越大，普朗特數(shù)記為Pr=v/α.平均Nu數(shù)和Ra之間存在著冪律關(guān)系：Nu=a(Ra)b，式中a=0.142和b=0.299.

圖2 封閉方腔示意圖Fig.2 Schematic diagram of closed square cavity

在開始階段，不論是否存在對流，枝晶生長速度都很大，如圖3所示，隨著枝晶不斷向界面前沿釋放溶質(zhì)，導(dǎo)致溶質(zhì)富集，各個尖端生長速度下降，擴散達到平衡時，速度趨于穩(wěn)定，穩(wěn)態(tài)的速度與LGK的解析值非常接近，因此，應(yīng)用LBM模型和算法計算自然對流是可行的.

圖3 PF-LBM模型與LGK模型理論值對比

2.2 自然對流對單晶粒生長的影響

將擇優(yōu)取向與水平夾角為0度的晶核放置在計算區(qū)域中心，設(shè)定半徑為R=10 mm，初始過冷度ΔT0=10 K、溶質(zhì)初始c0=0.45 wt.%，并以枝晶中心為原點建立平面直角坐標系，將兩條數(shù)軸分別置于水平位置與垂直位置，將枝晶分為4個象限，如圖4a所示.枝晶生長是一種動態(tài)過程，需要考慮自然對流對枝晶的影響.在相變過程中，因潛熱釋放和溶質(zhì)再分配會使液相中溶質(zhì)濃度和溫度分布不均勻，導(dǎo)致凝固場中產(chǎn)生溫度梯度、濃度梯度，使流體密度不均勻，從而在重力作用下引起流動，產(chǎn)生漩渦.枝晶生長初期，自然對流強度較弱，液體流動方向由方腔底部繞枝晶至頂部，由于枝晶的生長規(guī)律以及邊界條件的限制，流場出現(xiàn)左右對稱的兩個渦流，如圖4a1和圖4b1所示，此時枝晶尚不發(fā)達，枝晶整體較為對稱.隨著潛熱釋放以及溶質(zhì)排出，界面溫度和濃度升高，流動達到一定強度時，產(chǎn)生四個漩渦，四個垂直方向生長的主枝晶臂均逐漸變粗，枝晶對稱性遭到破壞，垂直向上的一次枝晶臂受到抑制，垂直向下的一次枝晶臂生長受到促進，水平方向一次枝晶臂與純擴散條件下的枝晶生長趨勢近似，如圖4a2和圖4b2所示.

圖4 自然對流條件下,不同時刻單晶粒相場、溶質(zhì)場演變 Fig.4 The evolution of single dendrite phase field and concentration field at different times under the condition of natural convection

枝晶臂尖端速度和濃度隨時間變化曲線如圖5所示.開始階段，晶粒處于一個低濃度、低溫度的情況下，容易獲得較大的生長動力，一次枝晶臂的尖端以較大速度開始生長.同時，隨著主枝晶臂向固液界面處延伸，枝晶排出的溶質(zhì)進行再分配，潛熱的釋放富集在固/液界面，導(dǎo)致溶質(zhì)成分增加，枝晶前沿液相濃度也隨之升高，并使固相溶質(zhì)濃度低于液相初始濃度，過冷度隨著液相濃度增加而降低，枝晶生長驅(qū)動力與開始時相比有明顯降低，生長速度也明顯下降.一段時間后，固相初始溫度與液相初始溫度基本達到平衡狀態(tài)，當t=0.52×10-3s時，生長驅(qū)動力基本保持在相對穩(wěn)定區(qū)間，枝晶尖端生長速度基本達到穩(wěn)定值.由圖5a可知，垂直向下的一次枝晶臂尖端生長速度最高，達到0.021×10-3m/s，垂直向上的一次枝晶臂尖端生長速度最低，為0.19×10-3m/s，水平方向一次枝晶臂生長速度相較于純擴散條件下的情況略高一些.在此過程中，液相流動將溶質(zhì)從方腔底部沖刷至頂部，垂直向下一次枝晶臂生長速度快，排出溶質(zhì)多，富集主要集中在枝晶前沿，即使液相流動會將溶質(zhì)從底部傳輸?shù)巾敳?，垂直向下枝晶臂生長快，前沿富集的溶質(zhì)還未傳送出去，由圖5b可知，垂直向下一次枝晶臂尖端穩(wěn)態(tài)濃度依舊明顯高于垂直向上一次枝晶臂.

其他物性參數(shù)保持不變，加入起伏結(jié)構(gòu)和能量起伏觀察枝晶生長過程.觀察發(fā)現(xiàn)，起伏結(jié)構(gòu)和能量起伏為晶核形成提供了能量，凝固界面在起伏的作用下失去穩(wěn)定性，形成二次或更高次枝晶臂.無起伏時，一次枝晶臂表面光滑，未出現(xiàn)側(cè)向分枝，加入起伏后，一次枝晶臂溶質(zhì)含量最低，與過冷液體界面之間穩(wěn)定性降低，一次枝晶臂上出現(xiàn)凸起，形成側(cè)向分枝，隨著凝固進行,二次枝晶臂變粗，間距增大.

圖5 單晶粒一次枝晶臂尖端速度，濃度對比圖Fig.5 Comparison diagram of single dendrite primary dendrite arm tip velocity and concentration

整個模擬過程中枝晶外部區(qū)域受自然對流影響顯著，渦流帶動溶質(zhì)擴散，導(dǎo)致上游區(qū)域枝晶臂尖端固/液界面溶質(zhì)濃度降低，釋放出的溶質(zhì)增多，濃度略高于下游區(qū)域枝晶臂尖端濃度，如圖6b和圖6d所示.由于整個枝晶曲率基本相同，可知上游區(qū)域枝晶臂總過冷度大于下游區(qū)域，過冷度較小處界面前沿排出的溶質(zhì)和熱量傳輸至遠離界面區(qū)域時間較長，上游區(qū)域枝晶生長驅(qū)動力更大，時間越久，枝晶的不對稱性越明顯，由圖6a和圖6c可知，3、4象限內(nèi)枝晶臂尖端生長速度大于1、2象限內(nèi)枝晶臂，其中垂直向下一次枝晶臂與3、4象限內(nèi)二次枝晶臂生長速度近似，垂直向上的一次枝晶臂與1、2象限內(nèi)二次枝晶臂生長速度近似.

圖6 單晶粒二次枝晶臂尖端生長速度、濃度對比圖Fig.6 Comparison diagram of growth speed and concentration at the tip of single grain secondary dendrite arm

2.3 自然對流情況下各向異性強度對單晶粒生長的影響

各向異性強度是指界面表面張力、界面厚度以及界面動力學的各向異性程度，記為γ.圖7為自然對流情況下t=0.1 s時，各向異性強度分別為0.3、0.05、0.07時枝晶生長形貌.在自然對流存在且各向異性強度較大時，枝晶沿著〈100〉方向生長，隨著界面前沿愈加不穩(wěn)定，枝晶周圍溫度梯度、熱擴散速度相繼減小，表面張力以及液膜厚度等相關(guān)參數(shù)使枝晶周圍擾動變明顯，枝晶尺寸變大，枝晶干變細長，內(nèi)部凹陷和尖端尖銳程度變明顯，二次枝晶臂間距增大，幅值減小，根部出現(xiàn)“頸縮”現(xiàn)象.自然對流增強了各向異性強度對枝晶形貌的影響，各向異性強度為0.05時出現(xiàn)明顯“頸縮”現(xiàn)象.

圖7 t=0.1 s各向異性強度不同時單晶粒相場、溶質(zhì)場示意圖

由于下游區(qū)域溶質(zhì)來不及擴散，大量溶質(zhì)富集在枝晶尖端，下游區(qū)域枝晶臂尖端凝固，析出溶質(zhì)濃度變高，遠離枝晶臂的液相區(qū)域溶質(zhì)濃度幾乎不受流動影響，溶質(zhì)濃度近似，直到溶質(zhì)擴散和對流遷移達到動態(tài)平衡時，枝晶尖端溶質(zhì)以及濃度也趨于穩(wěn)定，上游枝晶生長受到促進，生長速度較為平緩，下游枝晶受到抑制，生長速度有明顯下降趨勢，如圖8所示.各向異性強度在0.01～0.03時枝晶速率斜率最大，生長速度變化較大，取值在0.03～0.07時斜率趨于一致，生長速度基本相同，取值在0.07之后，枝晶形貌出現(xiàn)失真現(xiàn)象，與理論值的變化不符，如圖9所示.

圖8 各向異性強度不同時枝晶生長速度、濃度對比圖

圖9 t=0.1 s時各向異性強度與生長速度對比圖Fig.9 Comparison diagram of anisotropic strength and growth rate at t = 0.1 s

2.4 自然對流對多晶粒生長的影響

多晶粒生長過程與單晶粒相比，生長情況更復(fù)雜，不僅會受到自身生長動力學和自然對流的影響，而且還受到相鄰枝晶的制約和限制.為了更好地觀察自然對流對多晶粒的影響，將5個枝晶置于計算區(qū)域內(nèi)，其中1個放置在區(qū)域中心，其余4個對稱分布在中心枝晶周圍，各向異性強度γ=0.07，設(shè)置擇優(yōu)取向與水平夾角為0°，其他條件與單晶粒生長過程模擬條件相同.初始階段，自然對流不明顯，每個枝晶生長不受相鄰枝晶的影響.當凝固時間為0.03 s時，5個晶粒生長速度大致相同，分別在固/液界面處，形成溶質(zhì)富集邊界層，色卡顯示最高體積分數(shù)到達約為5.0%，此時5個晶粒生長不受影響，溶質(zhì)富集邊界沒有接觸，如圖10所示.當t=0.07 s時，枝晶生長至尖端相互接觸，各自排出的溶質(zhì)富集融合成更大的邊界層，中心枝晶周圍溶質(zhì)大量富集，不能有效地傳輸出去，生長受到抑制，外圍4個枝晶外側(cè)區(qū)域受自然對流影響明顯，不斷排出的溶質(zhì)凝固潛熱富集在枝晶尖端，體積分數(shù)約為4.0%，固相排出的溶質(zhì)被渦流帶到液相區(qū)域，枝晶生長速度較快.隨著凝固時間推移，枝晶臂生長成熟，中間空隙部分溶質(zhì)體積分數(shù)升高至5.2%.另外，與單晶粒在自然對流下生長過程相同，計算區(qū)域上游位置的枝晶明顯比下游位置的枝晶生長速度快，不同的是，單晶粒在自然對流的情況下產(chǎn)生左右對稱的四個渦流，而多晶粒生長時枝晶外側(cè)凹陷處出現(xiàn)大小不等的多個渦流，加速溶質(zhì)向遠處液相擴散，溶質(zhì)邊界層富集的溶質(zhì)在對流的沖刷下，導(dǎo)致上游枝晶的溶質(zhì)邊界層較薄，部分溶質(zhì)被沖走促進生長.

圖10 自然對流條件下多晶粒相場和溶質(zhì)場示意圖

將枝晶擇優(yōu)取向與水平夾角設(shè)置為30°，其他條件與無角度多晶粒生長過程模擬條件相同.由于重力垂直向下，自然對流使上游區(qū)域以及枝晶右側(cè)枝晶臂生長較快，較為發(fā)達，如圖11a和圖11b所示.隨著凝固進行，由于枝晶不同的生長取向及相對位置，枝晶互相制約對方生長，晶粒在生長方向與水平夾角為30°和120°上生長速度逐漸增大，溶質(zhì)在中心枝晶周圍無法擴散，體積分數(shù)高達5.2%，枝晶間內(nèi)部區(qū)域的枝晶臂主要受到相鄰枝晶的影響，枝晶外部區(qū)域受自然對流影響較為明顯，上游枝晶的主枝晶臂右側(cè)二次枝晶臂發(fā)達于左側(cè)，甚至出現(xiàn)了三次枝晶臂，如圖11c、d所示，自然對流流向由于枝晶生長方向影響，周圍漩渦發(fā)生改變，不再呈現(xiàn)左右對稱狀態(tài).

圖11 擇優(yōu)取向與水平夾角為30°時多晶粒生長規(guī)律示意圖 Fig.11 Schematic diagram of multi grain growth law with the angle of 30 degrees between preferred orientation and horizontal

2.5 自然對流情況下各向異性強度對多晶粒生長的影響

當γ=0.03時，由于自然對流，主枝晶臂a1生長得到強化，速度加快，枝晶臂變粗壯，二次枝晶臂a3受到a2以及其他相鄰枝晶臂擠壓，溶質(zhì)富集在枝晶尖端，導(dǎo)致a3生長驅(qū)動力變小，抑制其生長，如圖12a所示.隨著各向異性強度增大，枝晶在不同方向上的生長速度明顯不同，樹枝狀逐漸明顯，枝晶臂間距越來越小，枝晶臂越來越粗壯，具體表現(xiàn)為La1>Lb1>Lc1，La2>Lb2>Lc2，La3>Lb3>Lc3，如圖12b和圖12c所示.由于枝晶間相互抑制、擠壓，外圍枝晶排出的溶質(zhì)和潛熱富集，中心溶液濃度及溫度均升高，導(dǎo)致在各向異性強度為0.03時，出現(xiàn)明顯“頸縮”現(xiàn)象.

圖12 各向異性強度不同時擇優(yōu)取向與水平夾角為30°時多晶粒溶質(zhì)場示意圖 Fig.12 Schematic diagram of multi grain concentration field with the angle of 30 degrees between preferred orientation and horizontal at different anisotropic strength

3 結(jié)論

將相場法(PF)與格子玻爾茲曼(LBM)方法相結(jié)合，建立了單晶粒和多晶粒在自然對流作用下生長的PF-LBM模型，利用該模型再現(xiàn)枝晶凝固過程，計算封閉方腔自然對流問題，結(jié)論如下:

1) 用所建立的LBM模型模擬了經(jīng)典封閉方腔中自然對流的流態(tài)，PF-LBM模型模擬結(jié)果與LGK理論值吻合較好，驗證了計算自然對流的LBM模型正確性.

2) 在單晶粒生長過程中，自然對流將熱和溶質(zhì)從上游傳輸?shù)较掠危股嫌尉植窟^冷度高于下游，促進了上游區(qū)域枝晶臂的生長，抑制了下游區(qū)域枝晶臂的生長，加入能量起伏和結(jié)構(gòu)起伏后，枝晶產(chǎn)生的不對稱性更加明顯.隨著各向異性強度增大，枝晶尖端變尖銳，主枝晶臂變細，二次枝晶臂間隙變大.

3) 在多晶粒生長過程中，晶粒除了受到自然對流的影響外，還處于相互競爭和互相抑制的狀態(tài),隨著各向異性強度增加，晶粒在擇優(yōu)方向生長速度較快，上游區(qū)域的枝晶二次枝晶臂更發(fā)達，提前出現(xiàn)“頸縮”現(xiàn)象.