，，，?，，

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2 計(jì)算程序開(kāi)發(fā)

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

2.l前處理

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

圖l 每個(gè)球和壁面的存儲(chǔ)信息Fig.l Stored information of each pebble and wall

2.2 方程求解

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

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

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

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

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

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

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

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

2.3 后處理

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

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

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

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

3 模型驗(yàn)證

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

3.l一維驗(yàn)證

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

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

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

3.2三維驗(yàn)證

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

［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

中國(guó)科學(xué)院先導(dǎo)專項(xiàng)，項(xiàng)目編號(hào) XDA02050500

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

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