基于顆粒法的物質熱傳導建模與可視化仿真

陳曉峰，黃剛海，邱泓杰

（1.湖南科技大學 巖土工程穩定控制與健康監測湖南省重點實驗室，湖南 湘潭 411201；2.湖南科技大學 土木工程學院，湖南 湘潭 411201；3.中南大學 土木工程學院，湖南 長沙 410075；4.中南大學 資源與安全工程學院，湖南 長沙 410075）

1 研究綜述

顆粒材料與人們的日常生活密切相關，如土壤、砂石、糧食等。顆粒材料的熱傳導研究有著廣泛的應用前景。如在核廢料地質處置中，深埋地下核廢料中的核元素衰變會釋放大量熱，通過熱傳導使圍巖溫度升高，由熱應力引起的圍巖破裂可能導致核廢料泄漏，污染水資源；地熱具有清潔、可循環和安全等優點，是較好的可再生能源之一。在地熱開采中，將冷水注入熱巖，使巖石溫度下降，導致巖石發生熱破裂，這一問題也涉及巖石熱傳導。此外，地殼中礦床的形成規律研究、礦床勘探、二氧化碳的地質存儲和油氣開采等新領域都與顆粒材料熱傳導聯系緊密。因此，對顆粒的熱傳導研究有著重要的應用價值。

顆粒系統熱傳導問題研究分為兩類，即連續介質傳熱模型和離散顆粒傳熱模型。連續介質傳熱模型將系統的宏觀性質及宏觀傳遞過程與微觀結構、微觀特性相聯系，常用的有由I.Nozad 等[1]提出的有效介質法，用各向同性的有效介質代替顆粒系統，將材料的微觀結構特性簡化得到平均溫度分布，其不足之處在于忽略了顆粒系統的高度離散性；相反，基于微觀層次用離散模型研究顆粒材料的傳熱過程不存在此類問題。對于離散顆粒傳熱模型，S.Yagi等[2]最先提出了離散顆粒傳熱的概念，認為顆粒之間可通過5 種途徑傳熱。針對各傳熱途徑之間的差異，學者們提出了多種模型及相應的計算表達式。如M.Yovanovich[3]指出，接觸面間的導熱是顆粒系統傳熱的主要方式。Sun J.等[4]將兩顆粒通過接觸面的導熱問題近似為兩個半無限厚介質間的熱傳遞問題，但作了過多簡化。W.L.Vargas 等[5-6]研究發現，在理想情況下，接觸面上傳遞的熱量與顆粒間的法向應力存在一定的函數關系，在此基礎上發展了熱顆粒動力學法（thermal particle dynamics，TPD），但其在計算非固定顆粒時具有局限性。B.Chaudhuri[7]利用離散單元法（discrete element method，DEM），結合TPD，對轉鼓中顆粒在轉動情況下的傳熱問題進行了研究。Shi D.L.[8]結合DEM 與計算流體動力學（computational fluid dynamics，CFD），模擬了不同材料、不同顆粒大小情況下回轉窯中顆粒的導熱問題。Feng Y.T.等[9-10]以兩個顆粒間的接觸為基本單元，將顆粒體模型化為傳熱網絡，發展了顆粒材料中熱傳導問題的離散溫度單元法（discrete thermal element method，DTEM），并對較大規模的二維顆粒系統導熱問題進行了研究，但其計算傳熱需要引入兩個未知量。顆粒法（particle simulation method，PSM）將介質材料離散為一系列通過接觸連接在一起的顆粒結合體，可用于求解顆粒材料的熱傳導問題，已發展了成熟的顆粒法軟件PFC（particle flow code）[11]，現已被廣泛應用于顆粒材料特性及各類工程問題研究中。

與此同時，國內也有學者開展了關于顆粒熱傳導方面的研究。劉安源等[12]在分析顆粒碰撞傳熱機理的基礎上建立了一個新的顆粒碰撞傳熱理論模型，對流化床內床層與壁面的碰撞傳熱系數進行顆粒層次上的模擬。武錦濤等[13]綜合考慮顆粒接觸傳熱的多個過程，提出了顆粒接觸傳熱模型，并結合DEM模擬了移動床中顆粒與加熱面間的傳熱過程。李臘梅等[14]提出了一種全新的非連續介質熱傳導數值計算方法，并通過自主編制計算程序驗證了其正確性。于慶波等[15]采用CFD 軟件，對顆粒繞圓管傳熱現象進行了建模，并驗證了模型的正確性。王擎等[16]將離散元方法與顆粒熱傳導模擬結合，研究了回轉干餾爐內顆粒間的傳熱過程。姬國釗等[17]通過離散元方法分析了回轉窯內顆粒的運動和換熱過程。

