張文文，叢騰龍，田文喜，秋穗正，蘇光輝，謝永誠，蔣 興

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

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

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

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

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

1.1 耦合換熱模型

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

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

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

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

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

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

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

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

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

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

1.2 分布阻力模型

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

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

橫向流動阻力公式如下：

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

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

1.3 Boussinesq假設(shè)

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

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

1.4 湍流模型

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

1.5 數(shù)值計算流程

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

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

2 幾何模型與邊界條件

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

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

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

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

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

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

3 結(jié)果分析

3.1 IRWST側(cè)流動換熱計算結(jié)果

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)語

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

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

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

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

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

［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]薛若軍，鄧程程，彭敏俊.非能動余熱排出熱交換器數(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.