# Analysis of the Natural Convection Flow in the Upper Plenum of the MONJU Reactor with Trio_U

Authors: Ulrich Bieder; Gauthier Fauchet
Journal: Science and Technology of Nuclear Installations
 (2013)

---

## Abstract

The IAEA has coordinated a benchmark project on natural convection phenomena in the upper plenum of the MONJU reactor. JAEA has provided both detailed geometrical data of the plant and complete thermalhydraulic boundary conditions describing a pump trip transient, accomplished during the start-up experiments of the reactor. For the initial conditions of the pump trip transient, extensive sensitivity analyses have been made with the CFD code Trio_U. These calculations show a high sensitivity of the global flow pattern in the MONJU upper plenum depending on the initial order of the numerical scheme and the modelling of the geometrically complex upper core structure. During the pump trip, the formation of a thermal stratification within the plenum has been observed which persists for almost two hours. All calculations have shown a homogenization of the temperature in the plenum after about 15 minutes. A slight reduction of the mixing in the upper plenum could have been achieved by modifying the form of the flow holes in the inner barrel (fillets instead of sharp edges) in order to reduce their axial pressure loss.


---

## Body

## 1. Introduction

The IAEA has coordinated a research project (CRP) between 2008 and 2012 entitled “Benchmark Analysis of sodium Natural Convection in the upper plenum of the MONJU Reactor Vessel.” Eight research organizations from seven countries with an active program on sodium cooled fast reactors—namely, China, France, India, Japan, Republic of Korea, Russian Federation, and USA—contributed to this CRP. Japan Atomic Energy Agency (JAEA) has submitted to the CRP participants the data of sodium thermal stratification measurements in the MONJU reactor vessel upper plenum collected during a plant trip test conducted in December 1995. The benchmark partners have analysed this experiment by applying different codes and methodologies. The benchmark thus helped the members to improve their capability in the field of fast reactor in-vessel Sodium thermalhydraulics.Sodium cooled fast breeder reactors are under development for more than 50 years. Nevertheless, only very limited data are published to date which allow the validation of CFD codes in general and for natural and mixed convection phenomena in particular.Within the European Fast Breeder Reactor project, an experimental approach in the RAMONA facility has been largely used to study decay heat removal situations [1]. RAMONA was a 1/20 scale water model of the upper plenum with boundary conditions imposed at the core outlet and with active immersed coolers. Transient situations were tested and the interaction between cold sodium flowing downwards from the coolers and hot sodium flowing upwards from the core was analysed in detail. A larger water model at a scale 1/5 called NEPTUN was also used to estimate the scale effect on the results for a better extrapolation to the reactor [2].Similarity criteria are often used to extrapolate the results of water tests to sodium flow. Such similarity criteria for mixed convection have been evaluated at CEA in two geometrically identical experimental facilities; SUPERCAVNA for water flow and CORMORAN for sodium flow. The generic setup was a horizontal channel connected to a cavity, heated (or cooled) at one wall. The thermal stratification and thermal fluctuations have been measured in both facilities and additional velocity measurements were made in the water tests [3]. Sodium and water results have been compared to verify the similarity criteria, and the tests have been used to validate CFD calculations [4]. It is important to note that strong temperature gradients can lead to anisotropic turbulence. This effect can significantly influence the heat exchanges and cannot be modelled by linear eddy viscosity models [5].For the mixing of submerged jets encountered at the core outlet, two identical setups have been analysed in [6]: one for air and the other one for sodium. A basic configuration of two coaxial jets at different temperatures was used. The evaluation of the transfer of temperature fluctuations from the fluid to the structure requires experimental data in sodium, as the high thermal diffusivity of sodium plays an important role in the transfer function. Recent sodium data are provided by JAEA on a sodium test facility with jets mixing near a wall [7].Accompanying the construction of the Indian prototype fast breeder reactor, thermalhydraulic phenomena of internal heat exchangers have been analysed in detail by using CFD codes. On the primary side, the uniform distribution of the flow at the inlet window connected to the upper plenum has been evaluated in [8]. The flow and the temperature fields in the bundle also play an important role in the efficiency of the heat exchanger [8].Most of the phenomena expected to be important in the thermal stratification formation in the upper plenum of sodium cooled reactors have been analysed experimentally by separate effect tests. The only full scale, integral test published to date is the MONJU reactor pump trip used in the IAEA CRP. These data are used in this paper and are applied for CFD validation on reactor scale.


## 2. The MONJU Fast Breeder Reactor

### 2.1. The Upper Plenum of MONJU

The MONJU plant is a prototype nuclear power plant with sodium cooled fast breeder reactor. The upper plenum of the MONJU reactor is geometrically very complex. Figure1 is a view of the complete upper plenum from the support plate to the liquid level [9, 10]. The sodium is leaving the sub-assemblies (SA) of the core, is traversing the topologically very complex tube bundles of the upper core structure (UCS) to enter the open volume of the upper plenum. Then, the flow either crosses over the inner barrel (IB) or passes through the flow holes (FH) of the inner barrel to reach the three outlet nozzles (ON).Figure 1
Sketch of the MONJU upper plenum.The Fuel Handling Machine (FHM) with its hold down arm is also illustrated in Figure1, as well as the thermocouple plug TCP with its 38 thermocouples. These thermocouples are installed to measure the vertical temperature distribution in the upper plenum.The UCS is situated above the SA outlets and below the Upper Instrument Structure body (UIS). This region is shown in more detail in Figure2.Figure 2
Sketch of the MONJU upper core structure (UCS).The Flow Guide Tubes (FGTs) are a tube bundle maintained at its top by a perforated plate called Honeycomb Structure (HS). The Fingers (FS) are located above the FGT and the HS. Thermocouples are positioned at the FGT outlet centres, installed at the bottom ends of FS which are connected to the UIS. The main flow passes from the SA into the FGT to enter the Fingers region. 19 Control Rod Guide Tubes (CRGTs) are present inside both regions.


### 2.2. The MONJU Pump Trip Test

The purpose of the test was the confirmation of plant safety against turbine failure. At 40% nominal power, the main coolant pumps are tripped and the control rods were inserted. The decay heat is removed to the atmosphere through sodium-air heat exchangers. The temporal changes of the vertical distribution of the sodium temperature in the upper plenum was measured along the TCP as shown in Figure3. The reference height 0 m corresponds to the elevation of the liquid level.Figure 3
Measured temperature distribution on TCP.The thermal stratification front moved upward after the reactor scram and reached the top of the IB after about two hours.


