Cite this Article
Huang Yan, Sun Ji-Zhong, Cai Juan, Sun Zhen-Yue, Sang Chao-Feng, Wang De-Zhen. Numerical study of influence of J × B force on melt layer under conditions relevant to ITER ELMs. Chinese Physics B, 2019, 28(4): 045201  
Permissions
Numerical study of influence of J × B force on melt layer under conditions relevant to ITER ELMs
Huang Yan1, Sun Ji-Zhong2, †, Cai Juan3, Sun Zhen-Yue2, Sang Chao-Feng2, Wang De-Zhen2
School of Information Science and Engineering, Dalian Polytechnic University, Dalian 116034, China
Key Laboratory of Materials Modification by Laser, Ion and Electron Beams (Ministry of Education), School of Physics, Dalian University of Technology, Dalian 116024, China
School of Physics and Electronic Technology, Liaoning Normal University, Dalian 116021, China

 

† Corresponding author. E-mail: jsun@dlut.edu.cn

Project supported by the Scientific Research Foundation of Liaoning Province, China (Grant No. 2016J027), the Open Research Project of Key Laboratory of Materials Modification by Laser, Ion and Electron Beams, Ministry of Education (Grant No. KF1705), and the National Key Research and Development Program of China (Grant Nos. 2017YFA0402500, 2017YFE0300400, and 2017YFE0301200).

Abstract

The influence of the J × B force on the topographical modification of W targets during a type-I-like ELM in ITER has been studied numerically. A two-dimensional (2D) fluid dynamics model is employed by solving liquid hydrodynamic Navier-Stokes equation with the 2D heat conduction equation in addition to driving forces for surface topography, such as surface tension and pressure gradient, the J × B force is particularly addressed. The governing equations are solved with the finite volume method by adequate prediction of the moving solid-liquid interface. Numerical simulations are carried out for a range of type-I ELM characteristic parameters. Our results indicate that both the surface tension and the J × B force contributes to the melt motion of tungsten plates when the energy flux is under , the surface tension is a major driving force while the pressure gradient is negligible. Our results also indicate that the J × B force makes the small hills grow at different rates at both the crater edges under a type-I-like ELM heat load with a Gaussian power density profile.

PACS: 52.55.Fa;52.55.Rk;52.40.Hf;52.65.Kj
Keyword:ELMs;tungsten divertor plates;J × B force;melt motion
1. Introduction

The ELMy H-mode is considered to be the standard operating scenario for ITER.[1] Heat loads brought by type-I edge localized modes (ELMs) impact on the ITER divertor plates,[2] which may exceed the tungsten (W) melting threshold and form a melt layer. The melt layer is subjected to different driving forces, such as plasma pressure gradient, surface tension, and Lorentz force ( ),[3–6] etc., and it will then move. Melt layer propagation will increase surface roughness, which reduces W thermal conductivity, degrades its mechanic strength, initializes arcing,[7] and even causes droplet splashing. The melt motion will be the major mechanism of W divertor erosion in the next generation of tokamaks. Therefore, it is very important to take the melt motion into account to study the W thermal erosion caused by type-I ELMs. Currently, melting and the melt motion of W plates have been studied in tokamaks (e.g., JET and TEXTOR).[8] However, plasma–material interaction (PWI) is much stronger in ITER than in any existing fusion device and there is a big gap in our present understanding of PWI in ITER and in the present tokamaks. Therefore, extrapolations from current devices to ITER are uncertain, including transient peak thermal loads during ELMs.[9]

Theoretical study of the W erosion is growing in importance. Federici et al.[10–12] built a one-dimensional model to study the W divertor erosion under ELM-like heat loads. Hassanein et al.[13,14] evaluated the W divertor thermal erosion using the code HEIGHTS. However, the W divertor erosion owing to melt motion has not been included in these studies. Bazylev et al.[15–20] evaluated the thermal erosion of the W divertor plates under type-I ELMs using the code package MEMOS. To simplify the solution, they separated the heat conduction equation into two one-dimensional heat conduction equations: the first equation represents the heat transfer of W target in the vertical direction (through the liquid and the solid phases), and the second represents the heat transfer along the horizontal direction of the W target due to melt motion. To the best of our knowledge, no detailed study has reported on the W melt erosion under an energy flux of .

