2.1 Simplified multiphysical models for neutronics, temperature, and displacement

A previous work ^[5] has proved the assumption used in early works ^{[1, 2, 7, 9]}that in the Godiva I reactor, the spatial distribution of the fission rate is independent from the temporal variation. In one–group isotropic scattering theory, the point kinetic approximation for neutronics can be used. The point kinetic equation for the fission rate and the balance equation for the delayed neutron precursors are as follows: ^[10]

where ρ is the reactivity, β_eff is the effective delayed neutron fraction, Λ is the neutron generation time, ν_f is the average number of neutrons generated per fission, λ is the decay constant, F is the fission rate, C is the precursor density, t is the time, and the subscript i is the group number of the precursors.

The burst duration is typically ignorable compared to the heat balance time such that both the heat conduction in the core and the heat exchange through the surface can be neglected, i.e. the adiabatic approximation can be applied as follows:

where ${\rho }_{\text{V}}$ is the density, c_p is the specific heat capacity, T is the temperature rise, ${\epsilon }_{\text{f}}$ is the average energy released per fission, and r is the location in the radial coordinate system.

With the limitation in purely elastic behavior, i.e. the operation of the FBR is under a controlled condition, the displacement in the core can be calculated using the 1D elastic wave equation as follows: ^[10]

where the Lame constants are defined as follows:

Similar to the literature ^[9], the free boundary condition is applied on the core surface as follows:

where u is the displacement, α is the linear thermal expansion coefficient, v is the Poisson ratio, E is the Young’s modulus, and R is the core radius.

2.2 Numerical algorithms

To obtain the temporal variations and/or spatial distributions of the fission rate, temperature rise, and displacement, proper numerical algorithms should be adopted to calculate Eqs. (1)– (4).

In this study, the space–dependent terms in Eqs. (1)– (4) were discretized using the finite difference method. Considering that the variation in the fission rate will result in changes in the temperature rise and displacement, which in turns leads to variations in reactivity and fission rate, the iteration calculations of fission rate, temperature rise and displacement should be made at each time step. Therefore, the implicit scheme algorithm was applied on the time–dependent terms in Eqs. (1)– (4), i.e. the obtained fission rate from Eqs. (1) and (2) is provided as the input of the adiabatic Eq. (3) to calculate the temperature rise needed in the elastic wave equation (4). Then, the reactivity in Eqs. (1) and (2) is updated by using the developed displacement–reactivity feedback model and the displacement obtained in Eq. (4). Such an iteration carries on until reaching the convergence criterion and then moves to the coupled calculations at the next time step.

2.3 Traditional fission yield (temperature rise)–reactivity feedback models

For a burst without the inertia effect, the core expansion can accommodate the temperature rise, the distribution shape of the displacement remains the same as that in the static expansion caused by the same temperature rise, and only the magnitude of the displacement increases in proportion to the temperature rise ^{[1, 2]}, i.e.

where u₀ is the displacement for the average temperature rise $\overline{T}$ of 1 K in the static expansion.

For a spherical FBR, the relation between the reactivity change and displacement can be obtained by using the perturbation theory as follows:

where, $\Delta \rho$ is the reactivity change and W is the function of perturbation coefficient.

Substituting Eq.(7) into Eq.(8) and defining the constant coefficient $b_{\text{T}}$ we obtain the following:

The reactivity change is thus expected in a linear relation with the temperature rise (or the equivalent fission yield) and can be calculated using the Fuchs–Hansen model ^[9] as follows:

where, b_T is the temperature–reactivity feedback coefficient and b_f is the quench coefficient defined as the absolute value of the reactivity change per fission.

When the inertia effect is generated, a lag between the core expansion and temperature rise occurs. The core appears compressed as compared to the static expansion caused by the same temperature rise, i.e. the realistic displacement $u\left(r,t\right)$ is smaller than $u_0\left(r\right)\overline{T}\left(t\right)$. Thus, the application of Eq. (10) will result in an overestimation of the reactivity decrease and underestimations of the peak fission rate and fission yield. To account for the impact of the inertia effect, a widely used method is to establish the relation between the reactivity change and fission yield with the assumption of single–frequency core vibration (only the contribution of the fundamental mode is considered and the vibrational period is also the same as the fundamental frequency) such as depicted by the following: ^[10]

