1. 項(xiàng)目初印象與整體定位做生物傳熱仿真這幾年COMSOL一直是我主力工具之一。這個(gè)人體皮膚溫度場(chǎng)計(jì)算模型聽(tīng)起來(lái)好像就是幾層材料疊在一起算個(gè)溫度分布但真上手做一遍就會(huì)發(fā)現(xiàn)從幾何搭建、參數(shù)賦值到求解收斂每一步都有不少容易被忽略的細(xì)節(jié)。我最早接觸這個(gè)課題是想幫做低溫燙傷研究的同學(xué)建立一個(gè)皮膚受熱后的溫度響應(yīng)模型——實(shí)驗(yàn)不可能在人身上直接做而離體皮膚和活體皮膚的熱響應(yīng)又有明顯差異仿真就成了一個(gè)合理的替代手段。這篇內(nèi)容會(huì)把我從零搭建該模型的完整過(guò)程梳理一遍包括控制方程的選取邏輯、皮膚分層結(jié)構(gòu)的幾何建模細(xì)節(jié)、邊界條件的物理設(shè)定、網(wǎng)格劃分與求解器調(diào)試的實(shí)踐經(jīng)驗(yàn)以及模型后續(xù)可以擴(kuò)展的方向。不論你是剛接觸COMSOL的學(xué)生還是做醫(yī)工交叉項(xiàng)目、可穿戴設(shè)備熱設(shè)計(jì)、熱療方案評(píng)估的工程師這套建模思路都可以直接參考。模型本身不算復(fù)雜但涉及的關(guān)鍵概念——比如血液灌注項(xiàng)的物理含義、薄層的網(wǎng)格處理、瞬態(tài)求解的時(shí)間尺度匹配——都是這一領(lǐng)域里最容易踩坑的點(diǎn)。1.1 這個(gè)模型到底能做什么皮膚溫度場(chǎng)計(jì)算模型本質(zhì)上是在求解生物組織的傳熱方程得到皮膚表面和內(nèi)部各深度處的溫度分布。在工程和醫(yī)學(xué)場(chǎng)景里這個(gè)結(jié)果有多種用途。最基本的是評(píng)估皮膚在特定環(huán)境下的熱舒適性——比如判斷一塊加熱坐墊在多大功率下是溫暖而不是燙傷或者評(píng)估空調(diào)房里長(zhǎng)時(shí)間對(duì)著冷風(fēng)的體感。往醫(yī)學(xué)方向走它可以用于低溫燙傷風(fēng)險(xiǎn)評(píng)估、熱敷理療的溫度控制、射頻和激光治療的熱劑量估算。再往工業(yè)方向可穿戴設(shè)備、暖手寶、暖腳貼這類(lèi)貼身發(fā)熱產(chǎn)品的設(shè)計(jì)幾乎都繞不開(kāi)皮膚接觸面的溫度分布計(jì)算。我實(shí)際做過(guò)的一個(gè)場(chǎng)景是把一個(gè)恒溫加熱片貼在皮膚表面計(jì)算不同設(shè)定溫度下皮膚內(nèi)部溫度隨時(shí)間的變化。這種模擬如果只用實(shí)驗(yàn)試錯(cuò)成本高且不易覆蓋所有工況但仿真建模一次就可以把各種參數(shù)組合掃描一遍。模型跑通之后你會(huì)發(fā)現(xiàn)它還能反過(guò)來(lái)指導(dǎo)實(shí)驗(yàn)設(shè)計(jì)——比如告訴你傳感器應(yīng)該埋在哪個(gè)深度、多厚的組織才能測(cè)到穩(wěn)定信號(hào)。1.2 與普通固體傳熱的本質(zhì)區(qū)別很多做結(jié)構(gòu)傳熱的人第一次接觸生物傳熱容易把皮膚當(dāng)成一塊普通固體來(lái)算。這個(gè)思路能算出一個(gè)溫度場(chǎng)但結(jié)果和實(shí)測(cè)往往偏差很大根本原因是活體組織有一個(gè)普通材料沒(méi)有的機(jī)制——血液灌注。血液在體內(nèi)循環(huán)時(shí)會(huì)不斷把核心區(qū)的熱量帶到外周組織相當(dāng)于一個(gè)遍布全場(chǎng)的溫度調(diào)節(jié)熱源。當(dāng)你用熱水袋烤皮膚時(shí)不僅是皮膚本身在導(dǎo)熱血液也在不斷沖走局部熱量這個(gè)效應(yīng)直接決定了組織最終能達(dá)到多高的溫度。另一個(gè)明顯區(qū)別是邊界條件的復(fù)雜度。普通固體的外表面往往直接給定熱流或?qū)α鲹Q熱系數(shù)而人體皮膚表面有出汗蒸發(fā)、熱輻射、衣物隔熱等多重機(jī)制深部邊界又要錨定在核心體溫上。把這些耦合關(guān)系理清楚模型才會(huì)接近真實(shí)。COMSOL在這類(lèi)問(wèn)題上的優(yōu)勢(shì)恰恰在于可以靈活地添加物理場(chǎng)、修改方程源項(xiàng)、調(diào)整邊界條件而不需要像某些純編程方案那樣從零推導(dǎo)離散格式并自己處理網(wǎng)格。這也是我最終選擇COMSOL來(lái)做這個(gè)模型的主要原因。2. 控制方程與傳熱參數(shù)的設(shè)定邏輯建模第一步不是急著畫(huà)幾何而是把數(shù)學(xué)模型敲定。皮膚溫度場(chǎng)計(jì)算用的核心方程是Pennes生物傳熱方程。這個(gè)方程從1948年被提出到現(xiàn)在依然是絕大多數(shù)生物傳熱仿真包括熱療、冷凍、燒傷評(píng)估等的基礎(chǔ)框架。它不是唯一的選擇工程上也有用Wulff方程或更復(fù)雜的多層多孔介質(zhì)模型的但對(duì)皮膚溫度場(chǎng)這個(gè)尺度來(lái)說(shuō)Pennes方程在計(jì)算代價(jià)和物理準(zhǔn)確性之間是最平衡的。2.1 Pennes生物傳熱方程逐項(xiàng)拆解Pennes方程的基本形式是ρc·?T/?t ?·(k?T) ρb·cb·ωb·(Ta - T) Qmet如果你第一次看這個(gè)式子我建議把它拆成四個(gè)物理塊來(lái)理解。等式左邊是單位體積組織升溫所需的能量變化率。右邊第一項(xiàng)是熱傳導(dǎo)項(xiàng)和普通導(dǎo)熱方程一致描述熱量在組織間的擴(kuò)散。右邊的第二項(xiàng)是血液灌注項(xiàng)其中ρb和cb是血液的密度和比熱容ωb是血液灌注率Ta是動(dòng)脈血溫度T是組織局部溫度。這里有個(gè)生活化類(lèi)比想分享給大家——血液灌注項(xiàng)就像家里暖氣的管路管路里的熱水不停流動(dòng)把熱量帶到房間各處房間溫度越低于熱水溫度進(jìn)來(lái)補(bǔ)充的熱量就越多。第三項(xiàng)Qmet是代謝產(chǎn)熱表示組織自身新陳代謝產(chǎn)生的小量熱量安靜狀態(tài)下數(shù)值不大但在運(yùn)動(dòng)或炎癥狀態(tài)下會(huì)明顯上升。理解這個(gè)方程的關(guān)鍵是對(duì)血液灌注項(xiàng)的源項(xiàng)屬性有清晰認(rèn)識(shí)。它在方程里形式上不是時(shí)間或空間的微分項(xiàng)而是溫度T的代數(shù)項(xiàng)——所以它的物理效果是在每個(gè)計(jì)算點(diǎn)上都相當(dāng)于一個(gè)恒溫緩沖器局部溫度比血液溫度低得越多被加熱的速率越快反過(guò)來(lái)局部溫度比血液溫度高時(shí)它又變成散熱器。正因?yàn)檫@一項(xiàng)的存在皮膚在受到外部熱源加熱時(shí)溫度不會(huì)無(wú)限攀升而是會(huì)在某個(gè)位置達(dá)到平衡。實(shí)測(cè)數(shù)據(jù)也印證了這一點(diǎn)——同樣是45℃的接觸面活體組織因?yàn)檠汗嘧⒌拇嬖趦?nèi)部溫度峰值和壞死程度都比離體組織要低。2.2 各項(xiàng)參數(shù)的取值范圍與推導(dǎo)依據(jù)參數(shù)取值決定了模型的可靠性這也是實(shí)操過(guò)程中最容易出問(wèn)題的地方。皮膚各層的參數(shù)在專(zhuān)業(yè)文獻(xiàn)中是比較統(tǒng)一的但初學(xué)者極易在單位換算上栽跟頭。下面是我常用的一組參考值基本來(lái)自生物熱物理的經(jīng)典教科書(shū)和近年的相關(guān)論文參數(shù)表皮真皮皮下組織單位導(dǎo)熱系數(shù) k0.2350.3350.190W/(m·K)密度 ρ120012001000kg/m3比熱容 c359032202500J/(kg·K)血液灌注率 ωb0無(wú)血管0.00120.0007s?1代謝產(chǎn)熱 Qmet050080W/m3這里有兩個(gè)需要特別說(shuō)明的地方。第一個(gè)是表皮層一般設(shè)定為無(wú)血液灌注、無(wú)代謝產(chǎn)熱因?yàn)楸砥ぶ袥](méi)有血管分布這也是為什么表皮溫度對(duì)環(huán)境變化的響應(yīng)最靈敏、也最容易被燙傷。第二個(gè)是單位問(wèn)題——很多文獻(xiàn)給出的灌注率是以mL/(100g·min)為單位的要換算成SI單位制里的每秒s?1需要進(jìn)行兩步換算先把質(zhì)量換算成體積需要乘以組織密度再把分鐘換算成秒。我見(jiàn)過(guò)不止一個(gè)組因?yàn)橹苯幽梦墨I(xiàn)數(shù)值填入仿真軟件結(jié)果血液灌注項(xiàng)被放大了幾百倍整個(gè)溫度場(chǎng)看起來(lái)完全不合常理。血液的熱物性參數(shù)相對(duì)固定ρb 1060 kg/m3cb 3650 J/(kg·K)動(dòng)脈血溫度取體內(nèi)的36.7℃。動(dòng)脈血溫和核心體溫有細(xì)微差別但在皮膚淺層溫度場(chǎng)模型中直接取37℃或者36.7℃都行對(duì)結(jié)果影響很小。代謝產(chǎn)熱在靜息狀態(tài)下對(duì)深層組織溫度分布有一定影響但對(duì)皮膚表面溫度影響不大——不過(guò)如果你后面要擴(kuò)展到運(yùn)動(dòng)、炎癥組織或腫瘤熱療場(chǎng)景這兩個(gè)參數(shù)的取值就要認(rèn)真對(duì)待了。2.3 穩(wěn)態(tài)與瞬態(tài)模型的選擇依據(jù)很多人一開(kāi)始就會(huì)問(wèn)我應(yīng)該做穩(wěn)態(tài)還是瞬態(tài)答案取決于你的實(shí)際場(chǎng)景。如果只關(guān)心皮膚在恒定環(huán)境下足夠長(zhǎng)時(shí)間后的最終溫度分布——比如一個(gè)恒定功率的暖寶寶貼了半小時(shí)以上——穩(wěn)態(tài)模型就夠了計(jì)算量小調(diào)試方便。但如果你關(guān)心的是加熱前幾分鐘內(nèi)皮膚溫度的上升過(guò)程或者要評(píng)估多少秒內(nèi)會(huì)達(dá)到燙傷溫度這種問(wèn)題就必須用瞬態(tài)模型。我的經(jīng)驗(yàn)是先把穩(wěn)態(tài)模型跑通得到合理結(jié)果后再切換到瞬態(tài)研究。這樣做的好處在于穩(wěn)態(tài)模型中不存在時(shí)間步長(zhǎng)和初始條件的問(wèn)題一旦結(jié)果不對(duì)很可能是邊界條件或材料參數(shù)的問(wèn)題排查范圍小很多。穩(wěn)態(tài)驗(yàn)證通過(guò)后再加瞬態(tài)此時(shí)主要處理的就是初始條件和時(shí)間步長(zhǎng)了。這就好比先學(xué)會(huì)走再學(xué)跑避免一開(kāi)始就被一堆耦合因素繞暈。3. 幾何建模與材料賦值的實(shí)操記錄方程確定之后就要把皮膚結(jié)構(gòu)變成COMSOL里的幾何模型。皮膚在微觀上是很復(fù)雜的多層結(jié)構(gòu)但從傳熱角度看很多微觀細(xì)節(jié)比如毛囊、皮脂腺、毛細(xì)血管對(duì)宏觀溫度場(chǎng)影響很小可以忽略。我用的是一個(gè)帶軸對(duì)性的二維模型在足夠近似的前提下大幅降低計(jì)算量這是生物傳熱建模最常用的方法。3.1 皮膚分層結(jié)構(gòu)的幾何尺寸設(shè)定我要建的幾何是1cm × 1cm的皮膚截面從上到下分三層表皮厚度0.1mm、真皮厚度1.5mm、皮下組織厚度3mm。為什么這么取值這是有講究的。表皮在手掌和腳掌處可以到0.4mm以上但身體大部分區(qū)域只有0.05~0.1mm取0.1mm是一個(gè)典型均值。真皮一般在1~2mm之間取1.5mm。皮下組織的厚度個(gè)體差異最大身體不同部位可以從幾毫米到幾厘米不等我這里取3mm——原因在于我們關(guān)心的是皮膚淺層溫度分布皮下組織取得足夠厚就可以讓深部邊界的影響減弱到可接受范圍。COMSOL的操作并不復(fù)雜在模型向?qū)е羞x擇二維軸對(duì)稱(chēng)空間維度添加固體傳熱物理場(chǎng)接口然后從幾何節(jié)點(diǎn)創(chuàng)建三個(gè)矩形Rectangle并疊放起來(lái)對(duì)齊到同一坐標(biāo)軸。一個(gè)常見(jiàn)的錯(cuò)誤是把三個(gè)層建成分離的幾何體而忘記做并集操作導(dǎo)致層與層之間出現(xiàn)接觸熱阻溫度場(chǎng)會(huì)在界面處出現(xiàn)不自然的突變。正確做法是在形成聯(lián)合體Form Union時(shí)確保所有層被合并為一個(gè)域然后在材料節(jié)點(diǎn)按域分別賦予不同層的材料屬性。網(wǎng)格劃分的麻煩主要在表皮這一層。0.1mm的厚度相對(duì)于整個(gè)模型3~5mm的尺度很容易導(dǎo)致網(wǎng)格數(shù)量爆炸或長(zhǎng)寬比嚴(yán)重失調(diào)。如果你直接使用默認(rèn)的極細(xì)化網(wǎng)格COMSOL會(huì)在表皮層里塞進(jìn)幾萬(wàn)個(gè)極薄的三角形單元求解速度會(huì)非常慢。我的做法是對(duì)表皮層單獨(dú)使用邊界層網(wǎng)格Boundary Layer功能在表皮內(nèi)生成2~3層極度扁平的邊界層單元而在真皮和皮下組織中使用自由三角形網(wǎng)格。這樣既能保證薄層內(nèi)有足夠的節(jié)點(diǎn)分辨溫度梯度又能控制整體單元數(shù)量在可接受的范圍。3.2 材料參數(shù)在COMSOL中的配置技巧幾何畫(huà)好后就要給每一層指定材料參數(shù)。COMSOL的材料庫(kù)中本身沒(méi)有現(xiàn)成的皮膚材料需要自定義。我在材料節(jié)點(diǎn)下新建三個(gè)空材料分別命名Epidermis、Dermis、Subcutaneous然后在材料屬性中填入前述的k、ρ、c。這里要特別提醒COMSOL的默認(rèn)單位制是SI制導(dǎo)熱系數(shù)填0.235就是0.235 W/(m·K)密度填1200就是1200 kg/m3但血液灌注率這個(gè)參數(shù)沒(méi)有現(xiàn)成的屬性名你需要在固體傳熱物理場(chǎng)中添加一個(gè)熱源節(jié)點(diǎn)將熱源表達(dá)式寫(xiě)成0.0012 × 3650 × 1060 × (37[degC] - T)這就是把血液灌注項(xiàng)直接以熱源形式施加到真皮和皮下組織的域上。仔細(xì)看這個(gè)式子0.0012是灌注率3650是血液比熱容1060是血液密度三者相乘再乘以動(dòng)脈血溫度-局部溫度恰好就是Pennes方程中的血液灌注項(xiàng)。用這種方法等于把方程重新實(shí)現(xiàn)了一遍但好處是你不需要修改COMSOL的底層方程所有操作都可以在圖形界面里完成。有一個(gè)微妙但重要的細(xì)節(jié)COMSOL中的溫度變量T默認(rèn)是絕對(duì)溫度單位的物理量但當(dāng)你在表達(dá)式里寫(xiě)37[degC]時(shí)軟件會(huì)正確地將它表示為310.15K。不過(guò)表達(dá)式中的(T)也必須和它保持一致所以我建議用degC和K混寫(xiě)時(shí)要打心底里清楚自己在干什么。如果直接在表達(dá)式里寫(xiě)37-T那相當(dāng)于計(jì)算了37K與絕對(duì)溫度的差值結(jié)果會(huì)完全錯(cuò)誤。3.3 材料參數(shù)獲取的靠譜渠道很多讀者會(huì)問(wèn)你這些參數(shù)是哪里查的我個(gè)人的優(yōu)先級(jí)是首先查近五年的綜述文章特別是關(guān)于生物組織熱物性的綜述它們通常會(huì)匯總大量測(cè)量數(shù)據(jù)并給出推薦值其次查經(jīng)典教材比如生物傳熱學(xué)相關(guān)的標(biāo)準(zhǔn)參考書(shū)書(shū)中對(duì)皮膚各層參數(shù)給出了非常穩(wěn)定的數(shù)值最后才是單篇研究論文但單篇論文的參數(shù)往往針對(duì)特定部位或特定物種需要仔細(xì)甄別。需要特別注意的是人的個(gè)體差異。年齡、性別、身體部位、BMI都會(huì)導(dǎo)致皮膚厚度和熱物性出現(xiàn)明顯變化。比如老年人的表皮更薄、真皮層血液灌注下降同樣的熱刺激條件下他們的皮膚溫度會(huì)更高、燙傷風(fēng)險(xiǎn)也更大——這也是臨床上老年人對(duì)溫度不敏感卻更容易燙傷的原因之一。如果在模型里只使用一組平均人參數(shù)得到的結(jié)果對(duì)個(gè)體而言可能不完全準(zhǔn)確但是趨勢(shì)判斷和機(jī)制分析是足夠用的。4. 邊界條件、熱源與初始條件的實(shí)操配置模型的核心求解域是三層皮膚邊界條件決定了這個(gè)域與環(huán)境、體內(nèi)以及其他組織之間的熱交換方式。這一步做得不好哪怕前面參數(shù)再準(zhǔn)確溫度場(chǎng)也會(huì)嚴(yán)重失真。我把自己常用的邊界條件配置寫(xiě)下來(lái)每一項(xiàng)背后都說(shuō)清楚物理依據(jù)。4.1 皮膚表面與環(huán)境的對(duì)流換熱配置皮膚表面是最重要的邊界。它與空氣之間存在自然對(duì)流換熱這一項(xiàng)的數(shù)學(xué)表達(dá)是q h · (T_surface - T_ambient)其中h是對(duì)流換熱系數(shù)T_surface是皮膚表面溫度T_ambient是環(huán)境溫度。h的取值對(duì)結(jié)果十分敏感但恰恰是初學(xué)者最容易拍腦袋填的地方。靜止空氣中的自然對(duì)流系數(shù)一般在5~10 W/(m2·K)之間我通常取8有風(fēng)情況下這個(gè)值會(huì)大到20以上如果皮膚直接接觸固態(tài)熱源如暖寶寶、熱水袋則接觸傳熱的等效換熱系數(shù)可以達(dá)到幾百甚至不能單純用對(duì)流系數(shù)來(lái)描述。這里要給大家一個(gè)非常重要的實(shí)操建議在COMSOL中不要把對(duì)流邊界條件寫(xiě)成固定的熱流或固定溫度而是用對(duì)流熱通量邊界節(jié)點(diǎn)把h和環(huán)境溫度設(shè)置成參數(shù)。這樣做的直接好處是后面對(duì)h做參數(shù)掃描時(shí)完全不需要改動(dòng)邊界條件結(jié)構(gòu)。我第一版模型就是把h寫(xiě)成8然后直接硬編碼在邊界里后來(lái)要對(duì)比不同風(fēng)速下溫度場(chǎng)變化時(shí)不得不逐個(gè)手動(dòng)修改浪費(fèi)了大量時(shí)間。參數(shù)化配置是COMSOL使用中花幾分鐘省幾小時(shí)的典型操作。4.2 深部邊界與側(cè)邊界的處理深層邊界是另一個(gè)容易模糊的地方。人體核心區(qū)溫度基本恒定為37℃左右所以皮下組織最底層可以直接設(shè)置成一個(gè)固定溫度邊界即Dirichlet邊界條件溫度值為37℃。這背后的生理依據(jù)是當(dāng)皮下組織厚度大于幾個(gè)毫米時(shí)組織深處的溫度已經(jīng)非常接近核心體溫與外界環(huán)境的關(guān)系不大。這個(gè)邊界條件的設(shè)置方法是在固體傳熱物理場(chǎng)中添加溫度邊界節(jié)點(diǎn)選中皮下組織最底部的邊界設(shè)置T0為37℃即可。側(cè)邊界的處理有兩種思路。如果你建的是軸對(duì)稱(chēng)模型——比如把皮膚看成圓柱體——側(cè)邊默認(rèn)是軸對(duì)稱(chēng)邊界不需要額外設(shè)置。如果建的是二維平面模型側(cè)邊界應(yīng)當(dāng)設(shè)置為絕緣即零熱通量因?yàn)檫@個(gè)邊界代表的是皮膚中一個(gè)假想的對(duì)稱(chēng)截面兩側(cè)的溫度梯度應(yīng)該是一致的。這種絕緣邊界預(yù)設(shè)了側(cè)邊沒(méi)有熱量橫向流動(dòng)對(duì)于模擬一個(gè)均勻熱源貼在一大片皮膚中心區(qū)域的情況是合理的。但如果你的熱源很小且緊貼邊緣側(cè)邊界的設(shè)置就要改成對(duì)流或熱通量邊界否則熱量會(huì)被人為困在模型內(nèi)導(dǎo)致局部溫度偏高。我踩過(guò)這個(gè)坑——最初為了省計(jì)算量把模型尺寸縮小到3mm×3mm又沒(méi)有設(shè)置好側(cè)邊界結(jié)果發(fā)現(xiàn)同樣的邊界條件下溫度比實(shí)測(cè)高出近2℃后來(lái)把計(jì)算域擴(kuò)到10mm×10mm并設(shè)置絕緣邊界結(jié)果才和實(shí)測(cè)吻合。4.3 瞬態(tài)模擬的初始條件設(shè)置與時(shí)間尺度如果你要跑瞬態(tài)模擬初始條件是相對(duì)容易忽略的變量。這里需要想清楚一個(gè)物理問(wèn)題皮膚在未受外界熱源干擾時(shí)自身已經(jīng)有一個(gè)穩(wěn)定的溫度分布——從表面約32~33℃逐漸過(guò)渡到體內(nèi)37℃。這個(gè)初始溫度分布就是模型模擬皮膚從室溫狀態(tài)開(kāi)始被加熱過(guò)程的起點(diǎn)。獲取這個(gè)初始分布最可靠的辦法是先跑一個(gè)沒(méi)有加熱源的穩(wěn)態(tài)模擬得到皮膚的自然溫度梯度然后把這個(gè)穩(wěn)態(tài)解作為瞬態(tài)模擬的初始條件。在COMSOL里你可以在研究中先添加穩(wěn)態(tài)研究步驟求解一次再添加瞬態(tài)研究步驟并在瞬態(tài)步驟的設(shè)置中指定初始值不來(lái)自物理場(chǎng)默認(rèn)值而是來(lái)自上一個(gè)穩(wěn)態(tài)研究的解。這個(gè)做法在工程上叫做逐步求解或預(yù)加載初始條件它比直接給一個(gè)均勻的37℃初始值要精確得多尤其對(duì)表面溫度隨時(shí)間變化的早期階段影響顯著。時(shí)間尺度的把握也值得多說(shuō)一步。皮膚熱響應(yīng)有快有慢表皮層很薄熱響應(yīng)時(shí)間常數(shù)只有幾秒到幾十秒皮下組織厚時(shí)間常數(shù)可能到幾分鐘甚至更長(zhǎng)。如果時(shí)間步長(zhǎng)設(shè)得太大短時(shí)間的快速升溫過(guò)程會(huì)被嚴(yán)重平滑結(jié)果失真。我的經(jīng)驗(yàn)是如果模擬時(shí)長(zhǎng)為30分鐘左右設(shè)置步長(zhǎng)為1~5秒就足夠精確如果模擬幾個(gè)小時(shí)步長(zhǎng)可以放寬到10~30秒。求解器本身的自適應(yīng)時(shí)間步進(jìn)功能也可以打開(kāi)但初始的最小步長(zhǎng)還是建議手動(dòng)控制一下防止早期振蕩。5. 求解器設(shè)置、網(wǎng)格無(wú)關(guān)性與收斂性調(diào)試很多人在這一步會(huì)被求解器的各種選項(xiàng)勸退。實(shí)際上對(duì)于這種尺度的導(dǎo)熱模型COMSOL默認(rèn)的求解器設(shè)置已經(jīng)能應(yīng)付大部分場(chǎng)景真正需要手動(dòng)干預(yù)的是網(wǎng)格質(zhì)量和瞬態(tài)求解的時(shí)間步控制。以下是幾個(gè)可以完全照抄的經(jīng)驗(yàn)值配合一段時(shí)間對(duì)溫度的殘差監(jiān)控使用效果更佳。5.1 網(wǎng)格劃分策略與無(wú)關(guān)性驗(yàn)證網(wǎng)格無(wú)關(guān)性驗(yàn)證是仿真報(bào)告的必備環(huán)節(jié)。基本做法是至少劃分三套網(wǎng)格粗、中等、細(xì)觀察關(guān)注點(diǎn)比如皮膚表面中心溫度的數(shù)值變化。如果粗網(wǎng)格和細(xì)網(wǎng)格的同一物理量差值小于0.1℃就可以認(rèn)為結(jié)果已經(jīng)與網(wǎng)格無(wú)關(guān)。這里必須強(qiáng)調(diào)一個(gè)容易被忽略的原則網(wǎng)格無(wú)關(guān)性驗(yàn)證不是整模型一起加密而是要對(duì)溫差大的區(qū)域——尤其是表皮——單獨(dú)加密。只加密真皮而下層的網(wǎng)格表面溫度幾乎沒(méi)變化但如果你把表皮從兩層邊界層改成四層邊界層表面溫度的變化可能達(dá)到0.3℃甚至更多。我在實(shí)際建模中使用了COMSOL的邊界層功能做了三次加密測(cè)試。第一次表皮層只有1層邊界層計(jì)算很快但表面溫度比實(shí)驗(yàn)值低了0.8℃第二次改成了3層誤差縮到0.15℃以?xún)?nèi)第三次改成5層誤差幾乎不再變化。這個(gè)結(jié)果說(shuō)明3層邊界層已經(jīng)足夠。盲目加密到5層不會(huì)顯著提高精度只會(huì)讓求解時(shí)間從幾秒增加到幾十秒——不必為消除不了的那0.05℃誤差付出過(guò)大代價(jià)。5.2 穩(wěn)態(tài)與瞬態(tài)求解器的參數(shù)推薦對(duì)于穩(wěn)態(tài)問(wèn)題COMSOL默認(rèn)的穩(wěn)態(tài)求解器通常會(huì)自動(dòng)選擇直接求解器PARDISO或MUMPS就可以直接跑不需要額外改動(dòng)。有時(shí)候PARDISO報(bào)錯(cuò)內(nèi)存不足切換到MUMPS往往能解決問(wèn)題——兩者的數(shù)值穩(wěn)定性略有差異在對(duì)稱(chēng)問(wèn)題上MUMPS通常更穩(wěn)但內(nèi)存占用也稍高。如果求解仍然發(fā)散第一步檢查網(wǎng)格質(zhì)量第二步檢查材料參數(shù)是否出現(xiàn)負(fù)值或跨數(shù)量級(jí)的突變。三層皮膚的熱導(dǎo)率差異不算大但如果有人把表皮導(dǎo)熱系數(shù)從0.235填成了0.0235溫度場(chǎng)會(huì)在表皮內(nèi)出現(xiàn)不合理的巨大梯度。瞬態(tài)求解設(shè)置中有兩個(gè)關(guān)鍵參數(shù)初始步長(zhǎng)和最大步長(zhǎng)。COMSOL使用自適應(yīng)時(shí)間步進(jìn)時(shí)默認(rèn)值通?？梢詽M(mǎn)足一般需要但初始步長(zhǎng)我建議手動(dòng)指定不要依賴(lài)默認(rèn)的自由。因?yàn)槠つw表面在加熱瞬間會(huì)有一個(gè)很陡的溫度梯度變化如果初始步長(zhǎng)過(guò)大第一個(gè)時(shí)間步就會(huì)跳過(guò)關(guān)鍵的溫升過(guò)程。我的經(jīng)驗(yàn)值是初始步長(zhǎng)設(shè)為總模擬時(shí)長(zhǎng)的1/1000~1/500。比如模擬1800秒初始步長(zhǎng)設(shè)為1秒同時(shí)允許求解器在平滑階段自動(dòng)放大到20秒。5.3 調(diào)試現(xiàn)場(chǎng)一次溫度異常高的排查過(guò)程分享一個(gè)真實(shí)的調(diào)試案例。某次我把模型參數(shù)全部設(shè)置好穩(wěn)態(tài)求解后皮膚表面溫度居然高達(dá)42℃在常溫26℃環(huán)境下明顯不合理。第一反應(yīng)是檢查對(duì)流換熱系數(shù)——我以為自己填了8 W/(m2·K)結(jié)果發(fā)現(xiàn)界面上填的是8 W/(m2·K)沒(méi)錯(cuò)但環(huán)境溫度不小心設(shè)置成了45℃。所以第一輪就排除了邊界條件的問(wèn)題。第二懷疑是血液灌注熱源表達(dá)式寫(xiě)錯(cuò)。檢查后發(fā)現(xiàn)真皮和皮下組織的域選擇有誤皮下組織中填的是和真皮一樣的灌注率。皮下組織的灌注率通常比真皮低40%左右這個(gè)錯(cuò)誤會(huì)讓深部組織冷卻能力更強(qiáng)反而會(huì)降低表面溫度不會(huì)升高。這個(gè)因素影響方向不對(duì)繼續(xù)排查。第三輪才找到真正的問(wèn)題我設(shè)置的加熱熱源的功率密度是5000 W/m2但是加熱源的邊界被加在了皮膚表面的一個(gè)0.5mm寬的狹長(zhǎng)區(qū)域上而總面積只有1cm×1cm這樣一來(lái)等效熱流密度其實(shí)被放大了20倍——因?yàn)闊崃枯斎胧前疵娣e積分的同樣總量加在更小的面積上局部的熱流密度就會(huì)超高。把熱源的施加面積和功率密度統(tǒng)一為合理的值后表面溫度落到了34.5℃附近和預(yù)實(shí)驗(yàn)的紅外測(cè)溫?cái)?shù)據(jù)基本吻合。這個(gè)案例的教訓(xùn)是拿到異常結(jié)果時(shí)先從邊界條件找原因再?gòu)牟牧蠀?shù)找原因最后才懷疑求解器本身。6. 結(jié)果后處理與模型擴(kuò)展場(chǎng)景求解完成后后處理其實(shí)比大多數(shù)人想象得更重要。COMSOL自帶的繪圖和數(shù)據(jù)集提取功能已經(jīng)相當(dāng)強(qiáng)大關(guān)鍵是你要清楚自己需要什么樣的輸出。是定性看溫度云圖還是定量提取某條深度線(xiàn)上的溫度曲線(xiàn)又或者是輸出表面溫度隨時(shí)間變化的瞬態(tài)響應(yīng)曲線(xiàn)——不同的需求對(duì)應(yīng)不同的數(shù)據(jù)提取操作。6.1 溫度場(chǎng)可視化與關(guān)鍵數(shù)據(jù)提取對(duì)于二維軸對(duì)稱(chēng)模型最直接的結(jié)果展示是溫度分布云圖。用COMSOL的二維繪圖組功能選擇表面圖變量選T就可以得到彩色的溫度場(chǎng)分布。由于皮膚層很薄默認(rèn)的色彩映射下三層組織往往都擠在一起看不清細(xì)節(jié)。這里我建議選中高度表達(dá)式選項(xiàng)把幾何的Y軸按比例拉伸比如把表皮厚度從0.1mm在繪圖時(shí)顯示為1mm層間差異會(huì)一目了然。這個(gè)操作只改變顯示效果不修改幾何數(shù)據(jù)。定量提取數(shù)據(jù)的核心操作是繪制截線(xiàn)圖。在二維結(jié)果中你可以通過(guò)一維繪圖組-截線(xiàn)功能畫(huà)一條從體表到深部的垂直線(xiàn)然后輸出溫度沿深度的變化曲線(xiàn)。這條深度-溫度曲線(xiàn)是驗(yàn)證模型準(zhǔn)確性的最重要依據(jù)在正常室溫下它應(yīng)該呈現(xiàn)體表溫度約33℃、深層逐漸趨向37℃的單調(diào)遞增形態(tài)。如果曲線(xiàn)中出現(xiàn)了非單調(diào)的拐點(diǎn)或突變通常意味著材料參數(shù)或邊界條件的設(shè)置有誤。瞬態(tài)結(jié)果的提取則可以使用一維繪圖組-點(diǎn)圖功能選中體表上的一個(gè)點(diǎn)就可以輸出該點(diǎn)溫度隨時(shí)間的變化。比如加熱前30分鐘的表面溫度曲線(xiàn)可以用來(lái)評(píng)估燙傷風(fēng)險(xiǎn)——結(jié)合已知的43℃以上每上升1℃安全暴露時(shí)間減半的工程經(jīng)驗(yàn)?zāi)苤苯优袛喈a(chǎn)品的安全范圍。6.2 模型復(fù)用與多物理場(chǎng)擴(kuò)展思路基礎(chǔ)模型跑通之后擴(kuò)展空間非常大。我這里列幾個(gè)常見(jiàn)的擴(kuò)展方向方便你自己后續(xù)往下走。第一加上衣物或保溫層。很多實(shí)際問(wèn)題中皮膚表面并不是直接暴露在空氣里而是隔著一層布、一個(gè)加熱片、一層凝膠墊。在COMSOL中這些都可以作為額外的幾何層或殼層加入模型只需定義它們的導(dǎo)熱系數(shù)和厚度即可。比如評(píng)估暖手寶的熱性能時(shí)就可以在皮膚模型上疊加一層面料層再設(shè)定面料外表面的對(duì)流邊界條件。第二引入多物理場(chǎng)耦合。皮膚表面出汗時(shí)會(huì)蒸發(fā)散熱這個(gè)散熱通量取決于濕度差和蒸發(fā)潛熱。你可以用COMSOL的濕空氣傳熱模塊或手動(dòng)添加蒸發(fā)冷卻項(xiàng)從而評(píng)估運(yùn)動(dòng)狀態(tài)下人體的熱舒適性。此外血液灌注率會(huì)隨溫度變化而改變這也可以在熱源表達(dá)式中設(shè)置成與溫度相關(guān)的函數(shù)實(shí)現(xiàn)更真實(shí)的生理反饋。第三與熱損傷模型耦合。如果你有工程背景想進(jìn)一步做低溫燙傷或熱療劑量評(píng)估可以引入Arrhenius損傷積分計(jì)算組織在不同溫度-時(shí)間組合下的損傷累積程度。COMSOL中可以在固體傳熱之外再添加一個(gè)因變量方程把損傷度作為未知量與溫度場(chǎng)同步求解。這類(lèi)模型在醫(yī)學(xué)物理研究中很常見(jiàn)屬于比較經(jīng)典的應(yīng)用方向。6.3 參數(shù)掃描與結(jié)果全局分析COMSOL的參數(shù)掃描功能可以自動(dòng)化地完成修改參數(shù) - 求解 - 記錄結(jié)果的流程這在模型調(diào)試階段是極其高效的。以皮膚模型為例可以分別對(duì)加熱功率、環(huán)境溫度、對(duì)流換熱系數(shù)、血液灌注率這四個(gè)參數(shù)做掃描得到一系列溫度場(chǎng)結(jié)果然后匯總成不同條件下的表面溫度變化表。我自己的一個(gè)實(shí)際掃描案例是這樣把加熱片的表面溫度從40℃以1℃為步長(zhǎng)遞增到50℃對(duì)每個(gè)溫度都求解穩(wěn)態(tài)溫度場(chǎng)和1小時(shí)內(nèi)的瞬態(tài)變化。最終得到表面溫度-內(nèi)部峰值溫度-達(dá)到峰值的時(shí)間三張關(guān)系表直接交給產(chǎn)品團(tuán)隊(duì)評(píng)估設(shè)定多少度比較安全。如果沒(méi)有參數(shù)掃描功能這個(gè)工作量需要手動(dòng)重復(fù)11次建模求解非常痛苦。關(guān)于掃描結(jié)果的數(shù)據(jù)導(dǎo)出COMSOL通常可以生成一個(gè)表格你也可以用外部腳本比如通過(guò)編程導(dǎo)出結(jié)構(gòu)化數(shù)據(jù)深入分析。這種批量化思維不僅在COMSOL里有價(jià)值在任何仿真工作流里都是提升效率的關(guān)鍵。7. 常見(jiàn)問(wèn)題速查與個(gè)人經(jīng)驗(yàn)總結(jié)模型做完后我把這段時(shí)間積累的典型問(wèn)題和排查經(jīng)驗(yàn)整理成了一個(gè)速查表方便后續(xù)復(fù)用。這些問(wèn)題分布在整個(gè)建模鏈路的不同環(huán)節(jié)但只要掌握了排查邏輯成功率會(huì)有明顯提升。故障現(xiàn)象最可能原因排查順序與解法表面溫度偏高 / 溫度云圖整體過(guò)熱環(huán)境溫度設(shè)置錯(cuò)誤熱源施加面積過(guò)小先檢查邊界條件數(shù)值再檢查熱源邊界選擇表面溫度偏低對(duì)流換熱系數(shù)過(guò)大表皮層網(wǎng)格太少調(diào)低h增加表皮邊界層數(shù)層界面溫度突變不連續(xù)未做幾何并集層間材料未正確分配檢查形成聯(lián)合體并重新賦材料瞬態(tài)求解發(fā)散初始步長(zhǎng)過(guò)大初始條件不合理縮短初始步長(zhǎng)使用穩(wěn)態(tài)解作初始值血液灌注項(xiàng)無(wú)效果表達(dá)式未乘以密度/比熱容域選擇錯(cuò)誤核對(duì)完整表達(dá)式確認(rèn)域選定求解器內(nèi)存不足網(wǎng)格過(guò)細(xì)求解器選擇不當(dāng)換MUMPS求解器適度降低網(wǎng)格細(xì)度模型調(diào)試完成后我通常還會(huì)有最后一道驗(yàn)證工序把仿真得到的皮膚表面溫度與文獻(xiàn)中實(shí)測(cè)的紅外熱像圖數(shù)據(jù)對(duì)比。一般室溫26℃無(wú)風(fēng)環(huán)境下人體前臂皮膚表面溫度在32~34℃之間如果仿真結(jié)果落在這個(gè)范圍內(nèi)模型就基本可靠。這個(gè)驗(yàn)證錨點(diǎn)不需要復(fù)雜實(shí)驗(yàn)但非常有效。最后再分享一個(gè)我個(gè)人的小經(jīng)驗(yàn)在整個(gè)皮膚溫度場(chǎng)模型的開(kāi)發(fā)周期里花費(fèi)時(shí)間最長(zhǎng)的往往不是求解本身而是參數(shù)反復(fù)調(diào)整和結(jié)果驗(yàn)證。前期花時(shí)間把材料參數(shù)的單位、邊界條件的物理含義徹底搞清楚后面反而最省時(shí)間。仿真工具的本質(zhì)是幫我們快速驗(yàn)證物理假設(shè)而不是替我們決策——做模型之前先把問(wèn)題問(wèn)清楚會(huì)讓整個(gè)建模過(guò)程順暢很多。