The 3-D time-dependent transport code TORT-TD and its coupling with the 3D thermal-hydraulic code ATTICA3D for HTGR applications

Nuclear Engineering and Design 251 (2012) 173–180

Contents lists available at SciVerse ScienceDirect

Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes

The 3-D time-dependent transport code TORT-TD and its coupling with the 3D thermal-hydraulic code ATTICA3D for HTGR applications A. Seubert a,∗ , A. Sureda a , J. Bader b,c , J. Lapins b , M. Buck b , E. Laurien b a b c

Gesellschaft für Anlagen- und Reaktorsicherheit (GRS) mbH, Forschungszentrum, Boltzmannstraße 14, D-85748 Garching, Germany Institut für Kernenergetik und Energiesysteme (IKE), Universität Stuttgart, Pfaffenwaldring 31, D-70569 Stuttgart, Germany EnBW Kernkraft GmbH, Kernkraftwerk Philippsburg, Rheinschanzinsel, D-76661 Philippsburg, Germany

a r t i c l e

i n f o

Article history: Received 23 February 2011 Received in revised form 24 September 2011 Accepted 24 September 2011

a b s t r a c t The application of the time-dependent 3-D discrete-ordinates neutron transport code TORT-TD to HTGR of pebble bed type and its coupling with the 3-D thermal-hydraulic code ATTICA3D is described. TORTTD, developed at GRS, solves the time-dependent multi-group SN transport equation in both Cartesian and curvilinear geometry using a fully implicit time discretization scheme. To provide a detailed 3-D core model for pebble bed type HTGRs, TORT-TD has been coupled with the 3-D thermal-hydraulic code ATTICA3D which is developed at IKE, University of Stuttgart. ATTICA3D provides a 3-D modeling of heat conduction and convection using a porous media approach. In this paper, the TORT-TD/ATTICA3D coupling methodology is described. The PBMR-268 and PBMR-400 benchmark problems have been chosen as first test cases for TORT-TD and the coupled code TORT-TD/ATTICA3D. In TORT-TD, the control rods in the side reflector have been modeled both as grey curtain and spatially resolved. Results for first coupled steady-state and transient calculations of the PBMR-400 benchmark are shown. © 2011 Elsevier B.V. All rights reserved.

1. Introduction It is known that deterministic neutronics, thermal-hydraulics and transient analysis tools and methods available to design and analyse High-Temperature Gas-cooled Reactors (HTGR) have, in many cases, lagged behind the state of the art compared to other reactor technologies, e.g. light water reactors (LWR). For comprehensive safety studies of HTGR, full three dimensional coupled treatments of both neutron kinetics and thermal-hydraulics become more important to account for non-symmetric effects. As the neutron transport code TORT-TD and the thermal-hydraulics code ATTICA3D are both three dimensional steady-state and transient codes, the coupled code system TORT-TD/ATTICA3D may provide a basis for high-accuracy full three dimensional coupled neutronics and thermal-hydraulics simulation for HTGR of pebble bed type. 2. The time-dependent 3-D discrete ordinates transport code TORT-TD TORT-TD (Seubert et al., 2008) is a time-dependent 3-D multigroup discrete ordinates (SN ) neutron transport code developed at GRS. It is based on the DOORS steady-state neutron transport

∗ Corresponding author. Tel.: +49 89 32004 469. E-mail address: [email protected] (A. Seubert). 0029-5493/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.nucengdes.2011.09.067

code TORT (Rhoades and Childs, 1991; Rhoades and Simpson, 1991) and solves the steady-state and time-dependent multigroup transport equation with an arbitrary number of prompt and delayed neutron precursor groups in both Cartesian and cylindrical (r–ϑ–z) geometry. Unconditional stability in transient calculations is achieved using a fully implicit time discretization scheme. Scattering anisotropy is treated in terms of a Pl Legendre scattering cross section expansion. Computing time can be saved by extrapolating the angular fluxes to the next time step using the space-energy resolved inverse reactor period ωg (r). In order to reduce spatial homogenization errors at the level of pin cells (for LWR applications), Generalized Equivalence Theory (GET) has recently been implemented in TORT-TD in terms of pin cell discontinuity factors (PDF) (Seubert, 2010). TORT-TD has also been extended to account for time-dependent anisotropic distributed external sources (Seubert et al., 2009), an issue that becomes relevant for subcritical systems driven by external neutron sources. For treating thermal-hydraulic feedback, TORT-TD has been coupled with the GRS system code ATHLET (Lerchl and Austregesilo, 2003) and the subchannel code COBRA-TF (Christienne et al., 2010). Few-group macroscopic cross sections are fed into TORT-TD in terms of parameterized tabulated cross section libraries where the dependence of up to five thermal hydraulic parameters can be considered. Depending on the current thermal-hydraulic state in each spatial mesh, appropriate cross sections are interpolated by TORT-TD between given sampling points either linearly or using cubic spline polynomials. For handling buckling dependency of

