非能動(dòng)余熱排出換熱器換熱能力數(shù)值分析

張文文，叢騰龍，田文喜，秋穗正，蘇光輝，謝永誠(chéng)，蔣 興

（1.西安交通大學(xué) 能源與動(dòng)力工程學(xué)院，陜西 西安 710049；2.上海核工程研究設(shè)計(jì)院，上海 200233）

AP1000核電廠采用非能動(dòng)設(shè)計(jì)來(lái)提高核電廠安全性并簡(jiǎn)化核電廠設(shè)備。對(duì)于非喪失冷卻劑事故（LOCA），非能動(dòng)余熱排出系統(tǒng)（PRHRS）起著至關(guān)重要的作用。PRHRS主要包括非能動(dòng)余熱排出換熱器（PRHR－HX）與安全殼內(nèi)置換料水箱（IRWST）。在正常運(yùn)行階段，PRHR－HX出口管線上的閥門處于關(guān)閉狀態(tài)，當(dāng)安注信號(hào)發(fā)出后，蒸汽發(fā)生器被隔離，PRHRS投入運(yùn)行，一回路冷卻劑通過自然循環(huán)的方式將熱量帶至PRHR－HX。IRWST為衰變熱的熱阱，與PRHR－HX間存在兩種換熱方式：自然對(duì)流換熱與沸騰換熱。在PRHRS運(yùn)行初期，IRWST內(nèi)冷卻劑處于高度過冷狀態(tài)，單相自然對(duì)流換熱是其主要換熱機(jī)制。

PRHR－HX傳熱管數(shù)目龐大，采用直接建模對(duì)計(jì)算能力要求極高，并不適合瞬態(tài)分析。目前，針對(duì)PRHRS的數(shù)值模擬主要是通過減少傳熱管數(shù)來(lái)進(jìn)行模擬，如Strohecker［1］針對(duì)AP600試驗(yàn)臺(tái)架APEX，選取4根傳熱管分別對(duì)單相對(duì)流與過冷沸騰階段進(jìn)行了研究；薛若軍等［2］以 AP1000核電廠內(nèi) PRHR－HX 為研究對(duì)象，通過模型簡(jiǎn)化對(duì)其進(jìn)行了非穩(wěn)態(tài)特性的模擬。為了簡(jiǎn)化幾何建模，多孔介質(zhì)模型已成功應(yīng)用到相同復(fù)雜結(jié)構(gòu)的管殼式換熱器中，Prithiviraj等［3－4］基于多孔介質(zhì)模型研究了管殼式換熱器的流動(dòng)換熱特性；Zhang等［5－6］對(duì)冷凝器內(nèi)的冷凝現(xiàn)象進(jìn)行了準(zhǔn)三維模擬；Cong等［7］對(duì)蒸汽發(fā)生器二次側(cè)三維兩相流場(chǎng)進(jìn)行了穩(wěn)態(tài)分析。

基于上述工作，本文采用多孔介質(zhì)模型對(duì)PRHR－HX進(jìn)行全尺寸建模，并采用一回路側(cè)與IRWST側(cè)耦合換熱的方法對(duì)其非穩(wěn)態(tài)流動(dòng)換熱特性進(jìn)行分析計(jì)算。

1 數(shù)學(xué)物理模型

1.1 耦合換熱模型

傳熱管管壁兩側(cè)耦合換熱計(jì)算可克服采用給定傳熱管外壁熱邊界計(jì)算方法的不足，使得計(jì)算更接近實(shí)際物理現(xiàn)象，可更準(zhǔn)確地分析換熱器的熱量輸出能力。

傳熱管束區(qū)的熱流密度通過下式計(jì)算：

其中：Tp和Ts分別為一回路冷卻劑溫度與IRWST側(cè)冷卻劑溫度；hp為傳熱管內(nèi)壁換熱系數(shù)；hs為傳熱管外壁換熱系數(shù)；k為傳熱管熱導(dǎo)率；do為傳熱管外徑；di為傳熱管內(nèi)徑。

在多孔介質(zhì)模型的應(yīng)用過程中，換熱量以體積熱源的形式添加進(jìn)能量方程中，對(duì)于管束區(qū)，體積熱源計(jì)算公式如下：

其中，ASV為體積面積密度，即單位體積內(nèi)的傳熱管的換熱面積。