In this paper, we extend the previous works[21,22] to study the influence of force on the melt motion during the type-I ELMs under a set of plasma parameters similar to those of ITER. In the model, typical interaction forces (pressure, surface tension and force) responsible for melt layer motion are taken into account. Simulations are carried out for type-I ELMs under a range of characteristic parameters. In addition, an improved numerical method has been used to solve the governing equations.

2. Model description

In ITER, the ELMs heat loads are foreseen in a range of on a time scale of about 0.1 ms–1 ms, and inter-ELMs thermal loads are expected .[23,24] In this range of heat loads, the melt layer will be a thin boundary layer—far thinner than the target thickness. Therefore, it is a shallow water circulation problem, in which a velocity component parallel to the surface (a melt velocity averaged over the melt layer thickness) can be used to describe the melt motion.[25] The melt layer propagation is shown schematically in Fig. 1 under a Gaussian power density profile heat load. As shown in Fig. 1, the melt layer flows along the x-coordinate axis (the poloidal direction of ITER), and the z coordinate is orthogonal to the x axis (into the layer). The y axis is the opposite direction of the toroidal magnetic field. All of the key surface heat transfer processes (e.g., evaporation, melting, radiation, and melt motion) are considered in the present model.

Fig. 1. Sketch of melt layer and melt motion of the target.

Driving forces of the melt motion, such as surface tension, pressure, and Lorentz force due to external currents, are considered. Meanwhile, gravity is ignored in the present model because it is a small amount compared with the force and surface tension (our simulations indicate that the average value of the force at each point is about ten times larger than the gravity, and the average value of the surface tension is about ten times and even more than gravity). In addition, the model assumes the liquid to be incompressible. The equations of conservation of mass, momentum and energy can be written as follows: where p, c, , and T are the material density, specific heat capacity, velocity, and temperature, respectively. The temperature dependence of thermal conductivity k is described by a function ,[26] where a and b are empirical parameters,[27] Q is the sum of the volumetric heating, and p is the pressure across the target surface.

The boundary conditions at the liquid–vapor boundary are represented by equations as follows: Here equation (4) is the energy balance equation at the surface, where F pl is the surface heat load, L v is the latent heat of vaporization, V v is the evaporation front velocity, and T D is the unperturbed divertor target surface temperature. Equation (5) gives the force balance condition at the liquid–vapor interface, where the term on its left-hand side (LHS) is the friction force, the term on its right-hand side (RHS) is the surface tension gradient, α is the surface tension coefficient, and μ is the melt viscosity.

At the melting front ( ), the classic Stefan boundary condition[28] is set for the solution. Here, the subscript s refers to the solid, and subscript l refers to the liquid, L m is the latent heat of melting, and V m(t) is the melting front velocity.

The shallow water approximation allows us to model the propagation of melt layer using a modification of the St Venant equations, The second term on the LHS of Eq. (7) describes the velocity of melt layer thickness change due to melt motion (where is the melt velocity averaged over the melt layer thickness, h is the melt layer thickness). In Eq. (8), the second term on the RHS describes the friction due to mass addition of the melt layer, the third and fourth terms describe the effects of melt viscosity, the fifth term describes the effect of surface tension, and the last term describes the volumetric Lorentz force. Other parameters are the same as those defined in Refs. [16, 22, and 29].

In this model, we use the following 2D transient heat conduction equation for both liquid and solid phases, where V s is the total velocity of the melt surface, and equals to zero in the solid. Here, the new variable of a moving reference frame has been introduced to replace the old variable z of the fixed reference frame,[22] .

In this work, the numerical calculations were performed using the finite volume method to track the moving interface between solid and liquid phases. T, h, and are discretized on the staggered grid. Here, the origin of the coordinate frame is set at the plate edge. Since is different on the two sides of the separatrix strike point (SSP), the upwind scheme is adopted, respectively, for both sides. In this way, the difficulties associated with conventional approaches can be effectively surmounted. The temperature at the bottom side is fixed to T 0=350 K, representing the cooling effect, the initial temperature of the material is also fixed to T 0, and the magnetic field strength is set to B = 6 T.

