Simulation of Natural Convection of Helium in DEMO Cryostat Using OpenFOAM

(1)

705.1

Simulation of Natural Convection of Helium in DEMO Cryostat Using OpenFOAM

Rok Krpan, Matej Tekavčič, Boštjan Končar Jožef Stefan Institute

Reactor Engineering Division Jamova cesta 39 SI-1000, Ljubljana, Slovenia

rok.krpan@ijs.si, matej.tekavcic@ijs.si, bostjan.koncar@ijs.si ABSTRACT

The accident scenario considered in this study considers that the cryostat in DEMO fusion reactor is filled with helium gas. This enables the heat transfer with natural convection between the cryostat thermal shield, which remains actively cooled, and cryostat, which is on the outside in contact with the atmosphere at room temperature. This can cause a significant cool-down of the cryostat walls. The open source computational fluid dynamics (CFD) code OpenFOAM is used to simulate natural convection of helium gas and to assess the temperature distribution inside the compartment between the cryostat thermal shield and the cryostat, and the temperature of the inner cryostat wall. The results are also compared to previous results obtained with the CFD code ANSYS CFX.

1 INTRODUCTION

One of the hypothetical accident scenarios in the DEMO fusion reactor considers ingress of helium into the cryostat due to the break of the cryogenic cooling line, which is used to cool the superconducting magnets. The accident scenario used in this study assumes that the cryostat thermal shield (CTS) remains actively cooled (e.g. approximately at a constant temperature) and the cryostat is on the outside surrounded by air at room temperature [1]. For such scenario, the majority of the heat transfer occurs in the interspace between the CTS and the cryostat inner wall. The temperature difference between the CTS and the cryostat generates natural convection of helium in the interspace, which causes a significant local cool-down of the cryostat walls.

The cryostat structures in DEMO are directly connected to the bio-shield and the cryostat contraction is limited by its supports. The heat transfer and the temperature distributions on the cryostat walls are thus needed to design the supports, which must be able to withstand high temperature differences and forces caused by the cryostat contractions.

The steady-state simulations of the natural convection was performed using the open source computational fluid dynamics (CFD) code OpenFOAM [2]. A numerical model including a 20° sector of the isolated interspace between the CTS and the cryostat inner wall was developed. The aim of this study is to determine the temperature distribution on the cryostat walls and to compare the results obtained with the open-source CFD code to similar work performed with commercial CFD code ANSYS CFX 18.1 [3] by Tekavčič et al. [4].

(2)

Proceedings of the International Conference Nuclear Energy for New Europe, Portorož, Slovenia, September 9̶ 12, 2019

2 SIMULATION MODEL

Numerical simulations were performed with the open-source CFD code OpenFOAM, version v1606+ [2]. To simulate the steady state helium natural convection, the Reynolds Averaged Navier-Stokes (RANS) approach was used with the k-ω Shear Stress Transport (SST) turbulence model [5]. Root Mean Squared (RMS) residual values of less than 10^-7 were set as the convergence criteria. Because second order numerical schemes caused numerical instability and did not converge properly, first order numerical schemes were used.

2.1 Physical model

To simulate the natural convection of He gas in the isolated compartment between the actively cooled CTS and the cryostat wall, the buoyantSimpleFoam solver was used. This is a steady-state solver intended for buoyant, turbulent flows of compressible fluids. The momentum equation was solved:

𝜕𝜕

𝜕𝜕𝑥𝑥𝑗𝑗�𝜌𝜌𝑈𝑈𝑖𝑖𝑈𝑈𝑗𝑗�= − 𝜕𝜕𝜕𝜕

𝜕𝜕𝑥𝑥𝑖𝑖 + 𝜕𝜕

𝜕𝜕𝑥𝑥𝑗𝑗�(𝜇𝜇+𝜇𝜇𝑡𝑡)�𝜕𝜕𝑈𝑈_𝑖𝑖