## 2.1. The Upper Plenum of MONJU

The MONJU plant is a prototype nuclear power plant with sodium cooled fast breeder reactor. The upper plenum of the MONJU reactor is geometrically very complex. Figure1 is a view of the complete upper plenum from the support plate to the liquid level [9, 10]. The sodium is leaving the sub-assemblies (SA) of the core, is traversing the topologically very complex tube bundles of the upper core structure (UCS) to enter the open volume of the upper plenum. Then, the flow either crosses over the inner barrel (IB) or passes through the flow holes (FH) of the inner barrel to reach the three outlet nozzles (ON).Figure 1
Sketch of the MONJU upper plenum.The Fuel Handling Machine (FHM) with its hold down arm is also illustrated in Figure1, as well as the thermocouple plug TCP with its 38 thermocouples. These thermocouples are installed to measure the vertical temperature distribution in the upper plenum.The UCS is situated above the SA outlets and below the Upper Instrument Structure body (UIS). This region is shown in more detail in Figure2.Figure 2
Sketch of the MONJU upper core structure (UCS).The Flow Guide Tubes (FGTs) are a tube bundle maintained at its top by a perforated plate called Honeycomb Structure (HS). The Fingers (FS) are located above the FGT and the HS. Thermocouples are positioned at the FGT outlet centres, installed at the bottom ends of FS which are connected to the UIS. The main flow passes from the SA into the FGT to enter the Fingers region. 19 Control Rod Guide Tubes (CRGTs) are present inside both regions.


## 2.2. The MONJU Pump Trip Test

The purpose of the test was the confirmation of plant safety against turbine failure. At 40% nominal power, the main coolant pumps are tripped and the control rods were inserted. The decay heat is removed to the atmosphere through sodium-air heat exchangers. The temporal changes of the vertical distribution of the sodium temperature in the upper plenum was measured along the TCP as shown in Figure3. The reference height 0 m corresponds to the elevation of the liquid level.Figure 3
Measured temperature distribution on TCP.The thermal stratification front moved upward after the reactor scram and reached the top of the IB after about two hours.


## 3. The Numerical Model

### 3.1. The Simplified Geometry

A simplified CAD model has been developed at ANL [11] which consists in angular direction of 1/6 of the upper plenum of the MONJU reactor. The following modifications have been made in order to obtain a symmetric geometry.(i)
All of the SA outlets are taken as hexagonal.(ii)
The top of the SA outlets is placed on the levels 0.87 m, 0.9 m, and 0.93 m above the support plate to distinguish the different channels easily. In /10/, all SA outlets are located 0.93 m above the support plate.(iii)
The ON is rotated to respect the hexagonal core configuration (only half of the outlet nozzle is modelled).(iv)
Internal structures like FHM and TCP which are placed in a nonsymmetrical way in the upper plenum are not modelled.Within the UCS, only the Control Rod Guide Tubes (CRGTs) are realized in the CAD model. FGT and FS are not treated explicitly due to their small size and geometrical complexity. However, the external frontiers of the FGT and Fingers regions are realized with the intention of forming in that way two subregions which allow a specific macroscopic modelling by means of porous media with directional pressure loss correlations.


### 3.2. Modelling of the Upper Core Structure

The FGT and FS regions are densely packed tube bundles. The correlation (1) defining the friction coefficient Cf, which is used to predict the pressure loss within these regions, is valid for flow inside and outside of tube bundles with a circular pitch [12]
(1)Cf=aRe-bwithRe=UDν.
Re is the Reynolds number, U a characteristic velocity, D a characteristic diameter, and ν the kinematic viscosity. The correlation distinguishes the axial (subscript a) and the transverse (subscript t) direction of the velocity vector u→ with respect to the orientation to the tubes as well as the hydraulic diameter Dh and external diameter of tubes De. The values of the parameters of (1) are summarised in Table 1.Table 1
Parameters for the pressure loss correlation.
Direction
a
b
U
D
Axial
0.316
0.25
|
u
→
a
|
D
h
Transverse
4.03
0.27
|
u
→
t
|
D
eThe axial velocity and transversal velocity are calculated as a function of tube direction vectord→a which is simplified as the global vertical direction (0,0,1)T:
(2)u→a=(u→·d→a)d→a,u→t=u⃑-u→a.
The pressure loss of the tube bundle ∇→PTB is implemented in the Navier-Stokes equations as additional source term:
(3)∇→PTB=-ρ·Cfa|u→a|2·Dhu→a-ρ·Cft|u→t|2·Deu→t,
where ρ is the fluid density. It is important to verify that the implementation of the transversal pressure loss leads to angle-independent results for flow source in the centre of a symmetric rod bundle.The HS is modelled as perforated plate without any thickness. The singular pressure loss which is applied on the plate is predicted from a grid pressure loss correlation [12]. The directional pressure loss applied normal to the HS plate is defined as
(4)ΔPHS=12·ρ·K·|u→|·(u→·n→),
where K is the pressure loss coefficient of the HS and n→ the normal vector of the surface of the HS plate. The value of the pressure loss coefficient K is constant during the trip test and is equal to 60 [13].The tubes in the fingers region reduce the volume accessible by the sodium flow. A volumetric porosity of 0.83 is thus defined in the Navier-Stokes equations for the FS region to simulate the acceleration of the flow outside of the tubes due to mass conservation considerations.


### 3.3. Meshing

ANL has provided the CAD model of the simplified 60° geometry in IGES format. This model has been imported into the commercial mesh generator ICEMCFD. The UCS regions FGT and FS are defined as internal sub-domains, and the HS is implemented as an internal, semi permeable wall. Three pure tetrahedral meshes of 300.000 elements (coarse mesh) and 1.25 million elements (intermediate mesh) and 3.3 million elements (fine mesh) have been created. The coarse and intermediate meshes are shown in Figure4. A horizontal cut plane is visualized which goes through the lower line of the core barrel flow holes.Coarse and intermediate meshing in a horizontal cut plan.
(a)
(b)The coarse and fine meshes are created with a Delaunay procedure which leads to a better approximation of boundary layers on solid surfaces. The intermediate mesh is created with an Octree procedure which leads to a more homogeneous mesh inside of the plenum. Where possible, a similar thickness of the first wall near layer of cells has been defined for all meshes. They+ values on the plenum walls, the IB, and the UIS are large for the steady-state solution (order of 1000) and reach values around 100 during the transient. The flow holes are discretized in about 8 calculation points per diameter for the coarse mesh and in about 15 calculation points for the fine mesh.


