中間熱交換器建模仿真研究

張厚明，段天英，王 剛，劉 勇，熊文彬，劉兆陽

（1．國家核安全局華北核與輻射安全監(jiān)督站，北京100191；2．中國原子能科學(xué)研究院，北京102413；3．國家核安全局核與輻射安全中心，北京100082）

中間熱交換器是液態(tài)金屬冷卻快中子增殖堆（以下簡稱“快堆”）熱傳輸系統(tǒng)中隔離放射性的一回路載熱劑與非放射性二回路載熱劑的屏障，同時也是將堆芯釋熱傳遞給蒸汽發(fā)生器的重要設(shè)備。

對于快堆核電廠熱工水力、儀控系統(tǒng)的設(shè)計來講，中間熱交換器的建模及仿真研究可以為之提供設(shè)計依據(jù)。同時，可以將建模及仿真研究的成果用于設(shè)計工程仿真機、全范圍仿真機來為操作員提供培訓(xùn)。另外，仿真模型在高速計算機上運行，通過實時采集快堆核電廠數(shù)據(jù)提前計算出各個參數(shù)的運行趨勢，從而為操縱員提供技術(shù)支持。因此，對于上述三個方面，中間熱交換器建模及仿真研究都具有重要意義。本文給出一種計算迅速、精度較高、可以用于控制系統(tǒng)分析、實時計算的中間熱交換器仿真模型。

1 中間熱交換器建模簡介

早期實驗、原型快堆中間熱交換器類型、結(jié)構(gòu)形式較多［1，2］，例如：BN－350采用豎直U形傳熱管束置于水平矩形筒體內(nèi)；BR－5及Grand Quevilly French Test Station采用U形傳熱管束置于U形筒體內(nèi)；Fermi及FFTF采用直傳熱管置于豎直圓柱形筒體內(nèi)，采用可移動的封頭來補償熱膨脹效應(yīng)；Hallam Nuclear Power Facility采用直傳熱管，殼側(cè)有膨脹節(jié)來補償熱膨脹效應(yīng)。

示范、商用大型快堆，例如Super Phénix、CRBRP（未建造）、BN－800，均采用直換熱管、逆流管殼式中間熱交換器。二回路鈉均從頂部進入，通過中心下降管流到底部下腔室，然后流向反轉(zhuǎn)，向上流過換熱管束。一回路鈉均從上部進入，流經(jīng)傳熱管間的殼側(cè)，從底部流出，如圖1所示［3］。

因此，本文選取某個直換熱管、逆流殼式中間熱交換器作為研究對象。

描述中間熱交換器三維熱工水力的方程組非常復(fù)雜。由于動量、能量方程的相互耦合，并且，這些方程的非線性使得方程組很難求解。然而，即使在自然循環(huán)工況下，三維計算結(jié)果與一維計算結(jié)果也相差不大［4］，因此以往的研究中，大多假設(shè)：

圖1 CRBRP中間熱交換器概圖Fig．1 Sketch of CRBRP IHXs'

（a）用一根傳熱管代替?zhèn)鳠峁苁?/p>

（b）流體物性參數(shù)沿徑向保持不變；

（c）忽略載熱劑、管壁的軸向?qū)幔?/p>

（d）忽略中間熱交換器殼體的散熱。

在上述這些假設(shè)下，Gunby［5］采用“節(jié)距固定”方法將中間熱交換器的傳熱管劃分為10個等長節(jié)段，在每個節(jié)段上建立了能量守恒微分方程組，采用不同的差分來求解，研究了流量、入口溫度變化等各種工況下中間熱交換器的溫度的分布及變化。研究結(jié)果表明：對一、二回路不平衡流動（10∶1），流量大的一側(cè)其計算結(jié)果差別較大；并且，即使出口溫度計算結(jié)果較為可信，各種差分方法計算出的內(nèi)部溫度分布也均不正確；在要求出口溫度和內(nèi)部溫度分布計算結(jié)果同時可靠時，大擾動情況下，采用這種“節(jié)距固定”的方法，劃分節(jié)段數(shù)目較少時，各種差分方法均難以獲得正確的計算結(jié)果。

Brukx［6］針對Gunby方法中出現(xiàn)的內(nèi)部溫度分布不正確問題，給出了一種將能量守恒微分方程組在空間上離散的差分計算方法，計算結(jié)果較好，但這種方法中無法實時計算溫度變化，在控制系統(tǒng)設(shè)計，例如使用Matlab／Simulink軟件對控制系統(tǒng)進行設(shè)計中，諸多不便。