𝜕𝜕𝑥𝑥𝑗𝑗 +𝜕𝜕𝑈𝑈𝑗𝑗

𝜕𝜕𝑥𝑥𝑖𝑖�� − 𝑔𝑔𝑘𝑘𝑥𝑥𝑘𝑘 𝜕𝜕𝜌𝜌

𝜕𝜕𝑥𝑥𝑖𝑖 (1)

where ρ, U, p, μ, μt and g are density, velocity, pressure, dynamic and eddy viscosity, and gravitational acceleration, respectively. The transport of energy was solved using the enthalpy (h) equation:

𝜕𝜕

𝜕𝜕𝑥𝑥𝑗𝑗(𝜌𝜌𝑈𝑈_𝑖𝑖ℎ) = 𝜕𝜕

𝜕𝜕𝑥𝑥𝑗𝑗��𝜇𝜇 𝑃𝑃𝑃𝑃+ 𝜇𝜇_𝑡𝑡

𝑃𝑃𝑃𝑃𝑡𝑡� 𝜕𝜕ℎ

𝜕𝜕𝑥𝑥𝑗𝑗�+𝜌𝜌𝑔𝑔𝑘𝑘𝑈𝑈𝑘𝑘 (2)

where Pr and Prt are Prandtl number and turbulent Prandtl number, respectively. In buoyantSimpleFoam solver the default value of turbulent Prandtl is 0.85 and cannot be changed.

2.2 Computational mesh and boundary conditions

Tokamak and cryostat geometry used in this work are based on the 2015 DEMO baseline design with 18 toroidal sectors [6]. Figure 1 left shows a 20° axial symmetric wedge. The yellow highlighted part represents the closed compartment between CTS and cryostat wall, which is filled with He gas during the accident. Due to the symmetry it was cut in half and considered as a computational domain (Figure 1 right).

Numerical meshes used in this study were taken from Tekavčič et. al [4], where they were generated with ANSYS Workbench [3]. This mesh is composed of a combination of tetrahedral elements in the centre of the domain and prism elements in the wall boundary layer on CTS and cryostat walls. Another mesh was generated with the snappyHexMesh utility, which is included in OpenFOAM programme package. This second mesh consists entirely of hexahedral elements. In Table 1 the mesh parameters used in this study are listed. Figure 2 shows the results of the mesh convergence study, namely horizontal temperature profiles (left) and cryostat inner wall vertical temperature profiles (right) at p =10² Pa and heat transfer coefficient hw = 10 W/m²K. At the highest measuring position, the horizontal temperature profiles are the same, regardless of the mesh used, while at lower positions some discrepancies can be observed. On the other hand, there are some discrepancies in the cryostat inner wall temperature vertical profiles and the results do not converge towards the finest mesh. Because of the numerical artefacts generated by the meshes and the complex geometry, at some point there is smaller difference between the results given by the coarsest and the finest mesh (mesh 1 and 4), than by finer and the finest (mesh 3 and 4). Since there were no major discrepancies between the

(3)

results, and to reduce computational time, the mesh 1 was used for calculations performed with both CFD codes.

Figure 1: Reference tokamak geometry (left) and considered computational domain with boundary patches (right) [4][6].

Figure 2: Horizontal temperature profiles (left) and cryostat inner wall vertical temperature profiles (right) obtained with different meshes at p =10² Pa and hw = 10 W/m²K.

Figure 1 right shows prescribed boundaries. A constant temperature of 80 K was prescribed on CTS wall and adiabatic boundary condition was prescribed on ports. On cryostat

1

2

3

4 snappyHexMesh

80 120 160 200 240 280 320

14.5 15 15.5 16 16.5 17 17.5 18 18.5 19 19.5

T [K]

Z = 10 m

80 120 160 200 240 280 320

17.6 17.8 18 18.2 18.4 18.6 18.8 19 19.2 19.4

T [K]

Z = 0 m

