Analysis and design of fast flow liquid Li divertor for fusion nuclear science facility (FNSF) using coupled plasma boundary and LM MHD/heat transfer codes^*

Ne seeding is used to maintain the divertor heat flux lower than 10 MW m⁻² via heat dissipation via radiation. The D₂ puffing level is maintained at a constant $8\times 10^{21}$ atoms s⁻¹ while the Ne seeding condition is varied over a wide range. The simulation results presented in this subsection and section 3.2 with standalone SOLPS-ITER (i.e. without LM MHD/heat transfer code coupling), as the goal of Ne seeding is to first obtain fluxes and baseline conditions which meet the FNSF design requirements listed in table 1.

Figure 4(a) depicts the relationship between the upstream electron density and peak outer divertor plasma heat flux. As Ne seeding levels increase from 0 to $5\times 10^{20}$ atoms s⁻¹, there is a significant reduction in the peak divertor energy flux. Specifically, Ne seeding levels of 5$\times 10^{20}$ atoms s⁻¹ results in a peak outer divertor plasma energy flux of 1.75 MW m⁻². However, it is important to note that Ne seeding reduces the upstream electron and main ion density as well, due to the presence of highly radiating charge states of neon. For the highest Ne seeding level, the upstream electron density is reduced by nearly 50$\%$.

Figure 4(b) shows the dependence of the electron cooling rate on peak divertor heat flux for the divertor and core + SOL, with separate data point shapes to discriminate between the power radiated within the divertor volume and power radiated within the core + SOL volumes. The region definition used in the SOLPS simulation is depicted in figure 1, where blue, orange, yellow, and purple regions correspond to the core, SOL, inner divertor, and outer divertor, respectively. The electron cooling due to Ne dissipation processes (ionization, neutral radiation, line radiation, and bremsstrahlung radiation) is more effective in the divertor (facilitated by the seeding location). At elevated seeding levels, Ne also induces significant radiation in the core and SOL regions, contributing to a decrease in core plasma performance.

Figure 4(c) illustrates the dependency of the effective ion charge (Z $_\textrm{eff}$ = ($\sum_\alpha Z_\alpha^2 n_\alpha$)/n_e ) at the OMP and the entire SOL volume on the outer divertor peak heat flux. Ne concentration is increased as the Ne seeding rate is increased, as illustrated in figure 4(d). While Ne seeding effectively reduces the $q_{\perp,\textrm{max}}^\textrm{Odiv}$ to below 10 MW m⁻² through line radiation, this comes at the expense of an unacceptable increase in Z$_\textrm{eff}$.

In conclusion, simulations using only Ne starting from a D₂ puff level that gives $n_{e,\textrm{sep}}^\textrm{OMP}$ = 1× 10²⁰ m⁻³ fail to meet the FNSF requirements for upstream electron density and impurity concentration due to elevated core radiation and dilution. To address this, both the D₂ and Ne seeding levels are thoroughly scanned to identify an operating window. The D₂ puffing level is increased to counteract the core density drop, while a smaller amount of Ne is used to refine the operating parameters for FNSF, as discussed in the next subsection.

SOLPS is employed to identify the puffing and pumping levels that yield acceptable electron density at the OMP separatrix and the peak divertor energy flux. Through an extensive set of scans, it has been determined that a combination of Ne seeding (1$\times 10^{20}$ atoms s⁻¹) and D₂ puffing (4$\times 10^{22}$ atoms s⁻¹) meets all the FNSF demands on the divertor heat flux (refer to figure 5) and upstream density (refer to figure 6). Ne concentration at the OMP separatrix is found to be around 0.3$\%$. Ne seeding effectively mitigates plasma heat flux to the divertor through radiation while D₂ puffing counteracts the reduction in upstream density. More details can be found in [27].

Zoom In Zoom Out Reset image size

Figure 5 shows the radial profile of the total heat fluxes on the outer and inner divertor. Ne seeding significantly mitigates the plasma heat flux (from ∼25 to ∼4 MW m⁻²) to the outer divertor, resulting in the plasma, neutral, and radiation loads of similar magnitudes. In this work, the radiation heat load was computed in cylindrical approximation using the SOLPS internal tool (WLLD), however, the radiative heat load can also be computed with external tools, as referenced in [44, 45].