對(duì)于傳熱管管壁內(nèi)側(cè)強(qiáng)迫對(duì)流換熱系數(shù)采用Dittus－Bolter公式：

分別對(duì)水平管束區(qū)與豎直管束區(qū)進(jìn)行傳熱管管壁外側(cè)換熱系數(shù)的計(jì)算。對(duì)于水平管束區(qū)，管外單相自然對(duì)流使用了Langmuir關(guān)系式［8］，該關(guān)系式由無(wú)限空間內(nèi)水平單管換熱式發(fā)展而來(lái)，形式如下：

其中，瑞利數(shù)Ra計(jì)算過程中采用管外徑作為特征長(zhǎng)度。

對(duì)于豎直管束區(qū)管外單相自然對(duì)流，使用了 Churchill－Chu關(guān)系式［9］，形式如下：

其中，下部水平管束區(qū)的高度作為該式計(jì)算過程中的特征長(zhǎng)度。該關(guān)系式適用于格拉曉夫數(shù)大于1010的工況，而所有的計(jì)算工況均符合該條件。

1.2 分布阻力模型

本文計(jì)算中考慮的阻力由管束造成，并作為分布阻力添加進(jìn)相應(yīng)區(qū)域的各單元內(nèi)。管束區(qū)阻力的計(jì)算包括軸向流動(dòng)阻力及橫向流動(dòng)阻力，又因C型管的阻力為各向異性，因此在3個(gè)方向上分別添加以各方向速度分量計(jì)算的分布阻力。軸向流動(dòng)阻力由下式計(jì)算：

其中：G為質(zhì)量流量；Δpa為軸向壓降；l為軸向長(zhǎng)度；ρ為冷卻劑密度；de為流通通道當(dāng)量直徑；fa為軸向流動(dòng)摩擦系數(shù)，由 MacAdams公式［10］計(jì)算：

橫向流動(dòng)阻力公式如下：

其中：Δpc為長(zhǎng)度為l的橫向流動(dòng)壓降；fc為橫向流動(dòng)摩擦系數(shù)；vmax為管束間最窄流通截面的流速；pv為棒距。

與軸向流動(dòng)不同的是，采用棒距pv作為特征長(zhǎng)度，特征速度采用管束間最窄流通截面的流速vmax，橫向流動(dòng)阻力系數(shù)采用Grimison經(jīng)驗(yàn)關(guān)系式［11］：

1.3 Boussinesq假設(shè)

在計(jì)算過程中，IRWST被處理為封閉腔室，為滿足連續(xù)性方程，冷卻劑密度需假設(shè)為常數(shù)。因此，采用Boussinesq假設(shè)，在守恒方程求解過程中，密度為定值，而在求解動(dòng)量方程時(shí)，浮力項(xiàng)所包含的密度變化與溫度變化相關(guān)，即：

其中：ρ0為參考密度；T0為參考溫度；β為體積膨脹系數(shù)；g為重力加速度。

1.4 湍流模型

PRHR－HX／IRWST系統(tǒng)瞬態(tài)計(jì)算時(shí)間跨度較大，考慮到計(jì)算能力限制，在建模時(shí)采用粗網(wǎng)格處理方法，因此采用對(duì)網(wǎng)格要求不高的Spalart－Allmaras模型，該模型屬于以Boussinesq假設(shè)為前提的渦粘性模型，通過求解中間變量的輸運(yùn)方程獲得湍流運(yùn)動(dòng)的黏性系數(shù)［12］：

1.5 數(shù)值計(jì)算流程

采用商用CFD軟件FLUENT對(duì)模型進(jìn)行求解。為了得到實(shí)時(shí)的換熱量，分別準(zhǔn)確計(jì)算出傳熱管管壁兩側(cè)流場(chǎng)及溫度場(chǎng)以獲得其換熱系數(shù)，在每步迭代后進(jìn)行換熱量計(jì)算。圖1示出PRHR－HX數(shù)值計(jì)算的流程。

圖1 PRHR－HX數(shù)值計(jì)算流程Fig.1 Calculation procedure of PRHR－HX

2 幾何模型與邊界條件