段天英［7，8］根據(jù)能量守恒方程，在小擾動情況下推導(dǎo)出中間熱交換器出口溫度隨流量、入口溫度等參數(shù)變化的傳遞函數(shù)，并將傳遞函數(shù)進行修正，計算結(jié)果較為精確，這種方法可以用于擾動不大的控制系統(tǒng)設(shè)計中；但是，這種方法無法計算各種瞬態(tài)工況下中間熱交換器內(nèi)部的溫度變化，且不可以用于擾動較大的控制系統(tǒng)設(shè)計。

Tzanos［9］給出兩種中間熱交換器模型。第一種為“節(jié)距固定”方法：將傳熱管均分為145個節(jié)段，在每個節(jié)段上采用能量守恒方程進行描述，這樣獲得描述中間熱交換器模型的常微分方程組。由于早期計算機運算速度較慢，Tzanos給出第二種“節(jié)距可變”方法：將中間熱交換器按照等能量劃分為10～15個節(jié)段，采用Leibnitz方程將每個節(jié)段上的能量守恒方程偏微分方程組轉(zhuǎn)化為溫度導(dǎo)數(shù)與節(jié)段長度導(dǎo)數(shù)有關(guān)的常微分方程組。由于出口溫度變化與長度變化有關(guān)，“節(jié)距可變”方法使得節(jié)段平均溫度為節(jié)段出入口溫度的算術(shù)平均值這一假設(shè)在節(jié)段數(shù)目較少的情況下比“節(jié)距固定”的方法更為合理，于是，采用這種階段長度可變的方法可以大大減少常微分方程組維數(shù)和計算機機時。但是，這種“節(jié)段長度可變”的方法與“節(jié)距固定”方法的計算結(jié)果仍有一定差別。這種方法在EBR－Ⅱ的TEST－8A試驗［9，10］中得到驗證，最大計算誤差在6K左右。

綜上所述，在以往的對中間熱交換器建模的各種方法中，采用“節(jié)距固定”的方法來建立中間熱交換器模型，若取得較好的計算結(jié)果，瞬態(tài)過程中中間熱交換器內(nèi)溫度可信，不會出現(xiàn)數(shù)值不穩(wěn)定，那么節(jié)段數(shù)目的選取應(yīng)至少為100個左右［9］。本文給出150節(jié)段、適于Matlab／Simulink軟件中計算的“節(jié)距固定”模型。本模型建立過程中首先提出兩種穩(wěn)態(tài)計算方法，比較兩種方法計算的穩(wěn)態(tài)結(jié)果，然后給出一種瞬態(tài)模型，并將穩(wěn)態(tài)計算的結(jié)果作為瞬態(tài)計算的參數(shù)，對二回路鈉流量指數(shù)衰減的工況進行了計算，以下分別介紹穩(wěn)態(tài)、瞬態(tài)計算模型。

2 中間熱交換器穩(wěn)態(tài)計算方法

文獻［11］在假定穩(wěn)態(tài)情況下整個中間熱交換器內(nèi)鈉、換熱管材料的比熱容及換熱系數(shù)均不變，推導(dǎo)出沿高度方向溫度呈指數(shù)分布，這種計算方法較為粗糙，計算結(jié)果誤差較大。但本文在程序設(shè)計中將這種方法計算出的結(jié)果作為粗略估計值，從而可以減少穩(wěn)態(tài)計算的循環(huán)次數(shù)。

本文給出的兩種穩(wěn)態(tài)工況下的計算方法，將中間熱交換器傳熱管分為n個節(jié)段，如圖2所示。

圖2 中間熱交換器第i個節(jié)段Fig．2 Segment i in IHXs

穩(wěn)態(tài)工況下，中間熱交換器在節(jié)段i上的能量守恒方程為：

式（1）中，

其中，W為質(zhì)量流量，C為比熱容，U為傳熱系數(shù)，l為節(jié)段長度，T為溫度，λ為熱導(dǎo)率，a為一回路流道截面積，p為節(jié)距；角標p為一回路，m為管壁，s為二回路，in為管內(nèi)，out為管外。

根據(jù)式（1）采用包括：等節(jié)段長度、等換熱量兩種方法來劃分節(jié)段數(shù)目并編制程序，具體算法、循環(huán)語句如圖3所示。根據(jù)這兩種方法，分別調(diào)整α、β、δ、ε的值，可以獲得精度較高的計算結(jié)果，經(jīng)過計算兩種算法出的結(jié)果相差很小，如圖4所示。與實際值相比，所計算出的值絕對誤差誤差在3K以內(nèi)。

圖3 穩(wěn)態(tài)計算方法Fig．3 Steady state algorithms

圖4 穩(wěn)態(tài)工況下的計算結(jié)果Fig．4 Calculated result with steady state algorithms

