Simulating cyclic high heat flux loading of divertor cooling finger mock-up

Nuclear Engineering and Design 261 (2013) 317–325

Contents lists available at SciVerse ScienceDirect

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

Simulating cyclic high heat flux loading of divertor cooling finger mock-up Samo Koˇsmrlj ∗ , Boˇstjan Konˇcar Joˇzef Stefan Institute, Jamova Cesta 39, SI-1000 Ljubljana, Slovenia

h i g h l i g h t s • • • •

CFD simulations of high heat flux experiments. Comparison of temperature profiles from the experiments with the simulated profiles. Accuracy improved with new materials data and better approximation of heat flux. Transient analyses of cyclic heat loading using HTC.

a r t i c l e

i n f o

Article history: Received 3 February 2012 Received in revised form 6 February 2013 Accepted 7 March 2013

a b s t r a c t The cooling capability of the divertor mock-up under the real DEMO conditions can be tested on the highheat flux experiments. In this study the high heat flux (HHF) experiments, carried out at Efremov Institute are simulated using the computational fluid dynamics analysis. To accurately reproduce the experimental boundary conditions, the non-uniform distribution of the surface heat flux was taken into account and the heat losses were estimated. Thermal conductivity of the structure material was also found to be an important source of modelling error. The effects of the varied mass flow and inlet temperature boundary conditions on the wall heat transfer coefficient (HTC) distribution are analyzed by the means of steadystate simulation. Based on the obtained HTC, the transient analysis of cyclic surface heating was carried out. The predicted surface temperature cycles on the top of the cooling unit are validated against the experimental data. © 2013 Elsevier B.V. All rights reserved.

1. Introduction The design concept for future fusion reactor DEMO includes high heat flux (HHF) components that have to operate under extreme physical conditions. Although the exact mode of operation of DEMO reactor is not known yet, it is clear that the heat loads and neutron fluxes will be significantly higher than in what is currently the largest experimental fusion reactor under construction ITER (International Thermonuclear Experimental Reactor). Unlike ITER, DEMO is supposed to operate as a prototype power plant with higher energy conversion efficiency and should be capable of electricity generation, continuous operation and long lifetime. The divertor, as the highest thermally loaded HHF component in the fusion reactor takes care of cleaning of the plasma, which is a prerequisite for continuous reactor operation. About 15% of the total fusion power has to be removed by divertor, taking into account heat load peaks of over 10 MW/m2 and extensive volume heating by neutrons. Moreover, the total balance of the power plant depends on the efficient

∗ Corresponding author. Tel: +386 1 5885 260; fax: +386 1 5885 377. E-mail addresses: [email protected], [email protected] (S. Koˇsmrlj), [email protected] (B. Konˇcar). 0029-5493/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.nucengdes.2013.03.035

divertor operation, therefore several requirements have to be fulfilled, i.e. high heat removal capability and reasonably low pressure loss. Helium as a coolant and heat transfer agent is considered as the best option for DEMO divertor. It offers several advantages, such as chemical and neutronic inertness, and that the helium loop can be operated at much higher temperatures and lower pressures than for example water loop. Last but not least, the use of helium eliminates the risk of hydrogen generation, which cannot be excluded in the case of accident in water circuit, where hydrogen can be formed during the water-beryllium reaction. The main deficiency of the helium is its relatively low heat removal capability and a high pumping loss. However, the heat transfer can be improved by increasing the heat transfer surface (heat transfer promotors) or by employing turbulent impinging jets. Several advanced He-cooled divertor concepts for power plant application have been proposed in the past years as presented in the overview paper of Tillack et al. (2011). The cooling configurations considerably differ in size – from the largest with the characteristic length of the target plate of 1 m (Hermsmeyer and Malang, 2002), to the ARIES-CST-tube configuration with characteristic length in the order of 10 cm (Raffray et al., 2008) and eventually to the smallest cooling finger unit concept of KIT with the characteristic dimension of 1.5 cm (Norajitra et al., 2008a). The common feature of all

318

S. Koˇsmrlj, B. Konˇcar / Nuclear Engineering and Design 261 (2013) 317–325

