# Thermal Analysis for the Dense Granular Target of CIADS

Authors: Kang Chen; Yongwei Yang; Yucui Gao
Journal: Science and Technology of Nuclear Installations
 (2016)

---

## Abstract

For the China Initiative Accelerator Driven System (CIADS), the energy of the protons is 250 MeV, and the current intensity will reach 10 milliamperes. A new concept of a dense granular spallation target is proposed for which the tungsten granules are chosen as the target material. After being bombarded with the accelerated protons from the accelerator, the tungsten granules with high-temperature flow out of the subcritical reactor and the heat is removed by the heat exchanger. One key issue of the target is to remove the 2.5 MW heat deposition safely. Another one is the heat exchange between the target and the subcritical reactor. Based on the model of effective thermal conductivity, a new thermal code is developed in Matlab. The new code is used to calculate the temperature field of the target area near active zone and it is partly verified by commercial CFD code Fluent. The result shows that the peak temperature of the target zone is nearly 740°C and the reactor and the target are proved to be uncoupled in thermal process.


---

## Body

## 1. Introduction

China Initiative Accelerator Driven System (CIADS) is a project proposed by Chinese Academy of Science (CAS) [1]. It consists of three subsystems: accelerator, spallation target, and subcritical reactor (Figure 1). Among them, the spallation target is responsible for generating neutrons, which are used to drive the subcritical reactor [2]. For CIADS, a new concept of spallation target is proposed which is called dense granular target [3].Figure 1
Schematic of CIADS.In present design, the target material is a large collection of discrete tungsten (W) granules. The container for the granules is a hopper and W alloy grains flow into the spallation region under gravity from the upper annular duct where the beam pipe is located. Then a proton beam from the accelerator will interact with the flowing W grains below the pipe. The grains will pass through the spallation region quickly and discharge from the orifice of the hopper to avoid being melted down. What is more, a circuit, including a heat exchanger, grain filter, grain elevator, grain storage, and cover gas (helium) system, will guarantee the normal circulation of the working medium in the entire system.The proton beam from the accelerator with an energy of 250 MeV and a current of 10 mA reacts with the tungsten nuclei, and approximately 2.5 MW of heat is deposited in the target zone, which must be removed by circulating the granules. As the granules flow through the center channel of the subcritical reactor, if the target and the reactor are decoupled in heat transfer, the safety review for the spallation target can be simplified. Therefore, it is necessary to know the temperature field of the target. Furthermore, as the target and the reactor are mainly coupled by neutron near the active zone and the heat is mainly deposited there, we focus our attention on the target area near the active zone in present study (Figure2).Figure 2
Schematic of the granular target.Though granular material is one of the most common forms in existence, there has been no comprehensive theory that can reliably predict the behavior of such materials in technical devices [4]. The discrete element method (DEM) is now a popular method of simulating granular material but it needs a great amount of calculation to simulate the target and it can hardly calculate the temperature field. In this paper, a new thermal code based on the model of effective thermal conductivity is developed in Matlab which places smaller demands on calculation capability. The commercial CFD code Fluent is used to verify the new code. The result calculated from the new code and Fluent is nearly the same. Moreover, the result shows that the heat exchange between the target and the reactor can be neglected.


## 2. Computational Model

Concentrating on the target area near the active zone, the model is simplified to a cylinder with a diameter of 28 cm and a height of 80 cm (Figure3). The tungsten granules flow into the spallation region with a velocity of 0.5 m/s and a temperature of 250°C. A cylinder with a radius of 14 cm and a height of 5 mm is chosen as the control volume. By tracking it, the flow of the granules can be simulated. As to the control volume, it can be seen as stagnant helium-saturated porous beds of spheres and the heat transfer in the radial direction can be calculated based on the model of effective thermal conductivity. The heat transfer in the axial direction is taken into consideration through source correction.Figure 3
Schematic of the whole computational model.