### 3.4. Thermophysical Properties of Sodium

Liquid sodium is used as coolant in the MONJU reactor. The default values of the thermo-physical data used at CEA for 400°C are given in Table2.Table 2
Physical properties of sodium at 400°C.
Quantity
Value
Unit
Density
858
kg·m−3
Dynamic viscosity
2.81 * 10−4
kg·m−1·s−1
Therm. conductivity
69.7
W·m−1·K−1
Heat capacity
1284
J·kg−1·K−1
Therm. expansion coefficient
2.68 * 10−4
K−1In the temperature range of interest (380°C to 512°C), the fluid can be treated in first approximation as incompressible. This simplification allows the application of the Boussinesq hypothesis to simplify the Navier-Stokes equations. All physical properties are treated as temperature independent, and buoyancy effects are taken into account only via the gravitational acceleration.


### 3.5. The Trio_U Code and the Numerical Scheme

Trio_U [14] is a CFD code for strongly unsteady, low Mach number, turbulent flows. The code is especially designed for industrial CFD calculations on structured and nonstructured grids of several hundreds of millions of nodes. The platform-independent code, developed at the CEA, is based on an object-oriented, intrinsically parallel approach and is coded in C++. The flexible code structure allows the user to choose a suitable discretization method and to combine various appropriate physical models, including different treatments of turbulence. Several convection and time marching schemes as well as a wide range of boundary conditions are available. This flexibility is implemented for massively parallel computing without a significant reduction of the overall performance of the code.For unstructured, tetrahedral grids, a hybrid Finite Volume Element method (FVE) is applied. This method approximates a continuous problem by a discrete solution in the space of the finite elements by maintaining the balance notation of finite volumes. In Trio_U, the main unknowns as velocity and temperature are located in the centre of the faces of an element (P1NC). Thus, the number of control volumes for the momentum and scalar conservation is approximately two times the number of elements. The pressure is discretized in both the centre (P0) and the vertices (P1) of an element. The resulting staggered mesh arrangement improves the velocity/pressure coupling. An implicit velocity projection method is used to assure the mass conservation [15]. The fundamental numerical scheme of the calculations presented in this paper is summarised in Table 3. More detailed information on the numerical scheme and the solution method can be found in [14].Table 3
Fundamental numerical scheme used in the Trio_U calculations.
General
Dimension
3D calculation
Fluid
Sodium at 400°C
Mesh
Tetrahedral mesh
Discretization
P0/P1 for PP1NC for U,T,k,ε
Time scheme
1st order Euler explicit
Navier-Stokesequations
Convection
2nd order upwind
Diffusion
2nd order centred
Pressure solver
Cholesky method
Thermal effects
Boussinesq hypothesis
Wall law
Logarithmic wall law
Turbulence
RANS
Turbulencemodelling
Turbulence model
High Reynoldsk-ε
k, ε convection
2nd order upwind
k, ε diffusion
2nd order centred
Energy transportequation
Convection
2nd order upwind
Diffusion
2nd order centred
Wall law
Logarithmic wall law
Turbulence
Turbulent Prandtl NumberPrtTo reach a steady-state solution, a transient is calculated until all unknowns reach constant values. This procedure guarantees physically correct solutions at any time of the transient.


### 3.6. Turbulence Modelling

Turbulence is treated with the Reynolds averaging concept where the instantaneous velocityu→ is decomposed into a temporal mean value u-→ and a fluctuating value u→′. The set of Reynolds-averaged Navier-Stokes equations is closed by using the k-ε model. The high Reynolds form of the model is appropriated to fully developed turbulent flows and allows to some extent the presence of buoyancy effects [16]. In the Boussinesq hypothesis framework, the turbulent viscosity is linked to the turbulent kinetic energy k and the dissipation rate of the turbulent kinetic energy ε via
(5)νt=Cμk2ε.
Conservation equations are written for both the turbulent kinetic energy k and the turbulent dissipation rate ε: (6a)∂k∂t+u-→·∇→k=∇→·(νtσk∇→k)-ε+P+G,(6b)∂ε∂t+u-→·∇→ε=∇→·(νtσε∇→ε)-Cε2ε2k+Cε1Pεk+Cε1Cε3Gεk.The production of turbulence kinetic energy is calculated by(7)P=-ui′uj′¯∂u-i∂xjwith-ui′uj′¯=νt(∂u-i∂xj+∂u-j∂xi)-23kδij.Buoyancy effects for incompressible flows are treated by(8)G=-νtPrtβTg→·∇→T-.The empirical coefficients listed in Table4 are used.Table 4
Coefficients of the high Reynolds numberk-ε model.
C
μ
σ
k
σ
ε
C
ε
1
C
ε
2
C
ε
3
Pr
t
0.09
1.0
1.3
1.44
1.92
1.0
0.9To take into account the effect of thermal stratification, the following assumption is made [16]:(i)
for stable stratification with reduced buoyancy effects:Cε3=0ifG<0;(ii)
For unstable stratification with full buoyancy effects:Cε3=1ifG>0.This extension of the standard k-ε model has been tested for sodium flow by analysing various natural and mixed convection experiments [4].


### 3.7. Boundary Conditions

For each SA outlet, the time-depending mass flow rate and temperature are given in the data description report [17] for selected instances during the whole pump trip experiment. Continuous values are interpolated linearly between reported ones. The mass flow rates are converted to velocity inflow boundary conditions by applying on one hand the physical properties listed in Table 2 and on the other hand the free surface of an SA outlet (11573 mm2). A spatially uniform distribution of the temperature and the velocity is assumed for each SA outlet. The turbulent kinetic energy k and its dissipation rate ε are predicted from the assembly diameter and a fluctuation velocity which is assumed to be 10% of the mean SA outflow velocity.All walls are treated as adiabatic. Standard logarithmic wall functions [16] are applied for momentum equations. The free sodium surface and the two lateral surfaces of the 60° geometry are modelled with the symmetry hypothesis. A pressure outlet condition, associated with free out stream conditions for the velocity and the temperature is used as boundary condition for the outlet nozzle boundary.


## 3.1. The Simplified Geometry