IRWST內(nèi)包括1個(gè)PRHR－HX與2個(gè)ADS泄壓噴淋器。本文著重分析傳熱管束區(qū)及周邊的流動(dòng)換熱特性，不考慮在事故后期階段才投入使用的泄壓噴淋器。PRHR／IRWST系統(tǒng)結(jié)構(gòu)如圖2所示，忽略支撐板與外殼框架，IRWST為半銅錢結(jié)構(gòu)，傳熱管束區(qū)采用多孔介質(zhì)模型處理，只對(duì)外部輪廓進(jìn)行建模。一回路冷卻劑由上部入口封頭進(jìn)入，經(jīng)水平管束區(qū)熱段、豎直管束區(qū)以及水平管束區(qū)冷段后，由出口封頭流出。

圖2 計(jì)算區(qū)域結(jié)構(gòu)Fig.2 Structure of calculated domain

AP1000PRHR－HX共有689根傳熱管，分為29排29列，水平管束區(qū)縱橫棒距相同，而豎直管束區(qū)縱橫棒距相差一倍，其詳細(xì)結(jié)構(gòu)如圖3所示。為了進(jìn)行耦合換熱計(jì)算，將傳熱管束區(qū)沿一回路冷卻劑流動(dòng)方向劃分為20個(gè)控制體，如圖4所示，并假設(shè)每個(gè)控制體內(nèi)一回路側(cè)網(wǎng)格內(nèi)參數(shù)相同。

圖3 非能動(dòng)余熱排出換熱器結(jié)構(gòu)Fig.3 Structure of PRHR－HX

一回路側(cè)冷卻劑入口溫度與流量采用RELAP5對(duì)AP1000核電廠全廠斷電事故的分析結(jié)果，同時(shí)忽略一回路壓力降低對(duì)冷卻劑物性的影響。為了避免對(duì)水箱上部氣液交界面進(jìn)行建模，將IRWST頂部自由液面處理為零切應(yīng)力絕熱壁面，忽略上表面熱量散失。計(jì)算中，IRWST內(nèi)冷卻劑的初始溫度設(shè)為322K。

圖4 一回路側(cè)控制體劃分Fig.4 Division of control volume in primary loop side

3 結(jié)果分析

3.1 IRWST側(cè)流動(dòng)換熱計(jì)算結(jié)果

伴隨著熱量導(dǎo)入IRWST內(nèi)，冷卻劑受熱膨脹，在浮力的驅(qū)動(dòng)下，熱流體向上流動(dòng)，在壁面附近方向發(fā)生偏轉(zhuǎn)。圖5示出500s與3000s時(shí)IRWST內(nèi)的流線分布。由圖5可見：500s時(shí)自然循環(huán)尚處于啟動(dòng)階段，流體上升到上部區(qū)域后轉(zhuǎn)向向IRWST遠(yuǎn)端流動(dòng)，到達(dá)遠(yuǎn)端后從水箱下部返回到傳熱管束下部，形成豎直方向內(nèi)的自然循環(huán)；3000s時(shí)，流動(dòng)已經(jīng)基本穩(wěn)定，此時(shí)，在水箱上部形成水平渦，即流動(dòng)集中在上部，而豎直方向上的循環(huán)變得較弱。

圖5 IRWST內(nèi)的流線圖Fig.5 Stream line in IRWST

圖6示出不同時(shí)刻對(duì)稱面處的速度矢量分布。由圖6可看出，在500s時(shí)大流速區(qū)主要集中在豎直傳熱管束附近，而在3000s時(shí)上部水平管束區(qū)上方速度較大。這是由于在啟動(dòng)階段，管束區(qū)內(nèi)、外溫差較大，驅(qū)動(dòng)較強(qiáng)，同時(shí)豎直管束區(qū)通道阻力較小，流道長(zhǎng)，因此速度較大，可達(dá)2m／s左右。隨著循環(huán)建立，上部區(qū)域形成水平渦，折返的流體與水平管上升流體匯聚，故此處速度較大。同時(shí)可見，IRWST非管束主體區(qū)域的冷卻劑流動(dòng)性較差，下壁面附近冷卻劑近似停滯。