where α₀ is the reciprocal of the reactor period and τ is reciprocal of the fundamental angular frequency of the core vibration.

The total reactivity can be determined as follows:

where, ${\rho }_0$ is the initial reactivity.

2.4 Main parameters of the burst calculated using traditional reactivity feedback models

It is known that the performance of an FBR is primarily characterized by the peak fission rate and fission yield under controlled conditions. By neglecting the delayed neutron precursors in Eq. (1) (the burst time is very short as compared to the half–life such that the precursors nearly do not decay), analytical solutions for the fission rate F and fission yield Q can be calculated using the traditional reactivity feedback models in Section 2.3 as follows: ^{[1, 2, 10]}

where δ can be set as 0 or 1 when neglecting or considering the inertia effect and t_max is the time of peak fission rate defined as follows:

With the adiabatic approximation, the average temperature rise can be obtained as follows:

The displacement can be estimated as follows: ^[11]

where m is the material mass of the core.

By comparing Eqs. (10) and (11) it is clear that the impact of the inertia effect on the reactivity change is considered by modifying the quench coefficient $b_{\text{f}}$ with a constant factor of 1/(1+α₀²τ²), which indicates that the fission rate/yield, temperature rise, and displacement are (1+α₀²τ²) times those neglecting the inertia effect.

4.1 Validation of the static neutron transport component of the FBR–MPC code

The static neutron transport calculation was performed for the spherical model of the Godiva reactor and the calculated effective multiplication factor k_eff was 0.9978, very near the measured 1.000±0.001^[14] and 0.99629±0.00009 obtained using the JMCT code^[15]. Because of the lack of experimental data, the calculated neutron flux distribution was only compared to that obtained using JMCT as shown in Fig. 1. The two curves nearly overlap, indicating the validation of the static neutron transport component of the FBR–MPC code.

Fig. 1Full-screen View

Radial distribution of the normalized static neutron flux

For the Godiva I reactor studied, with the application of the point kinetic approximation for neutronics and the adiabatic approximation for temperature, the fission rate and temperature rise in the burst have the same spatial distribution shape as the static neutron flux and can be expressed as a function in the form of $\sin \left(r/R\text{'}\text{'}\right)/r$, where $R\text{'}\text{'}\text{=}10.65cm$ is the extrapolation radius determined by the static neutron flux distribution calculated using the FBR–MPC. The fitted function of $\sin \left(r/R\text{'}\text{'}\right)/r$ is also shown in Fig. 1 and is consistent with the static neutron fluxes obtained using the FBR–MPC and JMCT codes.

4.2 Main burst parameters

Bursts triggered with different initial stepwise reactivity (0.65, 1.08, 2.09, 3.29, 5.99, and 8.37¢ above supercritical) were simulated using the simplified multiphysics models described in Section 2.1 and the displacement–reactivity feedback model described in Section 3. The results were compared to the experimental data and the solutions described in Section 2.4 using traditional fission yield (temperature rise)–reactivity feedback models described in Section 2.3.

As can be seen from Tables 2 and 3, for all selected bursts, the peak fission rates calculated using the developed displacement–reactivity feedback model and Eq. (13) using the experimental quench coefficient $b_{\text{f}}=2.08E-17\text{\$/}f$ correspond very closely to the values of the experimental data. In addition, for bursts without initial reactivity (for instance ρ₀<2.09 ¢), the calculated fission yields are larger than the experimental yields whereas the temperature rises are smaller than the experimental ones. This can be attributed to the heat capacity of the core employed in this study being higher than the realistic capacity. A larger heat capacity causes a smaller temperature rise and core expansion, and its induced reactivity feedback effect are consequently smaller than the experimental values. As a result, the peak fission rate and fission yield tend to be overestimated.

