非能動余熱排出換熱器換熱能力數值分析

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

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

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

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

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

1 數學物理模型

1.1 耦合換熱模型

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

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

其中：Tp和Ts分別為一回路冷卻劑溫度與IRWST側冷卻劑溫度；hp為傳熱管內壁換熱系數；hs為傳熱管外壁換熱系數；k為傳熱管熱導率；do為傳熱管外徑；di為傳熱管內徑。

在多孔介質模型的應用過程中，換熱量以體積熱源的形式添加進能量方程中，對于管束區，體積熱源計算公式如下：

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

對于傳熱管管壁內側強迫對流換熱系數采用Dittus－Bolter公式：

分別對水平管束區與豎直管束區進行傳熱管管壁外側換熱系數的計算。對于水平管束區，管外單相自然對流使用了Langmuir關系式［8］，該關系式由無限空間內水平單管換熱式發展而來，形式如下：

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

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

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

1.2 分布阻力模型

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

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

橫向流動阻力公式如下：

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

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

1.3 Boussinesq假設

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

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

1.4 湍流模型

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

1.5 數值計算流程

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

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

2 幾何模型與邊界條件

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

圖2 計算區域結構Fig.2 Structure of calculated domain

AP1000PRHR－HX共有689根傳熱管，分為29排29列，水平管束區縱橫棒距相同，而豎直管束區縱橫棒距相差一倍，其詳細結構如圖3所示。為了進行耦合換熱計算，將傳熱管束區沿一回路冷卻劑流動方向劃分為20個控制體，如圖4所示，并假設每個控制體內一回路側網格內參數相同。

圖3 非能動余熱排出換熱器結構Fig.3 Structure of PRHR－HX

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

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

3 結果分析

3.1 IRWST側流動換熱計算結果

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

圖5 IRWST內的流線圖Fig.5 Stream line in IRWST

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

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

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

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

3.2 耦合換熱計算結果

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

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

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

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

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

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

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

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

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

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

4 結語

本文采用多孔介質模型，對AP1000核電廠PRHR－HX進行了數值模擬，得出如下結論。

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

2）系統投入運行后，IRWST內冷卻劑的流動主要發生在傳熱管束區及箱體上部區域，同時熱流體在上部區域的聚集造成了明顯的熱分層現象；

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

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

［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]薛若軍，鄧程程，彭敏俊.非能動余熱排出熱交換器數值模擬［J］.原子能科學技術，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.