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