concepts is that they use tungsten or tungsten alloy as a structural material. The development goes in the direction of smaller cooling units in order to minimize the thermal stress at the given heat flux load. In this study a modular divertor concept composed of numerous centimeter-scale cooling units, which are cooled by helium impinging jets is considered. Its development, carried out at Karlsruhe Institute of Technology (KIT) (Norajitra et al., 2007) is however still associated with several unknowns, especially from the point of view of cooling efficiency and constraints imposed by material properties. The divertor target plates should withstand peak heat loads of 10–15 MW/m2 . In addition, the divertor should survive around 1000 heat load cycles between the room temperature and higher limit operational temperature, which can exceed even 2000 ◦ C. To study and improve the heat removal capability of cooling fingers, experiments under the real DEMO operating conditions accompanied by numerical analyses are required. Therefore it is important to develop an accurate numerical model and compare it with the real experiments. The experimental campaign aiming to test the cooling finger mock-ups under extreme DEMO-like high heat flux conditions started several years ago with the first test series, carried out at Efremov Institute in St. Petersburg, Russia (Ovchinnikov et al., 2005; Norajitra et al., 2006). In 2008, a third series of experiments was conducted at Efremov Institute (Norajitra et al., 2008b), where an optimized cooling finger design with the rounded tile lower edge was used for the experimental mock-up. Thermo-mechanical simulations of the first experimental series have been performed on the simple adiabatic model with a uniform heat flux density (Kruessmann et al., 2007; Konˇcar et al., 2010; Simonovski et al., 2010; Widak and Norajitra, 2010). In the present paper the simulations of the third series of Efremov experiments are presented focusing on transient simulations of periodic heat flux loading accompanied with detailed representation of experimental boundary conditions and heat losses. The computational analyses were performed using the ANSYS CFX code (ANSYS CFX, 2009) and the results were validated on the experimental data.

2. Description of the experiment The combined helium loop and electron beam facility (TSEFEY facility) at Efremov Institute in St. Petersburg, Russia was designed for extremely high heat flux tests. The heart of the facility is a raster electron beam with a Gaussian intensity profile. The beam power was Pbeam = 14.5 kW and the beam diameter was 15 mm, expressed in full width at half maximum (FWHM). Connection between FWHM and the standard deviation , another common of measure of Gaussian function width, which is half the diameter √ the Gaussian function at inflection point, is FWHM = 2 2 ln 2. One-dimensional Gaussian function is shown in Fig. 1, together with its measures of width and the projected part of the function, which misses the target in case of raster pattern of the center of the beam covering only the cooling finger and its immediate neighborhood, as will be described in Section 2.1. Other parts of the facility, shown in Fig. 2, include helium cooling circuit for cooling the mock-up, shielding mask used to shield sides of the mock-up from unwanted heating, water cooling circuit for cooling the shielding mask and measuring equipment. The temperature of the heated mock-up surface was measured by an infrared (IR) camera and by a 2 color pyrometer. The raster pattern and the beam profile result in a non-uniform distribution of heat flux density on the cooling finger. Together with the fact that tungsten reflects some of the incident electron beam, this means it is difficult to estimate the absorbed power. Therefore in the experiment, only the power that was removed by the helium coolant was measured. Consequently, in the CFD analyses the absorbed power profile and heat losses had to be estimated

Fig. 1. Non-normalized Gaussian distribution function and its measures of width. The function has a width of FWHM = 15, which corresponds to  = 6.7, along with the projected extent of the cooling finger (dashed) and the part of the function, which gets deposited in case of having beam centered at the extreme distance from the axis of the cooling finger, as shown in Fig. 3 (hatched area).

and taken into account when applying boundary conditions. The absorbed power profile was calculated from the Gaussian shape of the beam and the raster pattern, which is shown in Fig. 3, along with the projected beam diameter (dashed line). As it can be seen in Fig. 1, even when the beam is at extreme right, a significant part of the beam is still deposited on the cooling finger (hatched area in Fig. 1). But due to the scanning pattern, that amount is lost because we do not continue moving the beam toward the right but rather turn the beam back to the left. This means that the central region of the cooling finger receives more heat flux than the outer regions and results in a non-uniform heat flux distribution. 2.1. Approximation of the heat flux density Commonly, for the purpose of CFD analyses, a uniform heat flux distribution over the top surface of the cooling finger is assumed (Kruessmann et al., 2007; Konˇcar et al., 2010). However, in the presented approach, we have modelled a realistic heat flux distribution in the experiment, taking into account the raster pattern of Gaussian shaped beam. The absorbed quasi steady state heat flux

Fig. 2. TSEFEY facility: rastered electron beam heats the cooling finger mock-up, which is cooled by a helium loop and shielded from unwanted additional heating on the sides by a water cooled mask. The surface temperature of the mock-up is measured by an IR camera and by an optical pyrometer. The absorbed power is calculated from the change of He temperature (Norajitra et al., 2007).

S. Koˇsmrlj, B. Konˇcar / Nuclear Engineering and Design 261 (2013) 317–325

319

Fig. 4. Heat flux distribution in case of 12 MW/m2 equivalent uniform absorbed heat flux. The absorbed heat flux ranges from 10.2 MW/m2 at the corner farthest from center to 13.3 MW/m2 in the center.

Fig. 3. Raster pattern of electron beam (red line) and its projected diameter at FWHM (dashed black line), together with the cooling finger extents (solid black line). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

distribution was calculated by superposition of the beam heat flux distributions over the top surface of the mock-up (Eq. (1)):

 1 exp 2 2 N N

f˙ (x, y) =



f (x, y), defined by beams f˙ (x, y) with a single distribution function  Eq. (5). We can see that the relative error made by approximation is below 3% at any point on the surface. (x, y) =

 2

−