### 2.1. Governing Equations

Because of the symmetric condition, a two-dimensional model is established. The height of the control volume is set as 5 mm (Figure4) and the time step size is set as 0.01 s to match with the velocity of the granular flow. Meanwhile, the number of time steps is set as 160 to match with the height of 80 cm. By means of the first boundary condition, the effect of the reactor on the heat transfer is taken into consideration to realize the coupling of target and reactor. The temperature field of the whole model can be obtained by putting the temperature fields of control volume at different time together depending on the relationship between time and height.Figure 4
Computational model of the control volume.In the control volume, based on the model of effective thermal conductivity, it is simplified to be a one-dimensional transient heat conduction problem. The governing equation can be seen as follows [5]:(1)ρc∂T∂t=1r∂∂rλr∂T∂r+S,0<r≤b,T=250,t=0,0<r≤b,T=Tt,r=b,∂T∂r=0,r=0,where 250°C is the designed inlet temperature, ρ is the average density, c is the average special heat capacity, λ is the effective thermal conductivity and will be explained in Section 2.2, and S is the heat source. The explicit time integration between t and t+Δt is performed in this equation. Then, the discrete form for the equation can be obtained as follows: (2)apTpt+Δt=aeTet+ap-ae-awTpt+awTwt+b,T=250,t=0,0<r≤b,Tn-1=TY,r=b,T0=T1,r=0,where(3)ap=ρcrpΔrΔt,aw=λrwΔr,ae=λreΔr,b=rpSΔr,and T(Y) changes from 300°C to 385°C [6], which is the designed temperature of the reactor.


### 2.2. Effective Thermal Conductivity

The control volume can be seen as stagnant helium-saturated porous beds of spheres and there are three main processes of heat transfer in the control volume: the radiant heat transfer in the gaps, the conductive heat transfer between the tungsten granules, and the conductive heat transfer of the helium. The present research converts the practical heat transfer in porous media into equivalent heat conduction of common solid materials which satisfies the heat conduction equations, and thus it becomes a one-dimensional transient heat conduction problem.Effective thermal conductivity is an important physical quantity describing the heat conduction in porous media. Based on a one-dimensional heat flow model for conduction through a packed bed of spherical particles, Zehner and Schlunder presented a correlation (Z-S formula) for the stagnant effective thermal conductivity [7] as follows:(4)λegλf=1-1-ε+21-ε1-λB1-λB1-λB2ln⁡1λB-B+12-B-11-λB,where(5)λ=λfλs,B=1.251-εε10/9.The equations had been verified against experimental data obtained by Kladias and Bandyopadhaya’s experiment [8] in which various combinations of solid and fluid are used to produce the porous media. The stagnant thermal conductivities measured by them agreed reasonably well with the predictions made by Zehner and Schlunder as long as the conductivity ratio is not too high.


### 2.3. The Correction of Heat Transfer in the Axial Direction

In the previous governing equation, only the heat transfer in the radial direction is taken into consideration. To consider the axial heat transfer, the heat source is modified in every step during iteration (Figure5).Figure 5
Procedure for the correction of axial heat transfer.The following equation is used to calculate the heat source which needs to be modified (Figure6):(6)q=-λΔTΔx∗Δh,q′=-λTm-Tm-1Δx∗Δh-λTm-Tm+1Δx∗Δh,Sn=S0+q′,where ΔT is the temperature difference between the same node in successive time, Δh is the height of the control volume, Δx is the distance of the same node in successive time, q′ is the heat transfer to be modified, Sn is the heat source after being modified, and S0 is the heat transfer calculated by neutronics.Figure 6
Schematic of the calculation of the axial heat transfer.


### 2.4. Heat Deposition and Material Properties