由計算結(jié)果可知：管內(nèi)熱阻不超過總熱阻的15%～20%，殼側(cè)的傳熱系數(shù)明顯小于管側(cè)，其熱阻約為總熱阻的30%～40%，管壁熱阻約為總熱阻的40%～55%。因此在設(shè)計中，進一步增加管內(nèi)流速不會大幅降低總熱阻；增加殼側(cè)的流速可以有效降低總熱阻，減小管壁厚度可以大大降低總熱阻。這與文獻［1，2］得出的結(jié)論是一致的。

鈉、管壁材料的物性參數(shù)來自文獻［1，11］；一、二回路的傳熱系數(shù)計算公式來自文獻［9］。

3 中間熱交換器瞬態(tài)模型

用來描述中間熱交換器瞬態(tài)變化較好的模型應(yīng)包括一、二回路入口流量、溫度變化帶來中間熱交換器內(nèi)部、出口溫度變化的數(shù)學(xué)模型，同時，應(yīng)選取一種計算迅速、易于使用、誤差小的計算方法。根據(jù)上述討論，描述中間熱交換器的能量守恒偏微分方程組為方程（3）的形式。

其中，密度ρ、比熱C、傳熱系數(shù)U，均為溫度的函數(shù)。

在每個節(jié)段上，將方程（3）離散為差分方程（4）：

對方程（4）這種差分方法進行了穩(wěn)定性分析，將方程（4）展開，并轉(zhuǎn)化為方程（5）的形式。

設(shè)3n－3維向量T＝（Tp，Tm，Ts），也即：

那么差分方程（5）可以寫成方程（6）的形式，簡寫為MT′＝NT，其中A－G均為系數(shù)矩陣。

根據(jù)差分方程組（4）依次建立每個節(jié)段上的模型，如圖5所示。并根據(jù)上述文獻及穩(wěn)態(tài)計算結(jié)果計算出每個節(jié)段上的物性參數(shù)與各微分方程組的系數(shù)。

圖5 第i個節(jié)段上的瞬態(tài)計算模型Fig．5 Transient state model of segment i in IHXs

通過本文中的瞬態(tài)模型中變量取值，帶入上述矩陣方程分析，在各種工況下增長矩陣M－1N的譜半徑始終小于1，符合Von－Neumann條件，故差分格式是穩(wěn)定的。

根據(jù)計算出來的每個微分方程組的系數(shù)，可得方程組的最大剛性比達106左右，對于這種剛性常微分方程組，采用ode15s方法求解，選取相對精度為10－4。

采用這種模型和算法，在其他參數(shù)均不變，對于為了便于比較，使輸出量趨于同一方向，輸入量分別選取100%功率工況下：一次側(cè)流量階躍降低5%、二次側(cè)流量階躍上升5%、一次側(cè)入口溫度階躍降低5%及二次側(cè)入口溫度階躍降低5%，對一、二次側(cè)出口溫度影響分別如圖6、圖7所示。

圖6 一回路出口溫度Fig．6 Temperature at primary outlet

圖7 二回路出口溫度Fig．7 Temperature at secondary outlet

計算結(jié)果表明：

（1）流通延遲

一次側(cè)入口溫度對其出口溫度影響有約6s的延遲、二次側(cè)（中間回路）入口溫度對其出口溫度有約2．9s的延遲，這與文獻［12］中相應(yīng)的分別有5．8s、2．9s延遲的結(jié)論一致的；

（2）一次側(cè)出口溫度

二次側(cè)入口溫度對一次側(cè)出口溫度影響較之其他3個參數(shù)的影響大。二次側(cè)入口溫度階躍降低后，一次側(cè)出口溫度迅速降低，經(jīng)過約80s穩(wěn)定在最終值附近。二次側(cè)入口溫度降低約29．2K，一次側(cè)出口溫度降低約23．5K，這與穩(wěn)態(tài)計算結(jié)果相差較小，這是因為，如圖5，在瞬態(tài)模型建立中將鈉溫度、流量對傳熱系數(shù)的影響加入到相應(yīng)的程序中。

圖8 一回路流量衰減工況計算結(jié)果Fig．8 Simulation result of loss of primary flow transient

一、二次側(cè)鈉流量對一次側(cè)鈉出口溫度的影響均較小，約5K，且動態(tài)響應(yīng)過程也類似，經(jīng)過約90s一次側(cè)出口鈉溫穩(wěn)定在最終值附近；一次側(cè)入口溫度對其出口溫度影響較之鈉流量的影響稍大，動態(tài)響應(yīng)過程也稍長一些，經(jīng)過約100s穩(wěn)定在最終值附近；