A simplified CAD model has been developed at ANL [11] which consists in angular direction of 1/6 of the upper plenum of the MONJU reactor. The following modifications have been made in order to obtain a symmetric geometry.(i)
All of the SA outlets are taken as hexagonal.(ii)
The top of the SA outlets is placed on the levels 0.87 m, 0.9 m, and 0.93 m above the support plate to distinguish the different channels easily. In /10/, all SA outlets are located 0.93 m above the support plate.(iii)
The ON is rotated to respect the hexagonal core configuration (only half of the outlet nozzle is modelled).(iv)
Internal structures like FHM and TCP which are placed in a nonsymmetrical way in the upper plenum are not modelled.Within the UCS, only the Control Rod Guide Tubes (CRGTs) are realized in the CAD model. FGT and FS are not treated explicitly due to their small size and geometrical complexity. However, the external frontiers of the FGT and Fingers regions are realized with the intention of forming in that way two subregions which allow a specific macroscopic modelling by means of porous media with directional pressure loss correlations.


## 3.2. Modelling of the Upper Core Structure

The FGT and FS regions are densely packed tube bundles. The correlation (1) defining the friction coefficient Cf, which is used to predict the pressure loss within these regions, is valid for flow inside and outside of tube bundles with a circular pitch [12]
(1)Cf=aRe-bwithRe=UDν.
Re is the Reynolds number, U a characteristic velocity, D a characteristic diameter, and ν the kinematic viscosity. The correlation distinguishes the axial (subscript a) and the transverse (subscript t) direction of the velocity vector u→ with respect to the orientation to the tubes as well as the hydraulic diameter Dh and external diameter of tubes De. The values of the parameters of (1) are summarised in Table 1.Table 1
Parameters for the pressure loss correlation.
Direction
a
b
U
D
Axial
0.316
0.25
|
u
→
a
|
D
h
Transverse
4.03
0.27
|
u
→
t
|
D
eThe axial velocity and transversal velocity are calculated as a function of tube direction vectord→a which is simplified as the global vertical direction (0,0,1)T:
(2)u→a=(u→·d→a)d→a,u→t=u⃑-u→a.
The pressure loss of the tube bundle ∇→PTB is implemented in the Navier-Stokes equations as additional source term:
(3)∇→PTB=-ρ·Cfa|u→a|2·Dhu→a-ρ·Cft|u→t|2·Deu→t,
where ρ is the fluid density. It is important to verify that the implementation of the transversal pressure loss leads to angle-independent results for flow source in the centre of a symmetric rod bundle.The HS is modelled as perforated plate without any thickness. The singular pressure loss which is applied on the plate is predicted from a grid pressure loss correlation [12]. The directional pressure loss applied normal to the HS plate is defined as
(4)ΔPHS=12·ρ·K·|u→|·(u→·n→),
where K is the pressure loss coefficient of the HS and n→ the normal vector of the surface of the HS plate. The value of the pressure loss coefficient K is constant during the trip test and is equal to 60 [13].The tubes in the fingers region reduce the volume accessible by the sodium flow. A volumetric porosity of 0.83 is thus defined in the Navier-Stokes equations for the FS region to simulate the acceleration of the flow outside of the tubes due to mass conservation considerations.


## 3.3. Meshing

ANL has provided the CAD model of the simplified 60° geometry in IGES format. This model has been imported into the commercial mesh generator ICEMCFD. The UCS regions FGT and FS are defined as internal sub-domains, and the HS is implemented as an internal, semi permeable wall. Three pure tetrahedral meshes of 300.000 elements (coarse mesh) and 1.25 million elements (intermediate mesh) and 3.3 million elements (fine mesh) have been created. The coarse and intermediate meshes are shown in Figure4. A horizontal cut plane is visualized which goes through the lower line of the core barrel flow holes.Coarse and intermediate meshing in a horizontal cut plan.
(a)
(b)The coarse and fine meshes are created with a Delaunay procedure which leads to a better approximation of boundary layers on solid surfaces. The intermediate mesh is created with an Octree procedure which leads to a more homogeneous mesh inside of the plenum. Where possible, a similar thickness of the first wall near layer of cells has been defined for all meshes. They+ values on the plenum walls, the IB, and the UIS are large for the steady-state solution (order of 1000) and reach values around 100 during the transient. The flow holes are discretized in about 8 calculation points per diameter for the coarse mesh and in about 15 calculation points for the fine mesh.


## 3.4. Thermophysical Properties of Sodium

Liquid sodium is used as coolant in the MONJU reactor. The default values of the thermo-physical data used at CEA for 400°C are given in Table2.Table 2
Physical properties of sodium at 400°C.
Quantity
Value
Unit
Density
858
kg·m−3
Dynamic viscosity
2.81 * 10−4
kg·m−1·s−1
Therm. conductivity
69.7
W·m−1·K−1
Heat capacity
1284
J·kg−1·K−1
Therm. expansion coefficient
2.68 * 10−4
K−1In the temperature range of interest (380°C to 512°C), the fluid can be treated in first approximation as incompressible. This simplification allows the application of the Boussinesq hypothesis to simplify the Navier-Stokes equations. All physical properties are treated as temperature independent, and buoyancy effects are taken into account only via the gravitational acceleration.


## 3.5. The Trio_U Code and the Numerical Scheme

Trio_U [14] is a CFD code for strongly unsteady, low Mach number, turbulent flows. The code is especially designed for industrial CFD calculations on structured and nonstructured grids of several hundreds of millions of nodes. The platform-independent code, developed at the CEA, is based on an object-oriented, intrinsically parallel approach and is coded in C++. The flexible code structure allows the user to choose a suitable discretization method and to combine various appropriate physical models, including different treatments of turbulence. Several convection and time marching schemes as well as a wide range of boundary conditions are available. This flexibility is implemented for massively parallel computing without a significant reduction of the overall performance of the code.For unstructured, tetrahedral grids, a hybrid Finite Volume Element method (FVE) is applied. This method approximates a continuous problem by a discrete solution in the space of the finite elements by maintaining the balance notation of finite volumes. In Trio_U, the main unknowns as velocity and temperature are located in the centre of the faces of an element (P1NC). Thus, the number of control volumes for the momentum and scalar conservation is approximately two times the number of elements. The pressure is discretized in both the centre (P0) and the vertices (P1) of an element. The resulting staggered mesh arrangement improves the velocity/pressure coupling. An implicit velocity projection method is used to assure the mass conservation [15]. The fundamental numerical scheme of the calculations presented in this paper is summarised in Table 3. More detailed information on the numerical scheme and the solution method can be found in [14].Table 3
Fundamental numerical scheme used in the Trio_U calculations.
General
Dimension
3D calculation
Fluid
Sodium at 400°C
Mesh
Tetrahedral mesh
Discretization
P0/P1 for PP1NC for U,T,k,ε
Time scheme
1st order Euler explicit
Navier-Stokesequations
Convection
2nd order upwind
Diffusion
2nd order centred
Pressure solver
Cholesky method
Thermal effects
Boussinesq hypothesis
Wall law
Logarithmic wall law
Turbulence
RANS
Turbulencemodelling
Turbulence model
High Reynoldsk-ε
k, ε convection
2nd order upwind
k, ε diffusion
2nd order centred
Energy transportequation
Convection
2nd order upwind
Diffusion
2nd order centred
Wall law
Logarithmic wall law
Turbulence
Turbulent Prandtl NumberPrtTo reach a steady-state solution, a transient is calculated until all unknowns reach constant values. This procedure guarantees physically correct solutions at any time of the transient.