3. Results and discussion
3.1. Melt surface evolution

We first look at the evolution of surface topology induced by an ELM heat load of with a Gaussian power density profile (our previous studies on power load on the divertor target of DIII-D,[30,31] EAST,[32,33] and HL-2A[34] by using SOLPS modeling show that the energy flux shape is very similar to Gaussian distribution. To analyze the influence of various factors on the divertor erosion straightly, the Gaussian power density profile was assumed with the same standard deviation on both sides of the SSP). Here a typical value of is taken to be .[16] The initial sample temperature is 1000 K (the steady-state temperature distribution in the PFC block of 1-cm thickness is about 1000 K under the steady-state inter-ELMs heat flux of , which is similar to that of ITER, and the predefined boundary conditions[10,21]). Recent theoretical and experimental works[35,36] indicate that the energy decay length at the SOL (λ q) is inversely proportional to the toroidal magnetic field strength (Bt), and predict for ITER.[36] Considering the flux expansion at the divertor target, it is assumed the energy deposition width during ELMs is 1 cm in the present work.

Figure 2 shows the effect of force on the temporal evolution of the melt surface profile. The result shows that both the surface tension and the force can induce melt motion, but the pressure gradient is negligible when the value of heat flux peak is less than . The W melt flow can lead to a rather deep erosion dent near the region of the SSP, and two humps at the dent edges during the ELM. The heights of the two W humps are the same when the force is neglected. In contrast, if the calculations are performed with force included, the melt W driven by the force can accumulate at the dent edge facing to the force, the forming humps are asymmetric at both the dent edges and the roughness of W plate surface turns more severe, as shown in Fig. 2. The hump height on the facing side increases by than without the force at t=0.8 ms, and at t=1.0 ms, respectively.

Fig. 2. The melt surface profile calculated under an ELM with the peak heat flux of , is . The initial sample temperature is 1000 K.

Table 1 gives the values of the hump height on the facing side at different irradiation time (from 0.75 ms to 1.0 ms in an interval of 0.05 ms) without and with the force. The increase in the amount of the hump height in 0.05 ms is much larger with the force than without the force. It can be seen that the hump height increases more rapidly as the irradiation lasts longer if the force is included. This mainly happens because the longer heat flux duration leads to a deeper melt layer and faster melt motion if the force is included. Garkusha et al.[37] calculated the height of the hump induced only by force is about . Their calculations were performed according to the experimental conditions (B, 1.4 T, maximum external current density, , current duration, 0.15 ms).

Table 1.

The height of peak hump on the facing side, and its increase in an interval of 0.05 ms without and with force.

.

The results obtained from the present work are in good agreement with those in the previous work,[37] considering there is a small difference in magnetic field strength and in current density (in the present paper, the current density varies between , which is similar to ITER parameters[8,38]). All of the results from both experiments and the present modelling indicate that the force makes the small hill on the facing side grow relatively faster, the force may aggravate the W plate erosion, whose extent depends on the temperature and properties of tungsten plates.

Figure 3 shows the temporal evolution of surface temperature. In the initial stage of melting, the difference in the surface temperature between without and with force is not noticeable because the significant melt flow has not yet arrived. After the beginning stage, as the melt flow generated by the surface tension and the force carries heat from the SSP toward its periphery on the facing side, the location of the maximal surface temperature moves a bit from the SSP to the facing side. When the heat flux duration progresses to 1.0 ms, the maximal surface temperature with the force is 1 K higher than without the force, and its location is away from the SSP. From these results, it can be seen that the flow induced by the force results in a higher melt surface temperature on the facing side (see Fig. 3).

Fig. 3. Temporal evolution of the surface temperature. Incident heat flux of , , , the initial sample temperature, 1000 K.

To further assess the effect of force on erosion of W plates, we prolong the calculation time to 1.8 ms. We find that the difference in the melt surface temperature between without and with the force becomes more apparent as the calculation time increases. Although the influence of the force on the melt motion is less than the surface tension (as shown in table 1), the accumulated effect of the force increases remarkably with the calculation time.