The peak energy flux on the inner and outer divertor is found to be lower than 10 MW m⁻² (refer to figure 5) while figure 6 shows the radial profile of electron density and $Z_\textrm{eff}$ at the OMP, with $n_{e,\textrm{sep}}^\textrm{OMP}$ found to be nearly $1 \times 10^{20}$ m⁻³ and $Z_\textrm{eff}$ also lower than 2.2, meeting all the requirements listed in table 1. A well-balanced combination of puffing is necessary to maintain plasma parameters in the upstream and in the divertor region. The combined Ne and D₂ puffing regime, identified through these scans, is utilized for coupling and detailed analysis involving liquid lithium, as discussed in section 4. The LM surface temperature is calculated based on the deposited energy flux (figure 5).

LM PFC components introduce new sources of plasma impurities, including physical sputtering, evaporation, and sputtering of loosely bound ad-atoms [3]. The emitted fluxes from an LM surface are calculated based on the equations and parameters outlined in [46]. Figure 7 shows the emitted Li fluxes from an LM surface as a function of T $_\textrm{surf}$. The emitted Li fluxes are calculated according to equation (4). Physical sputtering as considered here is independent of T $_\textrm{surf}$. In contrast, evaporation loss rapidly increases with the T $_\textrm{surf}$, while the ad-atom loss rate also rises with T $_\textrm{surf}$ until it largely saturates at temperatures ($\gt$500 ^∘C). The basic assumptions and simplifications of the work are as follows: (1) no chemical interactions of D and Ne with the LM surface, (2) pure Li emitted from the surface depending on the local plasma conditions, (3) D and Ne are fully recycled from the LM surface (Li is already 100$\%$ saturated by D), and targets, (4) Li is fully stuck on the wall and targets, (5) Li re-deposition is ignored, and (6) the sputtering yield is not dependent on the local plasma conditions.

Zoom In Zoom Out Reset image size

The total emitted Li flux from the LM surface can be expressed as follows:

$$\begin{equation} \Gamma_\textrm{Li}^\textrm{Loss} = \Gamma_\textrm{Evp} + \Gamma_\textrm{ad} + \Gamma_\textrm{Sput}, \end{equation} \tag{ 4 }$$

$$\begin{equation} \Gamma_\textrm{Evp} = \frac{\alpha P_\textrm{Li}} {\sqrt{2\pi m_\textrm{Li}kT_\textrm{surf}}}, \end{equation} \tag{ 5 }$$

$$\begin{equation} \Gamma_\textrm{Sput} = f_\textrm{neut} \sum_j \Gamma_j Y_{j, \textrm{Li}}^{\,{ps}}, \end{equation} \tag{ 6 }$$

$$\begin{eqnarray} &&\Gamma_\textrm{Ad} = f_\textrm{neut}^{\,\,\textrm{ad}}\frac {Y_\textrm{ad}/Y_{ps}} {1 + A \textrm{exp} \left(\frac {E_\textrm{eff}} {k\textrm{T}_\textrm{surf}}\right)} \sum_j \Gamma_j Y_{j, \textrm{Li}}^{\,\textrm{ad}} . \end{eqnarray} \tag{ 7 }$$

Where $\textrm{log} P_\textrm{Li} = a - b/\textrm{T}_\textrm{surf}$, a = 5.055, b = 8023, m $_\textrm{Li}$ is the mass of Li atom, T $_\textrm{surf}$ is the surface temperature (output of MHD code), α is the sticking coefficient of 1, f $_\textrm{neut}$ is the fraction emitted as neutrals and assumed to be 0.35 [46], $\Gamma_j$ is the incident main ion flux (D) on an LM, Y $^{ps}_{\textrm{D},\textrm{Li}}$ = 10⁻³, $f_\textrm{neut}^{\,\,\textrm{ad}}$ = 1, the fitting coefficient (A = 10⁻⁷) has been obtained for D⁺ impinging on liquid Li in PISCES-RF linear plasma device [47], E $_\textrm{eff}$ = 0.9 eV, $Y_{\textrm{D,Li}}^\textrm{ad}$ = 10⁻³, and $Y_{\,\,\textrm{ad}}/Y_{ps}$ is the ad-atom to physical sputtering yield. E $_\textrm{eff} = E_\textrm{ad} - E_\textrm{D}$, where E $_\textrm{ad}$ is the ad-atom surface binding energy while E $_\textrm{D}$ is the surface diffusion activation energy [47]. E $_\textrm{D}$ was measured to be 0.75 eV [48]. Significant uncertainties exist in ad-atom binding energy, ad-atom yield, and constant A, as discussed in section 4.2. Compared to [49] we account for uncertainties in the ad-atom yield resulting from physical sputtering. Although Li redeposition may alter the sputtering flux, it might not have a significant impact on the evaporation flux. As a result, we anticipate that Li redeposition may not significantly alter the operating window identified in this work. In the future, we will investigate the impact of Li redeposition.