174

A. Seubert et al. / Nuclear Engineering and Design 251 (2012) 173–180

LMW-Benchmark 300 TORT-TD

Power (MW)

250 QUABOX/CUBBOX 200 150 100 50 0,0

10,0

20,0

30,0

40,0

50,0

60,0

Time (sec) Fig. 1. TORT-TD solution (red line) of the LMW benchmark. The blue line denotes the QUABOX/CUBBOX reference solution. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

few-group cross sections, TORT-TD calculates the buckling over larger spatial regions, e.g. material zones, by evaluating the net leakages across the surfaces of the corresponding volumes. By implementing steady-state and transient Iodine–Xenon equations, TORT-TD has been prepared for the simulation of operational transients. TORT-TD, both stand-alone and coupled to ATHLET or COBRATF, has been applied to several, mainly LWR problems, and verified against transient benchmark problems. From the two examples shown here, the first one refers to the 3-D LMW benchmark (Langenbuch et al., 1977) which simulates two control rod groups moving in opposite directions in a simplified PWR, thereby initiating a delayed critical transient in the range of minutes. The total power evolution is shown in Fig. 1. The second example is the Dodds benchmark, a delayed supercritical transient in 2-D cylindrical geometry (ANL, 1985). Again, the TORT-TD solution agrees well with the reference solution regarding both the total power evolution (Fig. 2) and the radial and axial power density distributions (upper and lower panel of Fig. 3). 3. The 3-D thermal-hydraulic code ATTICA3D 3.1. Physical models and numerical methods The Advanced Thermal hydraulic Tool for In-vessel and Core Analysis in 3 Dimensions (ATTICA3D) is a three-dimensional tool

Dodds-Benchmark 2,8 2,6

for thermal hydraulic analysis of high temperature reactors and was initially introduced as TH3D (Hossain et al., 2008a). ATTICA3D supports both 3-D Cartesian and cylindrical geometry. Gas and solid phases are both considered as continua that occupy their respective volume fractions. For flow calculations, the porous medium approach is applied. This leaves out detailed description of the components; considered control volumes are subdivided into gas and solid fraction via the porosity parameter ε. Important flow parameters like the hydraulic diameter are input by the user. The coolant helium is a single-phase gas with variable density. For nominal operation, only one gas phase has to be accounted for, but provision is made for further gas components. Between solid and fluid phase, thermal non-equilibrium is assumed and the model is based on the following set of conservation equations: The time dependent energy equation is given by: (1 − ε)s cp,s

2,4

Relave power

Fig. 3. Radial (top) and axial (bottom) power density distributions for the Dodds benchmark obtained by TORT-TD (solid lines) at different instances in time compared to the ANL reference (dashed lines).

2,2 2 1,8 1,6

TORT-TD 3-D transport soluon

1,4

ANL Diffusion soluon

1,2 1 0,0

1,0

2,0

3,0

4,0

5,0

Time (sec) Fig. 2. TORT-TD power evolution (red) for the Dodds benchmark compared to the ANL reference (blue). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

∂Ts = (1 − ε)div(eff grad(Ts )) − q˙ conv + q˙ nu ∂t

(1)

