做相變仿真這行久了你會(huì)發(fā)現(xiàn)最繞不開的數(shù)學(xué)模型就是焓法。不管你是做金屬凝固、石蠟儲(chǔ)能、還是電池?zé)峁芾砝锏娜廴谙嘧儽砻嫔峡磫栴}是“材料隨著溫度變了物理狀態(tài)”但真正落到方程里難點(diǎn)全部集中在相變界面的處理上——界面在移動(dòng)、潛熱在釋放、流場跟著溫度走如果直接用傳統(tǒng)的比熱容法去硬解數(shù)值上有一萬種方式讓你懷疑人生。這篇文章我就把自己在多個(gè)項(xiàng)目里用焓法耦合溫度場和流場做相變模型的核心經(jīng)驗(yàn)拆開講從模型思路、數(shù)學(xué)表達(dá)、實(shí)操設(shè)置到問題排查盡量把每一步為什么這么做、有哪些坑、參數(shù)怎么選說透。適合正在被相變收斂問題折磨的仿真工程師也適合剛?cè)腴T想做材料相變過程模擬的同學(xué)參考。1. 相變模型為什么繞不開焓法先搞清楚問題本質(zhì)1.1 相變問題的難點(diǎn)在于移動(dòng)邊界先說個(gè)最簡單的物理圖景。一塊材料放在熱源旁邊溫度升高到熔點(diǎn)固體開始變成液體。這個(gè)過程里固液界面是不斷移動(dòng)的——外圈先熔化熔化的液體如果受熱不均還會(huì)流動(dòng)流動(dòng)又反過來參與傳熱界面就變成不規(guī)則的、隨時(shí)間變化的曲面。你讓仿真軟件去追蹤這個(gè)界面本質(zhì)上是在解一個(gè)移動(dòng)邊界問題。移動(dòng)邊界問題最讓人頭疼的地方在于界面位置不是事先知道的結(jié)果而是求解過程中由能量平衡自動(dòng)決定的。如果你用界面追蹤法比如水平集法在相變里的一些應(yīng)用每一步都要顯式地重構(gòu)界面、設(shè)置邊界條件、更新網(wǎng)格復(fù)雜幾何下很容易把網(wǎng)格扭成一團(tuán)而且兩相密度差明顯時(shí)界面附近的質(zhì)量和能量守恒很難同時(shí)保證。焓法不太一樣。它不直接追蹤界面而是把一個(gè)稱為“總焓”的物理量作為基本未知量。固相、液相、兩相區(qū)統(tǒng)一用一套控制方程描述相變潛熱被吸收進(jìn)焓值的變化里界面位置在求解之后從焓值反推出來不需要單獨(dú)處理。1.2 焓法和等效熱容法的本質(zhì)區(qū)別很多初學(xué)者一開始接觸的是等效熱容法也叫表觀熱容法。思路是把相變潛熱折算到相變溫度區(qū)間內(nèi)的比熱容上比如材料在固相時(shí)比熱是 2000 J/(kg·K)在液相時(shí)是 2100 J/(kg·K)相變潛熱是 200 kJ/kg相變區(qū)間是 40℃ 到 60℃那就把這段 20℃ 區(qū)間內(nèi)的比熱加成(200000 / 20) 10000 J/(kg·K)等效熱容一下子翻了五倍。這個(gè)方法好用是好用但它有個(gè)致命弱點(diǎn)等效熱容在相變區(qū)間形成了很高的尖峰如果網(wǎng)格不夠密、時(shí)間步長不夠小溫度場在穿越相變區(qū)間時(shí)會(huì)發(fā)生振蕩甚至出現(xiàn)溫度來回抖動(dòng)但不收斂的情況。而且尖峰區(qū)間內(nèi)的積分精度對(duì)網(wǎng)格極其敏感你在 A 方案里算得挺好把網(wǎng)格一加密等效熱容的積分反而出偏差了。焓法不這樣做。它的控制方程里直接把焓 h 作為待求變量潛熱以焓的連續(xù)變化形式存在不存在高聳的尖峰。從數(shù)值穩(wěn)定性上講焓法的非線性主要體現(xiàn)在焓到溫度的反解上處理起來要溫和得多。我個(gè)人的建議很直接只要你的模型里有明顯的相變現(xiàn)象且相變區(qū)間相對(duì)較小直接用焓法。等效熱容法只適合粗略估算穩(wěn)態(tài)溫度分布的情況動(dòng)態(tài)過程和時(shí)間相關(guān)的相變比例焓法的結(jié)果可靠得多。2. 焓法數(shù)學(xué)模型搭建幾個(gè)關(guān)鍵表達(dá)式必須吃透2.1 液相分?jǐn)?shù)的定義與插值方式焓法的核心變量是總焓 h它的物理定義是顯熱和潛熱之和。對(duì)于純物質(zhì)或者窄相變區(qū)間材料通常把焓寫成h ∫ c_p dT f_L · L其中 c_p 是比熱容L 是單位質(zhì)量潛熱f_L 是液相分?jǐn)?shù)。f_L 的取值從 0全固相到 1全液相連續(xù)變化。這個(gè) f_L 是整個(gè)模型里的靈魂變量溫度場、流場、材料參數(shù)都要掛在它身上。對(duì)于純物質(zhì)理想情況相變發(fā)生在單一溫度點(diǎn) T_mf_L 是一個(gè)階躍函數(shù)。但數(shù)值計(jì)算里沒法處理無窮陡的階躍必須給它一個(gè)小的過渡區(qū)間。風(fēng)扇常見做法是定義固相線溫度 T_s 和液相線溫度 T_l然后在線性插值假設(shè)下f_L (T - T_s) / (T_l - T_s)當(dāng) T_s ≤ T ≤ T_lf_L 0當(dāng) T T_sf_L 1當(dāng) T T_l這是工程中最常見的線性插值。如果相變不是線性的也可以用 smoother 的 S 形函數(shù)比如基于誤差函數(shù)或者正弦的函數(shù)補(bǔ)插但對(duì)于大多數(shù)工程模擬線性插值已經(jīng)能給出足夠準(zhǔn)確的結(jié)果而且數(shù)值穩(wěn)定性更好。這里有個(gè)細(xì)節(jié)值得注意相變區(qū)間的寬度不是隨便設(shè)的。區(qū)間設(shè)置太窄比如 T_l - T_s 只有 1℃焓法模型里熱容變化劇烈求解器每一步都要處理強(qiáng)非線性極易不收斂。區(qū)間設(shè)置太寬相變過程被人為拉長潛熱釋放被抹平均攤結(jié)果溫度曲線失真。通常建議把區(qū)間寬度設(shè)置為材料實(shí)際相變范圍的 1.5 到 2 倍或者至少大于二維/三維網(wǎng)格尺寸對(duì)應(yīng)的最小特征時(shí)間內(nèi)的溫度增量。2.2 焓值反解溫度的邏輯焓法里第二個(gè)關(guān)鍵點(diǎn)是如何從焓值反推出溫度。因?yàn)榭刂品匠探獬鰜淼氖庆?h但材料參數(shù)、換熱系數(shù)、浮力項(xiàng)都需要溫度 T所以必須建立 h 到 T 的逆關(guān)系。在分段線性假設(shè)下反解邏輯是這樣的如果 h h_s ∫_0^T_s c_p,s dT說明處于全固相T h / c_p,s如果 h h_l h_s L ∫_T_s^T_l c_p,l dT說明處于全液相T T_l (h - h_l) / c_p,l如果 h_s ≤ h ≤ h_l說明處于兩相區(qū)T T_s (h - h_s) / (h_l - h_s) · (T_l - T_s)這個(gè)分段邏輯在軟件里用條件表達(dá)式寫直接用 h 作為自變量。實(shí)現(xiàn)時(shí)的注意點(diǎn)當(dāng)材料在相變區(qū)間內(nèi)時(shí)它的實(shí)際比熱不是恒定的而是等效變化的。建議整個(gè)相變區(qū)間的焓差顯熱加潛熱用一個(gè)統(tǒng)一表達(dá)式計(jì)算避免分段點(diǎn)處的不連續(xù)導(dǎo)致殘差跳變。2.3 參數(shù)實(shí)例給一個(gè)可以抄作業(yè)的案例我在一個(gè)儲(chǔ)能項(xiàng)目里做過典型的石蠟相變模型參數(shù)直接拿來舉例點(diǎn)地?cái)?shù)據(jù)做了一定歸一化工程上符合典型石蠟特性參數(shù)數(shù)值說明固相密度860 kg/m3固液相密度差不大時(shí)可近似常數(shù)固相比熱2100 J/(kg·K)液相密度760 kg/m3如果考慮自然對(duì)流浮力這個(gè)差要體現(xiàn)液相比熱2400 J/(kg·K)導(dǎo)熱系數(shù)0.25 W/(m·K)固液導(dǎo)熱接近暫不分段潛熱180 kJ/kg典型石蠟級(jí)入能量密度固相線 T_s45℃液相線 T_l55℃實(shí)際石蠟熔化區(qū)間比這更寬我這里做的簡化熱膨脹系數(shù)0.001 /K用于自然對(duì)流浮力項(xiàng)動(dòng)力粘度0.005 Pa·s液態(tài)石蠟典型粘度隨溫度變化略照這個(gè)參數(shù)你算一下相變區(qū)間的焓差h_s 2100 × 45 94500 J/kgh_l 94500 180000 2400 × (55 - 45) 94500 180000 24000 298500 J/kg也就是在 45℃ 到 55℃ 之間每公斤石蠟要吸收大約 298.5 kJ 的能量其中潛熱占了 180 kJ比例非常高。這個(gè)數(shù)字直接決定了儲(chǔ)能時(shí)間和加熱速率你在后處理時(shí)看到某個(gè)區(qū)域內(nèi)溫度始終停在 45-55℃ 區(qū)間就是在相變儲(chǔ)能溫度曲線出現(xiàn)一個(gè)明顯的平臺(tái)段。有了這個(gè)參數(shù)表接下來搭建多物理場耦合就有了抓手。3. 溫度場和流場的耦合方式從單向到雙向的升級(jí)3.1 能量方程里的對(duì)流項(xiàng)別漏掉單純導(dǎo)熱模型只解能量方程溫度場和流場互不影響。但實(shí)際相變過程里液相一旦出現(xiàn)自然對(duì)流就會(huì)顯著改變熱量輸運(yùn)。這時(shí)候能量方程需要寫成包含對(duì)流項(xiàng)的形式ρ c_p ?T/?t ρ c_p u · ?T ?·(k ?T) Q_latent其中 Q_latent 是潛熱源項(xiàng)在焓法框架里可以做到溫度方程里。也可以用焓形式直接寫ρ ?h/?t ρ c_p u · ?T ?·(k ?T)注意對(duì)流項(xiàng)里用的是 c_p 乘溫度梯度不是焓梯度這在數(shù)學(xué)上更符合能量守恒的物理分解也是焓法實(shí)現(xiàn)時(shí)最容易出錯(cuò)的表達(dá)式。我在實(shí)操中的一個(gè)體會(huì)是如果材料固相和液相的比熱差別較大對(duì)流項(xiàng)里的 c_p 到底取固相值還是液相值會(huì)導(dǎo)致明顯的溫差。通常做法是使用與溫度相關(guān)的分段 c_p(T)借助液相分?jǐn)?shù)對(duì)固液兩相的值做插值c_p,eff c_p,s (1 - f_L) c_p,l f_L這個(gè)表達(dá)式看著簡單但它保證了能量守恒在固液界面處的一致性。直接把 c_p 設(shè)成常數(shù)會(huì)高估或低估融化速度跑出來的相變時(shí)間差個(gè) 10%-20% 很正常。3.2 自然對(duì)流的浮力驅(qū)動(dòng)Boussinesq 近似液態(tài)相變材料里的流動(dòng)絕大多數(shù)是自然對(duì)流由溫度引起的密度差驅(qū)動(dòng)。如果你直接寫流體密度隨溫度變化的全可壓縮 N-S 方程求解器會(huì)處理大量和相變無關(guān)的聲學(xué)尺度物理又慢又不穩(wěn)定。工程上最常用的處理是 Boussinesq 近似。思路是在動(dòng)量方程里只保留浮力項(xiàng)的密度變化把動(dòng)量方程寫成ρ0 ?u/?t ρ0 (u · ?)u -?p μ ?2u ρ0 β (T - T_ref) g其中 ρ0 是參考密度β 是熱膨脹系數(shù)T_ref 通常取初始溫度或平均溫度g 是重力矢量方向。浮力項(xiàng)只在液相區(qū)有物理意義固相區(qū)因?yàn)樗俣缺蛔枘醿鼋Y(jié)該項(xiàng)自然不起作用。這里有一個(gè)先決條件密度差不能太大一般要求 ρ_l / ρ_s 在 0.9-1.1 之間。石蠟這種材料固液相密度差在 10%-15%Boussinesq 近似還能湊合用如果要更嚴(yán)格就得走變密度流模型。金屬相變時(shí)固液相密度差也不大同樣適用。3.3 糊狀區(qū)阻尼保證固相區(qū)速度為 0 的物理手段這是相變耦合模型里最容易忽略、最影響收斂的關(guān)鍵項(xiàng)。試想液相分?jǐn)?shù) f_L 0 的區(qū)域材料還是固體不應(yīng)該有任何流動(dòng)速度。但如果只靠前面那些方程速度場會(huì)把整個(gè)計(jì)算域都當(dāng)流體去算固相區(qū)會(huì)出現(xiàn)泄漏速度物理上完全錯(cuò)誤。解決思路是在動(dòng)量方程里加一個(gè)“糊狀區(qū)阻力項(xiàng)”模擬固相區(qū)對(duì)流動(dòng)的巨大阻力。最常用的是基于 Darcy 定律的 Carman-Kozeny 表達(dá)式F -C (1 - f_L)2 / (f_L3 ε) · u這里的 C 是糊狀區(qū)常數(shù)工程上取 10^4 到 10^7 甚至更大 ε 是個(gè)很小的正數(shù)防止分母為零通常取 10^-3 量級(jí)。當(dāng) f_L 0全固相時(shí)這個(gè)阻力項(xiàng)趨向無窮大物理上被 C 限制在一個(gè)超大值速度被強(qiáng)推為 0。當(dāng) f_L 1全液相時(shí)阻力項(xiàng)為 0不影響流動(dòng)。有個(gè)實(shí)用經(jīng)驗(yàn)C 值不能無條件取最大。C 取 10^7 時(shí)動(dòng)量方程里這個(gè)阻力項(xiàng)的數(shù)量級(jí)非常大會(huì)破壞矩陣的條件數(shù)導(dǎo)致非線性求解器每步都很難收斂。我的習(xí)慣是先從 10^4 開始試算確認(rèn)固相區(qū)的速度泄漏量級(jí)已經(jīng)低于 10^-6 m/s 之后再逐步加加到 10^6 為止。大部分情況下 10^5 到 10^6 足夠。這個(gè)阻尼項(xiàng)還順便解決了相變界面的數(shù)值穩(wěn)定性問題。沒有它液相分?jǐn)?shù)在 0-1 之間的界面區(qū)域會(huì)產(chǎn)生虛假的對(duì)流嚴(yán)重時(shí)會(huì)在界面兩側(cè)形成振蕩的溫度漣漪畫面非常難看。3.4 耦合方式是單向還是雙向什么時(shí)候必須雙向如果你的相變材料是純導(dǎo)熱主導(dǎo)比如熱源功率很小、液相區(qū)流動(dòng)很弱那可以只做單向耦合先解溫度場根據(jù)溫度分布算出液相分?jǐn)?shù)再把液相分?jǐn)?shù)映射到材料參數(shù)粘度、熱容、導(dǎo)熱系數(shù)更新流場速度。這種解法和解耦求解速度快很多。但如果自然對(duì)流的 Rayleigh 數(shù)高熱源功率大流場反過來顯著影響溫度分布尤其你想準(zhǔn)確預(yù)測“融化前鋒的形狀”時(shí)單向耦合完全不夠。必須做雙向耦合也就是每個(gè)時(shí)間步內(nèi)同步求解動(dòng)量方程、連續(xù)性方程、能量方程和焓反解所有變量聯(lián)動(dòng)更新。焓法的非線性本身已經(jīng)夠強(qiáng)雙向耦合會(huì)讓非線性進(jìn)一步增強(qiáng)這時(shí)候求解策略就成了關(guān)鍵。我個(gè)人的取舍標(biāo)準(zhǔn)是看無量綱數(shù) Ra 是否超過臨界值Ra g β (T_hot - T_cold) H3 / (ν α)H 是特征尺度ν 是運(yùn)動(dòng)粘度α 是熱擴(kuò)散率。Ra 超過 10^5 時(shí)自然對(duì)流已經(jīng)是主要傳熱機(jī)制必須雙向耦合。低于 10^4 時(shí)單向耦合誤差可以接受。這個(gè)判斷標(biāo)準(zhǔn)我用了很多年基本靠譜。4. 實(shí)操流程從幾何搭建到后處理的完整執(zhí)行方案4.1 物理場選擇與接口組合我不強(qiáng)調(diào)具體軟件品牌工程圈都知道這類多物理場工具的通用流程但無論你用哪套平臺(tái)典型的相變-流-熱耦合模型物理場組合都是流體流動(dòng)接口求解 N-S 方程 連續(xù)性方程固體傳熱/流體傳熱接口用于能量方程自定義變量/方程用于計(jì)算液相分?jǐn)?shù) f_L、焓 h、溫度反解 T(h)這里有個(gè)自由度的選擇問題。用哪個(gè)溫度變量去求解決定了表達(dá)式寫法。我用過兩種方案。方案 A主變量直接用溫度 T焓法以“表觀熱容”加“潛熱源項(xiàng)”的形式實(shí)現(xiàn)。具體做法是能量方程里把熱容替換成熱焓對(duì)時(shí)間的偏導(dǎo)數(shù)項(xiàng)同時(shí)在瞬間添加一個(gè)潛熱釋放源項(xiàng)。優(yōu)點(diǎn)是與傳熱接口兼容度最高后處理直接看溫度。缺點(diǎn)是實(shí)質(zhì)上回到等效熱容路線對(duì)網(wǎng)格敏感。方案 B主變量改成焓 h 或使用熱容形式材料但這種寫法通常要到更高階的自定義弱形式才能實(shí)現(xiàn)通用接口下反而難操作。對(duì)于絕大多數(shù)工程場景我建議你使用表觀熱容形式的寫法但要把熱容表達(dá)式精心構(gòu)造避免直接用尖峰函數(shù)。具體做法定義等效熱容 C_eff(T)在整個(gè)溫區(qū)內(nèi)連續(xù)C_eff c_p,s (1 - f_L) c_p,l f_L L / (T_l - T_s) · df_L/dT由于 f_L 是線性插值時(shí)df_L/dT 在兩個(gè)相變邊界點(diǎn)是階躍的。為了讓數(shù)值穩(wěn)定我給 df_L/dT 用一個(gè)小平滑實(shí)際效果比理想階躍好得多。4.2 幾何建模與網(wǎng)格分區(qū)建立二維或軸對(duì)稱模型時(shí)我通常把計(jì)算域畫成矩形或圓柱腔體。關(guān)鍵網(wǎng)格處理原則相變前沿可能經(jīng)過的區(qū)域加密加密區(qū)的網(wǎng)格尺寸不超過總特征長度的 1/30壁面處添加邊界層網(wǎng)格至少 5 層第一層厚度根據(jù)預(yù)期對(duì)流強(qiáng)度設(shè)定確保 y 滿足低 Re 渦模擬要求如果使用高度函數(shù)初始化液相分?jǐn)?shù) 0初始時(shí)刻全固相逐步融化焓法和 Level Set 不同不需要界面附近高精度貼體網(wǎng)格因?yàn)榻缑媸请[式捕捉的只需要確保相變區(qū)間溫度變化對(duì)應(yīng)的焓值能被網(wǎng)格積分覆蓋。一個(gè)可量化的經(jīng)驗(yàn)在最大局部熱流密度下網(wǎng)格間距需要滿足使得一個(gè)時(shí)間步內(nèi)溫度變化小于 0.5×(T_l - T_s)。4.3 邊界條件設(shè)置實(shí)例以一個(gè)帶內(nèi)熱源的方形儲(chǔ)能單元為例左側(cè)壁面加熱溫度恒定 80℃右側(cè)壁面自然對(duì)流散熱對(duì)流系數(shù) 10 W/(m2·K)環(huán)境溫度 25℃上下壁面絕熱初始溫度 25℃全固相速度場初始為 0重力方向沿-y軸浮力驅(qū)動(dòng)液相向上流動(dòng)這些條件寫起來都不復(fù)雜但耦合設(shè)置里有兩個(gè)坑第一個(gè)坑是壓力參考點(diǎn)必須設(shè)置。純自然對(duì)流問題速度場是由浮力驅(qū)動(dòng)的壓力沒有唯一絕對(duì)參考不設(shè)參考點(diǎn)的話壓力方程奇異求解器會(huì)報(bào)錯(cuò)或長時(shí)間不收斂。我在所有涉及自然對(duì)流的模型里計(jì)算域角落固定一個(gè)壓力零點(diǎn)。第二個(gè)坑是壁面附近的浮力方向。如果加熱壁在左側(cè)液體受熱后沿左壁上升在頂部轉(zhuǎn)向右在右壁冷卻下降形成一個(gè)逆時(shí)針的大渦環(huán)。你不能把所有壁面都設(shè)為無滑移就完事還要檢查浮力項(xiàng)的方向設(shè)定是否正確。我第一次做的時(shí)候把重力項(xiàng)寫錯(cuò)方向算出來的渦環(huán)方向完全反了溫度場倒是沒太大差異但流場后處理一眼就看出來有問題。4.4 瞬態(tài)求解器參數(shù)與時(shí)間步策略瞬態(tài)求解設(shè)置里最影響結(jié)果是時(shí)間步進(jìn)方式用 BDF 2二階向后差分或等價(jià)的隱式時(shí)間格式避免顯示格式的時(shí)間步限制初始時(shí)間步設(shè)得非常小取 0.001 秒量級(jí)讓求解器從全固相平穩(wěn)啟動(dòng)開啟自適應(yīng)時(shí)間步進(jìn)步長上限控制在 1 秒左右下限 1e-5 秒非線性殘差容差設(shè)置為默認(rèn)值的 1/10比如 0.001保證界面推移的精度焓法模型的非線性非常強(qiáng)尤其是溫度穿越相變區(qū)間的瞬間物理特性突變。這種情況下直接全體方程全耦合求解求解器每步的雅可比矩陣更新成本很高但值得。我常用的策略是“先導(dǎo)熱后對(duì)流”兩步走第一步關(guān)閉流體求解只算純導(dǎo)熱的相變過程 0.5 秒得到初始溫度場第二步再開啟流體耦合求解。這個(gè)方法能快速跨過最困難的啟動(dòng)段節(jié)省大量時(shí)間并且顯著提高穩(wěn)定性。4.5 后處理重點(diǎn)指標(biāo)與可視化模型算完后處理要看哪些東西第一個(gè)是液相分?jǐn)?shù)云圖。它是最直觀的相變過程體現(xiàn)你能夠在不同時(shí)間切片上看到固相逐漸縮小的全過程輪廓。可以用 0.5 等值線代表界面位置看它的推進(jìn)路徑是否合理。第二個(gè)是溫度歷史曲線。選取幾個(gè)監(jiān)測點(diǎn)分別位于靠近熱源、中間區(qū)域、遠(yuǎn)離熱源位置畫出溫度隨時(shí)間變化。典型的相變曲線應(yīng)該出現(xiàn)斜率先升、中間變緩相變平臺(tái)、隨后斜率恢復(fù)的過程。平臺(tái)段的長度對(duì)應(yīng)潛熱吸收時(shí)間。第三個(gè)是速度場?？匆合鄥^(qū)的渦環(huán)結(jié)構(gòu)是否形成渦心位置大概在哪個(gè)區(qū)域最大流速數(shù)量級(jí)是否符合自然對(duì)流的常識(shí)通常幾毫米每秒到幾十毫米每秒。第四個(gè)是能量守恒校驗(yàn)。這個(gè)方法非常必要把整個(gè)計(jì)算域的焓隨時(shí)間的變化率與邊界熱流積分總量去做對(duì)比。我每次都要算這個(gè)總量差的曲線如果能量不平衡超過 5%說明模型或者求解參數(shù)有問題。曾經(jīng)有一次我漏掉了相變區(qū)間的熱量計(jì)算能量曲線直接缺口排查了半天才發(fā)現(xiàn)是液相分?jǐn)?shù)插值表達(dá)式寫錯(cuò)了一個(gè)邊界。5. 常見問題與排查技巧實(shí)錄5.1 溫度場鋸齒振蕩界面像鋸齒怎么辦鋸齒振蕩基本都出在液相分?jǐn)?shù)和潛熱的處理上。排查順序第一步檢查相變區(qū)間的寬度是不是太小。如果 T_l - T_s 只有 0.5℃網(wǎng)格稍粗就被溫度跳躍跨過了必然振蕩。把區(qū)間放大到 2-5℃如果振蕩消失就是這個(gè)問題。第二步檢查等效熱容是不是存在尖峰。理想線性插值產(chǎn)生的 df_L/dT 是常數(shù)不會(huì)尖峰但如果你用了非線性光滑函數(shù)相位變化集中在極窄區(qū)間內(nèi)等效熱容就有尖峰隱患。第三步檢查網(wǎng)格在相變前沿區(qū)域是否足夠細(xì)。界面處溫度梯度大網(wǎng)格太粗時(shí)界面推進(jìn)是一格一格跳的看起來像鋸齒。5.2 計(jì)算發(fā)散提示非線性求解器不收斂最典型的場景發(fā)生在材料剛開始熔化那幾步。這時(shí)候潛熱突然釋放溫度對(duì)焓的非線性映射變化劇烈牛頓迭代容易振蕩。排查技巧第一調(diào)小初始時(shí)間步長。把這個(gè)參數(shù)縮小一個(gè)數(shù)量級(jí)通常就能跨過啟動(dòng)難關(guān)。第二給液相分?jǐn)?shù)加“松動(dòng)”。英語里叫 under-relaxation也就是液分更新不用全量每一步只更新一部分f_L,new f_L_old α (f_L_calc - f_L_old)α 取 0.5-0.8。這個(gè)做法犧牲一點(diǎn)點(diǎn)時(shí)間精度但換來了極其穩(wěn)定的收斂過程。第三如果時(shí)間步長已經(jīng)非常小而仍然發(fā)散查看是不是糊狀區(qū)常數(shù) C 太大。把 C 從 10^6 降到 10^4 試試。有的模型里 C 過強(qiáng)會(huì)讓動(dòng)量方程局部條件數(shù)劇增。第四檢查初始條件是否自洽。初始溫度恰好卡在相變區(qū)間邊界上有時(shí)候會(huì)觸發(fā)初始?xì)埐顦O大。把初始溫度稍微調(diào)整 0.1℃或者初始給一個(gè)極小的液相分?jǐn)?shù) 1e-4也能繞過。5.3 速度場在固相區(qū)出現(xiàn)非零速度如果固相區(qū)速度不是零先檢查糊狀區(qū)阻尼項(xiàng)是不是寫進(jìn)了方程。這個(gè)大概率是漏寫了這個(gè)源項(xiàng)。不是漏寫的話檢查 C 值是不是太小增大一個(gè)數(shù)量級(jí)看看。另外注意 f_L 的數(shù)值精度。如果 f_L 在固相區(qū)因?yàn)閿?shù)值誤差變成 1e-5 而不是精確的 0那個(gè)小值被代入阻尼公式1-f_L2/(f_L3ε)雖然分母有 ε 保護(hù)但 C 值如果不夠大仍舊會(huì)出現(xiàn)微小泄漏。5.4 相變時(shí)間過短或過長和理論值明顯不符相變持續(xù)總時(shí)長對(duì)不上理論估算通常在兩個(gè)地方第一潛熱值是不是被正確“激活”了。有時(shí)候材料定義里 L 被填成 0潛熱完全不起作用融化時(shí)間比理論值少 40% 以上。排查方式是在能量守恒曲線里看潛熱貢獻(xiàn)段。第二相變區(qū)間寬度被設(shè)置得過寬。假設(shè)實(shí)際材料相變區(qū)間只有 1℃但你在模型里為了數(shù)值穩(wěn)定把它拉成 20℃潛熱在這個(gè)過程中均勻釋放雖然消耗的總能量一致但熱流分配方式和實(shí)際情況不一樣會(huì)導(dǎo)致融化過程持續(xù)時(shí)間偏長。這里要做一個(gè)平衡數(shù)值穩(wěn)定和物理保真之間的折中。5.5 快速排查清單癥狀首要懷疑點(diǎn)推薦措施溫度鋸齒振蕩相變區(qū)間太窄/網(wǎng)格太粗增大 T_l-T_s加密網(wǎng)格不收斂初始步長過大初始步長縮到 1e-4 再試流場凍結(jié)/泄漏糊狀區(qū)阻尼缺失檢查 C(1-f_L)2/(f_L3ε)項(xiàng)相變時(shí)間偏短潛熱沒接上檢查 L 的表達(dá)式鏈接溫度平臺(tái)不明顯等效熱容無尖峰但區(qū)間寬適當(dāng)壓縮相變區(qū)間5.6 一個(gè)非常隱蔽的坑熱物性參數(shù)插值固相和液相的熱物性不同尤其導(dǎo)熱系數(shù)可能差好幾倍。如果你只在材料屬性里用分段函數(shù)切換而沒有通過液相分?jǐn)?shù)做光滑過渡界面處的導(dǎo)熱系數(shù)突變會(huì)帶來嚴(yán)重的數(shù)值不連續(xù)性。建議所有依賴于相態(tài)的物性參數(shù)都統(tǒng)一寫成k_eff k_s (1 - f_L) k_l f_L這樣的形式。這不是額外的物理假設(shè)而是數(shù)值處理上的必需操作。我當(dāng)年做金屬相變時(shí)固態(tài)的導(dǎo)熱系數(shù)比液態(tài)高大約 60%最初沒做過渡處理溫度場界面處出現(xiàn)了無物理意義的高頻擾動(dòng)。改成這個(gè)直接插值后一切都穩(wěn)了。關(guān)于參數(shù)選擇我最后再啰嗦幾句相變仿真做到最后瓶頸往往不在軟件操作而在你對(duì)材料相變行為的物理理解有多深。焓法這個(gè)數(shù)學(xué)框架本身非常優(yōu)雅它把一個(gè)強(qiáng)非線性移動(dòng)邊界問題轉(zhuǎn)化為單域固定網(wǎng)格問題只要掌控了液相分?jǐn)?shù)、焓反解、糊狀區(qū)阻尼這三個(gè)核心環(huán)節(jié)絕大多數(shù)工程相變場景都能穩(wěn)定算出來。最后分享一個(gè)我個(gè)人的實(shí)踐經(jīng)驗(yàn)?zāi)玫饺魏蜗嘧兡Ｐ偷谝患虏皇羌敝{(diào)求解器而是先驗(yàn)算焓差和能量總量。在 Excel 里花十分鐘把固相升至熔點(diǎn)、潛熱吸收、液相繼續(xù)升溫這三段能量相加和總加熱功率乘以預(yù)期時(shí)間做個(gè)對(duì)比。如果這個(gè)粗算和仿真結(jié)果的量級(jí)對(duì)不上后面調(diào)什么都白費(fèi)。把這一步養(yǎng)成習(xí)慣之后你做相變模型的速度和一次成功率高得驚人。