，，，?，，

（1．環(huán)境保護部核與輻射安全中心，北京100082；2．清華大學(xué)核能與新能源技術(shù)研究院，北京100084）

高溫氣冷堆球床接觸導(dǎo)熱計算方法

李聰新1，任成2，許超1，?，溫麗晶1，劉宇生1

（1．環(huán)境保護部核與輻射安全中心，北京100082；2．清華大學(xué)核能與新能源技術(shù)研究院，北京100084）

接觸導(dǎo)熱是低溫條件下高溫堆球床堆芯傳熱的主要方式，接觸導(dǎo)熱的計算精度決定了低溫條件下高溫氣冷堆球床傳熱計算的精度。本文建立了類似熱離散單元法的球床接觸導(dǎo)熱的計算模型，通過計算流體力學(xué)方法進行模型中兩個球之間理想接觸的導(dǎo)熱量計算，給出了所編寫的大規(guī)模隨機堆積球床接觸導(dǎo)熱計算程序的詳細(xì)步驟。通過與計算流體力學(xué)方法對球床局部堆積結(jié)構(gòu)接觸導(dǎo)熱計算結(jié)果的比較，驗證了該計算程序的準(zhǔn)確性。該方法與球床局部結(jié)構(gòu)計算流體力學(xué)方法精度接近，并可用于顆粒尺度大規(guī)模球床計算。

高溫堆；球床；接觸導(dǎo)熱；顆粒尺度；理性接觸；熱離散單元法

高溫氣冷堆球床堆芯是由球形燃料元件隨機堆積而成的［1］，在低溫條件下，球床的主要傳熱方式為接觸導(dǎo)熱［2］，球床接觸導(dǎo)熱的計算精度決定了低溫條件下高溫氣冷堆球床傳熱計算的精度。目前，從顆粒尺度進行球床傳熱計算的主要方法為熱離散單元法［3］（TDEM，Thermal Discrete Element Method）。TDEM是在顆粒流動計算的離散單元法 （DEM，Discrete Element Method）基礎(chǔ)上增加了傳熱模型［4］，將每個球用一個平均溫度表示，從顆粒尺度進行球床傳熱計算［5］，從而得到每個球的溫度和傳熱量，既可以研究球床的局部傳熱特性［6］，又可以采用統(tǒng)計平均方法研究球床等效導(dǎo)熱系數(shù)［7］、熱流量等宏觀參數(shù)，是目前球床流動傳熱計算研究的熱點［8，9］。熱離散單元法的主要問題是缺乏嚴(yán)格的數(shù)學(xué)證明，模型中存在大量的工程假設(shè)和經(jīng)驗參數(shù)，所以熱離散單元法可以認(rèn)為是一種工程方法，而不是嚴(yán)格的理論計算方法。目前開源的DEM軟件有LIGGGHTS［10］和OpenFOAM［11］，商業(yè)軟件有EDEM和PFC3D等。在這些軟件中，顆粒間的接觸導(dǎo)熱量通常表示為：

ΔT為兩個接觸球的溫差，hc為接觸面?zhèn)鳠嵯禂?shù)，其表達式為：

λs1、λs2為兩個接觸顆粒各自的材料導(dǎo)熱系數(shù)；FN為重力；r?為兩個接觸顆粒的幾何平均半徑；E?為等效楊氏模量；括號中的部分表示顆粒間的接觸面積。從公式 （2）可以看出，接觸面?zhèn)鳠嵯禂?shù)表示為了導(dǎo)熱系數(shù)和接觸面積的一個簡單函數(shù)，而接觸面積則是按照彈性力學(xué)假設(shè)的一個經(jīng)驗關(guān)系式。該公式將傳熱和受力變形的經(jīng)驗關(guān)系式結(jié)合在一起，導(dǎo)致接觸傳熱系數(shù)的計算精度難以控制，用熱離散單元法進行的球床傳熱計算結(jié)果也將難以驗證。