The heat deposition is calculated by neutronics code and the proton beam has a radius of 5 cm with Gaussian distribution. Figure7 shows the energy distribution in the target zone. The maximum energy density reaches 6000 W/cm3, appearing near the inlet area.Figure 7
Heat deposition in the target zone.The material properties of tungsten and helium are listed in Table1.Table 1
The material properties of tungsten granules and helium.
Symbol Property Value d Diameter (mm) 1 ε Porosity 0.52 λ h Thermal conductivity of helium (W·m−1·K−1) 0.1426 ρ h Density of helium (Kg·m−3) 0.1785 c h Specific heat of helium (J·Kg−1·K−1) 3213.7 λ t Thermal conductivity of tungsten (W·m−1·K−1) 174 ρ t Density of tungsten (Kg·m−3) 19350 c t Specific heat of tungsten (J·Kg−1·K−1) 130The average special heat capacity can be calculated by the following equation:(7)c=Pt·ct+1-Pt·ch,where (8)Pt=mtm.Inserting ε, ρt, ct, ρh, and ch into the equation,(9)Pt=0.52∗193500.52∗19350+0.48∗0.1785≈1.c=130J·Kg-1·K-1The average density can be calculated by the following equation:(10)ρ=ε·ρt+1-ε·ρh.Insertingε, ρt, and ρh into the equation,(11)ρ=0.52∗19350+0.48∗0.1785=10062Kg·m-3.Insertingλh, λt, and ε in Z-S formula, the effective thermal conductivity of the spallation target can be calculated. The results are listed in Table 2.Table 2
Quantities calculated by Z-S formula.
Symbol Definition Value λ Conductivity ratio 8.2 e - 4 λ e Effective thermal conductivity (W·m−1·K−1) 1.176


## 2.1. Governing Equations

Because of the symmetric condition, a two-dimensional model is established. The height of the control volume is set as 5 mm (Figure4) and the time step size is set as 0.01 s to match with the velocity of the granular flow. Meanwhile, the number of time steps is set as 160 to match with the height of 80 cm. By means of the first boundary condition, the effect of the reactor on the heat transfer is taken into consideration to realize the coupling of target and reactor. The temperature field of the whole model can be obtained by putting the temperature fields of control volume at different time together depending on the relationship between time and height.Figure 4
Computational model of the control volume.In the control volume, based on the model of effective thermal conductivity, it is simplified to be a one-dimensional transient heat conduction problem. The governing equation can be seen as follows [5]:(1)ρc∂T∂t=1r∂∂rλr∂T∂r+S,0<r≤b,T=250,t=0,0<r≤b,T=Tt,r=b,∂T∂r=0,r=0,where 250°C is the designed inlet temperature, ρ is the average density, c is the average special heat capacity, λ is the effective thermal conductivity and will be explained in Section 2.2, and S is the heat source. The explicit time integration between t and t+Δt is performed in this equation. Then, the discrete form for the equation can be obtained as follows: (2)apTpt+Δt=aeTet+ap-ae-awTpt+awTwt+b,T=250,t=0,0<r≤b,Tn-1=TY,r=b,T0=T1,r=0,where(3)ap=ρcrpΔrΔt,aw=λrwΔr,ae=λreΔr,b=rpSΔr,and T(Y) changes from 300°C to 385°C [6], which is the designed temperature of the reactor.


## 2.2. Effective Thermal Conductivity

The control volume can be seen as stagnant helium-saturated porous beds of spheres and there are three main processes of heat transfer in the control volume: the radiant heat transfer in the gaps, the conductive heat transfer between the tungsten granules, and the conductive heat transfer of the helium. The present research converts the practical heat transfer in porous media into equivalent heat conduction of common solid materials which satisfies the heat conduction equations, and thus it becomes a one-dimensional transient heat conduction problem.Effective thermal conductivity is an important physical quantity describing the heat conduction in porous media. Based on a one-dimensional heat flow model for conduction through a packed bed of spherical particles, Zehner and Schlunder presented a correlation (Z-S formula) for the stagnant effective thermal conductivity [7] as follows:(4)λegλf=1-1-ε+21-ε1-λB1-λB1-λB2ln⁡1λB-B+12-B-11-λB,where(5)λ=λfλs,B=1.251-εε10/9.The equations had been verified against experimental data obtained by Kladias and Bandyopadhaya’s experiment [8] in which various combinations of solid and fluid are used to produce the porous media. The stagnant thermal conductivities measured by them agreed reasonably well with the predictions made by Zehner and Schlunder as long as the conductivity ratio is not too high.