where index s indicates the solid state, ε, , cp , eff , T, q˙ conv and q˙ nu denote porosity, density, heat capacity, effective conductivity, solid temperature, convective heat transferred from hot gas to solids, and the amount of generated nuclear power, respectively. The time dependent energy equation for the gas phase is solved for the flow region: ∂εh  h) = div(εg,eff grad(Tg )) + q˙ conv + div(εu ∂t

(2)

 denotes the velocity vector, where h is the specific gas enthalpy, u g,eff the effective heat conductivity of the gas, Tg the gas temperature and q˙ conv the heat transferred from the solid to the gas phase. The latter is calculated by q˙ conv = ˛(Ts − Tg ), where ˛ is the

A. Seubert et al. / Nuclear Engineering and Design 251 (2012) 173–180

(3)

where cn denotes the respective concentration of the gas component (if more than one is present, e.g. air ingression, n = 1, . . ., N gas components). A quasi-steady, simplified momentum equation is solved in the flow region. Simplifications are made since major friction and gravity body force are the dominating driving mechanisms, while other terms contribute only little, thereby saving computing time:

800 ATTICA3D

700

Benchmark upper lower

500 1100 1000

(4)

 , g denote the pressure, the friction force in the porous where p, R, medium, a parameter determining the algebraic sign, and the acceleration due to gravity.  is equal to +1 for upward flow and equal to −1 for downward flow. The set of conservation equations is transformed into a set of initial value ordinary differential equations. Spatial derivatives are replaced by discretisation for which the finite-volume method is applied. For the calculation of velocities and all other variables, a staggered grid approach is applied. Velocities are calculated at the interface of the control volumes while the other variables are calculated at the centre. For time integration, a fully implicit, time adaptive multi-step backward differentiation method (BDF) is used. The non-linear equations produced by BDF are subsequently solved for each time step by a modified Newton method together with sparse matrix techniques. In order to correctly capture the feedback from thermal hydraulics on neutronics, a heterogeneous fuel model that distinguishes between fuel kernel temperatures and surrounding graphite matrix temperatures is implemented in ATTICA3D. In this model the pebble can be subdivided in an arbitrary number of shells. While the outermost shell contains graphite only, the inner shells contain both graphite matrix and coated fuel kernels. For the pebble’s outermost shell the surface temperature serves as boundary condition. From there on the heat conduction is solved towards the pebble centre. In the heterogeneous approach a representative coated particle is determined in each shell. It is clear that the hottest spot will be at the very centre. For the processing of cross sections, these more detailed temperatures are employed.

800

700

600

500

400

300

200

100

0

Average radial fuel temperature 900

850

800

ATTICA3D Benchmark upper lower

750

700 100

110

120

130

140

150

160

170

180

190

Radius [cm] Average Axial Coolant Temperature in Core Region 1000

900

800

700

ATTICA3D Benchmark upper lower

600

500 1000

800

600

400

200

0

Axial Height [cm]

Average Axial Maximum Fuel Temperatures for Fuel Kernels and Inner Moderator Shell

3.2. Comparison of ATTICA3D to other codes

1100

Temperature [°C]

With ATTICA3D several calculations were performed for reasons of comparison to other code systems (Hossain et al., 2008a,b), e.g. the IAEA CRP-3 benchmark for the GT-MHR (Ball, 1997a,b; Kuzavkov, 1997) and the OECD/NEA/NSC PBMR benchmark (Reitsma et al., 2007). Therein, most participants applied particular versions of the well-recognized THERMIX/KONVEK code. Selected results for the PBMR-400 benchmark obtained with ATTICA3D are shown in Fig. 4. The thermal hydraulic exercise 2, a purely thermal hydraulic 2-dimensional exercise, was performed with ATTICA3D using the predefined power distribution provided by the benchmark. The obtained results for the PBMR-400 are compared to the published results (Reitsma et al., 2008), and demonstrate the computational capabilities of our new 3-D code. The distribution of temperatures coincides well for the axial distribution with only minor (∼5 K) deviations from the scatter band of participants’ results for the radial temperature distribution. Keeping in mind that all but one participant used their particular THERMIX versions, the agreement

900

Height [cm]

Temperature [°C]

 − εg εgrad(p) = −R

900

600

Temperature [°C]

∂(εcn )  cn ) = 0 + div(εu ∂t

Axial average fuel temperature

1000

Temperature [°C]

volumetric heat transfer coefficient. Values for ˛ for pebble beds are determined according to KTA standards (KTA, 1983). The time dependent, compressible mass conservation equation for n gas components is solved for the flow region:

175

900

ATTICA3D- fuel kernel ATTICA3D- moderator shell 700

Benchmark upper lower

500 1000

800

600

400

200

0

Height [cm] Fig. 4. PBMR-400 benchmark results (averaged values) for axial and radial fuel temperatures, axial coolant temperatures and axial fuel temperatures derived from innermost moderator shell temperatures (mixture of fuel kernels and graphite matrix) together with the maximum fuel kernel temperatures (from top to bottom) obtained with ATTICA3D.

176

A. Seubert et al. / Nuclear Engineering and Design 251 (2012) 173–180

is present if the 3-D core model covers the neutronics including the fuel temperature model and the fluid-dynamics of the core region, whereas the thermal-hydraulics code treats the whole thermal fluid-dynamics except for the core region (Fig. 5). For parallel coupling, the 3-D neutronics code – including the fuel temperature and the fluid-dynamics model – represents the reactor core. The thermal-hydraulics code simulates the thermal fluid-dynamics in the primary circuit and the core region in a simplified manner. The calculated boundary conditions of the thermal-hydraulics code are transferred in terms of time-dependent boundary conditions of the more detailed core calculation performed in parallel. The advantage of the internal coupling approach is that the fluid-dynamic equations can be integrated in closed form. This may be relevant for temperature and density transport processes. The coupled code system TORT-TD/ATTICA3D is represented by a single executable in which ATTICA3D acts as the main program and calls TORT-TD in terms of a subroutine whenever an update calculation of the power distribution is requested. For the data exchange between TORT-TD and ATTICA3D, already existing TORT-TD interface routines have been utilized in combination with the ATTICA3D mesh overlay feature that transfers 3-D distributions from its thermal-hydraulic mesh to a superimposed neutron-kinetics mesh and vice versa. This allows for efficient data transfer via direct memory access of array elements. 5. First 3-D test cases

Fig. 5. Types of coupling between a thermal-hydraulics and a 3-D neutron kinetics code. In TORT-TD/ATTICA3D, the topmost coupling approach is implemented. For the steady state calculation, ATTICA3D and TORT-TD are called repeatedly, followed by exchange of thermal-hydraulic and neutron kinetics data, until convergence of the 3-D temperature and power distributions are achieved. At the beginning of the iteration process, TORT-TD calculates for a given thermal-hydraulic initial distribution the corresponding power distribution that is transferred to ATTICA3D as first estimate.

of temperature distribution can be regarded as excellent for a code to code benchmark. It is noteworthy that the definition of the maximum fuel temperatures in the benchmark description is somewhat ambiguous. It was defined as the innermost shell temperature of a homogenized fuel/moderator mixture. But for failure of fuel particles it is more important to regard the hottest spots which will be found in the fuel kernels in the very centre of the pebble. To visualise the differences in maximum temperatures, innermost moderator shell are compared to the maximum kernel temperatures in Fig. 4. However, for further validation code to experiment benchmarks are crucial and will be carried out in the future. The ATTICA3D results shown here are subsequently taken as basis to establish a coupling with a 3-dimensional neutronics transport code as described in the following chapter. 4. Coupling ATTICA3D with TORT-TD For comprehensive safety studies of pebble bed type HTR, a full three dimensional treatment of both neutron kinetics and thermalhydraulics is mandatory. This is achieved by coupling TORT-TD with the 3-D thermal-hydraulic code ATTICA3D (Hossain et al., 2008a). The so-called internal coupling methodology has been adopted. The coupling is called internal if the thermal-hydraulics code models the entire thermal fluid-dynamics of the system including the core region. This is indicated in the top panel of Fig. 5. External coupling

In this paper, two pebble bed type HTR test cases derived from the PBMR-268 (Reitsma et al., 2005, 2006) and the PBMR-400 benchmark (Reitsma et al., 2007; Strydom et al., 2010) are considered. GRS has successfully participated in these benchmarks with the 2-D coupled code system DORT-TD/THERMIX (Pautz et al., 2008; Tyobeka et al., 2008). The applications reported in this paper, however, focus on 3-D models of both neutron kinetics and thermal hydraulics. 5.1. 3-D modeling of PBMR-268 in TORT-TD The PBMR-268 neutronics and transient benchmark problem has been designed with the objective to assess existent deterministic neutronics, thermal-hydraulics and transient analysis tools and methods available to design and analysis of the Pebble-Bed Modular Reactor (PBMR) as a High-Temperature Gas-cooled Reactor (HTGR) concept. For the PBMR-268 benchmark, several geometrical simplifications have been made to transform the core design into an equivalent 2-D r–z model. This in particular applies to the control rods in the side reflector that are modeled as a cylindrical skirt, also referred to as a grey curtain. In this paper, the PBMR-268 has been modeled in TORT-TD in two different ways: First according to the grey curtain model, i.e. with no azimuthal dependence, in order to verify TORT-TD with an existing DORT-TD solution. The second model is a full 3-D model with the 24 control rods in the side reflector spatially resolved. Both calculations are based on a fixed nuclear cross section set in 4 energy groups and P3 Legendre scattering cross section expansion order obtained with the MICROX-2 code (Mathews, 1997). A levelsymmetric S4 quadrature has been applied in TORT-TD. For the grey curtain model, an eigenvalue of keff = 0.9262 is obtained. This is in very good agreement with a corresponding DORT-TD calculation (keff = 0.9265) and thus regarded as a confirmation of the TORT-TD model. In order to spatially resolve the 24 control rods in the second TORT-TD model, the azimuthal domain was divided into 49 sectors as depicted in Fig. 6. With 56 radial and 66 axial meshes, this amounts to a total of 181,104 meshes. Since it is currently not

A. Seubert et al. / Nuclear Engineering and Design 251 (2012) 173–180

Fig. 6. Radial cross section of the TORT-TD model for the PBMR-268 with spatially resolved control rods.

possible in TORT-TD to model off-centred circular shapes in cylindrical coordinates, the circular control rods are approximated by azimuthal sector portions that preserve the control rod cross sectional areas. The total area covered by the control rod material is then significantly smaller compared to the grey curtain model. If the absorber cross sections remained unchanged, the effectiveness of the control rods would be by far too small and the corresponding eigenvalue (keff = 1.014) too large. This can roughly be compensated by scaling the grey curtain absorption cross sections by the ratio of the grey curtain area to the total area of the azimuthally resolved control rods. As expected, the resulting eigenvalue keff = 0.9169 is roughly comparable with the grey curtain model. Of course, an appropriate cross section set is necessary for an accurate full 3D neutronics description of the PBMR-268 with spatially resolved control rods. This has been done successfully with the MICROX-2 code. The resulting eigenvalue keff = 0.9256 agrees very well with the grey curtain model. 5.2. 3-D coupled TORT-TD/ATTICA3D steady-state analysis of the PBMR-400 benchmark The TORT-TD model for the PBMR-400 is derived from the PBMR-400 benchmark specification. In contrast to the specification, the PBMR-400 is represented in TORT-TD with a full 3-D model, i.e. with the possibility to spatially resolve the 24 control rods (RCS) in the side reflector. For verifying the PBMR-400 neutronics modeling in TORT-TD, stand-alone calculations applying the simplified cross section set (exercise 1) have been performed. The results, for which the radial and axial power profiles are shown in Figs. 7 and 8, have been obtained with level-symmetric S4 quadrature and compare very well with the solutions reported by the other benchmark participants. For the coupled steady-state and transient simulations, the parameterized nuclear cross section library provided with the benchmark specification has been used. It is a two energy group library that is parameterized with respect to the five parameters: fuel and moderator temperature, fast and thermal buckling and Xenon concentration. A software tool was written to calculate the missing self-scatter cross sections from the corresponding total and off-diagonal scattering cross sections and to transform the whole library into the format that can be processed by TORT-TD. In addition to the 5-D linear interpolation routine provided with the benchmark specification, the TORT-TD built-in interpolation with cubic spline polynomials was also applied. Again, a level-symmetric S4 quadrature has been used in TORT-TD.

177

Fig. 7. Axially integrated radial power density profile obtained with TORT-TD (red line) for the simplified cross section set compared to other benchmark participants’ solutions. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 8. Radially integrated axial power density profile obtained with TORT-TD (red line) for the simplified cross section set compared to other benchmark participants’ solutions. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

The ATTICA3D model is based on the simplified flow path depicted in Fig. 9 according to the benchmark specification. The main design parameters used are listed in Table 1. The thermal hydraulic model of the PBMR-400 as proposed by the benchmark description is depicted in Fig. 9. Apparently, a simplified model is applied. Helium enters at the inlet in the lower part of the side reflector with a mass flow rate of 192.7 kg/s, a temperature of 500 ◦ C and a pressure of 90 bar. It then flows upwards through the helium riser channels to the upper plenum which is located sideways of the core. When helium reaches the lower part of the reflector it flows through the flat porous bottom reflector to be gathered in the hot plenum. The hot plenum is also modeled as the gas outlet where gas is removed at a temperature of around 900 ◦ C. For the gas flow no by-pass in gaps or for rod cooling was assumed.

Table 1 Main design parameters used for the ATTICA3D model of the PBMR-400. Description

Unit

Thermal power Inlet temperature Outlet temperature Gas mass flow rate System pressure

MW ◦ C ◦ C kg/s bar

400 500 900 192.7 90

178

A. Seubert et al. / Nuclear Engineering and Design 251 (2012) 173–180

Radial power density profile

W/cm3 5,7 5,5 5,3 5,1 4,9 4,7 4,5 4,3 100

110

120

130

140

150

160

170

180

190

200

Radial distance from center (cm) DALTON-THERMIX PARCS/THERMIX-DIREKT DORT-TD/THERMIX (Transport)

PEBBED-FD/THERMIX-KONVEK TINTE PEBBED-ND/THERMIX-KONVEK

MARS-GCR/CAPP DORT-TD/THERMIX (Diffusion) TORT-TD/ATTICA3D het

Fig. 11. Axially integrated radial power density profile obtained with TORTTD/ATTICA3D (black line) for the steady state compared to other benchmark participants’ solutions.

Average axial fuel temperature profile

°C 950 900 850 800

Fig. 9. Simplified flow path and boundary conditions for the PBMR-400 according to benchmark description.

750 700

DALTON - THERMIX PEBBEDfd-THERMIX-KONVEK(94) MARS-GCR/CAPP PARCS / THERMIX-DIREKT TINTE PEBBEDnd-THERMIX-KONVEK(94) TORT-TD/ATTICA3D

650

The core is represented by an annulus that neglects cones of pebbles on top and bottom. The boundaries are modeled adiabatically at the top and the bottom, and are modeled isothermally (20 ◦ C) at the reactor cavity cooling system. So, heat can only be transferred out of the system via the side of the internal structure materials. This simplified model is subsequently applied to the coupled version of TORT-TD/ATTICA3D. Convergence of the coupled TORT-TD/ATTICA3D steady-state (exercise 3) calculation to less than 10−6 with respect to the maximum local relative change in power distribution between subsequent iterations was obtained after six iterations. The resulting eigenvalue of keff = 0.9862 is compatible with the values reported by the participants of the benchmark and differs by only 87 pcm from the 2-D DORT-TD/THERMIX solution (Tyobeka et al., 2008). Comparisons with the other benchmark participants’ steady state solutions have been made regarding both neutronics and thermalhydraulics quantities. The obtained axial and radial power density profiles are shown in Figs. 10 and 11, respectively. Good agreement

600 550 500 0

100

300

400

500

600

700

800

900

1000 1100 1200

Axial distance from top (cm) Fig. 12. Radially integrated axial fuel temperature profile obtained with TORTTD/ATTICA3D (black line) for the steady state compared to other benchmark participants’ solutions.

has also been achieved for the fuel and moderator temperature profiles as depicted in Figs. 12–15. 5.3. 3-D coupled TORT-TD/ATTICA3D simulation of the total control rod withdrawal transient of the PBMR-400 benchmark As first transient application of TORT-TD/ATTICA3D, the total control rod withdrawal (TCRW) transient specified in exercise 5a

Average axial moderator temperature profile

°C W/cm3 12

200

920

Axial power density profile

870 10

820

8

770

6

720 DALTON - THERMIX

4

PEBBEDfd-THERMIX-KONVEK(94)

670

MARS-GCR/CAPP

2

620

PARCS / THERMIX-DIREKT TINTE

0 0

200 DALTON-THERMIX PARCS/THERMIX-DIREKT DORT-TD/THERMIX (Transport)

400 600 800 Axial distance from top (cm) PEBBED-FD/THERMIX-KONVEK TINTE PEBBED-ND/THERMIX-KONVEK

1000

1200

570

PEBBEDnd-THERMIX-KONVEK(94) TORT-TD/ATTICA3D

520 MARS-GCR/CAPP DORT-TD/THERMIX (Diffusion) TORT-TD/ATTICA3D

Fig. 10. Radially integrated axial power density profile obtained with TORTTD/ATTICA3D (black line) for the steady state compared to other benchmark participants’ solutions.

0

100

200

300

400

500

600

700

800

900

1000

1100

1200

Axial distance from top (cm) Fig. 13. Radially integrated axial moderator temperature profile obtained with TORT-TD/ATTICA3D (black line) for the steady state compared to other benchmark participants’ solutions.

A. Seubert et al. / Nuclear Engineering and Design 251 (2012) 173–180

Total fission power

Average radial fuel temperature profile

°C 900

850

DALTON - THERMIX

800

PEBBEDfd-THERMIX-KONVEK(94)

880

179

MARS-GCR/CAPP

750

PARCS / THERMIX-DIREKT

860

700

TINTE

840

MW

PEBBEDnd-THERMIX-KONVEK(94) TORT-TD/ATTICA3D

820

650 600

KAERI

550

PBMR PURDUE

800

500

PSU-DORT TD Trans PSU-DORT TD Diffusion

450

780

TORT-TD/ATTICA3D

400

760 100

110

120

130

140

150

160

170

180

0

190

40

80

120

160

200

240

280

320

360

400

Time (sec)

Radial distance from center (cm) Fig. 14. Axially integrated radial fuel temperature profile obtained with TORTTD/ATTICA3D (black line) for the steady state compared to other benchmark participants’ solutions.

Fig. 16. Total fission power evolution for the TCRW transient obtained by TORTTD/ATTICA3D (black line) compared to other benchmark participants’ solutions.

Axial power offset has been considered. Starting with operating conditions defined by full power at equilibrium Xenon distribution with the control rods inserted 2.0 m below the bottom of the top reflector, i.e. 1.5 m alongside the pebble-bed, all control rods are completely withdrawn at constant velocity within 200 s. In ATTICA3D, the fuel temperature (Doppler) feedback was calculated by means of a detailed model where the Doppler temperatures are evaluated by explicitly considering the temperature distribution of a representative fuel kernel per individual spherical shell of the fuel element at respective location. The solution of the heat conduction equation is based on a quasi stationary approach, i.e. the steady state heat conduction equation is solved for the local power density distribution present at the current instance in time. The time step size was set to 0.1 s. Before the actual transient, a socalled zero transient starting from the steady state and lasting 2 s was performed in order to verify both codes’ stability when they change to their respective transient mode of operation. Large reactor periods indicate that the steady state calculation was well converged. The total fission power evolutions obtained by TORTTD/ATTICA3D during the TCRW transient is shown in Fig. 16. Obviously, the power evolutions are in reasonable agreement with other benchmark participants’ codes and solutions. The wellknown control rod cusping effect is also present in the TORT-TD

°C

0,85 0,80 0,75 0,70 0,65 0,60 KAERI

0,55

PBMR

0,50

PURDUE

0,45

PSU-DORT TD Trans

0,40

PSU-DORT TD Diffusion TORT-TD/ATTICA3D

0,35 0

40

80

120

280

320

360

400

solution, though to a smaller extent, i.e. with shorter periods due to the finer axial discretisation. The axial power offset has also been evaluated and is compared to the other benchmark participants’ solutions in Fig. 17. Regarding the temporal evolutions of maximum temperatures, they are also found to be within the range of other benchmark

Maximum fuel temperature 1350

DALTON - THERMIX

850

PEBBEDfd-THERMIX-KONVEK(94)

840

MARS-GCR/CAPP

1300 1250

PARCS / THERMIX-DIREKT TINTE

1200

PEBBEDnd-THERMIX-KONVEK(94) TORT-TD/ATTICA3D

°C

810

240

Fig. 17. Axial power offset evolution for the TCRW transient obtained by TORTTD/ATTICA3D (black line) compared to other benchmark participants’ solutions.

Average radial moderator temperature profile

820

200

Time (sec)

860

830

160

1150

800

1100

KAERI

790

1050

PBMR

780

1000

770

950

PURDUE

760 100

110

120

130

140

150

160

170

180

190

Radial distance from center (cm) Fig. 15. Axially integrated radial moderator temperature profile obtained with TORT-TD/ATTICA3D (black line) for the steady state compared to other benchmark participants’ solutions.

PSU-DORTTD Trans PSU-DORTTD Diffusion TORT-TD/ATTICA3D

900 0

40

80

120

160

200

240

280

320

360

400

Time (sec) Fig. 18. Maximum fuel temperature evolution for the TCRW transient obtained by TORT-TD/ATTICA3D (black line) compared to other benchmark participants’ solutions.

180

A. Seubert et al. / Nuclear Engineering and Design 251 (2012) 173–180

Acknowledgements

Maximum moderator temperature 1350

This work was supported by the German Federal Ministry of Economics and Technology due to a resolution of the German Bundestag. Additional support was provided by EnBW Kernkraft GmbH.

1300 1250

°C

1200 1150

References

1100

KAERI

1050

PBMR PURDUE

1000

PSU-DORT TD Trans PSU-DORT TD Diffusion

950

TORT-TD/ATTICA3D

900 0

40

80

120

160

200

240

280

320

360

400

Time (sec) Fig. 19. Maximum moderator temperature evolution for the TCRW transient obtained by TORT-TD/ATTICA3D (black line) compared to other benchmark participants’ solutions.

participants’ results. This is shown in Fig. 18 for the fuel temperature and in Fig. 19 for the moderator temperature.

6. Summary and outlook This paper describes the first application of the time-dependent 3-D discrete-ordinates neutron transport code TORT-TD to HTGR of pebble bed type and the coupling of TORT-TD with the 3-D thermal-hydraulic code ATTICA3D. For the coupled code system TORT-TD/ATTICA3D, the so-called internal coupling approach has been implemented. Therein, ATTICA3D simulates the whole 3D thermal-hydraulics including the reactor core and TORT-TD treats the neutron kinetics. TORT-TD/ATTICA3D is represented by a single executable. The 3-D distributions of temperatures from ATTICA3D and power density from TORT-TD are efficiently exchanged by direct memory access of array elements. For the data exchange, already existing interface routines in TORT-TD developed for coupling with ATHLET and COBRA-TF have been utilized in combination with an ATTICA3D interpolation routine that transfers 3-D data between different meshes. The PBMR-268 and PBMR-400 benchmarks have been considered for first applications of TORT-TD and TORT-TD/ATTICA3D to pebble bed HTGR. For the PBMR-268 benchmark, a 3-D model in TORT-TD has been prepared with azimuthally resolved control rods in the side reflector. A PBMR-268 grey curtain model is used to verify the cylindrical TORT-TD calculation method by comparing with a corresponding 2-D DORT-TD solution. For first coupled TORTTD/ATTICA3D applications, the steady state and the total control rod withdrawal transient as specified in the PBMR-400 benchmark have been simulated. The results are very promising and demonstrate that the coupled code system TORT-TD/ATTICA3D can be an important component of a future comprehensive 3-D code system for HTGR of pebble bed type.

ANL Benchmark Problem Book, 1985. ANL-Report ANL-7416, Suppl. 2 and 3. Ball, S.J., April 1997. CRP-3 Benchmark Problem Description for GT-MHR Pu Burner Accidents, Supplementary Information on the RCCS Design. Oak Ridge National Laboratory. Ball, S.J., February 1997. CRP-3 Benchmark Problem Description for GT-MHR Pu Burner Accidents. Oak Ridge National Laboratory. Christienne, M., Avramova, M., Périn, Y., Seubert, A., 2010. Coupled TORT-TD/CTF capability for high-fidelity LWR core calculations. Physor 2010, Pittsburgh, USA, May 9–14, 2010. Hossain, K., Buck, M., Ben Said, N., Bernnat, W., Lohnert, G., 2008a. Development of a fast 3D thermal-hydraulic tool for design and safety studies for HTRS. Nucl. Eng. Des. 238, 2976. Hossain, K., Buck, M., Bernnat, W., Lohnert, G., 2008b. TH3D, a three-dimensional thermal hydraulics tool, for design and safety analysis of HTRs. In: Proceedings of HTR2008, Washington DC, USA, September 28–October 1. KTA: Auslegung der Reaktorkerne von gasgekühlten Hochtemperaturreaktoren, Wärmeübergänge im Kugelhaufen. KTA 3102.2, 1983. Kuzavkov, N.G., December 1997. CRP-3 Benchmark Problem Description for GT-MHR Pu Burner Accidents, Supplementary Information on the Benchmark Problem. OKBM. Langenbuch, S., Maurer, W., Werner, W., 1977. Coarse mesh flux expansion method for the analysis of space-time effects in large light water reactor cores. Nucl. Sci. Eng. 63, 437. Lerchl, G., Austregesilo, H., 2003. ATHLET Mod2 Cycle A, User’s Manual, GRS. Mathews, D., 1997. An Improved Version of the MICROX-2 Code. Paul Scherrer Institute, Switzerland, PSI Bericht Nr. 97-11. Pautz, A., Tyobeka, B., Ivanov, K., 2008. Application of time-dependent neutron transport theory to high temperature reactors of pebble bed type. Physor 2008, Interlaken, Switzerland, September 14–19, 2008. Reitsma, F., Strydom, G., Sikik, E., 2005. PBMR-268 neutronics and transient benchmark problem, NRG-PBMR-PSU. Reitsma, F., et al., 2006. The PBMR steady-state and coupled kinetics core thermalhydraulics benchmark test problems. Nucl. Eng. Des. 236, 657. Reitsma, F., et al., 2007. PBMR coupled neutronics/thermal hydraulics transient benchmark the PBMR-400 core design, NEA/NSC/DOC(2007), Draft-V07. http://www.nea.fr/science/wprs/pbmr400. Reitsma, F., Han, J., Ivanov, K., Sartori, E., 2008. The OECD/NEA/NSC PBMR400 MW coupled neutronics thermal hydraulics transient benchmark—steady-state results and status. In: Proceedings of the International Conference of Physics of Reactors “Nuclear Power: A Sustainable Resource”, Interlaken, Switzerland, September 14–19. Rhoades, W.A., Childs, R.L., 1991. The TORT three dimensional discrete ordinates neutron/photon transport code. Nucl. Sci. Eng. 107, 397. Rhoades, W.A., Simpson, D.B., 1991. The TORT three dimensional discrete ordinates neutron/photon transport code (TORT version 3), ORNL/TM-13221. Seubert, A., 2010. Pin cell discontinuity factors in the transient discrete ordinates code TORT-TD. Physor 2010, Pittsburgh, PA, USA, May 9–14, 2010. Seubert, A., Velkov, K., Langenbuch, S., 2008. The time-dependent 3-D discrete ordinates code TORT-TD with thermal-hydraulic feedback by ATHLET models. Physor 2008, Interlaken, Switzerland, September 14–19, 2008. Seubert, A., Pautz, A., Becker, M., Dagan, R., 2009. Time-dependent anisotropic distributed source capability in transient 3-D discrete ordinates code TORT-TD. In: International Conference on Mathematics, Computational Methods & Reactor Physics (M&C 2009), Saratoga Springs, USA, May 3–7. Strydom, G., et al., 2010. The OECD/NEA/NSC PBMR 400 MW coupled neutronics thermal hydraulics transient benchmark: transient results. Physor 2010, Pittsburgh, PA, USA, May 9–14, 2010. Tyobeka, B., Pautz, A., Ivanov, K., 2008. DORT-TD/THERMIX solutions for the OECD/NEA/NSC PBMR400MW coupled neutronics thermal hydraulics transient benchmark. Physor 2008, Interlaken, Switzerland, September 14–19, 2008.