本文在接觸導(dǎo)熱模型中將受力變形假設(shè)和傳熱計算分離，只進行理想接觸條件下的傳熱計算，通過與計算流體力學(xué)方法對球床局部堆積結(jié)構(gòu)接觸導(dǎo)熱計算結(jié)果的比較，驗證了計算程序的準(zhǔn)確性。該方法可用于顆粒尺度大規(guī)模球床導(dǎo)熱計算，并得到與球床局部結(jié)構(gòu)計算流體力學(xué)方法接近的精度。在進行實際球床傳熱計算時，可以通過實際接觸面積與理想接觸面積的比值耦合力學(xué)參數(shù)。

l 接觸導(dǎo)熱計算模型

在低溫條件下，球表面的輻射傳熱可以忽略，則球之間通過接觸面的導(dǎo)熱將成為主要傳熱方式。將每一個球定義為一個單元，每個球的溫度用球的平均溫度Ti代表，定義流入球的熱流方向為正，相鄰兩個球之間的導(dǎo)熱量為：Ti為當(dāng)前球的平均溫度；Tj為與當(dāng)前球接觸的球的溫度；hci定義為接觸面?zhèn)鳠嵯禂?shù)。

根據(jù)能量守恒，傳入球內(nèi)的總熱量等于球的內(nèi)能變化：

其中，n為與當(dāng)前球接觸的球的個數(shù)。Cs為球的熱容，定義為球的質(zhì)量與比熱容的乘積：

當(dāng)處于穩(wěn)態(tài)時，球的凈熱流量為零，則：

對于靠近壁面的球，還需要相應(yīng)的邊界條件才能使方程組閉合。根據(jù)目前球床導(dǎo)熱計算的需要，只處理給定壁溫和絕熱兩類邊界條件。

當(dāng)球處于給定壁溫邊界時：hci為接觸半徑相同的球之間的接觸面?zhèn)鳠嵯禂?shù)，Tw為接觸壁面的溫度。即在同等接觸面積時，球和壁面之間的接觸面?zhèn)鳠嵯禂?shù)為兩個球之間的接觸面?zhèn)鳠嵯禂?shù)的兩倍。

當(dāng)處于絕熱邊界條件時，相當(dāng)于減少了一個熱流通道：

本文的接觸導(dǎo)熱模型與一般熱離散單元法的導(dǎo)熱模型類似，區(qū)別在于接觸面?zhèn)鳠嵯禂?shù)的處理上。熱離散單元法將傳熱和受力變形的經(jīng)驗關(guān)系式結(jié)合在一起，導(dǎo)致接觸傳熱系數(shù)計算精度難以控制。本文將受力變形假設(shè)和傳熱計算分離，進行單純的接觸導(dǎo)熱計算，受力變形關(guān)系由離散單元法在確定堆積結(jié)構(gòu)時確定。球之間的接觸導(dǎo)熱采用理想接觸面積計算，將由表面粗糙度引起的

Rc是半個球的接觸熱阻。研究兩個球之間接觸熱阻的文章很多，早期多采用解析方法和簡化假設(shè)得到近似解析解［12，13］，但不同解析解的差異較大，且都存在較大誤差。文獻［14］分析了導(dǎo)致差異和誤差的原因，并通過計算流體力學(xué)方法精細(xì)計算了不同接觸半徑的球之間的導(dǎo)熱量，得到了一個較為精確的接觸面?zhèn)鳠崦嫦禂?shù)表達式：

其中r為球的半徑，λs為球的導(dǎo)熱系數(shù)，Ks是一個由相對接觸半徑確定的形狀函數(shù)，其通過計算流體力學(xué)方法計算半徑為1m、導(dǎo)熱系數(shù)為1W·m-1·K-1、相對接觸半徑為γ的兩個半球在一定溫差下的導(dǎo)熱量得到：實際接觸面積減小用一個0和1之間的系數(shù)Kn加以考慮。