## 2.3. The Correction of Heat Transfer in the Axial Direction

In the previous governing equation, only the heat transfer in the radial direction is taken into consideration. To consider the axial heat transfer, the heat source is modified in every step during iteration (Figure5).Figure 5
Procedure for the correction of axial heat transfer.The following equation is used to calculate the heat source which needs to be modified (Figure6):(6)q=-λΔTΔx∗Δh,q′=-λTm-Tm-1Δx∗Δh-λTm-Tm+1Δx∗Δh,Sn=S0+q′,where ΔT is the temperature difference between the same node in successive time, Δh is the height of the control volume, Δx is the distance of the same node in successive time, q′ is the heat transfer to be modified, Sn is the heat source after being modified, and S0 is the heat transfer calculated by neutronics.Figure 6
Schematic of the calculation of the axial heat transfer.


## 2.4. Heat Deposition and Material Properties

The heat deposition is calculated by neutronics code and the proton beam has a radius of 5 cm with Gaussian distribution. Figure7 shows the energy distribution in the target zone. The maximum energy density reaches 6000 W/cm3, appearing near the inlet area.Figure 7
Heat deposition in the target zone.The material properties of tungsten and helium are listed in Table1.Table 1
The material properties of tungsten granules and helium.
Symbol Property Value d Diameter (mm) 1 ε Porosity 0.52 λ h Thermal conductivity of helium (W·m−1·K−1) 0.1426 ρ h Density of helium (Kg·m−3) 0.1785 c h Specific heat of helium (J·Kg−1·K−1) 3213.7 λ t Thermal conductivity of tungsten (W·m−1·K−1) 174 ρ t Density of tungsten (Kg·m−3) 19350 c t Specific heat of tungsten (J·Kg−1·K−1) 130The average special heat capacity can be calculated by the following equation:(7)c=Pt·ct+1-Pt·ch,where (8)Pt=mtm.Inserting ε, ρt, ct, ρh, and ch into the equation,(9)Pt=0.52∗193500.52∗19350+0.48∗0.1785≈1.c=130J·Kg-1·K-1The average density can be calculated by the following equation:(10)ρ=ε·ρt+1-ε·ρh.Insertingε, ρt, and ρh into the equation,(11)ρ=0.52∗19350+0.48∗0.1785=10062Kg·m-3.Insertingλh, λt, and ε in Z-S formula, the effective thermal conductivity of the spallation target can be calculated. The results are listed in Table 2.Table 2
Quantities calculated by Z-S formula.
Symbol Definition Value λ Conductivity ratio 8.2 e - 4 λ e Effective thermal conductivity (W·m−1·K−1) 1.176


## 3. Results and Discussion

To validate the solver in the new code, the commercial CFD code Fluent is used, in which the heat transfer in the axial direction can hardly be considered. And then, taking the axial heat transfer into consideration, the new thermal code developed in Matlab is used to calculate the temperature field of the target area near active zone. The effect of the axial heat transfer is also analyzed.


### 3.1. Code Validation