從國內外研究現狀看，熱傳導的研究常采用商業軟件模擬，獨立開發自主計算程序甚少。而對熱傳導研究方法而言，離散顆粒模型，因考慮了實際顆粒組成結構，能更好地反映實際過程，受到了更廣泛的關注。鑒于熱傳導研究現狀及離散顆粒模型的優點，本論文以顆粒法為傳熱算法，自主編寫C++代碼，開發熱傳導計算程序，對顆粒材料進行熱傳導建模及可視化仿真分析。

2 顆粒法基本原理

2.1 熱傳導基本方程

設顆粒系統中任一顆粒為單元i，溫度的變化規律滿足連續介質熱傳導方程：

式中：ρi為顆粒密度；Cv為材料比熱容；qi為熱流密度；qv為單位體積發熱率；Ti為顆粒溫度。

傅里葉熱傳導定律定義了熱流密度與溫度梯度間的關系：

式中：x為熱流方向；k為熱傳導系數。

常用的邊界條件分為3 類：

1）Dirichlet 邊界，指定邊界上的溫度值。

2）Neumann 邊界，指定邊界上的熱流密度值。

3）Robin 邊界，指定邊界上的對流換熱系數及邊界層溫度值。

本研究中將顆粒材料離散為一系列的熱存儲器和熱管道，熱量將通過熱存儲器之間的激活管道進行傳遞。熱管道的有效性與接觸有關，如果接觸過程中兩顆粒單元有重疊則認為連接兩顆粒間的熱管道有效，熱傳導接觸模型示意圖如圖1所示。

圖1 熱傳導接觸模型示意圖Fig.1 Schematic diagram of heat conduction contact model

從圖1所示模型可知，熱管道連接兩個熱存儲器，熱流僅在熱管中流動，兩熱存儲器轉移的熱量即為熱管的能量ΔQ，熱管有固定的熱阻η，將每個熱管看成單位長度上熱阻為η的一維物體，則有：

式中：ΔT為熱管道兩端熱存儲器的溫度差；L為熱管道的長度。

對顆粒單元i進行分析，可得：

式中：ηi-j與Li-j分別為連接單元i與單元j間熱管的熱阻與長度。

同理可得：

測量顆粒單元中熱傳導系數與熱管道中的熱阻抗系數有如下理論關系式[11]：

式中：p為顆粒材料的孔隙率；V為顆粒的體積。

2.2 熱傳導數值離散

給定初始條件下的溫度，可計算出每個熱管的能量變化量，從而熱源的溫度將改變，變化公式如下：

對式（10）中時間導數向前差分，可得：

式中Δt為時間步長。

2.3 熱傳導計算流程

基于上述顆粒法基本原理和熱傳導接觸模型，通過所有接觸對熱傳導的累積，可實現整個顆粒系統的熱傳導。計算流程如圖2所示，其中i_step為計算步標識，i_par為顆粒單元標識。每執行一次循環更新一次顆粒單元溫度，當達到時間總長時輸出結果，計算結束。

圖2 熱傳導計算流程圖Fig.2 Calculation flow chart

3 可視化模型建立

熱傳導的程序架構和軟件界面如圖3所示。程序主要分為計算內核模塊和可視化模塊兩部分。

圖3 熱傳導程序架構和軟件界面Fig.3 Program architecture and software interface

計算內核模塊將熱傳導相關數據和操作封裝成一個C++類——CHeatConduction 類，用于執行所有與圓盤顆粒熱傳導有關的運算。CHeatConduction 類的數據結構如圖4所示。

圖4 CHeatCondution 類數據結構Fig.4 Data structure of CHeatCondution

由圖4可知，該類用小類對象作為其成員變量，主要包括CDisk 類——圓盤類，每個CDisk 類對象表示圓盤顆粒系統中的一個圓盤；其成員變量主要包括圓盤溫度和圓盤溫度變化量，圓盤中心坐標和半徑；圓盤密度和質量以及與該圓盤接觸的相鄰圓盤數據等。pD 為CheatConduction 類對象數組，存儲了圓盤顆粒系統中所有圓盤的數據。

熱傳導計算主要包含如下操作。

Step 1 執行setInitialTemperature( )函數，設置初始溫度；

Step 2 執行searchContactPairs( ) 函數，訪問pD 數組，搜索系統中的所有圓盤-圓盤接觸對，根據接觸對數據求解出與傳熱有關的變量(如嵌入量等)；在計算中以上步驟只需執行一遍，這是因為本論文中顆粒固定不動，因此顆粒間的接觸方式及與傳熱有關的嵌入量等均保持不變。然后執行單位時步的計算，每一單位時步內的計算分為兩步：

Step 1 執行implementCalculation( )函數，求解每一組接觸對的溫度變化量，更新改變后的圓盤溫度值，并將求解結果保存至圓盤數組中；

Step 2 執行updateModelData( )更新顆粒溫度。