n

(x − xn )2 + (y − yn ) 2 2

.

(1)

Here xn and yn are coordinates corresponding to the positions of beam center, and  is the width of the Gaussian distribution function of the beam. The resulting absorbed heat flux distribution function f˙ (x, y)) can be approximated with a single Gaussian

 = 14.3 mm. function  f (x, y), which has a width of    2 2 1 x +y f˙ (x, y) ≈ exp − = f (x, y). 2 2 2 2

j0 AS

2 2



−

exp

x 2 + y2

2 2

(2)



.

f (x, y) f˙ (x, y) −  . f˙ (x, y)

(5)

2.2. Heat losses Cooling fingers, heated to temperatures up to 2000 K, lose some energy due to radiation losses. Radiation losses were calculated from the Stefan–Boltzmann law J ∗ = (T )T 4 ,

(6)

with total emissivity of tungsten modeled as

The approximation distribution function  f (x, y) is then renormalized, so that its integral over the surface of the mock-up is 1, and multiplied by the equivalent uniform heat flux to get the distribution of heat flux for nonuniform case. f (x, y) = Pabs (x, y) ≈ j0 AS

Fig. 4 shows distribution of heat flux, according to Eq. (1), and Fig. 5 shows the relative error , made by approximating the sum of

(3)



(T ) = 0.02 +

 −

5.0 × 10−5 K

4.0 × 10−11 K3





T 3.



T+

1.2 × 10−7 K2

 T2

(7)

The values of total emissivity used in interpolation were taken from (Omega Engineering Inc, 2011). Interpolation is valid in the range from 300 K to 2300 K. The resulting heat losses for selected simulated cases are shown in Table 1.

Here j0 is the absorbed power flux equivalent, which is the absorbed power flux in case of uniform power flux, S = 2.74 × 10−4 m2 is the surface of cooling finger and A = 5.11 is renormalization constant, f (x, y) over the surface of mock-up is 1. so that the integral of  The approximation of the sum of beams by a single heat flux density distribution can be justified by comparing the time, required for repeating the entire pattern and the timescale for the disturbance to travel across the cooling finger. The time required to repeat the entire pattern was 1.8 ms (Mueller et al., 2011). The timescale for a disturbance to travel across the cooling finger is given by the solution of a heat transport equation (Eq. (4)).



T (t, L) ∝ exp

−

t cp L2



 t

= exp −



.

(4)

Values for tungsten finger at 300 K are  ≈ 120 W/mK,  = 19, 250 kg/m3 , Cp ≈ 160 J/kg K and L = 8 mm, which gives the characteristic time  = 1.8 s, several orders of magnitude greater than the time required for entire raster pattern.

Fig. 5. Relative approximation error done by approximating the rastered beam with a single Gaussian distribution function (Eq. (5)).

320

S. Koˇsmrlj, B. Konˇcar / Nuclear Engineering and Design 261 (2013) 317–325

Table 1 Absolute and relative heat losses (total losses/absorbed power) of selected models. Case

Rad. top (W)

Rad. side of tile (W)

Rad. thimble (W)

Total losses (W)

Absorbed power (W)

Relative losses (%)

mfr 13.5 g/s, jabs 8 MW/m2 , THe 680 K mfr 13.5 g/s, jabs 12 MW/m2 , THe 680 K mfr 13.5 g/s, jabs 12 MW/m2 , THe 830 K

10.62 31.14 44.82

10.14 24.72 38.88

0.90 1.14 2.64

21.66 57.00 86.34

2193.60 3290.40 3290.40

1.0 1.7 2.6

the cell. For the solver used, values above 0.1 are considered to be acceptable. The mesh was oriented in such way that it followed the fluid flow direction and the mesh in the near-wall boundary layer was refined, so that the first node in the fluid domain was placed at the non-dimensional distance y+ ≈ 1 from the wall. The mesh is divided into fluid and solid domains (tile, which receives the heat flux from the plasma, thimble, which supports the tile and cartridge, which forms helium jets). The mesh, shown in Fig. 7, consisted of 1.5 million elements. 3.1. Governing equations Heat transfer in the solid domain is defined by the heat conduction equation:

∂T  ∂ T , = cp  ∂t ∂xi2 2

(8)

i

Fig. 6. Locations of characteristic points together with heat fluxes (absorbed heat flux, radiated heat flux and removed heat flux).

As a part of the analysis of heat losses, the temperature profiles of two models were compared: adiabatic model and the model with radiation losses. Temperatures, used for comparison, were taken at three characteristic points along the axis of the cooling finger. The cooling finger is thinnest at the axis and the heat flux density is the greatest there, therefore the temperatures are the highest on the axis of the cooling finger. T1 is located at fluid–thimble interface, T2 is located at the tile-thimble interface (brazing layer) and T3 is located at the top of the tile, indicating the highest temperature in the cooling finger. The characteristic points are shown in Fig. 6 and the temperatures at characteristic points for both models are shown in Table 2.