For bursts with large initial reactivity (such as ρ₀=8.37 ¢), because the burst width in the traditional models remains same as that upon neglecting the inertia effect (see Eq. (13)), the temporal integration of fission rate, i.e. the fission yield, is overestimated. It can thus be expected that the fission yield obtained using the displacement–reactivity feedback model will be closer to that obtained using the experimental data, because the smaller burst width and hence the smaller temporal integration of the fission rate caused by the inertia effect can be revealed. However, the calculated results deviate from the experimental data, as shown in Table 3. This deviation was also noticed in our previous work ^[5] and can be attributed to the neglection of the damping effect in Eq. (4). Notably, although the traditional models seem to provide more accurate results for fission yields for bursts with a large initial reactivity, they cannot reflect the realistic physical process of a burst because this coincidence is caused by the fact that the offsets of inertia and damping effects both lead to a small change in burst width. Moreover, the underestimation of fission yield results in a smaller temperature rise and surface displacement as compared to those of the experimental data and/or those obtained using Eqs. (16) and (17) with traditional models, as seen in Tabs. 4 and 5.

To obtain intuitive insight into the relation between the reactivity change and displacement, a new reactivity feedback coefficient b_s, defined as the absolute value of the reactivity change per unit change in surface displacement, is introduced, i.e. b_s=|Δρ|/u_R.

For bursts with ρ₀ ≤ 2.09 ¢, the results show that the inertia effect is not generated, i.e. the core expansion can accommodate the temperature rise, the shape of the displacement distribution remains the same as that in the case of static expansion, and only the magnitude of displacement changes. In such cases, b_s is equivalent to the quench coefficient b_f and they are in a linear relation of b_f=α_fsb_s, α_fs= 6.02E–17 cm/f for the Godiva I reactor. As shown in Fig. 2, b_s remains constant at 36.21$/cm, corresponding to b_f=2.18E–17$/f, very near the experimental value of 2.08E–17$/f ^[6].

Fig. 2Full-screen View

Temporal variations of b_s for various bursts (0 represents time of t_max)