在球床接觸導(dǎo)熱計算中，關(guān)鍵在于獲得兩個球之間的接觸面?zhèn)鳠嵯禂?shù)hci。接觸面?zhèn)鳠嵯禂?shù)的引入是為了計算方便，其定義為兩個半球之間接觸熱阻的倒數(shù)，即：

形狀函數(shù)Ks確定以后，理想接觸的兩球之間的導(dǎo)熱量可以通過公式 （10）確定。當(dāng)考慮實際物體的表面粗糙度時，實際導(dǎo)熱面積小于理想接觸面積，因此引入一個表示實際接觸面積和理想接觸面積的比例系數(shù)Kn，其取值范圍為0到1之間。因此，表示兩個球之間接觸面?zhèn)鳠嵯禂?shù)的完整公式為：

2 計算程序開發(fā)

在球床的接觸導(dǎo)熱計算中，將每個球的溫度用一個平均溫度表示，每個球的能量守恒可以確定一個線性方程，所有球的線性方程組成一個閉合的線性方程組。由于球床的隨機堆積狀態(tài)，每個球所接觸的球的個數(shù)和方向不同，導(dǎo)熱模型的計算過程相當(dāng)于一個非結(jié)構(gòu)化網(wǎng)格的固體導(dǎo)熱計算。球床的接觸導(dǎo)熱計算的完整計算程序包括前處理、方程求解、后處理三個部分，以下對本文所開發(fā)的程序的這三個部分進行簡要說明。

2.l前處理

前處理用于獲得球所接觸的球和壁面信息，相當(dāng)于流體力學(xué)計算中的網(wǎng)格信息。由于離散單元法程序只是給出了一定堆積結(jié)構(gòu)內(nèi)的球的位置坐標(biāo)，需要從球的坐標(biāo)中提取球間的接觸和邊界信息。以立方體隨機堆積為例，首先對球的坐標(biāo)沿某一方向排序 （可取主熱流方向），然后對每個球按排序結(jié)果編號，球的編號是求解計算的基礎(chǔ)，其與方程和節(jié)點相對應(yīng)，對于由一萬個球堆積成的球床，球的編號即為1～10 000。完成編號以后，對球床中任意兩個球的球心距離做一次比較，當(dāng)小于球直徑時即認(rèn)為兩個球接觸。在對球心距離進行比較的過程中，如圖1所示，需要存儲當(dāng)前球的接觸球總數(shù)、每個接觸球的編號和相對接觸半徑。最后進行一次壁面邊界信息搜索，比如對于球心Z方向坐標(biāo)小于一倍球半徑的球認(rèn)為與下壁面接觸，確定與邊界接觸的球的總個數(shù)以及每個邊界球的編號和相對接觸半徑。在立方體堆積球床中，有六個壁面，在環(huán)形球床中有四個壁面，壁面邊界信息中還應(yīng)存儲邊界條件信息，如溫度邊界條件或絕熱邊界條件。在獲得了每一個球和壁面的接觸信息之后，前處理過程完成。

圖l 每個球和壁面的存儲信息Fig.l Stored information of each pebble and wall

2.2 方程求解

方程求解的關(guān)鍵是通過球和壁面信息組成一個線性方程組，形式為：

A為系數(shù)矩陣，維度為N×N，N為球床中球的個數(shù)，x是N×1矩陣，即未知溫度，b方程中與未知溫度無關(guān)的常數(shù)，由邊界條件確定。由于上述方程可采用任何線性方程組的數(shù)值解法進行求解，因此計算的難點在于組成矩陣A和b。

根據(jù)前處理得到的球和壁面信息，可以通過兩個大的循環(huán)完成A和b的計算。初始時將A和b的全部元素值賦零。首先是對球的循環(huán)，對于每個球，根據(jù)公式 （12）確定當(dāng)前球i與鄰球j的接觸面?zhèn)鳠嵯禂?shù)，j為鄰球的編號。矩陣A的主對角線 （i，i）的值為其與所有鄰球的接觸面?zhèn)鳠嵯禂?shù)和的負(fù)值，矩陣A的 （i，j）元素的值為當(dāng)前球與鄰球j的接觸面?zhèn)鳠嵯禂?shù)。