Self-consistent calculations of the emitted Li fluxes from an LM surface require liquid MHD and heat transfer calculations to determine the T $_\textrm{surf}$ given the incident plasma flux. Higher T $_\textrm{surf}$ can result in higher emitted Li fluxes, which in turn reduce T $_\textrm{surf}$ due to a reduction in the plasma heat flux via Li dissipation processes. Therefore, the SOLPS is coupled self-consistently with an LM MHD/heat transfer code [31] to understand these feedback processes and to find operating conditions consistent with the reactor design scenarios.

The LM MHD/heat transfer code computes the thickness of the Li layer, the flow velocity, and the Li T $_\textrm{surf}$, using the surface energy fluxes (figure 5) computed by SOLPS as input data. A more detailed description of the MHD code is provided in the appendix A.1. A schematic view of the coupling process is shown in figure 8. During the coupling procedure, the inlet flow velocity is maintained constant, while the heat flux profile is computed using SOLPS based on the emitted Li fluxes. The procedure is iterated until convergence is achieved.

Zoom In Zoom Out Reset image size

Significant uncertainties exist in Li loss processes, such as physical sputtering yield (Y $^{ps}$) and ad-atom to sputtering yield (Y $_\textrm{ad}$/Y $_{ps}$). The aim of this section is to understand the sensitivity of total Li loss to these parameters and establish a physics model for an LM Li-divertor. Although higher Y $^{ps}$ significantly increases the physical and ad-atom processes, evaporation loss mainly dominates the total loss processes at higher surface temperature regions ($\gt$525 ^∘C). The dependency of Y $^{ps}$ on the incident energy has been calculated and reported in [50]. Y $^{ps}$ due to D⁺ incident on Li strongly depends on the incident energy, rapidly reducing with the incident energy. SOLPS results have been used in a single-pass coupling to the hPIC code, finding Y $^{ps}$ in the range of 0.001 [46], consistent with the reported data in [50]. Consequently, Y $^{ps}$ = 0.001 is used in this simulation.

Ad-atoms are mobile atoms excited from their bound state on the LM surface by the incident plasma flux but lack sufficient energy to be sputtered from the surface [47]. As a result, ad-atoms require less energy to be emitted from an LM surface. The uncertainties in ad-atom losses arise from factors such as ad-atom yield, surface binding energy, and constant A [3, 46, 47]. The ad-atoms to sputtering yield (Y $_\textrm{ad}/Y_{ps}$) covers a broad range; previous experiments and modeling suggest a variation from approximately ∼1–70 [47]. Figure 7 illustrates the dependency of uncertainties in ad-atom losses. When Y $_\textrm{ad}/Y_{ps}$ (referred to Y/Y in figure 7) is 1, the total emitted flux is primarily dominated by the physical sputtering and evaporation flux. However, at higher yields, ad-atoms contribute significantly to Li losses in the temperature range of 200 ^∘C–525 ^∘C, while evaporation becomes the predominant driver of Li fluxes at higher temperatures (i.e, $\gt$550 ^∘C).

Significant uncertainties also exist in ad-atom surface binding energy; PISCES-B data suggested E $_\textrm{eff}$ = 0.70 eV [47], leading to a calculated E $_\textrm{ad}$ to 1.45 eV (E $_\textrm{eff}$ = E $_\textrm{ad}$—E $_\textrm{D}$; E $_\textrm{D}$ was measured to be 0.75 eV [48]), which is lower than the surface binding energy of Li (1.67 eV). As ad-atom binding energy increases, ad-atom losses decrease (see figure 7). Ad-atom losses are mostly saturated at the lowest binding energy, implying that lower binding energy primarily affects Li losses at lower surface temperature regions (i.e. $\lt$400 ^∘C).

Previous PISCES-B data shows good agreement with the experimental data at A = 6.9$\times 10^{-7}$ [3, 47]. Ad-atom losses are reduced with an increase in the magnitude of the fitting coefficient A (see figure 7). This constant also affects Li loss processes at the lower surface temperature range.

The Y $_\textrm{ad}/Y_{ps}$ is varied over a wide range (1–100) to understand its impact on the plasma parameters and the sensitivity is studied by SOLPS (refer to figures 16 and 17 in the appendix). Simulation results indicate that Y $_\textrm{ad}/Y_{ps}$ uncertainties show a minor impact on the main plasma (see also figure 17); ad-atom uncertainties only affect Li loss processes at the lower surface temperature regions (i.e. $\lt$525 ^∘C) as shown in figure 16, while evaporation flux mostly provides higher Li fluxes and contaminates the main plasma by fuel dilution. Given that the impact of Y $_\textrm{ad}/Y_{ps}$ uncertainties is very small at higher surface temperatures (refer to figure 16), Y $_\textrm{ad}/Y_{ps}$ = 1 is used for the subsequent section. In this work, A = 10⁻⁷, E $_\textrm{eff}$ = 0.9 eV, and Y $_\textrm{ad}$/Y$_{ps}$ = 1 are used for detailed Li sourcing in section 4.3, and all results and discussion in this work are valid based on these assumptions.