Figure 4 shows the temporal evolution of the melt front location in the W plates. At t=0.8 ms, before significant flow of melt W arrives, the melt front locations are not noticeably different between without and with force. At t=1.0 ms, small difference can be seen. Compared with the case without force, the melt depth is a little shallower on the facing side and a little deeper on the other side of the SSP. The melt layer flow carries heat from the central part of the dent, increasing the surface temperature nearby, which siphons a part of the total heat flux penetrating into the bulk. The difference in the melt motion between the cases without and with force gives rise to the difference in the conduction flux penetrating the melt layer. When the duration time increases to 1.8 ms, the difference in the melt front location between without and with force is more evident.

Fig. 4. Temporal evolution of the melt boundary under peak heat flux of , and is . The black dot-dashed line is for the case with force.
3.2. dependence

Impurities in the W plates will change the surface tension.[16] In Eq. (8), is the gradient of surface tension coefficient with regarding to temperature, and the term in the RHS describes the effect of surface tension. Therefore it is necessary to evaluate the influence of different values of (the values of are taken from Ref. [16]) on the erosion of W plates when the force cannot be ignored. Figure 5 presents the melt surface profiles under different values of for a scenario of a sample with initial temperature of 1000 K, a heat flux of , and a heat flux duration of 0.8 ms. We can see that the hump height on the facing side of the SSP is larger for different values of . It is found that the humps become more asymmetric when the value of is smaller. The peak value of the hump height is higher than in the case without the force (increased by 87.2%) for of , and is higher (increased by 3.3%) for of . When the value of is small, the contribution driven by the surface tension drops in the total of driving force while the contribution of force relatively increases. The calculation results show that the presence of impurities in the melt layer which changes the value of surface tension will accordingly affect the effects of surface tension and force on the melt layer.

Fig. 5. Calculated melt surface profile under peak heat flux of for different coefficients of surface tension. The heat load duration is 0.8 ms.
3.3. Heat flux dependence

The heat deposition Q on the divertor is anticipated in the range about on a time of 0.1 ms–1 ms for ELMs for ITER. Figure 6 shows temporal evolution of the W melt surface profiles without and with the force under different peak heat fluxes of 2000, 2500, and , respectively. For them, the value of is ,[16] and the heat load duration is 0.8 ms. Although the chosen peak fluxes in the present work are in the ITER range, they are relatively moderate. The humps at both edges of the dent caused by such an ELM heat load with a Gaussian power density profile are the same at height when the force is neglected. In contrast, when the force is considered, the humps at both the dent edges grow at different rates—the hump on the facing side grows noticeably faster. When the force is taken into account, the maximal hump height increases rapidly from less than under heat flux of to about under heat flux of at t=0.8 ms. This mainly happens because as the heat flux goes up, the force increases, and the asymmetry becomes larger. In this work, only the Lorentz force caused by external ionic currents is considered. If the eddy and halo currents are also considered, then the force will have an even larger effect on the melt motion.

Fig. 6. W melt surface profile without and with force under different peak heat fluxes ( , , and heat load duration, 0.8 ms).
3.4. Surface roughness

Melt motion under the force can induce some surface instability, which may result in the splashing of W observed in TEXTOR,[17,39] causing significant local and core plasma contamination. The melt motion will also increase the hump on the facing side, where the erosion will be increased and the electric field will be enhanced, thus the heat flux and electron emission will increase, resulting in initiating the arc.[7] The surface topography gives the change of target surface under the force visually while the surface roughness introduced here can be calculated quantitatively to describe the change of target surface. In the following, we further assess the effect of force on the W divertor erosion quantitatively by evaluating the surface roughness. Surface roughness[40,41] has been defined as here h(i) is the actual height of each point on the plate surface (a total of 150 points with an fixed interval of have been taken between x=0.35 cm and 0.65 cm in the calculation), h is the average height of plate surface. The surface roughness varies with heat flux, the value of and duration, as shown in Fig. 7. We can see that the surface roughness is about zero when the incident heat flux is less than about , because the divertor plates just begin to melt under heat flux of , as shown in Fig. 7(a). As the heat flux goes up, the surface roughness becomes larger, and the difference in surface roughness becomes even larger between without and with the force at the same . From Fig. 7(b) we can see that as the absolute value of goes up, the target surface roughness becomes more severe, and the difference in surface roughness becomes smaller between without and with the force under the same heat flux. The bigger absolute value of means that there is a bigger surface tension gradient, which gradually overshadows the effect of force. The longer duration of the ELM, means that the surface roughness grows larger, and the difference in surface roughness between without and with force becomes more evident, as shown in Fig. 7(c).

