算)
1. 為什么會(huì)有這篇番外最近在啃MFiX的源代碼主線的顆粒流體力學(xué)的部分還沒啃完中途遇到一個(gè)不算起眼的小程序文件calc_thermo_des.f本來打算掃一眼就跳過去結(jié)果越看越覺得有意思。它不涉及復(fù)雜的算法推導(dǎo)也不涉及高深的數(shù)值格式但它是顆粒相能量計(jì)算里繞不開的一個(gè)環(huán)節(jié)。更重要的是它作為代碼里的“附帶題”把我之前零零散散看到的熱傳導(dǎo)、對(duì)流傳熱、輻射換熱這些概念全都串了起來。先說點(diǎn)題外話。讀源代碼這件事很多剛?cè)腴T的朋友容易陷入一個(gè)誤區(qū)把主程序、核心模塊啃得特別細(xì)遇到邊邊角角的子程序就直接跳。但實(shí)際模擬工程里真正讓結(jié)果對(duì)不上的往往就是這些邊邊角角的細(xì)節(jié)。比如顆粒溫度迭代發(fā)散、床層溫度場(chǎng)分布不對(duì)、熱流密度異常跳變排查到最后通常都在一些容易被忽略的輔助子程序里找到根因。我個(gè)人的習(xí)慣是主程序大框架看兩三遍即可反而要把每個(gè)被調(diào)用的子程序徹底弄清“輸入、輸出、內(nèi)部邏輯三步”因?yàn)榇a的骨架是恒定不變的血肉才是各項(xiàng)目差異的根源。這篇文章就以calc_thermo_des.f為起點(diǎn)講一講它在MFiX DES離散元法模塊里的定位、子程序的結(jié)構(gòu)邏輯、關(guān)鍵參數(shù)怎么理解以及在調(diào)試和后續(xù)修改中常見的坑。為了便于閱讀文章不會(huì)貼出完整源碼而是基于該文件的整體邏輯按行理解的方式拆解。如果你正打算閱讀MFiX的源代碼或者正在使用DEM方法做顆粒流動(dòng)復(fù)合傳熱方面的模擬這篇文章應(yīng)該能幫你節(jié)約不少時(shí)間。2. 定位這個(gè)文件在MFiX代碼體系中處于什么角色2.1 DES模塊整體結(jié)構(gòu)里的坐標(biāo)MFiX里面和離散元相關(guān)的內(nèi)容主要集中在des/目錄下。calc_thermo_des.f在DES模塊里的地位可以類比為“顆粒相能量計(jì)算的中轉(zhuǎn)站”。為什么要強(qiáng)調(diào)“中轉(zhuǎn)站”這個(gè)說法因?yàn)樗钦麄€(gè)顆粒相熱量計(jì)算中被調(diào)用的入口子程序但真正算東西的還在后面的具體函數(shù)里。我讀代碼時(shí)有個(gè)習(xí)慣拿到一個(gè)子程序先別急著看內(nèi)部實(shí)現(xiàn)而是先查調(diào)用關(guān)系。calc_thermo_des.f的調(diào)用情況比較單一通常在 DES 的時(shí)間步進(jìn)主循環(huán)中被顯式調(diào)用。也就是說在每個(gè)DEM時(shí)間步長(zhǎng)里程序都會(huì)調(diào)用一次這個(gè)子程序然后對(duì)每個(gè)顆粒進(jìn)行熱量計(jì)算。這個(gè)文件之所以命名為“calc_thermo_des”理解起來很簡(jiǎn)單calc 是計(jì)算的意思thermo 指熱力學(xué)/熱物理量des 標(biāo)記它屬于離散元模塊。組合起來它就是“計(jì)算離散元顆粒相熱物性相關(guān)量的子程序”。需要注意的一點(diǎn)這個(gè)文件在較新的MFiX版本中命名可能是calc_thermo_des_mod.f或類似帶_mod后綴的模塊化文件取決于版本。舊版本中以calc_thermo_des.f單獨(dú)存在的形式比較常見后續(xù)版本為了統(tǒng)一代碼結(jié)構(gòu)把很多子程序都模塊化了。讀老版本代碼和新版本代碼之間最大的坑就是文件名的變化但內(nèi)部邏輯基本是一脈相承的。2.2 它算的到底是什么既然叫 thermo那它計(jì)算的內(nèi)容自然和溫度、熱量相關(guān)。具體來說這個(gè)子程序主要負(fù)責(zé)的是顆粒相的能量方程組裝與求解準(zhǔn)備。DEM方法中每個(gè)顆粒都有自己獨(dú)立的溫度這個(gè)溫度受以下幾種傳熱機(jī)制影響顆粒與流體之間的對(duì)流傳熱顆粒與顆粒直接接觸時(shí)的導(dǎo)熱顆粒與壁面之間的導(dǎo)熱和對(duì)流輻射換熱化學(xué)反應(yīng)熱如涉熱反應(yīng)顆粒MFiX將熱量源項(xiàng)顯式地分門別類在計(jì)算顆粒溫度變化率DES_COND_T等變量之前先把相關(guān)的熱量量綱數(shù)組都正確組裝好。calc_thermo_des.f在這個(gè)組裝過程中承擔(dān)了“總體調(diào)度”的責(zé)任。它本身的邏輯并不復(fù)雜核心就三步對(duì)所有顆粒做初始化或清零操作判斷哪種傳熱機(jī)制需要被激活調(diào)用對(duì)應(yīng)的子程序?qū)⒂?jì)算出的熱流量分配到每個(gè)顆粒上理解了這個(gè)定位再看整個(gè)文件就會(huì)順暢很多。它不是物理模型的實(shí)現(xiàn)者而是調(diào)度員。2.3 為什么還要單獨(dú)一個(gè)“番外”很多主線的代碼閱讀文章會(huì)集中講力學(xué)的部分接觸力模型、阻尼、切向力、摩擦力以及顆粒在流場(chǎng)里的運(yùn)動(dòng)軌跡更新。熱學(xué)部分在DEM里往往不是初學(xué)者的第一站但工程應(yīng)用中幾乎躲不開。說個(gè)我實(shí)際參與過的例子。某個(gè)處理顆粒干燥的模擬項(xiàng)目剛開始只看流動(dòng)和混合不考慮顆粒溫度。后來發(fā)現(xiàn)顆粒在設(shè)備里停留時(shí)間較長(zhǎng)熱風(fēng)換熱對(duì)顆粒含水率有決定性影響就必須加上熱能計(jì)算。這個(gè)時(shí)候calc_thermo_des.f這類子程序就成了調(diào)試的主力對(duì)象因?yàn)轭w粒溫度算不對(duì)干燥速率就完全失真后面的組分輸運(yùn)全部白搭。所以這篇番外想換個(gè)視角不按照結(jié)構(gòu)力學(xué)的路線走而是單獨(dú)把熱能計(jì)算這條線拉出來捋清楚。這樣對(duì)于正在做干燥、燃燒、熱解、氣化這類主題的人會(huì)更有針對(duì)性。3. 參數(shù)逐個(gè)拆解進(jìn)入子程序之前先認(rèn)識(shí)這些變量3.1 子程序的入?yún)⑴c出參在MFiX的DES模塊中calc_thermo_des.f的入?yún)⑼ǔ０w粒數(shù)、顆粒直徑、位置坐標(biāo)、速度、溫度等基礎(chǔ)信息以及流體相的溫度、速度、體積分?jǐn)?shù)等。出參則主要體現(xiàn)在溫度變化率數(shù)組和熱流數(shù)組上它們會(huì)在之后的時(shí)間積分中被用來更新顆粒溫度。這里有個(gè)比較容易混淆的點(diǎn)MFiX中有大量的全局?jǐn)?shù)組和局部數(shù)組。calc_thermo_des.f內(nèi)部很多變量并不通過子程序參數(shù)列表傳入傳出直接在模塊變量層面使用。這一點(diǎn)和純函數(shù)式編程的思路差異很大但這也是老牌Fortran代碼的常見風(fēng)格。首次閱讀MFiX代碼的人最需要適應(yīng)的地方就在這里一個(gè)子程序的“真實(shí)輸入輸出”要比參數(shù)列表多很多必須對(duì)照模塊定義才能完整了解。為了方便理解我將該文件中涉及的核心數(shù)組做一個(gè)歸類說明。需要注意的是不同版本的變量名可能略有差異但含義是穩(wěn)定的。顆粒溫度數(shù)組存儲(chǔ)每個(gè)顆粒當(dāng)前時(shí)刻的溫度值。單位通常是K。熱量源項(xiàng)數(shù)組存儲(chǔ)每個(gè)顆粒在當(dāng)前時(shí)間步內(nèi)接收到的總熱量單位通常是W。傳熱系數(shù)相關(guān)數(shù)組存儲(chǔ)顆粒與流體之間的對(duì)流傳熱系數(shù)或者顆粒與顆粒接觸時(shí)的等效導(dǎo)熱系數(shù)單位W/(m2·K)或W/(m·K)。體積分?jǐn)?shù)數(shù)組流體網(wǎng)格中的顆粒體積分?jǐn)?shù)會(huì)影響傳熱面積和換熱系數(shù)的計(jì)算。時(shí)間步長(zhǎng)控制變量DES子循環(huán)中實(shí)際采用的時(shí)間步長(zhǎng)在熱流計(jì)算中用于后續(xù)顯式更新。這些數(shù)組的關(guān)系可以這樣理解顆粒溫度是“狀態(tài)量”傳熱系數(shù)是“傳遞能力”熱量源項(xiàng)是“結(jié)果量”時(shí)間步長(zhǎng)是“推進(jìn)尺度”。四者構(gòu)成完整的顯式溫度更新閉環(huán)。3.2 傳熱面積怎么來的顆粒與流體之間的對(duì)流傳熱工業(yè)上最常用的基礎(chǔ)公式是牛頓冷卻定律乘以傳熱面積。在DEM代碼里每個(gè)顆粒的傳熱面積并不等于球體真實(shí)表面積因?yàn)椴皇撬斜砻娑紖⑴c換熱比如顆粒被其他顆粒緊挨著遮擋的部分就不參與。calc_thermo_des.f里有一種常見的傳熱面積計(jì)算方式如果顆粒自由懸浮在流體中默認(rèn)傳熱面積等于球體外表面面積。如果顆粒與周圍顆粒形成接觸則會(huì)對(duì)傳熱面積做修正。實(shí)際代碼中這個(gè)修正可能并不在calc_thermo_des.f內(nèi)部完成可能由接觸檢測(cè)子程序算好之后傳遞給熱學(xué)計(jì)算模塊。但調(diào)度邏輯里必須存在可供判斷的接觸狀態(tài)。從代碼閱讀角度我建議先把“傳熱面積從哪來”這個(gè)問題抓住。因?yàn)楹罄m(xù)所有的熱流量計(jì)算結(jié)果都直接正比于傳熱面積面積算錯(cuò)一處顆粒溫度就整體偏移。在MFiX的DEM中單個(gè)顆粒默認(rèn)傳熱面積通常取顆粒表面積也就是 π*d_p^2。有些版本中的顆粒直徑是實(shí)際非球形等效直徑有些則直接使用顆粒直徑需要仔細(xì)看顆粒物性初始化部分的定義。3.3 努塞爾數(shù)和傳熱系數(shù)之間的橋你可能好奇對(duì)流傳熱系數(shù)本身是怎么確定的。這就涉及顆粒雷諾數(shù)和顆粒普朗特?cái)?shù)再借助關(guān)聯(lián)式計(jì)算努塞爾數(shù)。MFiX的DEM模塊中的顆粒對(duì)流傳熱大多沿用經(jīng)典的球形顆粒在流場(chǎng)中強(qiáng)迫對(duì)流傳熱關(guān)聯(lián)式。例如常見的里茲-馬歇爾型關(guān)聯(lián)式Nu 2 0.6 * Re_p^0.5 * Pr^0.33這個(gè)式子的物理含義是當(dāng)顆粒與流體的相對(duì)速度很小、顆粒雷諾數(shù)趨近于0時(shí)努塞爾數(shù)下限是2對(duì)應(yīng)靜止球體在靜止流體中的分子導(dǎo)熱極限第二項(xiàng)則代表強(qiáng)迫對(duì)流帶來的額外增強(qiáng)。從代碼角度看calc_thermo_des.f未必直接寫這個(gè)關(guān)聯(lián)式它可能調(diào)用一個(gè)獨(dú)立函數(shù)來計(jì)算顆粒努塞爾數(shù)。我讀該文件時(shí)比較關(guān)注的是以下幾點(diǎn)顆粒雷諾數(shù)使用的是顆粒直徑不是水力直徑相對(duì)速度取的是流體速度與顆粒速度的矢量差絕對(duì)值氣體在顆粒表面溫度下的物性而不是來流平均溫度下的物性在源文件中顆粒努塞爾數(shù)的計(jì)算一般發(fā)生在傳熱系數(shù)賦值之前。傳熱系數(shù) h 則按照 Nu * λ / d_p 來還原λ 是流體熱導(dǎo)率。這里有個(gè)容易被忽略的細(xì)節(jié)如果顆粒速度和流體速度完全相等相對(duì)速度為零那么 Nu 2傳熱系數(shù)退化為純導(dǎo)熱模式。這個(gè)退化極限很重要因?yàn)樗跀?shù)值上提供了穩(wěn)定性下限。4. 內(nèi)部邏輯逐步拆解calc_thermo_des.f的代碼流程4.1 第一步初始化或清零進(jìn)入子程序后首先進(jìn)行的操作通常是對(duì)所有顆粒的傳熱相關(guān)累加量做清零或初始化。這一步在很多代碼里看起來枯燥但在數(shù)值穩(wěn)定性上卻至關(guān)重要。如果不做清零上一個(gè)時(shí)間步的殘留熱流量會(huì)被帶到下一個(gè)時(shí)間步導(dǎo)致顆粒溫度持續(xù)累積誤差最終溫度漂移、發(fā)散甚至產(chǎn)生非物理的負(fù)溫度。在舊版本代碼中這一部分可能是直接的循環(huán)賦值DO NP 1, MAX_PIP Q_SOURCE(NP) ZERO DES_COND_T(NP) ZERO END DO新版本代碼中一般會(huì)引入數(shù)組整體初始化的寫法但邏輯本質(zhì)相同。我的建議是讀源代碼時(shí)先確認(rèn)清零循環(huán)覆蓋了哪些數(shù)組。如果后續(xù)新增了自定義傳熱數(shù)組卻忘了添加清零操作那么顆粒溫度會(huì)在幾十萬步時(shí)間推進(jìn)后慢慢失控。這個(gè)坑我自己就踩過排查了幾個(gè)小時(shí)最后發(fā)現(xiàn)只是初始化數(shù)組少了一行。4.2 第二步判斷是否需要進(jìn)行傳熱計(jì)算MFiX在設(shè)計(jì)上并不是每個(gè)網(wǎng)格、每個(gè)顆粒都必然執(zhí)行全量傳熱計(jì)算它會(huì)通過邏輯變量和開關(guān)來控制哪些物理過程處于激活狀態(tài)。這部分的判斷常見有以下幾類流體相計(jì)算是否開啟傳熱顆粒相溫度是否被納入求解輻射模型是否被激活化學(xué)反應(yīng)是否釋放或吸收熱量顆粒與壁面之間是否允許換熱從調(diào)度邏輯角度看這些開關(guān)好比是控制閘門決定熱量要不要流入某個(gè)通道。我之前遇到一個(gè)情況模擬結(jié)果中顆粒溫度始終不變化排查到最后發(fā)現(xiàn)是顆粒能量方程開關(guān)沒有開啟連熱量計(jì)算都沒有進(jìn)入。這類問題發(fā)生的原因通常不是因?yàn)樵创a邏輯復(fù)雜而是因?yàn)檩斎胛募械腄ES_ENERGY_EQ或者ENERGY_EQ標(biāo)志位設(shè)置不正確導(dǎo)致子程序被調(diào)用但內(nèi)部提前返回。讀源代碼前先看輸入開關(guān)往往事半功倍。4.3 第三步類型判斷與分支選擇MFiX中支持多種傳熱模型。calc_thermo_des.f內(nèi)部最核心的分支結(jié)構(gòu)是根據(jù)顆粒與顆粒的接觸狀態(tài)以及周圍流場(chǎng)狀態(tài)決定采用哪條傳熱計(jì)算路徑。這個(gè)分支結(jié)構(gòu)是我認(rèn)為該子程序中最值得細(xì)讀的部分。我先梳理幾種典型路徑。路徑A自由懸浮顆粒只和流體換熱不涉及顆粒間接觸導(dǎo)熱。這種情形下熱量?jī)H通過對(duì)流輸運(yùn)進(jìn)入顆粒。路徑B顆粒與顆粒相互接觸存在接觸導(dǎo)熱。此時(shí)需要額外考慮兩顆粒接觸點(diǎn)附近的導(dǎo)熱效應(yīng)通常會(huì)涉及接觸半徑、等效熱導(dǎo)率等參數(shù)。路徑C顆粒與壁面接觸涉及壁面導(dǎo)熱。處理方式與顆粒間接觸導(dǎo)熱相似但壁面的熱物性參數(shù)與顆粒不同。路徑D同時(shí)存在上述多種傳熱機(jī)制需要疊加計(jì)算。在實(shí)際工程模擬中一個(gè)時(shí)間步內(nèi)單個(gè)顆?？赡芡瑫r(shí)與多個(gè)顆粒接觸并同時(shí)處于流體對(duì)流之中所以路徑D才是最常見的。因此分支代碼里往往并不是互斥的 if-elseif 結(jié)構(gòu)而是并行的多個(gè)獨(dú)立 if 判斷每條判斷都在給熱量源項(xiàng)數(shù)組做累加。這里插一句經(jīng)驗(yàn)把接觸導(dǎo)熱與對(duì)流傳熱耦合起來的時(shí)候最怕的是時(shí)間步長(zhǎng)過大導(dǎo)致顯式更新不穩(wěn)定。熱擴(kuò)散的特征時(shí)間通常比顆粒碰撞的特征時(shí)間短如果時(shí)間步長(zhǎng)按力學(xué)標(biāo)準(zhǔn)取溫度更新很容易振蕩。所以后續(xù)在時(shí)間積分環(huán)節(jié)常常需要對(duì)能量方程另做時(shí)間步長(zhǎng)限制代碼中一般會(huì)有單獨(dú)的熱時(shí)間步長(zhǎng)控制。只不過這個(gè)控制不一定在calc_thermo_des.f中可能位于主DES時(shí)間步進(jìn)循環(huán)里。4.4 第四步調(diào)用底層傳熱計(jì)算函數(shù)calc_thermo_des.f自身承擔(dān)的“粗活”是調(diào)度和累加真正的“細(xì)活”在各底層函數(shù)或子程序中完成。常見的被調(diào)用對(duì)象包括顆粒對(duì)流傳熱系數(shù)計(jì)算函數(shù)顆粒間接觸導(dǎo)熱的計(jì)算子程序輻射換熱的計(jì)算子程序壁面換熱計(jì)算子程序這些底層的函數(shù)中有的可以直接獨(dú)立測(cè)試。我的建議是閱讀完calc_thermo_des.f后再單獨(dú)去跟入其中一個(gè)底層函數(shù)比如顆粒間導(dǎo)熱計(jì)算函數(shù)看看它在不同接觸半徑下的行為是否合理。曾有同行跟我交流過說他在讀這段代碼時(shí)覺得函數(shù)跳來跳去很麻煩線索容易斷。我的習(xí)慣是手繪調(diào)用樹理清楚calc_thermo_des.f一共調(diào)用了幾層每層函數(shù)的輸入是什么輸出是什么。理完這個(gè)調(diào)用樹整條顆粒相熱量計(jì)算鏈路就裝進(jìn)腦子了。之后不管遇到多復(fù)雜的傳熱bug都能快速定位到具體某一層。5. 關(guān)鍵模型細(xì)節(jié)不能跳過的物理量含義5.1 顆粒間接觸導(dǎo)熱的等效處理接觸導(dǎo)熱在DEM代碼里是一個(gè)比較微妙的問題。宏觀直覺上兩個(gè)顆粒相互接觸熱量會(huì)通過接觸點(diǎn)傳導(dǎo)。但真實(shí)顆粒表面粗糙接觸并不是理想球面與球面的點(diǎn)接觸而是包含若干微凸體的小面積接觸。DEM為了工程可行通常使用等效接觸半徑和等效熱導(dǎo)率的簡(jiǎn)化模型。MFiX在這部分的常見處理是根據(jù)兩個(gè)顆粒的接觸重疊量和楊氏模量等參數(shù)計(jì)算出接觸半徑再用顆粒材料的導(dǎo)熱系數(shù)和接觸半徑構(gòu)造接觸熱阻。從數(shù)值實(shí)現(xiàn)上看接觸導(dǎo)熱量往往正比于溫度差乘以一個(gè)與接觸幾何相關(guān)的系數(shù)。這個(gè)系數(shù)形式類似于Q_cond H_c * (T_j - T_i)其中 H_c 為接觸傳熱系數(shù)單位與熱導(dǎo)乘以特征長(zhǎng)度的形式有關(guān)。不同文獻(xiàn)給出的模型表達(dá)式差異較大所以源代碼中具體采用哪種形式只能逐個(gè)版本確認(rèn)。在閱讀這部分代碼時(shí)我建議重點(diǎn)關(guān)注接觸半徑的計(jì)算方法。有些版本使用簡(jiǎn)化的幾何關(guān)系有些版本則從法向接觸力反推接觸面積。兩者在結(jié)果上的差異不容忽視尤其在顆粒堆積比較密實(shí)的床層中顆粒間接觸點(diǎn)多接觸導(dǎo)熱路徑復(fù)雜對(duì)溫度分布的影響非常顯著。5.2 壁面導(dǎo)熱的特殊之處MFiX DEM中壁面通常被視為特殊的顆粒或半無限大固體。calc_thermo_des.f在該情形下的邏輯往往在顆粒循環(huán)中額外判斷顆粒是否與某壁面單元重疊然后單獨(dú)計(jì)算壁面導(dǎo)熱量。壁面導(dǎo)熱和顆粒間導(dǎo)熱的差異主要體現(xiàn)在熱物性參數(shù)賦值上。壁面的導(dǎo)熱系數(shù)可能是常數(shù)也可能是溫度的函數(shù)壁面溫度可能是固定壁溫也可能是給定熱流邊界條件。這些都會(huì)改變熱流表達(dá)式的形式。我自己調(diào)試時(shí)有一個(gè)教訓(xùn)壁面溫度如果在代碼中直接使用初始溫度而不是每步更新熱量累積就會(huì)失真。這類問題不易發(fā)現(xiàn)因?yàn)槌蹩礈囟葓?chǎng)變化趨勢(shì)是合理的但定量對(duì)照實(shí)驗(yàn)結(jié)果時(shí)差異明顯。所以讀源碼時(shí)務(wù)必確認(rèn)壁面溫度的更新位置。5.3 輻射換熱的簡(jiǎn)化顆粒相在高溫場(chǎng)景下的輻射換熱不能忽略比如燃燒、熱解、氣化。但DEM中顆粒數(shù)量巨大直接做逐顆粒的輻射角系數(shù)計(jì)算幾乎不現(xiàn)實(shí)所以MFiX中往往采用各種簡(jiǎn)化輻射模型。常見的簡(jiǎn)化方式包括將顆粒視為灰體引入等效輻射溫度場(chǎng)忽略顆粒間的輻射遮擋直接按經(jīng)驗(yàn)系數(shù)修正把輻射熱流作為體積源項(xiàng)在流體相中處理calc_thermo_des.f在輻射處理環(huán)節(jié)通常只是調(diào)用模型函數(shù)真正的實(shí)現(xiàn)放在對(duì)應(yīng)的輻射子程序中。我從個(gè)人使用角度看如果模擬溫度不超過700K左右輻射換熱的貢獻(xiàn)相對(duì)有限可以暫緩但如果涉及高溫燃燒輻射權(quán)重非常大必須確認(rèn)代碼中輻射源項(xiàng)的量綱和符號(hào)否則可能把凈熱流方向搞反。量綱和符號(hào)問題是源代碼閱讀里最容易被忽視、卻也最容易出大錯(cuò)的細(xì)節(jié)。熱量傳入顆粒為正傳出為負(fù)這個(gè)約定看似簡(jiǎn)單但在多個(gè)子程序多層調(diào)用后很容易在某個(gè)環(huán)節(jié)乘一個(gè)負(fù)號(hào)導(dǎo)致顆粒溫度不升反降。6. 實(shí)操過程中最常見的問題與排查建議6.1 顆粒溫度整體偏高或偏低出現(xiàn)這種情況首先檢查對(duì)流傳熱系數(shù)模型和顆粒面積的定義。一般來說整體溫度偏高說明熱量輸入偏多要么傳熱面積偏大要么熱流密度偏大。整體溫度偏低則反過來。我實(shí)踐中的排查順序是這樣的先檢查顆粒初始化溫度是否合理檢查傳熱面積是否采用了表面積而不是投影面積檢查努塞爾數(shù)關(guān)聯(lián)式中的物性取值溫度是否正確檢查時(shí)間步長(zhǎng)是否過大導(dǎo)致顯式更新的數(shù)值耗散多數(shù)情況下問題出在物性參數(shù)的溫度取值基準(zhǔn)不統(tǒng)一。比如氣體熱導(dǎo)率如果取的是入口溫度下的值而實(shí)際顆粒表面溫度遠(yuǎn)高于入口溫度那么兩者差異會(huì)導(dǎo)致計(jì)算出的傳熱系數(shù)和真實(shí)偏差很大。6.2 顆粒溫度出現(xiàn)不正常的局部尖峰局部尖峰通常和接觸導(dǎo)熱有關(guān)。當(dāng)兩個(gè)顆粒發(fā)生碰撞接觸半徑若被高估接觸導(dǎo)熱量會(huì)在一瞬間顯著增大表現(xiàn)為局部溫度尖峰。還有一個(gè)容易被忽略的原因是接觸判斷的邏輯利用了接觸重疊量而重疊量在碰撞力模型中可能因?yàn)闀r(shí)間步長(zhǎng)過大出現(xiàn)振蕩進(jìn)而導(dǎo)致接觸導(dǎo)熱熱流也振蕩。處理思路是降低DEM時(shí)間步長(zhǎng)觀察尖峰是否減弱如果減弱說明是接觸力-導(dǎo)熱耦合導(dǎo)致的數(shù)值問題。如果降低步長(zhǎng)無效則檢查接觸半徑計(jì)算中是否合理限幅。6.3 顆粒溫度比流體溫度還平穩(wěn)如果顆粒溫度基本不變而流體溫度在變化檢查顆粒與流體之間是否真的啟用了對(duì)流傳熱。還有可能是因?yàn)轭w粒熱容設(shè)置得非常大比如誤將顆粒體積熱容寫成了比熱容導(dǎo)致熱量輸入對(duì)溫度變化貢獻(xiàn)極小。顆粒熱容的量綱在MFiX中需要特別注意。有些模型輸入的是體積熱容J/m3·K有些輸入的是質(zhì)量熱容J/kg·K弄混之后溫度響應(yīng)會(huì)差好幾個(gè)數(shù)量級(jí)。我建議在讀代碼時(shí)順手把這個(gè)量綱鏈條推導(dǎo)一遍熱流量單位W時(shí)間步長(zhǎng)單位s熱量單位J溫度變化量單位K溫度變化量 熱量 / (顆粒質(zhì)量 × 顆粒比熱容)。如果發(fā)現(xiàn)計(jì)算結(jié)果和這個(gè)關(guān)系對(duì)不上那一定是代碼中定義的某個(gè)熱容或溫度相關(guān)變量理解有誤。6.4 傳熱開關(guān)全開了但就是不算熱這類問題多半是邏輯開關(guān)的命名理解有偏差。比如MFiX中某些開關(guān)是默認(rèn)關(guān)閉需要在輸入文件中顯式開啟而另一些開關(guān)默認(rèn)是開啟的。還有一部分開關(guān)之間存在從屬關(guān)系比如總體能量方程開關(guān)不開即使顆粒能量方程開關(guān)開啟也無效。排查這類問題我通常會(huì)在calc_thermo_des.f的入口處加一個(gè)臨時(shí)輸出打印當(dāng)前步數(shù)、顆粒數(shù)、開關(guān)狀態(tài)然后觀察日志。讀取Fortran代碼時(shí)在關(guān)鍵位置臨時(shí)插入WRITE(*,*)是非常常用的調(diào)試手法雖然原始項(xiàng)目代碼中不推薦留這類輸出但臨時(shí)用一下確實(shí)高效。7. 讀這個(gè)文件時(shí)值得順手做的備忘7.1 版本差異備忘MFiX發(fā)展過程中熱學(xué)計(jì)算模塊經(jīng)歷了較大的結(jié)構(gòu)重構(gòu)。老版本里calc_thermo_des.f是獨(dú)立子程序文件新版本里可能被整合進(jìn)模塊中作為subroutine保留。文件名后綴_mod不是簡(jiǎn)單改名而是從子程序轉(zhuǎn)成了模塊化封裝調(diào)用方式也發(fā)生了變化。如果你正在讀的版本中找不到calc_thermo_des.f不要意外搜索calc_thermo_des或者calc_thermo這樣的關(guān)鍵詞大概率能找到對(duì)應(yīng)的模塊。源碼閱讀工具的全局搜索在這里比手動(dòng)翻目錄高效得多。7.2 數(shù)組維度問題MFiX中的顆粒數(shù)組維度通常不是固定常數(shù)而是動(dòng)態(tài)分配。這意味著循環(huán)邊界條件往往和顆粒最大編號(hào)相關(guān)而不是和實(shí)際當(dāng)前顆粒數(shù)相關(guān)。在calc_thermo_des.f中循環(huán)邊界如果寫的是MAX_PIP那就會(huì)遍歷所有可能的顆粒槽位包括尚未激活的顆粒。未激活顆粒的物理量可能包含垃圾值因此代碼內(nèi)部幾乎都通過IF (SOLIDS ...)或者顆粒存在性判斷來保護(hù)防止對(duì)空顆粒做熱量累加。閱讀時(shí)如果發(fā)現(xiàn)循環(huán)體里面沒有存在性判斷就要警惕。這是少數(shù)會(huì)導(dǎo)致結(jié)果莫名包含NaN或無窮大的原因。7.3 并行與串行差異MFiX支持并行計(jì)算calc_thermo_des.f的邏輯在并行情況下執(zhí)行時(shí)可能會(huì)涉及顆粒分布式存儲(chǔ)即每個(gè)進(jìn)程只負(fù)責(zé)自己區(qū)域內(nèi)的顆粒。這樣一來相鄰進(jìn)程的邊界顆粒換熱需要依賴通訊而不是簡(jiǎn)單的本地循環(huán)。如果你只是在串行環(huán)境下調(diào)試則不需要過度關(guān)注并行通訊邏輯。但如果要做大規(guī)模并行模擬就必須搞清楚顆粒熱量更新時(shí)的邊界交換策略否則顆粒溫度在跨界時(shí)會(huì)出現(xiàn)跳變。關(guān)于并行版本的代碼我現(xiàn)在還沒有完全吃透后續(xù)如果繼續(xù)深入讀完MPI通訊相關(guān)模塊再另外成文分享。8. 從代碼閱讀到實(shí)際修改的經(jīng)驗(yàn)小結(jié)讀源代碼最大的收獲不是記住了某個(gè)公式是怎么寫的而是建立起一種“物理量在代碼中如何流轉(zhuǎn)”的直覺。算例出問題時(shí)能快速判斷是物性參數(shù)問題、邊界條件問題還是模型選擇問題。calc_thermo_des.f這個(gè)文件雖然簡(jiǎn)單但它的調(diào)度邏輯幾乎涵蓋了DEM熱能計(jì)算的全部入口。把它徹底弄懂之后再去看具體的傳熱模型細(xì)節(jié)就會(huì)發(fā)現(xiàn)不過是在這個(gè)框架里填入不同的計(jì)算表達(dá)式而已。下面是我在為一個(gè)流化床干燥算例修改代碼時(shí)走過的一小段真實(shí)流程供參考。第一步確認(rèn)模擬目標(biāo)是顆粒含水率變化而不是單純溫度變化。因此需要在能量計(jì)算之外再引入組分輸運(yùn)和蒸發(fā)模型但蒸發(fā)模型的源項(xiàng)必須和calc_thermo_des.f產(chǎn)生的熱量源項(xiàng)做疊加。第二步在calc_thermo_des.f中找到顆粒熱量累加的循環(huán)在循環(huán)末尾單獨(dú)加上一個(gè)蒸發(fā)潛熱的熱量扣除項(xiàng)。這個(gè)扣除項(xiàng)是負(fù)的因?yàn)檎舭l(fā)會(huì)帶走熱量。第三步調(diào)試中發(fā)現(xiàn)顆粒溫度在水分快蒸干時(shí)出現(xiàn)輕微振蕩。原因是最小含濕率判斷和蒸發(fā)源項(xiàng)切換之間的遲滯。解決辦法是添加滯回區(qū)間讓熱源項(xiàng)平滑過渡。這段經(jīng)歷想說明的是讀源代碼不是目的達(dá)成自己的物理模型目標(biāo)才是目的。但如果不讀源代碼連往哪里加蒸發(fā)源項(xiàng)都無從下手。9. 給新手的一些更實(shí)際的閱讀方法9.1 先跑通一個(gè)算例再讀代碼硬啃源代碼尤其是沒有運(yùn)行過案例的情況下很容易迷失在變量和循環(huán)里。我建議先找一個(gè)帶傳熱過程的標(biāo)準(zhǔn)案例MFiX自帶有許多測(cè)試案例跑通一次再帶著“模擬結(jié)果是怎么算出來”的問題去讀calc_thermo_des.f。這時(shí)候變量不再是紙上的符號(hào)而是和實(shí)際模擬曲線聯(lián)系起來的量。比如跑一個(gè)顆粒靜止床層的導(dǎo)熱案例觀察顆粒溫度從初值趨向壁面溫度的過程。再到代碼里找哪個(gè)變量控制了這個(gè)導(dǎo)熱行為會(huì)很直觀。9.2 打日志比看邏輯快Fortran代碼的閱讀有一個(gè)“笨辦法”在關(guān)鍵子程序入口和出口加上臨時(shí)輸出打印顆粒編號(hào)、顆粒溫度、累計(jì)熱流、接觸狀態(tài)等信息。雖然每次運(yùn)行日志量很大但對(duì)于單個(gè)子程序的邏輯驗(yàn)證這種方式遠(yuǎn)比反復(fù)推敲代碼邏輯快。我以前的調(diào)試習(xí)慣是加一個(gè)條件輸出只打印特定顆粒編號(hào)的信息比如IF (NP 100) THEN WRITE(*,*) NP100 T, DES_TEMP_NEW(NP), Q, Q_SOURCE(NP) END IF這樣就能鎖定單個(gè)顆粒觀察它在整個(gè)時(shí)間步進(jìn)中的溫度和熱流變化。這個(gè)方法對(duì)排查顆粒溫度異常特別有效。9.3 善用差異對(duì)比MFiX支持同一算例用不同模型開關(guān)跑出不同結(jié)果。直接對(duì)比這兩種結(jié)果的溫度場(chǎng)差異再回到代碼里找對(duì)應(yīng)條件分支比憑空猜測(cè)分支邏輯更直接。我之前對(duì)比過開啟接觸導(dǎo)熱和關(guān)閉接觸導(dǎo)熱兩種設(shè)置下的床層溫度分布發(fā)現(xiàn)差異集中出現(xiàn)在堆積區(qū)說明接觸導(dǎo)熱模型在堆積區(qū)域起主導(dǎo)作用。據(jù)此再反推代碼理解就快多了。10. 最后想說的源文件雖小串聯(lián)的知識(shí)面卻不小calc_thermo_des.f在MFiX整個(gè)項(xiàng)目里算是一個(gè)規(guī)模不大的子程序。但它串聯(lián)起了顆粒相運(yùn)動(dòng)、流體與顆粒對(duì)流傳熱、顆粒間接觸導(dǎo)熱、壁面換熱和輻射換熱這些模型。如果只從代碼規(guī)模上判斷它不重要就很容易錯(cuò)過DEN模塊熱能計(jì)算體系的關(guān)鍵入口。我現(xiàn)在讀到這個(gè)階段越來越覺得讀源代碼像在拆一個(gè)復(fù)雜的機(jī)械設(shè)備看起來只有一臺(tái)機(jī)器但每個(gè)零件都有自己的角色拆得越多對(duì)整個(gè)系統(tǒng)的理解越完整。calc_thermo_des.f就是其中一個(gè)不大但絕對(duì)不能說無關(guān)緊要的零件。如果你也在讀MFiX源代碼建議以后遇到這個(gè)文件時(shí)不要跳過花一兩個(gè)小時(shí)把它的調(diào)用關(guān)系和變量流轉(zhuǎn)弄清楚。之后再看任何和顆粒傳熱相關(guān)的代碼都會(huì)有一種豁然開朗的感覺。 MFiX源代碼解析：calc_thermo_des.f與顆粒相能量計(jì)算
1. 為什么會(huì)有這篇番外最近在啃MFiX的源代碼主線的顆粒流體力學(xué)的部分還沒啃完中途遇到一個(gè)不算起眼的小程序文件calc_thermo_des.f本來打算掃一眼就跳過去結(jié)果越看越覺得有意思。它不涉及復(fù)雜的算法推導(dǎo)也不涉及高深的數(shù)值格式但它是顆粒相能量計(jì)算里繞不開的一個(gè)環(huán)節(jié)。更重要的是它作為代碼里的“附帶題”把我之前零零散散看到的熱傳導(dǎo)、對(duì)流傳熱、輻射換熱這些概念全都串了起來。先說點(diǎn)題外話。讀源代碼這件事很多剛?cè)腴T的朋友容易陷入一個(gè)誤區(qū)把主程序、核心模塊啃得特別細(xì)遇到邊邊角角的子程序就直接跳。但實(shí)際模擬工程里真正讓結(jié)果對(duì)不上的往往就是這些邊邊角角的細(xì)節(jié)。比如顆粒溫度迭代發(fā)散、床層溫度場(chǎng)分布不對(duì)、熱流密度異常跳變排查到最后通常都在一些容易被忽略的輔助子程序里找到根因。我個(gè)人的習(xí)慣是主程序大框架看兩三遍即可反而要把每個(gè)被調(diào)用的子程序徹底弄清“輸入、輸出、內(nèi)部邏輯三步”因?yàn)榇a的骨架是恒定不變的血肉才是各項(xiàng)目差異的根源。這篇文章就以calc_thermo_des.f為起點(diǎn)講一講它在MFiX DES離散元法模塊里的定位、子程序的結(jié)構(gòu)邏輯、關(guān)鍵參數(shù)怎么理解以及在調(diào)試和后續(xù)修改中常見的坑。為了便于閱讀文章不會(huì)貼出完整源碼而是基于該文件的整體邏輯按行理解的方式拆解。如果你正打算閱讀MFiX的源代碼或者正在使用DEM方法做顆粒流動(dòng)復(fù)合傳熱方面的模擬這篇文章應(yīng)該能幫你節(jié)約不少時(shí)間。2. 定位這個(gè)文件在MFiX代碼體系中處于什么角色2.1 DES模塊整體結(jié)構(gòu)里的坐標(biāo)MFiX里面和離散元相關(guān)的內(nèi)容主要集中在des/目錄下。calc_thermo_des.f在DES模塊里的地位可以類比為“顆粒相能量計(jì)算的中轉(zhuǎn)站”。為什么要強(qiáng)調(diào)“中轉(zhuǎn)站”這個(gè)說法因?yàn)樗钦麄€(gè)顆粒相熱量計(jì)算中被調(diào)用的入口子程序但真正算東西的還在后面的具體函數(shù)里。我讀代碼時(shí)有個(gè)習(xí)慣拿到一個(gè)子程序先別急著看內(nèi)部實(shí)現(xiàn)而是先查調(diào)用關(guān)系。calc_thermo_des.f的調(diào)用情況比較單一通常在 DES 的時(shí)間步進(jìn)主循環(huán)中被顯式調(diào)用。也就是說在每個(gè)DEM時(shí)間步長(zhǎng)里程序都會(huì)調(diào)用一次這個(gè)子程序然后對(duì)每個(gè)顆粒進(jìn)行熱量計(jì)算。這個(gè)文件之所以命名為“calc_thermo_des”理解起來很簡(jiǎn)單calc 是計(jì)算的意思thermo 指熱力學(xué)/熱物理量des 標(biāo)記它屬于離散元模塊。組合起來它就是“計(jì)算離散元顆粒相熱物性相關(guān)量的子程序”。需要注意的一點(diǎn)這個(gè)文件在較新的MFiX版本中命名可能是calc_thermo_des_mod.f或類似帶_mod后綴的模塊化文件取決于版本。舊版本中以calc_thermo_des.f單獨(dú)存在的形式比較常見后續(xù)版本為了統(tǒng)一代碼結(jié)構(gòu)把很多子程序都模塊化了。讀老版本代碼和新版本代碼之間最大的坑就是文件名的變化但內(nèi)部邏輯基本是一脈相承的。2.2 它算的到底是什么既然叫 thermo那它計(jì)算的內(nèi)容自然和溫度、熱量相關(guān)。具體來說這個(gè)子程序主要負(fù)責(zé)的是顆粒相的能量方程組裝與求解準(zhǔn)備。DEM方法中每個(gè)顆粒都有自己獨(dú)立的溫度這個(gè)溫度受以下幾種傳熱機(jī)制影響顆粒與流體之間的對(duì)流傳熱顆粒與顆粒直接接觸時(shí)的導(dǎo)熱顆粒與壁面之間的導(dǎo)熱和對(duì)流輻射換熱化學(xué)反應(yīng)熱如涉熱反應(yīng)顆粒MFiX將熱量源項(xiàng)顯式地分門別類在計(jì)算顆粒溫度變化率DES_COND_T等變量之前先把相關(guān)的熱量量綱數(shù)組都正確組裝好。calc_thermo_des.f在這個(gè)組裝過程中承擔(dān)了“總體調(diào)度”的責(zé)任。它本身的邏輯并不復(fù)雜核心就三步對(duì)所有顆粒做初始化或清零操作判斷哪種傳熱機(jī)制需要被激活調(diào)用對(duì)應(yīng)的子程序?qū)⒂?jì)算出的熱流量分配到每個(gè)顆粒上理解了這個(gè)定位再看整個(gè)文件就會(huì)順暢很多。它不是物理模型的實(shí)現(xiàn)者而是調(diào)度員。2.3 為什么還要單獨(dú)一個(gè)“番外”很多主線的代碼閱讀文章會(huì)集中講力學(xué)的部分接觸力模型、阻尼、切向力、摩擦力以及顆粒在流場(chǎng)里的運(yùn)動(dòng)軌跡更新。熱學(xué)部分在DEM里往往不是初學(xué)者的第一站但工程應(yīng)用中幾乎躲不開。說個(gè)我實(shí)際參與過的例子。某個(gè)處理顆粒干燥的模擬項(xiàng)目剛開始只看流動(dòng)和混合不考慮顆粒溫度。后來發(fā)現(xiàn)顆粒在設(shè)備里停留時(shí)間較長(zhǎng)熱風(fēng)換熱對(duì)顆粒含水率有決定性影響就必須加上熱能計(jì)算。這個(gè)時(shí)候calc_thermo_des.f這類子程序就成了調(diào)試的主力對(duì)象因?yàn)轭w粒溫度算不對(duì)干燥速率就完全失真后面的組分輸運(yùn)全部白搭。所以這篇番外想換個(gè)視角不按照結(jié)構(gòu)力學(xué)的路線走而是單獨(dú)把熱能計(jì)算這條線拉出來捋清楚。這樣對(duì)于正在做干燥、燃燒、熱解、氣化這類主題的人會(huì)更有針對(duì)性。3. 參數(shù)逐個(gè)拆解進(jìn)入子程序之前先認(rèn)識(shí)這些變量3.1 子程序的入?yún)⑴c出參在MFiX的DES模塊中calc_thermo_des.f的入?yún)⑼ǔ０w粒數(shù)、顆粒直徑、位置坐標(biāo)、速度、溫度等基礎(chǔ)信息以及流體相的溫度、速度、體積分?jǐn)?shù)等。出參則主要體現(xiàn)在溫度變化率數(shù)組和熱流數(shù)組上它們會(huì)在之后的時(shí)間積分中被用來更新顆粒溫度。這里有個(gè)比較容易混淆的點(diǎn)MFiX中有大量的全局?jǐn)?shù)組和局部數(shù)組。calc_thermo_des.f內(nèi)部很多變量并不通過子程序參數(shù)列表傳入傳出直接在模塊變量層面使用。這一點(diǎn)和純函數(shù)式編程的思路差異很大但這也是老牌Fortran代碼的常見風(fēng)格。首次閱讀MFiX代碼的人最需要適應(yīng)的地方就在這里一個(gè)子程序的“真實(shí)輸入輸出”要比參數(shù)列表多很多必須對(duì)照模塊定義才能完整了解。為了方便理解我將該文件中涉及的核心數(shù)組做一個(gè)歸類說明。需要注意的是不同版本的變量名可能略有差異但含義是穩(wěn)定的。顆粒溫度數(shù)組存儲(chǔ)每個(gè)顆粒當(dāng)前時(shí)刻的溫度值。單位通常是K。熱量源項(xiàng)數(shù)組存儲(chǔ)每個(gè)顆粒在當(dāng)前時(shí)間步內(nèi)接收到的總熱量單位通常是W。傳熱系數(shù)相關(guān)數(shù)組存儲(chǔ)顆粒與流體之間的對(duì)流傳熱系數(shù)或者顆粒與顆粒接觸時(shí)的等效導(dǎo)熱系數(shù)單位W/(m2·K)或W/(m·K)。體積分?jǐn)?shù)數(shù)組流體網(wǎng)格中的顆粒體積分?jǐn)?shù)會(huì)影響傳熱面積和換熱系數(shù)的計(jì)算。時(shí)間步長(zhǎng)控制變量DES子循環(huán)中實(shí)際采用的時(shí)間步長(zhǎng)在熱流計(jì)算中用于后續(xù)顯式更新。這些數(shù)組的關(guān)系可以這樣理解顆粒溫度是“狀態(tài)量”傳熱系數(shù)是“傳遞能力”熱量源項(xiàng)是“結(jié)果量”時(shí)間步長(zhǎng)是“推進(jìn)尺度”。四者構(gòu)成完整的顯式溫度更新閉環(huán)。3.2 傳熱面積怎么來的顆粒與流體之間的對(duì)流傳熱工業(yè)上最常用的基礎(chǔ)公式是牛頓冷卻定律乘以傳熱面積。在DEM代碼里每個(gè)顆粒的傳熱面積并不等于球體真實(shí)表面積因?yàn)椴皇撬斜砻娑紖⑴c換熱比如顆粒被其他顆粒緊挨著遮擋的部分就不參與。calc_thermo_des.f里有一種常見的傳熱面積計(jì)算方式如果顆粒自由懸浮在流體中默認(rèn)傳熱面積等于球體外表面面積。如果顆粒與周圍顆粒形成接觸則會(huì)對(duì)傳熱面積做修正。實(shí)際代碼中這個(gè)修正可能并不在calc_thermo_des.f內(nèi)部完成可能由接觸檢測(cè)子程序算好之后傳遞給熱學(xué)計(jì)算模塊。但調(diào)度邏輯里必須存在可供判斷的接觸狀態(tài)。從代碼閱讀角度我建議先把“傳熱面積從哪來”這個(gè)問題抓住。因?yàn)楹罄m(xù)所有的熱流量計(jì)算結(jié)果都直接正比于傳熱面積面積算錯(cuò)一處顆粒溫度就整體偏移。在MFiX的DEM中單個(gè)顆粒默認(rèn)傳熱面積通常取顆粒表面積也就是 π*d_p^2。有些版本中的顆粒直徑是實(shí)際非球形等效直徑有些則直接使用顆粒直徑需要仔細(xì)看顆粒物性初始化部分的定義。3.3 努塞爾數(shù)和傳熱系數(shù)之間的橋你可能好奇對(duì)流傳熱系數(shù)本身是怎么確定的。這就涉及顆粒雷諾數(shù)和顆粒普朗特?cái)?shù)再借助關(guān)聯(lián)式計(jì)算努塞爾數(shù)。MFiX的DEM模塊中的顆粒對(duì)流傳熱大多沿用經(jīng)典的球形顆粒在流場(chǎng)中強(qiáng)迫對(duì)流傳熱關(guān)聯(lián)式。例如常見的里茲-馬歇爾型關(guān)聯(lián)式Nu 2 0.6 * Re_p^0.5 * Pr^0.33這個(gè)式子的物理含義是當(dāng)顆粒與流體的相對(duì)速度很小、顆粒雷諾數(shù)趨近于0時(shí)努塞爾數(shù)下限是2對(duì)應(yīng)靜止球體在靜止流體中的分子導(dǎo)熱極限第二項(xiàng)則代表強(qiáng)迫對(duì)流帶來的額外增強(qiáng)。從代碼角度看calc_thermo_des.f未必直接寫這個(gè)關(guān)聯(lián)式它可能調(diào)用一個(gè)獨(dú)立函數(shù)來計(jì)算顆粒努塞爾數(shù)。我讀該文件時(shí)比較關(guān)注的是以下幾點(diǎn)顆粒雷諾數(shù)使用的是顆粒直徑不是水力直徑相對(duì)速度取的是流體速度與顆粒速度的矢量差絕對(duì)值氣體在顆粒表面溫度下的物性而不是來流平均溫度下的物性在源文件中顆粒努塞爾數(shù)的計(jì)算一般發(fā)生在傳熱系數(shù)賦值之前。傳熱系數(shù) h 則按照 Nu * λ / d_p 來還原λ 是流體熱導(dǎo)率。這里有個(gè)容易被忽略的細(xì)節(jié)如果顆粒速度和流體速度完全相等相對(duì)速度為零那么 Nu 2傳熱系數(shù)退化為純導(dǎo)熱模式。這個(gè)退化極限很重要因?yàn)樗跀?shù)值上提供了穩(wěn)定性下限。4. 內(nèi)部邏輯逐步拆解calc_thermo_des.f的代碼流程4.1 第一步初始化或清零進(jìn)入子程序后首先進(jìn)行的操作通常是對(duì)所有顆粒的傳熱相關(guān)累加量做清零或初始化。這一步在很多代碼里看起來枯燥但在數(shù)值穩(wěn)定性上卻至關(guān)重要。如果不做清零上一個(gè)時(shí)間步的殘留熱流量會(huì)被帶到下一個(gè)時(shí)間步導(dǎo)致顆粒溫度持續(xù)累積誤差最終溫度漂移、發(fā)散甚至產(chǎn)生非物理的負(fù)溫度。在舊版本代碼中這一部分可能是直接的循環(huán)賦值DO NP 1, MAX_PIP Q_SOURCE(NP) ZERO DES_COND_T(NP) ZERO END DO新版本代碼中一般會(huì)引入數(shù)組整體初始化的寫法但邏輯本質(zhì)相同。我的建議是讀源代碼時(shí)先確認(rèn)清零循環(huán)覆蓋了哪些數(shù)組。如果后續(xù)新增了自定義傳熱數(shù)組卻忘了添加清零操作那么顆粒溫度會(huì)在幾十萬步時(shí)間推進(jìn)后慢慢失控。這個(gè)坑我自己就踩過排查了幾個(gè)小時(shí)最后發(fā)現(xiàn)只是初始化數(shù)組少了一行。4.2 第二步判斷是否需要進(jìn)行傳熱計(jì)算MFiX在設(shè)計(jì)上并不是每個(gè)網(wǎng)格、每個(gè)顆粒都必然執(zhí)行全量傳熱計(jì)算它會(huì)通過邏輯變量和開關(guān)來控制哪些物理過程處于激活狀態(tài)。這部分的判斷常見有以下幾類流體相計(jì)算是否開啟傳熱顆粒相溫度是否被納入求解輻射模型是否被激活化學(xué)反應(yīng)是否釋放或吸收熱量顆粒與壁面之間是否允許換熱從調(diào)度邏輯角度看這些開關(guān)好比是控制閘門決定熱量要不要流入某個(gè)通道。我之前遇到一個(gè)情況模擬結(jié)果中顆粒溫度始終不變化排查到最后發(fā)現(xiàn)是顆粒能量方程開關(guān)沒有開啟連熱量計(jì)算都沒有進(jìn)入。這類問題發(fā)生的原因通常不是因?yàn)樵创a邏輯復(fù)雜而是因?yàn)檩斎胛募械腄ES_ENERGY_EQ或者ENERGY_EQ標(biāo)志位設(shè)置不正確導(dǎo)致子程序被調(diào)用但內(nèi)部提前返回。讀源代碼前先看輸入開關(guān)往往事半功倍。4.3 第三步類型判斷與分支選擇MFiX中支持多種傳熱模型。calc_thermo_des.f內(nèi)部最核心的分支結(jié)構(gòu)是根據(jù)顆粒與顆粒的接觸狀態(tài)以及周圍流場(chǎng)狀態(tài)決定采用哪條傳熱計(jì)算路徑。這個(gè)分支結(jié)構(gòu)是我認(rèn)為該子程序中最值得細(xì)讀的部分。我先梳理幾種典型路徑。路徑A自由懸浮顆粒只和流體換熱不涉及顆粒間接觸導(dǎo)熱。這種情形下熱量?jī)H通過對(duì)流輸運(yùn)進(jìn)入顆粒。路徑B顆粒與顆粒相互接觸存在接觸導(dǎo)熱。此時(shí)需要額外考慮兩顆粒接觸點(diǎn)附近的導(dǎo)熱效應(yīng)通常會(huì)涉及接觸半徑、等效熱導(dǎo)率等參數(shù)。路徑C顆粒與壁面接觸涉及壁面導(dǎo)熱。處理方式與顆粒間接觸導(dǎo)熱相似但壁面的熱物性參數(shù)與顆粒不同。路徑D同時(shí)存在上述多種傳熱機(jī)制需要疊加計(jì)算。在實(shí)際工程模擬中一個(gè)時(shí)間步內(nèi)單個(gè)顆?？赡芡瑫r(shí)與多個(gè)顆粒接觸并同時(shí)處于流體對(duì)流之中所以路徑D才是最常見的。因此分支代碼里往往并不是互斥的 if-elseif 結(jié)構(gòu)而是并行的多個(gè)獨(dú)立 if 判斷每條判斷都在給熱量源項(xiàng)數(shù)組做累加。這里插一句經(jīng)驗(yàn)把接觸導(dǎo)熱與對(duì)流傳熱耦合起來的時(shí)候最怕的是時(shí)間步長(zhǎng)過大導(dǎo)致顯式更新不穩(wěn)定。熱擴(kuò)散的特征時(shí)間通常比顆粒碰撞的特征時(shí)間短如果時(shí)間步長(zhǎng)按力學(xué)標(biāo)準(zhǔn)取溫度更新很容易振蕩。所以后續(xù)在時(shí)間積分環(huán)節(jié)常常需要對(duì)能量方程另做時(shí)間步長(zhǎng)限制代碼中一般會(huì)有單獨(dú)的熱時(shí)間步長(zhǎng)控制。只不過這個(gè)控制不一定在calc_thermo_des.f中可能位于主DES時(shí)間步進(jìn)循環(huán)里。4.4 第四步調(diào)用底層傳熱計(jì)算函數(shù)calc_thermo_des.f自身承擔(dān)的“粗活”是調(diào)度和累加真正的“細(xì)活”在各底層函數(shù)或子程序中完成。常見的被調(diào)用對(duì)象包括顆粒對(duì)流傳熱系數(shù)計(jì)算函數(shù)顆粒間接觸導(dǎo)熱的計(jì)算子程序輻射換熱的計(jì)算子程序壁面換熱計(jì)算子程序這些底層的函數(shù)中有的可以直接獨(dú)立測(cè)試。我的建議是閱讀完calc_thermo_des.f后再單獨(dú)去跟入其中一個(gè)底層函數(shù)比如顆粒間導(dǎo)熱計(jì)算函數(shù)看看它在不同接觸半徑下的行為是否合理。曾有同行跟我交流過說他在讀這段代碼時(shí)覺得函數(shù)跳來跳去很麻煩線索容易斷。我的習(xí)慣是手繪調(diào)用樹理清楚calc_thermo_des.f一共調(diào)用了幾層每層函數(shù)的輸入是什么輸出是什么。理完這個(gè)調(diào)用樹整條顆粒相熱量計(jì)算鏈路就裝進(jìn)腦子了。之后不管遇到多復(fù)雜的傳熱bug都能快速定位到具體某一層。5. 關(guān)鍵模型細(xì)節(jié)不能跳過的物理量含義5.1 顆粒間接觸導(dǎo)熱的等效處理接觸導(dǎo)熱在DEM代碼里是一個(gè)比較微妙的問題。宏觀直覺上兩個(gè)顆粒相互接觸熱量會(huì)通過接觸點(diǎn)傳導(dǎo)。但真實(shí)顆粒表面粗糙接觸并不是理想球面與球面的點(diǎn)接觸而是包含若干微凸體的小面積接觸。DEM為了工程可行通常使用等效接觸半徑和等效熱導(dǎo)率的簡(jiǎn)化模型。MFiX在這部分的常見處理是根據(jù)兩個(gè)顆粒的接觸重疊量和楊氏模量等參數(shù)計(jì)算出接觸半徑再用顆粒材料的導(dǎo)熱系數(shù)和接觸半徑構(gòu)造接觸熱阻。從數(shù)值實(shí)現(xiàn)上看接觸導(dǎo)熱量往往正比于溫度差乘以一個(gè)與接觸幾何相關(guān)的系數(shù)。這個(gè)系數(shù)形式類似于Q_cond H_c * (T_j - T_i)其中 H_c 為接觸傳熱系數(shù)單位與熱導(dǎo)乘以特征長(zhǎng)度的形式有關(guān)。不同文獻(xiàn)給出的模型表達(dá)式差異較大所以源代碼中具體采用哪種形式只能逐個(gè)版本確認(rèn)。在閱讀這部分代碼時(shí)我建議重點(diǎn)關(guān)注接觸半徑的計(jì)算方法。有些版本使用簡(jiǎn)化的幾何關(guān)系有些版本則從法向接觸力反推接觸面積。兩者在結(jié)果上的差異不容忽視尤其在顆粒堆積比較密實(shí)的床層中顆粒間接觸點(diǎn)多接觸導(dǎo)熱路徑復(fù)雜對(duì)溫度分布的影響非常顯著。5.2 壁面導(dǎo)熱的特殊之處MFiX DEM中壁面通常被視為特殊的顆粒或半無限大固體。calc_thermo_des.f在該情形下的邏輯往往在顆粒循環(huán)中額外判斷顆粒是否與某壁面單元重疊然后單獨(dú)計(jì)算壁面導(dǎo)熱量。壁面導(dǎo)熱和顆粒間導(dǎo)熱的差異主要體現(xiàn)在熱物性參數(shù)賦值上。壁面的導(dǎo)熱系數(shù)可能是常數(shù)也可能是溫度的函數(shù)壁面溫度可能是固定壁溫也可能是給定熱流邊界條件。這些都會(huì)改變熱流表達(dá)式的形式。我自己調(diào)試時(shí)有一個(gè)教訓(xùn)壁面溫度如果在代碼中直接使用初始溫度而不是每步更新熱量累積就會(huì)失真。這類問題不易發(fā)現(xiàn)因?yàn)槌蹩礈囟葓?chǎng)變化趨勢(shì)是合理的但定量對(duì)照實(shí)驗(yàn)結(jié)果時(shí)差異明顯。所以讀源碼時(shí)務(wù)必確認(rèn)壁面溫度的更新位置。5.3 輻射換熱的簡(jiǎn)化顆粒相在高溫場(chǎng)景下的輻射換熱不能忽略比如燃燒、熱解、氣化。但DEM中顆粒數(shù)量巨大直接做逐顆粒的輻射角系數(shù)計(jì)算幾乎不現(xiàn)實(shí)所以MFiX中往往采用各種簡(jiǎn)化輻射模型。常見的簡(jiǎn)化方式包括將顆粒視為灰體引入等效輻射溫度場(chǎng)忽略顆粒間的輻射遮擋直接按經(jīng)驗(yàn)系數(shù)修正把輻射熱流作為體積源項(xiàng)在流體相中處理calc_thermo_des.f在輻射處理環(huán)節(jié)通常只是調(diào)用模型函數(shù)真正的實(shí)現(xiàn)放在對(duì)應(yīng)的輻射子程序中。我從個(gè)人使用角度看如果模擬溫度不超過700K左右輻射換熱的貢獻(xiàn)相對(duì)有限可以暫緩但如果涉及高溫燃燒輻射權(quán)重非常大必須確認(rèn)代碼中輻射源項(xiàng)的量綱和符號(hào)否則可能把凈熱流方向搞反。量綱和符號(hào)問題是源代碼閱讀里最容易被忽視、卻也最容易出大錯(cuò)的細(xì)節(jié)。熱量傳入顆粒為正傳出為負(fù)這個(gè)約定看似簡(jiǎn)單但在多個(gè)子程序多層調(diào)用后很容易在某個(gè)環(huán)節(jié)乘一個(gè)負(fù)號(hào)導(dǎo)致顆粒溫度不升反降。6. 實(shí)操過程中最常見的問題與排查建議6.1 顆粒溫度整體偏高或偏低出現(xiàn)這種情況首先檢查對(duì)流傳熱系數(shù)模型和顆粒面積的定義。一般來說整體溫度偏高說明熱量輸入偏多要么傳熱面積偏大要么熱流密度偏大。整體溫度偏低則反過來。我實(shí)踐中的排查順序是這樣的先檢查顆粒初始化溫度是否合理檢查傳熱面積是否采用了表面積而不是投影面積檢查努塞爾數(shù)關(guān)聯(lián)式中的物性取值溫度是否正確檢查時(shí)間步長(zhǎng)是否過大導(dǎo)致顯式更新的數(shù)值耗散多數(shù)情況下問題出在物性參數(shù)的溫度取值基準(zhǔn)不統(tǒng)一。比如氣體熱導(dǎo)率如果取的是入口溫度下的值而實(shí)際顆粒表面溫度遠(yuǎn)高于入口溫度那么兩者差異會(huì)導(dǎo)致計(jì)算出的傳熱系數(shù)和真實(shí)偏差很大。6.2 顆粒溫度出現(xiàn)不正常的局部尖峰局部尖峰通常和接觸導(dǎo)熱有關(guān)。當(dāng)兩個(gè)顆粒發(fā)生碰撞接觸半徑若被高估接觸導(dǎo)熱量會(huì)在一瞬間顯著增大表現(xiàn)為局部溫度尖峰。還有一個(gè)容易被忽略的原因是接觸判斷的邏輯利用了接觸重疊量而重疊量在碰撞力模型中可能因?yàn)闀r(shí)間步長(zhǎng)過大出現(xiàn)振蕩進(jìn)而導(dǎo)致接觸導(dǎo)熱熱流也振蕩。處理思路是降低DEM時(shí)間步長(zhǎng)觀察尖峰是否減弱如果減弱說明是接觸力-導(dǎo)熱耦合導(dǎo)致的數(shù)值問題。如果降低步長(zhǎng)無效則檢查接觸半徑計(jì)算中是否合理限幅。6.3 顆粒溫度比流體溫度還平穩(wěn)如果顆粒溫度基本不變而流體溫度在變化檢查顆粒與流體之間是否真的啟用了對(duì)流傳熱。還有可能是因?yàn)轭w粒熱容設(shè)置得非常大比如誤將顆粒體積熱容寫成了比熱容導(dǎo)致熱量輸入對(duì)溫度變化貢獻(xiàn)極小。顆粒熱容的量綱在MFiX中需要特別注意。有些模型輸入的是體積熱容J/m3·K有些輸入的是質(zhì)量熱容J/kg·K弄混之后溫度響應(yīng)會(huì)差好幾個(gè)數(shù)量級(jí)。我建議在讀代碼時(shí)順手把這個(gè)量綱鏈條推導(dǎo)一遍熱流量單位W時(shí)間步長(zhǎng)單位s熱量單位J溫度變化量單位K溫度變化量 熱量 / (顆粒質(zhì)量 × 顆粒比熱容)。如果發(fā)現(xiàn)計(jì)算結(jié)果和這個(gè)關(guān)系對(duì)不上那一定是代碼中定義的某個(gè)熱容或溫度相關(guān)變量理解有誤。6.4 傳熱開關(guān)全開了但就是不算熱這類問題多半是邏輯開關(guān)的命名理解有偏差。比如MFiX中某些開關(guān)是默認(rèn)關(guān)閉需要在輸入文件中顯式開啟而另一些開關(guān)默認(rèn)是開啟的。還有一部分開關(guān)之間存在從屬關(guān)系比如總體能量方程開關(guān)不開即使顆粒能量方程開關(guān)開啟也無效。排查這類問題我通常會(huì)在calc_thermo_des.f的入口處加一個(gè)臨時(shí)輸出打印當(dāng)前步數(shù)、顆粒數(shù)、開關(guān)狀態(tài)然后觀察日志。讀取Fortran代碼時(shí)在關(guān)鍵位置臨時(shí)插入WRITE(*,*)是非常常用的調(diào)試手法雖然原始項(xiàng)目代碼中不推薦留這類輸出但臨時(shí)用一下確實(shí)高效。7. 讀這個(gè)文件時(shí)值得順手做的備忘7.1 版本差異備忘MFiX發(fā)展過程中熱學(xué)計(jì)算模塊經(jīng)歷了較大的結(jié)構(gòu)重構(gòu)。老版本里calc_thermo_des.f是獨(dú)立子程序文件新版本里可能被整合進(jìn)模塊中作為subroutine保留。文件名后綴_mod不是簡(jiǎn)單改名而是從子程序轉(zhuǎn)成了模塊化封裝調(diào)用方式也發(fā)生了變化。如果你正在讀的版本中找不到calc_thermo_des.f不要意外搜索calc_thermo_des或者calc_thermo這樣的關(guān)鍵詞大概率能找到對(duì)應(yīng)的模塊。源碼閱讀工具的全局搜索在這里比手動(dòng)翻目錄高效得多。7.2 數(shù)組維度問題MFiX中的顆粒數(shù)組維度通常不是固定常數(shù)而是動(dòng)態(tài)分配。這意味著循環(huán)邊界條件往往和顆粒最大編號(hào)相關(guān)而不是和實(shí)際當(dāng)前顆粒數(shù)相關(guān)。在calc_thermo_des.f中循環(huán)邊界如果寫的是MAX_PIP那就會(huì)遍歷所有可能的顆粒槽位包括尚未激活的顆粒。未激活顆粒的物理量可能包含垃圾值因此代碼內(nèi)部幾乎都通過IF (SOLIDS ...)或者顆粒存在性判斷來保護(hù)防止對(duì)空顆粒做熱量累加。閱讀時(shí)如果發(fā)現(xiàn)循環(huán)體里面沒有存在性判斷就要警惕。這是少數(shù)會(huì)導(dǎo)致結(jié)果莫名包含NaN或無窮大的原因。7.3 并行與串行差異MFiX支持并行計(jì)算calc_thermo_des.f的邏輯在并行情況下執(zhí)行時(shí)可能會(huì)涉及顆粒分布式存儲(chǔ)即每個(gè)進(jìn)程只負(fù)責(zé)自己區(qū)域內(nèi)的顆粒。這樣一來相鄰進(jìn)程的邊界顆粒換熱需要依賴通訊而不是簡(jiǎn)單的本地循環(huán)。如果你只是在串行環(huán)境下調(diào)試則不需要過度關(guān)注并行通訊邏輯。但如果要做大規(guī)模并行模擬就必須搞清楚顆粒熱量更新時(shí)的邊界交換策略否則顆粒溫度在跨界時(shí)會(huì)出現(xiàn)跳變。關(guān)于并行版本的代碼我現(xiàn)在還沒有完全吃透后續(xù)如果繼續(xù)深入讀完MPI通訊相關(guān)模塊再另外成文分享。8. 從代碼閱讀到實(shí)際修改的經(jīng)驗(yàn)小結(jié)讀源代碼最大的收獲不是記住了某個(gè)公式是怎么寫的而是建立起一種“物理量在代碼中如何流轉(zhuǎn)”的直覺。算例出問題時(shí)能快速判斷是物性參數(shù)問題、邊界條件問題還是模型選擇問題。calc_thermo_des.f這個(gè)文件雖然簡(jiǎn)單但它的調(diào)度邏輯幾乎涵蓋了DEM熱能計(jì)算的全部入口。把它徹底弄懂之后再去看具體的傳熱模型細(xì)節(jié)就會(huì)發(fā)現(xiàn)不過是在這個(gè)框架里填入不同的計(jì)算表達(dá)式而已。下面是我在為一個(gè)流化床干燥算例修改代碼時(shí)走過的一小段真實(shí)流程供參考。第一步確認(rèn)模擬目標(biāo)是顆粒含水率變化而不是單純溫度變化。因此需要在能量計(jì)算之外再引入組分輸運(yùn)和蒸發(fā)模型但蒸發(fā)模型的源項(xiàng)必須和calc_thermo_des.f產(chǎn)生的熱量源項(xiàng)做疊加。第二步在calc_thermo_des.f中找到顆粒熱量累加的循環(huán)在循環(huán)末尾單獨(dú)加上一個(gè)蒸發(fā)潛熱的熱量扣除項(xiàng)。這個(gè)扣除項(xiàng)是負(fù)的因?yàn)檎舭l(fā)會(huì)帶走熱量。第三步調(diào)試中發(fā)現(xiàn)顆粒溫度在水分快蒸干時(shí)出現(xiàn)輕微振蕩。原因是最小含濕率判斷和蒸發(fā)源項(xiàng)切換之間的遲滯。解決辦法是添加滯回區(qū)間讓熱源項(xiàng)平滑過渡。這段經(jīng)歷想說明的是讀源代碼不是目的達(dá)成自己的物理模型目標(biāo)才是目的。但如果不讀源代碼連往哪里加蒸發(fā)源項(xiàng)都無從下手。9. 給新手的一些更實(shí)際的閱讀方法9.1 先跑通一個(gè)算例再讀代碼硬啃源代碼尤其是沒有運(yùn)行過案例的情況下很容易迷失在變量和循環(huán)里。我建議先找一個(gè)帶傳熱過程的標(biāo)準(zhǔn)案例MFiX自帶有許多測(cè)試案例跑通一次再帶著“模擬結(jié)果是怎么算出來”的問題去讀calc_thermo_des.f。這時(shí)候變量不再是紙上的符號(hào)而是和實(shí)際模擬曲線聯(lián)系起來的量。比如跑一個(gè)顆粒靜止床層的導(dǎo)熱案例觀察顆粒溫度從初值趨向壁面溫度的過程。再到代碼里找哪個(gè)變量控制了這個(gè)導(dǎo)熱行為會(huì)很直觀。9.2 打日志比看邏輯快Fortran代碼的閱讀有一個(gè)“笨辦法”在關(guān)鍵子程序入口和出口加上臨時(shí)輸出打印顆粒編號(hào)、顆粒溫度、累計(jì)熱流、接觸狀態(tài)等信息。雖然每次運(yùn)行日志量很大但對(duì)于單個(gè)子程序的邏輯驗(yàn)證這種方式遠(yuǎn)比反復(fù)推敲代碼邏輯快。我以前的調(diào)試習(xí)慣是加一個(gè)條件輸出只打印特定顆粒編號(hào)的信息比如IF (NP 100) THEN WRITE(*,*) NP100 T, DES_TEMP_NEW(NP), Q, Q_SOURCE(NP) END IF這樣就能鎖定單個(gè)顆粒觀察它在整個(gè)時(shí)間步進(jìn)中的溫度和熱流變化。這個(gè)方法對(duì)排查顆粒溫度異常特別有效。9.3 善用差異對(duì)比MFiX支持同一算例用不同模型開關(guān)跑出不同結(jié)果。直接對(duì)比這兩種結(jié)果的溫度場(chǎng)差異再回到代碼里找對(duì)應(yīng)條件分支比憑空猜測(cè)分支邏輯更直接。我之前對(duì)比過開啟接觸導(dǎo)熱和關(guān)閉接觸導(dǎo)熱兩種設(shè)置下的床層溫度分布發(fā)現(xiàn)差異集中出現(xiàn)在堆積區(qū)說明接觸導(dǎo)熱模型在堆積區(qū)域起主導(dǎo)作用。據(jù)此再反推代碼理解就快多了。10. 最后想說的源文件雖小串聯(lián)的知識(shí)面卻不小calc_thermo_des.f在MFiX整個(gè)項(xiàng)目里算是一個(gè)規(guī)模不大的子程序。但它串聯(lián)起了顆粒相運(yùn)動(dòng)、流體與顆粒對(duì)流傳熱、顆粒間接觸導(dǎo)熱、壁面換熱和輻射換熱這些模型。如果只從代碼規(guī)模上判斷它不重要就很容易錯(cuò)過DEN模塊熱能計(jì)算體系的關(guān)鍵入口。我現(xiàn)在讀到這個(gè)階段越來越覺得讀源代碼像在拆一個(gè)復(fù)雜的機(jī)械設(shè)備看起來只有一臺(tái)機(jī)器但每個(gè)零件都有自己的角色拆得越多對(duì)整個(gè)系統(tǒng)的理解越完整。calc_thermo_des.f就是其中一個(gè)不大但絕對(duì)不能說無關(guān)緊要的零件。如果你也在讀MFiX源代碼建議以后遇到這個(gè)文件時(shí)不要跳過花一兩個(gè)小時(shí)把它的調(diào)用關(guān)系和變量流轉(zhuǎn)弄清楚。之后再看任何和顆粒傳熱相關(guān)的代碼都會(huì)有一種豁然開朗的感覺。