Investigation of heat flux deposition on divertor target on the Large Helical Device with EMC3-EIRENE modelling

The measured divertor heat flux profiles are compared to the EMC3-EIRENE simulations for two different times of an LHD discharge, corresponding to higher and lower edge temperatures. The relation between the three-dimensional magnetic field structure and the heat flux distributions on the divertor has been analysed. The modelled heat flux for the lower plasma temperature case has a better agreement with the experimental result obtained by the Langmuir probes, which shows a qualitative reproduction of the experimental profile shape. However, the heat flux distribution for the high plasma temperature case shows a different behaviour between the simulation results and the experimental measurements. The detailed analysis of the heat flux distribution for the higher temperature case which has a larger discrepancy has been performed, both quantitatively and qualitatively. The radiation of the eroded impurity from divertor target plates has a minor effect on the heat flux distribution. Non-uniform cross-field transport coefficients are used in the simulations and its impact on the heat flux distributions is discussed for the case of the high plasma temperature.


Introduction
The investigation of the heat flux distribution on the divertor target is considered as one of the most critical issues in fusion devices, which has a strong implication for the divertor performance [1]. For the tokamaks with axisymmetric structure, the heat flux is primarily transported along the parallel magnetic field lines in the scrape-off layer to the divertor. The studies about the edge plasma behaviour and heat flux deposition on the divertor target have been performed by SOLPS code for ASDEX Upgrade [2,3] and DIII-D [4].
While for the stellarators and tokamaks with resonant magnetic perturbation (RMP) application, the heat transport is more complicated due to the existence of the magnetic field stochasticity in the edge plasma [5][6][7][8][9][10][11][12][13]. In particular, the effect of the perpendicular heat transport becomes more remarkable due to the long connection length in the stochastic magnetic field configuration. Hence, the distribution of the heat flux deposition on the divertor target is determined by the three-dimensional (3D) magnetic configuration and the balance of parallel and perpendicular transport. Moreover, the 3D footprint of the heat flux would have an impact on the erosion/deposition balance on the divertor target [14]. Therefore, studies of heat flux distribution in 3D magnetic field are important to obtain a better understanding of the underlying mechanisms of edge plasma transport and the lifetime of divertor target.
In the LHD experiments, the heat flux deposition on the divertor target has been investigated by using the Langmuir probes, thermocouples and infrared camera [15,16]. It is found that the heat flux distributions are strongly associated with the 3D magnetic configuration and operational regime. The 3D modelling of the heat flux deposition on the divertor target of LHD has been performed by the edge transport code EMC3-EIRENE [17,18]. In the modelling, it is revealed that the cross-field transport has an impact on the heat flux distribution on the divertor target [15,16]. However, there are still some uncertainties remaining due to the technical difficulties in the experiments and simulations. In this study, the electron density and temperature profiles at the mid-plane obtained by the Thomson scattering system in LHD are used to determine the cross-field transport coefficients of background plasma used by EMC3-EIRENE modelling. On the other hand, the EMC3-EIRENE code has been upgraded to be able to deal with spatially non-uniform cross-field transport coefficients [19] and the new code version has been benchmarked with analytical models [20]. Furthermore, the development of the computational grid for the divertor leg regions of LHD has been achieved [21], which results in a good description of the plasma transport behaviour in the divertor leg regions connecting to the divertor target. Based on these improvements, the investigations of the heat flux distributions on the divertor target in LHD experiments have been motivated and revisited by the EMC3-EIRENE modelling. The EMC3-EIRENE code has a very good flexibility for treating arbitrary magnetic geometry such as helical devices and nonaxisymmetric tokamaks with RMP fields, which has been widely used for the 3D edge plasma modelling [22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38].
In section 2, the magnetic field structure of LHD is introduced. In section 3, the results of EMC3-EIRENE modelling are presented and benchmarked against the experimental results. The impact of non-uniform cross-field transport for different scenarios is studied. The discussion on the relation between the parallel and perpendicular heat transport is attempted for the high plasma temperature case. Finally, the results are summarized in section 5. The magnetic field geometry of LHD is created by a pair of helical coils twisted with poloidal and toroidal pitch numbers of l=2 and n=10, respectively. Thus, the magnetic field structure has 10 field periods in the toroidal direction. The plasma distributions at the toroidal angles of j =    0 , 18 and 36 in LHD have the up-down symmetry structures. Figure 1(b) shows the poloidal cross section of LHD at the toroidal angle of j =  18 . The stochastic magnetic field geometry is generated by the helical coil system due to the overlapping of magnetic island chains. Two X-points are created by the helical coils of LHD as shown in figure 1(b), and these X-points rotate poloidally when they move toroidally. Beyond these X-points, four divertor legs induced by the stretch of the flux tubes are connected to the divertor target plates. In order to facilitate discussion, the divertor target plates as shown in figure 1(b) are denominated as Open 1-4, respectively. The entire edge region outside the last closed flux surface (LCFS) is termed as the stochastic layer in this study. The vacuum vessel wall and divertor target plates in LHD are made of stainless steel and graphite, respectively. The magnetic field structure of the divertor leg is wider at the inboard side compared to the outboard side for j =  18 in figure 1(b). This configuration leads to in-out asymmetric distributions of background plasma along the divertor legs, and thus the recycling flux at the inboard and outboard sides.

