張厚明，段天英，王 剛，劉 勇，熊文彬，劉兆陽(yáng)

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

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

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

1 中間熱交換器建模簡(jiǎn)介

早期實(shí)驗(yàn)、原型快堆中間熱交換器類(lèi)型、結(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)，采用可移動(dòng)的封頭來(lái)補(bǔ)償熱膨脹效應(yīng)；Hallam Nuclear Power Facility采用直傳熱管，殼側(cè)有膨脹節(jié)來(lái)補(bǔ)償熱膨脹效應(yīng)。

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

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

描述中間熱交換器三維熱工水力的方程組非常復(fù)雜。由于動(dòng)量、能量方程的相互耦合，并且，這些方程的非線性使得方程組很難求解。然而，即使在自然循環(huán)工況下，三維計(jì)算結(jié)果與一維計(jì)算結(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個(gè)等長(zhǎng)節(jié)段，在每個(gè)節(jié)段上建立了能量守恒微分方程組，采用不同的差分來(lái)求解，研究了流量、入口溫度變化等各種工況下中間熱交換器的溫度的分布及變化。研究結(jié)果表明：對(duì)一、二回路不平衡流動(dòng)（10∶1），流量大的一側(cè)其計(jì)算結(jié)果差別較大；并且，即使出口溫度計(jì)算結(jié)果較為可信，各種差分方法計(jì)算出的內(nèi)部溫度分布也均不正確；在要求出口溫度和內(nèi)部溫度分布計(jì)算結(jié)果同時(shí)可靠時(shí)，大擾動(dòng)情況下，采用這種“節(jié)距固定”的方法，劃分節(jié)段數(shù)目較少時(shí)，各種差分方法均難以獲得正確的計(jì)算結(jié)果。

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

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

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

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

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

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

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

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

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

式（1）中，

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

計(jì)算結(jié)果表明：

（1）流通延遲

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

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

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

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

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

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

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

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

4 結(jié)論及展望

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

未來(lái)將在中國(guó)實(shí)驗(yàn)快堆（CEFR）的試驗(yàn)中對(duì)本模型進(jìn)行驗(yàn)證，或者在其他公開(kāi)文獻(xiàn)中獲得數(shù)據(jù)對(duì)其進(jìn)行驗(yàn)證。

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

［1］ 蘇著亭，葉長(zhǎng)源，閻鳳文，等．鈉冷快增殖堆［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］ 段天英．中國(guó)實(shí)驗(yàn)快堆控制系統(tǒng)的仿真研究［D］．中國(guó)原子能科學(xué)研究院博士學(xué)位論文，1999．

［8］ 段天英．熱交換器仿真模型的修正［J］．中國(guó)試驗(yàn)快堆年報(bào)，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ì)熱物性計(jì)算程序及數(shù)據(jù)手冊(cè)［M］．北京：原子能出版社，1990：92－99．

［12］ 中國(guó)實(shí)驗(yàn)快堆最終安全分析報(bào)告［R］．中國(guó)原子能科學(xué)研究院，2008．