其中，可視化模塊被封裝成了另一個C++類——COpenGL 類，該類從Qt 預定義的QOpenGLWidget類和QOpenGLFunctions 類派生，并重寫其內嵌虛函數，用于渲染程序界面的可視化視窗。CHeatConduction 類和COpenGL 類通過指針相互訪問，COpenGL類可讀取CHeatConduction類的數據（如顆粒的溫度、位置坐標及顆粒大小等）并將模型實時顯示在計算機屏幕上。

4 模型驗證與結果分析

通過分析一維薄板、二維對稱圓環的熱傳導特性對其進行正確性驗證。本研究算例所用參數如下：熱傳導系數k=1.6 W/(m·℃)；比熱容Cv=0.2 J/(kg·℃)；密度ρ=1 000 kg/m3；初始溫度T0=0 ℃。

4.1 一維薄板熱傳導

一維平面薄板熱傳導示意圖如圖5所示。由圖5a 可知，平面薄板的長度L=20 m，寬度W=0.2 m。采用100 個圓盤搭建一維薄板熱傳導模型。模型左右兩側設定不同的常溫度邊界，分別為T1=100 ℃，T2=0 ℃。

圖5 一維薄板熱傳導模型與變化規律Fig.5 One-dimensional thin plate heat conduction

該薄板在固定溫度邊界下逐漸達熱平衡狀態。根據文獻[18]在“常溫度邊界下薄板的瞬態熱響應”中的推導，薄板各處的溫度T(Z,t)為

式中：t為時間；Z為距離模型左側距離。

圖5b 展示了自主程序算出的薄板各處溫度及對應的解析解，可看出數值解與理論解高度吻合。

4.2 二維對稱圓環熱傳導

二維對稱圓環熱傳導示意圖如圖6所示。

圖6a 可知，圓環的內徑Ra=2 m，外徑Rb=5 m，該模型由6 596 個半徑為0.05 m 的圓盤組成。在模型內、外側分別設定100 ℃，0 ℃的常溫度邊界。

由圖6a 中的圓環模型，當熱流穩定時各處溫度T(r)的解析解[18]為

式中r為監測點到圓心的距離。

圖6b 為自主程序計算出的圓環各處溫度及對應的理論解，從中可以看出數值解與理論解高度吻合。

圖6 二維對稱圓環熱傳導模型與變化規律Fig.6 Two-dimensional symmetrical ring heat conduction

從以上一維薄板及對稱圓環兩個算例的結果可以看出：自主程序計算正確，且精度較高。

4.3 影響因素分析

4.3.1 顆粒粒徑

圖7為粒徑不同的兩個顆粒模型，在模型左右兩側分別設定100 ℃和0 ℃常溫度邊界進行熱傳導模擬。對兩個模型分別監測，取監測點的位置為距離模型左側4,8,16 m 處，所得結果如圖8所示。從圖8中的曲線圖可以看出：顆粒半徑對傳熱的影響非常小，大顆粒模型傳熱略快，但顆粒系統最終的穩定傳熱狀態不受顆粒粒徑的影響。

圖7 兩種不同粒徑的顆粒模型Fig.7 Two particle models with different particle sizes

圖8 顆粒粒徑對傳熱的影響Fig.8 Influence of particle size on heat transfer

4.3.2 堆積方式

本研究建立了3 種不同的顆粒堆積模型，如圖9所示。在堆積方式不同的3 種顆粒模型左右兩側分別設定100 ℃和0 ℃常溫度邊界，進行熱傳導模擬。對3 種模型取不同的位置點分別進行監測，所得結果如圖10 所示。從圖10 可以看出：堆積方式對傳熱的影響非常小，方形堆積傳熱略快，但顆粒系統最終的穩定傳熱狀態不受顆粒堆積方式影響。

圖9 3 種堆積方式的顆粒模型Fig.9 Disks model formed by three stacking methods

圖10 堆積方式對傳熱的影響Fig.10 Influence of stacking mode accumulation methods

5 結論

本文討論了模擬顆粒材料熱傳導的顆粒法接觸模型、接觸熱交換計算方法及可視化實現算法的基本實施步驟。在此基礎上，通過分析一維薄板、二維對稱圓環的熱傳導特性，驗證了模型的正確性。然后分析了顆粒粒徑、堆積方式對傳熱的影響，得出如下結論：1）在其他參數一致時，顆粒半徑對傳熱的影響非常小，大顆粒模型傳熱略快，但顆粒系統最終的穩定傳熱狀態不受顆粒粒徑影響；2）顆粒堆積方式對傳熱的影響較小，方形堆積傳熱略快，但顆粒系統最終的穩定傳熱狀態不受顆粒堆積方式影響。

以上研究結果表明，所提方法能較好地模擬顆粒的熱傳導，但本研究也有一定局限性，僅對一維薄板、二維圓環進行了研究，對三維立體顆粒物質的傳熱還有待進一步研究。