其次是對壁面的循環(huán)，給定溫度和絕熱邊界條件的處理方法不同。對于給定溫度的邊界條件，每個壁面球按公式 （7）確定其與壁面的接觸面?zhèn)鳠嵯禂?shù)。對于編號為j的壁面球：

hwj為當(dāng)前球和壁面的接觸面?zhèn)鳠嵯禂?shù)，Tj為壁面的溫度。A（j，j）n-1為在此次循環(huán)之前確定的系數(shù)矩陣中主對角線上元素的值，b（j）n-1為當(dāng)前循環(huán)之前b（j）的值。

對于給定的絕熱邊界條件，其相當(dāng)于對該接觸點對球能量守恒方程無貢獻，直接忽略該循環(huán)即可。

在完成了上述兩個循環(huán)之后，系數(shù)矩陣A和右端項b即計算完成。在一般流體力學(xué)計算中，系數(shù)矩陣的非零元素在二維情況下為五對角矩陣，在三維情況下為七對角矩陣。與一般流體力學(xué)固體導(dǎo)熱計算的系數(shù)矩陣不同，由于球床中球接觸的隨機性，系數(shù)矩陣也具有一定的隨機性。以前文模擬的立方體隨機堆積球床為例，系數(shù)矩陣A（左上角的部分）的非零元素如圖2所示，形成一個寬度約為600的隨機對對角線，每一行中非零元素為2～13個，由于壁面效應(yīng)造成的配位數(shù)變化也可以在一定程度上反應(yīng)在帶狀邊界寬度上。由于系數(shù)矩陣A為嚴(yán)格對角占優(yōu)矩陣，因此方程組是穩(wěn)定可解的。為提高計算效率，可采用一些大型稀疏矩陣算法，如multifrontal算法［15］。

圖2 系數(shù)矩陣A中非零元素位置Fig.2 Positions of non-zero elements in the coefficient matrix A

2.3 后處理

后處理主要是根據(jù)計算得到的每個球的溫度和熱流量進行相關(guān)物理參數(shù)的計算和統(tǒng)計，如球床溫度分布、邊界熱流量、球床等效導(dǎo)熱系數(shù)等。由于球床內(nèi)每個球的溫度都已經(jīng)確定，這些物理參數(shù)不過是一個求和、平均的統(tǒng)計過程。

以立方體內(nèi)的球在上下兩側(cè)給定溫度、四周絕熱條件下的熱流量和等效導(dǎo)熱系數(shù)為例。給定溫度的壁面導(dǎo)熱量的計算與第一類邊界條件設(shè)置算法類似，每個與壁面接觸的球與壁面之間傳遞的熱量為：

其中hwi為球與壁面的接觸面?zhèn)鳠嵯禂?shù)，Tw為設(shè)定的壁面溫度，Twi為計算得到的與壁面接觸的球的溫度。通過該壁面的總的熱流量為：

其中ΔT為球床兩次溫度差，ΔL為球床的高度，A為球床的橫截面積。

3 模型驗證

n為前處理階段得到的與該壁面接觸的球個數(shù)。球床的等效導(dǎo)熱系數(shù)為：

3.l一維驗證

將10個直徑為60 mm，導(dǎo)熱系數(shù)為100 W· m-1·K-1的石墨球等間距連成一排，球四周絕熱，左右兩側(cè)所接觸壁面的溫度分別為0℃和100℃，通過改變球心距離計算在不同相對接觸半徑下通過單個球的熱流量。球心距離在12～58 mm之間以1 mm為步長遞增，在58．1～59．9 mm之間以0．1 mm為步長遞增。