Neglecting the correction of the axial heat transfer, the new code and Fluent were used to calculate the temperature field of the target because in Fluent the axial heat transfer can hardly be taken into consideration. In Fluent and the new code, the same physics model is adopted.In Fluent, the granular flow is also considered to be solid because the momentum equation for the fluid in Fluent is not suitable. The same control volume whose height is 5 mm and width is 14 cm is chosen by tracking which flow of the tungsten granules is simulated. The mesh of the control volume can be seen in Figure8. The two-dimensional transient solver is chosen. The time step size is set as 0.01 s and the number of time steps is 160 s. The initial temperature is set as 250°C. The upper and lower boundaries are set as adiabatic. The left boundary is set as the axis and the right boundary adopts the first boundary condition. The temperature of the right boundary changes from 380°C to 300°C. Moreover, the heat deposition in Section 2.4 is imported to Fluent with the help of UDF.Figure 8
Mesh of the control volume.The result shows that the temperature field calculated by these two codes is almost the same (Figure9). The temperature changes greatly near the inlet area and tends to keep unchanged in the following flow because the heat source concentrates near the inlet. Moreover, the peak temperature is 740.5°C calculated by the new code and 739.5°C calculated by Fluent. The result validates the governing equations and the discretization methods in the new code.Figure 9
Temperature field of the target area (neglecting the axial heat transfer). ((a) is calculated by the new code and (b) is calculated by Fluent.)
(a) (b)The average outlet temperature calculated by these two codes is 314°C. For the control volume, the average temperatures at the very beginning and the eighth second are 250°C and 314°C, respectively. To verify the heat source, the heat removed by the target is calculated by the following equation:(12)ΔQ=c∗m∗ΔT=c∗ρ∗π∗r2∗h∗ΔT.Inserting values into the equation,(13)ΔQ=130∗10062∗3.14∗0.142∗0.005∗314-250=25761J,P=ΔQΔt=257610.01=2.58MW.The result shows that the heat removed in the present model is 2.58 MW which is similar to the design index: 2.5 MW.


### 3.2. Temperature Field of the Target Area near the Active Zone

Considering the correction of the axial heat transfer, the new code is used to calculate temperature field of the target zone again. Figure10 shows the new temperature field. The peak temperature is also 740.5°C and the average outlet temperature is also 314°C. Table 3 shows the iterative process of the axial heat transfer correction. By 4 times of iteration, the results satisfy the demand. The results reveal that the axial heat transfer has a little effect on the temperature field. It is easy to understand because the temperature difference between successive time is small, as well as the effective thermal conductivity.Table 3
Iterative process of the axial heat transfer correction.
Iteration number e e ∗ δ Peak temperature (°C) 1  1.0 e - 08 740.4897 2 1.5748 e - 05 1.0 e - 08 740.4912 3 7.2384 e - 08 1.0 e - 08 740.4913 4 2.8562 e - 11 1.0 e - 08 740.4913 e e   ∗ = ∑ i n ( T i k - T i k - 1 ) / ∑ i n T i k and δ is set as 1.0e-08.Figure 10
Temperature field of the target area (considering the axial heat transfer).Figure11 shows the temperature distribution in radial direction at different distance. The high-temperature zone lies in the middle of the model and the heat transferred from the reactor is limited to the boundary area, which is mainly due to the low effective thermal conductivity. And thus, it reveals that the reactor and the target are uncoupled in thermal process.Figure 11
Temperature distribution in radial direction at 20 cm, 40 cm, 60 cm, and 80 cm.


## 3.1. Code Validation