where , cp and  are the heat conduction coefficient, specific heat capacity at constant pressure and the density of the solid, respectively. For fluid equations, the Einstein summation convention has been used (if an index variable appears more than once in a single term, that implies we are summing over all possible values of the index variable). Using this formulation, the fluid motion and heat transfer in the fluid domain are calculated by the continuity equation

∂ ∂(uj ) = 0, + ∂t ∂xj

(9)

3. CFD steady state analyses After the distribution of heating power and the heat losses are determined, we can start setting up our model. First we must create computational mesh, which allows us to solve the heat transfer and fluid motion equations. For creating the computational mesh, hexahedral meshing was used. One-sixth cutout was simulated, taking into account the axial symmetry of the mock-up. The mesh was created with a minimum quality for a single mesh element of 0.35, which is a weighted average of determinant, angle and warpage of Table 2 Comparison of temperature at characteristic points for the model with rod thermal conductivity and mfr = 13.5 g/s, Tin = 830 K and jabs = 12 MW/m2 . Case

T1 [K]

T2 [K]

T3 [K]

Adiabatic model Model with radiation losses

1194 1184.1

1327.9 1314.2

1975.1 1946.5

Fig. 7. Computational mesh of the one sixth cutout: Helium domain is colored blue, WL10 domains (the cartridge has holes which form helium jets, thimble supports the tungsten tile) are colored brown and tungsten domain is colored light green. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

S. Koˇsmrlj, B. Konˇcar / Nuclear Engineering and Design 261 (2013) 317–325

• Specific dissipation rate:

the momentum equation (in absence of external forces)



∂(ui ) ∂(ui uj ) ∂p ∂ =− +

eff + ∂xj ∂t ∂xi ∂xj and the total energy equation

∂(htot ) ∂p ∂ (ui htot ) ∂ = − + ∂t ∂t ∂xi ∂xi



+

∂ ui eff ∂xi





∂ui ∂uj − ∂xj ∂xj



∂ui ∂uj + ∂xi ∂xj

∂T  ∂xi



,

(10)

+ 2(1 − F1 )



(11)

Here xi is the spatial coordinate, ui is component of the velocity field, htot is the total enthalpy, related to the static enthalpy h(T, p) by htot = h +

ui ui , 2

(12)

and p is effective pressure, defined as p = p +

2 2 ∂ui . k + eff 3 3 ∂xi

(13)

Here k is turbulence kinetic energy. A Boussinesque approximation is used to define turbulent fluctuations in the fluid. Hence, the effective viscosity in momentum equation can be defined as a sum of molecular ( ) and the so-called turbulent viscosity ( t ):

eff = + t .

(14)

In the present approach, the turbulent viscosity is modeled by a Shear-Stress Transport (SST) model (Menter, 1994). The formulation of SST model has several advantages over the classical k −  turbulence models, since its formulation enables that it can be used for free shear and boundary flows, as is the case at impinging jets. SST model combines the commonly used k −  and k − ω models. Turbulence kinematic viscosity is linked to the turbulence kinetic energy k and specific dissipation rate ω via the relation, characteristic for SST model: T =

  ∂(uj ω) ∂(ω) ∂

t ∂ω ω =

+ + ˛ P − ˇω2 + ω ∂xj k ∂t ∂xj ∂xj



.

t a1 k . =  max(a1 ω, SF2 )

(15)



2Sij Sij , where Sij are components of the mean S is defined as S = strain rate tensor. SST model is characterized by the above expression for kinematic viscosity. Values for k and ω come from the equations for the turbulence kinetic energy and specific dissipation rate: • Turbulence kinetic energy:

  ∂(k) ∂(uj k)

t ∂k ∂ =

+ + P − ˇ kω. + k ∂xj ∂t ∂xj ∂xj

321

(16)

ω ∂k ∂ω . ω ∂xj ∂xj

(17)