80 120 160 200 240 280 320

15 15.5 16 16.5 17 17.5 18 18.5 19 19.5

T [K]

X [m]

Z = -10 m

-10 -5 0 5 10

288 290 292 294 296

Z [m]

T [K]

Ports

CTS Wall

Cryostat

Symmetry planes

(4)

walls, heat transfer coefficient (hw) and outside air temperature (Tair = 293 K) were set. On all boundaries (except the symmetry planes), wall functions were prescribed for turbulence quantities and the no-slip boundary condition was prescribed for velocity.

Table 1: Mesh generation software, total number of mesh cells, maximal cell size, number of prism layer, and maximal calculated non-dimensional first cell to wall distance.

Software Mesh N Max. size [m] Prism layers # y+

ANSYS Workbench

1 78,233 0.5

11

0.22

2 148,161 0.3 0.16

3 380,737 0.2 0.14

4 2,856,790 0.1 0.13

snappyHexMesh 297,857 0.5 0.49

Several cases with different initial and boundary conditions were considered. This conditions were already assessed and discussed in Tekavčič et. al [4]. Helium reference pressure values inside the computational domain were 10² Pa, 10³ Pa, 10⁴ Pa and 10⁵ Pa. Heat transfer coefficients prescribed on cryostat walls were 1 W/m²K and 10 W/m²K. The lowest value is conservative and represents excessively poor heat transfer from the outer air.

3 RESULTS AND DISCUSSION

Horizontal profiles are extracted from simulations at positions showed in Figure 1 right.

Red lines are located at z = 10 m, z = 0 m and at z = -10 m.

Figure 3 shows volume-average temperature inside computational domain (top), the surface-average temperature of the cryostat wall (middle) and the lowest temperature of the cryostat wall (bottom) for different pressures and heat transfer coefficients. It can be observed that all temperatures decrease with increasing p and decreasing hw. Higher p means higher mass of helium released in the compartment, which increases the ability to transfer the heat. Lower hw between the cryostat outer wall and the surrounding ambient air decreases the heat influx into the compartment. The highest temperature differences from the normal operation thus appear at the highest p (p = 10⁵ Pa) and the lowest hw (hw = 1 W/m²K), where the cryostat inner wall temperature locally drops to 110.2 K, which is approximately 180 K lower than during normal operation. Figure 4 shows He temperature distribution inside the computational domain and at cryostat wall at p = 10⁵ Pa and hw = 1 W/m²K (left) and hw = 10 W/m²K (right). It can be observed that warmer He gas accumulates on the top of the domain. Similar is with cryostat inner wall temperature where colder regions are located below the half height. Besides this the warmer regions are located also above the ports, where the natural convection is obstructed.

Figure 5 shows cryostat inner wall temperature profiles along the selected line on the vertical wall (red line on the left side of Figure 5) for different p and hw. As shown on the graphs (right side of Figure 5), the temperature slightly increases with height. The temperature values increase greatly at higher reference pressure p or hw. Further, in Figure 5 the results obtained with OpenFOAM (dashed line) are also compared to the previous results obtained by the commercial code ANSYS CFX (full line) [4]. Some smaller discrepancies can be observed between the both results, mostly due to differences in turbulence model; the k-ω SST turbulence model implemented in OpenFOAM does not have additional terms for turbulence production and dissipation due to buoyancy, which are by default included in ANSYS CFX. The highest discrepancies (but in absolute terms still small) between the results can be observed at the highest p (p = 10⁵ Pa) and at the highest hw (hw = 10 W/m²K). Figure 6 shows code-to-code comparison of turbulence kinetic energy (left) and vertical velocity component (right) along the

(5)

three horizontal lines (shown on the left side of Figure 5) at p = 10⁵ Pa and hw = 10 W/m²K.

The differences in modelled turbulence kinetic energy are significant especially at the highest position. Some minor discrepancies can be noted also in the vertical velocity values (Figure 6 right).