Neglecting the correction of the axial heat transfer, the new code and Fluent were used to calculate the temperature field of the target because in Fluent the axial heat transfer can hardly be taken into consideration. In Fluent and the new code, the same physics model is adopted.In Fluent, the granular flow is also considered to be solid because the momentum equation for the fluid in Fluent is not suitable. The same control volume whose height is 5 mm and width is 14 cm is chosen by tracking which flow of the tungsten granules is simulated. The mesh of the control volume can be seen in Figure8. The two-dimensional transient solver is chosen. The time step size is set as 0.01 s and the number of time steps is 160 s. The initial temperature is set as 250°C. The upper and lower boundaries are set as adiabatic. The left boundary is set as the axis and the right boundary adopts the first boundary condition. The temperature of the right boundary changes from 380°C to 300°C. Moreover, the heat deposition in Section 2.4 is imported to Fluent with the help of UDF.Figure 8
Mesh of the control volume.The result shows that the temperature field calculated by these two codes is almost the same (Figure9). The temperature changes greatly near the inlet area and tends to keep unchanged in the following flow because the heat source concentrates near the inlet. Moreover, the peak temperature is 740.5°C calculated by the new code and 739.5°C calculated by Fluent. The result validates the governing equations and the discretization methods in the new code.Figure 9
Temperature field of the target area (neglecting the axial heat transfer). ((a) is calculated by the new code and (b) is calculated by Fluent.)
(a) (b)The average outlet temperature calculated by these two codes is 314°C. For the control volume, the average temperatures at the very beginning and the eighth second are 250°C and 314°C, respectively. To verify the heat source, the heat removed by the target is calculated by the following equation:(12)ΔQ=c∗m∗ΔT=c∗ρ∗π∗r2∗h∗ΔT.Inserting values into the equation,(13)ΔQ=130∗10062∗3.14∗0.142∗0.005∗314-250=25761J,P=ΔQΔt=257610.01=2.58MW.The result shows that the heat removed in the present model is 2.58 MW which is similar to the design index: 2.5 MW.


## 3.2. Temperature Field of the Target Area near the Active Zone

Considering the correction of the axial heat transfer, the new code is used to calculate temperature field of the target zone again. Figure10 shows the new temperature field. The peak temperature is also 740.5°C and the average outlet temperature is also 314°C. Table 3 shows the iterative process of the axial heat transfer correction. By 4 times of iteration, the results satisfy the demand. The results reveal that the axial heat transfer has a little effect on the temperature field. It is easy to understand because the temperature difference between successive time is small, as well as the effective thermal conductivity.Table 3
Iterative process of the axial heat transfer correction.
Iteration number e e ∗ δ Peak temperature (°C) 1  1.0 e - 08 740.4897 2 1.5748 e - 05 1.0 e - 08 740.4912 3 7.2384 e - 08 1.0 e - 08 740.4913 4 2.8562 e - 11 1.0 e - 08 740.4913 e e   ∗ = ∑ i n ( T i k - T i k - 1 ) / ∑ i n T i k and δ is set as 1.0e-08.Figure 10
Temperature field of the target area (considering the axial heat transfer).Figure11 shows the temperature distribution in radial direction at different distance. The high-temperature zone lies in the middle of the model and the heat transferred from the reactor is limited to the boundary area, which is mainly due to the low effective thermal conductivity. And thus, it reveals that the reactor and the target are uncoupled in thermal process.Figure 11
Temperature distribution in radial direction at 20 cm, 40 cm, 60 cm, and 80 cm.


## 4. Conclusions and Perspective

The dense granular target is a new concept for a high-power spallation target. In present work, a new thermal code based on the effective thermal conductivity model is developed to simulate its temperature field.(i)
Neglecting the axial heat transfer correction, the solver of the new code is validated by Fluent.(ii)
The target can remove 2.58 MW heat deposition with the peak temperature 740°C.(iii)
The axial heat transfer has a little effect on the temperature field.(iv)
The heat deposition in the target can hardly be transferred to the subcritical reactor and the heat transferred from the reactor is limited to the boundary area of the target, which means that the heat exchange between the target and the reactor can be neglected.More work needs to be done in the following studies:(i)
Experiment could be carried out to measure the stagnant thermal conductivity of helium-saturated porous beds of tungsten granules. Zhao and Yao [9] have done some similar experiments to measure the stagnant thermal conductivity of porous beds, which shows that the experiment may not be very complex or expensive.(ii)
The whole target can be taken into consideration, not only the area near active zone.



---

*Source: 5158610-2016-11-07.xml*