As ρ₀ increases, a lag between the core expansion and temperature rise occurs because of the faster increase in the fission rate/power, the curve for the distribution shape of displacement is under those of bursts with ρ₀ ≤ 2.09 ¢, and the core appears to be compressed, as shown in Fig. 3. Thus, the application of Eq. (10) will underestimate both the peak fission rate and fission yield because the core expansion and its induced reactivity decrease are overestimated. To account for the inertia effect, early works modified the quench coefficient $b_{\text{f}}$ such as in Eq. (11), whereas the modified quench coefficient $b_{\text{f}}^{\text{"}}$ remains constant, indicating that the shape of the displacement distribution remains the same as that in the bursts without an inertia effect, as in Eq. (17). This is in conflict with the fact that when the inertia effect is generated, the shape of the displacement distribution varies in the burst, as shown in Fig. 4. As a result, the impact of the variation in the distribution shape of the displacement on neutronics cannot be shown. Fortunately, this difficulty can be overcome using the developed displacement–reactivity feedback model.

Fig. 3Full-screen View

Radial distribution of normalized displacement during increasing periods of bursts

Fig. 4Full-screen View

Radial distributions of normalized displacement for a burst with ρ₀=8.37 ¢

As shown in Fig. 2, for bursts with ρ₀>2.09 ¢, at the initiation of the burst, because the fission power is very small, the temperature rise and its induced core expansion can be neglected. Thus, b_s remains practically constant. As the fission rate intensely increases, the temperature rise and core expansion vary sharply, leading to a rapid increase in b_s. For a burst with ρ₀ ≥ 3.29 ¢, oscillations caused by the core vibration, noticed in other experiments ^{[6, 7]} and our previous work ^[5] but which cannot be reproduced using the traditional models, can be observed. Notably, for a burst with ρ₀=3.29¢, b_s finally does not tend to attain a constant value for the burst without an inertia effect but oscillates slightly with time because of the small inertia effect.

Moreover, using the displacement–reactivity feedback model, the more rapidly increasing velocity of the fission rate at the initiation of the burst, which has been noticed in the CFBR–II reactor ^[7] but cannot be reproduced using Eq. (13), can be observed. This is because the realistic reactivity feedback coefficient $b_{\text{s}}$ is smaller than $b_{\text{s}}^{\text{"}}$ (corresponding to $b_{\text{f}}^{\text{'}}$ in Eq. (11), forming a straight line below that for the burst with ρ₀ ≤ 2.09 ¢ but still above the initial part of $b_{\text{s}}$ in Fig. 2). In addition, during the drop period of the fission rate (t>0 in Fig. 2), the larger b_s indicates a faster changing ratio of fission rate than during the increasing period (t<0), leading to an asymmetry of the fission rate curve and smaller burst width (also see Fig. 6).

Fig. 6Full-screen View

Temporal variations in the fission rate of various bursts

4.3 Temporal variations

Because of the lack of detailed experimental data, only the calculated temporal variations in reactivity, fission rate, average temperature rise, and surface displacement are presented and analyzed. The bursts with ρ₀ ≤ 1.08 ¢ are not considered because in such cases these parameters are small (see Tabs. 2–5) and their temporal variations are very slow.

As shown in Figs. 5 and 6, for a burst without an inertia effect (ρ₀=2.09 ¢), the calculated results are similar to those obtained using Eqs. (10) and (13), i.e. the reactivity decreases from an initial value of ρ₀ to –ρ₀. This change is terminated during the drop period of the fission rate; the fission rate nearly increases/decreases in an exponential manner and its curve is symmetric about the peak.

Fig. 5Full-screen View

Temporal variations in the reactivity of various bursts

For bursts with the inertia effect, the reactivity does not monotonically decrease to –ρ₀ as expected in early theoretical works ^[1] but a temporal oscillation is generated. For the bursts with ρ₀=5.99 ¢ and 8.37 ¢, the changing ranges of reactivity are –6.00±1.20 ¢ and –9.74±5.51 ¢, respectively, in the same period of core vibration (see Figs. 6 and 8). As the initial reactivity ρ₀ increases, the delay between the core expansion and temperature rise during the initial period of the burst becomes more significant (the smaller $b_{\text{s}}$ in Fig. 2), which subsequently results in the more prominent core vibration and hence the larger oscillation amplitude of reactivity. In addition, the average reactivity is smaller than –ρ₀, which indicates that the FBR system will finally be in a deeper subcritical state rather than the recognized –ρ₀; this tendency is consistent with the experimental data of the CFBR–II reactor ^[7]. Moreover, because the reactor core changes from expansion to compression during the drop period of the burst (see Fig. 8), the reactivity changes from a decrease to an increase and the decrease ratio of the fission rate thus becomes smaller and the fission rate curve is distorted, as shown in Fig. 6. This tendency was first noticed in the literature ^[5] and was reconfirmed in this study.

Fig. 8Full-screen View

Temporal variations in surface displacement for various bursts

The temporal variations in average temperature rises are shown in Fig. 7. The temperature sharply increases only when the fission rate (or the corresponding fission power) is sufficiently high (F > 1E17/s). At the beginning of the plateau following the burst, the fission power (the corresponding fission rate is in the magnitude of 1E16/s, see Fig. 6) is much larger than the heat loss ratio through the core surface given the heat transfer and thermal radiation effects. Therefore, the temperature rises at a very slow rate. Notably, if the reactor core is not scrammed, for bursts with a prominent inertia effect, the increasing temperature will lead to a stronger core vibration and the stress also increases and may exceed the yield strength, resulting in material failure.

Fig. 7Full-screen View

Temporal variations in the average temperature rise of various bursts

For a burst without an inertia effect (ρ₀=2.09 ¢), the core expansion can accommodate the temperature rise and the displacement at an arbitrary position of the core increases in the same manner as expected using Eq. (7). The surface displacement monotonically increases and tends to a constant value such as that in the static expansion. Whereas for a burst with an inertia effect, the core finally converts to an undamped free vibration, and the vibration amplitude increases for bursts with a larger initial reactivity, as shown in Fig. 8. The calculated vibration period is 60 μs, consistent with that obtained performing an accurate 3D dynamic coupled multiphysical simulation using the FBR–MPC code ^[5]. Because of the underestimation of fission yield for a burst with ρ₀=8.37 ¢ using the developed displacement–reactivity feedback model, the calculated maximal equivalent stress is ~121 MPa on the core surface, smaller than that obtained using the dynamic simulation of the FBR–MPC code.