## 3.6. Turbulence Modelling

Turbulence is treated with the Reynolds averaging concept where the instantaneous velocityu→ is decomposed into a temporal mean value u-→ and a fluctuating value u→′. The set of Reynolds-averaged Navier-Stokes equations is closed by using the k-ε model. The high Reynolds form of the model is appropriated to fully developed turbulent flows and allows to some extent the presence of buoyancy effects [16]. In the Boussinesq hypothesis framework, the turbulent viscosity is linked to the turbulent kinetic energy k and the dissipation rate of the turbulent kinetic energy ε via
(5)νt=Cμk2ε.
Conservation equations are written for both the turbulent kinetic energy k and the turbulent dissipation rate ε: (6a)∂k∂t+u-→·∇→k=∇→·(νtσk∇→k)-ε+P+G,(6b)∂ε∂t+u-→·∇→ε=∇→·(νtσε∇→ε)-Cε2ε2k+Cε1Pεk+Cε1Cε3Gεk.The production of turbulence kinetic energy is calculated by(7)P=-ui′uj′¯∂u-i∂xjwith-ui′uj′¯=νt(∂u-i∂xj+∂u-j∂xi)-23kδij.Buoyancy effects for incompressible flows are treated by(8)G=-νtPrtβTg→·∇→T-.The empirical coefficients listed in Table4 are used.Table 4
Coefficients of the high Reynolds numberk-ε model.
C
μ
σ
k
σ
ε
C
ε
1
C
ε
2
C
ε
3
Pr
t
0.09
1.0
1.3
1.44
1.92
1.0
0.9To take into account the effect of thermal stratification, the following assumption is made [16]:(i)
for stable stratification with reduced buoyancy effects:Cε3=0ifG<0;(ii)
For unstable stratification with full buoyancy effects:Cε3=1ifG>0.This extension of the standard k-ε model has been tested for sodium flow by analysing various natural and mixed convection experiments [4].


## 3.7. Boundary Conditions

For each SA outlet, the time-depending mass flow rate and temperature are given in the data description report [17] for selected instances during the whole pump trip experiment. Continuous values are interpolated linearly between reported ones. The mass flow rates are converted to velocity inflow boundary conditions by applying on one hand the physical properties listed in Table 2 and on the other hand the free surface of an SA outlet (11573 mm2). A spatially uniform distribution of the temperature and the velocity is assumed for each SA outlet. The turbulent kinetic energy k and its dissipation rate ε are predicted from the assembly diameter and a fluctuation velocity which is assumed to be 10% of the mean SA outflow velocity.All walls are treated as adiabatic. Standard logarithmic wall functions [16] are applied for momentum equations. The free sodium surface and the two lateral surfaces of the 60° geometry are modelled with the symmetry hypothesis. A pressure outlet condition, associated with free out stream conditions for the velocity and the temperature is used as boundary condition for the outlet nozzle boundary.


## 4. Analysis of the Steady-State Condition

Before the pump trip at 40% power, the temperature field and the velocity flow were well established in the upper plenum. A first calculation was performed in order to obtain this initial condition. The boundary conditions for the calculation are the conditions at the beginning of the test [17]. Starting from an initially guessed velocity and temperature field, a transient calculation was carried out for constant boundary conditions to reach a steady-state solution. The flow field is considered as established in the plenum when all 120 temperature samples placed axially along the thermocouple plug (TCP) have reached a constant value.A reduced Froude number Fr of the jet leaving the upper core structure is calculated from(9)Fr=|u-→||g→|·Δρ/ρ·d.
Using the mean velocity of the jet (|u-→|≈0.5 m/s), the density difference between hot jet and cooler environment (Δρ/ρ≈3·10-3), and the mean diameter of the jet (d≈0.5 m), the reduced Froude number is in the order of 4. This reduced Froude number characterises mixed convection what indicates that both buoyancy forces and the inertia of the jet are in the same order of magnitude.


### 4.1. Momentum Dominated Solution

A three-step procedure was applied to ensure convergence of the solution on both the meshing and the numerical scheme. In this procedure, solutions converged on coarser meshes were interpolated on finer meshes to define the new initial condition for the subsequent calculation on the finer mesh.(i)
Initializing with a reposing (u→=0) isothermal flow at T=418°C, a preliminary calculation was performed on the coarse mesh using a 1st order convection scheme. This calculation led to solution S1c.(ii)
Initializing with the solutionS1c, a reference calculation was performed on the intermediate mesh using a 2nd order convection scheme. This calculation led to solution S1i.(iii)
Initializing with the solutionS1i, a convergence test calculation was performed on the fine mesh using a 2nd order convection scheme. This calculation led to solution S1f.Identical solutionsS1i and S1f indicate a mesh independency of the solution which is based on a 2nd order numerical scheme.In Figure5(a), both the temperature distribution and the flow field are given for the symmetry plane going through the outlet nozzle. The solution S1i is shown. The reference height 0 m corresponds to the elevation of the support plate. The momentum of the flow leaving the upper core structures is sufficiently high to allow the jet to traverse the upper plenum, to impinge on the outer barrel, and to cross over the outer barrel to finally reach the outlet nozzle. We call this solution “momentum dominated solution.”(a) Momentum dominated solution (S1i). (b) Buoyancy dominated solution (S2i). (c) Mixed convection solution (S3i).
(a)
(b)
(c)