Figure 9 presents the calculated LM T $_\textrm{surf}$ and corresponding total emitted Li fluxes from an LM surface using equation (4). While maintaining a constant LM thickness of 5 mm, the inlet flow velocity was adjusted to achieve different peak T $_\textrm{surf}$ values in order to identify the upper limit of T $_\textrm{surf}$. The emitted fluxes from an LM surface increase with increasing surface temperature. Components of the emitted Li fluxes are shown in the appendix (refer to figure 18). For increasing surface temperature, the evaporation flux is strongly enhanced due to its strong dependency on the surface temperature (see figure 18).

Zoom In Zoom Out Reset image size

Figure 10 shows the dependency of the radiative power losses due to Li dissipation processes as a function of the peak LM surface temperature. The radiation power loss due to Li dissipation processes shows a weak dependency at the lower surface temperature due to lower emitted Li fluxes. However, the radiation loss is rapidly increased after 525 ^∘C (2D profile is shown in the appendix, figure 19), which indicates the contribution of evaporation processes, as evaporation mostly dominates Li loss processes after 525 ^∘C (see figure 18). Therefore, Li only reduces the divertor heat flux at the higher T $_\textrm{surf}$ range ($\lt$525 ^∘C) via radiation. Increasing T $_\textrm{surf}$ to 600 ^∘C, results in a radiation power loss of 15 MW.

Zoom In Zoom Out Reset image size

In this paper, for the sake of simplicity, the peak surface temperature (T $_\textrm{surf}^\textrm{max}$) represents the 1st iteration temperature (i.e. figure 3), while the simulation results in section 4.3 correspond to the converged condition results. Figure 11(a) illustrates the relationship between the final and initial condition T $_\textrm{surf}^\textrm{max}$ to understand the extent of temperature change from the initial condition due to Li dissipation processes, thereby illustrating the impact of the coupling feedback process. On the other hand, figure 3 depicts the relationship between T $_\textrm{surf}^\textrm{max}$ at a given Li flow and surface heat flux (single pass MHD code output). The main motivation of the coupling is to understand the dynamics and overall impact of the feedback process. During this process, the emitted Li fluxes reduce plasma heat flux through Li dissipation processes, particularly at low flow rates where low Li flow results in higher T $_\textrm{surf}$ and increased evaporation. Consequently, the final iterations yield a lower peak surface temperature compared to the initial conditions (refer to figure 11). This occurs because T $_\textrm{surf}$ depends on both the flow velocity and incident heat flux, demonstrating the influence of the coupling process.

Zoom In Zoom Out Reset image size

The higher emitted Li fluxes reduce plasma heat flux via Li dissipation processes at the higher surface temperature region. As a result, the converged case provides a slightly lower peak surface temperature compared to the initial (1st) iteration (see figure 11). We have identified the acceptable Li flow ($\gt$2 m s⁻¹), surface temperature ($\lt$525 ^∘C), and Li flux ($\lt \phi_\textrm{Li}$ ∼1×10²³ atom s⁻¹, corresponds to 1.118 gm s⁻¹) in terms of upstream Li concentration.

Li concentration at the OMP increases with rising surface temperature (refer to figure 12), and the Li concentration is found to be in the acceptable range when the surface temperature is kept below 525 ^∘C. This indicates an upper limit of Li sourcing range, and beyond that sourcing range, Li contaminates the upstream main ion density via fuel dilution. Within the assumptions of this model, an acceptable operating window for an LM Li divertor has been identified.

Zoom In Zoom Out Reset image size

Higher Li sourcing significantly reduces main ion density due to fuel dilution, while the electron density remains almost similar (the ions are fully stripped), as shown in figure 13. Fuel dilution (1 - $C_\textrm{Li}Z$) results in a 12$\%$ reduction in fusion power by C$_\textrm{Li} \sim$ 2$\%$. $Z_\textrm{eff}$ is increased slightly with increasing Li sourcing.

Zoom In Zoom Out Reset image size