Detailed description of the blending functions F1 and F2 and of the constants, used in the above equations, can be found in Menter (1994). The complete set of equations for fluid and solid domains is solved for steady-state conditions only, using the ANSYS-CFX code (ANSYS CFX, 2009). In previous studies of Draksler and Konˇcar (2011, 2010) turbulence models and heat transfer predictions were validated on a standard benchmark of a simple jet impingement test case (Baughn and Shimizu, 1989) and also on the multiple jet impingement test case (Geers et al., 2004). The results show reasonably good predictions of flow and heat transfer phenomena. Furthermore, the results of other authors comparing different turbulence models confirmed that the SST model performs well in the case of multiple impinging jets (Rao et al., 2009; Zu et al., 2009). 3.2. Thermal conductivity of structural materials Diverter materials are subjected to extreme thermal conditions. Therefore it is very important to accurately describe the material properties. Recently, new measurements of pure tungsten and 1 wt% La2 O3 doped tungsten (WL10) thermal conductivities became available (Rieth et al., 2010). Sacrificial top layer of cooling finger (tile) is made of pure tungsten and WL10 is the material used for cooling finger structure (thimble). The thermal conductivities reported in (Rieth et al., 2010) differ significantly from the frequently used reference conductivities, given in ITER materials handbook (ITER, 2001). Higher thermal conductivities mean lower temperature gradients in the structure and consequently lower temperature in the cooling finger mock-ups. Lower temperature gradients also mean lower thermal stress in the cooling fingers. The effect of different conductivities on temperature distribution inside a cooling finger along the axis is shown in Table 3. According to (Rieth et al., 2010), the thermal conductivities depend on the manufacturing process (plate, rod). The highest thermal conductivity and consequently the lowest temperature for the WL10 material is obtained for the thimble made from rod material. We can see that the different thermal conductivities mean significant differences in temperature, especially at the top of the finger. The differences can be almost 100 K, or up to 5%, at top of tile, and over 30 K at the tile–thimble interface. Tile and thimble were brazed together at the temperature of 1050 ◦ C (1323 K) and it is very important that this temperature limit is not exceeded at the tile-thimble interface. Otherwise the brazing layer starts melting and consequently, the heat transfer through the cooling finger

