有段時間沒管過這類問題了看到標(biāo)題里那句結(jié)果總不理想求助大佬我是真能體會那種憋屈感。COMSOL做納秒脈沖激光燒蝕模型看著不難傳熱加移動網(wǎng)格點幾下就能跑可出來的溫度場就是不對勁不是中心溫度炸到幾百萬開爾文就是燒蝕坑邊緣鋸齒狀更糟的是網(wǎng)格跑著跑著直接畸形報錯。我當(dāng)年在這個坑里蹲了至少一個月后來把整個模型從物理到數(shù)值扒了一遍才理順。這篇東西我就按自己的排查順序來寫先拆物理過程再講移動網(wǎng)格然后是熱源建模最后給一套能直接抄作業(yè)的實操流程和排查清單。目標(biāo)是讓你跑出來的溫度場既對得上理論量級也能呈現(xiàn)合理的燒蝕形貌。1. 先別急著調(diào)參數(shù)把納秒激光燒蝕的物理過程拆清楚很多人在模型里折騰半天發(fā)現(xiàn)結(jié)果不對問題根源其實不在軟件操作而在物理圖景沒建對。納秒脈沖激光和材料的相互作用跟連續(xù)激光、飛秒激光完全是兩碼事。搞清這一點后面的每一項設(shè)置你才知道自己到底在模擬什么東西。1.1 納秒脈沖激光燒蝕到底發(fā)生了什么納秒脈沖的脈寬通常在1到100納秒范圍這個時間尺度決定了激光能量在材料內(nèi)部的沉積和擴散方式。拿金屬舉例熱擴散深度可以用一個簡單式子估算l_th sqrt(4 * α * τ)α是材料的熱擴散系數(shù)典型金屬大概在1e-4平方米每秒量級τ是脈沖寬度。如果脈寬取10納秒算下來熱擴散深度大約是2微米左右。這個數(shù)字很關(guān)鍵它直接決定你網(wǎng)格在表面附近要怎么加密——單元尺寸最好在熱擴散深度的五分之一到十分之一以下不然溫度梯度根本解析不出來。納秒脈沖的能量沉積過程光先被材料表面極淺的一層吸收金屬的光吸收深度通常是幾十納米量級然后熱量靠熱傳導(dǎo)往深處走。脈寬只有幾十納秒熱傳導(dǎo)來不及把能量送到很遠(yuǎn)的地方所以燒蝕區(qū)主要集中在光斑半徑和熱擴散深度限定的范圍內(nèi)。材料表面溫度急劇升高達(dá)到沸點甚至超過沸點后表層材料發(fā)生汽化蒸汽帶走了大量能量表面開始后退這就是燒蝕。這里有個很多人忽略的點燒蝕不僅僅是溫度超過了某個閾值那么簡單。納秒激光燒蝕在物理上是由蒸發(fā)動力學(xué)控制的表面溫度不是到了沸點就停在沸點而是可以繼續(xù)升高蒸發(fā)速率隨溫度按指數(shù)規(guī)律猛增。所以你在模型里不能靠溫度到了某個值就刪掉網(wǎng)格這種粗暴方式而是要跟蒸發(fā)相關(guān)的熱流損失和表面后退速度掛上鉤。1.2 移動網(wǎng)格在模型里到底扮演什么角色移動網(wǎng)格Moving Mesh在COMSOL里實際是基于ALE方法也就是任意拉格朗日-歐拉方法。這套方法的通俗理解是網(wǎng)格節(jié)點可以在空間里動但動的速度不必等于材料質(zhì)點的速度這樣可以保留清晰的材料界面同時避免網(wǎng)格像純拉格朗日方法那樣劇烈扭曲。在激光燒蝕這個場景里燒蝕表面在往里退如果你不處理這個邊界那幾何就一直保持原始形狀算到后面表面溫度早就超出物理范圍了。移動網(wǎng)格的職責(zé)就是跟隨燒蝕邊界往回退同時保持網(wǎng)格質(zhì)量、避免單元翻轉(zhuǎn)。要注意的是COMSOL里有兩個容易搞混的接口一個是移動網(wǎng)格Moving Mesh對應(yīng)物理場接口名是ale另一個是變形幾何Deformed Geometry對應(yīng)接口名是dg。做燒蝕模擬絕大多數(shù)情況下用移動網(wǎng)格更合適。區(qū)別在于移動網(wǎng)格里材料域還是那個材料域只是網(wǎng)格在空間里變形適應(yīng)新的邊界位置變形幾何則更偏向幾何本身被重新定義常用于制造過程仿真或形狀優(yōu)化。燒蝕模擬里材料是不斷被去除的但你不需要真的把材料域挖掉只需要讓燒蝕邊界動態(tài)后退移動網(wǎng)格夠用了。1.3 結(jié)果不理想先分清是哪一類不理想排查問題最忌諱的就是漫無目的地亂試。溫度場不理想先說清楚到底哪里不理想我把常見的情況歸成幾類溫度分布的形狀不對比如應(yīng)該是高斯鼓包狀的等溫線結(jié)果變成了扁平餅狀或者出現(xiàn)奇怪的環(huán)形振蕩。溫度量級離譜比如中心溫度算出來幾百萬開爾文或者算完整個工件還是室溫沒反應(yīng)。燒蝕坑的形狀不對比如應(yīng)該是近似高斯截面的凹坑結(jié)果坑壁角度怪異、邊緣翹起或者深度明顯偏淺偏深。計算中途報錯基本都是網(wǎng)格畸變、雅可比矩陣非正定、在某個時間步不收斂這一類。多脈沖場景下溫度場不隨激光移動而移動掃過的路徑上有殘留的熱積累異常。不同的問題對應(yīng)完全不同的根源。溫度量級離譜大概率是熱源單位、吸收率或時間步長的問題形狀不對多半是熱源空間分布或者網(wǎng)格分辨率的問題網(wǎng)格畸變是移動網(wǎng)格設(shè)置的問題。下面幾個章節(jié)我會把每個方向的具體細(xì)節(jié)展開講你對照自己的現(xiàn)象去查就行。2. 移動網(wǎng)格設(shè)置80%的結(jié)果不理想都出在這里移動網(wǎng)格是燒蝕模擬里最容易出問題的一環(huán)也是最勸退的一環(huán)。很多人的溫度場不對根本原因其實是網(wǎng)格動得不對導(dǎo)致傳熱計算的幾何坐標(biāo)亂掉了。這一章我把移動網(wǎng)格的選型、邊界配置和網(wǎng)格質(zhì)量控制拆開講。2.1 移動網(wǎng)格的數(shù)學(xué)基礎(chǔ)與幾何配置ALE方法的核心思想是網(wǎng)格的位移場由你指定COMSOL默認(rèn)在計算域內(nèi)部用平滑算法來分配網(wǎng)格位移。當(dāng)你給燒蝕表面指定了一個向內(nèi)的位移或速度移動網(wǎng)格接口就會用Laplace平滑或Winslow平滑把邊界運動的影響擴散到整個域內(nèi)部越靠近運動邊界網(wǎng)格變形越大遠(yuǎn)處逐漸衰減。在COMSOL的設(shè)置中幾何配置要先保證移動網(wǎng)格接口的變形域覆蓋所有可能發(fā)生幾何變形的區(qū)域。如果你的模型里燒蝕只發(fā)生在光斑附近的局部區(qū)域你可以把變形域縮小到激光作用區(qū)周圍這樣能大幅降低網(wǎng)格畸變風(fēng)險計算量也小很多。遠(yuǎn)處不需要變形的區(qū)域直接設(shè)成固定網(wǎng)格別讓它們跟著動。實際模型里我習(xí)慣這樣配置幾何整個材料域是一個長方形二維軸對稱或者長方體三維激光從上方入射。頂部表面是燒蝕面設(shè)置為變形邊界底面和側(cè)面設(shè)置為固定邊界如果模型中間有對稱軸對稱軸也要做相應(yīng)處理。變形域單獨劃出一個局部區(qū)域比如以光斑為中心、半徑3到5倍光斑半徑、深度3到5倍熱擴散深度的子域這個子域內(nèi)網(wǎng)格可以自由變形外面的區(qū)域網(wǎng)格保持不動。這里有個很多人容易踩的坑變形域范圍太小網(wǎng)格變形集中在很窄的條帶里很快單元就擰成麻花了變形域范圍太大邊界移動的位移被平滑掉大半燒蝕表面的實際內(nèi)退量跟計算值對不上溫度場自然不準(zhǔn)。經(jīng)驗上先取3倍熱擴散深度作為變形域厚度比較穩(wěn)妥然后根據(jù)網(wǎng)格質(zhì)量再做調(diào)整。2.2 哪個邊界該動、哪個不該動邊界配置是移動網(wǎng)格的核心操作也是出錯最密集的地方。燒蝕表面要動這是肯定的但一動就有講究你給它指定的是指定網(wǎng)格位移還是指定網(wǎng)格速度。指定網(wǎng)格位移Specified Mesh Displacement適合每個時間步內(nèi)位移量明確的情況。燒蝕表面每個點的位移等于燒蝕速率對時間的積分。這個方式的好處是直觀、穩(wěn)定只要不出現(xiàn)大規(guī)模網(wǎng)格翻轉(zhuǎn)基本不會出幺蛾子。壞處是一旦積分誤差累積位移量會漂移。指定網(wǎng)格速度Specified Mesh Velocity適合燒蝕速率作為表面溫度函數(shù)直接計算的情況。你把蒸發(fā)模型算出來的表面后退速度直接給到邊界法向網(wǎng)格每步按這個速度運動。這個方式更貼近物理但對求解器的時間步長和容差要求更高速度值有波動時網(wǎng)格容易抖。我建議新手先用指定網(wǎng)格速度因為它跟蒸發(fā)模型的耦合路徑最短邏輯最清晰。具體操作上在移動網(wǎng)格接口里給燒蝕邊界加一個指定網(wǎng)格速度節(jié)點速度方向設(shè)為邊界法向大小由蒸發(fā)模型給出。注意COMSOL的網(wǎng)格速度方向需要使用邊界法向變量比如-t_MN.nx或者類似的法向表達(dá)式別自己寫成分量形式否則在弧形邊界上速度方向會亂。其余邊界的處理原則是能不動的都別動。底面和側(cè)壁直接設(shè)為固定邊界對稱面如果沿法向沒有位移也是固定。溫度場求解的邊界條件跟網(wǎng)格邊界條件是兩套系統(tǒng)不要混淆。網(wǎng)格邊界管的是幾何邊界的運動物理邊界管的是熱流和溫度的約束二者需要在同一幾何邊界上共存但設(shè)置是獨立的。2.3 網(wǎng)格變形控制與單元質(zhì)量移動網(wǎng)格的網(wǎng)格質(zhì)量控制決定了你的計算能跑多遠(yuǎn)。就算物理模型全對網(wǎng)格一旦翻轉(zhuǎn)扭曲求解器分分鐘報錯。COMSOL的移動網(wǎng)格接口里有網(wǎng)格平滑設(shè)置默認(rèn)是Laplace平滑對大的邊界位移適應(yīng)性較差容易出現(xiàn)單元被壓扁的情況。我實際用下來Winslow平滑在大變形問題里表現(xiàn)更好因為它考慮了各向異性的網(wǎng)格變形如果燒蝕深度很大還可以試試超彈性平滑。超彈性的思路是把網(wǎng)格當(dāng)成一個彈性體用應(yīng)變能最小化來分配變形對極端變形的容忍度更高但每次迭代的計算量也更大。這里還要補充一個預(yù)先做好網(wǎng)格的準(zhǔn)備在燒蝕表面附近也就是光斑作用區(qū)域網(wǎng)格要加密而且要有意識地把網(wǎng)格設(shè)計成過渡平緩的樣式。如果粗糙網(wǎng)格和加密網(wǎng)格之間尺寸跳變太劇烈變形后的過渡區(qū)很容易出現(xiàn)負(fù)體積單元。表面用邊界層網(wǎng)格是好習(xí)慣但邊界層的層數(shù)和厚度要跟移動網(wǎng)格兼容——邊界層太密太細(xì)變形時最外層單元會最先扭曲。我常用的網(wǎng)格策略是表面最大單元尺寸取光斑半徑的二十分之一左右深度方向用邊界層首層厚度取熱擴散深度的十分之一層數(shù)以8到10層為基準(zhǔn)增長率1.2到1.3。二維軸對稱模型網(wǎng)格數(shù)一般控制在2萬到5萬三維模型控制在30萬到80萬既保證精度又不會算到天荒地老。3. 納秒脈沖激光熱源建模結(jié)果偏差的頭號來源如果移動網(wǎng)格設(shè)置都正常溫度場還是不對勁那問題十有八九在熱源上。熱源建模的細(xì)節(jié)相當(dāng)多我單獨列一章講因為這里一個單位錯誤、一個吸收率參數(shù)不對溫度場能差出好幾個數(shù)量級。3.1 熱源項的函數(shù)表達(dá)與關(guān)鍵參數(shù)納秒脈沖激光在材料內(nèi)部的熱源常用體熱源形式。假設(shè)激光沿z軸入射光斑在x-y平面上是高斯分布光束沿x方向以速度v移動那么熱源表達(dá)式可以寫成Q(x,y,z,t) (1 - R) * I0 * exp(-2 * ((x - v*t)^2 y^2) / ω0^2) * exp(-α_abs * z) * pulse(t)各項含義R是材料表面反射率I0是激光峰值功率密度單位W/m2ω0是高斯光斑的束腰半徑α_abs是材料對該波長激光的吸收系數(shù)單位1/mpulse(t)是脈沖的時間包絡(luò)函數(shù)。很多人的單位問題就出在這里。COMSOL默認(rèn)用國際單位制長度單位是米時間單位是秒。如果你的光斑半徑習(xí)慣用微米、功率密度習(xí)慣用W/cm2換算錯一位就是100倍甚至10000倍的誤差。我見過一個模型I0按W/cm2填進去而COMSOL按W/m2處理中心溫度直接爆到上百萬開爾文。建議在全局參數(shù)里把單位都顯式換算好比如I0 1e12[W/m^2]這種寫法讓COMSOL自己檢查單位一致性。反射率R這個參數(shù)很傷人。金屬在室溫下對近紅外激光的反射率通常很高鋁大概0.9左右銅更高鋼大概0.6到0.7。但材料溫度升高后反射率會下降特別是接近熔點和沸點時吸收率顯著上升。如果模型里用固定反射率算出來的能量沉積往往偏低或偏高。嚴(yán)格做法是把反射率設(shè)成溫度的函數(shù)比如用一個分段函數(shù)在室溫到熔點的區(qū)間內(nèi)線性插值。如果要先跑通模型可以先用一個折中的固定值比如鋼取0.35到0.4的有效吸收率跟實驗對一下再修正。3.2 脈沖時序的處理方式納秒脈沖的時間包絡(luò)有幾種常見寫法矩形脈沖、高斯脈沖、或者更真實的平頂-高斯混合模式。矩形脈沖最簡單直接用(t - n*T_period)落在脈寬τ內(nèi)時為1否則為0。但真實納秒激光器輸出的脈沖形狀更接近高斯或平頂高斯特別是調(diào)Q激光器出來的脈沖時間波形大致是高斯形。在COMSOL里實現(xiàn)多脈沖序列最簡單的辦法是用解析函數(shù)加取模運算。比如脈沖周期為T_period單脈沖寬度為τ_pulse那么第n個脈沖的起始時刻是n*T_period。定義一個變量t_pulse mod(t, T_period)然后脈沖包絡(luò)pulse(t)等于在t_pulse τ_pulse時為1矩形或高斯形狀。高斯包絡(luò)可以寫成pulse(t) exp(-((t_pulse - τ_pulse/2) / (τ_pulse/4))^2)只要t_pulse τ_pulse成立。這個寫法要注意mod函數(shù)在脈沖開始瞬間會有跳變?nèi)绻蠼馄鲿r間步長太大容易漏掉脈沖起始沿導(dǎo)致能量缺失。這時候可以用COMSOL的事件接口來精確控制脈沖的開啟和關(guān)閉在脈沖開始和結(jié)束時刻強制求解器插入時間步。多脈沖模擬如果結(jié)果總覺得能量不夠檢查一下是不是時間步長跨過了脈沖沿。3.3 材料對激光的吸收與相變潛熱的處理體熱源表達(dá)式里吸收系數(shù)α_abs對不同材料差異很大。金屬對納秒激光的吸收深度非常淺典型值在10到50納米這遠(yuǎn)小于熱擴散深度。如果把真實吸收深度放進體熱源你需要在表面劃分厚度納微米級的網(wǎng)格才能解析這個能量沉積層計算代價極高。實際工程處理中很多文獻(xiàn)直接用表面熱通量來近似等效成熱流邊界條件。這樣做的物理依據(jù)是吸收深度遠(yuǎn)小于熱擴散深度時能量沉積可以看成表面事件對溫度場的差別在微米尺度之外可忽略。如果你一定要用體熱源可以把吸收深度人為放寬到跟網(wǎng)格尺度相當(dāng)?shù)某潭鹊ＷC總吸收能量不變。也就是說吸收系數(shù)α_abs取1微米量級I0相應(yīng)調(diào)大使得積分吸收的總功率跟物理情況一致。這種做法是數(shù)值上的合理近似不是物理造假前提是你能評估誤差范圍。相變潛熱這一塊納秒激光燒蝕過程中材料會經(jīng)歷熔化、汽化、甚至等離子體化。但嚴(yán)格模擬所有這些相變在COMSOL里非常復(fù)雜。最實用的做法是把熔化潛熱考慮到材料的熱容項里用一個等效熱容Cp_eff Cp L_m * delta(T - T_m)來處理固液相變汽化潛熱則通過蒸發(fā)模型的熱流損失來體現(xiàn)而不是顯式引入一個相變溫度區(qū)間。汽化熱流損失可以寫成q_evap m_dot * L_v其中m_dot是蒸發(fā)質(zhì)量流率L_v是汽化潛熱。這個損失項加到燒蝕表面的熱邊界條件里同時m_dot通過Hertz-Knudsen方程由表面溫度決定。這樣一來溫度越高蒸發(fā)越快帶走的熱量越多自然形成一種負(fù)反饋防止溫度無限制飆升。4. 一套可以直接上手的實操流程二維軸對稱示例前面講了一大堆原理和坑這章給一套可以照做的完整流程。我以一個二維軸對稱模型的單脈沖納秒激光燒蝕金屬靶材為例這套配置在學(xué)校工作站和普通臺式機上都能跑得動大概半小時到兩小時能出結(jié)果。4.1 幾何、材料與物理場選擇打開COMSOL新建模型向?qū)нx擇二維軸對稱添加兩個物理場接口固體傳熱ht和移動網(wǎng)格ale。研究選瞬態(tài)。幾何直接畫一個矩形代表材料截面比如半徑200微米、厚度100微米。這里維度要說明一下二維軸對稱模型的對稱軸是z軸r是徑向坐標(biāo)激光沿z軸方向入射光斑中心在對稱軸上。如果你的激光是移動的二維軸對稱就不適用了需要切到二維平面模型或者三維模型。移動激光的場景我放到后面講單脈沖靜止光斑先用軸對稱把核心邏輯跑通。材料參數(shù)給一個典型金屬的值以304不銹鋼為例導(dǎo)熱系數(shù)15 W/(m·K)密度7900 kg/m3比熱500 J/(kg·K)熔點1700 K沸點3100 K熔化潛熱2.6e5 J/kg汽化潛熱6.0e6 J/kg。熱擴散系數(shù)算下來約3.8e-6 m2/s10納秒脈沖熱擴散深度約0.4微米。注意這里導(dǎo)熱系數(shù)我給了室溫值。真實情況導(dǎo)熱系數(shù)隨溫度升高會下降不銹鋼從室溫到高溫導(dǎo)熱系數(shù)變化很明顯。如果做定量對比建議設(shè)置一個線性溫度相關(guān)導(dǎo)熱系數(shù)比如從15降到25 W/(m·K)再降到更高溫時的28左右這種趨勢具體數(shù)值按材料手冊插值。4.2 移動網(wǎng)格與熱源耦合的具體操作移動網(wǎng)格接口中選擇整個材料域作為變形域。燒蝕表面z0邊界添加指定網(wǎng)格速度節(jié)點設(shè)置法向速度為蒸發(fā)模型計算出的表面后退速度。表面后退速度可以從Hertz-Knudsen方程出發(fā)用飽和蒸氣壓和溫度的關(guān)系來推導(dǎo)v_surface m_dot / ρ其中 m_dot P_v(T) * sqrt(M / (2πRT))。P_v用Clausius-Clapeyron近似表達(dá)P_v(T) P_atm * exp((L_v * M / R) * (1/T_boil - 1/T))把這一串簡化后在COMSOL的變量定義里寫好就可以。擔(dān)心表達(dá)式太復(fù)雜的話第一版模型可以先測試一個簡化版本表面溫度超過沸點時給一個恒定的表面后退速度比如1到5 m/s。這個數(shù)量級跟納秒激光燒蝕的典型后退速度吻合。先跑通再慢慢把蒸發(fā)模型完善進去。熱源設(shè)置這一步如果是單脈沖軸對稱模型可以用一個熱通量邊界條件施加在燒蝕表面這個方式最簡單。熱通量表達(dá)式為q_laser I0 * exp(-2*r^2/ω0^2) * pulse(t)配合燒蝕表面同時放一個熱通量節(jié)點減去蒸發(fā)冷卻項再加上輻射和對流損失。蒸發(fā)冷卻項就是-m_dot * L_v輻射項用表面發(fā)射率算對流項可以用經(jīng)驗對流系數(shù)但納秒時間尺度對流換熱對結(jié)果影響很小可以忽略。COMSOL里有個隱形的坑當(dāng)移動網(wǎng)格啟用時固體傳熱默認(rèn)按歐拉描述還是ALE描述處理取決于你是否勾選了移動網(wǎng)格坐標(biāo)選項。在固體傳熱接口的物理場設(shè)置里找到坐標(biāo)系的運動或者移動網(wǎng)格選項必須打開這樣控制方程里會自動帶上網(wǎng)格運動帶來的對流項。很多人網(wǎng)格動了但傳熱方程沒感知到幾何在動溫度場當(dāng)然不對。4.3 求解器與時間步長的設(shè)置納秒脈沖激光燒蝕計算的時間跨度分布很寬單脈沖過程在納秒尺度但多脈沖之間可能是微秒甚至毫秒間隔。求解器設(shè)置要分情況。單脈沖模擬時間范圍可以從0到幾個微秒比如0到2微秒。關(guān)鍵是要保證脈沖期間每個納秒都有好幾個時間步。COMSOL用BDF向后差分時間步進默認(rèn)是自適應(yīng)步長理論上會自己加密但實際中如果容差放松了BDF會得出一個偏大的時間步導(dǎo)致脈沖能量注入不連續(xù)。我在研究設(shè)置里會把時間步序列直接給定比如range(0, 0.2e-9, 20e-9) range(20e-9, 1e-8, 200e-9) range(220e-9, 5e-8, 2e-6)這樣前20納秒用0.2納秒的步長抓脈沖后面用50納秒甚至更大的步長算冷卻。這個方法很笨但很穩(wěn)。求解器里的相對容差建議設(shè)嚴(yán)一些1e-4起最好到1e-5。納秒激光的溫度梯度非常陡絕對容差和相對容差太松溫度峰值會被明顯平滑掉。我見過有人用默認(rèn)的0.01跑出來溫度場分布總是偏圓潤怎么調(diào)熱源都調(diào)不出尖銳的高斯峰最后發(fā)現(xiàn)就是容差太松。4.4 后處理溫度場結(jié)果的提取與評估計算完先別急著截圖有幾樣?xùn)|西一定要提取出來看第一是表面中心點的溫度時間曲線。這個曲線能直接告訴你峰值溫度是不是合理。對金屬而言納秒脈沖的峰值表面溫度通常在幾千開爾文量級如果遠(yuǎn)低于沸點說明能量不夠或者吸收率太低如果遠(yuǎn)超上萬K說明能量輸入過量或者蒸發(fā)冷卻沒起作用。第二是燒蝕深度隨時間的變化。在移動網(wǎng)格接口里可以直接提取燒蝕表面的位移路徑上取軸向?qū)ΨQ點繪制位移時間曲線。這個值跟激光參數(shù)和材料參數(shù)之間的量級關(guān)系可以做個估算如果平均燒蝕深度是幾微米量級那跟文獻(xiàn)數(shù)據(jù)差不多可以對上。第三是瞬時溫度場云圖。建議在脈沖峰值時刻、脈沖結(jié)束后10納秒、100納秒這幾個時間點各存一幀。對比這幾幀可以看到熱擴散是否正常如果等溫線在脈沖結(jié)束后還在劇烈向外擴說明熱擴散系數(shù)設(shè)置偏大如果幾乎不擴散說明導(dǎo)熱系數(shù)設(shè)置有問題。5. 溫度場不理想對照這份排查清單這一章我把自己踩過的坑和幫別人看模型時的高頻問題整理成一份排查清單。你跑出不對的溫度場時按順序過一遍大部分問題能定位。5.1 溫度場分布形狀不對如果你期待高斯形溫度鼓包出來的溫度場卻是平頭狀或者環(huán)形先檢查你的熱源表達(dá)式里的空間分布函數(shù)。常見錯誤是exp(-2*r^2/ω0^2)里的光斑半徑ω0取了直徑值或者單位是微米沒有換算成米。這類錯誤的表現(xiàn)一般是光斑區(qū)域偏大或者偏小等溫線間距跟預(yù)期不一致。另一個常見原因是對流項丟失。前面提過固體傳熱接口里要勾選移動網(wǎng)格坐標(biāo)選項如果沒勾選在邊界移動的區(qū)域溫度方程會出錯等溫線會扭曲變形。這種變形在燒蝕坑附近最明顯溫度場會出現(xiàn)不對稱或者被拉伸的奇怪形狀。如果溫度分布出現(xiàn)振蕩的環(huán)狀條紋那基本是網(wǎng)格分辨率不夠溫度梯度沒法在兩個相鄰網(wǎng)格單元之間被平滑表達(dá)。解決方法是加密表面邊界層網(wǎng)格同時把相對容差收緊。5.2 溫度量級離譜溫度峰值過高一般逃不出這幾個原因激光峰值功率密度I0的單位換算錯誤最常見W/cm2和W/m2差了4個數(shù)量級。反射率設(shè)成0實際材料吸收了全部入射能量。時間步長太大脈寬時間內(nèi)只算了一個步能量在一步之內(nèi)全部注入沒有熱傳導(dǎo)的時間。蒸發(fā)冷卻沒設(shè)或者蒸發(fā)模型沒起作用表面溫度失去了最重要的散熱通道。溫度峰值過低同樣好查I0數(shù)量級算錯光斑半徑設(shè)大了能量密度攤薄吸收系數(shù)偏低脈沖時間包絡(luò)的積分不等于單脈沖能量。一個實用的自檢方法用單脈沖能量密度單位J/cm2來估算。如果你知道激光器的單脈沖能量和光斑面積峰值功率密度I0跟單脈沖能量密度的關(guān)系是I0 ≈ F / τ矩形脈沖或者I0 2F/(τ*sqrt(π))高斯脈沖。先把I0量級算對再填進COMSOL。5.3 網(wǎng)格畸變導(dǎo)致計算終止網(wǎng)格畸變是移動網(wǎng)格模擬的頭號殺手。報錯信息通常是找不到一致的初始值或者在為某些網(wǎng)格單元求解時發(fā)生錯誤。出現(xiàn)這種問題排查順序是先看報錯發(fā)生前最后幾步的網(wǎng)格質(zhì)量在結(jié)果里創(chuàng)建一個單元質(zhì)量繪圖可視化最后成功時間步的網(wǎng)格。如果發(fā)現(xiàn)某些單元的偏斜度接近0或者出現(xiàn)負(fù)值那就是這些單元導(dǎo)致計算崩潰。檢查變形域是否太小變形全擠在幾個單元里。把變形域擴大讓網(wǎng)格位移有足夠的空間攤開。修改平滑設(shè)置。Laplace換成Winslow或超彈性對比網(wǎng)格質(zhì)量的變化。增大最小單元質(zhì)量指標(biāo)或者開啟自動重新劃分網(wǎng)格如果版本支持。這個功能允許在網(wǎng)格質(zhì)量惡化到一定程度時自動重新生成網(wǎng)格把舊解插值到新網(wǎng)格上繼續(xù)計算。時間步長縮小防止單步網(wǎng)格位移過大。新手可以給網(wǎng)格速度乘一個松弛因子比如實際速度的0.8倍先跑通了再逐步逼近真實值。5.4 燒蝕坑的形貌異常燒蝕坑太淺、太深或者邊緣形狀不合理原因不完全在網(wǎng)格往往在蒸發(fā)模型和熱源耦合。坑太淺說明表面后退速度偏小可以檢查蒸發(fā)模型里飽和蒸氣壓的溫度依賴是否太弱或者表面后退速度被低估了?？犹钫f明速度偏大可能是汽化潛熱L_v取小了蒸發(fā)帶走的熱量偏少溫度降不下來?？颖谛螤钊绻腹獍叩母咚狗植歼吘壞芰勘旧砭偷驼?yīng)該是平滑的圓坑。如果坑邊緣出現(xiàn)翹起或者向內(nèi)翻大概率是表面網(wǎng)格運動方向出了問題。指定網(wǎng)格速度時務(wù)必確認(rèn)方向是材料內(nèi)部法向別設(shè)反了成了往外鼓。5.5 移動激光場景溫度場跟不上光斑位置最后提一下移動激光的情況。標(biāo)題里專門提到移動的納秒脈沖激光這一類模型溫度場的典型問題就是熱積累位置和光斑位置對不上。移動激光需要在空間坐標(biāo)里把光束位置的時變性寫進去。如果激光沿x方向以速度v移動熱源表達(dá)式里的空間分布函數(shù)要寫成(x - v*t)而不是x。COMSOL里可以在全局參數(shù)設(shè)速度然后用變量引用。這一步不難但很多人容易在某一個維度上漏掉時間項導(dǎo)致溫度場中心跟光斑錯位。還有一個細(xì)節(jié)是網(wǎng)格變形方向。移動激光的燒蝕坑是沿著掃描方向的溝槽網(wǎng)格變形不再只是軸向的而是有橫向移動的趨勢。如果變形域設(shè)得太小網(wǎng)格畸變會更快出現(xiàn)。建議移動掃描場景下變形域的寬度取光斑半徑的4到5倍深度取燒蝕深度的3倍以上同時把網(wǎng)格速度平滑一下避免因為光束移動和燒蝕速度疊加出局部的大變形。6. 我自己踩過的幾個坑順便說點建議最后隨便聊幾句實操中的體會。有個細(xì)節(jié)我一直覺得值得提醒做這類多物理場耦合模擬千萬別一開始就追求全細(xì)節(jié)。第一版模型盡量簡單粗暴比如固定反射率、矩形脈沖、簡化蒸發(fā)模型跑通出結(jié)果再逐項增加復(fù)雜度。我見過太多人第一步就把反射率做成溫度函數(shù)、飽和蒸氣壓公式寫了幾十行、等離子體效應(yīng)都塞進去結(jié)果模型報錯都找不到是哪個公式的問題。仿真也是工程分步走永遠(yuǎn)比一步到位靠譜。另外一個技巧是惡意測試故意把激光功率密度調(diào)大10倍看溫度場是不是跟著飆升、燒蝕深度是不是明顯加深。如果響應(yīng)方向不對說明模型里某項耦合順序錯了這比對著結(jié)果圖紙排查快得多。還有就是要養(yǎng)成看一眼網(wǎng)格質(zhì)量的習(xí)慣。很多結(jié)果不理想其實在計算過程中就已經(jīng)埋了雷只是到了后面才爆。跑完第一步就生成單元質(zhì)量圖看看比跑完整個模型再回頭找原因節(jié)省幾小時起步。納秒脈沖激光燒蝕的COMSOL模擬難就難在時間尺度、空間尺度和多物理場耦合這三個維度同時起作用。溫度場不理想往根上說要么是物理過程沒表達(dá)對要么是數(shù)值參數(shù)沒配合好。按這篇的順序把移動網(wǎng)格和熱源建模重新過一遍你對模型的控制力會上一個臺階。等溫度場做對了后面再往上加材料相變、流體流動甚至等離子體模型都只是錦上添花。