Setup of EMC3-EIRENE simulations
For simulation of the edge magnetic configuration of LHD, the whole torus is reduced to one field period, i.e., j=0°-36°(=360°/10) with the assumption of a complete periodicity in the toroidal direction. Further, the simulation domain can be reduced to a toroidal range of j=0°-18°, i.e., onehalf of a field period, assuming that another one-half field period is reproduced by stellarator symmetry (=36°/2). Therefore, the physical quantities in the domain of j = 18°-36°can be obtained with the mapping of , where f represents physical quantity in the domain of j=0°-18°, and j ¢ ¢ ¢ ¢ f R Z , , , represent physical quantity and coordinates in the domain of j=18°-36°. Here, an up-down symmetry is imposed at the boundaries of j =    0 , 18 and 36 in order to ensure consistency of the mapping. The computational mesh has been divided into two blocks of j =     --0 9 and 9 18 in order to ensure the up-down symmetry property at the boundaries of j =   0 and 18 . The reversible field line mapping technique is used for tracing individual MC particle trajectories at the toroidal block boundary of j =  9 [39,40]. The detailed information of the computational grid for LHD is introduced in [21].
The input parameters for EMC3-EIRENE modelling are specified according to the experimental measurements during LHD exposure. Two scenarios for different times during the discharge (shot #65769) are simulated in this work, which represent the cases of the high and low edge plasma temperatures at the LCFS. The input power is 10 MW and the upstream density n LCFS is set to 1.8×10 19 m −3 for the high plasma temperature case. For the low plasma temperature case, the input power is 4 MW and the upstream density n LCFS is set to 2.2×10 19 m −3 . The decay lengths of 3, 5 and 30 cm are set at the outer radial boundaries for particle, energy and momentum transport, respectively. The toroidal magnetic field is 2.75 T and the magnetic axis is located at R ax =3.6 m. The Bohm sheath boundary condition is used at the divertor target plates. The ion and electron sheath heat transmission coefficients are assumed to be 2.5 and 4.5, respectively. The sputtering coefficient is assumed to be 0.005 in the simulation. The erosion distribution of carbon impurity is determined by the incident plasma flux and sputtering coefficient in EMC3-EIRENE. At the present version of EMC3, the sputtering coefficient is a free parameter fixed in the simulation, hence the distribution of the released carbon impurity is proportional to the particle flux deposition pattern. The cross-field particle and energy transport coefficients of background plasma,D and c^, are determined by fitting the electron density and temperature profiles measured by the Thomson scattering system in LHD. Figure 2 shows the plasma temperature and density distributions measured at the midplane of LHD and modelled by EMC3-EIRENE code. A good matching between the experiment and modelling is obtained by usingD =0.6 m 2 s −1 and c^=0.6 m 2 s −1 for the high plasma temperature case and usingD =0.4 m 2 s −1 and c^=0.2 m 2 s −1 for the low plasma temperature case, respectively. Here, it should be noted that the cross-field transport coefficients of background plasma are constant in space. The impact of the non-uniform cross-field transport will be analysed below. The above-mentioned parameters are used as the default values for the following EMC3-EIRENE simulations.
3.2. The heat flux distribution with the constant cross-field transport coefficients Figure 3 shows the 2D distributions of the heat flux deposition on the divertor target plates for the high and low plasma temperatures, respectively. The definition of the poloidal angle in LHD starts from the inboard midplane and increases to upper divertor. The black cross points in figure 3 show the start and end positions of the Langmuir probe array. The toroidal range of the Langmuir probe array is from  17.0 to  15.6 . The label '0 mm' at  17.0 indicates the original point of the horizontal axis in figure 4. The red arrow indicates the direction of the embedded Langmuir probes, which is also the positive direction of the horizontal axis in figure 4. The main heat flux is deposited at the divertor target plates (Open 2 and 3) for the high and low plasma temperatures, which is the inboard side of the torus as shown in figure 1(b). This is because more flux tubes are connected to the inboard divertor target as mentioned above, which act as the heat transport channels in the stochastic layer. The input power for the high plasma temperature is 10 MW, while it is 4 MW for the low plasma temperature. This leads to a higher heat flux deposition on the divertor target plates of Open 2 and 3 for the high plasma temperature as shown in figure 3. The power lost through the outermost boundary is usually small in the simulation, which is guaranteed by properly choosing the domain size in grid construction.
The respective integrated powers deposition on the divertor are about 6.8 and 2.6 MW for the high and low plasma temperature cases. For the high plasma temperature, the integrated power is about 5.1 MW on the inboard divertor   fluxes on the in-and out-board divertor targets are about 6.7 and 4.7 MW m −2 for the high plasma temperature; while the respective peak heat fluxes on the in-and out-board divertor targets are around 3.4 and 1.9 MW m −2 for the low plasma temperature. In the modelling, the hydrogen radiation powers are about 2.3% and 3.6% of the individual input powers for the high and low plasma temperature cases, respectively. The lost powers by the plasma and neutral collisions are about 29.5% and 31.9% of the individual input powers for the high and low plasma temperature cases, respectively. Figure 4 displays the experimental and simulated heat flux and connection length distributions for the high and low plasma temperature cases in the modelling. The horizontal axis is along the line where the Langmuir probe array is embedded on the inboard divertor plate (Open 2). It can be seen that the peak values of heat flux are about 6.0 and 3.0 MW m −2 for the high and low plasma temperatures, respectively. There are three locations at about 50, 70, and 100 mm where there are larger connection lengths compared to the other neighbour positions. The larger connection length leads to a longer particle transit time and further a higher plasma density compared to the shorter connection length, which results in a higher heat flux at the long connection length as shown in figure 4. The red lines in figure 4 are the experimental data of the heat flux distribution for the high and low plasma temperatures, which are measured by the Langmuir probe array at the divertor target of Open 2. For the low plasma temperature, the simulated and experimental distributions show a better agreement compared to the high plasma temperature. For the high plasma temperature, the simulated peak value of the heat flux located at 70 mm is 3 times higher than the experimental result. Further, the heat flux in the modelling at the positions between 80 and 100 mm is much higher than that in the experiment in figure 4(a).
The investigation of the intrinsic carbon impurity eroded from the divertor target plates is carried out in order to assess the impact of the impurity radiative losses on the heat flux distribution for the high plasma temperature case. Figure 5 shows the heat flux distributions for different sputtering coefficients for the high plasma temperature. It is seen that the heat flux distributions do not change much for different sputtering coefficients. The radiation power of the eroded impurity increases 10 times when the sputtering coefficient increases from 0.5% to 5.0% for the high plasma temperature. However, the impurity radiation power is about 4.0% of the input power for the sputtering coefficient of 5.0%, which is so small that the heat flux distribution does not change too much as shown in figure 5. It can be seen that even for a high sputtering coefficient of 5.0%, the concentration of carbon impurity in the plasma is still not high enough to radiate the power and to effectively mitigate the heat flux in the region with a large discrepancy. In the following subsection, the large difference in the heat flux between the simulation and experiment for the high plasma temperature is studied based on the impact of the non-uniform transport.
3.3. The impact of the non-uniform cross-field transport on the heat flux distribution 3.3.1. Zone-dependent non-uniform cross-field transport.
Since the magnetic shear inside the divertor leg region of LHD is strong, the computational mesh for plasma simulation is divided into an edge region and four divertor leg regions at each toroidal angle [21]. Hence, the cross-field transport coefficients in the edge region and divertor leg regions can be specified separately in the modelling. Figure 6 shows the experimental and simulated heat flux distributions on the inboard divertor plate with different non-uniform cross-field coefficients c^in the divertor leg regions for the high plasma temperature. Here, the cross-field coefficients in the edge region is maintained same as the above constant cross-field coefficients (c^=0.6 m 2 s −1 ) in order to guarantee that the plasma parameters have a good agreement with the Thomson scattering data at the mid-plane. In addition, the test  simulations are carried out to check the effect of the crossfield coefficient of particle transportD on the heat flux distribution. However, the influence ofD is small compared to c^because the heat conduction is dominant for the high plasma temperature case. Hence, we focus on the studies on the variation of c^in the following simulations. In figure 6, it is seen that the heat flux distribution has a small change for c =6.0 m 2 s −1 compared to c^=0.6 m 2 s −1 . When the cî ncreases to 60.0 m 2 s −1 , the heat flux reduces obviously between 70 and 100 mm, almost by a factor of about 2. However, the peak value of heat flux at 70 mm is still two times higher than the experimental result. The heat flux between 40 and 60 mm also decreases and becomes even lower than the experimental result.

B-dependent non-uniform cross-field transport.
On the other hand, the dependence of the cross-field transport coefficient on the magnetic field strength is also studied. cross-field transport coefficient, which is dependent on the magnetic field strength. Here, c j and B j represent the crossfield transport coefficient and magnetic field strength for each physical cell of EMC3-EIRENE modelling, and B is the average value of B j at each toroidal angle. The exponent of f is chosen to be    1, 2 and 3 to study the impact of B-field dependent cross-field transport coefficient. It should be noted that the fixed value of c is changed for different f to make sure that the plasma parameters have a good agreement with the Thomson scattering data at the mid-plane. Figure 7(b) presents the non-uniform cross-field transport coefficient c j for the case of f=3 and c=2.0 m 2 s −1 at the poloidal cross section of j =  18 . It is seen that the values of c j are larger at the divertor leg regions compared to the edge region. Figure 8 shows the experimental and simulated heat flux distributions on the inboard divertor plate with different factor f for the high plasma temperature. The peak values of the simulated heat flux at 70 and 90 mm reduce slightly when the factor f increases from 1 to 3 in figure 8(a). This is due to the enhancement of the perpendicular transport coefficients in the divertor leg regions as shown in figure 7(b). The tendency is similar for changing the factor f from −3 to −1 as shown in figure 8(b). The peak values of the heat flux at 70 and 90 mm decrease slightly due to the increase of the cross-field transport coefficients. However, the distributions of heat flux for different cases still have a large difference compared to the experimental result for the case of the high plasma temperature.

T-dependent non-uniform cross-field transport.
The dependence of the cross-field transport coefficient on the plasma temperature is studied in figure 9, which presents the experimental and simulated heat flux distributions on the inboard divertor plate for the high plasma temperature case. The above result for the uniform cross-field coefficient (c =0.6 m 2 s −1 ) is also shown for comparison. The formula of c c = ( ) j T T f j is used to calculate the non-uniform cross-field transport coefficient, where T j is the plasma temperature for each physical cell and T is the average value of T j at each toroidal angle. Here, the result shown in figure 9 is just for the case of f=1 and c=0.2 m 2 s −1 . The scans for different factors f are also performed, but the simulated plasma parameters have no agreement with the Thomson scattering data at the mid-plane. It can be seen that the simulated heat flux for the non-uniform transport is still much larger than the experimental result in figure 9. Moreover, the simulation result is even higher than the case with the uniform transport coefficients since the small cross-field coefficients are used in the divertor leg regions for the non-uniform transport.

Discussion
Based on the above simulations, it is indicated that the parallel heat transport is stronger compared to the perpendicular heat transport for the high plasma temperature case. In order to assess the perpendicular scale of the heat transport on the divertor target plate, an analytic study of the perpendicular transport scale has been performed during the parallel transport time scale based on the available simulation data. This can quantify the scale of cross-field transport needed to be to compete with the parallel transport and to significantly reduce where L c is connection length along the magnetic field line. The flux tube tracing method from the divertor target to the upstream is used to obtain the upstream plasma parameters by Kmag code [36,41]. The flux tube tracing of L c =10.0 m from the line along the Langmuir probe position to the upstream has been conducted to make sure more than one poloidal turn for all flux tubes. According to the upstream plasma temperature and density, the ratio of the parallel conductive heat transport to the perpendicular conductive heat transport 0.5 where k e i , is heat conductivity coefficient [42]. The electron parallel conductive heat transport is used to approximate the parallel heat transport due to the high heat conductivity coefficient. Based on this ratio, one can estimate the perpendicular energy transport scale as shown in figure 10, which presents the distributions of the heat flux and the perpendicular energy transport length along the line where the Langmuir probe array is embedded for the high plasma temperature. Here, the default value of c^=0.6 m 2 s −1 is used for the calculation of the perpendicular energy transport length. It can be seen that the magnitude of the perpendicular energy transport scale is about 1 mm for the high heat flux region between 60 and 100 mm in figure 10. Especially for the peak flux position, the perpendicular energy transport length is about 0.5 mm.   Here we briefly revisit the analyses in the preceding sections to see the impact of the perpendicular heat transport on the simulated heat flux distribution. For the zone-dependent non-uniform cross-field transport, the perpendicular transport coefficients increase 10 and 100 times compared to the default value of c^=0.6 m 2 s −1 at the divertor leg regions. The perpendicular transport scales increase by a factor of about 3 and 10 for c^=6.0 and 60.0 m 2 s −1 , respectively. According to the results in figure 10, one can estimate that the individual perpendicular energy transport lengths are 3 and 10 mm at most due to shorter parallel transport length at the divertor leg regions (<10.0 m). This leads to a slight change of the heat flux distribution for c =6.0 m 2 s −1 , while a relatively large variation for c =60.0 m 2 s −1 . For the B-dependent non-uniform cross-field transport, the maximum perpendicular transport is about c =5.0 m 2 s −1 by using the factor f=3 for the inboard divertor regions as shown in figure 7(b). This results in a similar change of the heat flux distribution in figure 8(a) as that in figure 6 for the case of c^=6.0 m 2 s −1 . Hence, the change of the heat flux distribution is also marginal. In addition, the use of the negative factor f from −1 to −3 leads to a suppression of the perpendicular energy transport in the divertor leg regions, which even results in a higher peak value of heat flux as shown in figure 8(b). For the T-dependent nonuniform cross-field transport, the tendency of the simulation results in figure 9 are similar to the results in figure 8(b) due to the lower perpendicular energy transport coefficients used in the divertor leg regions, which leads to a larger discrepancy between the simulation and experiment in figure 9.

Summary
The heat flux distributions on the divertor target plates have been investigated with the 3D edge transport code EMC3-EIRENE in comparison with the experimental measurements obtained by the Langmuir probes in LHD. The distribution of 3D magnetic field structure plays an important role in the edge heat transport in the stochastic layer, which leads to a higher heat flux deposition on the inboard divertor target plates. The modelled heat flux distribution shows a better agreement with the experimental result at the low plasma temperature. However, the heat flux distribution at the high plasma temperature has a different behaviour between the simulation results and the experimental measurements.
The enhanced sputtering coefficient leads to a higher impurity radiation, but this has a negligible effect on the heat flux distribution for the high plasma temperature. A detailed analysis of the non-uniform cross-field transport coefficients has been performed to study the impact of the balance of parallel and perpendicular heat transport on the heat flux deposition on the divertor target. The strong increase of the cross-field transport in the divertor leg regions can reduce the heat flux, but the simulated distribution is still different from the experimental results. The dependence of the non-uniform cross-field transport coefficient on the magnetic field structure and plasma temperature is investigated. However, the change of the heat flux distribution for different scenarios is minor for the high plasma temperature.
For the present version of the EMC3 code, parallel transport is purely classical without kinetic corrections. In view of the complex 3D topology, where flux tubes with completely different L c strongly interact with each other via perpendicular transport, it is difficult to make a reasonable estimation of kinetic effects on the heat flux deposition analysis. In addition, the power loss through the last radial surface rises with increasing the cross-field transport coefficients for fixed temperature decay lengths as boundary conditions. However, even for the high cross-field transport coefficient c =60 m 2 s −1 , the power loss is not high enough to affect the simulation results.
There are some possible reasons leading to the discrepancy between the modelling results and experimental measurements for the high plasma temperature case. First, the equilibrium magnetic field structure can be changed in the experiments due to the local current flowing somewhere in the edge regions. However, the unperturbed magnetic field structure is used for the current modelling. Second, the heat load distributions can also be affected due to the energetic particles, which are not taken into account in the fluid model. The impacts of the magnetic perturbations and the energetic particles can be studied by the magnetohydrodynamic models. The corresponding investigations are ongoing to provide more accurate information. The present analysis should be revisited in the future to reveal the underlying mechanisms causing the discrepancy between simulations and experiments at the high plasma temperature.