### 4.2. Buoyancy Dominated Solution

A two-step procedure was applied to ensure convergence of the solution on the meshing.(i)
Initializing with a reposing (u→=0) isothermal flow at T=418°C, the reference calculation was performed on the intermediate mesh using a 2nd order convection scheme. This calculation led to solution S2i.(ii)
Initializing with the solutionS2i, a convergence test calculation was performed on the fine mesh using a 2nd order convection scheme. This calculation led to solution S2f.Identical solutionsS2i and S2f indicate a mesh independency of the solution which is based on a 2nd order numerical scheme.In Figure5(b), both the temperature distribution and the flow field are given for the symmetry plane going through the outlet nozzle. The solution S2i is shown. The reference height 0 m corresponds to the elevation of the support plate. The buoyancy of the flow leaving the upper core structures prevents the jet from traversing the upper plenum as in solution S1. The flow leaving the upper core structure rises along the upper instrument structure and crosses over the outer barrel to reach the outlet nozzle. We call this solution “buoyancy dominated solution.”It has been tested if it is possible to change the solution from the buoyancy dominated flow to the momentum driven flow by changing the numerical scheme. Initializing with the solutionS2f, a test calculation was performed on the coarse mesh using a 1st order convection scheme. This calculation maintained the solution S2 after convergence.The coldest sodium leaving the core region has a temperature of 418°C. When initialising the upper plenum with this temperature, buoyancy forces act instantly on the hotter jet. These buoyancy forces are minimized by initialising the upper plenum at the temperature of the hottest sodium (513°C) leaving the core region. Starting from a reposing (u→=0) isothermal flow at T=513°C, an additional calculation was performed on the intermediate mesh using a 2nd order convection scheme. This calculation also led to solution S2i.


### 4.3. Mixed Convection Solution

The possible influence of the upper core structure (UCS) on both the thermal stratification and the flow pattern within the upper plenum is not a priori evident. To show the possible effect of the UCS on the flow, an additional sensitivity calculation without UCS has been performed.Initializing with the solutionS2i (with UCS), a test calculation without UCS was performed on an intermediate mesh using a 2nd order convection scheme. This calculation led to solution S3i.Convergence tests on further mesh refinements have not been made. The solutionS3i is shown in Figure 5(c). The temperature distribution and the flow field are given for the symmetry plane going through the outlet nozzle. The reference height 0 m corresponds to the elevation of the support plate. The buoyancy of the flow leaving the upper core structures is sufficiently high to prevent the jet from traversing the upper plenum as in solution S1. However, the inertia is high enough to also prevent the jet from simply rising along the Upper Instrument Structure as in solution S2. After spreading out in the upper plenum, the flow crosses over the outer barrel and reaches the outlet nozzle. We call this solution “mixed convection solution.”


### 4.4. Lessons Learned from the Steady-State Calculations

Using different modelling options, three different steady solutions have been achieved which can be distinguished regarding the driving physical process:(i)
a momentum dominated solution S1;(ii)
a buoyancy dominated solution S2;(iii)
a mixed convection solution S3.The calculated temperature profiles along the TCP are compared in Figure6 to the experimental values. The reference height 0 m corresponds to the elevation of the support plate.Figure 6
Comparison of the axial temperature profiles.From these temperature profiles it seems not possible to identify the flow pattern present in the upper plenum before the pump trip. JAEA has estimated from many scaled tests using sodium and water that the momentum dominated solution ought to be close to the real behaviour. According to JAEA, this assumption is supported by the two measured temperature changes in the middle part of the plug, 3 m and 3.75 m above the vessel bottom. However, the momentum dominated solution S1 presented in Figure6 did not reproduce this temperature change, whereas the buoyancy dominated solution S2 does roughly.Various sensitivity calculations have been performed to assure that the presence of the two solutions is not related to incomplete convergence on the mesh or on the numerical scheme. The effects of the sensitivity studies are summarized in Table5.Table 5
Influence of varied parameters on the solution.
Objective
Initialisation
Varied parameter
Solution
Convergence of solution S1
Reposing flow;T=418°C
Coarse mesh; 1st order convection scheme
S1c
Solution S1c;
Intermediate mesh; 2nd order convection scheme
S1i
Solution S1i;
Fine mesh; 2nd order convection scheme
S1f
Convergence of solution S2
Reposing flow;T=418°C
Intermediate mesh; 2nd order convection scheme
S2i
Solution S2i;
Fine mesh; 2nd order convection scheme
S2f
Solution S2f;
Coarse mesh; 1st order convection scheme
S2c
Reposing flow;T=513°C
Intermediate mesh; 1st order convection scheme
S2i
Effect of UCS
Solution S2i;  without UCS
Intermediate mesh; 2nd order convection scheme
S3iThe following conclusions of the sensitivity calculations can be made.(i)
A bifurcation into momentum (S1) and buoyancy (S2) dominated solutions is observed. Each solution has been verified to be meshing independent.(ii)
The momentum dominated solution seems to be attained by using initially coarse meshes and low-order convection scheme (initializing with a reposing fluid).(iii)
The buoyancy dominated solution seems to be attained by using initially finer meshes and higher order convection schemes (initializing with a reposing fluid).(iv)
Once a steady-state solution is achieved (either S1 or S2), this solution cannot be altered by changing the numerical scheme and/or the mesh refinement (all solutions were converged on the fine mesh and 2nd order numerical scheme).(v)
The modelling of the UCS can have a strong influence on the flow pattern (solution S3).A dependency of the flow pattern on the experimental management has been observed at CEA in sodium mixed convection experiments. In a small scale experiment comparable to the MONJU pump trip experiment, a momentum dominated flow similar to the solution S1 has been observed for mixed convection conditions after reducing the inlet flow velocity from an initially higher value to the target value. A buoyancy dominated flow similar to S2 has been observed for identical mixed convection conditions after increasing the flow velocity from an initially lower value to the target value. It cannot be excluded that the flow pattern in the MONJU plenum before the pump trip is as sensitive to the plant management before experiment as the CEA experiment was to the initial conditions. Today, almost 20 years after the experiment, it seems hardly possible to know from numerical simulations the flow pattern at the beginning of the experiment without any doubt.


## 4.1. Momentum Dominated Solution