Table 3 Effect of different thermal conductivities (plate and rod material from (Rieth et al., 2010) and data from ITER handbook (ITER, 2001), on temperature at various locations along the axis of cooling finger (see Fig. 6) for several cases with different boundary conditions: mass flow rate of helium (mfr), absorbed heat flux equivalent (jabs ) and inlet helium temperature (THe ). Rod

Plate

ITER handbook

Case

T1 [K]

T2 [K]

T3 [K]

T1 [K]

T2 [K]

T3 [K]

T1 [K]

T2 [K]

T3 [K]

1. mfr = 12 g, jabs = 8 MW/m2 , THe = 680 K 2. mfr = 13.5 g, jabs = 8 MW/m2 , THe = 680 K 3. mfr = 13.5 g, jabs = 8 MW/m2 , THe = 830 K 4. mfr = 13.5 g, jabs = 12 MW/m2 , THe = 680 K 5. mfr = 13.5 g, jabs = 12 MW/m2 , THe = 830 K

949.9 933.2 1068.1 1056.0 1184.1

1032.2 1014.5 1153.1 1173.0 1306.3

1438.8 1420.1 1566.5 1798.0 1946.5

950.4 933.8 1070.6 1056.2 1185.4

1045.1 1027.3 1169.5 1200.4 1335.9

1450.8 1432.0 1585.0 1823.9 1969.6

951.2 934.7 1070.5 1057.5 1186.1

1050.2 1032.6 1172.0 1206.1 1339.3

1487.2 1467.6 1623.5 1892.7 2039.7

322

S. Koˇsmrlj, B. Konˇcar / Nuclear Engineering and Design 261 (2013) 317–325

Table 4 Mesh independence study. Mesh

Number of mesh elements in the fluid domain

T1 [K]

T2 [K]

Coarse Medium Fine

439,308 953,864 1,649,984

1304.9 1289.5 1285.9

1177.5 1160.4 1153.6

HTCavg [W/m2 K]

18,029 19,240 19,819

deteriorates. This finally leads to complete failure of the cooling finger. 3.3. Mesh refinement study In order to obtain mesh independent results, the mesh refinement study was carried out for three different meshes of the fluid domain. In the study, the temperatures in points T1 and T2 described in Fig. 6 were compared. Additionally, the average heat transfer coefficient (HTCavg ) on the fluid-solid interface was calculated. The mesh was refined in particular in the gap region close to the fluid–thimle interface, where the highest velocity and temperature gradients can be expected. The three meshes are named as “coarse”, “medium” and “fine”. Number of mesh elements in the fluid domain and comparison of the results on the different meshes is shown in Table 4. For better representation of the temperature distribution on different meshes, the fluid temperature on the thimble–helium interface is shown in Fig. 8. In Fig. 8, the quotient l/l0 represents the non-dimensional distance from the finger central vertical axis (see Fig. 6). As may be observed, the temperature distributions are almost the same for the medium and fine mesh. Results from Table 4 and Fig. 8 show that satisfactory mesh independence can be obtained already for the medium mesh, therefore all further analyses were performed on this mesh. 3.4. Temperature profile on the surface of the tile After estimating the heat losses, it is possible to estimate the power that was absorbed on the top of the tile from the power that was removed by helium. The absorbed power is used as a boundary condition on top of the tile. With all the boundary conditions defined, we consider the computational model complete and ready for the simulations. Fig. 9 shows the simulated and experimental temperature distributions along the line, shown in Fig. 10 for the case with mass flow rate 13 g/s, inlet temperature 830 K and absorbed power equivalent of 10 MW/m2 .

Fig. 8. Temperature along the line, projected to the thimble–helium interface for different meshes.

Fig. 9. Comparison of experimental and simulated temperature profile for the uniform and Gaussian shape of incident beam and with added reflected Gaussian beam (Ovchinnikov et al., 2008).

The experimental profile shows asymmetry. This asymmetry can be explained as follows: while the experimental data was being recorded, the heating beam was still on for a fraction of exposure time and a part of it reflected into camera. This would cause the camera to measure higher temperatures in the region where the beam was focused. To account for this, the heat flux into the camera was recalculated from the measured temperature, according to the Stefan–Boltzmann law (Eq. (6)). The heat flux, radiated into the camera from the cooling finger as a result of the temperature distribution calculated in the simulation, was subtracted and the remaining difference was approximated by a Gaussian shaped beam with the same width as the electron beam, using the least squares method. The beam which fits best to the unaccounted heat flux corresponds to a 830 W beam which is centered on the measurement line about 4 mm away of the central axis of the cooling finger, or to a 14.5 kW electron beam with its focus located at 25 mm from the measurement line. If a 14.5 kW electron beam with FWHM of 15 mm was directed straight into the camera, the camera would record temperatures up to 5500 K.

Fig. 10. Temperature on the top of cooling finger (hexahedral structure in the middle) and shielding mask (circular structure surrounding the cooling finger). Also shown is the line along which the temperature, shown in Fig. 9, was measured. It is likely, that the hotter region at bottom left is a result of a reflected beam (Ovchinnikov et al., 2008).

S. Koˇsmrlj, B. Konˇcar / Nuclear Engineering and Design 261 (2013) 317–325

323

To evaluate, which of the calculated profiles fits the experimental results better, we can use the Root Mean Square (RMS) difference defined in Eq. (18).

 ıTRMS =

(Tsim (x) − Texp (x))2 dx



Texp (x)2 dx

.

(18)

Here Tsim (x) denotes the temperature, calculated from the sum of the heat flux from the temperature, which was calculated in the simulation and the best fit (using the least squares method) of the Gaussian shape beam to the remaining heat flux. Texp is the measured temperature from the actual experiment and x is the coordinate along the measurement line. The ıTRMS for the uniform absorbed heat flux is 0.74% and for the Gaussian beam it is 0.63%. Another possible explanation is that the brazing layer failed and as a result, the thermal conductivity at the fault is decreased, which would result in higher temperatures above the failure of the brazing layer. Additionally, the simulated data is affected by approximation error, though the relative error in temperature, calculated in simulation is proportional to the relative error in heat flux. This was estimated to be at most 5 K for 12 MW/m2 of equivalent absorbed heat flux, which was estimated by comparing two absorbed power density distributions with slightly different heat fluxes. Experimental data was recorded with an IR camera with a ± 5% accuracy.

3.5. Heat transfer coefficient calculation Simulating turbulent flows with CFD programs is computationally very expensive. In order to capture the details of the turbulent flow with sufficient precision, a very fine mesh must be used in the region of helium jets. This means that in order for the courant number to stay sufficiently low, very short timestep for a transient simulation must be used (≈0.1 ␮s). For this reason, the fluid flow was only simulated in steady state analyses to extract the heat transfer coefficient (HTC) distribution on the fluid-solid interface. HTC distributions were then used in transient analyses as a boundary condition on the bottom side of the thimble. The effect of varying boundary conditions, i.e. mass flow rate of the coolant, inlet temperature of the coolant, total absorbed power and distribution of the absorbed power has been estimated (Koˇsmrlj and Konˇcar, 2011) and it was found that only varying the fluid parameters, i.e. mass flow rate and inlet temperature of the coolant significantly changed the HTC distribution. As it can be seen from Table 3 by comparing the cases with the same mass flow rate and the same absorbed heat flux equivalent but with different inlet temperatures (cases 2 and 3 or cases 4 and 5), we are able to estimate the effect of the inlet temperatures on the temperature at the fluid–thimble interface. For example, let us compare the cases 4 and 5 with the plate thermal conductivity. Temperature at T1 (see Section 2.2) is 1056.2 K for the case with the inlet temperature of 680 K and 1185.4 K for the case with the inlet temperature of 830 K. Difference between temperatures at T1 is 129.2 K. Difference in inlet temperatures is 150 K, which means that the case with higher inlet temperature has slightly higher HTC (at the region near the axis of the cooling finger), and the difference in HTC contributes −0.13 K to the temperature difference at the fluid–thimble interface for every 1 K difference in inlet temperature. As it will be discussed in Section 4, the inlet temperature during the simulated cyclic heat flux experiment varied for about 40 K, and can therefore contribute about 5 K to the maximum error of the interface temperature, if a constant HTC distribution is used instead of using different HTC distribution, depending on inlet temperature, for every time step in transient simulation.

Fig. 11. Time dependence of heat flow rate (red), transient average of heat flow rate (green) and heat flow rate from the steady state simulation (blue). It can be seen that the heat flow rate can be considered constant already in a few 100 ␮s. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

4. Transient analyses After all the above preparations of computational model, we are finally ready to perform transient analyses. Transient analyses were done on the tile and thimble regions of the cooling finger structure only, where the heat conduction equation (Eq. (8)) was solved. This greatly reduces the mesh size, because the fluid domain – especially in the near-wall region – needs to be meshed with a very fine mesh in order to be able to resolve the velocity and temperature variations in the fluid, which occur on very small spatial scales. Also, in the fluid region the number of numeric operations is much higher because of the complexity and number of equations solved. Altogether, this means that one computational iteration for the cooling finger with the fluid domain (medium mesh) takes about 30 times more time than the computational iteration for the cooling finger without the fluid domain. The time scale of the fluctuations in the heat flow rate from the solid domain into fluid domain, due to the unsteady helium flow, is on the order of microseconds. This is calculated in a short transient simulation with the fluid domain included and is shown in Fig. 11. Transient calculation of only 100 ␮s of the simulated time, where both the heat transfer and fluid motion equations were solved in parallel on 24 high performance Intel Xeon 5670 processors, took about two weeks of real time. On the other hand, transient simulation, by orders of magnitude longer (610 s), where solely the heat transfer in solid domain was solved, took only two days of CPU time. The timescale for heat transfer fluctuations (≈100 ␮s) is much smaller than the time scale of temperature fluctuations in the solid domain (≈2 s), described by Eq. (4). Therefore it is safe to assume that the average heat transfer on the solid–fluid interface does not depend on the stochastic fluctuations in the helium flow. HTC distribution was therefore assumed to be constant throughout the entire simulated case, neglecting the effect of varying inlet temperature of the helium, as justified in Section 3.5. In transient analyses, the mass flow rate of helium and the shape of absorbed heat flux distribution function were considered to be constant, while the total absorbed power and the helium inlet temperature were varied. The soft ramp case was investigated with the ramp cycle as follows: 25 s for the heating power to reach its maximum from zero, 25 s at full power, 25 s to reach zero from maximum and 25 s at zero. The reported experimental data however showed significant discrepancy between the reported heating power and the calculated power, removed by the helium. This indicates, that there were

324

S. Koˇsmrlj, B. Konˇcar / Nuclear Engineering and Design 261 (2013) 317–325

the estimated absorbed heat flux was 10.1 MW/m2 . To account for the possible additional losses, another cycle was simulated with the absorbed power increased by 10%. Fig. 12 shows the time dependency of removed heat flux equivalent together with inlet temperature of helium. Figs. 13 and 14 show comparison of experimental and simulated maximum temperature on the top of tile for the case where no additional heat losses except for the losses due to thermal radiation were assumed (Fig. 13), and for the case where additional losses have also been taken into account (Fig. 14). Also note that IR camera, used to record the temperature, has a lower operating limit, hence the mismatch at lower temperatures. 5. Conclusions

Fig. 12. Heat flux equivalent and cooling temperature in a simulated cycle.

The obtained results show reasonable agreement with the experimental data. With a good computational model we are able to accurately simulate the conditions inside the cooling finger, which provides a unique insight into the limitations of a cooling finger design and allows us to improve the parameters for relatively low computational cost. This means making development and optimization of the cooling finger design much easier and faster. The computational models have been gradually improved over time and can nowadays provide a good approximation of actual experimental conditions and therefore one can feel safe to use the results of computational analyses to support the design process. Acknowledgements The authors gratefully acknowledge P. Norajitra from Karlsruhe Institute of Technology for valuable discussions. The authors also acknowledge the financial support by the Slovenian Research Agency and the support by the European Commission within the framework of European Fusion Development Agreement (EFDA) and the Slovenian Fusion Association (SFA).

Fig. 13. Comparison of experimental and simulated maximum temperature on top of the tile.

additional heat losses. For example, the reported absorbed heat flux in the simulated cycle was around 11.0 MW/m2 (Ovchinnikov et al., 2008), but the heat flux, calculated from the outlet helium temperature was around 9.8 MW/m2 . Together with the radiated power,

Fig. 14. Comparison of experimental and simulated maximum temperature on top of the tile, where heat flux in simulation was increased by 10% to account for the discrepancy between the sum of removed power and calculated radiated power and the reported absorbed power.

References ANSYS CFX, 2009. Release 12.1, ANSYS CFX-Solver Theory Guide, ANSYS. Baughn, J.W., Shimizu, S., 1989. Heat transfer measurements from a surface with uniform heat flux and an impinging jet. J. Heat Transfer 111, 1096–1098. Draksler, M., Konˇcar, B., 2011. Analysis of heat transfer and flow characteristics in turbulent impinging jet. Nucl. Eng. Des. 241 (4), 1248–1254. Draksler, M., Konˇcar, B., 2010. Heat transfer and jet interaction for different arrays of impinging jets. In: Proceedings of Nuclear Energy for New Europe, NENE2010, September 6–9, Portoroz, Slovenia. ´ K., 2004. Experimental investigation of Geers, L.F.G., Tummers, M.J., Hanjalic, impinging jet arrays. Exp. Fluids 36, 946–958. Hermsmeyer, S., Malang, S., 2002. Gas-cooled high performance divertor for a power plant. Fusion Eng. Des. 61–62, 197–202. ITER Material Properties Handbook, 2001. Konˇcar, B., Norajitra, P., Oblak, K., 2010. Effect of nozzle sizes on jet impingement heat transfer in He-cooled divertor. Appl. Therm. Eng. 30, 697–705. Koˇsmrlj, S., Konˇcar, B., 2011. CFD Analyses of EFREMOV Institute Cooling Finger Experiments, IJS-DP-10812, Joˇzef Stefan Institute. Kruessmann, R., Messemer, G., Norajitra, P., Weggen, J., Zinn, K., 2007. Validation of computational fluid dynamics (CFD) tools for the development of a heliumcooled divertor. Fusion Eng. Des. 82, 2812–2816. Menter, F.R., 1994. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA J. 32, 1598–1605. Mueller, G., Gervash, A., Norajitra, P., 2011. EB loading pattern for HHF testing of mock-ups, personal communication. Norajitra, P., Abdel-Khalik, S., Giancarli, L.M., Ihli, T., Janeschitz, G., Malang, S., Mazul, I.V., Sardain, P., 2008a. Divertor conceptual designs for a fusion power plant. Fusion Eng. Des 83, 893–902. Norajitra, P., Gervash, A., Giniyatulin, R., Ihli, T., Krauss, W., Kruessmann, R., Kuznetsov, V., Makhankov, A., Mazul, I., Ovchinnikov, I., 2006. He-cooled divertor for DEMO: experimental verification of the conceptual modular design. Fusion Eng. Des. 81, 341–346. Norajitra, P., Giniyatulin, R., Ihli, T., Janeschitz, G., Krauss, W., Kruessmann, R., Kuznetsov, V., Mazul, I., Widak, V., Ovchinnikov, I., Ruprecht, R., Zeep, B., 2007. He-cooled divertor development for DEMO. Fusion Eng. Des. 82, 2740–2744.

S. Koˇsmrlj, B. Konˇcar / Nuclear Engineering and Design 261 (2013) 317–325 Norajitra, P., Krauss, W., Reiser, J., Widak, V., 2008b. He-cooled divertor development: the 3rd high-heat-flux test series of HEMJ 1-finger mock-ups and manufacturing and thermo-hydraulic tests of 9-finger module. Report for TASK of the EFDA Technology Programme, FZK. Raffray, A.R., El-Guebaly, L., Ihli, T., Malang, S., Wang, X., The ARIES-CS Team, 2008. Engineering design and analysis of the ARIES-CS power plant. Fusion Sci. Technol. 54, 725–746. Rao, G.A., Kitron-Belinkov, M., Levy, Y., 2009. Numerical analysis of a multiple jet impingement system. ASME Paper No. GT-2009-59719. Omega Engineering Inc One Omega Drive, Stamford, Connecticut, USA, Table of Total Emissivity. Available from: http://www.omega.com/ temperature/z/pdf/z088-089.pdf (accessed 12.09.11). Ovchinnikov, I., Giniyatulin, R., Ihli, T., Janeschitz, G., Komarov, A., Kruessmann, R., Kuznetsov, V., Mikhailov, S., Norajitra, P., Smirnov, V., 2005. Experimental study of DEMO helium cooled divertor target mock-ups to estimate their thermal and pumping efficiencies. Fusion Eng. Des. 73, 181– 186.

325

Ovchinnikov, I., Giniyatulin, R., Norajitra, P., 2008. Thermo-hydraulic and surface temperature data of the 3rd series of HHF 1-finger mock-up tests, personal communication. Rieth, M., Armstrong, D., Dafferner, B., Heger, S., Hoffmann, A., Hoffmann, M.D., Jaentsch, U., Kuebel, C., Materna-Morris, E., Reiser, J., Rohde, M., Scherer, T., Widak, V., Zimmermann, H., 2010. Tungsten as a structural divertor material. Adv. Sci. Technol. 73, 11–21. Simonovski, I., Konˇcar, B., Cizelj, L., 2010. Thermo-mechanical analysis of a DEMO divertor cooling finger under the EFREMOV test conditions. Fusion Eng. Des. 85, 130–137. Tillack, M.S., Raffray, A.R., Wang, X.R., Malang, S., Abdel-Khalik, S., Yoda, M., Youchison, M., 2011. Recent US activities on advanced He-cooled W-alloy divertor concepts for fusion power plants. Fusion Eng. Des. 86, 71–98. Widak, V., Norajitra, P., 2010. Optimization of He-cooled divertor cooling fingers using a CAD-FEM method. Fusion Eng. Des. 84, 1973–1978. Zu, Y.Q., Yan, Y.Y., Maltson, J.D., 2009. CFD prediction for multi-jet impingement heat transfer. ASME Paper No. GT-2009-59488.