采用計算流體力學(xué)軟件進行計算， （如CFX、Fluent等）計算結(jié)果作為比較對象，在不同接觸半徑下，球床接觸導(dǎo)熱模型和計算流體力學(xué)軟件得到的通過每球的熱流量如圖3所示，在計算范圍內(nèi)。兩者計算得到的熱流量的相對誤差最大值不超過2%，導(dǎo)熱模型和計算流體力學(xué)方法的計算結(jié)果應(yīng)保持一致，但可以看出，二者之間仍然存在一定誤差。誤差產(chǎn)生的原因是球接觸熱阻形狀函數(shù)是通過數(shù)值方法得到的，其計算結(jié)果與網(wǎng)格疏密有關(guān)，在計算形狀函數(shù)時，單一球的等效網(wǎng)格數(shù)目為千萬量級，遠(yuǎn)遠(yuǎn)大于10個球計算時每個球的網(wǎng)格數(shù)目。若可以進一步增加網(wǎng)格數(shù)量，兩種計算方法的計算結(jié)果應(yīng)趨于一致。

圖3 不同接觸半徑下通過單個球的熱流量比較Fig.3 Comparison of heat flow through one pebble with different contact radius

3.2三維驗證

為了驗證球場接觸導(dǎo)熱計算方法可以應(yīng)用于三維堆積結(jié)構(gòu)中，本文仍然采用與計算流體力學(xué)軟件直接計算結(jié)果作比較的方法。采用60 mm、導(dǎo)熱系數(shù)為100 W·m-1·K-1的石墨球通過立方體 （SC）堆積而成，X、Y、Z方向球的個數(shù)分別為6、5、4。球床六個表面同時施加不同的溫度邊界條件，其中X方向為0℃和400℃，Y方向為100℃和500℃，Z方向為200℃和600℃，比較球內(nèi)溫度為三維分布時接觸導(dǎo)熱模型與計算流體力學(xué)軟件計算得到的球的溫度和熱流量。由于120個球的計算流體力學(xué)軟件固體導(dǎo)熱計算耗時較長，因此不再進行每一接觸半徑下的比較。由于隨著相對接觸半徑減小，同樣的導(dǎo)熱量計算精度所需的網(wǎng)格數(shù)據(jù)迅速增加，因此只比較球心距離為58 mm，實際接觸半徑為7．68 mm，相對接觸半徑為0．256的情況，此時所需網(wǎng)格數(shù)目尚在目前16G內(nèi)存服務(wù)器的計算能力范圍內(nèi)。

由于接觸導(dǎo)熱模型基于一維球接觸熱阻得到，其在三維條件下是否適用是模型能否用于隨機球床導(dǎo)熱計算的關(guān)鍵。因此，上述邊界條件的設(shè)置是為了使球床中的溫度處于完全的三維分布狀態(tài)。隨機取出堆積區(qū)域中間的兩個球，計算流體力學(xué)軟件計算得到的等溫面和熱流線如圖4所示。由于接觸點的溫度各不相同，球內(nèi)部的溫度場已經(jīng)完全處于三維分布狀態(tài)，從隨機取出的兩個球看，其等溫面已經(jīng)不存在相似之處。

圖4 球等溫面和熱流線的計算流體利力學(xué)軟件計算結(jié)果Fig.4 The isothermic surface and Heat flux lines of CFD software Calculated result

采用計算流體力學(xué)方法和球床接觸導(dǎo)熱方法計算得到的溫度分布如圖5所示，從整體看，兩種方法得到的計算結(jié)果都反映了所設(shè)置的邊界條件特征，溫度最小值出現(xiàn)在X=0的平面，溫度最大值出現(xiàn)在最上層。由于堆積結(jié)構(gòu)的六個表面設(shè)置的溫度各不相同，球床中的球溫度均顯示出緩慢漸變的特點。導(dǎo)熱模型相當(dāng)于計算了球的平均溫度，不能反映球在接觸點附近的溫度變化，例如在X=0的平面上設(shè)置的溫度為0℃，可以從計算流體力學(xué)軟件計算結(jié)果中明顯看出X=0的平面上球接觸點的溫度均為0℃，而導(dǎo)熱模型計算結(jié)果則為求平均后的溫度。