A three-step procedure was applied to ensure convergence of the solution on both the meshing and the numerical scheme. In this procedure, solutions converged on coarser meshes were interpolated on finer meshes to define the new initial condition for the subsequent calculation on the finer mesh.(i)
Initializing with a reposing (u→=0) isothermal flow at T=418°C, a preliminary calculation was performed on the coarse mesh using a 1st order convection scheme. This calculation led to solution S1c.(ii)
Initializing with the solutionS1c, a reference calculation was performed on the intermediate mesh using a 2nd order convection scheme. This calculation led to solution S1i.(iii)
Initializing with the solutionS1i, a convergence test calculation was performed on the fine mesh using a 2nd order convection scheme. This calculation led to solution S1f.Identical solutionsS1i and S1f indicate a mesh independency of the solution which is based on a 2nd order numerical scheme.In Figure5(a), both the temperature distribution and the flow field are given for the symmetry plane going through the outlet nozzle. The solution S1i is shown. The reference height 0 m corresponds to the elevation of the support plate. The momentum of the flow leaving the upper core structures is sufficiently high to allow the jet to traverse the upper plenum, to impinge on the outer barrel, and to cross over the outer barrel to finally reach the outlet nozzle. We call this solution “momentum dominated solution.”(a) Momentum dominated solution (S1i). (b) Buoyancy dominated solution (S2i). (c) Mixed convection solution (S3i).
(a)
(b)
(c)


## 4.2. Buoyancy Dominated Solution

A two-step procedure was applied to ensure convergence of the solution on the meshing.(i)
Initializing with a reposing (u→=0) isothermal flow at T=418°C, the reference calculation was performed on the intermediate mesh using a 2nd order convection scheme. This calculation led to solution S2i.(ii)
Initializing with the solutionS2i, a convergence test calculation was performed on the fine mesh using a 2nd order convection scheme. This calculation led to solution S2f.Identical solutionsS2i and S2f indicate a mesh independency of the solution which is based on a 2nd order numerical scheme.In Figure5(b), both the temperature distribution and the flow field are given for the symmetry plane going through the outlet nozzle. The solution S2i is shown. The reference height 0 m corresponds to the elevation of the support plate. The buoyancy of the flow leaving the upper core structures prevents the jet from traversing the upper plenum as in solution S1. The flow leaving the upper core structure rises along the upper instrument structure and crosses over the outer barrel to reach the outlet nozzle. We call this solution “buoyancy dominated solution.”It has been tested if it is possible to change the solution from the buoyancy dominated flow to the momentum driven flow by changing the numerical scheme. Initializing with the solutionS2f, a test calculation was performed on the coarse mesh using a 1st order convection scheme. This calculation maintained the solution S2 after convergence.The coldest sodium leaving the core region has a temperature of 418°C. When initialising the upper plenum with this temperature, buoyancy forces act instantly on the hotter jet. These buoyancy forces are minimized by initialising the upper plenum at the temperature of the hottest sodium (513°C) leaving the core region. Starting from a reposing (u→=0) isothermal flow at T=513°C, an additional calculation was performed on the intermediate mesh using a 2nd order convection scheme. This calculation also led to solution S2i.


## 4.3. Mixed Convection Solution

The possible influence of the upper core structure (UCS) on both the thermal stratification and the flow pattern within the upper plenum is not a priori evident. To show the possible effect of the UCS on the flow, an additional sensitivity calculation without UCS has been performed.Initializing with the solutionS2i (with UCS), a test calculation without UCS was performed on an intermediate mesh using a 2nd order convection scheme. This calculation led to solution S3i.Convergence tests on further mesh refinements have not been made. The solutionS3i is shown in Figure 5(c). The temperature distribution and the flow field are given for the symmetry plane going through the outlet nozzle. The reference height 0 m corresponds to the elevation of the support plate. The buoyancy of the flow leaving the upper core structures is sufficiently high to prevent the jet from traversing the upper plenum as in solution S1. However, the inertia is high enough to also prevent the jet from simply rising along the Upper Instrument Structure as in solution S2. After spreading out in the upper plenum, the flow crosses over the outer barrel and reaches the outlet nozzle. We call this solution “mixed convection solution.”


## 4.4. Lessons Learned from the Steady-State Calculations

Using different modelling options, three different steady solutions have been achieved which can be distinguished regarding the driving physical process:(i)
a momentum dominated solution S1;(ii)
a buoyancy dominated solution S2;(iii)
a mixed convection solution S3.The calculated temperature profiles along the TCP are compared in Figure6 to the experimental values. The reference height 0 m corresponds to the elevation of the support plate.Figure 6
Comparison of the axial temperature profiles.From these temperature profiles it seems not possible to identify the flow pattern present in the upper plenum before the pump trip. JAEA has estimated from many scaled tests using sodium and water that the momentum dominated solution ought to be close to the real behaviour. According to JAEA, this assumption is supported by the two measured temperature changes in the middle part of the plug, 3 m and 3.75 m above the vessel bottom. However, the momentum dominated solution S1 presented in Figure6 did not reproduce this temperature change, whereas the buoyancy dominated solution S2 does roughly.Various sensitivity calculations have been performed to assure that the presence of the two solutions is not related to incomplete convergence on the mesh or on the numerical scheme. The effects of the sensitivity studies are summarized in Table5.Table 5
Influence of varied parameters on the solution.
Objective
Initialisation
Varied parameter
Solution
Convergence of solution S1
Reposing flow;T=418°C
Coarse mesh; 1st order convection scheme
S1c
Solution S1c;
Intermediate mesh; 2nd order convection scheme
S1i
Solution S1i;
Fine mesh; 2nd order convection scheme
S1f
Convergence of solution S2
Reposing flow;T=418°C
Intermediate mesh; 2nd order convection scheme
S2i
Solution S2i;
Fine mesh; 2nd order convection scheme
S2f
Solution S2f;
Coarse mesh; 1st order convection scheme
S2c
Reposing flow;T=513°C
Intermediate mesh; 1st order convection scheme
S2i
Effect of UCS
Solution S2i;  without UCS
Intermediate mesh; 2nd order convection scheme
S3iThe following conclusions of the sensitivity calculations can be made.(i)
A bifurcation into momentum (S1) and buoyancy (S2) dominated solutions is observed. Each solution has been verified to be meshing independent.(ii)
The momentum dominated solution seems to be attained by using initially coarse meshes and low-order convection scheme (initializing with a reposing fluid).(iii)
The buoyancy dominated solution seems to be attained by using initially finer meshes and higher order convection schemes (initializing with a reposing fluid).(iv)
Once a steady-state solution is achieved (either S1 or S2), this solution cannot be altered by changing the numerical scheme and/or the mesh refinement (all solutions were converged on the fine mesh and 2nd order numerical scheme).(v)
The modelling of the UCS can have a strong influence on the flow pattern (solution S3).A dependency of the flow pattern on the experimental management has been observed at CEA in sodium mixed convection experiments. In a small scale experiment comparable to the MONJU pump trip experiment, a momentum dominated flow similar to the solution S1 has been observed for mixed convection conditions after reducing the inlet flow velocity from an initially higher value to the target value. A buoyancy dominated flow similar to S2 has been observed for identical mixed convection conditions after increasing the flow velocity from an initially lower value to the target value. It cannot be excluded that the flow pattern in the MONJU plenum before the pump trip is as sensitive to the plant management before experiment as the CEA experiment was to the initial conditions. Today, almost 20 years after the experiment, it seems hardly possible to know from numerical simulations the flow pattern at the beginning of the experiment without any doubt.


