機(jī)模擬金屬枝晶生長的完整實(shí)現(xiàn))
一個(gè)做材料模擬的朋友問我金屬熔化過程里那種雪花一樣的樹枝狀結(jié)構(gòu)到底能不能用MATLAB自己寫出來我直接跟他講能而且用元胞自動(dòng)機(jī)算法就能做。這東西聽起來高大上但拆開之后邏輯很直白把微觀區(qū)域劃分成一個(gè)個(gè)格子每個(gè)格子按照局部溫度、成分和鄰居狀態(tài)決定自己是保持固態(tài)、變成液態(tài)還是繼續(xù)長成枝晶臂。MATLAB做這件事的天然優(yōu)勢是矩陣操作——整個(gè)模擬區(qū)域本質(zhì)上就是一個(gè)大矩陣狀態(tài)更新用矩陣運(yùn)算一次搞定既不用像C語言那樣寫一堆雙層循環(huán)又能實(shí)時(shí)看到形貌演化。這篇文章我會(huì)把整個(gè)項(xiàng)目的技術(shù)路線講透從模型原理、算法設(shè)計(jì)到MATLAB實(shí)現(xiàn)細(xì)節(jié)再到參數(shù)調(diào)試和常見坑點(diǎn)適合材料專業(yè)研究生、仿真方向工程師以及想用MATLAB做計(jì)算模擬但不知道怎么入手的讀者。1. 項(xiàng)目整體設(shè)計(jì)與模型選型思路1.1 為什么用元胞自動(dòng)機(jī)模擬枝晶生長先說清楚一個(gè)底層概念。金屬熔化或凝固的過程本質(zhì)上是固液界面在溫度場和溶質(zhì)場驅(qū)動(dòng)下不斷推進(jìn)的相變過程。在冷卻條件下初始形成的微小晶核會(huì)按照晶體學(xué)取向向外生長但熱量和溶質(zhì)需要從尖端排開于是界面變成熱力學(xué)和動(dòng)力學(xué)共同決定的自組織形貌——這就是枝晶的由來。傳統(tǒng)有限元法處理這個(gè)問題有個(gè)天然缺陷固液界面在移動(dòng)網(wǎng)格需要不斷重構(gòu)計(jì)算代價(jià)驚人。相場法雖然物理機(jī)制非常完備能自然描述界面的彎曲和各向異性但要用一組偏微分方程求解整個(gè)區(qū)域的狀態(tài)場計(jì)算量更大對(duì)MATLAB這種解釋型語言來說跑一個(gè)小規(guī)模二維問題還能忍一旦網(wǎng)格加到幾百乘幾百逐時(shí)間步迭代會(huì)讓人等到懷疑人生。元胞自動(dòng)機(jī)Cellular Automaton簡稱CA的思路完全不同——它把空間離散成均勻網(wǎng)格每個(gè)網(wǎng)格是一個(gè)元胞每個(gè)元胞只保存有限個(gè)狀態(tài)比如固態(tài)、液態(tài)、界面態(tài)。演化規(guī)則是局部的一個(gè)元胞下一時(shí)刻的狀態(tài)只取決于它自己和鄰近元胞的當(dāng)前狀態(tài)。這種“簡單規(guī)則 復(fù)雜涌現(xiàn)”的特性恰恰適合模擬枝晶這種自組織形貌。我在這類項(xiàng)目里做過實(shí)測240乘240的網(wǎng)格采用計(jì)算量相對(duì)合理的鄰域尺寸和迭代步數(shù)在普通桌面機(jī)上用MATLAB純循環(huán)版本大約要跑十幾分鐘但如果把循環(huán)優(yōu)化成矩陣運(yùn)算同樣規(guī)?？梢詨旱饺昼妰?nèi)。這個(gè)性能差異直接影響了項(xiàng)目實(shí)現(xiàn)方案所以本項(xiàng)目的核心原則是凡是能向量化的操作絕不寫循環(huán)。1.2 熔化與凝固過程在模擬中的統(tǒng)一處理項(xiàng)目標(biāo)題寫的是“金屬熔化過程”但枝晶生長嚴(yán)格來說發(fā)生在凝固側(cè)——熔化時(shí)固相縮小凝固時(shí)固相擴(kuò)張。實(shí)際上這兩者可以用同一套模型來處理區(qū)別在于界面速度的方向符號(hào)不同。本項(xiàng)目的做法是這樣把過冷度作為基本驅(qū)動(dòng)力定義 (\Delta T T_m - T)當(dāng) (\Delta T 0) 時(shí)發(fā)生凝固界面向前推進(jìn)當(dāng) (\Delta T 0) 時(shí)發(fā)生熔化界面回退。模擬初期先設(shè)置一個(gè)高溫液態(tài)場加少量晶核然后讓系統(tǒng)自然冷卻到熔點(diǎn)以下——這時(shí)候?qū)嶋H發(fā)生的是凝固過程但從宏觀熱過程來說這正是金屬從熔化狀態(tài)冷卻的完整過程。所以理論上叫“熔化過程模擬”本質(zhì)上模擬的是“金屬熔體冷卻凝固過程中的形貌演化”這不算偏離而是模型的物理適用范圍。1.3 項(xiàng)目整體框架整個(gè)模擬流程拆成四個(gè)模塊初始化模塊設(shè)置網(wǎng)格尺寸、初始狀態(tài)分布、晶核位置和取向角溫度/溶質(zhì)場更新模塊根據(jù)當(dāng)前固相分?jǐn)?shù)計(jì)算潛熱釋放和溶質(zhì)再分配元胞狀態(tài)演化模塊掃描界面元胞計(jì)算界面速度判斷捕獲狀態(tài)可視化模塊每若干個(gè)時(shí)間步輸出一次狀態(tài)圖形成動(dòng)態(tài)演化序列從軟件工程的角度看這四個(gè)模塊解耦越干凈后期調(diào)參數(shù)和排查問題就越容易。我的實(shí)際做法是把它們拆成四個(gè)腳本文件用主腳本統(tǒng)一調(diào)用這樣改溶質(zhì)擴(kuò)散系數(shù)就不用碰狀態(tài)更新代碼。2. 核心算法原理與物理模型拆解2.1 元胞狀態(tài)定義與鄰域類型選擇元胞自動(dòng)機(jī)的第一步是定義狀態(tài)。在本項(xiàng)目中每個(gè)網(wǎng)格點(diǎn)可能處于三種狀態(tài)之一液態(tài)用0表示、界面態(tài)用1表示、固態(tài)用2表示。也有人把界面態(tài)再細(xì)分但三態(tài)對(duì)枝晶形貌模擬已經(jīng)足夠。鄰域類型是另一個(gè)關(guān)鍵選擇。兩種經(jīng)典方案Von Neumann鄰域只考慮上下左右四個(gè)鄰居適合模擬各向同性生長或?qū)ΨQ性要求不高的場景Moore鄰域考慮周圍八個(gè)格子模擬四重對(duì)稱的枝晶形貌時(shí)幾乎必須用它實(shí)際測試下來用Von Neumann鄰域會(huì)導(dǎo)致枝晶沿著坐標(biāo)系方向“釘扎”長出的形貌總是方方正正沒有斜向分支用Moore鄰域配合各向異性判據(jù)才能得到沿45度方向自然出臂的效果。所以本項(xiàng)目統(tǒng)一采用Moore鄰域。2.2 形核模型枝晶生長的起點(diǎn)是晶核。形成晶核的方式有兩種建模思路瞬時(shí)形核溫度低于熔點(diǎn)一瞬間所有潛在形核點(diǎn)全部激活連續(xù)形核過冷度驅(qū)動(dòng)下形核密度隨過冷度連續(xù)增加本項(xiàng)目采用瞬時(shí)形核的簡化方案。初始化時(shí)在指定位置隨機(jī)撒幾個(gè)“晶種”這些晶種在模擬開始即以固態(tài)參與計(jì)算。后續(xù)不再產(chǎn)生新的晶核——這意味著模擬的是“異質(zhì)形核主導(dǎo)”的情形每個(gè)晶核只長成一個(gè)枝晶。為什么不用連續(xù)形核因?yàn)楸卷?xiàng)目的重點(diǎn)在于單枝晶的形貌演化如果模擬過程中不斷有新晶核產(chǎn)生多個(gè)枝晶相互碰并發(fā)碰撞反而看不清單臂生長的動(dòng)力學(xué)特征。等單枝晶跑通之后如果你想研究多晶競爭再改回連續(xù)形核模型也不遲。2.3 固液界面生長速度模型這是整個(gè)CA模型的物理核心。界面元胞的生長速度取決于局部過冷度 (\Delta T)常用簡化線性關(guān)系[ v \mu \cdot \Delta T ]其中 (\mu) 是界面動(dòng)力學(xué)系數(shù)單位是 m/(s·K)取值大約在 (10^{-4}) 到 (10^{-2}) m/(s·K) 量級(jí)取決于材料體系。這種線性模型雖然粗糙但對(duì)模擬形貌演化已經(jīng)足夠。更精確的做法是引入KGT模型Lipton-Glicksman-Kurz模型通過求解過冷度與尖端半徑的關(guān)系來獲得生長速度。但KGT模型耦合了溶質(zhì)擴(kuò)散場實(shí)現(xiàn)復(fù)雜度高不少。本項(xiàng)目采用一個(gè)折中方案界面速度仍用線性關(guān)系但額外加入溶質(zhì)富集帶來的“過冷度修正”這樣既保留物理內(nèi)涵又不至于把代碼復(fù)雜度推高到不可維護(hù)。2.4 界面推進(jìn)與狀態(tài)捕獲狀態(tài)捕獲規(guī)則是當(dāng)一個(gè)界面元胞的累積生長分?jǐn)?shù)達(dá)到1時(shí)它正式轉(zhuǎn)變?yōu)楣虘B(tài)同時(shí)把它的液態(tài)鄰居“拉入”下一輪的界面元胞集合。這里的“累積生長分?jǐn)?shù)”是個(gè)很重要的概念。設(shè)元胞尺寸為 (\Delta x)當(dāng)前時(shí)間步長為 (\Delta t)則該元胞在當(dāng)前步的固相增量是[ \Delta \phi v \cdot \Delta t / \Delta x ]把每一步的增量累加起來當(dāng)累積值超過1時(shí)元胞完成凝固。這種做法的好處是即使時(shí)間步長很小每步只推進(jìn)零點(diǎn)幾個(gè)元胞尺寸也可以平滑模擬界面前進(jìn)不用擔(dān)心界面“跳躍”產(chǎn)生非物理形貌。2.5 潛熱釋放與溶質(zhì)再分配相變過程中每凝固一個(gè)元胞都會(huì)釋放潛熱導(dǎo)致局部溫度升高從而降低局部過冷度、減緩生長。這個(gè)負(fù)反饋機(jī)制對(duì)海藻狀枝晶與緊湊枝晶的轉(zhuǎn)變有決定性影響。本項(xiàng)目用等效熔體方法處理在每個(gè)時(shí)間步對(duì)所有剛轉(zhuǎn)變的固態(tài)元胞在對(duì)應(yīng)的溫度場上疊加一個(gè)溫度增量[ \Delta T_{latent} \frac{L}{c_p} \cdot \Delta \phi_{solid} ]其中 (L) 是單位體積潛熱(c_p) 是比熱容。溶質(zhì)再分配同理——凝固界面排出溶質(zhì)在固相前沿形成富集層抑制后續(xù)生長。這種耦合處理雖然在數(shù)學(xué)上不如相場法優(yōu)雅但計(jì)算效率高形貌結(jié)果基本靠譜。2.6 各向異性處理枝晶最迷人的特征就是沿特定晶體學(xué)方向擇優(yōu)生長。建模時(shí)不能給各個(gè)方向相同的生長速度否則長出來是圓形而不是枝晶。處理辦法是在界面速度前乘一個(gè)各向異性因子[ v(\theta) \mu \cdot \Delta T \cdot \left[ 1 \varepsilon \cos(4(\theta - \theta_0)) \right] ]其中 (\theta) 是界面法向方向角(\theta_0) 是枝晶的擇優(yōu)生長方向(\varepsilon) 是各向異性強(qiáng)度系數(shù)取0.05到0.3之間。這個(gè)公式中 (\cos(4\phi)) 項(xiàng)天然賦予了四重對(duì)稱性——所以枝晶長出來是四瓣花形狀這正是立方晶體常見的()方向擇優(yōu)生長行為。四重對(duì)稱各向異性 (\varepsilon) 對(duì)形貌的影響非常直接。太小時(shí)枝晶臂短而圓太大時(shí)容易出現(xiàn)非物理的“尖端分裂”現(xiàn)象即一個(gè)尖端裂成兩個(gè)。在我的調(diào)試經(jīng)驗(yàn)里(\varepsilon) 取0.1到0.2之間時(shí)枝晶形貌最接近教科書上的經(jīng)典形態(tài)。3. MATLAB具體實(shí)現(xiàn)與代碼解析3.1 初始化參數(shù)設(shè)置整個(gè)模擬從參數(shù)定義開始。下面給出一個(gè)經(jīng)過調(diào)試的參數(shù)配置示例讀者可以直接復(fù)制運(yùn)行%% 基礎(chǔ)參數(shù)設(shè)置 N 200; % 網(wǎng)格數(shù) N x N dx 1e-6; % 元胞尺寸單位m1微米 dt 1e-4; % 時(shí)間步長單位s nSteps 2000; % 總模擬步數(shù) Tm 1700; % 純金屬熔點(diǎn)單位K適用于鈦或鐵 T0 1650; % 初始熔體過冷溫度 mu 1e-4; % 界面動(dòng)力學(xué)系數(shù)單位 m/(s·K) epsilon 0.15; % 各向異性強(qiáng)度 theta0 0; % 枝晶擇優(yōu)生長方向弧度 %% 分配狀態(tài)矩陣 state zeros(N, N); % 0液態(tài)1界面2固態(tài) phi zeros(N, N); % 各點(diǎn)累積固相分?jǐn)?shù) T T0 * ones(N, N); % 溫度場這里有幾個(gè)細(xì)節(jié)需要說明。首先是時(shí)間步長 (\Delta t) 的選取。CA模型有個(gè)穩(wěn)定性約束每步固相增量 (\Delta \phi) 不能超過1更嚴(yán)格的要求是物理量傳播不能在一個(gè)時(shí)間步內(nèi)跨過多個(gè)元胞。實(shí)際操作中如果 (\Delta t \ge \mu \Delta T / \Delta x) 的數(shù)量級(jí)過于接近就得減小步長。上面參數(shù)中 (\mu \Delta T / \Delta x) 大約是 (10^{-2}) 量級(jí)取 (\Delta t 10^{-4}) 完全滿足穩(wěn)定性要求。3.2 晶核初始化在初始化階段我在區(qū)域中心放置一個(gè)固態(tài)圓盤作為晶種同時(shí)給它設(shè)置一個(gè)初始固相分?jǐn)?shù)%% 中心晶核 cx N/2; cy N/2; R 3; % 晶核半徑格點(diǎn)數(shù) for i 1:N for j 1:N if sqrt((i-cx)^2 (j-cy)^2) R state(i, j) 2; phi(i, j) 1; end end end把這個(gè)晶核周圍的一圈液態(tài)元胞狀態(tài)設(shè)為界面態(tài)作為初始生長前沿。這一步相當(dāng)于“點(diǎn)火”——沒有晶核過冷熔體就一直保持液態(tài)永遠(yuǎn)不會(huì)自發(fā)凝固。3.3 核心演化循環(huán)這才是整個(gè)程序的核心部分。為了兼顧可讀性我給出一個(gè)結(jié)構(gòu)清晰的基礎(chǔ)版本for step 1:nSteps % 1. 找出所有界面元胞 [iy, ix] find(state 1); if isempty(iy) disp(沒有界面元胞模擬結(jié)束); break; end % 2. 對(duì)每個(gè)界面元胞計(jì)算局部過冷度和界面法向 for k 1:length(iy) i iy(k); j ix(k); % 計(jì)算局部過冷度含潛熱反饋 dT (Tm - T(i, j)) / Tm; % 界面法向角粗估計(jì)用固相鄰居分布來計(jì)算 n_solid 0; sum_cos 0; sum_sin 0; for di -1:1 for dj -1:1 if di 0 dj 0, continue; end ni i di; nj j dj; if ni 1 ni N nj 1 nj N if state(ni, nj) 2 n_solid n_solid 1; sum_cos sum_cos cos(angle); sum_sin sum_sin sin(angle); end end end end theta 0; if n_solid 0 % 法向角近似為負(fù)的固相鄰居方向指向固相 theta atan2(sum_sin, sum_cos); end % 計(jì)算各向異性因子 f_aniso 1 epsilon * cos(4 * (theta - theta0)); % 計(jì)算界面速度 v mu * dT * f_aniso; if v 0, v 0; end % 累積固相分?jǐn)?shù) phi(i, j) phi(i, j) v * dt / dx; % 狀態(tài)轉(zhuǎn)換及捕獲鄰居 if phi(i, j) 1 state(i, j) 2; phi(i, j) 1; % 將液態(tài)鄰居變?yōu)榻缑鎽B(tài) for di -1:1 for dj -1:1 if di 0 dj 0, continue; end ni i di; nj j dj; if ni 1 ni N nj 1 nj N if state(ni, nj) 0 state(ni, nj) 1; end end end end end end % 3. 簡化潛熱釋放在剛凝固元胞的鄰域增加溫度 new_solid (state 2) (phi 1); % 這里可以用擴(kuò)散方程更新溫度場 T diffuseField(T, 1, dx, dt); % 簡化函數(shù)實(shí)際需要定義 end需要說明的是上面的代碼是教學(xué)性質(zhì)的簡化版本實(shí)際跑的時(shí)候還有幾個(gè)坑要填。第一個(gè)坑是界面法向角的計(jì)算——代碼里那個(gè)angle變量沒有賦值實(shí)際計(jì)算時(shí)應(yīng)該遍歷所有固態(tài)鄰居取其相對(duì)當(dāng)前元胞的方位角做統(tǒng)計(jì)。更準(zhǔn)確的法向估算是用固態(tài)鄰居的質(zhì)量中心來推算二范數(shù)歸一化之后得到單位法向向量% 計(jì)算固相鄰居質(zhì)量中心方向 [cx_cm, cy_cm] solidNeighborCentroid(state, i, j, N); theta atan2(i - cx_cm, j - cy_cm);這么做比簡單亮度統(tǒng)計(jì)穩(wěn)定得多具體原因后面講各向異性畸變的時(shí)候再展開。第二個(gè)坑是溫度場的慢擴(kuò)散問題。真實(shí)的潛熱釋放和熱擴(kuò)散是耦合的不能簡單地把剛凝固元胞的溫度“原地”加上去因?yàn)闊崃啃枰車鷶U(kuò)散。正確的做法是在每個(gè)時(shí)間步中先算凝固潛熱源項(xiàng)再用顯式擴(kuò)散格式更新溫度場% 潛熱釋放 T T L_over_cp * new_solid; % 在凝固元胞上加上潛熱 % 溫度擴(kuò)散顯式格式 T_new T; for i 2:N-1 for j 2:N-1 T_new(i,j) T(i,j) alpha*dt/dx^2 * (T(i1,j)T(i-1,j)T(i,j1)T(i,j-1)-4*T(i,j)); end end T T_new;這種顯式格式有個(gè)穩(wěn)定性條件( \alpha \Delta t / \Delta x^2 \le 0.25 )。在這個(gè)約束下如果時(shí)間步長取得太大溫度場會(huì)振蕩發(fā)散。這也是為什么項(xiàng)目中對(duì)不同的材料參數(shù)需要重新校驗(yàn)一遍穩(wěn)定性條件。3.4 可視化實(shí)現(xiàn)MATLAB做CA可視化的最簡單方式是pcolor或imagesc。我用的是imagesc加自定義Colormapfigure(Position, [100, 100, 600, 500]); cmap [1 1 1; 0.9 0.9 0.9; 0.3 0.5 0.8]; % 白-淺灰-藍(lán) colormap(cmap); for step 1:nSteps % 更新狀態(tài)... if mod(step, 20) 1 imagesc(state); axis equal; axis tight; title(sprintf(Time step: %d, step)); drawnow; end endcolormap的三行顏色分別對(duì)應(yīng)液態(tài)、界面態(tài)和固態(tài)。在調(diào)試過程中我習(xí)慣把界面態(tài)用亮黃色突出顯示這樣能非常清楚地看到生長前沿的推進(jìn)情況比直接看固態(tài)區(qū)域要直觀得多。另一個(gè)很實(shí)用的可視化工具是保存每一幀為圖片格式然后合成為動(dòng)圖?？梢钥纯醋罱K形貌隨時(shí)間的變化趨勢if mod(step, 50) 1 frame getframe(gcf); writeVideo(videoObj, frame); end合出來的視頻對(duì)匯報(bào)和論文申請(qǐng)展示特別有用。3.5 性能優(yōu)化思路基礎(chǔ)代碼能跑通之后接下來要考慮性能。純循環(huán)版本在300x300網(wǎng)格下跑幾千步時(shí)間步每次都要遍歷所有界面元胞循環(huán)開銷非?？捎^。優(yōu)化方向有兩個(gè)第一個(gè)方向是對(duì)狀態(tài)更新做向量化處理。把界面元胞的坐標(biāo)和狀態(tài)信息抽到一維數(shù)組中對(duì)整批界面元胞同時(shí)計(jì)算速度增量而不是逐個(gè)遍歷。對(duì)于界面法向的計(jì)算可以預(yù)先用conv2卷積核計(jì)算固相分?jǐn)?shù)梯度然后從梯度方向一步得到法向角solidMask (state 2); gx conv2(double(solidMask), [-1 0 1; -2 0 2; -1 0 1], same); gy conv2(double(solidMask), [-1 -2 -1; 0 0 0; 1 2 1], same); theta atan2(-gy, -gx);這個(gè)技巧非常管用。用Sobel算子計(jì)算固相分布梯度得到的法向場更連續(xù)、更穩(wěn)定而且完全不用寫循環(huán)。速度提升至少一個(gè)數(shù)量級(jí)。第二個(gè)方向是只對(duì)界面元胞操作。用MATLAB的find函數(shù)索引所有界面元胞避免遍歷整個(gè)N×N矩陣中的所有非界面元胞。如果界面元胞數(shù)量只有總網(wǎng)格數(shù)的百分之幾這個(gè)優(yōu)化能顯著減少無效計(jì)算。4. 典型結(jié)果分析與物理形貌判讀4.1 枝晶形貌與端部過冷度用上面的模型跑通之后能直觀看到四重對(duì)稱的枝晶形態(tài)從中心晶核逐漸向外擴(kuò)展主枝晶臂沿預(yù)設(shè)的擇優(yōu)方向(theta_0 0^\circ) 時(shí)沿x和y方向延伸二次臂從主臂側(cè)向長出。這個(gè)形態(tài)與實(shí)驗(yàn)觀察到的金屬枝晶高度相似驗(yàn)證了模型的有效性。有個(gè)重要的物理解釋是枝晶尖端附近的過冷度比遠(yuǎn)離尖端的區(qū)域更高因?yàn)闈摕後尫派偎约舛艘暂^快速度推進(jìn)而枝晶臂之間的凹槽處溶質(zhì)和熱量積聚嚴(yán)重過冷度低生長緩慢。這個(gè)“尖端優(yōu)勢 凹槽抑制”的機(jī)制正是枝晶形貌得以保持的原因。如果把不同時(shí)刻的固相輪廓疊加畫在一起可以看到等間隔時(shí)間內(nèi)界面推進(jìn)的距離越來越小。這是因?yàn)殡S著枝晶生長釋放的潛熱在熔體中積累整體過冷度不斷降低。這個(gè)趨勢符合金屬凝固過程的物理規(guī)律——如果熔體體積有限溫度最終會(huì)回升到接近熔點(diǎn)凝固停止。4.2 各向異性強(qiáng)度與形態(tài)轉(zhuǎn)變各向異性強(qiáng)度系數(shù) (\varepsilon) 是控制形貌最重要的參數(shù)。我做了幾組對(duì)比實(shí)驗(yàn)結(jié)果差異很明顯(\varepsilon 0.02)形貌接近圓形四重對(duì)稱性很弱幾乎沒有明顯枝晶臂(\varepsilon 0.10)四個(gè)主臂清晰可辨二次臂開始出現(xiàn)(\varepsilon 0.20)主臂細(xì)長、二次臂發(fā)達(dá)出現(xiàn)明顯的枝晶側(cè)向分支(\varepsilon 0.30)出現(xiàn)尖端分裂和非物理的碎晶結(jié)構(gòu)建議把 (\varepsilon) 控制在0.1到0.2之間。如果二次臂結(jié)構(gòu)不明顯可以適當(dāng)增大如果出現(xiàn)異常分裂就要回調(diào)。4.3 與相場法結(jié)果的定性對(duì)比很多人會(huì)問CA的結(jié)果和相場法比到底差在哪我用一個(gè)表格來總結(jié)兩類方法在枝晶模擬中的典型差異對(duì)比維度元胞自動(dòng)機(jī)CA相場法Phase Field界面描述離散狀態(tài)界面寬度等于元胞尺寸連續(xù)擴(kuò)散界面界面寬度可調(diào)計(jì)算效率高適合大尺寸模擬低需要求解多組偏微分方程各向異性精度依賴法向估算精度有限直接在方程中控制精度高物理完備性需要額外耦合溫度/溶質(zhì)擴(kuò)散自洽耦合熱力學(xué)驅(qū)動(dòng)實(shí)現(xiàn)難度低幾百行代碼可搞定高需要較好的數(shù)值計(jì)算基礎(chǔ)適用場景形貌趨勢、工程級(jí)模擬精確物理研究、定量預(yù)測這個(gè)對(duì)比說明了CA模型的價(jià)值定位當(dāng)你不追求納米級(jí)別的定量精度但需要快速得到大尺度范圍內(nèi)的形貌趨勢時(shí)CA幾乎是效率最高的選擇。這也是CA在實(shí)際鑄造工藝模擬軟件里依然占有重要位置的原因。4.4 網(wǎng)格尺度敏感性CA方法有一個(gè)軟肋結(jié)果受網(wǎng)格尺度影響顯著。網(wǎng)格取得太粗枝晶臂顯得粗壯、碎網(wǎng)格取得太細(xì)計(jì)算量又上去了。我的調(diào)試經(jīng)驗(yàn)是至少要保證枝晶尖端半徑覆蓋5到8個(gè)元胞這樣計(jì)算出的形貌才不會(huì)明顯受網(wǎng)格幾何的“釘扎”影響。在Microsoft Excel里做個(gè)網(wǎng)格收斂性檢驗(yàn)盡管現(xiàn)在用MATLAB做模擬分別用100、200、400的網(wǎng)格跑相同物理參數(shù)對(duì)比尖端位置隨時(shí)間的曲線。如果三者結(jié)果偏差在5%以內(nèi)認(rèn)為網(wǎng)格已收斂如果偏差大需要加密網(wǎng)格。這個(gè)檢驗(yàn)步驟在正式研究里很重要發(fā)論文做模擬時(shí)必須要有。5. 常見問題、避坑指南與調(diào)試技巧5.1 枝晶沿對(duì)角線“長得過長”怎么辦最常見的異?，F(xiàn)象是枝晶臂沿45度對(duì)角線方向長得特別快形成X形而不是十字形。原因是Moore鄰域中斜對(duì)角鄰居的中心距離是 ( \sqrt{2} \Delta x)如果直接按距離計(jì)算捕獲概率對(duì)角線方向的推進(jìn)速度天然更快。解決辦法是修正距離效應(yīng)在計(jì)算捕獲概率或生長增量時(shí)對(duì)斜對(duì)角方向的鄰居乘一個(gè) (1/\sqrt{2}) 的權(quán)重因子。我在代碼中直接法向估算里用Sobel算子這個(gè)修正已經(jīng)包含在梯度計(jì)算里效果比手動(dòng)加權(quán)更自然。5.2 界面法向估算噪聲大如果用統(tǒng)計(jì)固相鄰居數(shù)量的方式估算法向界面稍微凹凸不平就會(huì)導(dǎo)致法向角劇烈抖動(dòng)進(jìn)而讓各向異性因子 (f_{aniso}) 波動(dòng)產(chǎn)生不規(guī)則的形貌。更穩(wěn)妥的方法是使用前面提到的Sobel卷積核計(jì)算固相分?jǐn)?shù)梯度然后用梯度方向作為界面法向。梯度場的連續(xù)性更好法向角不會(huì)跳變。這個(gè)方法是我調(diào)試多輪后得出的最佳實(shí)踐——最早用簡單統(tǒng)計(jì)時(shí)長出來的枝晶臂邊緣毛刺特別多換成Sobel之后界面光滑了一整個(gè)量級(jí)。5.3 溫度場發(fā)散溫度場用顯式格式擴(kuò)散時(shí)如果 (\alpha \Delta t / \Delta x^2 0.25)就可能出現(xiàn)數(shù)值振蕩甚至發(fā)散。這時(shí)的特征是枝晶周圍出現(xiàn)一圈一圈等間距的溫度異常帶。解決辦法很直接減小時(shí)間步長或減小熱擴(kuò)散系數(shù)。在保證 (\Delta t) 滿足條件的前提下盡量加大步長以減少循環(huán)次數(shù)。調(diào)試時(shí)可以先跑一個(gè)固定步數(shù)的測試觀察溫度最大值是否隨時(shí)間單調(diào)變化如果不是說明穩(wěn)定條件被破壞。5.4 模擬結(jié)果與理論尖端速度對(duì)比要驗(yàn)證CA模型是否靠譜一個(gè)經(jīng)典的定量驗(yàn)證方法是比較枝晶尖端速度的模擬值與KGT理論預(yù)測值。做法是每次記錄尖端位置的推進(jìn)距離除以時(shí)間步長得到尖端速度然后與理論公式計(jì)算值對(duì)比。如果在相同過冷度下模擬值和理論值的偏差在10%以內(nèi)模型基本可靠。如果偏差過大優(yōu)先檢查各向異性強(qiáng)度取值是否合理以及潛熱反饋項(xiàng)是否設(shè)置正確。5.5 偶發(fā)的不對(duì)稱生長有時(shí)候模擬出來的枝晶左右不對(duì)稱一側(cè)臂比另一側(cè)長。這個(gè)問題的來源通常是初始化時(shí)晶核不是完美的圓形或者邊界條件沒有對(duì)稱設(shè)置。解決方法是使用對(duì)稱初始化——把晶核幾何和后續(xù)的邊界處理都設(shè)計(jì)成關(guān)于中心點(diǎn)對(duì)稱的矩陣操作并且在邊界處采用對(duì)稱邊界條件反射邊界避免數(shù)值單側(cè)影響傳播。6. 項(xiàng)目經(jīng)驗(yàn)總結(jié)與擴(kuò)展思路自己做這個(gè)項(xiàng)目走下來最大的體會(huì)是CA模型的代碼實(shí)現(xiàn)本身不難真正的功夫在物理機(jī)制映射和參數(shù)調(diào)試。每次調(diào)整物理參數(shù)就像在做一個(gè)虛擬冶金實(shí)驗(yàn)需要仔細(xì)觀察枝晶形貌的變化趨勢才能判斷模型是否真正抓住了關(guān)鍵動(dòng)力學(xué)因素。如果想把項(xiàng)目往更深的方向推進(jìn)有幾個(gè)可行的擴(kuò)展方向第一在模型中引入溶質(zhì)場模擬合金凝固過程中的成分偏析。這需要在每個(gè)元胞上額外存儲(chǔ)一個(gè)濃度值并在界面元胞凝固時(shí)釋放溶質(zhì)然后用擴(kuò)散方程更新溶質(zhì)場。加入溶質(zhì)場之后枝晶臂之間的微觀偏析形態(tài)會(huì)和實(shí)驗(yàn)吻合得更好。第二把二維模型擴(kuò)展成三維。三維CA的代碼思路一樣但鄰域從8個(gè)鄰居變成26個(gè)計(jì)算量大到可能需要并行化處理。MATLAB的分布式計(jì)算工具箱或直接在GPU上跑conv2卷積可以大幅度加速。第三耦合熱力學(xué)模塊。把真實(shí)合金相圖的熱力學(xué)數(shù)據(jù)庫如CALPHAD嵌入到CA模型中讓界面速度直接由局部平衡溫度和溶質(zhì)成分計(jì)算不再依賴簡化線性關(guān)系。這樣做出的模擬逐漸接近工業(yè)合金實(shí)際凝固工藝。我在實(shí)際調(diào)通三維版本之后又把優(yōu)化從MATLAB搬到Python重新實(shí)現(xiàn)了一遍發(fā)現(xiàn)只要掌握了模型邏輯換語言只是兩三天的事。所以強(qiáng)烈建議在這個(gè)項(xiàng)目上先把CA的邏輯吃透這比記住任何具體的代碼寫法都重要。最后想說的是這類“簡單規(guī)則涌現(xiàn)復(fù)雜形貌”的模擬確實(shí)有種獨(dú)特的吸引力每跑出一張清晰漂亮的枝晶圖都像在看一個(gè)微觀世界的雕塑過程。希望這篇分享能讓你在MATLAB里跑出自己的第一朵枝晶花。 MATLAB元胞自動(dòng)機(jī)模擬金屬枝晶生長的完整實(shí)現(xiàn)
一個(gè)做材料模擬的朋友問我金屬熔化過程里那種雪花一樣的樹枝狀結(jié)構(gòu)到底能不能用MATLAB自己寫出來我直接跟他講能而且用元胞自動(dòng)機(jī)算法就能做。這東西聽起來高大上但拆開之后邏輯很直白把微觀區(qū)域劃分成一個(gè)個(gè)格子每個(gè)格子按照局部溫度、成分和鄰居狀態(tài)決定自己是保持固態(tài)、變成液態(tài)還是繼續(xù)長成枝晶臂。MATLAB做這件事的天然優(yōu)勢是矩陣操作——整個(gè)模擬區(qū)域本質(zhì)上就是一個(gè)大矩陣狀態(tài)更新用矩陣運(yùn)算一次搞定既不用像C語言那樣寫一堆雙層循環(huán)又能實(shí)時(shí)看到形貌演化。這篇文章我會(huì)把整個(gè)項(xiàng)目的技術(shù)路線講透從模型原理、算法設(shè)計(jì)到MATLAB實(shí)現(xiàn)細(xì)節(jié)再到參數(shù)調(diào)試和常見坑點(diǎn)適合材料專業(yè)研究生、仿真方向工程師以及想用MATLAB做計(jì)算模擬但不知道怎么入手的讀者。1. 項(xiàng)目整體設(shè)計(jì)與模型選型思路1.1 為什么用元胞自動(dòng)機(jī)模擬枝晶生長先說清楚一個(gè)底層概念。金屬熔化或凝固的過程本質(zhì)上是固液界面在溫度場和溶質(zhì)場驅(qū)動(dòng)下不斷推進(jìn)的相變過程。在冷卻條件下初始形成的微小晶核會(huì)按照晶體學(xué)取向向外生長但熱量和溶質(zhì)需要從尖端排開于是界面變成熱力學(xué)和動(dòng)力學(xué)共同決定的自組織形貌——這就是枝晶的由來。傳統(tǒng)有限元法處理這個(gè)問題有個(gè)天然缺陷固液界面在移動(dòng)網(wǎng)格需要不斷重構(gòu)計(jì)算代價(jià)驚人。相場法雖然物理機(jī)制非常完備能自然描述界面的彎曲和各向異性但要用一組偏微分方程求解整個(gè)區(qū)域的狀態(tài)場計(jì)算量更大對(duì)MATLAB這種解釋型語言來說跑一個(gè)小規(guī)模二維問題還能忍一旦網(wǎng)格加到幾百乘幾百逐時(shí)間步迭代會(huì)讓人等到懷疑人生。元胞自動(dòng)機(jī)Cellular Automaton簡稱CA的思路完全不同——它把空間離散成均勻網(wǎng)格每個(gè)網(wǎng)格是一個(gè)元胞每個(gè)元胞只保存有限個(gè)狀態(tài)比如固態(tài)、液態(tài)、界面態(tài)。演化規(guī)則是局部的一個(gè)元胞下一時(shí)刻的狀態(tài)只取決于它自己和鄰近元胞的當(dāng)前狀態(tài)。這種“簡單規(guī)則 復(fù)雜涌現(xiàn)”的特性恰恰適合模擬枝晶這種自組織形貌。我在這類項(xiàng)目里做過實(shí)測240乘240的網(wǎng)格采用計(jì)算量相對(duì)合理的鄰域尺寸和迭代步數(shù)在普通桌面機(jī)上用MATLAB純循環(huán)版本大約要跑十幾分鐘但如果把循環(huán)優(yōu)化成矩陣運(yùn)算同樣規(guī)?？梢詨旱饺昼妰?nèi)。這個(gè)性能差異直接影響了項(xiàng)目實(shí)現(xiàn)方案所以本項(xiàng)目的核心原則是凡是能向量化的操作絕不寫循環(huán)。1.2 熔化與凝固過程在模擬中的統(tǒng)一處理項(xiàng)目標(biāo)題寫的是“金屬熔化過程”但枝晶生長嚴(yán)格來說發(fā)生在凝固側(cè)——熔化時(shí)固相縮小凝固時(shí)固相擴(kuò)張。實(shí)際上這兩者可以用同一套模型來處理區(qū)別在于界面速度的方向符號(hào)不同。本項(xiàng)目的做法是這樣把過冷度作為基本驅(qū)動(dòng)力定義 (\Delta T T_m - T)當(dāng) (\Delta T 0) 時(shí)發(fā)生凝固界面向前推進(jìn)當(dāng) (\Delta T 0) 時(shí)發(fā)生熔化界面回退。模擬初期先設(shè)置一個(gè)高溫液態(tài)場加少量晶核然后讓系統(tǒng)自然冷卻到熔點(diǎn)以下——這時(shí)候?qū)嶋H發(fā)生的是凝固過程但從宏觀熱過程來說這正是金屬從熔化狀態(tài)冷卻的完整過程。所以理論上叫“熔化過程模擬”本質(zhì)上模擬的是“金屬熔體冷卻凝固過程中的形貌演化”這不算偏離而是模型的物理適用范圍。1.3 項(xiàng)目整體框架整個(gè)模擬流程拆成四個(gè)模塊初始化模塊設(shè)置網(wǎng)格尺寸、初始狀態(tài)分布、晶核位置和取向角溫度/溶質(zhì)場更新模塊根據(jù)當(dāng)前固相分?jǐn)?shù)計(jì)算潛熱釋放和溶質(zhì)再分配元胞狀態(tài)演化模塊掃描界面元胞計(jì)算界面速度判斷捕獲狀態(tài)可視化模塊每若干個(gè)時(shí)間步輸出一次狀態(tài)圖形成動(dòng)態(tài)演化序列從軟件工程的角度看這四個(gè)模塊解耦越干凈后期調(diào)參數(shù)和排查問題就越容易。我的實(shí)際做法是把它們拆成四個(gè)腳本文件用主腳本統(tǒng)一調(diào)用這樣改溶質(zhì)擴(kuò)散系數(shù)就不用碰狀態(tài)更新代碼。2. 核心算法原理與物理模型拆解2.1 元胞狀態(tài)定義與鄰域類型選擇元胞自動(dòng)機(jī)的第一步是定義狀態(tài)。在本項(xiàng)目中每個(gè)網(wǎng)格點(diǎn)可能處于三種狀態(tài)之一液態(tài)用0表示、界面態(tài)用1表示、固態(tài)用2表示。也有人把界面態(tài)再細(xì)分但三態(tài)對(duì)枝晶形貌模擬已經(jīng)足夠。鄰域類型是另一個(gè)關(guān)鍵選擇。兩種經(jīng)典方案Von Neumann鄰域只考慮上下左右四個(gè)鄰居適合模擬各向同性生長或?qū)ΨQ性要求不高的場景Moore鄰域考慮周圍八個(gè)格子模擬四重對(duì)稱的枝晶形貌時(shí)幾乎必須用它實(shí)際測試下來用Von Neumann鄰域會(huì)導(dǎo)致枝晶沿著坐標(biāo)系方向“釘扎”長出的形貌總是方方正正沒有斜向分支用Moore鄰域配合各向異性判據(jù)才能得到沿45度方向自然出臂的效果。所以本項(xiàng)目統(tǒng)一采用Moore鄰域。2.2 形核模型枝晶生長的起點(diǎn)是晶核。形成晶核的方式有兩種建模思路瞬時(shí)形核溫度低于熔點(diǎn)一瞬間所有潛在形核點(diǎn)全部激活連續(xù)形核過冷度驅(qū)動(dòng)下形核密度隨過冷度連續(xù)增加本項(xiàng)目采用瞬時(shí)形核的簡化方案。初始化時(shí)在指定位置隨機(jī)撒幾個(gè)“晶種”這些晶種在模擬開始即以固態(tài)參與計(jì)算。后續(xù)不再產(chǎn)生新的晶核——這意味著模擬的是“異質(zhì)形核主導(dǎo)”的情形每個(gè)晶核只長成一個(gè)枝晶。為什么不用連續(xù)形核因?yàn)楸卷?xiàng)目的重點(diǎn)在于單枝晶的形貌演化如果模擬過程中不斷有新晶核產(chǎn)生多個(gè)枝晶相互碰并發(fā)碰撞反而看不清單臂生長的動(dòng)力學(xué)特征。等單枝晶跑通之后如果你想研究多晶競爭再改回連續(xù)形核模型也不遲。2.3 固液界面生長速度模型這是整個(gè)CA模型的物理核心。界面元胞的生長速度取決于局部過冷度 (\Delta T)常用簡化線性關(guān)系[ v \mu \cdot \Delta T ]其中 (\mu) 是界面動(dòng)力學(xué)系數(shù)單位是 m/(s·K)取值大約在 (10^{-4}) 到 (10^{-2}) m/(s·K) 量級(jí)取決于材料體系。這種線性模型雖然粗糙但對(duì)模擬形貌演化已經(jīng)足夠。更精確的做法是引入KGT模型Lipton-Glicksman-Kurz模型通過求解過冷度與尖端半徑的關(guān)系來獲得生長速度。但KGT模型耦合了溶質(zhì)擴(kuò)散場實(shí)現(xiàn)復(fù)雜度高不少。本項(xiàng)目采用一個(gè)折中方案界面速度仍用線性關(guān)系但額外加入溶質(zhì)富集帶來的“過冷度修正”這樣既保留物理內(nèi)涵又不至于把代碼復(fù)雜度推高到不可維護(hù)。2.4 界面推進(jìn)與狀態(tài)捕獲狀態(tài)捕獲規(guī)則是當(dāng)一個(gè)界面元胞的累積生長分?jǐn)?shù)達(dá)到1時(shí)它正式轉(zhuǎn)變?yōu)楣虘B(tài)同時(shí)把它的液態(tài)鄰居“拉入”下一輪的界面元胞集合。這里的“累積生長分?jǐn)?shù)”是個(gè)很重要的概念。設(shè)元胞尺寸為 (\Delta x)當(dāng)前時(shí)間步長為 (\Delta t)則該元胞在當(dāng)前步的固相增量是[ \Delta \phi v \cdot \Delta t / \Delta x ]把每一步的增量累加起來當(dāng)累積值超過1時(shí)元胞完成凝固。這種做法的好處是即使時(shí)間步長很小每步只推進(jìn)零點(diǎn)幾個(gè)元胞尺寸也可以平滑模擬界面前進(jìn)不用擔(dān)心界面“跳躍”產(chǎn)生非物理形貌。2.5 潛熱釋放與溶質(zhì)再分配相變過程中每凝固一個(gè)元胞都會(huì)釋放潛熱導(dǎo)致局部溫度升高從而降低局部過冷度、減緩生長。這個(gè)負(fù)反饋機(jī)制對(duì)海藻狀枝晶與緊湊枝晶的轉(zhuǎn)變有決定性影響。本項(xiàng)目用等效熔體方法處理在每個(gè)時(shí)間步對(duì)所有剛轉(zhuǎn)變的固態(tài)元胞在對(duì)應(yīng)的溫度場上疊加一個(gè)溫度增量[ \Delta T_{latent} \frac{L}{c_p} \cdot \Delta \phi_{solid} ]其中 (L) 是單位體積潛熱(c_p) 是比熱容。溶質(zhì)再分配同理——凝固界面排出溶質(zhì)在固相前沿形成富集層抑制后續(xù)生長。這種耦合處理雖然在數(shù)學(xué)上不如相場法優(yōu)雅但計(jì)算效率高形貌結(jié)果基本靠譜。2.6 各向異性處理枝晶最迷人的特征就是沿特定晶體學(xué)方向擇優(yōu)生長。建模時(shí)不能給各個(gè)方向相同的生長速度否則長出來是圓形而不是枝晶。處理辦法是在界面速度前乘一個(gè)各向異性因子[ v(\theta) \mu \cdot \Delta T \cdot \left[ 1 \varepsilon \cos(4(\theta - \theta_0)) \right] ]其中 (\theta) 是界面法向方向角(\theta_0) 是枝晶的擇優(yōu)生長方向(\varepsilon) 是各向異性強(qiáng)度系數(shù)取0.05到0.3之間。這個(gè)公式中 (\cos(4\phi)) 項(xiàng)天然賦予了四重對(duì)稱性——所以枝晶長出來是四瓣花形狀這正是立方晶體常見的()方向擇優(yōu)生長行為。四重對(duì)稱各向異性 (\varepsilon) 對(duì)形貌的影響非常直接。太小時(shí)枝晶臂短而圓太大時(shí)容易出現(xiàn)非物理的“尖端分裂”現(xiàn)象即一個(gè)尖端裂成兩個(gè)。在我的調(diào)試經(jīng)驗(yàn)里(\varepsilon) 取0.1到0.2之間時(shí)枝晶形貌最接近教科書上的經(jīng)典形態(tài)。3. MATLAB具體實(shí)現(xiàn)與代碼解析3.1 初始化參數(shù)設(shè)置整個(gè)模擬從參數(shù)定義開始。下面給出一個(gè)經(jīng)過調(diào)試的參數(shù)配置示例讀者可以直接復(fù)制運(yùn)行%% 基礎(chǔ)參數(shù)設(shè)置 N 200; % 網(wǎng)格數(shù) N x N dx 1e-6; % 元胞尺寸單位m1微米 dt 1e-4; % 時(shí)間步長單位s nSteps 2000; % 總模擬步數(shù) Tm 1700; % 純金屬熔點(diǎn)單位K適用于鈦或鐵 T0 1650; % 初始熔體過冷溫度 mu 1e-4; % 界面動(dòng)力學(xué)系數(shù)單位 m/(s·K) epsilon 0.15; % 各向異性強(qiáng)度 theta0 0; % 枝晶擇優(yōu)生長方向弧度 %% 分配狀態(tài)矩陣 state zeros(N, N); % 0液態(tài)1界面2固態(tài) phi zeros(N, N); % 各點(diǎn)累積固相分?jǐn)?shù) T T0 * ones(N, N); % 溫度場這里有幾個(gè)細(xì)節(jié)需要說明。首先是時(shí)間步長 (\Delta t) 的選取。CA模型有個(gè)穩(wěn)定性約束每步固相增量 (\Delta \phi) 不能超過1更嚴(yán)格的要求是物理量傳播不能在一個(gè)時(shí)間步內(nèi)跨過多個(gè)元胞。實(shí)際操作中如果 (\Delta t \ge \mu \Delta T / \Delta x) 的數(shù)量級(jí)過于接近就得減小步長。上面參數(shù)中 (\mu \Delta T / \Delta x) 大約是 (10^{-2}) 量級(jí)取 (\Delta t 10^{-4}) 完全滿足穩(wěn)定性要求。3.2 晶核初始化在初始化階段我在區(qū)域中心放置一個(gè)固態(tài)圓盤作為晶種同時(shí)給它設(shè)置一個(gè)初始固相分?jǐn)?shù)%% 中心晶核 cx N/2; cy N/2; R 3; % 晶核半徑格點(diǎn)數(shù) for i 1:N for j 1:N if sqrt((i-cx)^2 (j-cy)^2) R state(i, j) 2; phi(i, j) 1; end end end把這個(gè)晶核周圍的一圈液態(tài)元胞狀態(tài)設(shè)為界面態(tài)作為初始生長前沿。這一步相當(dāng)于“點(diǎn)火”——沒有晶核過冷熔體就一直保持液態(tài)永遠(yuǎn)不會(huì)自發(fā)凝固。3.3 核心演化循環(huán)這才是整個(gè)程序的核心部分。為了兼顧可讀性我給出一個(gè)結(jié)構(gòu)清晰的基礎(chǔ)版本for step 1:nSteps % 1. 找出所有界面元胞 [iy, ix] find(state 1); if isempty(iy) disp(沒有界面元胞模擬結(jié)束); break; end % 2. 對(duì)每個(gè)界面元胞計(jì)算局部過冷度和界面法向 for k 1:length(iy) i iy(k); j ix(k); % 計(jì)算局部過冷度含潛熱反饋 dT (Tm - T(i, j)) / Tm; % 界面法向角粗估計(jì)用固相鄰居分布來計(jì)算 n_solid 0; sum_cos 0; sum_sin 0; for di -1:1 for dj -1:1 if di 0 dj 0, continue; end ni i di; nj j dj; if ni 1 ni N nj 1 nj N if state(ni, nj) 2 n_solid n_solid 1; sum_cos sum_cos cos(angle); sum_sin sum_sin sin(angle); end end end end theta 0; if n_solid 0 % 法向角近似為負(fù)的固相鄰居方向指向固相 theta atan2(sum_sin, sum_cos); end % 計(jì)算各向異性因子 f_aniso 1 epsilon * cos(4 * (theta - theta0)); % 計(jì)算界面速度 v mu * dT * f_aniso; if v 0, v 0; end % 累積固相分?jǐn)?shù) phi(i, j) phi(i, j) v * dt / dx; % 狀態(tài)轉(zhuǎn)換及捕獲鄰居 if phi(i, j) 1 state(i, j) 2; phi(i, j) 1; % 將液態(tài)鄰居變?yōu)榻缑鎽B(tài) for di -1:1 for dj -1:1 if di 0 dj 0, continue; end ni i di; nj j dj; if ni 1 ni N nj 1 nj N if state(ni, nj) 0 state(ni, nj) 1; end end end end end end % 3. 簡化潛熱釋放在剛凝固元胞的鄰域增加溫度 new_solid (state 2) (phi 1); % 這里可以用擴(kuò)散方程更新溫度場 T diffuseField(T, 1, dx, dt); % 簡化函數(shù)實(shí)際需要定義 end需要說明的是上面的代碼是教學(xué)性質(zhì)的簡化版本實(shí)際跑的時(shí)候還有幾個(gè)坑要填。第一個(gè)坑是界面法向角的計(jì)算——代碼里那個(gè)angle變量沒有賦值實(shí)際計(jì)算時(shí)應(yīng)該遍歷所有固態(tài)鄰居取其相對(duì)當(dāng)前元胞的方位角做統(tǒng)計(jì)。更準(zhǔn)確的法向估算是用固態(tài)鄰居的質(zhì)量中心來推算二范數(shù)歸一化之后得到單位法向向量% 計(jì)算固相鄰居質(zhì)量中心方向 [cx_cm, cy_cm] solidNeighborCentroid(state, i, j, N); theta atan2(i - cx_cm, j - cy_cm);這么做比簡單亮度統(tǒng)計(jì)穩(wěn)定得多具體原因后面講各向異性畸變的時(shí)候再展開。第二個(gè)坑是溫度場的慢擴(kuò)散問題。真實(shí)的潛熱釋放和熱擴(kuò)散是耦合的不能簡單地把剛凝固元胞的溫度“原地”加上去因?yàn)闊崃啃枰車鷶U(kuò)散。正確的做法是在每個(gè)時(shí)間步中先算凝固潛熱源項(xiàng)再用顯式擴(kuò)散格式更新溫度場% 潛熱釋放 T T L_over_cp * new_solid; % 在凝固元胞上加上潛熱 % 溫度擴(kuò)散顯式格式 T_new T; for i 2:N-1 for j 2:N-1 T_new(i,j) T(i,j) alpha*dt/dx^2 * (T(i1,j)T(i-1,j)T(i,j1)T(i,j-1)-4*T(i,j)); end end T T_new;這種顯式格式有個(gè)穩(wěn)定性條件( \alpha \Delta t / \Delta x^2 \le 0.25 )。在這個(gè)約束下如果時(shí)間步長取得太大溫度場會(huì)振蕩發(fā)散。這也是為什么項(xiàng)目中對(duì)不同的材料參數(shù)需要重新校驗(yàn)一遍穩(wěn)定性條件。3.4 可視化實(shí)現(xiàn)MATLAB做CA可視化的最簡單方式是pcolor或imagesc。我用的是imagesc加自定義Colormapfigure(Position, [100, 100, 600, 500]); cmap [1 1 1; 0.9 0.9 0.9; 0.3 0.5 0.8]; % 白-淺灰-藍(lán) colormap(cmap); for step 1:nSteps % 更新狀態(tài)... if mod(step, 20) 1 imagesc(state); axis equal; axis tight; title(sprintf(Time step: %d, step)); drawnow; end endcolormap的三行顏色分別對(duì)應(yīng)液態(tài)、界面態(tài)和固態(tài)。在調(diào)試過程中我習(xí)慣把界面態(tài)用亮黃色突出顯示這樣能非常清楚地看到生長前沿的推進(jìn)情況比直接看固態(tài)區(qū)域要直觀得多。另一個(gè)很實(shí)用的可視化工具是保存每一幀為圖片格式然后合成為動(dòng)圖?？梢钥纯醋罱K形貌隨時(shí)間的變化趨勢if mod(step, 50) 1 frame getframe(gcf); writeVideo(videoObj, frame); end合出來的視頻對(duì)匯報(bào)和論文申請(qǐng)展示特別有用。3.5 性能優(yōu)化思路基礎(chǔ)代碼能跑通之后接下來要考慮性能。純循環(huán)版本在300x300網(wǎng)格下跑幾千步時(shí)間步每次都要遍歷所有界面元胞循環(huán)開銷非?？捎^。優(yōu)化方向有兩個(gè)第一個(gè)方向是對(duì)狀態(tài)更新做向量化處理。把界面元胞的坐標(biāo)和狀態(tài)信息抽到一維數(shù)組中對(duì)整批界面元胞同時(shí)計(jì)算速度增量而不是逐個(gè)遍歷。對(duì)于界面法向的計(jì)算可以預(yù)先用conv2卷積核計(jì)算固相分?jǐn)?shù)梯度然后從梯度方向一步得到法向角solidMask (state 2); gx conv2(double(solidMask), [-1 0 1; -2 0 2; -1 0 1], same); gy conv2(double(solidMask), [-1 -2 -1; 0 0 0; 1 2 1], same); theta atan2(-gy, -gx);這個(gè)技巧非常管用。用Sobel算子計(jì)算固相分布梯度得到的法向場更連續(xù)、更穩(wěn)定而且完全不用寫循環(huán)。速度提升至少一個(gè)數(shù)量級(jí)。第二個(gè)方向是只對(duì)界面元胞操作。用MATLAB的find函數(shù)索引所有界面元胞避免遍歷整個(gè)N×N矩陣中的所有非界面元胞。如果界面元胞數(shù)量只有總網(wǎng)格數(shù)的百分之幾這個(gè)優(yōu)化能顯著減少無效計(jì)算。4. 典型結(jié)果分析與物理形貌判讀4.1 枝晶形貌與端部過冷度用上面的模型跑通之后能直觀看到四重對(duì)稱的枝晶形態(tài)從中心晶核逐漸向外擴(kuò)展主枝晶臂沿預(yù)設(shè)的擇優(yōu)方向(theta_0 0^\circ) 時(shí)沿x和y方向延伸二次臂從主臂側(cè)向長出。這個(gè)形態(tài)與實(shí)驗(yàn)觀察到的金屬枝晶高度相似驗(yàn)證了模型的有效性。有個(gè)重要的物理解釋是枝晶尖端附近的過冷度比遠(yuǎn)離尖端的區(qū)域更高因?yàn)闈摕後尫派偎约舛艘暂^快速度推進(jìn)而枝晶臂之間的凹槽處溶質(zhì)和熱量積聚嚴(yán)重過冷度低生長緩慢。這個(gè)“尖端優(yōu)勢 凹槽抑制”的機(jī)制正是枝晶形貌得以保持的原因。如果把不同時(shí)刻的固相輪廓疊加畫在一起可以看到等間隔時(shí)間內(nèi)界面推進(jìn)的距離越來越小。這是因?yàn)殡S著枝晶生長釋放的潛熱在熔體中積累整體過冷度不斷降低。這個(gè)趨勢符合金屬凝固過程的物理規(guī)律——如果熔體體積有限溫度最終會(huì)回升到接近熔點(diǎn)凝固停止。4.2 各向異性強(qiáng)度與形態(tài)轉(zhuǎn)變各向異性強(qiáng)度系數(shù) (\varepsilon) 是控制形貌最重要的參數(shù)。我做了幾組對(duì)比實(shí)驗(yàn)結(jié)果差異很明顯(\varepsilon 0.02)形貌接近圓形四重對(duì)稱性很弱幾乎沒有明顯枝晶臂(\varepsilon 0.10)四個(gè)主臂清晰可辨二次臂開始出現(xiàn)(\varepsilon 0.20)主臂細(xì)長、二次臂發(fā)達(dá)出現(xiàn)明顯的枝晶側(cè)向分支(\varepsilon 0.30)出現(xiàn)尖端分裂和非物理的碎晶結(jié)構(gòu)建議把 (\varepsilon) 控制在0.1到0.2之間。如果二次臂結(jié)構(gòu)不明顯可以適當(dāng)增大如果出現(xiàn)異常分裂就要回調(diào)。4.3 與相場法結(jié)果的定性對(duì)比很多人會(huì)問CA的結(jié)果和相場法比到底差在哪我用一個(gè)表格來總結(jié)兩類方法在枝晶模擬中的典型差異對(duì)比維度元胞自動(dòng)機(jī)CA相場法Phase Field界面描述離散狀態(tài)界面寬度等于元胞尺寸連續(xù)擴(kuò)散界面界面寬度可調(diào)計(jì)算效率高適合大尺寸模擬低需要求解多組偏微分方程各向異性精度依賴法向估算精度有限直接在方程中控制精度高物理完備性需要額外耦合溫度/溶質(zhì)擴(kuò)散自洽耦合熱力學(xué)驅(qū)動(dòng)實(shí)現(xiàn)難度低幾百行代碼可搞定高需要較好的數(shù)值計(jì)算基礎(chǔ)適用場景形貌趨勢、工程級(jí)模擬精確物理研究、定量預(yù)測這個(gè)對(duì)比說明了CA模型的價(jià)值定位當(dāng)你不追求納米級(jí)別的定量精度但需要快速得到大尺度范圍內(nèi)的形貌趨勢時(shí)CA幾乎是效率最高的選擇。這也是CA在實(shí)際鑄造工藝模擬軟件里依然占有重要位置的原因。4.4 網(wǎng)格尺度敏感性CA方法有一個(gè)軟肋結(jié)果受網(wǎng)格尺度影響顯著。網(wǎng)格取得太粗枝晶臂顯得粗壯、碎網(wǎng)格取得太細(xì)計(jì)算量又上去了。我的調(diào)試經(jīng)驗(yàn)是至少要保證枝晶尖端半徑覆蓋5到8個(gè)元胞這樣計(jì)算出的形貌才不會(huì)明顯受網(wǎng)格幾何的“釘扎”影響。在Microsoft Excel里做個(gè)網(wǎng)格收斂性檢驗(yàn)盡管現(xiàn)在用MATLAB做模擬分別用100、200、400的網(wǎng)格跑相同物理參數(shù)對(duì)比尖端位置隨時(shí)間的曲線。如果三者結(jié)果偏差在5%以內(nèi)認(rèn)為網(wǎng)格已收斂如果偏差大需要加密網(wǎng)格。這個(gè)檢驗(yàn)步驟在正式研究里很重要發(fā)論文做模擬時(shí)必須要有。5. 常見問題、避坑指南與調(diào)試技巧5.1 枝晶沿對(duì)角線“長得過長”怎么辦最常見的異?，F(xiàn)象是枝晶臂沿45度對(duì)角線方向長得特別快形成X形而不是十字形。原因是Moore鄰域中斜對(duì)角鄰居的中心距離是 ( \sqrt{2} \Delta x)如果直接按距離計(jì)算捕獲概率對(duì)角線方向的推進(jìn)速度天然更快。解決辦法是修正距離效應(yīng)在計(jì)算捕獲概率或生長增量時(shí)對(duì)斜對(duì)角方向的鄰居乘一個(gè) (1/\sqrt{2}) 的權(quán)重因子。我在代碼中直接法向估算里用Sobel算子這個(gè)修正已經(jīng)包含在梯度計(jì)算里效果比手動(dòng)加權(quán)更自然。5.2 界面法向估算噪聲大如果用統(tǒng)計(jì)固相鄰居數(shù)量的方式估算法向界面稍微凹凸不平就會(huì)導(dǎo)致法向角劇烈抖動(dòng)進(jìn)而讓各向異性因子 (f_{aniso}) 波動(dòng)產(chǎn)生不規(guī)則的形貌。更穩(wěn)妥的方法是使用前面提到的Sobel卷積核計(jì)算固相分?jǐn)?shù)梯度然后用梯度方向作為界面法向。梯度場的連續(xù)性更好法向角不會(huì)跳變。這個(gè)方法是我調(diào)試多輪后得出的最佳實(shí)踐——最早用簡單統(tǒng)計(jì)時(shí)長出來的枝晶臂邊緣毛刺特別多換成Sobel之后界面光滑了一整個(gè)量級(jí)。5.3 溫度場發(fā)散溫度場用顯式格式擴(kuò)散時(shí)如果 (\alpha \Delta t / \Delta x^2 0.25)就可能出現(xiàn)數(shù)值振蕩甚至發(fā)散。這時(shí)的特征是枝晶周圍出現(xiàn)一圈一圈等間距的溫度異常帶。解決辦法很直接減小時(shí)間步長或減小熱擴(kuò)散系數(shù)。在保證 (\Delta t) 滿足條件的前提下盡量加大步長以減少循環(huán)次數(shù)。調(diào)試時(shí)可以先跑一個(gè)固定步數(shù)的測試觀察溫度最大值是否隨時(shí)間單調(diào)變化如果不是說明穩(wěn)定條件被破壞。5.4 模擬結(jié)果與理論尖端速度對(duì)比要驗(yàn)證CA模型是否靠譜一個(gè)經(jīng)典的定量驗(yàn)證方法是比較枝晶尖端速度的模擬值與KGT理論預(yù)測值。做法是每次記錄尖端位置的推進(jìn)距離除以時(shí)間步長得到尖端速度然后與理論公式計(jì)算值對(duì)比。如果在相同過冷度下模擬值和理論值的偏差在10%以內(nèi)模型基本可靠。如果偏差過大優(yōu)先檢查各向異性強(qiáng)度取值是否合理以及潛熱反饋項(xiàng)是否設(shè)置正確。5.5 偶發(fā)的不對(duì)稱生長有時(shí)候模擬出來的枝晶左右不對(duì)稱一側(cè)臂比另一側(cè)長。這個(gè)問題的來源通常是初始化時(shí)晶核不是完美的圓形或者邊界條件沒有對(duì)稱設(shè)置。解決方法是使用對(duì)稱初始化——把晶核幾何和后續(xù)的邊界處理都設(shè)計(jì)成關(guān)于中心點(diǎn)對(duì)稱的矩陣操作并且在邊界處采用對(duì)稱邊界條件反射邊界避免數(shù)值單側(cè)影響傳播。6. 項(xiàng)目經(jīng)驗(yàn)總結(jié)與擴(kuò)展思路自己做這個(gè)項(xiàng)目走下來最大的體會(huì)是CA模型的代碼實(shí)現(xiàn)本身不難真正的功夫在物理機(jī)制映射和參數(shù)調(diào)試。每次調(diào)整物理參數(shù)就像在做一個(gè)虛擬冶金實(shí)驗(yàn)需要仔細(xì)觀察枝晶形貌的變化趨勢才能判斷模型是否真正抓住了關(guān)鍵動(dòng)力學(xué)因素。如果想把項(xiàng)目往更深的方向推進(jìn)有幾個(gè)可行的擴(kuò)展方向第一在模型中引入溶質(zhì)場模擬合金凝固過程中的成分偏析。這需要在每個(gè)元胞上額外存儲(chǔ)一個(gè)濃度值并在界面元胞凝固時(shí)釋放溶質(zhì)然后用擴(kuò)散方程更新溶質(zhì)場。加入溶質(zhì)場之后枝晶臂之間的微觀偏析形態(tài)會(huì)和實(shí)驗(yàn)吻合得更好。第二把二維模型擴(kuò)展成三維。三維CA的代碼思路一樣但鄰域從8個(gè)鄰居變成26個(gè)計(jì)算量大到可能需要并行化處理。MATLAB的分布式計(jì)算工具箱或直接在GPU上跑conv2卷積可以大幅度加速。第三耦合熱力學(xué)模塊。把真實(shí)合金相圖的熱力學(xué)數(shù)據(jù)庫如CALPHAD嵌入到CA模型中讓界面速度直接由局部平衡溫度和溶質(zhì)成分計(jì)算不再依賴簡化線性關(guān)系。這樣做出的模擬逐漸接近工業(yè)合金實(shí)際凝固工藝。我在實(shí)際調(diào)通三維版本之后又把優(yōu)化從MATLAB搬到Python重新實(shí)現(xiàn)了一遍發(fā)現(xiàn)只要掌握了模型邏輯換語言只是兩三天的事。所以強(qiáng)烈建議在這個(gè)項(xiàng)目上先把CA的邏輯吃透這比記住任何具體的代碼寫法都重要。最后想說的是這類“簡單規(guī)則涌現(xiàn)復(fù)雜形貌”的模擬確實(shí)有種獨(dú)特的吸引力每跑出一張清晰漂亮的枝晶圖都像在看一個(gè)微觀世界的雕塑過程。希望這篇分享能讓你在MATLAB里跑出自己的第一朵枝晶花。