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