## 5. Analysis of the Pump Trip Experiment

The intermediate mesh, 2nd order discretization schemes, and the 1st order explicit Euler time marching scheme are applied to analyse the transient of the pump trip experiment. The Courant-Friedrich-Levy stability condition is respected during the whole transient (CFL = 1). The calculation is initialized with solutionS2i. This calculation serves as reference calculation.The thermal stratification in the symmetry plane, 2 minutes and 10 minutes after the pump trip, is given in Figures7(a) and 7(b), respectively. The flow leaving the reactor core is colder than the sodium in the upper plenum and is thus accumulating in the lower part of the plenum. In accordance with the experiment, the Trio_U analysis of the pump trip transient shows a fairly well prediction of the thermal stratification formation in the first 10 minutes of the pump trip. However, about 15 minutes after the pump trip the stratification disappears and a homogeneous temperature distribution is calculated for the upper plenum. In the experiment, a homogeneous temperature distribution was measured but only after about 120 minutes (see Figure 3).(a) The flow field in the symmetry plane att=2 min. (b) The flow field in the symmetry plane at t=10 min.
(a)
(b)During the 3rd CRP meeting, an underestimation of the flow passing through the flow holes (FH) of the inner barrel has been identified by the benchmark participants as a potential source of uncertainty regarding long-term stratification formation. The exact form of the holes is not known. JAEA has provided circular holes with sharp edges in the CAD file [10]. The pressure drop in a sharp edge circular orifices can be more than 50% higher compared to that of a circular orifice of the same cross section but with smoothed edges [13]. It is important to note that the steady-state solution is not affected by this geometrical uncertainty. The pressure drop in the flow holes is so high that the main flow crosses above the inner barrel whatever the shape of the FH is (see Figure 5(a) to Figure 5(c)).The reference transient calculation with FH with sharp edges has been repeated with FH with rounded edges. During the transient, the overall flow patterns of sharp edge and round edge calculations look somewhat different [18]. Nevertheless, the resulting thermal stratifications seem not to be influenced significantly by the flow field. This can be seen from the axial temperature profiles along the thermocouple plug which are shown in Figures 8(a) and 8(b). The Figures show the vertical thermal stratification 2 minutes and 10 minutes after the pump trip. The effect of the shape of the flow holes on the calculated temperature profiles is clearly visible. However, the calculated stratification front moves upward faster than that measured in the experiment.(a) Axial temperature profiles att=2 min. (b) Axial temperature profiles at t=10 min.
(a)
(b)Sofu [19] has also considered inner barrel flow holes with perfectly round edges. The numerical results obtained with this geometrical model were found to be in excellent agreement with the experimental data, confirming the importance of the flow through the FH on the thermal stratification patterns in the MONJU upper plenum. Separate effect studies on the meshing of inner barrel flow holes with sharp edges and round edges are underway for flow conditions which are representative of the pump trip. Special attention is taken to get a solution for the MONJU pump trip which is independent on the mesh refinement within the holes.Later then about 10 minutes after the pump trip, an axially linear temperature distribution is present in the upper plenum at elevations higher than about 2 m above the support plate (see Figure3). The used pressure discretization (P0/P1) cannot balance a constant temperature gradient by a parabolic pressure course without creating spurious velocity modes. Thus, it is possible that the overestimated mixing in the upper part of the upper plenum later then, about 10 minutes after the pump trip, is related to the presence of such spurious velocity modes. More detailed investigations based on an improved discretisation of the pressure are under way.


## 6. Conclusion

In the framework of an IAEA CRP, benchmark analyses of sodium convection in the upper plenum of the MONJU reactor vessel have been performed at CEA by using the CFD code Trio_U. In angular direction only 1/6 of the MONJU upper plenum has been taken into account in the modelling. This reduced domain was discretized in up to 3.3 million tetrahedral elements. A high Reynolds number turbulence model (k-ε) was used to account for turbulent mixing with additional terms for treating buoyancy effects.For the steady-state initial condition before the pump trip transient, a bifurcation of the solutions into a “momentum dominated solution” and a “buoyancy dominated solution” was observed. Such a bifurcation has also been found experimentally at CEA in mixed convection experiments. Initializing the coarsest mesh with a reposing fluid and using a 1st order convection scheme, the momentum dominated solution was attained. Initializing the intermediate mesh with a reposing fluid and using a 2nd order convection scheme, the buoyancy dominated solution was attained. After successively refining the mesh and using 2nd order numerical schemes, each of the two solutions can be considered as meshing independent, because each solution was identical on the two finest grids. The experimentally found temperature stratification seems not to be a good indicator to distinguish the flow fields of the two solutions. For a more profound understanding, further numerical analysis is necessary. This concerns the solution procedure (iterative methods as SIMPLE instead of a transient calculation as presented) and may be the use of more sophisticated turbulence modelling approaches.The pump trip experiment has shown the formation of a thermal stratification within the plenum. The calculations with the intermediate mesh predict well the stratification formation in the first 10 minutes of the pump trip. In the experiment, the stratification persists for more than two hours, whereas all calculations of Trio_U have shown a homogenization of the temperature in the plenum after about 15 to 20 minutes (significant overestimation of the internal mixing). Further analysis is under way to better understand the overestimation of the mixing in the upper plenum. Special attention will be turned on balancing in reposing fluids linear temperature gradients by the pressure.



---

*Source: 987197-2013-05-30.xml*