圖5 三維球床溫度分布比較Fig.5 Comparison of 3D temperature distribution of pebble bed

為了定量比較導(dǎo)熱模型和計算流體力學(xué)軟件計算得到的球溫度，取Z方向最上層的30個球，編號如圖5（b）所示。對計算流體力學(xué)軟件計算結(jié)果采取體積平均的方式計算每個球的平均溫度，與導(dǎo)熱模型得到的球的溫度作比較，結(jié)果如圖6所示。可見，兩種方法得到的球溫度幾乎相同，最大相對誤差不超過1．7%。從圖6可以看到兩種方法計算得到的溫度都存在5個周期變化，每個周期有6個溫度點，正好與球的排列一致。這30個球的溫度變化幅度較大，這是由于所設(shè)置的邊界條件劇烈變化所致，可見采用接觸導(dǎo)熱模型并不限于小溫差假設(shè)。

圖6 Z方向最上層球溫度比較Fig.6 Temperature comparison of the top layer pebble in Z direction

在上述堆積方式中，六個邊界上的總熱流量見表1，其中熱流以流入球內(nèi)的方向為正。在計算流體力學(xué)軟件計算中，邊界總熱流量采用對每個接觸面的熱流密度的積分得到；在球床接觸導(dǎo)熱計算方法中，采用公式 （16）計算，即對與相應(yīng)壁面接觸的每個球的熱流量求和。從表1可以看出，邊界熱流量最大相對誤差為4．7%，考慮到每個球內(nèi)部的三維溫度分布和非常大的溫度梯度，這一精度完全可以接受。

通過接觸導(dǎo)熱模型與計算流體力學(xué)方法在熱流量和溫度計算方面的比較，可見基于一維熱阻推導(dǎo)的接觸導(dǎo)熱模型在三維球床導(dǎo)熱計算中同樣適用。接觸導(dǎo)熱模型可以進行大規(guī)模球床導(dǎo)熱的計算，與計算流體力學(xué)方法的相對誤差控制在5%左右。這一精度是目前熱離散單元法難以達到的。由于現(xiàn)有的計算流體力學(xué)軟件無法進行大規(guī)模的隨機堆積球床的導(dǎo)熱計算，因此不再進行大規(guī)模球床的接觸導(dǎo)熱驗證。

表l 邊界總熱流量比較Table l Comparison of boundary total heat flux

4 結(jié)論

本文建立了類似熱離散單元法的球床接觸導(dǎo)熱的計算方法，將導(dǎo)熱計算和力學(xué)經(jīng)驗關(guān)系式區(qū)分離，相當(dāng)于一種理想條件下球床導(dǎo)熱的高精度計算方法。通過計算流體力學(xué)方法進行了模型中兩個球之間理想接觸的導(dǎo)熱量計算，提高了接觸面?zhèn)鳠嵯禂?shù)計算的凈多。本文給出了所編寫的大規(guī)模隨機堆積球床接觸導(dǎo)熱計算程序的詳細(xì)步驟。通過與計算流體力學(xué)方法對球床局部堆積結(jié)構(gòu)接觸導(dǎo)熱計算結(jié)果的比較，驗證了本文計算程序的準(zhǔn)確性。通過引入一個表示實際接觸面積和球心距離確定的理想接觸面積的比例系數(shù)，可以進行實際球床接觸導(dǎo)熱的工程計算。

［1］Zhang Z，Wu Z，Sun Y，et al．Design aspects of the Chinese modular high-temperature gas-cooled reactor HTR-PM［J］．Nuclear Engineering and Design，2006，236（5）：485-490．

［2］Nishino K，Yamashita S，Torii K．Thermal contact conductance under low applied load in a vacuum environment［J］．Experimental thermal and fluid science，1995，10（2）：258-271．

［3］張品．基于DEM方法的散體顆粒熱傳導(dǎo)模擬 ［D］．湖南：中南大學(xué)，2011．

［4］Nguyen V D，Cogn E C，Guessasma M，et al．Discrete modeling of granular flow with thermal transfer：application to the discharge of silos［J］．Applied thermal engineering，2009，29（8）：1846-1853．