圖7示出不同時(shí)刻IRWST內(nèi)冷卻劑的溫度分布剖面圖，計(jì)算結(jié)果非常直觀地表現(xiàn)了熱分層現(xiàn)象的形成及發(fā)展。熱分層現(xiàn)象的形成主要與兩個(gè)因素有關(guān)，熱流體在上部區(qū)域的積聚與IRWST內(nèi)主體區(qū)域冷卻劑較小的流動(dòng)性。前者使得熱邊界層在箱內(nèi)上部區(qū)域形成并逐漸向底部區(qū)域發(fā)展，后者則降低了各溫度層間的相互攪混，使得熱分層較為穩(wěn)定。同時(shí)，從圖7中還可看出，最高溫度出現(xiàn)于上部水平管束區(qū)，而箱體底部的“死水區(qū)”的溫度則基本維持在初始溫度水平。

圖6 IRWST對(duì)稱面處的速度矢量分布Fig.6 Velocity vector distribution on symmetry plane in IRWST

圖7 不同時(shí)刻IRWST內(nèi)冷卻劑溫度分布Fig.7 Coolant temperature distribution of IRWST at different time

3.2 耦合換熱計(jì)算結(jié)果

圖8、9分別示出傳熱管內(nèi)壁及外壁對(duì)流換熱系數(shù)的分布。因管內(nèi)冷卻劑溫度沿一回路流動(dòng)方向逐漸降低，造成冷卻劑密度增大而黏性減小，導(dǎo)致管內(nèi)換熱系數(shù)逐漸降低。冷卻劑溫度隨時(shí)間的增加而降低，同樣導(dǎo)致?lián)Q熱系數(shù)降低。由圖8可見，管內(nèi)壁換熱系數(shù)基本處于9～15kW／（m2·K）范圍內(nèi)，變化相對(duì)較小。

圖8 一回路側(cè)沿冷卻劑流動(dòng)方向的換熱系數(shù)Fig.8 Heat transfer coefficient of primary loop side along flow direction of primary fluid

對(duì)于管外對(duì)流換熱，與管內(nèi)對(duì)流換熱相比影響因素較多，管束的排列方式及流動(dòng)方向均對(duì)其產(chǎn)生影響。由圖9可看出，管外壁換熱系數(shù)較小且分布較不均勻，最大值為2.5kW／（m2·K），最小值僅為1kW／（m2·K）左右。豎直管束區(qū)的換熱系數(shù)最大，盡管該區(qū)域的冷卻劑流動(dòng)方式為換熱性較差的平行于管束流動(dòng)，因該區(qū)域長(zhǎng)通道及低阻力的特點(diǎn)，使得其速度遠(yuǎn)大于水平管束區(qū)，進(jìn)而換熱系數(shù)也高于橫掠管束流動(dòng)的水平管束區(qū)。

圖9 IRWST側(cè)沿一回路冷卻劑流動(dòng)方向的換熱系數(shù)Fig.9 Heat transfer coefficient of IRWST side along flow direction of primary fluid

圖10示出一回路側(cè)冷卻劑不同時(shí)刻的溫度分布。由圖10可見，一回路側(cè)入出口溫降由500s時(shí)的130K減小至5000s時(shí)110K。水平管束區(qū)熱段及豎直管束區(qū)上半部各部分溫度隨時(shí)間的增加而降低，而豎直管束下半部與水平管束區(qū)冷段溫度隨時(shí)間的增加反而升高。產(chǎn)生這種現(xiàn)象的原因主要是由于上游換熱能力的降低使得一回路側(cè)冷卻劑的沿程溫降減小所致。

圖10 一回路側(cè)沿冷卻劑流動(dòng)方向的平均溫度分布Fig.10 Average temperature distribution of primary loop side along flow direction of primary fluid

圖11示出管束區(qū)域不同時(shí)刻的能量源項(xiàng)分布。換熱量的大小主要與兩側(cè)流體溫差及熱阻相關(guān)，其中，分布極不均勻的管壁外側(cè)換熱系數(shù)是傳熱過程的主要熱阻。但因豎直管束區(qū)的傳熱管排列稀疏，單位體積有效換熱面積較小，故其能量源項(xiàng)分布與圖9所示的換熱系數(shù)分布并不相同。

圖11 一回路側(cè)沿冷卻劑流動(dòng)方向的能量源項(xiàng)分布Fig.11 Energy source distribution of primary loop side along flow direction of primary fluid

對(duì)各部分換熱量求和即得到換熱器總換熱量隨時(shí)間的變化趨勢(shì)，如圖12所示。在前5000s內(nèi)，PRHR－HX的熱負(fù)荷由150MW 下降至約80MW。一回路側(cè)入口溫度及流量降低，IRWST內(nèi)冷卻劑溫度升高，使得兩側(cè)溫差減小，導(dǎo)致?lián)Q熱器換熱能力下降。