Figure 12 shows the radial profile of Li concentration at the OMP. Lower surface temperature is correlated with lower Li concentration, while higher surface temperature significantly increases Li concentration at the OMP. Li neutrals are confined near the divertor surfaces while Li³⁺ is distributed over the whole plasma volume at high sourcing levels (see figure 20), reaching up to approximately ∼5% in the core.

Higher Li sourcing leads to plasma detachment by enhancing both energy and momentum losses (indicated by the pressure drop in figure 21). The radial profile of electron temperature at the outer divertor as a function of T $_\textrm{surf}$ is shown in figure 22. Combined puffing of D₂ ($4\times 10^{-22}$ atoms s⁻¹) and Ne ($1\times 10^{-20}$ atoms s⁻¹) is used. The impact of these gases has already reduced T $_\textit{e}$ to lower than 1 eV at the separatrix at the divertor, while T_e increases in the far SOL. The poloidal profile of T_e along the 3rd flux tube outside of the separatrix is shown in figure 23 to understand electron cooling due to Li in the upstream region. The description of 3rd flux tube can be found in [27]. It is shown that Li electron cooling effects can reduce T_e in the divertor region while showing a minor impact at the OMP (refer to figure 23).

Total plasma pressure ($n_eT_e+n_iT_i+0.5mu^2$) is strongly reduced at higher T $_\textrm{surf}$ on the outer divertor. On the other hand, a slight increase in plasma pressure is shown at the inner divertor when T $_\textrm{surf}$ is increased up to 550 $^{\circ} \textrm{C}$, as electron density is increased due to Li ionization processes. A rollover is shown in electron density, as higher Li sourcing reduces the electron temperature, increasing recombination processes. From the above discussion, Li LM T $_\textrm{surf}$ should be limited to 525 $^{\circ} \textrm{C}$ to avoid contamination of the main plasma by Li dissipation processes.

The mathematical model employed in this study for the analysis of open-surface LM MHD flows in a magnetic field along an inclined solid substrate to the horizon, with or without lateral walls including heat transfer is based on the boundary layer (thin-shear-layer) approximation as proposed for this type of LM MHD flows in [52]. In this work, the full 3D fluid flow is reduced to a 2D form using the standard boundary layer approximation [53]. In the case of a segmented design (S1 in equation (9)), the original 3D MHD flow equations are simplified to a quasi-two-dimensional (Q2D) form by integrating them in the direction of the applied (toroidal) magnetic field (z-direction) between the two lateral walls. This simplification is made under the assumption of a Hartmann-type velocity profile in the z-direction. The resulting Q2D flow equations, also known as the SM82 model, were first introduced in [54] for a nonconducting duct (quoted from the [31]). In [55], an attempt was undertaken to expand this model to an electrically conducting duct in the name of the 'averaged flow model'. This averaged flow model is employed in the present study for open-surface MHD flows.

In the present model, the lateral walls are treated as either thin electrically conducting or insulating walls, the choice depending on the wall conductance ratio c_w = ($\sigma_w t_w)/\sigma b$, where index 'w' denotes the wall, t_w is the wall thickness, b is the toroidal half-width of each segment, and σ and σ_w are the electrical conductivity of the liquid and that of the wall. The design goal is to make c_w as small as possible by using appropriate materials with low electrical conductivity. The MHD flow problem for the segmented divertor is characterized by the dimensionless Hartmann number Ha_tor , which is built using the toroidal magnetic field B_tor , dimension b and the physical properties (electrical conductivity σ, kinematic viscosity ν, and density ρ, refer to table 2) of the liquid:

Table 2. Li properties used in the MHD/heat transfer code.

µ (m2 s−1)	ρ (kg m−3)	cp (J (kg·k)−1)	σ (Ω−1 m−1)	m$_{Li}$ (kg m−3)	Lv (J mol−1)	k (W (m·K)−1)
$8.5\times 10^{-7}$	485	4200	$1.1\times 10^{-6}$	$1.17\times 10^{-26}$	145 920	49.6

$$\begin{equation} Ha_{tor} = B_{tor} b \sqrt{\frac{\sigma}{\nu\rho}}. \end{equation} \tag{ 8 }$$

Two other basic parameters, the hydrodynamic Reynolds and Froude numbers, are built as $R_e = (U_{in} h_{in})/\nu$ and $F_r = (U_{in}^2)/(gh_{in})$, where U_in is the velocity at the inlet, h_in is the flow thickness at the inlet, and g is the acceleration due to gravity. The flow aspect ratio is denoted by $\beta: \beta = h_{in}/b$. The electrical conductivity of the substrate is anticipated to exert minimal influence on the MHD flow, given that the induced electric currents form a closed circuit within the liquid through the side boundary layer at the substrate (the seeback effect is neglected in the present work). Importantly, the electrical resistance of this boundary layer is much smaller than that of the substrate. Consequently, there has been no effort to include the impact of the conducting substrate in the model. Following the integration, the problem for the sectioned divertor is reduced to 2D equations in the ($x, y$) plane. These new equations incorporate the influence of the flow resisting the electromagnetic Lorentz force, represented by a dedicated term F_em on the right-hand side of equation (9).

