基于PF-LBM模型的自然對流影響枝晶生長模擬

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

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

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

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

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

1 數學模型

1.1 枝晶生長相場模型

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

(1)

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

(2)

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

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

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

(3)

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

(4)

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

1.2 LBM流場模型

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

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

(5)

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

(6)

式中：ωi為權函數，i取0～8，cs為格子聲速.

圖1 晶格結構Fig.1 Lattice structure

LBM用于計算受對流和擴散影響的濃度場和溫度場，其分布函數演化方程與流體分布函數類似，如下式所示[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)為濃度場分布函數，hi(x,t)為溫度場分布函數，Gi、Hi分別為枝晶生長過程中濃度梯度、溫度梯度造成的溶質源和溫度源項.

1.3 數值模擬計算

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

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

2 模擬結果與分析

2.1 LBM模擬自然對流驗證

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

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

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

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

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

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

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

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

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

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

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

圖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)定，枝晶周圍溫度梯度、熱擴散速度相繼減小，表面張力以及液膜厚度等相關參數使枝晶周圍擾動變明顯，枝晶尺寸變大，枝晶干變細長，內部凹陷和尖端尖銳程度變明顯，二次枝晶臂間距增大，幅值減小，根部出現“頸縮”現象.自然對流增強了各向異性強度對枝晶形貌的影響，各向異性強度為0.05時出現明顯“頸縮”現象.

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

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

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

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

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

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

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

將枝晶擇優(yōu)取向與水平夾角設置為30°，其他條件與無角度多晶粒生長過程模擬條件相同.由于重力垂直向下，自然對流使上游區(qū)域以及枝晶右側枝晶臂生長較快，較為發(fā)達，如圖11a和圖11b所示.隨著凝固進行，由于枝晶不同的生長取向及相對位置，枝晶互相制約對方生長，晶粒在生長方向與水平夾角為30°和120°上生長速度逐漸增大，溶質在中心枝晶周圍無法擴散，體積分數高達5.2%，枝晶間內部區(qū)域的枝晶臂主要受到相鄰枝晶的影響，枝晶外部區(qū)域受自然對流影響較為明顯，上游枝晶的主枝晶臂右側二次枝晶臂發(fā)達于左側，甚至出現了三次枝晶臂，如圖11c、d所示，自然對流流向由于枝晶生長方向影響，周圍漩渦發(fā)生改變，不再呈現左右對稱狀態(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以及其他相鄰枝晶臂擠壓，溶質富集在枝晶尖端，導致a3生長驅動力變小，抑制其生長，如圖12a所示.隨著各向異性強度增大，枝晶在不同方向上的生長速度明顯不同，樹枝狀逐漸明顯，枝晶臂間距越來越小，枝晶臂越來越粗壯，具體表現為La1>Lb1>Lc1，La2>Lb2>Lc2，La3>Lb3>Lc3，如圖12b和圖12c所示.由于枝晶間相互抑制、擠壓，外圍枝晶排出的溶質和潛熱富集，中心溶液濃度及溫度均升高，導致在各向異性強度為0.03時，出現明顯“頸縮”現象.

圖12 各向異性強度不同時擇優(yōu)取向與水平夾角為30°時多晶粒溶質場示意圖 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 結論

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

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

2) 在單晶粒生長過程中，自然對流將熱和溶質從上游傳輸到下游，使上游局部過冷度高于下游，促進了上游區(qū)域枝晶臂的生長，抑制了下游區(qū)域枝晶臂的生長，加入能量起伏和結構起伏后，枝晶產生的不對稱性更加明顯.隨著各向異性強度增大，枝晶尖端變尖銳，主枝晶臂變細，二次枝晶臂間隙變大.

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