（3）二次側(cè)（中間回路）出口溫度

一次側(cè)入口溫度對二次側(cè)（中間回路）出口溫度影響較之其他3個參數(shù)的影響大。一次側(cè)入口溫度階躍降低后，二次側(cè)出口溫度迅速降低，經(jīng)過約90s穩(wěn)定在最終值附近。一次側(cè)入口溫度降低約39．4K，二次側(cè)出口溫度降低約36．5K，這與穩(wěn)態(tài)計算結(jié)果相差很小。一、二次側(cè)鈉流量和對二次側(cè)出口溫度影響均較小，并且過程較為緩慢、平滑，經(jīng)過約90s穩(wěn)定在最終值附近。

另外，在一回路鈉流量以時間常數(shù)為15s指數(shù)衰減的工況下進行計算，給出了一、二次側(cè)溫度的25s、50s、75s時的計算結(jié)果，如圖8所示，其中各組曲線中溫度較高的為一次側(cè)、較低的為二次側(cè)。

4 結(jié)論及展望

本文使用Matlab／Simulink軟件建立的中間熱交換器模型，兩種穩(wěn)態(tài)計算方法計算出的符合性很好，與實際差別較小。瞬態(tài)模型計算迅速、數(shù)值穩(wěn)定。與以往的中間熱交換器建模方法相比，本文給出的兩種穩(wěn)態(tài)計算模型考慮了各物性參數(shù)變化對計算結(jié)果造成的影響。瞬態(tài)計算模型考慮了流量及鈉溫變化對換熱系數(shù)造成的影響，適用于大擾動的計算，例如主泵墮轉(zhuǎn)流量衰減。同時，瞬態(tài)模型采用Simulink模塊庫中的模塊建模，從而避開了具體算法程序的編寫，可以方便地對模型進行修改。因此本模型可用于中間熱交換器的熱工水力設(shè)計以及快堆核電廠的控制系統(tǒng)設(shè)計中。

未來將在中國實驗快堆（CEFR）的試驗中對本模型進行驗證，或者在其他公開文獻中獲得數(shù)據(jù)對其進行驗證。

同時，將進一步對中間熱交換器建模仿真進行研究，從一維向三維熱工、水力耦合模型擴展，從正常工況向自然循環(huán)工況擴展，從而獲取一種更為精確的模型。

［1］ 蘇著亭，葉長源，閻鳳文，等．鈉冷快增殖堆［M］．北京：原子能出版社，1991：332，361－362．

［2］ Tang Y S，Coffield Jr R D，Merkley R A．Thermal Analysis of Liquid Metal Fast Breeder Reactors［M］．ISBN 0894480111．La Grange Park，Illinois，USA：American Nuclear Society，1978．319－332，369－373．

［3］ Devlin R W，Bresnahan J D．FFTF and CRBRP intermediate heat exchanger design，testing，and fabrication［C］．Seminar on LMFBR components，Canoga Park，CA，USA，1977

［4］ CheungA C，Kao T T，Cho S M，et al．A Comparison of Between DEMO 1－D and the COMMIX 3－D Analysis of the CRBRP－IHX Under a Natural Circulation Event［C］．21st ASME／AIChE National Heat Transfer Conference．Seattle，Washington，USA，1983．

［5］ GunbyA L．Intermediate Heat Exchanger Modeling for FFTF Simulation［R］．AEC Research ＆Development Report of Battle Northwest Library．BNWL－1367，1970．

［6］ Brukx J F L M．CSDT Simulation Model for Prediction of IHX Transient Behavior Using S／360CSMP［D］．Delft，Netherlands：Delft University of Technology，1974．

［7］ 段天英．中國實驗快堆控制系統(tǒng)的仿真研究［D］．中國原子能科學(xué)研究院博士學(xué)位論文，1999．

［8］ 段天英．熱交換器仿真模型的修正［J］．中國試驗快堆年報，2000：5．

［9］ Constantine P，Tzanos．Liquid－Metal Fast Reactor Breeder Beactor Intermediate Heat Exchanger Transient Modeling for Faster than Real－Time Analysis［J］．Nuclear Technology，1987，76（3）：337－351．

［10］ Constantine P，Tzanos．Validation of an Intermediate Heat Exchanger Model for Real Time Analysis［C］．American Nuclear Society and Atomic Industrial Forum Joint Meeting．Washington DC，USA，1986．

［11］ 居懷明，徐元輝，李懷萱．載熱質(zhì)熱物性計算程序及數(shù)據(jù)手冊［M］．北京：原子能出版社，1990：92－99．

［12］ 中國實驗快堆最終安全分析報告［R］．中國原子能科學(xué)研究院，2008．