［5］Bu C，Liu D，Chen X，et al．Modeling and Coupling Particle Scale Heat Transfer with DEM through Heat Transfer Mechanisms［J］．Numerical Heat Transfer，Part A：Applications，2013，64（1）：56-71．

［6］Malone K F，Xu B H．Particle-scale simulation of heat transfer in liquid-fluidised beds［J］．Powder Technology，2008，184（2）：189-204．

［7］Zhang H W，Zhou Q，Xing H L，et al．A DEM study on the effective thermal conductivity of granular assemblies［J］．Powder Technology，2011，205（1）：172-183．

［8］Nguyen V，Benhabib K，Marie C，et al．Heat Transfer in a Granular Media Modeled by a Coupled DEM-Finite Difference Method：Application to Fluidized Bed Processes［J］．Procedia Engineering，2012，42：895-904．

［9］Brosh T，Levy A．Modeling of Heat Transfer in Pneumatic Conveyer Using a Combined DEM-CFD Numerical Code［J］．Drying Technology，2010，28（2）：155-164．

［10］Vedachalam V，Virdee D．Discrete Element Modelling Of Granular Snow Particles Using LIGGGHTS［D］．Ph．D Dissertation．Edinburgh Parallel Computing Centre．The University of Edinburgh，UK，2011．

［11］Jacobsen N G，F(xiàn)uhrman D R，F(xiàn)redsoe J．A wave generation toolbox for the open-source CFD library：OpenFoam?［J］．International Journal for Numerical Methods in Fluids，2012，70（9）：1073-1088．

［12］Carslaw H S，Jaeger J C．Conduction of heat in solids［J］．Oxford：Clarendon Press，1959，2nd ed．，1959，1．

［13］Yovanovich M M．Thermal contact resistance across elastically deformed spheres［J］．Journal of Spacecraft and Rockets，1967，4（1）：119-122．

［14］李聰新．高溫氣冷堆球床等效導(dǎo)熱系數(shù)研究 ［D］．清華大學(xué)，2014．

［15］Davis T A，Duff I S．A combined unifrontal/multifrontal meth

od for unsymmetric sparse matrices［N］．ACM Transactions on Mathematical Software（TOMS），1999，25（1）：1-20．

Calculating the Contact Heat Conduction in a High Temperature Gas-cooled Reactor

LI Congxin1，REN Cheng2，XU Chao1?，LIU Yusheng1，ZHANG Pan1
（1．Nuclear and Radiation Safety Center，MEP，Beijing 100082，China；2．Institute of Nuclear and New Energy Technology，Tsinghua University，Beijing，100084，China）

At low temperatures，heat transfer within the core of the high temperature gas-cooled reactor（HTGR）is dominated by contact heat conduction，the prediction of which thereby determining how accurate the heat transfer calculation for any HTGR could be．On the basis of the thermal discrete element method，this paper establishes a model to calculate the contact heat conduction occurring in the pebble bed．In doing so，we first resolve the conduction flux between two ideally-contacting balls via CFD approach，and then develop the detailed procedure of our code which calculates the contact heat conduction for large-scale，randomly-packing pebble bed．Our code has been validated by evaluations against predictions obtained by CFD approach for the contact heat conduction within local structure of pebble bed．Our method enjoys accuracy comparable to CFD simulation of local pebble bed structure，and it is further applicable to large-scale pebble bed simulation in the particle level．

HTGR；pebble bed；contact heat conduction；particle scale；ideal contact；thermal discrete element method

TL331

：A

：1672-5360（2016）04-0052-07

2016-06-24

2016-08-03

中國科學(xué)院先導(dǎo)專項，項目編號 XDA02050500

李聰新 （1984—），男，河北衡水人，工程師，核科學(xué)與工程，現(xiàn)主要從事核安全相關(guān)熱工水力試驗研究工作

?通訊作者：許 超，E-mail：seanwillian＠126．cn