In the toroidally symmetric divertor design (S2 in equation (9)), the flow is 2D. The flow-opposing electromagnetic force arises due to interaction between the wall-normal magnetic field $B_{w-n}$ with the induced toroidal electric current $j_{tor} = \sigma UB_{w-n}$. The associated Hartmann number is based on the wall-normal magnetic field and the inlet flow thickness: $Ha_{(w-n)}$ = $B_{w-n}h_{in} \sqrt{\sigma/\nu\rho}$. The mathematical formulation for the MHD flow is similar between the two divertor designs except for the electromagnetic force F_em . All equations are written in the dimensionless form by using appropriate scales:

$$\begin{align} \frac{\partial U}{\partial t} + U\frac{\partial U}{\partial x}+ V \frac{\partial U}{\partial y} & = \frac{1}{F_r} \left[sin\left(\alpha\right)- \frac{\partial h}{\partial x}\textrm{cos}\left(\alpha\right)\right]\nonumber\\ &\quad + \nu \frac{\partial^2 U}{\partial x^2} + F_{em}, \end{align} \tag{ 9 }$$

$$\begin{align} F_{em} = \left\{ \begin{array}{l} - \beta^2 \left[\frac{Ha_{tor}}{\textrm{Re}} + \frac{c_w}{1+c_w} \frac{Ha^2_{tor}}{\textrm{Re}}\right]U \quad\text{for} \quad S1\\ - \frac{Ha_{w-n}^2}{\textrm{Re}}U \quad \text{for} \quad S2, \end{array} \right. .\end{align}$$

The symbol S1 represents the segmented divertor, while S2 represents the symmetric divertor.

$$\begin{equation} \frac{\partial U}{\partial x}+ \frac{\partial V}{\partial y} = 0, \end{equation} \tag{ 10 }$$

$$\begin{equation} \frac{\partial h}{\partial t}+ U_h \frac{\partial h}{\partial x} = V_h - V_1 \left(J_\textrm{Li}\right). \end{equation} \tag{ 11 }$$

Where h is the thickness of the LM. In this model, equation (9) is employed to compute the axial velocity component U. The continuity equation is utilized to compute the velocity component, V, normal to the substrate. Equation (13), known as the kinematic free surface condition, is applied to track the free surface using the computed velocities U_h and V_h at the free surface. The inclination angle of the substrate with respect to the horizon is denoted as α. Additional term V₁ ($J_\textrm{Li}$) in equation (11) that has a unit of velocity, stands for the evaporation of Li from the free surface. It is linearly proportional to the Li evaporation rate $J_\textrm{Li}$, which in turn depends on the Li temperature at the interface. In the dimensional form, this additional velocity term can be written as:

$$\begin{equation} V_1 = \frac{J_\textrm{Li} m_\textrm{Li}}{\rho}. \end{equation} \tag{ 12 }$$

Where ρ represents the density, as listed in table 2. It should be noted that the thin shear layer approximation, employed in the derivation of the equations, limits the applicability of the model to flows with moderate changes in flow thickness, i.e. when the condition $V\lt \lt U$ applies. Significant variations of the flow such that the velocity component V is comparable with the component U are outside the applicability range of the model. In particular, flow regimes, such as a 'hydraulic jump' cannot be predicted with the model. In the computations, the flow cases that demonstrate significant changes in the flow thickness with the axial distance x usually do not converge even when fine meshes and small time increments are used. The flow cases with no convergence are likely those experiencing the hydraulic jump.

The flow equations (9)–(11) need to be solved simultaneously with the energy equation (13), which is also written in the dimensionless form using the thin-shear-layer approximation. This equation is used to obtain the temperature distribution $T(x,y)$:

$$\begin{equation} \frac{\partial T}{\partial t}+ U \frac{\partial T}{\partial x} + V \frac{\partial T}{\partial y} = \frac{1}{Pe} \frac{\partial^2 T}{\ \partial y^2} + \frac{1}{Pe} q^{^{^{\prime\prime\prime}}}. \end{equation} \tag{ 13 }$$

In equation (13), Pe for Peclet number (Re × Pr), Re for Reynolds number, and Pr for Prandtl number, and $q{^{^{\prime\prime}}{}^{^{\prime}}}$ is the volumetric heating. The Prandtl number is Pr = $\nu\rho C_p/\kappa$, where C_p is the specific heat and κ is the thermal conductivity of the liquid. Li properties used in the code are listed in the table 2. The boundary conditions at the interface between the substrate and Li at y = 0 are the no-slip condition and convective He cooling. In this study, the thermal boundary condition is reduced to zero heat flux assuming high Li velocities. This implies that T diffusion is neglected, as the convective transport in the flow direction vastly outweighs diffusion. Consequently, only a negligibly small fraction of heat is removed with He:

$$\begin{equation} y = 0: U = 0, V = 0, \partial T/ \partial y = 0. \end{equation} \tag{ 14 }$$

At the free surface y = h, the tangential stress is zero and the heat flux includes the applied heat flux from plasma $q_0^{^{^{\prime\prime}}}$ minus the evaporation heat flux $q_1^{^{^{\prime\prime}}}$ $(J_\textrm{Li})$:

$$\begin{equation} y = h: \partial U/ \partial y = 0, \partial T/ \partial y = q_0^{^{^{\prime\prime}}} - q_1^{^{^{\prime\prime}}} \left(J_\textrm{Li}\right). \end{equation} \tag{ 15 }$$

Based on the mathematical model described by equations (9)–(15), the numerical code, detailed in [4], is employed to compute the developing LM MHD flows and the temperature field in the flowing liquid metal. The code achieves a steady solution through a time-advancement process, initializing with a slug-type flow and constant flow thickness. As quoted in [4], all equations are implicitly approximated using the finite-volume approach on a nonuniform mesh that clusters grid points near the solid substrate and free surface. The grid clustering is conducted using the stretching transformation, specifically designed for boundary-layer-type problems. In the case of a uniform mesh, the finite-volume scheme yields a second-order approximation. The momentum equation follows a marching type approach concerning the axial coordinate. In each time-step, the Blottner-type technique, developed specifically for marching problems [56], is employed to solve the equations. Typically, a single computation takes a few minutes on a personal computer using a mesh of 2001 by 201 points before reaching a steady state.

The EIRENE model has been improved compared to Lore et al [46] to incorporate Li sourcing in SOLPS via sputtering processes. The Li sputtering yield due to the incident main ion energy flux and angle is included in the TRIM data. As a result, EIRENE can calculate the physical sputtering yield due to the incident main ion (D⁺). However, Li loss processes also rely on the LM surface temperature, as discussed in section 4. Consequently, we have calculated an effective sputtering yield, which encompasses the total Li processes. For this purpose, we use the case with T$_\textrm{surf}^\textrm{max}$ = 450 ^∘C as a test case. The effective sputter yield is calculated as follows:

$$\begin{align} Y_\textrm{eff}^{\,j} = \frac {\sum_{i = 1}^n {\Gamma_{\textrm{Li},i}^\textrm{Loss} A_{i}}} {\sum_{i = 1}^n {\Gamma_{\textrm{D}^+}^i A_i}}. \end{align} \tag{ 16 }$$

Where j represents the number of segments (EIRENE surfaces corresponding to the B2.5 one) that split the target in the B2.5 mesh (here, j = 7 for both the outer and inner divertor), i is the index scanning the radial direction in the B2.5 mesh at the target, A is the area of the cells, n is the number of radial cells, $\Gamma_{\textrm{Li},i}^\textrm{Loss}$ is the emitted Li fluxes from an LM surface (equation (4)), and $\Gamma_{\textrm{D}^+}$ is the main ion flux.

The target plates are divided into multiple segments to achieve a better resolution depicting the sputtered Li fluxes. The area of each segment of the outer divertor surface and gas puffing surface is illustrated in figure 14(a). A relatively smaller area (segment 4) is allocated for the sputter model due to the rapid change of the main ion flux in these regions. The effective sputtering yield is then applied to the corresponding surface to source Li-neutral atoms from the LM surfaces. The energy of the emitted neutral is assumed to be at room temperature (i.e. 0.0259 eV). We have also investigated the impact of the emitted neutral energy and found a negligible impact, as the ionization potential of Li neutral is small, resulting in Li-neutral atoms being mostly ionized immediately by the plasma. However, the energy of Li-neutral atoms may influence ion energy via charge-exchange processes, which are not included in the present work. We will address this issue in the future.

Zoom In Zoom Out Reset image size

The effective sputtering yield is adjusted to yield similar emitted Li fluxes between this model and the gas puffing model (i.e. section 4). These two sourcing models are compared to validate the improved EIRENE model for future simulations. The total emitted Li for the gas puffing model is 9.64 × 10²¹ atom s⁻¹ while the sputter model provides a total emitted flux of 9.68 × 10²¹ atom s⁻¹, resulting in a relative error of 0.417$\%$.

The emitted Li fluxes from each surface via sputtering are depicted in figure 14(b). The Li-sputtering model reflects SOLPS output, while the magnitude of gas puffing, defined as an input, is shown for the gas puffing model. Both models exhibit similar levels of agreement, but the sputtering model is deemed more accurate and realistic. This is attributed to its use of the incident particle flux to sputter Li from the target plate, accounting for changes in the main ion flux during iterations. Therefore, we have tested the sensitivity of iteration on the effective sputtering yield (figure 15). Each iteration runs for 0.5 ms, and the output of each iteration is used to calculate $Y_\textrm{eff}$. The magnitude of $Y_\textrm{eff}$ varies across segments, as the product of $Y_\textrm{eff}$ and main ion flux is scaled to $\Gamma_\textrm{Li}^\textrm{Loss}$. Simulation results demonstrate that $Y_\textrm{eff}$ remains consistent across iterations, indicating the robustness and validity of the model. Higher Li sourcing has the potential to alter the main ion flux as a result of Li dissipation processes, introducing the possibility of different scenarios. While this aspect is beyond the scope of the present work, it will be explored in detail in a future paper.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The LM MHD/heat transfer code computes the thickness of the Li layer, the velocity, and the Li surface temperature, using the surface heat fluxes (figure 5) computed by SOLPS as input data. The calculated LM surface temperature and total emitted Li fluxes from an LM surface based on the corresponding SOLPS heat flux are shown in figure 9. The evaporation flux is strongly increased with increasing the surface temperature, as evaporation rates strongly depend on the surface temperature. The LM thickness was kept constant at 5 mm while the inlet flow velocity was adjusted to achieve different peak T$_\textrm{surf}$ to identify the tolerable T$_\textrm{surf}$. Table 3 shows the MHD simulation parameters, higher flow speed provides lower T$_\textrm{surf}$ due to the cooling. Li flow at the inner divertor is shown to be comparatively slower due to lower energy flux, while the outer divertor Li flow is shown to be higher to keep the same LM surface temperature.

Table 3. MHD code simulation parameters showing the relationship between inlet Li flow velocity and T$_\textrm{surf}^\textrm{max}$ based on the incident heat flux (figure 5).

T$_\textrm{surf}^\textrm{max} (^{\circ}\textrm C)$	U$_\textrm{O-div}$ (m s−1)	U$_\textrm{I-div}$ (m s−1)
450	7.5	1.50
500	3.0	0.53
525	2.0	0.40
545	1.5	0.37
565	1.0	0.35
600	0.5	0.275

Table 4. Li ions and neutrals are returned to the wall and targets depending on the sourcing rate: at lower sourcing rates, mostly Li ions are redeposited on the target, while at higher sourcing rates, Li neutrals are returned to the walls, targets, and pump. Here, the returned percentage is indicated by the ratio ($\phi_\textrm{Li}^\textrm{returned}/\phi_\textrm{Li}^\textrm{source}$).

		$\phi_\textrm{Li-ions}^\textrm{returned}$ ($\%$)
T$_\textrm{surf}^\textrm{max} (^{\circ} \textrm{C})$	$\phi_\textrm{Li}^\textrm{source}$ ($\times10^{23}$s−1)	φtargets	φwalls	$\phi_\textrm{Li-neutrals}^\textrm{returned}$ ($\%$)
450	0.271	93	1	5.9
500	0.094	83	1.5	14.9
525	2.32	73	1.3	25.28
550	3.35	72	1	25.7
565	4.71	71	1	27.5
600	13.66	51.1	0.67	47.3

The iteration process continued until the solution converged, meaning that the newly computed heat flux based on the emitted Li fluxes is applied as an input to the MHD code, and vice-versa (see figure 25). The emitted Li fluxes reduce plasma heat flux via Li dissipation processes at the higher surface temperature region. Consequently, the 2nd iteration and subsequent iterations provide a slightly lower peak surface temperature compared to the 1st iteration (see figure 11), as T$_\textrm{surf}$ depends on the flow velocity and incident heat flux. Furthermore, the surface temperature (i.e. emitted Li fluxes) becomes saturated in the higher iterations (see figure 26), indicating the convergence of the coupling process. More specifically, the convergence of the coupling processes is confirmed by checking both code outputs for each iteration, and the relative changes are found to be very small as shown in figure 25.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size