Fig. 7. Surface roughness of W plates versus heat flux (a), value of (b), irradiation duration (c). The initial sample temperature is 1000 K.
4. Conclusions

A two-dimensional (2D) fluid dynamics model has been employed in this work by solving the liquid hydrodynamic Navier-Stokes equation together with the 2D heat conduction equation to assess the thermal erosion and topographical change of divertor target plates during a type-I like ELM in ITER with a Gaussian power density profile heat load. Apart from the driving forces, such as surface tension and pressure gradient, this paper particularly addresses the force and its effect on the melt layer. The present simulation results indicate: (i) the force makes the humps grow with different rates at the crater edges, and the difference in topographical modification between without and with the force becomes noticeable when the W plates is exposed to a heat flux higher than for 0.8 ms with equal to , (ii) the small absolute value of can enhance the effect of force on the melt motion due to the reduction of the surface tension gradient, (iii) the force makes the surface roughness increase faster as the heat flux increases, which can make the erosion of the divertor target more inhomogeneous, (iv) the melt motion induced by the force during a type-I like ELM in ITER will increase the target erosion. These can all have a great impact on the lifetime of the divertor target. Moreover, the eroded and splashing W impurity may contaminate the core plasma and even cause the disruption, which is a big threat to the ITER operation.

Acknowledgment

