A specialized code for operation transient analysis and its application in fluoride salt-cooled high temperature reactors

A specialized code for operation transient analysis and its application in fluoride salt-cooled high temperature reactors 1 2 first-author 1 2 1 1 1 corresp zouyang@sinap.ac.cn; zouyang@sinap.ac.cn; zouyang@sinap.ac.cn 1 corresp Corresponding author E-mail address: xuhongjie@sinap.ac.cn; Corresponding author E-mail address: xuhongjie@sinap.ac.cn; xuhongjie@sinap.ac.cn Shanghai Institute of Applied Physics, CAS, Shanghai, 201800, China Shanghai Institute of Applied Physics, CAS, Shanghai, 201800, China University of Chinese Academy of Sciences, Beijing 100049, China University of Chinese Academy of Sciences, Beijing 100049, China Fluoride salt-cooled High-temperature Reactors (FHRs) includemany attractive features, such as high temperature, large heat capacity, low pressureand strong inherent safety. Transient characteristics of FHR areparticularly important for evaluating its operation performance.Thus, a specialized code OCFHR (Operation and Control analysis code of FHR)isused to studyan experimental FHR’s operation behaviors. The geometric modeling of OCFHR is based on one-dimensional lumped parameter method,and some simplifications are taken into consideration during simulation due to the existence of complex structures such as pebble bed, intermediate heat exchanger (IHX),air radiator (AR) and multiply channels. Apoint neutron kinetics modelis developed and neutronphysics calculation is neededto provide some key inputs including axial power density distribution, reactivity coefficients and parameters about delayed neutron precursors. For analyzing the operational performance, five disturbed transients are simulated, involvingreactivity step insertion, variations of coolant mass flow rate of primaryloop andintermediate loop, adjustment of air inlet temperature,and mass flow rate of air-cooling system. Simulation results indicate that inherent self-stabilityof FHRrestrains severe consequencesunderabove transients, and some dynamic featuresareobserved, such as large negative temperature feedbacks, remarkable thermal inertia, and high response delay. FHR Simulation Pebble bed Transient analysis; 1 Introduction 1 Fluoride salt-cooled high-temperature reactors (FHRs),which are notable for thecoated particle fuel and fluoride salt coolant, provide outstanding priority in electricity generation and process heat applications.Its concept prototypewas developed by Oak Ridge National Laboratory (ORNL), Sandia National Laboratory (SNL) and the University of California at Berkeley (UCB) in 2001–2003, assimilating technical advantages of Molten Salt Reactor (MSR) and High Temperature Reactor (HTR)[ 1 1 ]. Chinese FHR research startedin2012 by TMSR (Thorium Molten Salt Reactor) Group of Chinese Academy of Sciences, with aconceptual design offluoride salt-cooled high-temperature experimental reactor (FHER)[ 2 2 ], as shown in Fig.1 Fig.1 . Molten salt and particle fuel affectgreatly FHRs’ transient performance. The coolant in primary loop is FLiBe(2LiF-BeF 2 ) which has low neutron absorption cross-section, large heat capacity (~2400J·kg −1 ·°C −1 ), and a wide range of liquid-phase temperature (450°C–1400°C), providing an abundant temperature safety margin [ 3 3 ]. It allows FHRs to increase average power density, reducemaximum pressure tolerance of material and enhance spatial compactness. The fuel element consists of a fuel core (either in plate type or spherical type) and several shielding materials ( Fig.2 Fig.2 ). The fuel core is mainly made up of tristructural-isotropic (TRISO), offering several hundred degrees safety margin comparing to rated operational temperature and leaving little possibility for fission products to leak out [ 4 4 ]. In this paper, spherical pebble type fuel element is considered. Understanding of transient performances is significant in FHR design, asthe large heat capacity of salt may extend transient response time and the multilayered structure of fuel element means non-negligible temperature variation delay. Moreover, the core, pipelines and other components also restrain thermal transfer speed as it costs no less than 70s for coolant to flow around the main loop, postponing the process of transients when disturbances occurs. Therefore, it is quite necessary to study FHRs’ performances under some typical operation disturbances. Analysis of transient operation has been done experimentallyor theoretically. Chen used MATLAB platform to analyze the transient performances of an FHR coupled with a helium Brayton power cycle [ 5 5 ]. Li simulated some kinds of disturbed transients’ influence on TMSR-SF1, with the RELAP5 codes,to summarize its characteristics[ 6 6 ]. Lv developed a computer code to evaluate the direct reactor auxiliary cooling system (DRACS) [ 7 7 ]. Ablay used a modeling and control approach to study the power conversion performance of advanced nuclear plants with gas turbines [ 8 8 ]. In this paper, the OCFHR (Operation and Control analysis code of FHR) is introduced in detail, including theequations, simplifications andequivalences of complex structures, and verification of calculation models.Five disturbed transients are simulated, involvingreactivity step insertion, andchanges in mass flow rate of coolant ofthree loops and in air inlet temperature of air-cooling system.The simulation resultsof theentire FHR systemare compared with thosefromRELAP5 code.Additional work will be performed to study operation control logic and power conversion[ 9 9 ]. 2 Method 1 2.1. Algorithm 2 CFHRis implemented by fourth order Runge-Kutta(R-K) method for numerical solution of differential equations. This algorithm is an explicit numerical difference scheme, and it can construct arbitrary high order calculation method [ 10 10 ]. Besides, Gear algorithm, a variable step length algorithm, is programmed in comparison to R-K method on computation speed and accuracy. 2.2. Simplification 2 The FHER system consists of the pebble bed, cavities, pipes, pumps, heat exchangers, load etc . The simplification of FHR system in OCFHR is based on lumped parameter method, which axially divides the object structure into several independent nodes, each representing segmental geometrical and physical parameters. Differential equationsformaterial properties and heat transfer, and conservation equations, are developed within a node. Every component model in OCFHR is labeled by a series of numbers, each denoting a node. Whenan additional structure is added in, so will the label number be inserted into the array. Hence, the whole system is simplified under this modularization thinking,which brings in significant progress in modeling speed [ 11 11 ]. Fig.3 Fig.3 shows the node diagram. The primary coolant heated from the core flows along the upper cavity and the hot leg to the inlet of IHX. Then the cooled salt reaches to the inlet of downcomer and lower cavity, eventually returns back to the core. Positions E and D areinlet and outlet of IHX cooling side, while F and G areinlet and outlet of SHX heating side. Fig.4 Fig.4 represents simplifications of the pebble bed. Relevant parameters of FHER are givenin Table 1 Table 1 . 1 1 Parameters 1 1 Values 1 1 Fuel pebble diameter(cm) 1 1 6.00 1 1 Fuel pebble number 1 1 10855 1 1 Graphite shell thickness(cm) 1 1 1.00 1 1 TRISO particles in one pebble 1 1 10000 1 1 Diameters of TRISOshells(outwards)(μm) 1 1 250/430/500/580 1 1 Active zone cylinder diameter(m) 1 1 1.35 1 1 Active zone cylinder height(m) 1 1 1.30 1 1 Active zone volume(m 3 ) 1 1 2.76 1 1 Upper reflector volume(m 3 ) 1 1 0.72 1 1 Lower reflector volume (m 3 ) 1 1 0.63 1 1 Upper plenum volume(m 3 ) 1 1 0.64 1 1 Lower plenum volume(m 3 ) 1 1 0.27 1 1 Diameter of pipe(m) 1 1 0.126 1 1 Length of hot leg(m) 1 1 10.0 1 1 Length of cold leg (m) 1 1 10.0 1 1 Reflector diameter(m) 1 1 2.85 1 1 Reflector height(m) 1 1 3.00 1 1 Total length of secondary loop(m) 1 1 20.0 1 1 FLiBe volume in IHX(m 3 ) 1 1 0.64 1 1 FLiNaK volume in IHX(m 3 ) 1 1 1.35 1 1 FLiNaK volume in SHX(m 3 ) 1 1 0.70 It is impossible to model the complicated bed structurewithout compromise because the core contains more than ten thousand pebbles and each pebble has thousands of TRISO particles. In OCFHR,the bed is divided into two areas: cylindrical graphite shells and cylindrical fuel mixing zone. The former structure strengthens the pebble while the latter contains TRISO particles and graphite moderator. Some parameters areacquired from neutron physics programsuch as temperature reactivity coefficients, decay constants and portions of neutron precursors.Heat transfer calculation of the core uses the mature porous medium Wakao correlation due to the impractical calculation of heat transfer between pebble bed and salt coolant. For heat transfer calculation between thebed and coolant, the Wakao correlation is often used [ 12 12 ]. <math display="block" id="M1"><mrow><mi>N</mi><mi>u</mi><mo>=</mo><mn>2.0</mn><mo>+</mo><mn>1.1</mn><mi>R</mi><msubsup><mi>e</mi><mi>D</mi><mrow><mn>0.6</mn></mrow></msubsup><mi>P</mi><msup><mi>r</mi><mrow><mn>1</mn><mo>/</mo><mn>3</mn></mrow></msup><mo>,</mo></mrow></math> where Nu is the Nusselt number, Re is the Reynolds number, Pr is the Prandtl number, and D is the equivalent diameter. Eq.(1) is applicable for 15 < Re < 8500[ 13 13 ]. Before calculation, some assumptions are adopted as follows: 1) The whole system is under perfect thermal insulation and all heat losses to the environment are neglected. 2) Consider the heat transfer between pebble and salt is merely heat convection. 3) Neglect the deformation of bed geometrycaused byvibration and flow of fluid. 4) Consider that average temperature of the fuel mixing zone equals to the surface temperature of fuelmixing zone in a pebble which is at the center of node while the node occupies the highest power distribution factor. 5) The heat fluxes of TRISO particles withinsame pebble are same. 2.3. Mathematic Model 2 The mass conservation in one fluid node conforms to that mass flow rate at the outlet equals to that at the inlet. The second item in Eq.(1) can be rewritten as one direction due to that the model is one dimensional. <math display="block" id="M2"><mrow><mfrac><mrow><mo>∂</mo><mi>ρ</mi></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mo>∇</mo><mo>·</mo><mrow><mo>(</mo><mrow><mi>ρ</mi><mi>V</mi></mrow><mo>)</mo></mrow><mo>=</mo><mn>0</mn><mo>,</mo></mrow></math> <math display="block" id="M3"><mrow><mfrac><mrow><mo>∂</mo><mi>ρ</mi></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>∂</mo><mi>ρ</mi><mi>u</mi></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac><mo>=</mo><mn>0</mn><mo>,</mo></mrow></math> where V is velocity vector; ρ is the density of coolant[kg/m 3 ]; and u is the average flow rate of coolant on x-dimension[m/s]. The coolant is driven by pump, which offers the pressure headto compensate the friction pressure drop, form pressure drop and gravitation pressure drop. The flow rate of coolant is determined by the difference between the pump pressure head and resistance pressure drop, <math display="block" id="M4"><mrow><mfrac><mrow><mo>∂</mo><mi>W</mi></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mo>∂</mo><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac><mrow><mo>(</mo><mrow><mfrac><mrow><msup><mi>W</mi><mn>2</mn></msup></mrow><mrow><mi>ρ</mi><mi>A</mi></mrow></mfrac></mrow><mo>)</mo></mrow><mo>+</mo><mfrac><mrow><mi>f</mi><msup><mi>W</mi><mn>2</mn></msup></mrow><mrow><mn>2</mn><mi>ρ</mi><mi>D</mi><mi>A</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><mi>K</mi><msup><mi>W</mi><mn>2</mn></msup></mrow><mrow><mn>2</mn><mi>ρ</mi><mi>A</mi></mrow></mfrac><mo>=</mo><mfrac><mrow><mo>∂</mo><mi>P</mi></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac><mo>,</mo></mrow></math> where W is mass flow rate of coolant [kg/s], f is the factor of friction pressure drop, K is the factor of form pressure drop, and P is pressure head of pump.Therefore, the integral form of Eq.(4) can be written as, <math display="block" id="M5"><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>m</mi></msubsup><mrow><msub><mrow><mrow><mo>(</mo><mrow><mfrac><mi>L</mi><mi>A</mi></mfrac></mrow><mo>)</mo></mrow></mrow><mi>i</mi></msub></mrow></mstyle><mfrac><mrow><mtext>d</mtext><mi>W</mi></mrow><mrow><mtext>d</mtext><mi>t</mi></mrow></mfrac><mo>+</mo><mi>Δ</mi><msub><mi>P</mi><mtext>f</mtext></msub><mo>+</mo><mi>Δ</mi><msub><mi>P</mi><mtext>K</mtext></msub><mo>+</mo><mi>Δ</mi><mi>P</mi><mo>=</mo><mn>0</mn><mo>,</mo></mrow></math> where m is the number of pipes or other sections, ΔP f is the friction pressure drop,and ΔP K is the form pressure drop. Fora TRISOparticle,thermal conductivityformula is adopted to calculate temperature gradient from its outer pyrolytic carbon shell to fuel particle.Based on the simplifications and assumptionsabove, the temperature of center fuel zone in the hottest TRISO can be calculated as, <math display="block" id="M6"><mrow><msub><mi>c</mi><mrow><mtext>pf</mtext></mrow></msub><msub><mi>m</mi><mtext>f</mtext></msub><mfrac><mrow><mo>∂</mo><msub><mi>T</mi><mtext>f</mtext></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>=</mo><mo>−</mo><mfrac><mrow><mrow><mo>(</mo><mrow><msub><mi>T</mi><mtext>f</mtext></msub><mo>−</mo><msub><mi>T</mi><mtext>g</mtext></msub></mrow><mo>)</mo></mrow></mrow><mrow><msub><mi>R</mi><mrow><mtext>fg</mtext></mrow></msub></mrow></mfrac><mo>+</mo><msub><mi>P</mi><mtext>f</mtext></msub><mi>F</mi><mo>.</mo></mrow></math> The heat transfer form in pebble graphite shell can be calculated by: <math display="block" id="M7"><mrow><msub><mi>c</mi><mrow><mtext>pg</mtext></mrow></msub><msub><mi>m</mi><mtext>g</mtext></msub><mfrac><mrow><mo>∂</mo><msub><mi>T</mi><mtext>g</mtext></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>=</mo><mfrac><mrow><mrow><mo>(</mo><mrow><msub><mi>T</mi><mtext>f</mtext></msub><mo>−</mo><msub><mi>T</mi><mtext>g</mtext></msub></mrow><mo>)</mo></mrow></mrow><mrow><msub><mi>R</mi><mrow><mtext>fg</mtext></mrow></msub></mrow></mfrac><mo>−</mo><mrow><mo>(</mo><mrow><msub><mi>T</mi><mtext>g</mtext></msub><mo>−</mo><msub><mi>T</mi><mtext>s</mtext></msub></mrow><mo>)</mo></mrow><msub><mi>h</mi><mrow><mtext>gs</mtext></mrow></msub><msub><mi>A</mi><mrow><mtext>gs</mtext></mrow></msub><mo>−</mo><msub><mi>ε</mi><mtext>g</mtext></msub><msub><mi>A</mi><mrow><mtext>ge</mtext></mrow></msub><mi>σ</mi><mrow><mo>(</mo><mrow><msubsup><mi>T</mi><mtext>g</mtext><mn>4</mn></msubsup><mo>−</mo><msubsup><mi>T</mi><mrow><mtext>ge</mtext></mrow><mn>4</mn></msubsup></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where P f is the reactor power[W]; F is the power axial distribution factor; c pf m f is the specific heat capacity of the fuel mixing zone [J/°C]; T f is the fuel mixing zone temperature [°C]; T s is the primary coolant temperature [°C]; T ge is the graphite reflector temperature [°C]; T g is the external graphite temperature in a node [°C]; c pg m g is the specific heat capacity of the node [J/°C]; h gs is equivalent heat transfer coefficient between graphite and molten salt [W·m −2 ·°C −1 ]; A gs is the heat transfer area between graphite zone and external salt; A ge is the radiation area between graphite zone and graphite reflector; ε g is the effective radiation coefficient between external graphite zone and reflector; and R fg is the thermal resistance between internal fuel mixing zone and external graphite. The subscripts f, g, s and ge denote the mixing zone, graphite shell, primary coolant and reflector, respectively. Fig 5 Fig 5 shows the heat transfer process in the core. In aprimary coolant node of core, <math display="block" id="M8"><mrow><msub><mi>c</mi><mrow><mtext>ps</mtext></mrow></msub><msub><mi>m</mi><mtext>s</mtext></msub><mfrac><mrow><mo>∂</mo><msub><mi>T</mi><mtext>s</mtext></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>=</mo><mrow><mo>(</mo><mrow><msub><mi>T</mi><mtext>g</mtext></msub><mo>−</mo><msub><mi>T</mi><mtext>s</mtext></msub></mrow><mo>)</mo></mrow><msub><mi>h</mi><mrow><mtext>gs</mtext></mrow></msub><msub><mi>A</mi><mrow><mtext>gs</mtext></mrow></msub><mo>−</mo><mrow><mo>(</mo><mrow><msub><mi>T</mi><mtext>s</mtext></msub><mo>−</mo><msub><mi>T</mi><mrow><mtext>ge</mtext></mrow></msub></mrow><mo>)</mo></mrow><msub><mi>h</mi><mrow><mtext>ge</mtext></mrow></msub><msub><mi>A</mi><mrow><mtext>ge</mtext></mrow></msub><mo>−</mo><msub><mi>c</mi><mrow><mtext>ps</mtext></mrow></msub><msub><mover accent="true"><mi>m</mi><mo>˙</mo></mover><mtext>s</mtext></msub><mrow><mo>(</mo><mrow><msub><mi>T</mi><mrow><mtext>out</mtext></mrow></msub><mo>−</mo><msub><mi>T</mi><mrow><mtext>in</mtext></mrow></msub></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where m s is the coolant mass in one node [kg]; c ps is the coolant specific heat [J·kg −1 ·°C −1 ]; <math id="M9"><mrow><msub><mover accent="true"><mi>m</mi><mo>˙</mo></mover><mi>s</mi></msub></mrow></math> is the coolant mass flow rate [kg/s]; and T out and T in are molten salt outlet and inlet temperatures, respectively,in one node. Inagraphite reflectornode, <math display="block" id="M10"><mrow><msub><mi>c</mi><mrow><mtext>pge</mtext></mrow></msub><msub><mi>m</mi><mrow><mtext>ge</mtext></mrow></msub><mfrac><mrow><mo>∂</mo><msub><mi>T</mi><mrow><mtext>ge</mtext></mrow></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>=</mo><msub><mi>ε</mi><mrow><mtext>ge</mtext></mrow></msub><msub><mi>A</mi><mrow><mtext>ge</mtext></mrow></msub><mi>σ</mi><mrow><mo>(</mo><mrow><msubsup><mi>T</mi><mtext>g</mtext><mn>4</mn></msubsup><mo>−</mo><msubsup><mi>T</mi><mrow><mtext>ge</mtext></mrow><mn>4</mn></msubsup></mrow><mo>)</mo></mrow><mo>+</mo><mrow><mo>(</mo><mrow><msub><mi>T</mi><mtext>s</mtext></msub><mo>−</mo><msub><mi>T</mi><mrow><mtext>ge</mtext></mrow></msub></mrow><mo>)</mo></mrow><msub><mi>h</mi><mrow><mtext>ge</mtext></mrow></msub><msub><mi>A</mi><mrow><mtext>ge</mtext></mrow></msub><mo>−</mo><mrow><mo>(</mo><mrow><msub><mi>T</mi><mrow><mtext>ge</mtext></mrow></msub><mo>−</mo><msub><mi>T</mi><mrow><mtext>sd</mtext></mrow></msub></mrow><mo>)</mo></mrow><msub><mi>h</mi><mrow><mtext>gd</mtext></mrow></msub><msub><mi>A</mi><mrow><mtext>gd</mtext></mrow></msub><mo>,</mo></mrow></math> where c pge is the graphite specific heat [J·kg −1 ·°C −1 ]; m ge is the graphite reflector mass of one node; T ge is the graphite temperature in one node [°C]; T sd is the temperature of salt in one node of downcomer [°C]; h ge is the equivalent heat transfer coefficient between salt and reflector [W·m −2 ·°C −1 ]; h gd is the equivalent heat transfer coefficient between salt in the downcomer and reflector [W·m −2 ·°C −1 ]; and A gd is the effective heat transfer area. In a coolant node of downcomer, <math display="block" id="M11"><mrow><msub><mi>c</mi><mrow><mtext>ps</mtext></mrow></msub><msub><mi>m</mi><mrow><mtext>sd</mtext></mrow></msub><mfrac><mrow><mo>∂</mo><msub><mi>T</mi><mrow><mtext>sd</mtext></mrow></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>=</mo><mrow><mo>(</mo><mrow><msub><mi>T</mi><mrow><mtext>ge</mtext></mrow></msub><mo>−</mo><msub><mi>T</mi><mrow><mtext>sd</mtext></mrow></msub></mrow><mo>)</mo></mrow><msub><mi>h</mi><mrow><mtext>gd</mtext></mrow></msub><msub><mi>A</mi><mrow><mtext>gd</mtext></mrow></msub><mo>−</mo><msub><mi>c</mi><mrow><mtext>ps</mtext></mrow></msub><msub><mover accent="true"><mi>m</mi><mo>˙</mo></mover><mtext>s</mtext></msub><mrow><mo>(</mo><mrow><msub><mi>T</mi><mrow><mtext>out</mtext></mrow></msub><mo>−</mo><msub><mi>T</mi><mrow><mtext>in</mtext></mrow></msub></mrow><mo>)</mo></mrow><mo>−</mo><msub><mover accent="true"><mi>P</mi><mo>˙</mo></mover><mn>0</mn></msub><mo>,</mo></mrow></math> where m sd is the coolant mass in one node [kg]; <math id="M12"><mrow><msub><mover accent="true"><mi>m</mi><mo>˙</mo></mover><mtext>s</mtext></msub></mrow></math> is the coolant mass flow rate [kg/s]; and <math id="M13"><mrow><msub><mover accent="true"><mi>P</mi><mo>˙</mo></mover><mn>0</mn></msub></mrow></math> is the heat release rate of the reactor vessel, which is set as zero under the first assumption. For the pressure drop in the core, the Ergun equation, which is suitable for porous medium, can be used[ 14 14 ]: <math display="block" id="M14"><mrow><mfrac><mrow><mi>Δ</mi><msub><mi>P</mi><mtext>c</mtext></msub></mrow><mi>L</mi></mfrac><mo>=</mo><mo>−</mo><mfrac><mrow><mn>150</mn><mi>μ</mi></mrow><mrow><msubsup><mi>D</mi><mtext>p</mtext><mn>2</mn></msubsup></mrow></mfrac><mfrac><mrow><msup><mrow><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><mi>ε</mi></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow><mrow><msup><mi>ε</mi><mn>3</mn></msup></mrow></mfrac><msub><mi>u</mi><mi>i</mi></msub><mo>−</mo><mfrac><mrow><mn>1.75</mn><mi>ρ</mi></mrow><mrow><msub><mi>D</mi><mi>p</mi></msub></mrow></mfrac><mfrac><mrow><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><mi>ε</mi></mrow><mo>)</mo></mrow></mrow><mrow><msup><mi>ε</mi><mn>3</mn></msup></mrow></mfrac><mrow><mo>|</mo><mi>u</mi><mo>|</mo></mrow><msub><mi>u</mi><mi>i</mi></msub><mo>,</mo></mrow></math> where, Δ P C is the core pressure drop, ε is the volume fraction, μ is the viscosity, D p is the pebble diameter, and L is the height of pebble bed. Nuclear power calculation is based on the fifteen-group point kinetic equations: <math display="block" id="M15"><mrow><mfrac><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>d</mi></mstyle><mi>n</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>d</mi></mstyle><mi>t</mi></mrow></mfrac><mo>=</mo><mfrac><mrow><mi>ρ</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>−</mo><mi>β</mi></mrow><mi>l</mi></mfrac><mi>n</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>+</mo><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mrow><mn>15</mn></mrow></munderover><msub><mi>λ</mi><mi>i</mi></msub><msub><mi>C</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mfrac><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>d</mi></mstyle><msub><mi>C</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>d</mi></mstyle><mi>t</mi></mrow></mfrac><mo>=</mo><mfrac><mrow><msub><mi>β</mi><mi>i</mi></msub></mrow><mi>l</mi></mfrac><mi>n</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>−</mo><msub><mi>λ</mi><mi>i</mi></msub><msub><mi>C</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>,</mo></mrow></math> where n is the neutron density [n/m 3 ]; C is the concentration of neutron precursor; β is the total delayed neutron fraction; β i is the i th delayed neutron fraction; λ is the decay constant [s −1 ]; l is the mean neutron generation time [s]. β , λ and l are gained from neutron physics calculation while t he external source of neutron is neglected.Among the 15delayed-neutron precursor groups, nine groups are added in due to the photo-neutron emission. This enlarges neutron lifetime and increases the effects of delayed neutrons[ 15 15 ]. Negative reactivity feedback modelcontains the effects of fuel, coolant, moderator and reflector. Fuel reactivity coefficient is negative due to the Doppler effect, and that of primary coolant also shows negativedue to its neutron property. Besides, the control rod effect is an important way to control reactivity in most FHR designs. <math display="block" id="M16"><mrow><mi>ρ</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>=</mo><msub><mi>ρ</mi><mn>0</mn></msub><mo>+</mo><msub><mi>ρ</mi><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi><mi>d</mi></mstyle></mrow></msub><mo>+</mo><msub><mi>α</mi><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>f</mi><mi>g</mi></mstyle></mrow></msub><mstyle mathvariant="bold" mathsize="normal"><mi>Δ</mi></mstyle><msub><mover accent="true"><mi>T</mi><mo>¯</mo></mover><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>f</mi><mi>g</mi></mstyle></mrow></msub><mo>+</mo><msub><mi>α</mi><mstyle mathvariant="bold" mathsize="normal"><mi>s</mi></mstyle></msub><mstyle mathvariant="bold" mathsize="normal"><mi>Δ</mi></mstyle><msub><mover accent="true"><mi>T</mi><mo>¯</mo></mover><mstyle mathvariant="bold" mathsize="normal"><mi>s</mi></mstyle></msub><mo>+</mo><msub><mi>α</mi><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>g</mi><mi>e</mi></mstyle></mrow></msub><mstyle mathvariant="bold" mathsize="normal"><mi>Δ</mi></mstyle><msub><mover accent="true"><mi>T</mi><mo>¯</mo></mover><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>g</mi><mi>e</mi></mstyle></mrow></msub><mo>,</mo></mrow></math> where ρ 0 is the initial reactivity; ρ rd is the control rod reactivity; α fg is the mixing zone reactivity coefficient; α S is the primary coolant reactivity coefficient; α ge is the graphite reflector reactivity coefficient; <math id="M17"><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>Δ</mi></mstyle><msub><mover accent="true"><mi>T</mi><mo>¯</mo></mover><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>f</mi><mi>g</mi></mstyle></mrow></msub></mrow></math> , <math id="M18"><mrow><mo> </mo><mstyle mathvariant="bold" mathsize="normal"><mi>Δ</mi></mstyle><msub><mover accent="true"><mi>T</mi><mo>¯</mo></mover><mstyle mathvariant="bold" mathsize="normal"><mi>s</mi></mstyle></msub></mrow></math> and <math id="M19"><mrow><mo> </mo><mstyle mathvariant="bold" mathsize="normal"><mi>Δ</mi></mstyle><msub><mover accent="true"><mi>T</mi><mo>¯</mo></mover><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>g</mi><mi>e</mi></mstyle></mrow></msub></mrow></math> are mean temperature change values [°C] of the mixing zone, primary coolant and reflector, respectively. The IHX modelis the smooth shell-and tube type heat exchanger (STHX).Compared with spirally-fluted tube, its pressure drop and fabricating cost is lower, while heat transfer efficiency is weakened by half. Fig.6 Fig.6 shows a cross-section of the IHX. The IHX operates under a FLiBe-to-FLiNaK environment while the primary coolant flows on the tube side and the secondary coolant flows on the shell side. Detailed parameters of IHX are shown in Table 2 Table 2 . 1 1 Parameters 1 1 Values 1 1 IHX primary side inlet temperature(°C) 1 1 700.0 1 1 IHX primary side outlet temperature(°C) 1 1 672.0 1 1 IHX secondary side inlet temperature(°C) 1 1 610.0 1 1 IHX secondary side outlet temperature(°C) 1 1 630.0 1 1 IHX tube number 1 1 1210 1 1 IHX tube thickness(mm) 1 1 1.65.0 1 1 IHX effective tube length(m) 1 1 2.8.0 1 1 IHX equivalent heat transfer coefficient(  W·m −2 ·K −1   ) 1 1 1200.0 1 1 SHX secondary side inlet temperature(°C) 1 1 40.0 1 1 SHX secondary side outlet temperature(°C) 1 1 151.0 1 1 SHX tube number 1 1 2016 1 1 SHX tube thickness(mm) 1 1 1.7 1 1 SHX effective tube length(m) 1 1 3.3 1 1 SHX equivalent heat transfer coefficient(W·m −2 ·K −1 ) 1 1 108.0 For the tube side, <math display="block" id="M20"><mrow><msub><mi>c</mi><mrow><mtext>ps</mtext></mrow></msub><msub><mi>m</mi><mrow><mtext>ep</mtext></mrow></msub><mfrac><mrow><mo>∂</mo><msub><mi>T</mi><mrow><mtext>ep</mtext></mrow></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>=</mo><mo>−</mo><mrow><mo>(</mo><mrow><msub><mi>T</mi><mrow><mtext>ep</mtext></mrow></msub><mo>−</mo><msub><mi>T</mi><mrow><mtext>ed</mtext></mrow></msub></mrow><mo>)</mo></mrow><msub><mi>h</mi><mrow><mtext>ep</mtext></mrow></msub><msub><mi>A</mi><mrow><mtext>ep</mtext></mrow></msub><mo>+</mo><msub><mi>c</mi><mrow><mtext>ps</mtext></mrow></msub><msub><mover accent="true"><mi>m</mi><mo>˙</mo></mover><mtext>s</mtext></msub><mrow><mo>(</mo><mrow><msub><mi>T</mi><mrow><mtext>inep</mtext></mrow></msub><mo>−</mo><msub><mi>T</mi><mrow><mtext>outep</mtext></mrow></msub></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where c ps is the internal salt specific heat [J·kg −1 ·°C −1 ]; m ep is internal coolant mass in a node [kg]; T ep is internal coolant temperature [°C]; T ed is the intermediate wall temperature [°C]; h ep is the equivalent heat transfer coefficient between salt in the tube and wall [W·m −2 ·°C −1 ]; A p is the effective heat transfer area; and T inep and T outep are the tube side coolant temperatures at the tube inlet and outlet, respectively. For the shell side, <math display="block" id="M21"><mrow><msub><mtext>c</mtext><mrow><mtext>ss</mtext></mrow></msub><msub><mi>m</mi><mrow><mtext>es</mtext></mrow></msub><mfrac><mrow><mo>∂</mo><msub><mi>T</mi><mrow><mtext>es</mtext></mrow></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>=</mo><mrow><mo>(</mo><mrow><msub><mi>T</mi><mrow><mtext>ed</mtext></mrow></msub><mo>−</mo><msub><mi>T</mi><mrow><mtext>es</mtext></mrow></msub></mrow><mo>)</mo></mrow><msub><mi>h</mi><mrow><mtext>es</mtext></mrow></msub><msub><mi>A</mi><mrow><mtext>es</mtext></mrow></msub><mo>+</mo><msub><mi>c</mi><mrow><mtext>ss</mtext></mrow></msub><msub><mover accent="true"><mi>m</mi><mo>˙</mo></mover><mrow><mtext>ss</mtext></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>T</mi><mrow><mtext>inss</mtext></mrow></msub><mo>−</mo><msub><mi>T</mi><mrow><mtext>outss</mtext></mrow></msub></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where c ss is the external salt specific heat [J·kg −1 ·°C −1 ]; m es is external coolant mass in a node [kg]; T es is coolant temperature [°C]; h es is the equivalent heat transfer coefficient between salt in the tube and wall [W·m −2 ·°C −1 ]; A es is the effective heat transfer area; and T inss and T outss are the shell side coolant temperatures at the shell inlet and outlet, respectively. The heat transfer between tube side and shell side of IHX and SHX is heat convection. Based on the simplification of heat exchanger, pipe bundles are approximated to a pipe of thesame length, equivalent diameter and heat exchanger area. The correlations used in this study are listed as follows[ 16 16 ]. For the tube side <math display="block" id="M22"><mrow><mi>N</mi><mi>u</mi><mo>=</mo><mn>0.027</mn><mi>R</mi><msup><mi>e</mi><mrow><mn>0.8</mn></mrow></msup><mi>P</mi><msup><mi>r</mi><mrow><mn>1</mn><mo>/</mo><mn>3</mn></mrow></msup><msup><mrow><mrow><mo>(</mo><mrow><mfrac><mi>μ</mi><mrow><msub><mi>μ</mi><mtext>w</mtext></msub></mrow></mfrac></mrow><mo>)</mo></mrow></mrow><mrow><mn>0.14</mn></mrow></msup><mo>,</mo></mrow></math> For the shell side <math display="block" id="M23"><mrow><mi>N</mi><mi>u</mi><mo>=</mo><mn>0.3</mn><mo>+</mo><mfrac><mrow><mn>0.62</mn><mi>R</mi><msup><mi>e</mi><mrow><mn>0.5</mn></mrow></msup><mi>P</mi><msup><mi>r</mi><mrow><mn>1</mn><mo>/</mo><mn>3</mn></mrow></msup></mrow><mrow><msup><mrow><mrow><mo>[</mo><mrow><mn>1</mn><mo>+</mo><msup><mrow><mrow><mo>(</mo><mrow><mn>0.4</mn><mo>/</mo><mi>P</mi><mi>r</mi></mrow><mo>)</mo></mrow></mrow><mrow><mn>2</mn><mo>/</mo><mn>3</mn></mrow></msup></mrow><mo>]</mo></mrow></mrow><mrow><mn>1</mn><mo>/</mo><mn>4</mn></mrow></msup></mrow></mfrac><msup><mrow><mrow><mo>[</mo><mrow><mn>1</mn><mo>+</mo><msup><mrow><mrow><mo>(</mo><mrow><mfrac><mrow><mi>R</mi><mi>e</mi></mrow><mrow><mn>28200</mn></mrow></mfrac></mrow><mo>)</mo></mrow></mrow><mrow><mn>5</mn><mo>/</mo><mn>8</mn></mrow></msup></mrow><mo>]</mo></mrow></mrow><mrow><mn>0.8</mn></mrow></msup><mo>,</mo></mrow></math> where μ and μ w are coolant viscosities calculated by coolant temperature and wall temperature, respectively. Currently, FHER adoptsan air fanas the load. In OCFHR, the AR model also uses an air fan as the heat sinkfor the hot secondary coolant flowingin the SHX from IHXcold side outlet. The type of SHX model is similar to IHX, with FLiNaK flowing in the tube and the air flowing in the shell. The inlet air is at 40°Cunder normal conditions. 3 Results and Discussion 1 3.1. Steady-state operation 2 Only when the model reaches its rated steady state, will disturbed transients be introduced.The steady-state operation parameters,node temperatures and delayed neutron group parameters are givenin Tables 3 Tables 3 , 4 4 and 5 5 , respectively. * TCR , temperature coefficient of reactivity 1 1 Parameters 1 1 Values 1 1 Total thermal power(MW) 1 1 10.0 1 1 Primary coolant inlet temperature(°C) 1 1 672.0 1 1 Primarycoolant outlet temperature(°C) 1 1 700.0 1 1 Primary coolant mass flow rate(kg·s −1 ) 1 1 150.0 1 1 Pressure drop of pebble bed(Pa) 1 1 2200 1 1 Fuel mixing zone TCR (pcm·°C −1 ) 1 1 −5.69 1 1 Salt TCR(pcm·°C −1 ) 1 1 −2.27 1 1 Graphite TCR(pcm·°C −1 ) 1 1 0.74 1 1 Secondary coolant mass flow rate(kg·s −1 ) 1 1 260.0 1 1 Secondary coolant inlet temperature(°C) 1 1 610.0 1 1 Secondary coolant outlet temperature(°C) 1 1 630.0 1 1 Air fan mass flow rate(kg·s −1 ) 1 1 51.0 1 1 Inlet of downcomer 1 1 672 1 1 Inlet of IHX hot side 1 1 700 1 1 Outlet of IHX hot side 1 1 672 1 1 Outlet of IHX cold side 1 1 630 1 1 Inlet of IHX cold side 1 1 610 1 1 Inlet of SHXhot side 1 1 630 1 1 Outlet of SHX hot side 1 1 610 1 1 Outlet of SHX cold side 1 1 227 1 1 Inlet of SHX cold side 1 1 40 1 1 No. 1 1 β 1 1 λ 1 1 1 1 1 2.11E-4 1 1 1.24E-2 1 1 2 1 1 1.395E-3 1 1 3.05E-2 1 1 3 1 1 1.25E-3 1 1 1.11E-1 1 1 4 1 1 2.514E-3 1 1 3.01E-1 1 1 5 1 1 7.35E-4 1 1 1.14 1 1 6 1 1 2.684E-4 1 1 3.02 1 1 7 1 1 5.512E-7 1 1 6.24E-7 1 1 8 1 1 3.67E-7 1 1 2.48E-6 1 1 9 1 1 2.514E-6 1 1 1.59E-5 1 1 10 1 1 3.094E-5 1 1 6.2E-5 1 1 11 1 1 3.481E-6 1 1 2.67E-4 1 1 12 1 1 3.559E-5 1 1 7.42E-4 1 1 13 1 1 1.789E-5 1 1 3.6E-3 1 1 14 1 1 3.54E-5 1 1 8.85E-3 1 1 15 1 1 2.0E-5 1 1 2.26E-2 3.2. Model solver verification 2 OCFHR employs a 15-group point kinetic model to calculate nuclear power. To verify the model’s reliability, the isolated point-reactor model is compared with thatbased on Gear and NumericalDifferentiationFormulas (NDFs) algorithmsby inserting a negative step reactivity of −100pcm.Also, a single-group point kinetic model based on the R-K method is compared with the analytical solutionby insertinga positive step reactivity of 250pcmto test its accuracy. The time step is 0.01s and the results are shown in Fig. 7 Fig. 7 . In Fig. 7(a) Fig. 7(a) , the numerical solution of OCFHR agree well with those of Gear and NDFs. Although different algorithms have different error handling methods which may cause small error, these results provethat R-K method is viable. From Fig. 7(b) Fig. 7(b) , R-K numerical solution matches well with single-group analytical solutionand the truncation error can be neglected[ 17 17 ]. The thermal hydraulic calculation is based on lumped parameter method. A simple case ( Fig.8 Fig.8 and Table 6 Table 6 ) is simulated and compared with RELAP5:FLiBe salt flows through a heated pipe ora cooled pipe.The start temperature of the salt and walls is 600°C. Then it begins to change as the heating power and cooling power isadded to pipe walls. 1 1 Parameters 1 1 Values 1 1 Pipe length (m) 1 1 50 1 1 Pipe thickness (m) 1 1 0.05 1 1 Pipe cross section area (m 2 ) 1 1 1.0 1 1 Pipe diameter (m) 1 1 1.128 1 1 Salt inlet temperature (°C) 1 1 600 1 1 Mass flow rate (kg/s) 1 1 100 1 1 Heater power density  (W/m 2 )  1 1 50000 1 1 Cooler power density (W/m 2 ) 1 1 −50000 1 1 Radial node number of heater and cooler 1 1 10 1 1 Axial node number of heater and cooler 1 1 20 1 1 Calculation time range (s) 1 1 10000 Figure 9 Figure 9 showsthat the temperatures of walls near the heater and the coolerresponse quickly while the temperature of cold salt drops at first. Nodes 1 and 2 arethe first and last wall node right next to the fluid,respectively.The maximum wall temperature difference between OCFHR and RELAP5 is 1.1°Cwith a calculation deviation of0.11%, caused probably by different calculation methods.The errors of salt and wall temperatures in Fig.10 Fig.10 prove the validity of thermal hydraulic calculation methodand its reliability tosimulate transient behavior of the FHER. 3.3. Transient simulation 2 OCFHR is programmed for simulating and analyzingnormal operation transient behavior of FHR. Before simulation, the types of operation transient should be clarified. The disturbances simulated are based on the possibilities of happening during normal operation.According to experimental and engineering experience, fourtypical transients are listed in Table7. Table7. 1 1 Title 1 1 Reason 1 1 Value 1 1 Reactivity step insertion 1 1 Maloperation of control rod 1 1 ±100(pcm) 1 1 Coolant mass flow rate step change 1 1 Maloperation of pump  /short electricity cut of pump 1 1 ±  (0–10%) 1 1 Air inlet   temperature change 1 1 Ambient temperature change 1 1 ±10(°C) 1 1 Air inlet flow rate change 1 1 Misoperation of air radiator 1 1 ±(0–20%) The dynamic model must reach a rated steady state in the first place. Then, the disturbed transients shown above will be introduced into the system.Considering the reactivity inserted and mass flow rate changes in the transients are oflarger impacts than other disturbances, and all safety control modules are forbidden during the simulations, this work can help to figure outthe operation stability and security of FHER. 3.3.1 100pcmpositive reactivity step insertion 3 A reactivity step insertion accident mayhappen during normal operation, ifcontrol rod is extracted or inserted by accident, or even by a blackout. When reactivity is step inserted into the reactor, the nuclear power may fluctuate strongly in a short time. Figure 11(a) Figure 11(a) shows tha tafter a step insertion of positive reactivity, the reactor power reaches the maximum power at about 20sbefore it is overwhelmed bythe negative temperature reactivity feedback mainly produced from the fuel and coolant. At about 800s, the power reaches a steady state of 10.36MW, beingabout 0.36MW higher than the origin. Fig.11(b) Fig.11(b) shows the coolant temperature change after the reactivity step insertion.Following the rising nuclear power, the temperature of the pebble center increases quicklytothe peak temperature of 955°C, which is much less than the design critical temperature. After nearly 65s, the coolant temperature at core inlet begins to increase while the interval is roughlythe same as the circulation time of primary loop. The secondary coolant temperature atIHX outlet begins to increase at about 22s and the coolant at IHX inlet increases at 30s, whichis also similar to the circulation time of secondary loop.After nearly 600s, the temperatures of systemstabilize at new operating points.The average temperatures of fuel mixing zone and graphite increase by 14°C and 12°C. The temperatures of primary loop at IHX inlet and outlet increase by about10°C and 11°C, while those in the secondary loopare 12.2°C and 12.6°C, respectively. 3.3.2 10% of primary coolant mass flow rate step increases 3 This accidenthappens due to the shifting of pump operation state or a breakdown, and 10% increase of mass flow rate is relatively large for normal operation. This disturbance can decrease fuel temperature for a time, and then fuel temperature increases as positive feedback is inserted. Fig.12(a) Fig.12(a) shows that when the primary coolant mass flow rate risesdueto incorrect operation of the pump, the heat exchange in core becomes intensive. Positive reactivity is added into the core due to the negative reactivity feedback of fuel and coolant. After that, negative reactivity is induced back to counteract the previous effect as temperature continues to go up. At about 800s, a stable value is reached and the power is nearly 0.1MW higher than its initial value. Fig. 12(b) Fig. 12(b) shows fluid temperature curve in this process. It is observed that the primary coolant temperature decreases first due to the extra flow rate, then increases going after the rising power. However,the temperature of fluid and fuel merely changes during the case. Still, after nearly 300s, all temperatures stabilizeapparently. The average temperatures of fuel zone and graphite shell decrease by only 0.1°C and 0.6°C. The core inlet temperature ascends by 2°C and the core outlet temperature decreases by 0.5°C. Besides, the secondary coolant temperature increases 3.3°C at IHX inlet and 3.5°C at IHX outlet, approximately. It can be concluded that the change of fluid flow rate matters little on the fluid temperature. 3.3.3 10% of secondary coolant mass flow rate step increases 3 As the secondary coolant mass flow rate changes, the nuclear power needs a relatively longer time to response. This is due to the large heat transfer delay of the system and the large thermal inertia of the molten salt. Fig. 13(a) Fig. 13(a) shows that at about 40s, the power begins to response according to the positive reactivity inserted by fuel and coolant temperature changes. At nearly 600s, the power becomes stable at 10.046MW. From Fig. 13(b) Fig. 13(b) , at the start of this simulation, the temperature of secondary coolant changes quickly but not apparently, and the primary coolant temperature varies inconspicuously. FLiBe temperatures at core inlet and outlet decrease by 0.34°C and 0.4°C, respectively. The FLiNaK temperatures at IHX inlet and outlet increase by 2.4°C and 0.6°C, respectively.The characteristics of molten salt make the change of mass flow have little influence on the nuclear power and primary coolant operation condition. When this incorrect operation occurs, FHR can deal well with the case without any other intervening measures. 3.3.4 20% of air mass flow rate step increases 3 Reactor’s load cannot be shifted synchronously with the nuclear power, and sometimes step change may occur in its operation point. This case is used in the load following operation mode. Fig.14(a) Fig.14(a) shows that the power increases slowly from 50s and reaches steady-state at about 1000s. The maximum of power is about 11.05MW. Fig. 14(b) Fig. 14(b) shows that temperatureof theaverage fuel zone increases by 2°C while the shell average temperature decreases by 4.5°C. The temperatures of primary coolant in core inlet and outlet begin to decrease at nearly 40s and reach to a stable value at about 600s, which are about 8°Clower and 6°C lower than the initial values respectively. On the contrary, the secondary coolant temperature responds quickly and it reaches stable state at about 600s. The secondary coolant temperatures in IHX inlet and outlet decrease by18.3°C and 16°C, respectively. From Fig.14 Fig.14 , wecan know that the nuclear power is connected to the fan’s flow rate, but not a positive correlation relationship, since the fan’s air temperature changes with the nuclear power. 3.3.5 10°C of air generator inlet temperatureincreases 3 Due to the use of fan as a PCU, a new disturbance of air inlet temperature is inserted into FHER system if no equipment exists to control the inlet air temperature, and temperature change of 10°C is large enough for surroundings. We assume that the air inlet temperature increases from 40°C to 50°C in 10 min to research the transient variation of the FHER. In Fig. 15(a) Fig. 15(a) , the nuclear power decreases at about 50s and stabilizes at its minimum value of 9.82 MW after 1500s.This is because the increasing air inlet temperature can impede heat transfer. The temperature of secondary coolant increases slowly, and the temperatures of fuel, shell and primary coolant do not change apparently. Both of the coolant temperatures become stable after about 800s. Besides, the temperature of fuel zone decreases by 0.5°C and the temperature of graphite shell increases by 1°C. The temperatures of primary coolant in core inlet and outlet increase by 1.8°C and 1.1°C, respectively. The secondary coolant temperatures in IHX inlet and outlet increase by 3.6°C and 3.2°C, respectively. These results argue the necessity to control the variation range and rate of inlet temperature of air coolant. 3.3.6 Variations of FHER during transient disturbances 3 We summarizethe simulation results above in Tables 8 Tables 8 and 9 9 , so to analyze influence of the disturbance and to compare with other disturbances. In the first rowof the tables, No.1–5 represent the five transient states. The ‘+’ and ‘−’ signs denote the change direction of variables. The response time isthe time period fromthe initial time of transient disturbances insertion to apparent change of target variable. 1 1 Parameters 1 1 0 1 1 1 1 1 2 1 1 3 1 1 4 1 1 5 1 1 Power(MW) 1 1 10 1 1 10.36 1 1 10.07 1 1 10.05 1 1 10.96 1 1 9.81 1 1 Max deviation from initial Power(MW) 1 1 ― 1 1 +3.76 1 1 +0.28 1 1 +0.06 1 1 +0.96 1 1 −0.19 1 1 Fuel zone average temperature (°C) 1 1 777.1 1 1 791.2 1 1 777.0 1 1 777.2 1 1 779.0 1 1 776.7 1 1 Max deviation from initialfuel zone average temperature (°C) 1 1 ― 1 1 +14.1 1 1 −0.1 1 1 +0.12 1 1 +1.9 1 1 −0.4 1 1 Graphite average temperature(°C) 1 1 710.4 1 1 722.1 1 1 709.8 1 1 710.2 1 1 706.0 1 1 711.3 1 1 Max deviation from initialgraphite average temperature(°C) 1 1 ― 1 1 +11.7 1 1 −0.6 1 1 −1.9 1 1 −4.4 1 1 +0.9 1 1 Primary coolant temperature at IHX inlet(°C) 1 1 700.0 1 1 711.5 1 1 699.5 1 1 699.7 1 1 694.4 1 1 701.1 1 1 Max deviation from initial primary coolant temperature at IHX inlet(°C) 1 1 ― 1 1 +11.5 1 1 −0.5 1 1 −0.34 1 1 −5.6 1 1 +1.1 1 1 Primary coolant temperature at IHX outlet(°C) 1 1 672.4 1 1 683.1 1 1 674.3 1 1 672.04 1 1 664.1 1 1 674.1 1 1 Max deviation from initial primary coolant temperature at IHX outlet(°C) 1 1 ― 1 1 +10.7 1 1 +1.9 1 1 −0.4 1 1 −8.3 1 1 +1.8 1 1 Secondary coolant temperature at IHX inlet(°C) 1 1 610.1 1 1 622.3 1 1 613.4 1 1 612.5 1 1 591.8 1 1 613.7 1 1 Max deviation from initialsecondary coolant temperature at IHX inlet(°C) 1 1 ― 1 1 +12.2 1 1 +3.3 1 1 +2.4 1 1 −18.3 1 1 +3.6 1 1 Secondary coolant temperature at IHX outlet(°C) 1 1 630.2 1 1 642.8 1 1 633.7 1 1 630.8 1 1 614.1 1 1 633.4 1 1 Max deviation from initial secondary coolant temperature at IHX outlet(°C) 1 1 ― 1 1 +12.6 1 1 +3.5 1 1 +0.6 1 1 −16.1 1 1 +3.2 * P.C. , primary coolant; S.C. , secondary coolant 1 1 Transient states 1 1 1 1 1 2 1 1 3 1 1 4 1 1 5 1 1 Power(MW) 1 1 ― 1 1 5.5 1 1 22.5 1 1 40.0 1 1 60.0 1 1 P.C. at core inlet 1 1 68.0 1 1 23.0 1 1 23.0 1 1 30.0 1 1 168.0 1 1 P.C. at coreoutlet 1 1 5.5 1 1 1.0 1 1 40.0 1 1 46.0 1 1 183.0 1 1 S.C. at IHX inlet 1 1 30.5 1 1 25.0 1 1 0.8 1 1 1.0 1 1 41.0 1 1 S.C. at IHX outlet 1 1 22.0 1 1 11.5 1 1 0.5 1 1 8.0 1 1 48.0 The transient of reactivity insertion affect FHER greatly. Itincreases thenuclear power sharply while the coolant responses slowly in temperature, but the system can return to normal level benefited from the negative temperature reactivity feedbackand large heat capacity.For other disturbed transients, the power and coolantcurvesbarely change observably and the heat transfer process of coolant in core and IHX is considerably slow. The results alsoindicate that the system possesses sufficient heat inertia and an abundant inherent self-stabilization to reach steady state, even with an insertion of disturbed transient. 4 Conclusion 1 OCFHR code is developedto analyze FHRs transient behaviors under normal operation. By comparing with other algorithm models and point reactor analytical solution, the results provethat the algorithm adopted in OCFHR is suitable for solving stiff equations. According to the five disturbed transient modes, the results fully verify the FHRs’ inherent characteristics, such as high temperature, low pressure, and large heat capacity, negative reactivity coefficients, which result in slow response and unobviousvariations during transient operation. For the reactivity step inserted mode, thenuclearpower changes promptly while the temperature of molten salt responses quite slowly. The results indicate FHRs negative reactivity feedback plays a vital role to stabilize power. For the latter flow rate shifting modes, the power’s responsesto these changes have certain delays because of salt’s large heat capacity and long circulation looping time. Besides, the IHX and SHX equipment canbuffer rapid thermal dynamic impact and increase delay-time, which can help to increase operator’s response time. Concerning the power conversion unit, the fan model adopted is not appropriate and a power cycle unit shall be developed for generation application. The research of OCFHR is to complete thecontrol and protection modulefor detailedFHR transient behaviors and suitable operation strategy. References C. Forsberg , The advanced high-temperature reactor: high-temperature fuel, liquid salt coolant, and liquid-metal reactor plant . Prog. Nucl. Energy 47 ( 1-4 ), 32 - 43 ( 2005 ). DOI: 10.1016/j.pnucene.2005.05.002 http://doi.org/10.1016/j.pnucene.2005.05.002 The advanced high-temperature reactor: high-temperature fuel, liquid salt coolant, and liquid-metal reactor plant TMSR center , SINAP . Pre-conceptual Design of a 2MW Pebble Bed Fluoride Salt Coolant High Temperature Test Reactor . ( 2012 ). Pre-conceptual Design of a 2MW Pebble Bed Fluoride Salt Coolant High Temperature Test Reactor D.F. Williams , Assessment of CandidateMolten Salt Coolants for theAdvanced High-Temperature Reactor (AHTR). ORNL-TM-2006-12 U.S. Department of energy ( 2006 ). DOI: 10.2172/885975. http://doi.org/10.2172/885975. Assessment of CandidateMolten Salt Coolants for theAdvanced High-Temperature Reactor (AHTR). ORNL-TM-2006-12 C.H. Andrades , A.T. Cisneros , J.K. Choi , et al ., Technical Description of the “Mark 1” Pebble-Bed Fluoride-salt-Cooled High-Temperature Reactor (PB-FHR) Power Plant. UCBTH-14-002 . Department of Nuclear Engineering, University of California, Berkeley ( 2014 ). DOI: 10.13140/RG.2.1.4915.7524 http://doi.org/10.13140/RG.2.1.4915.7524 . Technical Description of the “Mark 1” Pebble-Bed Fluoride-salt-Cooled High-Temperature Reactor (PB-FHR) Power Plant. UCBTH-14-002 M. Chen . Transient analysis of an FHR coupled to a helium Brayton power cycle . Prog. Nucl. Energy 83 , 283 - 293 ( 2015 ). DOI: 10.1016/j.pnucene.2015.02.015 http://doi.org/10.1016/j.pnucene.2015.02.015 . Transient analysis of an FHR coupled to a helium Brayton power cycle M. Li , J. Zhang , Y. Zou et al ., Disturbed Transient Analysis with Stable Operation Mode of TMSR-SF1. NURETH-16, Hyatt Regency Chicago, America, August 30 - September 4, on CD-ROM (2015) http://www.ans.org/store/item-700399/ http://www.ans.org/store/item-700399/ Disturbed Transient Analysis with Stable Operation Mode of TMSR-SF1. NURETH-16, Hyatt Regency Chicago, America, August 30 - September 4, on CD-ROM (2015) Q. Lv , H.C. Lin , I.H. Kim et al ., DRACS thermal performance evaluation for FHR . Annuals of Nuclear Energy 77 , 115 - 128 ( 2015 ). DOI: 10.1016/j.anucene.2014.10.032 http://doi.org/10.1016/j.anucene.2014.10.032 . DRACS thermal performance evaluation for FHR G. Ablay . A modeling and control approach to advanced nuclear power plants withgas turbines . Energy Conversion and Management 76 , 899 - 909 ( 2013 ). DOI: 10.1016/j.enconman.2013.08.048 http://doi.org/10.1016/j.enconman.2013.08.048 A modeling and control approach to advanced nuclear power plants withgas turbines C. Forsberg , L.W. Hu , P. Peterson et al . Fluoride-salt-cooled high-Temperature Reactor(FHR) for Power and Process Heat. MIT-ANP-TR-157 ( 2014 ). Fluoride-salt-cooled high-Temperature Reactor(FHR) for Power and Process Heat. MIT-ANP-TR-157 X. Huang . Research of Runge-Kutta method . Heilongjiang Sci. Tec . 23 , 23 - 27 ( 2012 ). DOI: 10.3969/j.issn.1673-1328.2012.23.063 http://doi.org/10.3969/j.issn.1673-1328.2012.23.063 (In Chinese) Research of Runge-Kutta method M.R. Galvin . System Model of a Natural Circulation Integral Test Facility . Oregon State University ( 2009 ). System Model of a Natural Circulation Integral Test Facility W. vanAntwerpen et al ., A review of correlations to model the packing structure and effective thermalconductivity in packed beds of mono-sized spherical particles . Nuclear Engineering and Design 240 , 1803 - 1818 ( 2010 ). DOI: 10.1016/j.nucengdes.2010.03.009 http://doi.org/10.1016/j.nucengdes.2010.03.009 A review of correlations to model the packing structure and effective thermalconductivity in packed beds of mono-sized spherical particles N. Wakao et al ., Effect of Fluid Dispersion coefficients on Particle-to-Fluid heat transfer coefficients in Packed Beds . Chemical Engineering Sci . 34 , 325 - 336 ( 1978 ). DOI: 10.1016/0009-2509(79)85064-2 http://doi.org/10.1016/0009-2509(79)85064-2 . Effect of Fluid Dispersion coefficients on Particle-to-Fluid heat transfer coefficients in Packed Beds S. Ergun . Fluid Flow Through Packed Columns . AIChE J , 18 ( 2 ): 361 - 371 ( 1972 ). Fluid Flow Through Packed Columns S.S. Lomakin , Yu. A. Nechaev . Transient processes and the measurement of reactivity of a reactor containing Beryllium . Translated from AtomnayaEnergiya . Vol 18 . No. 1 , January, 33-40 (1965) . DOI: 10.1007/BF01116353 http://doi.org/10.1007/BF01116353 . Transient processes and the measurement of reactivity of a reactor containing Beryllium S. Yang et al . (In Chinese) Heat Transfer . BeiJing Higher Education Press : 243 - 263 ,( 2006 ). (In Chinese) Heat Transfer M. Lin et al ., Extension of RELAP5 core power calculation model . Nuclear Power Engineering 28 , 16 - 19 ( 2007 ). DOI: 10.3969/j.issn.0258-0926.2007.06.005 http://doi.org/10.3969/j.issn.0258-0926.2007.06.005 Extension of RELAP5 core power calculation model 10.1007/s41365-017-0268-x **This work wassupported by the Strategic Priority Program of Chinese Academy of Science (XDA02010200). 2016-07-16 2016-11-15 2017-04-01 2017-07-26 2017-08-01 2021-04-11 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2017 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-017-0268-x http://dx.doi.org/10.1007/s41365-017-0268-x . 2017 Nuclear Science and Techniques 08 28