2.1. Algorithm

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]. Besides, Gear algorithm, a variable step length algorithm, is programmed in comparison to R-K method on computation speed and accuracy.

2.2. Simplification

The FHER system consists of the pebble bed, cavities, pipes, pumps, heat exchangers, loadetc. 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]. 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 represents simplifications of the pebble bed. Relevant parameters of FHER are givenin Table 1.

Fig. 3Full-screen View

Node diagram of FHER

Fig. 4Full-screen View

(Color online) Simplifications of Pebble bed

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].

where Nu is the Nusselt number, Reis the Reynolds number, Pr is the Prandtl number, and D is the equivalent diameter. Eq.(1) is applicable for 15<Re<8500[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

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.

where V is velocity vector; ρis the density of coolant[kg/m³]; anduis 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,

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, andP is pressure head of pump.Therefore, the integral form of Eq.(4) can be written as,

where m is the number of pipes or other sections, ΔP_f is the friction pressure drop,and ΔP_Kis 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,

The heat transfer form in pebble graphite shell can be calculated by:

where P_f is the reactor power[W];F is the power axial distribution factor; c_pfm_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_gis the external graphite temperature in a node [°C]; c_pgm_g is the specific heat capacity of the node [J/°C]; h_gsis equivalent heat transfer coefficient between graphite and molten salt [W·m⁻²·°C⁻¹]; 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 shows the heat transfer process in the core.

Fig. 5Full-screen View

(Color online) Heat transfer process in the core

In aprimary coolant node of core,

where m_s is the coolant mass in one node [kg];c_psis the coolant specific heat [J·kg⁻¹·°C⁻¹]; ${\dot{m}}_s$ is the coolant mass flow rate [kg/s]; andT_out andT_in are molten salt outlet and inlet temperatures, respectively,in one node.

Inagraphite reflectornode,

where c_pge is the graphite specific heat [J·kg⁻¹·°C⁻¹]; 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⁻²·°C⁻¹];h_gd is the equivalent heat transfer coefficient between salt in the downcomer and reflector [W·m⁻²·°C⁻¹]; andA_gd is the effective heat transfer area.

In a coolant node of downcomer,

where m_sd is the coolant mass in one node [kg]; ${\dot{m}}_{\text{s}}$ is the coolant mass flow rate [kg/s]; and ${\dot{P}}_0$ 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]:

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:

where n is the neutron density [n/m³]; 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⁻¹]; l is the mean neutron generation time [s]. β, λ and l are gained from neutron physics calculation while the 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].

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.

whereρ₀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; $\Delta {\overline{T}}_{fg}$, $\Delta {\overline{T}}_s$ and $\Delta {\overline{T}}_{ge}$ 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 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.

Fig. 6.Full-screen View

Schematic diagram of IHX model

For the tube side,

where c_ps is the internal salt specific heat [J·kg⁻¹·°C⁻¹]; m_ep is internal coolant mass in a node [kg];T_ep is internal coolant temperature [°C];T_edis the intermediate wall temperature [°C]; h_ep is the equivalent heat transfer coefficient between salt in the tube and wall [W·m⁻²·°C⁻¹];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,

where c_ss is the external salt specific heat [J·kg⁻¹·°C⁻¹]; 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⁻²·°C⁻¹]; 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].

For the tube side

For the shell side

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.1. Steady-state operation

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, 4 and 5, respectively.

3.2. Model solver verification

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.Full-screen View

Reactor power level behavior following negative (a) and positive (b) stepreactivity insertion based on different algorithm models.

In 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), R-K numerical solution matches well with single-group analytical solutionand the truncation error can be neglected[17].

The thermal hydraulic calculation is based on lumped parameter method. A simple case (Fig.8 and 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.

Fig. 8Full-screen View

(Color online) Loop structure of the case

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 prove the validity of thermal hydraulic calculation methodand its reliability tosimulate transient behavior of the FHER.

Fig. 9Full-screen View

(Color online) Simulation results of temperaturesof (a) the heater and cooler and (b) the salt, forthe case shown in Fig.8 and Table 6.

Fig. 10Full-screen View

(Color online) Relative temperature errors of (a) the heater and cooler and (b) the salts

3.3. Transient simulation

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 inTable7.

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

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) 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) 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.

Fig. 11Full-screen View

(Color online) Reactor power evolution (a)and coolant temperatures(b) after reactivity step insertion.

3.3.2 10% of primary coolant mass flow rate step increases

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) 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) 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.

Fig. 12Full-screen View

(Color online) Reactor power evolution (a) and coolant temperatures (b) after mass flow change.

3.3.3 10% of secondary coolant mass flow rate step increases

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) 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), 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.

Fig. 13Full-screen View

(Color online) Reactor power evolution (a) and coolant temperatures (b) after mass flow change.

3.3.4 20% of air mass flow rate step increases

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) 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) 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, 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.

Fig. 14Full-screen View

(Color online) Reactor power evolution (a) and coolant temperature (b) after load change

3.3.5 10°C of air generator inlet temperatureincreases

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), 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.

Fig. 15Full-screen View

(Color online) Reactor power evolution (a) and coolant temperatures (b) after air temperature change.