This work was supported by the cluster of PSI & APD group (http://plasmaphys.dlut.edu.cn) and Supercomputing Center of Dalian University of Technology.

Reference
[1] Qian X Y Peng X B Wang L Song Y T Ye M Y Zhang J W Li W X Zhu C C 2016 Nucl. Fusion 56 026010
[2] Gao J M Liu Y Li W Cui Z Y Zhou Y Huang Y A Ji X Q 2010 Chin. Phys. 19 115201
[3] Lingertat J Tabasso A Ali-Arshad S Alper B van Belle P Borrass K Clement S Coad J P Monk R 1997 J. Nucl. Mater. 241�?43 402
[4] Xue L Duan X R Zheng G Y Liu Y Q Yan S L Dokuka V N Lukash V E Khayrutdinov R R 2016 Chin. Phys. Lett. 33 055201
[5] Qiu Q L Xiao B J Guo Y Liu L Wang Y H 2017 Chin. Phys. 26 065205
[6] Pitts R A Bardinb S Bazylevc B van den Berg M A Bunting P Carpentier-Chouchana S Coenenf J W Corre Y Dejarnac R Escourbiac F Gaspar J Gunng J P Hirai T Hong S H Horacek J Iglesias D Kommh M Krieger K Lasnier C Matthews G F Morgan T W Panayotis S Pestchanyi S Podolnik A Nygren R E Rudakov D L De Temmermana G Vondracek P Watkins J G 2017 Nucl. Mater. Energ. 12 60
[7] Federici G Skinner C H Brooks J N Coad J P Grisolia C Haasz A A Hassanein A Philipps V Pitcher C S Roth J Wampler W R Whyte D G 2001 Nucl. Fusion 41 1967
[8] Sergienko G Bazylev B Huber A Kreter A Litnovsky A Rubel M Philipps V Pospieszczyk A Mertens Ph Samm U Schweer B Schmitz O Tokar M Team Textor 2007 J. Nucl. Mater 363�?65 96
[9] Federici G Andrew P Barabaschi P Brooks J Doerner R Geier A Herrmann A Janeschitz G Krieger K Kukushkin A 2003 J. Nucl. Mater. 313�?16 11
[10] Federici G Loarte A Strohmayer G 2003 Plasma Phys. Control. Fusion 45 1523
[11] Raffray A R Federici G 1997 J. Nucl. Mater. 244 85
[12] Federici G Raffray A R 1997 J. Nucl. Mater. 244 101
[13] Hassanein A Konkashbaev I 2000 Fusion Eng. 51�?2 681
[14] Sizyuk V Hassanein A 2015 Phys. Plasmas 22 013301
[15] Wurz H Bazylev B Landman I Pestchanyi S Gross S 2001 Fusion Eng. 56�?7 397
[16] Bazylev B Wuerz H 2002 J. Nucl. Mater. 307�?11 69
[17] Coenen J W Bazylev B Brezinsek S Philipps V Hirai T Kreter A Linke J Sergienko G Pospieszczyk A Tanabe T Ueda Y Samm U R-Team T E X T O 2011 J. Nucl. Mater. 415 S78
[18] Bazylev B N Janeschitz G Landman I S Pestchanyi S E 2005 Fusion Eng. 75�?9 407
[19] Bazylev B N Janeschitz G Landman I S Loarte A Pestchanyi S E 2007 J. Nucl. Mater. 363�?65 1011
[20] Igitkhanov Y Bazylev B 2014 IEEE Trans. Plasma Sci. 42 2284
[21] Huang Y Sun J Z Sang C F Ding F Wang D Z 2014 Acta Phys. Sin. 63 035204 (in Chinese) https://doi:10.7498/aps.63.035204
[22] Huang Y Sun J Z Hu W P Sang C F Wang D Z 2016 Fusion Eng. 102 28
[23] Miloshevsky G V Hassanein A 2010 Nucl. Fusion 50 115005
[24] Loarte A Saibene G Sartori R Campbell D Becoulet M Horton L Eich T Herrmann A Matthews G Asakura N Chankin A Leonard A Porter G Federici G Janeschitz G Shimada M Sugihara M 2003 Plasma Phys. Control. Fusion 45 1549
[25] Jiang C B Zhang Y L Ding Z P 2007 Comput. Fluid Mech. Beijing China Power Press 211 in Chinese
[26] Carslaw H W Jaeger J C 1959 Conduction of Heat in Solids 2 Clarendon Oxford 89 91
[27] Behrisch R 2010 J. Surf. Invest-X-Ray + 4 549
[28] Udaykumar H S Shyy W 1995 Int. J. Heat Mass Tran. 38 2057
[29] Semak V V Damkroger B Kempka S 1999 J. Phys. D: Appl. Phys. 32 1819
[30] Sang C F Guo H Y Stangeby P C Lao L Taylor T S 2017 Nucl. Fusion 57 056043
[31] Sang C F Stangeby P C Guo H Y Leonard A W Covele B Lao L L Moser A L Thomas D M 2017 Plasma Phys. Control. Fusion 59 025009
[32] Sang C F Du H L Zuo G Z Bonnin X Sun J Z Wang L Wang D Z 2016 Nucl. Fusion 56 106018
[33] Sang C F Ding R Bonnin X Wang L Wang D Z Team E A S T 2018 Phys. Plasma 25 072511
[34] Sang C F Wang Z H Min Xu Wang Q Wang D Z 2018 Fusion Eng. 136 1041
[35] Goldston R J 2012 Nucl. Fusion 52 013009
[36] Eich T Leonard A W Pitts R A Fundamenski W Goldston R J Gray T K Herrmann A Kirk A Kallenbach A Kardaun O Kukushkin A S LaBombard B Maingi R Makowski M A Scarabosio A Sieglin B Terry J Thornton A A S D E X Upgrade Team and Contributorsa A 2013 Nucl. Fusion 53 093031
[37] Garkusha I E Bazylev B N Bandura A N Byrka O V Chebotarev V V Landman I S Kulik N V Makhlaj V A Petrov Yu V Solyakov D G Tereshin V I 2007 J. Nucl. Mater. 363�?65 1021
[38] Bazyle B N Janeschitz G Landman I S Pestchanyi S E 2005 J. Nucl. Mater. 337�?39 766
[39] Coenen J W Philipps V Brezinsek S Bazylev B Kreter A Hirai T Laengner M Tanabe T Ueda Y Samm U TEXTOR Team 2011 Nucl. Fusion 51 083008
[40] Elsholz F Scholl E Scharfenorth C Seewald G Eichler H J Rosenfeld A 2005 J. Appl. Phys. 98 103516
[41] Elsholz F Scholl E Rosenfeld A 2004 Appl. Phys. Lett. 84 4167