圖12 換熱器總換熱量隨時(shí)間的變化Fig.12 Total heat transfer of heat exchanger vs time

4 結(jié)語(yǔ)

本文采用多孔介質(zhì)模型，對(duì)AP1000核電廠PRHR－HX進(jìn)行了數(shù)值模擬，得出如下結(jié)論。

1）采用多孔介質(zhì)模型模擬PRHR－HX，相比采用常規(guī)CFD網(wǎng)格的方法，大幅節(jié)省了幾何建模時(shí)間與計(jì)算量，同時(shí)可獲得精度較高的計(jì)算結(jié)果；

2）系統(tǒng)投入運(yùn)行后，IRWST內(nèi)冷卻劑的流動(dòng)主要發(fā)生在傳熱管束區(qū)及箱體上部區(qū)域，同時(shí)熱流體在上部區(qū)域的聚集造成了明顯的熱分層現(xiàn)象；

3）傳熱管管壁外側(cè)換熱系數(shù)要小于一回路側(cè)，是換熱器工作的主要熱阻，且分布較不均勻，換熱量計(jì)算的正確與否取決于該換熱系數(shù)計(jì)算的準(zhǔn)確性；

4）隨著時(shí)間的增加，一回路自然循環(huán)能力的降低導(dǎo)致入口流量的下降，以及換料水箱內(nèi)冷卻劑溫度升高使得傳熱管兩側(cè)溫差降低，其自然對(duì)流強(qiáng)度下降，造成換熱器總換熱量逐漸減小。

［1]STROHECKER M F.Investigation of the IRWST flow patterns during a simulated station blackout experiment on the OSU APEX facility［D］.USA：Oregon State University，1998.

［2]薛若軍，鄧程程，彭敏俊.非能動(dòng)余熱排出熱交換器數(shù)值模擬［J］.原子能科學(xué)技術(shù)，2010，44（4）：429－435.XUE Ruojun，DENG Chengcheng，PENG Minjun.Numerical simulation of passive residual heat removal heat exchanger［J］.Atomic EnergyScience and Technology，2010，44（4）：429－435（in Chinese）.

［3]PRITHIVIRAJ M，ANDREWS M.Three dimensional numerical simulation of shell－and－tube heat exchangers，Ⅰ：Foundation and fluid mechanics［J］.Numer Heat Transfer：Part A Appl，1998，33（8）：799－816.

［4]PRITHIVIRAJ M， ANDREWS M. Threedimensional numerical simulation of shell－andtube heat exchangers，Ⅱ：Heat transfer［J］.Numer Heat Transfer：Part A Appl，1998，33（8）：817－828.

［5]ZHANG C，BOKIL A.A quasi－three－dimensional approach to simulate the two－phase fluid flow and heat transfer in condensers［J］.International Journal of Heat and Mass Transfer，1997，40（10）：3537－3546.

［6]ZHANG C.Numerical modeling using aquasithree－dimensional procedure for large power plant condensers［J］.Journal of Heat Transfer，1994，116（1）：180－188.

［7]CONG Tenglong，TIAN Wenxi，SU Guanghui，et al.Three－dimensional study on steady thermohydraulics characteristics in secondary side of steam generator［J］.Progress in Nuclear Energy，2014，70（1）：188－198.

［8]YONOMOTO T，KUKITA Y，SCHULTZ R R.Heat transfer analysis of the passive residual heat removal system in ROSA／AP600experiments［J］.Nuclear Technology，1998，124（1）：18－30.

［9]CHURCHILL S W，CHU H S.Correlating equations for laminar and turbulent free convection from a vertical plate［J］.International Journal of Heat and Mass Transfer，1975，18（11）：1323－1329.

［10］MacADAMS W H.Heat transmission［M］.New York，US：McGraw Hill，1954.

［11］GRIMISON E.Correlation and utilization of new data on flow resistance and heat transfer for cross flow of gases over tube banks［J］.Trans ASME，1937，59：583－594.

［12］SPALART P，ALLMARAS S.A one－equation turbulence model for aerodynamic flows，Technical report AIAA－92－0439［R］.US：American Institute of Aeronautics and Astronautics，1992.