Figure 3: Volume-average temperature in the compartment (top), surface-average cryostat wall temperature (middle) and lowest local cryostat temperature (bottom) for different initial

and boundary conditions.

Figure 4: Temperature distribution inside compartment at p = 10⁵ Pa and hw = 1 W/m²K (left) and hw = 10 W/m²K (right).

h w = 1 W/m ²K h w = 10 W/m ²K

80

120 160 200

100 1000 10000 100000

T [K]

Average T in compartment

120 160 200 240 280 320

100 1000 10000 100000

T [K]

Average cryostat wall T

80 120 160 200 240 280 320

100 1000 10000 100000

T [K]

p [Pa]

Lowest cryostat wall T

hw = 1 W/m²K hw = 10 W/m²K

(6)

Figure 5: Cryostat inner wall vertical temperature profiles for different p and different hw.

Figure 6: Turbulence kinetic energy (left) and vertical velocity horizontal profiles (right) at p

= 10⁵ Pa and hw = 10 W/m²K

Figure 7 and Figure 8 show horizontal temperature profiles inside the compartment obtained with both CFD codes at different p and hw. It can be observed that the results obtained by the OpenFOAM match the ANSYS CFX results relatively well. At the lowest p a smaller temperature gradient can be observed near the walls, with higher temperature increase in the middle of the compartment, while for higher p only a minor temperature increase can be observed in the middle but higher temperature gradient occurs near the walls. At the highest measuring position the temperature profile towards the warmer wall slightly decreases, which happens due to the geometry of the compartment. The strait near the ports obstructs the flow

ANSYS CFX 10²Pa 10³Pa

OpenFOAM 10⁴Pa 10⁵Pa

-10 -5 0 5 10

100 150 200 250 300

Z [m]

T [K]

h w= 1 W/m ²K

-10 -5 0 5 10

200 225 250 275 300

Z [m]

T [K]

h w= 10 W/m ²K

ANSYS CFX OpenFOAM

0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16

14.5 15 15.5 16 16.5 17 17.5 18 18.5 19 19.5

TKE [m2/s2 ]

Z = 10 m

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3

14.5 15 15.5 16 16.5 17 17.5 18 18.5 19 19.5

Uz [m/s]

Z = 10 m

0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.24

17.6 17.8 18 18.2 18.4 18.6 18.8 19 19.2 19.4

TKE [m2/s2 ]

Z = 0 m

-1.5 -1 -0.5 0 0.5 1 1.5 2 2.5

17.6 17.8 18 18.2 18.4 18.6 18.8 19 19.2 19.4

Uz [m/s]

Z = 0 m

0 0.02 0.04 0.06 0.08 0.1 0.12 0.14

15 15.5 16 16.5 17 17.5 18 18.5 19 19.5

TKE [m2/s2 ]

X [m]

Z = -10 m

-3 -2.5-2 -1.5-1 -0.5 0 0.5 1 1.5 2 2.5

15 15.5 16 16.5 17 17.5 18 18.5 19 19.5

Uz [m/s]

X [m]

Z = -10 m

(7)

from lower regions so strongly, that a new recirculation zone is generated above the port, where the compartment widens again.

Figure 7: Horizontal temperature profiles at hw = 1 W/m²K.

Figure 8: Horizontal temperature profiles at hw = 10 W/m²K.

80 120 160 200 240 280

14.5 15 15.5 16 16.5 17 17.5 18 18.5 19 19.5

T [K]

Z = 10 m

80 100 120 140 160 180 200

14.5 15 15.5 16 16.5 17 17.5 18 18.5 19 19.5 Z = 10 m

80 120 160 200 240 280

17.6 18 18.4 18.8 19.2

T [K]

Z = 0 m

80 100 120 140 160 180

17.6 18 18.4 18.8 19.2 Z = 0 m

80 120 160 200 240 280

15 15.5 16 16.5 17 17.5 18 18.5 19 19.5

