- Email: [email protected]
Development and validation of a coupled neutron diffusion-thermal hydraulic calculation procedure for fast reactor applications
Annals of Nuclear Energy 139 (2020) 107243
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
Development and validation of a coupled neutron diffusion-thermal hydraulic calculation procedure for fast reactor applications Xuebei Zhang a, Qin Zeng b, Hongli Chen a a b
School of Physics, University of Science & Technology of China, Hefei 230027, China School of Electricity, South China University of Technology, Guangzhou 510641, China
a r t i c l e
i n f o
Article history: Received 29 June 2019 Received in revised form 1 December 2019 Accepted 4 December 2019
Keywords: Neutron diffusion and thermal-hydraulics coupling UDF and UDS functions 5 5 PWR assembly M2LFR-1000 hot assembly
a b s t r a c t The neutron diffusion equation is defined based on the User Defined Function (UDF) and the User Defined Scalar (UDS) functions of the FLUENT. The neutron diffusion equation is solved iteratively by using the solver of the FLUENT with the Finite Volume Method (FVM). At the same time, the mass, momentum and energy equations are solved iteratively. At each iteration, the power distribution (flux distribution) obtained by the iteration of the neutron diffusion equation is transferred to the thermal-hydraulics calculation and is used as the heat source term. At the same time, the temperature distribution obtained from the thermal-hydraulics calculation is transferred to the neutron diffusion calculation and the macroscopic cross sections of the materials are corrected to realize the coupling calculation of the neutron diffusion and the thermal-hydraulics under the same solver of the FLUENT without needing to develop the interface program and the computational cost is saved. 2D-TWIGL benchmark problem is calculated by the FLUENT solver to verify the feasibility of this method to solve neutron diffusion equation. Through the modeling and calculation of the 5 5 PWR assembly model, the calculation results are compared with the results of other programs to verify the feasibility of the coupling method and the correctness of data transfer. Then this coupling method is applied to calculate the hot assembly of a modular lead-cooled fast reactor (M2LFR-1000) to verify that the thermal-hydraulics characteristics (the maximum fuel temperature and the maximum cladding outer surface temperature) are all within the corresponding thermal-hydraulics design limits. Ó 2019 Elsevier Ltd. All rights reserved.
1. Introduction Traditionally, the best estimation procedure is generally used in reactor design and reactor safety analysis. With the improvement of computer performance and the development of parallel computing, the high confidence simulation of reactor has been paid more attention in the research of reactor design, scheme optimization and safety analysis. Only by considering the multi-physical feedback in reactor simulation, can high confidence simulation be realized. And the neutron diffusion and thermal-hydraulics coupling calculation is an important part of multi-physics coupling calculation (Kulesza et al., 2014; Cacuci, 2004). The actual operation of the reactor is a process of neutrons and thermal reciprocal feedback. Temperature coefficient (fuel temperature coefficient and moderator temperature coefficient) is an important factor for reactivity control in normal operation of reactor, Xie (2004). To achieve accurate calculation of reactor operation and transient conditions, the
E-mail address: [email protected] (X. Zhang) https://doi.org/10.1016/j.anucene.2019.107243 0306-4549/Ó 2019 Elsevier Ltd. All rights reserved.
effects of fuel temperature, moderator temperature and density on local neutron flux and system reactivity must be considered. Computational Fluid Dynamics (CFD) program FLUENT can realize the fine-mesh simulation, Tashakor et al. (2012) of reactor core and fuel assembly, Liu and Cheng (2009), Jiménez et al. (2010), Zhao et al. (2017), Jiménez et al. (2008) by coupling the mass continuity equation, momentum equation and energy conservation equation. The UDS (User Defined Scalar) in the FLUENT can solve a kind of diffusion equation by using the solver in Fluent. It has been widely used in multiphase flow coupling calculation and flow-field and electric field coupling calculation. Wang et al (2008) used the UDS to add the water diffusion equation in air and solid phase to FLUENT. And the hydrodynamic parameters of heat and mass transfer between two phases were added to the UDF of FLUENT to simulate the complex gas-solid multiphase process of batch fluidized bed drying. Donoso-GarcíaL and HenríquezVargas (2015) used the two-dimensional numerical simulation method to simulate the turbulent state of the adiabatic regenerative porous medium burner coupled with thermoelectric components. The time and volume averaged transport equation and the
2
X. Zhang et al. / Annals of Nuclear Energy 139 (2020) 107243
two order turbulence model were adopted. The FLUENT was used to simulate the burner through its UDF (user-defined function) and UDS (user-defined scalar) interface to obtain additional terms involving turbulence and thermal energy. The flow field and electric field were calculated considering the effects of inlet velocity and composition, thermal conductivity of porous media and thermal insulation materials on the burner. Liu et al. (2014) established a three-dimensional (3D) unsteady mathematical model of alumina ball regenerator, and solved it by commercial computational fluid dynamics (CFD) software FLUENT based on the porous medium hypothesis. The standard k-e turbulence model was combined with standard wall function to simulate gas flow and the momentum equation was modified to consider the effect of porous media on fluid flow. The user defined function (UDF) program was programmed in C language and connected with FLUENT. The user-defined scalar (UDS) transport equation of solid energy conservation was defined. And the heat transfer and thermo-physical properties between gas and solid phases were calculated. Jang and Arastoopour (2014) used ANSYS/FLUENT computational fluid dynamics (CFD) program to simulate the gas-solid two-phase flow pattern, the mixing and drying process of drug particles in three different scales of bubbling fluidized bed dryers. The capacity of water transfer and simulation of drying process was calculated. The user-defined scalar transfer equation (UDS) was added to FLUENT to simulate the flow pattern and heat and mass transfer process of drug drying process based on bubbling fluidized bed. Based on the UDF (User Defined Function) and UDS (User Defined Scalar) functions of FLUENT, Xi’an Jiaotong University developed the TASNAM program (Zhang et al., 2009), which is mainly used for numerical calculation of neutron diffusion in molten salt reactor. Based on the user interface programming function of commercial software CFX, Naval Engineering University added three-dimensional space-time neutron dynamics model, coupled with CFD thermal-hydraulics, and simulated the local threedimensional flow behavior and three-dimensional physical characteristics of PWR under steady state Gui et al. (2010). In this paper, the neutron diffusion equation is defined, the modeling tool (GAMBIT) is used to generate meshes and the solver in FLUENT is used to solve the neutron diffusion equation iteratively based on UDF and UDS functions of FLUENT. The thermal-hydraulics calculation is carried out at the same time. In each iteration, thermal power is transferred to thermalhydraulics calculation by solving neutron diffusion equation. The temperature obtained by thermal-hydraulics calculation is transferred to the neutron diffusion calculation, and the macroscopic cross sections of the materials are modified until the convergence of the iterative calculation of neutron diffusion equation and thermal-hydraulics iteration is achieved. The iterative flow chart of coupling calculation is shown by Fig. 1. In order to verify the correctness of the coupling calculation method and data transfer, the 5x5 PWR assembly model, Klas et al. (2014) is modeled and calculated, and the results are compared with other coupling programs. Then the hot assembly of the M2LFR-1000, Chen et al. (2018) is modeled and calculated. And the neutron flux, temperature distribution and the thermal-hydraulics characteristics (the maximum fuel temperature and the maximum cladding outer surface temperature) on the steady state has good agreement with the results calculated by sub-channel program (KMC-SUB), Li et al. (2017). In this paper, the calculation method and mathematical model are introduced in the Section 2. Section 3 describes the calculation of 2D-TWIGL, Ahmad and Abolhasan (2016), Ge et al. (2015) benchmark problem by the FLUENT solver. And the section 4 describes the coupling calculation of 5x5 PWR assembly. The calculation method is applied to the hot assembly of the M2LFR-1000 in Section 5. Section 6 summarizes the general conclusion.
Fig. 1. The flow chart of coupling calculation.
2. Calculation method and mathematical model 2.1. The UDS module of the FLUENT The UDS module of the FLUENT can define a kind of equation and solve it by the inner solver with the finite volume method. The form is shown in formula (1):
@ q/ @ @/ ðw/ C Þ ¼ Su þ @t @xj @xj
ð1Þ
The definitions of equations in the FLUENT are shown in Table 1.
3
X. Zhang et al. / Annals of Nuclear Energy 139 (2020) 107243
where q is density, Cu is the effective diffusion coefficient for the scalar u, Su is the source term ofuper unit volume. Applying Eq. (5) to each control volume, the discretization equation on each given cell is:
Table 1 Definition of equation in the UDS of FLUENT. Name
Expression
Definition
Corresponding functions in UDS
Unsteady-state term
@ q/ @t
DEFINE_UDS_UNSTEADY
Convection term
@ @xj
Discrete form of unsteady state term Flux(/)
Diffusion term
C @@x/2
w/
Diffusivity (CðTÞ)
2
j
Boundary condition
!
1 @/g ð r ;tÞ @t
vg
P P g 0 –g s;g 0 !g
þ
ND P i¼1
DEFINE_UDS_FLUX
!
Specified Value
@/ ) Flux(@x j
Specified Flux
diffusion
equation
is
shown
of u at the cell face, Af is the surface area vector, which means that ! its direction is normal to the surface and Af is the area of the sur-
in
P! ! r Dg r/g ð r ; tÞ ¼ ð r ; tÞþ
g 0 ¼1
m
P f ;g 0
face, ruf is the gradient of u at the face f and V is the cell volume. The equations given above can be applied to multidimensional, unstructured meshes composed of arbitrary polyhedral in FLUENT. For the steady state neutron diffusion equation (3), applying equation (6) to each control volume, the discretization equation on each given cell is:
P! P P ! ! ! Dg r/g ð r Þ Af ¼ ð ð r Þ þ /g0 ð r Þ
N P
r;g
G P ! /g0 ð r ; tÞ þ vg ð1 bÞ
! /g0 ð r ; tÞ
ð2Þ
v
þvg k g
eff
The first term on the left side of the equation (2) is unsteady term, the second term on the left side is diffusion term, and the term on the right side of the equation is the source term. For steady state calculation, the unsteady term in the equation is neglected. And the equation is shown in equation (3): g 0 –g s;g 0 !g
r;g
þvg k
eff
g 0 ¼1
m
f ;g 0
G P g 0 ¼1
m
P f ;g 0
ð7Þ
! /g0 ð r ÞÞV; g ¼ 1; :::; G
2.3. Thermal-hydraulics calculation model The energy equation of coolant region is shown in equation (8):
ðqcoolant cPcoolant U rTÞ ¼ bcoolant U rPþ
ð8Þ
r ðkcoolant rTÞ
P! P P ! ! r Dg r/g ð r Þ ¼ ð r Þ þ /g0 ð r Þ P
g 0 –g s;g 0 !g
r;g
f
vg;i ki C i ð! r ; tÞ; g ¼ 1; :::; G
G vg P
ð6Þ
where N is the number of the faces enclosing a cell; uf is the value
DEFINE_DIFFUSIVITY
Value(/)
The transient neutron equation (2):
N N X X ! ! ! @ qu q f u f U f Af ¼ Cu ruf Af þSu V Vþ @t f f
ð3Þ
! /g0 ð r Þ; g ¼ 1; :::; G
For the equation (3), /g ðrÞ represents the neutron flux, unit cm2s1, Drepresents the neutron diffusion coefficient, unit cm, P f represents the macroscopic fission cross section, the unit is cm1, P g 0 !g represents the macroscopic transfer cross section, the P unit is cm1, r represents the macroscopic removal cross section, the unit is cm1, m represents the average number of neutrons emitted per fission, vg represents the fission spectrum. The effective multiplication factor is calculated by equation (4): G R P P ðnÞ keff
V
¼ ½
i¼1
G R P P V
i¼1
m
m
ðnÞ f ;i /i dV
ð4Þ
ðn1Þ ðn1Þ dV=keff f ;i /i
2.2. Numerical method Finite Volume Method (FVM), Ge et al. (2015) is widely used in CFD methodology to discretize governing equations and is adopted by almost all the popular CFD softwares, ANSYS (2013). FLUENT converts a general transport equation to an algebraic equation using a control-volume-based technique which consists of integrating the general transport equation on each discrete control volume. For a general scalaru, the integral form of a transport equation on a control volume V can be illustrated as follows:
Z
V
@ qu dV þ @t
I
!
!
qu U d A ¼
I
!
Cu r u d A þ
Z
V
Su dV
ð5Þ
Fig. 2. 2D-TWIGL 1/4 core model.
Table 2 Section parameter 3 of 2D-TWIGL. Region
Group
Dg (cm1)
Ra,g (cm1)
tRf,g (cm1)
R1?2 (cm1)
1
1 2 1 2 1 2
1.4 0.4 1.4 0.4 1.3 0.5
0.01 0.15 0.01 0.15 0.008 0.05
0.007 0.2 0.007 0.2 0.003 0.006
0.01
2 3
0.01 0.01
4
X. Zhang et al. / Annals of Nuclear Energy 139 (2020) 107243
(2018) developed by the author’s research group calculates one set of group-wise neutronics parameters including the groupwise macroscopic cross sections, the diffusion coefficients (Dg), the neutron fission yields (mg ) and the fission spectrum (vg ) at 400–1600 K of fuel, air gap, cladding and coolant. The continuous macroscopic cross sections at 400–1600 K of materials can be obtained through interpolation, Bi (2008). Then, the macroscopic cross-section distribution functions of temperature simply shown by equation (13) are added to FLUENT solver through UDF, and the macroscopic cross-section of each grid is updated by reading the new temperature distribution of each grid after each iteration.
X
¼ f ðTÞ
ð13Þ
3. Validation of the FLUENT solver for neutron diffusion Fig. 3. Meshes of the 2D-TWIGL by the solver of FLUENT.
U is the value of the coolant velocity, bcoolant is the temperature dependence of the enthalpy The equation of heat conduction in fuel area is shown in equation (9):
r kfuel rT ¼ q qðrÞ ¼ c
G X X g¼1
ðrÞ/ðrÞ
ð9Þ ð10Þ
In this section, calculation results for 2D-TWIGL benchmark problem are presented and compared with reference values to prove that the FLUENT solver is feasible to solve the neutron diffusion based on the UDS and UDF function. The meshes adopted in this paper are generated by the general mesh generation tool Gambit. The mesh independent solutions are obtained but not presented here.
3.1. 2D-TWIGL seed blanket problem
f ;g
c is the release energy of each fission The heat conduction equation of the cladding and gap is shown in equation (11), (12) respectively:
r kcladding rT ¼ 0
ð11Þ
r kgas rT ¼ 0
ð12Þ
2.4. Material macroscopic cross section library Under the condition of knowing nuclear density and considering the energy group emerging and the influence of temperature on material density, the nuclear database program, Zhou et al.
The 2D-TWIGL benchmark problem, Ahmad and Abolhasan (2016), Ge et al. (2015) is a simplified neutron kinetics model with two neutron energy groups and one delayed neutron precursor family. A steady state calculation and two transient calculations are included in this problem. The reactor core is a 160 cm square consisting of three regions including (1) perturbed seed region containing primary fissile materials with time-dependent properties in the transient situation; (2) unperturbed seed regions containing primary fissile materials with constant properties in the transient situation; (3) a blanket region also containing fissile materials surrounding the whole core. Due to the symmetry, one quadrant of the reactor is modeled for the calculation, as displayed in Fig. 2. The group constants are given in Table 2.
Fig.4. Neutron flux distributions along the diagonal line.
X. Zhang et al. / Annals of Nuclear Energy 139 (2020) 107243
5
Fig. 5. Normalized power diagram for 2D-TWIGL assemblies.
The maximum mesh size is 0.1 cm in x & y directions and meshes are produced by GAMBIT. The number of mesh is about 1280000 and a part of the meshes is shown by Fig. 3.
3.2. Result analysis The results obtained by the FLUENT solver are in good agreement with the reference values. The effective multiplication factor calculated under steady state is 0.913306, which is very close to the reference value 0.913214. The relative error is 10pcm. Fig. 4 shows the fast and thermal neutron flux distributions on the diagonal line of the calculation domain respectively. Fig. 5 shows the comparison between the normalized power and reference value of the calculated assemblies (R represents the reference value, C represents the calculated value by FLUENT solver, e represents the relative error). We can find that the errors mainly occur on the boundary and the interface of different materials. In theory,
there are strong material heterogeneities on the interface areas. And the change of neutron flux on the boundary and interface areas is severe. When the neutron diffusion equation is discretised, the integration of diffusion term is carried out based on the linear assumption of spatial neutron flux. Therefore, the meshes on the boundary and interface areas should be refined to satisfy the linear assumption and reduce the error. So further work is needed to refine the mesh on the boundary and interface areas.
4. Validation of the coupling method 4.1. Model introduction In this paper, the neutron diffusion and thermal-hydraulics coupling calculation method based on UDF and UDS functions of the FLUENT is used to calculate the 5x5 PWR assembly model. The thermo-physical properties of the materials are given by Table 3,
6
X. Zhang et al. / Annals of Nuclear Energy 139 (2020) 107243
the model structure and parameters are given by Fig. 6 and Table 4, and the change of coolant density and specific heat with temperature is given by equation (14) and equation (15). 4.2. Modeling and mesh generation The gambit is used to model and mesh the fuel assembly model of 5x5 PWR. The radial mesh of fuel rod is shown in Fig. 7. The radial mesh of fuel assembly is shown in Fig. 8. The same meshes are used for neutron diffusion and thermal-hydraulics calculation. The axial mesh size is 0.1 m and the total mesh number of the assembly model is 3.10E6. The mesh quality checking tool in the Gambit is used to check the assembly meshes. As shown in Fig. 9, there are 99.32% of the meshes whose EquiSize Skew ranges from 0 to 0.4.
qðTÞ ¼
cP ðTÞ ¼
8 > < > : 8 > < > :
Table 4 Size and material parameters of the assembly. Fuel pin radius Cladding inner radius Cladding outer radius Pitch Fuel height Bottom reflector height Top reflector height Fuel Coolant Gap Cladding Power Energy group Energy boundary
1006; T < 277:6 4:12 108 T 4 þ 6:74 105 T 3 4:31 T 2 T 2 þ 11:8T 154; 277:6 < T < 600
4.1 mm 4.2 mm 4.8 mm 12.5 mm 3m 0.2 m 0.2 m UOX(2% ,4% enrichment) Water, 1000 ppm boron Helium, 0.1 MPa Zircaloy-2 12.5 MW 2 0.625 eV
kg=m3
ð14Þ
657:4; T > 600 4570; T < 500 5:0 105 T 4 0:11T 3 þ 85:2T 2 3:0 104 T þ 4:03 106 ; 500 < T < 600 ðJ=kg:K Þ
ð15Þ
6778; T > 600
4.3. Boundary conditions
4.4. Calculation results and analysis
The coupled calculation boundary conditions are shown in Table 5.
The reference value of effective multiplication factor of the module is 1.17109, Klas et al. (2014) and the effective multiplication factor calculated in this paper is 1.17100. Fig. 10 shows the axial power density distribution of fuel rod center with 2% and 4% enrichment when the neutron diffusion and thermal- hydraulics calculations converge. The blue and black points are the results of other coupling programs. The red and green lines are the results of this paper. The maximum power density of the 4% enrichment fuel rod center is about 4.25E8W/m3, and the maximum power density of the 2% enrichment fuel rod center is about 2.50E8W/ m3 in this paper. And the reference maximum power density of the 4% enrichment fuel rod center is about 4.178E8W/m3, and the reference maximum power density of the 2% enrichment fuel rod center is about 2.45E8W/m3.
Table 3 Thermal properties of the assembly materials. Materials
Density (g/cm3)
Thermal conductivity (W/m.K)
Specific heat capacity (J/kg.K)
Fuel (UO2) Cladding (Zircaloy-2) Coolant (Water) Gap (Helium)
10.3 6.5
3.0 11.0
310 330
0.0001625
0.53 0.152
5193
Viscosity (Pa.s)
0.00009177
Fig. 6. Radial and axial Geometric structure of the assembly.
Fig. 7. Radial mesh of the 5x5 assembly fuel rod.
7
X. Zhang et al. / Annals of Nuclear Energy 139 (2020) 107243 Table 5 Boundary conditions of the coupling calculation of the 5 5 assembly. Field
Boundary
Type
Value
Temperature (T)
Inlet Outlet Side walls Inlet
Constant value Zero gradient Symmetry Extrapolation boundary Extrapolation boundary Symmetry Zero gradient Constant value Constant value Zero gradient
540 K
Neutron flux (/)
Outlet
Pressure (P) Velocity (U)
Fig. 8. Radial mesh of the 5x5 assembly.
It is can be seen from Fig. 10 that the calculation deviation is mainly at the inlet and outlet of the assembly model. The reference value of the power density increases slightly at the inlet and outlet of the assembly, which is mainly due to the influence of the upper and lower reflectors, which make some neutrons be reflected into the fuel area, then cause a slight increase of the power density. In this paper, the influence of the upper and lower reflectors is neglected due to considering the convenience of modeling and meshing by the Gambit. Fig. 11 gives the temperature distribution of the fuel pellet outer surface, Fig. 12 gives the temperature distribution of the fuel cladding inner surface, Fig. 13 gives the temperature distribution of the fuel cladding outer surface and Fig. 14 gives the temperature distribution of the assembly inlet and outlet. Fig. 15 shows the axial coolant temperature distribution of the adjacent fuel rod center and diagonal fuel rod center, Fig. 16 shows the temperature distribution along the X axis direction (z = 0.0 m, y = 0.0 m, z = 0.0 m is symmetry axis). And they are all compared with reference values calculated by other coupling programs and
Side walls Inlet Outlet Inlet Outlet
Gradient on boundary Gradient on boundary
15.5 MPa (0,0,3)m/s
in good agreement with them. The maximum temperature of the 4% enrichment fuel rod center is 1506.97 K, and the maximum temperature of the 2% enrichment fuel rod center is 1066.42 K. And the reference maximum temperature of the 4% enrichment fuel rod center is 1502.22 K, and the reference maximum temperature of the 2% enrichment fuel rod center is 1047.60 K. Through the coupling calculation of the 5x5 PWR assembly, it is proved that this coupling method is feasible and the data transfer is correct. 5. Application of the coupling method The coupling calculation of 5 5 PWR assembly model proves that it is feasible to realize the coupling calculation of neutron diffusion and thermal-hydraulics by utilizing UDF and UDS functions of FLUENT. Now the M2LFR-1000 hot assembly model is calculated by this coupling method. 5.1. Model description M2LFR-1000 is a modular lead-cooled fast reactor. The structure of components and fuel rods is given by Fig. 17 and Fig. 18 respectively. The fuel rods in the core fuel assembly are arranged in a
Fig. 9. mesh quality check.
8
X. Zhang et al. / Annals of Nuclear Energy 139 (2020) 107243
Power density/(W/m3)
500000000
-2
400000000 300000000 200000000 100000000 0 -1.5
-1
-0.5
0
0.5
1
1.5
2
z/m the calculated axial power density distribution of the 2% enrichment fuel rods the calculated axial power density distribution of the 4% enrichment fuel rods the reference axial power density distribution of the 2% enrichment fuel rods the reference axial power density distribution of the 4% enrichment fuel rods Fig. 10. Axial power density distribution the fuel rod center.
Fig. 11. Temperature of the fuel pellet outer surface (unit: K).
Fig. 12. Temperature of fuel cladding inner surface (unit: K).
regular triangular matrix, and the bundles are hexagonal. The bundles are wrapped in the assembly box with a thickness of 4 mm. The distance between the fuel rods is 14 mm. Each fuel assembly contains 169 fuel rods. The inner margin of the fuel assembly box is 185 mm, the outer margin is 193 mm, and the component center distance is 198 mm. The core of the fuel pellet has a central hole with a diameter of 1.9 mm, which can reduce the core temperature and improve the core safety margin under the condition of the same linear power density of the fuel rod. The outer diameter of the fuel pellets is 8.6 mm, and the MOX fuel is added with a small amount of MA or without MA. There is a gap of 0.15 mm between fuel pellet and cladding, which is filled with He of ~ 0.5 MPa. The gas pressure is higher than the operating pressure of the primary circulation. This can improve
the gap heat conduction between fuel pellet and cladding and provide inert environment. On the other hand, it can prevent cladding from contacting with pellet due to external pressure and creep collapse. The cladding damage can also be tested. The fuel cladding is FMS T91 with thickness of 0.55 mm, the diameter of the whole fuel rod is 10.0 mm, and the length of the active zone of the fuel rod is 1000 mm. FMS T91 with good comprehensive performance is selected as the core structure and cladding material. The nuclide composition of FMS T91 is shown in Table 6, David (2015). And the coolant is Pb. The temperature of core inlet coolant is set as 673.15 K and the temperature of core outlet coolant is set as 753.15 K. Under all operational conditions including design basis accidents (DBAs), the maximum fuel pellet temperature should be lower than
9
X. Zhang et al. / Annals of Nuclear Energy 139 (2020) 107243
5.2. Modeling and mesh generation Considering the computational cost, the 1/6 assembly is selected to be modeled and calculated. The assembly is the hot assembly and the power is 3.6 MW. The M2LFR-1000 assembly is modeled and meshed by the Gambit. The axial mesh size is 0.05 m, and the total mesh number is 3.40E6. The radial mesh of the fuel rod and the assembly are shown by Fig. 19 and Fig. 20 respectively. The meshes are checked by Gambit’s grid quality checking tool. There are 96.67% of meshes whose EquiSize Skew ranges from 0 to 0.4. 5.3. Boundary conditions and properties of the materials, Popov et al. (2000) The boundary conditions are showed by Table 7. 5.4. Calculation results and analysis
Fig. 13. Temperature of fuel cladding outer surface (unit: K).
Fig. 21 and Fig. 22 show the unnormalized fast neutron flux distribution on the outlet and the unnormalized thermal neutron flux distribution on the outlet respectively. Because of the fission in the fuel region, the fast neutron is mainly in the fuel region and the thermal neutron is mainly in the coolant region. And great changing of the fuel temperature which makes the macroscopic cross sections of fuel region change greatly leads to apparent changing of the fast neutron flux and thermal neutron flux distribution in this region. However, for the coolant region, the operation temperature is 673 K-753 K and the macroscopic cross sections of Pb hardly change. So the fast neutron flux and thermal neutron flux distribution on the coolant change a little. Fig. 23 shows the fuel pellet centerline temperature distribution along Z axis direction and there is the peak temperature when Z = 0.6 m. And the maximum fuel pellet centerline temperature deviation is 21 K which occurs on the outlet. Fig. 24 shows the assembly temperature distribution along Y axis direction when Z = 0.6 m. Fig. 25 shows the coolant temperature distribution along Z axis direction (average value on X-Y plane). Fig. 26 displays the cladding outer surface temperature distribution along Z axis direction. It is obvious that the central hole makes the central fuel pellet temperature distribution flat, decreases the peak temperature of fuel pellet and improve the safety margin. And there is the maximum fuel temperature in the central fuel rod of the hot assembly. The coolant outlet temperature is 755.11 K. The maximum fuel temperature is 1643.41 K and the maximum cladding outer surface temperature is 773.83 K calculated by coupling calculation. And thermal conductivity of the MOX (W/m.K)
kMOX ¼ 1:0=ð2:85 0:03 þ 0:035 þ T ð0:286 0:715 0:03ÞÞ þ ð6400=T 2:5 Þ eð16:35=TÞ Fig. 14. Temperature of assembly inlet and outlet (unit: K).
densityMOX
Density of the MOX (kg/m )
8 < ð10970 þ 490 0:05Þ ð9:9734 101 þ T ð9:802 106 þ T ð2:705 1010 þ T 4:391 1013 ÞÞÞ3 ; T < 923 ¼ 3 : ð10970 þ 490 0:05Þ ð9:9672 101 þ T ð1:179 105 þ T ð2:429 1009 þ T 1:219 1012 ÞÞÞ ; T > 923
2946.15 K, Carbajo et al. (2001). Under normal conditions, the maximum cladding temperature should be lower than 823.15 K, Chen (2015) with sufficient safety margin.
ð16Þ
3
ð17Þ
Thermal conductivity of the T91 (W/m.K)
kT91 ¼ 32:196 6:534 eðT=546:792Þ
ð18Þ
10
X. Zhang et al. / Annals of Nuclear Energy 139 (2020) 107243
Temperature(K)
570 560 550 540 530
z/m The calculated axial coolant temperature distribution of adjacent fuel rods center The calculated axial coolant temperature distribution of diagonal fuel rods center The reference axial coolant temperature distribution of diagonal fuel rods center The reference axial coolant temperature distribution of adjacent fuel rods center Fig.15. Axial coolant temperature distribution of the contiguous fuel rod center and diagonal fuel rod center.
1600 1400
Temperature/K
1200 1000 800 600 400 200 0 -0.01
0
0.01
0.02
0.03
0.04
0.05
0.06
X/m The reference temperature distribution along the X direction at z=0.0m and y=0.0m The calculated temperature distribution along the X direction at z=0.0m and y=0.0m Fig.16. Temperature distribution in the X axis direction (z = 0.0 m, y = 0.0 m).
densityPb ¼ 11441 1:2795 T
Specific heat capacity of the T91 (J/kg.K)
C PT91 ¼ 404:336 þ T ð0:881 þ T ð2:04 103 þ T 6
2:6755 10 ÞÞ
Specific heat capacity of the Pb (J/kg.K)
ð19Þ ð20Þ 3
Density of the Pb (kg/m )
C PPb ¼ 176:2 4:923 102 T þ 1:544 105 T T 1:524 106 =ðT TÞ
Thermal conductivity of the Pb (W/m.K)
kPb ¼ 9:2 þ 0:011 T
ð21Þ
ð22Þ
Viscosity of the Pb (Pa.s)
v iscosityPb ¼ 4:55 104 e1069=T
ð23Þ
11
X. Zhang et al. / Annals of Nuclear Energy 139 (2020) 107243
Fig. 19. Radial mesh of the M2LFR-1000 fuel rod.
Fig. 17. Fuel assembly cross section (unit: mm).
Fig. 18. Fuel rod cross section.
the maximum fuel temperature is 1647.17 K and the maximum cladding outer surface temperature is 777.28 K calculated by subchannel code (KMC-SUB) which are all within the corresponding thermal-hydraulics design limits. Fig. 20. Radial mesh of the M2LFR-1000 assembly.
6. Conclusion In this study, based on the UDF and UDS function of the FLUENT, the neutron diffusion equation is defined. There is no requirement to develop the interface program of the coupling calculation. Then the assembly is fine modeled and the solver in the FLUENT is used to solve the neutron diffusion equation. The thermal-hydraulics calculation is carried out at the same time. Therefore the coupling
calculation between neutron diffusion and thermal-hydraulics is achieved on the same solver of the FLUENT. In order to achieve the convenient data transfer, the neutron diffusion and thermalhydraulics calculation use the same meshes. The feasibility of the FLUENT solver based on the Finite Volume Method (FVM) to calculate the neutron diffusion is verified
Table 6 The nuclide composition of the FMS T91. Nuclide
Nucleon density (b1cm1)
Nuclide
Nucleon density (b1cm1)
C Cr Ni Mn Mo Si Nb
3.8900E04 7.8690E03 1.5948E04 3.8300E04 4.6270E04 5.8230E04 4.0200E05
N P S Cu V Al Fe
1.6600E04 3.0200E05 7.2845E06 7.3600E05 1.9700E04 6.9300E05 7.4232E02
12
X. Zhang et al. / Annals of Nuclear Energy 139 (2020) 107243 Table 7 Boundary conditions of the coupling calculation of the M2LFR-1000 assembly. Field
Boundary
Type
Temperature (T)
Side walls Inlet Outlet Inlet Side walls Outlet Inlet
Symmetry Constant value Extrapolation boundary Extrapolation boundary Symmetry Constant value Constant value
Neutron flux (/)
Pressure (P) Velocity (U)
Fig. 21. Fast neutron flux on the outlet (unit: n*s/cm*cm).
through calculating the 2D-TWIGL benchmark problem. The correctness of this coupling method is verified through modeling and calculating the 5 5 PWR assembly Then the coupling method is applied to a modular lead-cooled fast reactor (M2LFR-1000). Through modeling and calculating the hot assembly of the M2LFR-1000: (1) The neutron flux and temperature distribution of the hot assembly are obtained on the steady state. The fast neutron is mainly in the fuel region and the thermal neutron is mainly in the coolant region. And the fast neutron flux and thermal neutron flux distribution on the coolant region change a little. But they change a lot on the fuel region.
Value 673 K Gradient on boundary Gradient on boundary 101 KPa (0,0,1.66)m/s
Fig. 22. Thermal neutron flux on the outlet (unit: n*s/cm*cm).
(2) It is obvious that the central hole makes the central fuel pellet temperature distribution flat, decreases the peak temperature of fuel pellet and improves the safety margin at the same power density. (3) The maximum fuel temperature and the maximum cladding outer surface temperature are obtained by the coupling calculation and are compared with the reference values calculated by sub-channel code (KMC-SUB). The error of maximum fuel temperature is 1.25 K and error of maximum cladding outer surface temperature is 2.68 K. The coolant outlet temperature is 755.11 K which is very close to the design value (753.15 K). And these thermal-hydraulics characteristics are all within the corresponding thermalhydraulics design limits.
13
X. Zhang et al. / Annals of Nuclear Energy 139 (2020) 107243
1700
Temperature/K
1600 1500 1400 1300 1200 1100 1000 900 0
0.2
0.4
0.6
0.8
1
Z/m The reference fuel pellet centerline temperature(KMC-SUB) The calculated fuel pellet centerline temperature Fig. 23. The fuel pellet centerline temperature distribution along Z axis direction.
Fig. 24. Temperature distribution along the Y axis direction (Z = 0.6 m, X = 0.0 m).
1.2
14
X. Zhang et al. / Annals of Nuclear Energy 139 (2020) 107243
Fig. 25. The coolant temperature distribution along Z axis direction (average value on X-Y plane).
Fig. 26. The cladding outer surface temperature distribution along Z axis direction.
X. Zhang et al. / Annals of Nuclear Energy 139 (2020) 107243
Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgment This work is supported by National Natural Science Foundation of China (Grant No. 11705195). Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi.org/10.1016/j.anucene.2019.107243. References Kulesza, J.A., Franceschini, F., Evans, T.M., et al., 2014. Overview of the consortium for the advanced simulation of light water reactor (CASL). CASL-U-2014-0099000. Cacuci, D., 2004. European platform for nuclear reactor simulation (NURESIM). Integrated Project NUCTECH-2004-3. 4. 3. 1-1. EURATOM Research and Training Program on Nuclear Energy. Xie, Z.S., 2004. Physical Analysis of Nuclear Reactor. Xi’an: Xi’an Jiao Tong University press. Tashakor, S., Salehi, A.A., Jahanfarnia, G., Abbaspour, A., Fard, Tehrani, 2012. Neutronic analysis of HPLWR fuel assembly cluster. Ann. Nucl. Energy 50, 38– 43. Liu, X.J., Cheng, X., 2009. Thermal-hydraulic and neutron-physical characteristics of a new SCWR fuel assembly. Ann. Nucl. Energy 36, 28–36. Jiménez, J. et al., 2010. A domain decomposition methodology for pin by pin coupled neutronic and thermal–hydraulic analyses in COBAYA3. Nucl. Eng. Des. 240, 313–320. Zhao, C., et al., 2016. Coupled neutronics and thermal–hydraulics analysis of annular fuel assembly for SCWR. In: Jiang H. (Eds.), Proceedings of The 20th Pacific Basin Nuclear Conference. PBNC 2016. Springer, Singapore. Jiménez, J., et al., 2008. Comparative analysis of neutronics/thermal-hydraulics multi-scale coupling for LWR analysis. In: International Conference on the Physics of Reactors ‘‘Nuclear Power: A Sustainable Resource” Casino-Kursaal Conference Center, Interlaken, Switzerland, September 14–19.
15
Wang, H.G., et al., 2008. Investigation of batch fluidized-bed drying by mathematical modeling, CFD simulation and ECT measurement. Wiley J., 54, 427–444 doi: 10.1002/aic.11406. Donoso-GarcíaL, P., Henríquez-Vargas, L., 2015. Numerical study of turbulent porous media combustion coupled with thermoelectric generation in a recuperative reactor. Energy 93, 1189–1198. https://doi.org/10.1016/j. energy.2015.09.123. Liu, Y. et al., 2014. Three-dimensional analysis of gas flow and heat transfer in a regenerator with alumina balls. Appl. Thermal Eng. 69, 113–122. https://doi. org/10.1016/j.applthermaleng.2014.04.058. Jang, J., Arastoopour, H., 2014. CFD simulation of a pharmaceutical bubbling bed drying process at three different scales. Powder Technol. 263, 14–25. https:// doi.org/10.1016/j.powtec.2014.04.054. Zhang, D.L. et al., 2009. Development of a steady state analysis code for a molten salt reactor. Ann. Nucl. Energy 36, 590–603. Gui, X.W. et al., 2010. Study on coupling of local three-dimension flow model based on CFD method and space-time neutron kinetics model. Chin. J. Nucl. Sci. Eng. 3, 216–222. Klas, J. et al., 2014. Fine-mesh deterministic modeling of PWR fuel assemblies: proof-of-principle of coupled neutronic/thermal–hydraulic calculations. Ann. Nucl. Energy 68, 247–256. Chen, H.L. et al., 2018. Preliminary design of a medium-power modular lead-cooled fast reactor with the application of optimization methods. Int. J. Energy Res. 42, 3643–3657. Li, S.Z. et al., 2017. Development of a sub-channel thermal hydraulic analysis code and its application to lead cooled fast reactor. Appl. Therm. Eng. 117, 443–451. Pirouzmand, Ahmad, Nabavi, Abolhasan, 2016. Simulation of nuclear reactor dynamics equations using reconfigurable computing. Prog. Nucl. Energy 89, 197–203. Ge, J. et al., 2015. Steady and transient solutions of neutronics problems based on finite volume method(FVM) with a CFD code. Prog. Nucl. Energy 85, 366–374. ANSYS, 2013. ANSYS FLUENT Theory Guide, Release 15.0. Zhou, X.B. et al., 2018. Development and preliminary test of data library ANDL-ADS for accelerator-driven systems. Nucl. Tech. 41 (03), 65–70. Bi, G.W., 2008. Interpolation Method Development for Temperature Based Neutron Cross-sections. Tsinghua University, Beijing. David Jaluvka, 2015. Development of a Core Management Tool for the MYRRHA Irradiation Research Facility [D]. KU Leuven. Carbajo, J.J., Yoder, G.L., Popov, S.G., Ivanov, V.K., 2001. A review of the thermophysical properties of MOX and UO2 fuels. J. Nucl. Mater. 299 (3), 181–198. Chen, Z., 2015. Thermal-hydraulics design and safety analysis of a 100MWth small natural circulation lead cooled fast reactor SNCLFR-100. University of Science and Technology of China. Popov, S.G., Carbajo, J.J., Ivanov, V.K., Yoder, G.L., 2000. Thermophysical properties of MOX and UO2 fuels including the effects of irradiation. ORNL Report TM-2000/ 351.
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
Development and validation of a coupled neutron diffusion-thermal hydraulic calculation procedure for fast reactor applications Xuebei Zhang a, Qin Zeng b, Hongli Chen a a b
School of Physics, University of Science & Technology of China, Hefei 230027, China School of Electricity, South China University of Technology, Guangzhou 510641, China
a r t i c l e
i n f o
Article history: Received 29 June 2019 Received in revised form 1 December 2019 Accepted 4 December 2019
Keywords: Neutron diffusion and thermal-hydraulics coupling UDF and UDS functions 5 5 PWR assembly M2LFR-1000 hot assembly
a b s t r a c t The neutron diffusion equation is defined based on the User Defined Function (UDF) and the User Defined Scalar (UDS) functions of the FLUENT. The neutron diffusion equation is solved iteratively by using the solver of the FLUENT with the Finite Volume Method (FVM). At the same time, the mass, momentum and energy equations are solved iteratively. At each iteration, the power distribution (flux distribution) obtained by the iteration of the neutron diffusion equation is transferred to the thermal-hydraulics calculation and is used as the heat source term. At the same time, the temperature distribution obtained from the thermal-hydraulics calculation is transferred to the neutron diffusion calculation and the macroscopic cross sections of the materials are corrected to realize the coupling calculation of the neutron diffusion and the thermal-hydraulics under the same solver of the FLUENT without needing to develop the interface program and the computational cost is saved. 2D-TWIGL benchmark problem is calculated by the FLUENT solver to verify the feasibility of this method to solve neutron diffusion equation. Through the modeling and calculation of the 5 5 PWR assembly model, the calculation results are compared with the results of other programs to verify the feasibility of the coupling method and the correctness of data transfer. Then this coupling method is applied to calculate the hot assembly of a modular lead-cooled fast reactor (M2LFR-1000) to verify that the thermal-hydraulics characteristics (the maximum fuel temperature and the maximum cladding outer surface temperature) are all within the corresponding thermal-hydraulics design limits. Ó 2019 Elsevier Ltd. All rights reserved.
1. Introduction Traditionally, the best estimation procedure is generally used in reactor design and reactor safety analysis. With the improvement of computer performance and the development of parallel computing, the high confidence simulation of reactor has been paid more attention in the research of reactor design, scheme optimization and safety analysis. Only by considering the multi-physical feedback in reactor simulation, can high confidence simulation be realized. And the neutron diffusion and thermal-hydraulics coupling calculation is an important part of multi-physics coupling calculation (Kulesza et al., 2014; Cacuci, 2004). The actual operation of the reactor is a process of neutrons and thermal reciprocal feedback. Temperature coefficient (fuel temperature coefficient and moderator temperature coefficient) is an important factor for reactivity control in normal operation of reactor, Xie (2004). To achieve accurate calculation of reactor operation and transient conditions, the
E-mail address: [email protected] (X. Zhang) https://doi.org/10.1016/j.anucene.2019.107243 0306-4549/Ó 2019 Elsevier Ltd. All rights reserved.
effects of fuel temperature, moderator temperature and density on local neutron flux and system reactivity must be considered. Computational Fluid Dynamics (CFD) program FLUENT can realize the fine-mesh simulation, Tashakor et al. (2012) of reactor core and fuel assembly, Liu and Cheng (2009), Jiménez et al. (2010), Zhao et al. (2017), Jiménez et al. (2008) by coupling the mass continuity equation, momentum equation and energy conservation equation. The UDS (User Defined Scalar) in the FLUENT can solve a kind of diffusion equation by using the solver in Fluent. It has been widely used in multiphase flow coupling calculation and flow-field and electric field coupling calculation. Wang et al (2008) used the UDS to add the water diffusion equation in air and solid phase to FLUENT. And the hydrodynamic parameters of heat and mass transfer between two phases were added to the UDF of FLUENT to simulate the complex gas-solid multiphase process of batch fluidized bed drying. Donoso-GarcíaL and HenríquezVargas (2015) used the two-dimensional numerical simulation method to simulate the turbulent state of the adiabatic regenerative porous medium burner coupled with thermoelectric components. The time and volume averaged transport equation and the
2
X. Zhang et al. / Annals of Nuclear Energy 139 (2020) 107243
two order turbulence model were adopted. The FLUENT was used to simulate the burner through its UDF (user-defined function) and UDS (user-defined scalar) interface to obtain additional terms involving turbulence and thermal energy. The flow field and electric field were calculated considering the effects of inlet velocity and composition, thermal conductivity of porous media and thermal insulation materials on the burner. Liu et al. (2014) established a three-dimensional (3D) unsteady mathematical model of alumina ball regenerator, and solved it by commercial computational fluid dynamics (CFD) software FLUENT based on the porous medium hypothesis. The standard k-e turbulence model was combined with standard wall function to simulate gas flow and the momentum equation was modified to consider the effect of porous media on fluid flow. The user defined function (UDF) program was programmed in C language and connected with FLUENT. The user-defined scalar (UDS) transport equation of solid energy conservation was defined. And the heat transfer and thermo-physical properties between gas and solid phases were calculated. Jang and Arastoopour (2014) used ANSYS/FLUENT computational fluid dynamics (CFD) program to simulate the gas-solid two-phase flow pattern, the mixing and drying process of drug particles in three different scales of bubbling fluidized bed dryers. The capacity of water transfer and simulation of drying process was calculated. The user-defined scalar transfer equation (UDS) was added to FLUENT to simulate the flow pattern and heat and mass transfer process of drug drying process based on bubbling fluidized bed. Based on the UDF (User Defined Function) and UDS (User Defined Scalar) functions of FLUENT, Xi’an Jiaotong University developed the TASNAM program (Zhang et al., 2009), which is mainly used for numerical calculation of neutron diffusion in molten salt reactor. Based on the user interface programming function of commercial software CFX, Naval Engineering University added three-dimensional space-time neutron dynamics model, coupled with CFD thermal-hydraulics, and simulated the local threedimensional flow behavior and three-dimensional physical characteristics of PWR under steady state Gui et al. (2010). In this paper, the neutron diffusion equation is defined, the modeling tool (GAMBIT) is used to generate meshes and the solver in FLUENT is used to solve the neutron diffusion equation iteratively based on UDF and UDS functions of FLUENT. The thermal-hydraulics calculation is carried out at the same time. In each iteration, thermal power is transferred to thermalhydraulics calculation by solving neutron diffusion equation. The temperature obtained by thermal-hydraulics calculation is transferred to the neutron diffusion calculation, and the macroscopic cross sections of the materials are modified until the convergence of the iterative calculation of neutron diffusion equation and thermal-hydraulics iteration is achieved. The iterative flow chart of coupling calculation is shown by Fig. 1. In order to verify the correctness of the coupling calculation method and data transfer, the 5x5 PWR assembly model, Klas et al. (2014) is modeled and calculated, and the results are compared with other coupling programs. Then the hot assembly of the M2LFR-1000, Chen et al. (2018) is modeled and calculated. And the neutron flux, temperature distribution and the thermal-hydraulics characteristics (the maximum fuel temperature and the maximum cladding outer surface temperature) on the steady state has good agreement with the results calculated by sub-channel program (KMC-SUB), Li et al. (2017). In this paper, the calculation method and mathematical model are introduced in the Section 2. Section 3 describes the calculation of 2D-TWIGL, Ahmad and Abolhasan (2016), Ge et al. (2015) benchmark problem by the FLUENT solver. And the section 4 describes the coupling calculation of 5x5 PWR assembly. The calculation method is applied to the hot assembly of the M2LFR-1000 in Section 5. Section 6 summarizes the general conclusion.
Fig. 1. The flow chart of coupling calculation.
2. Calculation method and mathematical model 2.1. The UDS module of the FLUENT The UDS module of the FLUENT can define a kind of equation and solve it by the inner solver with the finite volume method. The form is shown in formula (1):
@ q/ @ @/ ðw/ C Þ ¼ Su þ @t @xj @xj
ð1Þ
The definitions of equations in the FLUENT are shown in Table 1.
3
X. Zhang et al. / Annals of Nuclear Energy 139 (2020) 107243
where q is density, Cu is the effective diffusion coefficient for the scalar u, Su is the source term ofuper unit volume. Applying Eq. (5) to each control volume, the discretization equation on each given cell is:
Table 1 Definition of equation in the UDS of FLUENT. Name
Expression
Definition
Corresponding functions in UDS
Unsteady-state term
@ q/ @t
DEFINE_UDS_UNSTEADY
Convection term
@ @xj
Discrete form of unsteady state term Flux(/)
Diffusion term
C @@x/2
w/
Diffusivity (CðTÞ)
2
j
Boundary condition
!
1 @/g ð r ;tÞ @t
vg
P P g 0 –g s;g 0 !g
þ
ND P i¼1
DEFINE_UDS_FLUX
!
Specified Value
@/ ) Flux(@x j
Specified Flux
diffusion
equation
is
shown
of u at the cell face, Af is the surface area vector, which means that ! its direction is normal to the surface and Af is the area of the sur-
in
P! ! r Dg r/g ð r ; tÞ ¼ ð r ; tÞþ
g 0 ¼1
m
P f ;g 0
face, ruf is the gradient of u at the face f and V is the cell volume. The equations given above can be applied to multidimensional, unstructured meshes composed of arbitrary polyhedral in FLUENT. For the steady state neutron diffusion equation (3), applying equation (6) to each control volume, the discretization equation on each given cell is:
P! P P ! ! ! Dg r/g ð r Þ Af ¼ ð ð r Þ þ /g0 ð r Þ
N P
r;g
G P ! /g0 ð r ; tÞ þ vg ð1 bÞ
! /g0 ð r ; tÞ
ð2Þ
v
þvg k g
eff
The first term on the left side of the equation (2) is unsteady term, the second term on the left side is diffusion term, and the term on the right side of the equation is the source term. For steady state calculation, the unsteady term in the equation is neglected. And the equation is shown in equation (3): g 0 –g s;g 0 !g
r;g
þvg k
eff
g 0 ¼1
m
f ;g 0
G P g 0 ¼1
m
P f ;g 0
ð7Þ
! /g0 ð r ÞÞV; g ¼ 1; :::; G
2.3. Thermal-hydraulics calculation model The energy equation of coolant region is shown in equation (8):
ðqcoolant cPcoolant U rTÞ ¼ bcoolant U rPþ
ð8Þ
r ðkcoolant rTÞ
P! P P ! ! r Dg r/g ð r Þ ¼ ð r Þ þ /g0 ð r Þ P
g 0 –g s;g 0 !g
r;g
f
vg;i ki C i ð! r ; tÞ; g ¼ 1; :::; G
G vg P
ð6Þ
where N is the number of the faces enclosing a cell; uf is the value
DEFINE_DIFFUSIVITY
Value(/)
The transient neutron equation (2):
N N X X ! ! ! @ qu q f u f U f Af ¼ Cu ruf Af þSu V Vþ @t f f
ð3Þ
! /g0 ð r Þ; g ¼ 1; :::; G
For the equation (3), /g ðrÞ represents the neutron flux, unit cm2s1, Drepresents the neutron diffusion coefficient, unit cm, P f represents the macroscopic fission cross section, the unit is cm1, P g 0 !g represents the macroscopic transfer cross section, the P unit is cm1, r represents the macroscopic removal cross section, the unit is cm1, m represents the average number of neutrons emitted per fission, vg represents the fission spectrum. The effective multiplication factor is calculated by equation (4): G R P P ðnÞ keff
V
¼ ½
i¼1
G R P P V
i¼1
m
m
ðnÞ f ;i /i dV
ð4Þ
ðn1Þ ðn1Þ dV=keff f ;i /i
2.2. Numerical method Finite Volume Method (FVM), Ge et al. (2015) is widely used in CFD methodology to discretize governing equations and is adopted by almost all the popular CFD softwares, ANSYS (2013). FLUENT converts a general transport equation to an algebraic equation using a control-volume-based technique which consists of integrating the general transport equation on each discrete control volume. For a general scalaru, the integral form of a transport equation on a control volume V can be illustrated as follows:
Z
V
@ qu dV þ @t
I
!
!
qu U d A ¼
I
!
Cu r u d A þ
Z
V
Su dV
ð5Þ
Fig. 2. 2D-TWIGL 1/4 core model.
Table 2 Section parameter 3 of 2D-TWIGL. Region
Group
Dg (cm1)
Ra,g (cm1)
tRf,g (cm1)
R1?2 (cm1)
1
1 2 1 2 1 2
1.4 0.4 1.4 0.4 1.3 0.5
0.01 0.15 0.01 0.15 0.008 0.05
0.007 0.2 0.007 0.2 0.003 0.006
0.01
2 3
0.01 0.01
4
X. Zhang et al. / Annals of Nuclear Energy 139 (2020) 107243
(2018) developed by the author’s research group calculates one set of group-wise neutronics parameters including the groupwise macroscopic cross sections, the diffusion coefficients (Dg), the neutron fission yields (mg ) and the fission spectrum (vg ) at 400–1600 K of fuel, air gap, cladding and coolant. The continuous macroscopic cross sections at 400–1600 K of materials can be obtained through interpolation, Bi (2008). Then, the macroscopic cross-section distribution functions of temperature simply shown by equation (13) are added to FLUENT solver through UDF, and the macroscopic cross-section of each grid is updated by reading the new temperature distribution of each grid after each iteration.
X
¼ f ðTÞ
ð13Þ
3. Validation of the FLUENT solver for neutron diffusion Fig. 3. Meshes of the 2D-TWIGL by the solver of FLUENT.
U is the value of the coolant velocity, bcoolant is the temperature dependence of the enthalpy The equation of heat conduction in fuel area is shown in equation (9):
r kfuel rT ¼ q qðrÞ ¼ c
G X X g¼1
ðrÞ/ðrÞ
ð9Þ ð10Þ
In this section, calculation results for 2D-TWIGL benchmark problem are presented and compared with reference values to prove that the FLUENT solver is feasible to solve the neutron diffusion based on the UDS and UDF function. The meshes adopted in this paper are generated by the general mesh generation tool Gambit. The mesh independent solutions are obtained but not presented here.
3.1. 2D-TWIGL seed blanket problem
f ;g
c is the release energy of each fission The heat conduction equation of the cladding and gap is shown in equation (11), (12) respectively:
r kcladding rT ¼ 0
ð11Þ
r kgas rT ¼ 0
ð12Þ
2.4. Material macroscopic cross section library Under the condition of knowing nuclear density and considering the energy group emerging and the influence of temperature on material density, the nuclear database program, Zhou et al.
The 2D-TWIGL benchmark problem, Ahmad and Abolhasan (2016), Ge et al. (2015) is a simplified neutron kinetics model with two neutron energy groups and one delayed neutron precursor family. A steady state calculation and two transient calculations are included in this problem. The reactor core is a 160 cm square consisting of three regions including (1) perturbed seed region containing primary fissile materials with time-dependent properties in the transient situation; (2) unperturbed seed regions containing primary fissile materials with constant properties in the transient situation; (3) a blanket region also containing fissile materials surrounding the whole core. Due to the symmetry, one quadrant of the reactor is modeled for the calculation, as displayed in Fig. 2. The group constants are given in Table 2.
Fig.4. Neutron flux distributions along the diagonal line.
X. Zhang et al. / Annals of Nuclear Energy 139 (2020) 107243
5
Fig. 5. Normalized power diagram for 2D-TWIGL assemblies.
The maximum mesh size is 0.1 cm in x & y directions and meshes are produced by GAMBIT. The number of mesh is about 1280000 and a part of the meshes is shown by Fig. 3.
3.2. Result analysis The results obtained by the FLUENT solver are in good agreement with the reference values. The effective multiplication factor calculated under steady state is 0.913306, which is very close to the reference value 0.913214. The relative error is 10pcm. Fig. 4 shows the fast and thermal neutron flux distributions on the diagonal line of the calculation domain respectively. Fig. 5 shows the comparison between the normalized power and reference value of the calculated assemblies (R represents the reference value, C represents the calculated value by FLUENT solver, e represents the relative error). We can find that the errors mainly occur on the boundary and the interface of different materials. In theory,
there are strong material heterogeneities on the interface areas. And the change of neutron flux on the boundary and interface areas is severe. When the neutron diffusion equation is discretised, the integration of diffusion term is carried out based on the linear assumption of spatial neutron flux. Therefore, the meshes on the boundary and interface areas should be refined to satisfy the linear assumption and reduce the error. So further work is needed to refine the mesh on the boundary and interface areas.
4. Validation of the coupling method 4.1. Model introduction In this paper, the neutron diffusion and thermal-hydraulics coupling calculation method based on UDF and UDS functions of the FLUENT is used to calculate the 5x5 PWR assembly model. The thermo-physical properties of the materials are given by Table 3,
6
X. Zhang et al. / Annals of Nuclear Energy 139 (2020) 107243
the model structure and parameters are given by Fig. 6 and Table 4, and the change of coolant density and specific heat with temperature is given by equation (14) and equation (15). 4.2. Modeling and mesh generation The gambit is used to model and mesh the fuel assembly model of 5x5 PWR. The radial mesh of fuel rod is shown in Fig. 7. The radial mesh of fuel assembly is shown in Fig. 8. The same meshes are used for neutron diffusion and thermal-hydraulics calculation. The axial mesh size is 0.1 m and the total mesh number of the assembly model is 3.10E6. The mesh quality checking tool in the Gambit is used to check the assembly meshes. As shown in Fig. 9, there are 99.32% of the meshes whose EquiSize Skew ranges from 0 to 0.4.
qðTÞ ¼
cP ðTÞ ¼
8 > < > : 8 > < > :
Table 4 Size and material parameters of the assembly. Fuel pin radius Cladding inner radius Cladding outer radius Pitch Fuel height Bottom reflector height Top reflector height Fuel Coolant Gap Cladding Power Energy group Energy boundary
1006; T < 277:6 4:12 108 T 4 þ 6:74 105 T 3 4:31 T 2 T 2 þ 11:8T 154; 277:6 < T < 600
4.1 mm 4.2 mm 4.8 mm 12.5 mm 3m 0.2 m 0.2 m UOX(2% ,4% enrichment) Water, 1000 ppm boron Helium, 0.1 MPa Zircaloy-2 12.5 MW 2 0.625 eV
kg=m3
ð14Þ
657:4; T > 600 4570; T < 500 5:0 105 T 4 0:11T 3 þ 85:2T 2 3:0 104 T þ 4:03 106 ; 500 < T < 600 ðJ=kg:K Þ
ð15Þ
6778; T > 600
4.3. Boundary conditions
4.4. Calculation results and analysis
The coupled calculation boundary conditions are shown in Table 5.
The reference value of effective multiplication factor of the module is 1.17109, Klas et al. (2014) and the effective multiplication factor calculated in this paper is 1.17100. Fig. 10 shows the axial power density distribution of fuel rod center with 2% and 4% enrichment when the neutron diffusion and thermal- hydraulics calculations converge. The blue and black points are the results of other coupling programs. The red and green lines are the results of this paper. The maximum power density of the 4% enrichment fuel rod center is about 4.25E8W/m3, and the maximum power density of the 2% enrichment fuel rod center is about 2.50E8W/ m3 in this paper. And the reference maximum power density of the 4% enrichment fuel rod center is about 4.178E8W/m3, and the reference maximum power density of the 2% enrichment fuel rod center is about 2.45E8W/m3.
Table 3 Thermal properties of the assembly materials. Materials
Density (g/cm3)
Thermal conductivity (W/m.K)
Specific heat capacity (J/kg.K)
Fuel (UO2) Cladding (Zircaloy-2) Coolant (Water) Gap (Helium)
10.3 6.5
3.0 11.0
310 330
0.0001625
0.53 0.152
5193
Viscosity (Pa.s)
0.00009177
Fig. 6. Radial and axial Geometric structure of the assembly.
Fig. 7. Radial mesh of the 5x5 assembly fuel rod.
7
X. Zhang et al. / Annals of Nuclear Energy 139 (2020) 107243 Table 5 Boundary conditions of the coupling calculation of the 5 5 assembly. Field
Boundary
Type
Value
Temperature (T)
Inlet Outlet Side walls Inlet
Constant value Zero gradient Symmetry Extrapolation boundary Extrapolation boundary Symmetry Zero gradient Constant value Constant value Zero gradient
540 K
Neutron flux (/)
Outlet
Pressure (P) Velocity (U)
Fig. 8. Radial mesh of the 5x5 assembly.
It is can be seen from Fig. 10 that the calculation deviation is mainly at the inlet and outlet of the assembly model. The reference value of the power density increases slightly at the inlet and outlet of the assembly, which is mainly due to the influence of the upper and lower reflectors, which make some neutrons be reflected into the fuel area, then cause a slight increase of the power density. In this paper, the influence of the upper and lower reflectors is neglected due to considering the convenience of modeling and meshing by the Gambit. Fig. 11 gives the temperature distribution of the fuel pellet outer surface, Fig. 12 gives the temperature distribution of the fuel cladding inner surface, Fig. 13 gives the temperature distribution of the fuel cladding outer surface and Fig. 14 gives the temperature distribution of the assembly inlet and outlet. Fig. 15 shows the axial coolant temperature distribution of the adjacent fuel rod center and diagonal fuel rod center, Fig. 16 shows the temperature distribution along the X axis direction (z = 0.0 m, y = 0.0 m, z = 0.0 m is symmetry axis). And they are all compared with reference values calculated by other coupling programs and
Side walls Inlet Outlet Inlet Outlet
Gradient on boundary Gradient on boundary
15.5 MPa (0,0,3)m/s
in good agreement with them. The maximum temperature of the 4% enrichment fuel rod center is 1506.97 K, and the maximum temperature of the 2% enrichment fuel rod center is 1066.42 K. And the reference maximum temperature of the 4% enrichment fuel rod center is 1502.22 K, and the reference maximum temperature of the 2% enrichment fuel rod center is 1047.60 K. Through the coupling calculation of the 5x5 PWR assembly, it is proved that this coupling method is feasible and the data transfer is correct. 5. Application of the coupling method The coupling calculation of 5 5 PWR assembly model proves that it is feasible to realize the coupling calculation of neutron diffusion and thermal-hydraulics by utilizing UDF and UDS functions of FLUENT. Now the M2LFR-1000 hot assembly model is calculated by this coupling method. 5.1. Model description M2LFR-1000 is a modular lead-cooled fast reactor. The structure of components and fuel rods is given by Fig. 17 and Fig. 18 respectively. The fuel rods in the core fuel assembly are arranged in a
Fig. 9. mesh quality check.
8
X. Zhang et al. / Annals of Nuclear Energy 139 (2020) 107243
Power density/(W/m3)
500000000
-2
400000000 300000000 200000000 100000000 0 -1.5
-1
-0.5
0
0.5
1
1.5
2
z/m the calculated axial power density distribution of the 2% enrichment fuel rods the calculated axial power density distribution of the 4% enrichment fuel rods the reference axial power density distribution of the 2% enrichment fuel rods the reference axial power density distribution of the 4% enrichment fuel rods Fig. 10. Axial power density distribution the fuel rod center.
Fig. 11. Temperature of the fuel pellet outer surface (unit: K).
Fig. 12. Temperature of fuel cladding inner surface (unit: K).
regular triangular matrix, and the bundles are hexagonal. The bundles are wrapped in the assembly box with a thickness of 4 mm. The distance between the fuel rods is 14 mm. Each fuel assembly contains 169 fuel rods. The inner margin of the fuel assembly box is 185 mm, the outer margin is 193 mm, and the component center distance is 198 mm. The core of the fuel pellet has a central hole with a diameter of 1.9 mm, which can reduce the core temperature and improve the core safety margin under the condition of the same linear power density of the fuel rod. The outer diameter of the fuel pellets is 8.6 mm, and the MOX fuel is added with a small amount of MA or without MA. There is a gap of 0.15 mm between fuel pellet and cladding, which is filled with He of ~ 0.5 MPa. The gas pressure is higher than the operating pressure of the primary circulation. This can improve
the gap heat conduction between fuel pellet and cladding and provide inert environment. On the other hand, it can prevent cladding from contacting with pellet due to external pressure and creep collapse. The cladding damage can also be tested. The fuel cladding is FMS T91 with thickness of 0.55 mm, the diameter of the whole fuel rod is 10.0 mm, and the length of the active zone of the fuel rod is 1000 mm. FMS T91 with good comprehensive performance is selected as the core structure and cladding material. The nuclide composition of FMS T91 is shown in Table 6, David (2015). And the coolant is Pb. The temperature of core inlet coolant is set as 673.15 K and the temperature of core outlet coolant is set as 753.15 K. Under all operational conditions including design basis accidents (DBAs), the maximum fuel pellet temperature should be lower than
9
X. Zhang et al. / Annals of Nuclear Energy 139 (2020) 107243
5.2. Modeling and mesh generation Considering the computational cost, the 1/6 assembly is selected to be modeled and calculated. The assembly is the hot assembly and the power is 3.6 MW. The M2LFR-1000 assembly is modeled and meshed by the Gambit. The axial mesh size is 0.05 m, and the total mesh number is 3.40E6. The radial mesh of the fuel rod and the assembly are shown by Fig. 19 and Fig. 20 respectively. The meshes are checked by Gambit’s grid quality checking tool. There are 96.67% of meshes whose EquiSize Skew ranges from 0 to 0.4. 5.3. Boundary conditions and properties of the materials, Popov et al. (2000) The boundary conditions are showed by Table 7. 5.4. Calculation results and analysis
Fig. 13. Temperature of fuel cladding outer surface (unit: K).
Fig. 21 and Fig. 22 show the unnormalized fast neutron flux distribution on the outlet and the unnormalized thermal neutron flux distribution on the outlet respectively. Because of the fission in the fuel region, the fast neutron is mainly in the fuel region and the thermal neutron is mainly in the coolant region. And great changing of the fuel temperature which makes the macroscopic cross sections of fuel region change greatly leads to apparent changing of the fast neutron flux and thermal neutron flux distribution in this region. However, for the coolant region, the operation temperature is 673 K-753 K and the macroscopic cross sections of Pb hardly change. So the fast neutron flux and thermal neutron flux distribution on the coolant change a little. Fig. 23 shows the fuel pellet centerline temperature distribution along Z axis direction and there is the peak temperature when Z = 0.6 m. And the maximum fuel pellet centerline temperature deviation is 21 K which occurs on the outlet. Fig. 24 shows the assembly temperature distribution along Y axis direction when Z = 0.6 m. Fig. 25 shows the coolant temperature distribution along Z axis direction (average value on X-Y plane). Fig. 26 displays the cladding outer surface temperature distribution along Z axis direction. It is obvious that the central hole makes the central fuel pellet temperature distribution flat, decreases the peak temperature of fuel pellet and improve the safety margin. And there is the maximum fuel temperature in the central fuel rod of the hot assembly. The coolant outlet temperature is 755.11 K. The maximum fuel temperature is 1643.41 K and the maximum cladding outer surface temperature is 773.83 K calculated by coupling calculation. And thermal conductivity of the MOX (W/m.K)
kMOX ¼ 1:0=ð2:85 0:03 þ 0:035 þ T ð0:286 0:715 0:03ÞÞ þ ð6400=T 2:5 Þ eð16:35=TÞ Fig. 14. Temperature of assembly inlet and outlet (unit: K).
densityMOX
Density of the MOX (kg/m )
8 < ð10970 þ 490 0:05Þ ð9:9734 101 þ T ð9:802 106 þ T ð2:705 1010 þ T 4:391 1013 ÞÞÞ3 ; T < 923 ¼ 3 : ð10970 þ 490 0:05Þ ð9:9672 101 þ T ð1:179 105 þ T ð2:429 1009 þ T 1:219 1012 ÞÞÞ ; T > 923
2946.15 K, Carbajo et al. (2001). Under normal conditions, the maximum cladding temperature should be lower than 823.15 K, Chen (2015) with sufficient safety margin.
ð16Þ
3
ð17Þ
Thermal conductivity of the T91 (W/m.K)
kT91 ¼ 32:196 6:534 eðT=546:792Þ
ð18Þ
10
X. Zhang et al. / Annals of Nuclear Energy 139 (2020) 107243
Temperature(K)
570 560 550 540 530
z/m The calculated axial coolant temperature distribution of adjacent fuel rods center The calculated axial coolant temperature distribution of diagonal fuel rods center The reference axial coolant temperature distribution of diagonal fuel rods center The reference axial coolant temperature distribution of adjacent fuel rods center Fig.15. Axial coolant temperature distribution of the contiguous fuel rod center and diagonal fuel rod center.
1600 1400
Temperature/K
1200 1000 800 600 400 200 0 -0.01
0
0.01
0.02
0.03
0.04
0.05
0.06
X/m The reference temperature distribution along the X direction at z=0.0m and y=0.0m The calculated temperature distribution along the X direction at z=0.0m and y=0.0m Fig.16. Temperature distribution in the X axis direction (z = 0.0 m, y = 0.0 m).
densityPb ¼ 11441 1:2795 T
Specific heat capacity of the T91 (J/kg.K)
C PT91 ¼ 404:336 þ T ð0:881 þ T ð2:04 103 þ T 6
2:6755 10 ÞÞ
Specific heat capacity of the Pb (J/kg.K)
ð19Þ ð20Þ 3
Density of the Pb (kg/m )
C PPb ¼ 176:2 4:923 102 T þ 1:544 105 T T 1:524 106 =ðT TÞ
Thermal conductivity of the Pb (W/m.K)
kPb ¼ 9:2 þ 0:011 T
ð21Þ
ð22Þ
Viscosity of the Pb (Pa.s)
v iscosityPb ¼ 4:55 104 e1069=T
ð23Þ
11
X. Zhang et al. / Annals of Nuclear Energy 139 (2020) 107243
Fig. 19. Radial mesh of the M2LFR-1000 fuel rod.
Fig. 17. Fuel assembly cross section (unit: mm).
Fig. 18. Fuel rod cross section.
the maximum fuel temperature is 1647.17 K and the maximum cladding outer surface temperature is 777.28 K calculated by subchannel code (KMC-SUB) which are all within the corresponding thermal-hydraulics design limits. Fig. 20. Radial mesh of the M2LFR-1000 assembly.
6. Conclusion In this study, based on the UDF and UDS function of the FLUENT, the neutron diffusion equation is defined. There is no requirement to develop the interface program of the coupling calculation. Then the assembly is fine modeled and the solver in the FLUENT is used to solve the neutron diffusion equation. The thermal-hydraulics calculation is carried out at the same time. Therefore the coupling
calculation between neutron diffusion and thermal-hydraulics is achieved on the same solver of the FLUENT. In order to achieve the convenient data transfer, the neutron diffusion and thermalhydraulics calculation use the same meshes. The feasibility of the FLUENT solver based on the Finite Volume Method (FVM) to calculate the neutron diffusion is verified
Table 6 The nuclide composition of the FMS T91. Nuclide
Nucleon density (b1cm1)
Nuclide
Nucleon density (b1cm1)
C Cr Ni Mn Mo Si Nb
3.8900E04 7.8690E03 1.5948E04 3.8300E04 4.6270E04 5.8230E04 4.0200E05
N P S Cu V Al Fe
1.6600E04 3.0200E05 7.2845E06 7.3600E05 1.9700E04 6.9300E05 7.4232E02
12
X. Zhang et al. / Annals of Nuclear Energy 139 (2020) 107243 Table 7 Boundary conditions of the coupling calculation of the M2LFR-1000 assembly. Field
Boundary
Type
Temperature (T)
Side walls Inlet Outlet Inlet Side walls Outlet Inlet
Symmetry Constant value Extrapolation boundary Extrapolation boundary Symmetry Constant value Constant value
Neutron flux (/)
Pressure (P) Velocity (U)
Fig. 21. Fast neutron flux on the outlet (unit: n*s/cm*cm).
through calculating the 2D-TWIGL benchmark problem. The correctness of this coupling method is verified through modeling and calculating the 5 5 PWR assembly Then the coupling method is applied to a modular lead-cooled fast reactor (M2LFR-1000). Through modeling and calculating the hot assembly of the M2LFR-1000: (1) The neutron flux and temperature distribution of the hot assembly are obtained on the steady state. The fast neutron is mainly in the fuel region and the thermal neutron is mainly in the coolant region. And the fast neutron flux and thermal neutron flux distribution on the coolant region change a little. But they change a lot on the fuel region.
Value 673 K Gradient on boundary Gradient on boundary 101 KPa (0,0,1.66)m/s
Fig. 22. Thermal neutron flux on the outlet (unit: n*s/cm*cm).
(2) It is obvious that the central hole makes the central fuel pellet temperature distribution flat, decreases the peak temperature of fuel pellet and improves the safety margin at the same power density. (3) The maximum fuel temperature and the maximum cladding outer surface temperature are obtained by the coupling calculation and are compared with the reference values calculated by sub-channel code (KMC-SUB). The error of maximum fuel temperature is 1.25 K and error of maximum cladding outer surface temperature is 2.68 K. The coolant outlet temperature is 755.11 K which is very close to the design value (753.15 K). And these thermal-hydraulics characteristics are all within the corresponding thermalhydraulics design limits.
13
X. Zhang et al. / Annals of Nuclear Energy 139 (2020) 107243
1700
Temperature/K
1600 1500 1400 1300 1200 1100 1000 900 0
0.2
0.4
0.6
0.8
1
Z/m The reference fuel pellet centerline temperature(KMC-SUB) The calculated fuel pellet centerline temperature Fig. 23. The fuel pellet centerline temperature distribution along Z axis direction.
Fig. 24. Temperature distribution along the Y axis direction (Z = 0.6 m, X = 0.0 m).
1.2
14
X. Zhang et al. / Annals of Nuclear Energy 139 (2020) 107243
Fig. 25. The coolant temperature distribution along Z axis direction (average value on X-Y plane).
Fig. 26. The cladding outer surface temperature distribution along Z axis direction.
X. Zhang et al. / Annals of Nuclear Energy 139 (2020) 107243
Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgment This work is supported by National Natural Science Foundation of China (Grant No. 11705195). Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi.org/10.1016/j.anucene.2019.107243. References Kulesza, J.A., Franceschini, F., Evans, T.M., et al., 2014. Overview of the consortium for the advanced simulation of light water reactor (CASL). CASL-U-2014-0099000. Cacuci, D., 2004. European platform for nuclear reactor simulation (NURESIM). Integrated Project NUCTECH-2004-3. 4. 3. 1-1. EURATOM Research and Training Program on Nuclear Energy. Xie, Z.S., 2004. Physical Analysis of Nuclear Reactor. Xi’an: Xi’an Jiao Tong University press. Tashakor, S., Salehi, A.A., Jahanfarnia, G., Abbaspour, A., Fard, Tehrani, 2012. Neutronic analysis of HPLWR fuel assembly cluster. Ann. Nucl. Energy 50, 38– 43. Liu, X.J., Cheng, X., 2009. Thermal-hydraulic and neutron-physical characteristics of a new SCWR fuel assembly. Ann. Nucl. Energy 36, 28–36. Jiménez, J. et al., 2010. A domain decomposition methodology for pin by pin coupled neutronic and thermal–hydraulic analyses in COBAYA3. Nucl. Eng. Des. 240, 313–320. Zhao, C., et al., 2016. Coupled neutronics and thermal–hydraulics analysis of annular fuel assembly for SCWR. In: Jiang H. (Eds.), Proceedings of The 20th Pacific Basin Nuclear Conference. PBNC 2016. Springer, Singapore. Jiménez, J., et al., 2008. Comparative analysis of neutronics/thermal-hydraulics multi-scale coupling for LWR analysis. In: International Conference on the Physics of Reactors ‘‘Nuclear Power: A Sustainable Resource” Casino-Kursaal Conference Center, Interlaken, Switzerland, September 14–19.
15
Wang, H.G., et al., 2008. Investigation of batch fluidized-bed drying by mathematical modeling, CFD simulation and ECT measurement. Wiley J., 54, 427–444 doi: 10.1002/aic.11406. Donoso-GarcíaL, P., Henríquez-Vargas, L., 2015. Numerical study of turbulent porous media combustion coupled with thermoelectric generation in a recuperative reactor. Energy 93, 1189–1198. https://doi.org/10.1016/j. energy.2015.09.123. Liu, Y. et al., 2014. Three-dimensional analysis of gas flow and heat transfer in a regenerator with alumina balls. Appl. Thermal Eng. 69, 113–122. https://doi. org/10.1016/j.applthermaleng.2014.04.058. Jang, J., Arastoopour, H., 2014. CFD simulation of a pharmaceutical bubbling bed drying process at three different scales. Powder Technol. 263, 14–25. https:// doi.org/10.1016/j.powtec.2014.04.054. Zhang, D.L. et al., 2009. Development of a steady state analysis code for a molten salt reactor. Ann. Nucl. Energy 36, 590–603. Gui, X.W. et al., 2010. Study on coupling of local three-dimension flow model based on CFD method and space-time neutron kinetics model. Chin. J. Nucl. Sci. Eng. 3, 216–222. Klas, J. et al., 2014. Fine-mesh deterministic modeling of PWR fuel assemblies: proof-of-principle of coupled neutronic/thermal–hydraulic calculations. Ann. Nucl. Energy 68, 247–256. Chen, H.L. et al., 2018. Preliminary design of a medium-power modular lead-cooled fast reactor with the application of optimization methods. Int. J. Energy Res. 42, 3643–3657. Li, S.Z. et al., 2017. Development of a sub-channel thermal hydraulic analysis code and its application to lead cooled fast reactor. Appl. Therm. Eng. 117, 443–451. Pirouzmand, Ahmad, Nabavi, Abolhasan, 2016. Simulation of nuclear reactor dynamics equations using reconfigurable computing. Prog. Nucl. Energy 89, 197–203. Ge, J. et al., 2015. Steady and transient solutions of neutronics problems based on finite volume method(FVM) with a CFD code. Prog. Nucl. Energy 85, 366–374. ANSYS, 2013. ANSYS FLUENT Theory Guide, Release 15.0. Zhou, X.B. et al., 2018. Development and preliminary test of data library ANDL-ADS for accelerator-driven systems. Nucl. Tech. 41 (03), 65–70. Bi, G.W., 2008. Interpolation Method Development for Temperature Based Neutron Cross-sections. Tsinghua University, Beijing. David Jaluvka, 2015. Development of a Core Management Tool for the MYRRHA Irradiation Research Facility [D]. KU Leuven. Carbajo, J.J., Yoder, G.L., Popov, S.G., Ivanov, V.K., 2001. A review of the thermophysical properties of MOX and UO2 fuels. J. Nucl. Mater. 299 (3), 181–198. Chen, Z., 2015. Thermal-hydraulics design and safety analysis of a 100MWth small natural circulation lead cooled fast reactor SNCLFR-100. University of Science and Technology of China. Popov, S.G., Carbajo, J.J., Ivanov, V.K., Yoder, G.L., 2000. Thermophysical properties of MOX and UO2 fuels including the effects of irradiation. ORNL Report TM-2000/ 351.