T [K]

X [m]

Z = -10 m

80 100 120 140 160 180

15 15.5 16 16.5 17 17.5 18 18.5 19 19.5 X [m]

Z = -10 m

80 120 160 200 240 280 320

14.5 15 15.5 16 16.5 17 17.5 18 18.5 19 19.5

T [K]

Z = 10 m

80 120 160 200 240 280

14.5 15 15.5 16 16.5 17 17.5 18 18.5 19 19.5 Z = 10 m

80 120 160 200 240 280 320

17.6 18 18.4 18.8 19.2

T [K]

Z = 0 m

80 120 160 200 240 280

17.6 18 18.4 18.8 19.2 Z = 0 m

80 120 160 200 240 280 320

15 15.5 16 16.5 17 17.5 18 18.5 19 19.5

T [K]

X [m]

Z = -10 m

80 120 160 200 240 280

15 15.5 16 16.5 17 17.5 18 18.5 19 19.5 X [m]

Z = -10 m

hw = 1 W/m²K

hw = 10 W/m²K

ANSYS CFX OpenFOAM 10² Pa 10³ Pa 10⁴ Pa 10⁵ Pa

(8)

4 CONCLUSION

During a cryogenic cooling line rupture, the cryostat is filled with helium gas and the conditions for heat transfer due to natural convection are established in the interspace between the CTS and the cryostat inner wall. The heat transfer between the CTS, which in the scenario used for this study remains actively cooled at 80 K, and the cryostat causes a significant local cool-down and a non-homogeneous temperature distribution of the cryostat walls.

The results obtained with OpenFOAM CFD code were compared to ANSYS CFX 18.1 CFD code and the discrepancies in the results mainly due to the differences in the implemented turbulence models were not significant.

In the worst scenario considered in this study is the case with the highest possible helium pressure inside the cryostat (10⁵ Pa) and the lowest heat transfer coefficient between cryostat wall and outer air (1 W/m²K). In this case the local cryostat temperature decreases down to 110 K, which is 180 K colder than during the normal operation. Otherwise, for pressures up to 10³ Pa and heat transfer coefficient of 10 W/m²K a much smaller temperature drop (by approximately 30 K) is observed. To minimize the cooling of the cryostat wall during the helium ingress accident, the heat transfer between cryostat outer wall and the outside should be maximized.

ACKNOWLEDGMENTS

The present work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014 - 2018 and 2019-2020 under grant agreement No. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. The financial support from the Slovenian Research Agency (research project J2-9209 and core research program P2-0405) is also gratefully acknowledged.

REFERENCES

[1] C. Bachmann, Task specifications - CFX assessment of cryostat air/helium ingress, PMI- 2.1, May 2018.

[2] OpenCFD Ltd., “OpenFOAM: The open source CFD toolbox”,

http://www.openfoam.com/, view: 16.08.2016

[3] ANSYS Inc., ANSYS ® Academic Research, Release 18.1, Help System, Solver Theory, Multiphase Flow Theory, 2017.

[4] M. Tekavčič, B. Končar, M. Draksler, C. Bachmann, “Simulation of the Incident Helium Ingress into the Cryostat of the DEMO Fusion Reactor”, In: I. Jenčič (Ed.), Proceedings of 27th International Conference Nuclear Energy for New Europe – NENE 2018, Portorož, Slovenia, September 10-13, Nuclear Society of Slovenia, 2018, pp. 616-1- 616- 10.

[5] F. R. Menter, T. Esch, "Elements of Industrial Heat Transfer Prediction”, 16th Brazilian Congress of Mechanical Engineering (COBEM), 2001.

[6] Baseline CAD Configuration Model: 2D3FBF, May 2015.

[7] B. E. Launder, D. B. Spalding, “The Numerical Computation of Turbulent Flows”, Computer Methods in Applied Mechanics and Engineering, 3, 1974, pp. 269-289.