Dynamic modelling of local fuel inventory and desorption in the whole tokamak vacuum vessel for auto-consistent plasma-wall interaction simulations

Contents lists available at ScienceDirect Nuclear Materials and Energy journal homepage: www.elsevier.com/locate/nme Dynamic modelling of local fuel inventory and desorption in the whole tokamak vacuum vessel for auto-consistent plasma-wall interaction simulations J. Denis ⁎,a , J. Bucalossi a , G. Ciraolo a , E.A. Hodille 1,b , B. Pégourié a , H. Bufferand a , C. Grisolia a , T. Loarer a , Y. Marandet b , E. Serre c , JET Contributors 2 a CEA, IRFM, F-13108 Saint-Paul-lez-Durance, France b CNRS, Aix-Marseille Univ., PIIM, Marseille 13397 Cedex 20, France c Aix-Marseille Univ., CNRS, Centrale Marseille, Marseille, M2P2, France ARTICLE INFO Keywords: Plasma-wall interaction simulation Recycling Dynamics of hydrogen isotopes thermal desorption Dynamic retention Edge plasma physics ABSTRACT An extension of the SolEdge2D-EIRENE code package, named D-WEE, has been developed to add the dynamics of thermal desorption of hydrogen isotopes from the surface of plasma facing materials. To achieve this purpose, DWEE models hydrogen isotopes implantation, transport and retention in those materials. Before launching autoconsistent simulation (with feedback of D-WEE on SolEdge2D-EIRENE), D-WEE has to be initialised to ensure a realistic wall behaviour in terms of dynamics (pumping or fuelling areas) and fuel content. A methodology based on modelling is introduced to perform such initialisation. A synthetic plasma pulse is built from consecutive SolEdge2D-EIRENE simulations. This synthetic pulse is used as a plasma background for the D-WEE module. A sequence of plasma pulses is simulated with D-WEE to model a tokamak operation. This simulation enables to extract at a desired time during a pulse the local fuel inventory and the local desorption flux density which could be used as initial condition for coupled plasma-wall simulations. To assess the relevance of the dynamic retention behaviour obtained in the simulation, a confrontation to post-pulse experimental pressure measurement is performed. Such confrontation reveals a qualitative agreement between the temporal pressure drop obtained in the simulation and the one observed experimentally. The simulated dynamic retention during the consecutive pulses is also studied. 1. Introduction In fusion reactors, once the plasma is established, particle recycling on Plasma Facing Components (PFCs) represents the main particle source for the plasma. In steady-state conditions, the ratio between the gas throughput and the total recycling flux is estimated around 10% in present day tokamaks and is expected around 1% in ITER [1]. Hence, a complete understanding of the recycling phenomenon (reflection and molecular desorption of fuel) is important to ensure a reliable plasma density control, which has proven to be a critical issue for long-term operation of fusion reactors [2]. The recycling strongly depends on plasma facing materials. This phenomenon was highlighted at JET with the change from a carbon wall to the beryllium (Be) tungsten (W) ITERlike Wall (ILW). The change in the recycling process between both wall configurations has been pointed out as a possible explanation for the degraded confinement in H-mode observed in JET-ILW [3]. However, its experimental study remains challenging. Modelling can help the understanding of such process. SolEdge2D-EIRENE [4], with its ability to simulate plasma-wall interaction in complex wall geometry, is well suited for such modelling effort. Reflection is already handled in EIRENE through tabulated TRIM results [4]. An extension has been developed to model the dynamic thermal desorption of implanted fuel particles (Desorption from Wall ElemEnts, D-WEE) to complete the description of the recycling process. To achieve this purpose, D-WEE simulates hydrogen implantation, transport and retention in materials using Reaction-Diffusion equations. However, before launching coupled simulations, one has to define a realistic initial wall state (local fuel inventory and local desorption flux https://doi.org/10.1016/j.nme.2019.03.019 Received 14 August 2018; Received in revised form 19 March 2019; Accepted 19 March 2019 ⁎ Corresponding author. E-mail address: [email protected] (J. Denis). 1 Present working address: Department of Physics, University of Helsinki, P. O. Box 43, FI-00014, Finland. 2 See the author list of “Overview of the JET results in support to ITER” by X. Litaudon et al., Nuclear Fusion 57 (2017) 102001. Nuclear Materials and Energy 19 (2019) 550–557 Available online 04 June 2019 2352-1791/ © 2019 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/BY-NC-ND/4.0/). T density), ideally representative of the experimental one during a pulse. It is fundamental to determine this local wall state as it sets its behaviour w.r.t. the fuel (pumping or fuelling areas) and the available fuel reservoir before transient events. Such wall state cannot be evaluated experimentally. Indeed, post-mortem analysis are mostly performed weeks after tile removals, and therefore are representative of long-term fuel retention. The fuel uptake by the wall during pulses and its subsequent release between pulses, often referred to as dynamic retention, cannot be assessed by this type of analysis. In JET-ILW diverted discharges, a ratio dynamic retention to long-term retention of 10 have been measured [5]. This indicates the prevailing role of dynamic retention in this machine (which is mainly due to its short duty cycle). Only gas balance analysis gives a global idea of this dynamic retention [6]. Still the local wall state cannot be estimated from this kind of experimental studies. In the present contribution, a methodology to define this local wall state is presented. A synthetic plasma discharge is built from successive SolEdge2D-EIRENE simulations. This synthetic pulse is used as a plasma background for D-WEE. This code is used in standalone mode (independently of SolEdge2D-EIRENE) to investigate the wall dynamics during a sequence of plasma pulses. From that type of simulations, the wall state at a desired time during a pulse can be calculated and then be used as an initial condition for coupled SolEdge2D-EIRENE – D-WEE simulations (not shown in this contribution). Such methodology is applied to the JET tokamak in Section 3 after a brief introduction of the DWEE module in Section 2. Especially, the relevance of the dynamic retention behaviour obtained in the simulation is assessed in Section 3.3.2 through confrontation to post-pulse experimental pressure measurement and the simulated dynamic retention during the consecutive pulses is studied in Section 3.3.3. 2. D-WEE: a dynamic thermal desorption module for SolEdge2DEIRENE In this section, the D-WEE module developed for the SolEdge2DEIRENE code package is introduced. SolEdge2D-EIRENE is able to give the plasma parameters (density, particle flux density and temperature) along the complex Vacuum Vessel (VV) wall. In the following, the s symbol denotes the curvilinear coordinate along this wall (cf. Fig. 1 for the definition of this coordinate in the case of the JET tokamak), while the x symbol denotes the depth of the material (transverse to s ). 2.1. Desorption model: MHIMS The desorption model used in D-WEE is based on the code MHIMS (Migration of Hydrogen Isotopes in MaterialS) [7]. This code uses the Reaction-Diffusion system of equations (R-D equations) [8] to describe the transport and trapping of Hydrogen Isotopes (HI) in metals. In this model, two HI populations are considered: the mobile HI, which can diffuse in the metal lattice, and the HI which are trapped in the lattice defects (vacancies, grain boundaries, etc.). This system is made of a diffusion equation for the mobile particles which is coupled to equations describing the kinetic of reaction between the mobile particles and the available traps of kind i +HI Trap HI mit,i : = + + = + = + n s x t t x D T n x n tS S n s x t tT n T n n n ( , , ) ( ) ( , , ) ( ) ( ) ( ) m m i N t i ext iext at t i dt i t i t i m i t i 1 , , , , , , trap (1) where nm and nt i, are the mobile particle number density and the trapped particle number density in the trap of kind i respectively [ m3 ], T is the wall temperature [ K ], =D T D E k T( ) exp( /( )) D B0 is the diffusion coefficient of the HI in the considered material [ m .s 2 1 ] with activation energy E[eV] D , kB is the Boltzmann constant [ eV.K 1 ], =E k Texp( /( )) dt i dt i B, 0 , is the detrapping frequency of the trap of kind i [ s1 ] with activation energy E[eV] dt i, , =T D T n( ) ( )/( ) t i solute,2 is the trapping rate coefficient [ m .s 3 1 ], is the distance between two HI solute site [ m ], nsolute the HI solute site number density [ m3 ], Ntrap the number of traps in the material and n s x( , ) i is the trap number density [ m3 ]. The terms + Sext i and Sext at are the volume sources due to implantation of ions and atoms respectively: =S s x t s t f E s t s t( , , ) ( , ) ( ( , ), ( , )) ext j imp j ximp j imp j (2) where s t( , ) imp j is the particle implantation flux density and fx is a form function given by the SRIM code [9]. This form function depends on the mean particle impact energy Eimp j [eV] and the mean angle of incidence imp j [ ∘ ]. The coupling with SolEdge2D-EIRENE is made through these three terms ( imp j , Eimp j and imp j ). The system (1) is solved considering a Dirichlet boundary condition at the front surface of the material ( =n s 0 t( , , ) 0 m , the recombination and desorption steps are assumed to be immediate) and either a Dirichlet or a Neumann boundary condition at the rear surface (at =x L , respectively =n s L t( , , ) 0 m or = n x s L t( , , ) 0 m ) depending on the considered material. The feedback to SolEdge2D-EIRENE is performed through the term Rm , that one calls molecular recycling coefficient: = + + R s t s t s t s t ( , ) ( , ) ( , ) ( , ) m out imp i imp at (3) where out is the desorption flux density [ m .s 2 1 ], defined as the diffusive flux density at the material surface. This coefficient should not be confused with the total recycling coefficient as fuel reflection of both ions and atoms should also be considered in its calculation. Moreover, it should be noted that the flux density of molecules desorbing from the wall is half the value of out . 2.2. Thermal model: WEE-temp Diffusion and detrapping in Reaction-Diffusion model are thermally activated processes. Therefore, such model has to be coupled with a thermal model to estimate the temperature depth profile in each PFC. In D-WEE, such calculation is performed by the Wall ElemEnts temperature code (WEE-temp), which solves the 1D heat equation through the so-called Duhamel integral [10]: = +T s x t T s x s T s x t td( , , ) ( , , 0) ( , ) ( , , ) t net step 0 (4) where T s x 0( , , ) is the initial PFC temperature [K] , net is the net heat flux density [ W.m 2 ] given by SolEdge2D-EIRENE (taking into account the plasma and radiation contributions), and Tstep is the PFC temperature response to a unit heat flux density step [ K.(W.m ) 2 1 ]. The convolution integral in Eq. (4) is calculated using Laplace transform. The step response is calculated using the quadrupole method [11] which enables the modelling of actively-cooled and inertiallycooled PFCs with different material layers [12]. This technique enables the calculation of the temperature in the zone of interest (defined by the depth simulated in MHIMS, ≤ 1 mm) while taking into account the boundary condition at the PFC back surface (at a several mm depth). 3. Modelling of dynamic desorption and retention during JET discharges In this section, the methodology to define the local wall state (i.e. local fuel inventory and local desorption flux density) is presented and applied to the JET tokamak. The SolEdge2D-EIRENE code is used to simulate the different steady-state plasma phases of a given JET pulse. A synthetic plasma pulse is built from these simulations, which are used as consecutive fixed plasma backgrounds for the D-WEE module. A sequence of plasma pulses can then be simulated with D-WEE. For this simulation, the initial deuterium density in the wall was taken equal to J. Denis, et al. Nuclear Materials and Energy 19 (2019) 550–557 551 zero. At this stage, four successive pulses followed by 30 min resting time were simulated. 3.1. Experiment The present study focuses on the JET pulse number 89044 (#JPN89044) ( = =B I2.4 T, 2 MA t p , pure D plasma). The time evolution of the plasma pulse parameters is plotted in Fig. 2. The plasma is run in diverted configuration with inner and outer strike-points located on the top of vertical tile 3 and on the horizontal tile 5 Stack C respectively (cf. Fig. 1). Two steady-state plasma phases are identified, defined by constant magnetic equilibrium (X-point configuration and strike-points location), particle injection rate, input power and plasma density: (1) H-mode phase: from 48.3 s to 50.3 s. 15 MW of input power ( =P12 MW NBI ), gas puff injection rate = ×Q4.1 10 D.s inj 21 1 from the divertor base inner ring, NBI injection rate = ×Q1.5 10 D.s NBI 21 1 and upstream line-integrated electron density = ×n1.2 10 m e 20 2 . In the following, only the inter-ELM plasma conditions are considered. (2) L-mode phase: from 52 s to 59 s. 2.9 MW of input power, gas puff injection rate = ×Q4.1 10 D.s inj 21 1 from the tokamak midplane and = ×Q6.9 10 D.s inj 21 1 from the top of the machine and upstream line-integrated electron density = ×n9.2 10 m e 19 2 . One can wonder whether these two steady-state phases are representative of the plasma-wall interaction of the complete pulse. The answer to this question is complex and is probably negative. However, albeit both of them only represent 35% of the pulse duration (65% of the duration of the plasma current flat-top phase), almost 60% of the total divertor ion fluence measured by Langmuir probes is deposited during this two phases. If one only considered the flat-top phase, this percentage increases to 72%. Due to the impossibility to simulate the plasma ramp-up and ramp-down phases with SolEdge2D-EIRENE, the simulation of these two steady-state plasma phases can be considered as the best possible way to estimate the overall plasma-wall interaction of this pulse. 3.2. Simulation set-up 3.2.1. SolEdge2D-EIRENE plasma backgrounds The two identified plasma phases are simulated with the SolEdge2D-EIRENE code. Pure D plasma is considered in these simulations. The cross-field transport coefficients for particle and energy are automatically set up to fit experimental reference upstream profiles of density and temperatures (taken at the JET outer mid-plane). These reference profiles are extracted from coherently time averaged high resolution Thomson Scattering measurements during each phase of #JPN89044 (only inter-ELM profiles are considered for the H-mode phase since the modelling of the ELM transients is not conceivable due to its computational cost). They are also automatically shifted with respect to the magnetic separatrix in order to force the Scrape-Off Layer (SOL) input power ( P SOL ) to match the experimental one ( =P P P SOL in rad core, with P in the input power and P rad core, the radiative power in the plasma core, =P10 MW SOL for phase 1, 2.5 MW for phase 2). The experimental fuelling rates and locations (cf. Section 3.1) are respected and the NBI fuelling is considered during the H-mode phase. A Be first-wall (from 1 → 2 clockwise in Fig. 1) and a W divertor (from 2 → 1) are considered for particle reflection properties. A Be surface is also assumed on the top of HFGC tile and tile 1, in agreement with postmortem analysis that revealed the presence of a thick (up to 40 µm ) Be co-deposit layers in this zone of JET-ILW [13]. A molecular recycling coefficient ( Rm ) of unity is considered everywhere in the wall (under the assumption that plasma conditions are weakly dependent on Rm ), except in the two divertor throats (where the pumping ducts are located). A specified value is set there in order to force the pumping flux, Qpump , to ensure a steady particle balance ( + =Q Q Q inj NBI pump ) in the simulations. The SolEdge2D-EIRENE simulation results for the H-mode phase are shown in Fig. 3, where the distribution along the wall of the required Fig. 1. JET poloidal cross-section with the different PFCs considered in the following simulation. Points 1 and 2 indicate the boundary between the two regions of the vacuum vessel: from 1 → 2 clockwise the first-wall, from 2 → 1 the divertor. The curvilinear coordinate along the wall, s , is also defined. Fig. 2. Time evolution of pulse parameters for #JPN89044: (a) plasma current I p , (b) gas puff injection rate Q inj , (c) input and radiated powers, (d) line-integrated electron density n e , (e) inner strike point vertical position ZSIL and outer strike point horizontal position RSOL (cf. Fig. 1 for the JET coordinates). Two steady-state phases, labelled 1 and 2, are identified. They are characterized by constant pulse parameters. J. Denis, et al. Nuclear Materials and Energy 19 (2019) 550–557 552 D-WEE inputs is plotted. The implanted particle flux density presents two peak in the strike-points location ( ×2 10 m .s 22 2 1 ) with high impact energy ( ∼ 600 eV). A third ion flux density peak is observed in the inner side of the divertor where the Be co-deposit is located (high recycling region). Overall, the implanted particle flux density for both ions and atoms remains around 10 m .s 21 2 1 and 10 m .s 20 2 1 in the divertor and the first-wall region respectively, with impact energy of ∼ 100 eV. The mean angle of incidence are not plotted. It remains between 60 ∘ and 70 ∘ for the ions in all the vacuum vessel, while its value is around 40 ∘ for the atoms. Eventually, the heat flux density remains below 0.1 MW.m 2 in all the VV, with excursion in the strikepoints vicinity up to 5.8 MW.m 2 in the outer strike-point. The SolEdge2D-EIRENE results for phase 2 are not shown in this contribution. For the first-wall, results are quite similar to those of phase 1, except for the impact energy which decreases by a factor of 2. In the divertor region, the plasma is detached. The divertor implanted ion flux density decreases by a factor of ∼ 2 while the one for atoms increases by one order of magnitude, becoming prevailing with e.g. an outer strike-point value of ×4.5 10 m .s 22 2 1 . The impact energy remains below 10 eV and the heat flux density decreases up to 0.1 MW.m 2 . In conclusion, the total implanted particle flux ( Qimp tot, , including + D and D), calculated through integration of the implanted particle flux densities (assuming toroidal symmetry of the plasma-wall interaction), are ×1.6 10 D.s 23 1 and ×3.6 10 D.s 23 1 for phase 1 and phase 2 respectively. It represents an increase of the implantation by a factor of ∼ 2.3 between the two phases. However, if only ions are considered, the total implanted ion flux decreases by a factor of ∼ 0.74 between the two phases ( ×+ 1.2 10 D .s 23 1 vs. ×+ 8.9 10 D .s 22 1 , the difference with Qimp tot, representing the atomic contribution). These results show that a non-negligible part of the implanted D - the majority in the case of phase 2 - comes from D atoms. 3.2.2. D-WEE simulation set up JET-ILW tokamak is composed of inertially-cooled PFCs. In WEEtemp, such PFCs are simulated like actively-cooled PFCs [12] with a very low heat convection coefficient ( ×1 10 W.m .K 3 2 1 ). This enables the model to mimic the cooling of the PFCs by radiation and by heat conduction into the tokamak structure between the pulses, while keeping their inertial behaviour between plasma exposure. In the following simulation, 19 different PFCs have been considered, which are distributed in 22 different zones in the VV (cf. Fig. 1). Each PFC is defined by its material layers (massive W, W coating, CFC, Be etc.). The value of the material thermal properties (thermal conductivity, specific heat and material density) have been considered at the closest temperature to the base temperature of the different regions of the VV (i.e. at 200 ∘ C for the first-wall, at 100 ∘ C for the divertor). For each PFC, the main material layer thickness is set to an effective value, i.e. =e V S/ eff tile tile , where V tile and Stile are respectively the volume and the wetted area of the tile (taken from CAD drawings). The values of eeff range from 2 to 5 cm. This procedure enables the different zones where the heat flux density is uniform to experience a good inter-pulse temperature increase plasma after plasma. But such procedure is not working in regions where the heat flux density is peaked. In order to avoid inter-pulse over-heating, the two PFCs hosting the strike-points, i.e. tile 3 for inner strike-point and tile 5 Stack C for outer strike-point, have seen their effective thickness increased to 21 cm for the CFC layer and 23 cm for the W layer respectively. Eventually, the base temperature of the different PFCs is taken from thermocouples measurement at the beginning of the session for the divertor region, leading to temperatures ranging between 50 and 70 ∘ C. Such measurements are not available for the first-wall PFCs. They are considered at the minimum working temperature of JET-ILW first-wall (i.e. 200 ∘ C) at the start of the simulation. Concerning MHIMS, the main difficulty is to define the distribution Fig. 3. SolEdge2D-EIRENE simulation results for plasma phase 1 (H-mode, inter-ELM): distribution along the wall of the implanted particle flux density ( imp ) and the impact energy ( Eimp ) for both ions and atoms as well as the net heat flux density ( net ). These quantities are required as inputs for the D-WEE module. The definition of the curvilinear coordinate s can be found in Fig. 1. J. Denis, et al. Nuclear Materials and Energy 19 (2019) 550–557 553 of materials at the wall surface. Indeed, after two campaigns, JET-ILW exhibits a complex pattern of surface materials, with a mix of the original materials (massive W, W coating, massive Be, etc.) in the net erosion zones and co-deposit layers in remote areas [14]. Moreover, the lack of suitable retention models (eg. for H-Be interaction) and/or parameterisation of the R-D equations for the aforementioned materials hamper the modelling of the ILW. As an initial simplified approach, a full massive W wall configuration have been considered in the simulation. The trapping parameters (trap densities, depth profiles and detrapping energies) have been parameterised for polycrystalline W in a previous study [7]. Three kinds of defects have been identified: two intrinsic traps and an ion-induced trap. The latter one represents defects that are created by the simultaneous presence of impurities (C, O, etc.) with a large concentration of hydrogen isotopes in the W lattice during plasma exposure. This trap presents a high density in the implantation zone and a weaker density in the bulk of the material. The bulk part extent of its depth profile is linked to diffusion of impurities in the depth of the material. The higher the W working temperature is, the deeper the impurities diffuse and the deeper type 3 traps are created [15]. In the following, it was assumed that the divertor had experienced a working temperature of 1000 ∘ C during plasma operation while the first-wall was maintained at its base temperature of 200 ∘ C. This difference of working temperature leads to a deeper trap 3 density profile in the divertor than in the first-wall. The considered trap parameters are sum up in Table 1. They are considered to be representative of a W wall which has experienced hours of operation. Moreover, the detrapping pre-exponential factor 0 is taken to be equal to the Debye frequency, i.e. 10 s 13 1 . The solute sites for D in W are the tetrahedral sites which entails that =116 pm and = ×n6 solute W with W the W number density [7]. The diffusion coefficient from density functional theory (DFT) is taken in [16]: = ×D T k T( ) 1.9 10 exp( 0.2/( )) B 7 m .s 2 1 . This diffusion coefficient was calculated for hydrogen. To take into account the isotope effect, the diffusion pre-exponential factor is divided by 2 , the square root of the mass ratio between deuterium and hydrogen. Eventually, a Dirichlet boundary condition at the rear surface is considered (cf. Section 2.1) as the simulated depth of the material (0.1 mm) is high enough to ensure that no mobile deuterium reaches this surface in the time period of the simulation. In the following, an arbitrary differentiation of the traps is used, which is defined to highlight their dynamic retention capability in the following simulation. For the first-wall, traps 1 and 2 are gathered and labelled as low energy traps (with E1 eV dt i, ) while trap 3 is labelled as high energy trap. For the divertor, the low energy trap refers to trap 1 while the high energy traps (with E1 eV dt i, ) refer to traps 2 and 3. 3.3. Results In the following, the wall dynamics during the four pulses and in between pulses is studied. No temperature evolution is shown, as the majority of the PFCs remains around its base temperature even during plasma. Only the two strike-point tiles exhibit a large increase of temperature during the plasma H-mode phase ( + °450 C for the outer strike-point, + °280 C for the inner one). 3.3.1. Wall outgassing and inventory in between pulses The total D outgassing flux Qout tot, [ at.s 1 ] and total D wall inventory Nwall tot, [ at ] are obtained by integration of the desorption flux density and of the particle areal density respectively (assuming toroidal symmetry of the plasma-wall interaction). The time evolution of these quantities is shown in Fig. 4. The total outgassing flux undergoes a steep increase (5 orders of magnitude) during each pulse (due to plasma exposure), followed by a sharp decrease at the plasma shutdown. This decrease is then weaker up to the next pulse. The outgassing flux from different zones of the VV is also displayed. The first-wall PFCs are clearly dominating the outgassing in between pulses. The wall inventory time evolution is also plotted for the first-wall and for the divertor, with the differentiation for the low and high energy traps define in Section 3.2.2. Most of the total wall inventory is found in the high energy traps (cf. Fig. 4b), mainly in trap 3 due to the combination of its high detrapping energy with its high number density at the material surface. However, the high energy trap inventories do not share the same temporal evolution in between pulses: the divertor inventory tends to increase slightly (+2% in the 30 min after the 4 th pulse) while the first-wall one decreases (−10% in the same time period). This can be explained by the different base temperatures of the two VV regions. Indeed, the first-wall remains at 200 ∘ C leading to detrapping from trap 3 in between pulses (while trap 1 and 2 do not seem to retain deuterium during the pulses, cf. Fig. 4c). In the divertor, the base temperature is between 50 and 70 ∘ C. Therefore, Table 1 Trap parameters used in the simulation. The trap densities are given in at.% of W. For intrinsic traps, the trap densities are expected to be constant in all the depth of the material and equal to the values presented in this table. For the ioninduced traps, the densities presented here are the maximum densities up to the depths reported in the table. Divertor First-wall Trap 1 (intrinsic) =E0.85 eV dt,1 =n0.13 at.% 1 Trap 2 (intrinsic) =E1.00 eV dt,2 =n0.035 at.% 2 Trap 3 (ion-induced) =E1.5 eV dt,3 =n1 at.% 3 a up to 10.8 µm =n1 at.% a 3 up to 2.7 µm =n15 at.% 3 b up to 30 nm =n15 at.% 3 b up to 20 nm Fig. 4. Time evolution of simulated D outgassing flux Qout (a) and D wall inventory Nwall (b-c). Each peak of outgassing indicates a plasma pulse. The outgassing fluxes from different zones of the vacuum vessel are displayed. The first-wall and divertor traps inventory are arbitrarily differentiated according to their detrapping energies (cf. Section 3.2.2). J. Denis, et al. Nuclear Materials and Energy 19 (2019) 550–557 554 trapping is highly efficient in the high energy traps (i.e. traps 2 and 3), even between the pulses. The low energy trap also retains deuterium during plasma exposure, but releases its population in between pulses. These detrapped particles diffuse in the material, partly up to the surface (leading to outgassing), but also in the bulk material where they fill the high energy traps. This phenomenon has been exposed in [7] and will be highlighted in the following section. However, this effect is only transitory and will decrease plasma after plasma as the bulk high energy traps get more and more saturated in the depth of the material. 3.3.2. Short and long term outgassing after the last plasma pulse After a plasma discharge, outgassing of particles from the wall occurs and is the signature of dynamic retention of D in the wall during the pulse. In JET-ILW, such outgassing leads almost to a complete recovery of the injected particles during the pulse [5]. This flux can be evaluated by means of pressure measurements in the vacuum vessel after the discharge. In the following, a method is proposed to compare the time evolution of the simulated VV pressure with experimental pressure measurement. First, the general particle balance equation for the tokamak VV has to be introduced: + + = d dt N N N Q Q( ) [at.s ] ions atoms wall tot inj tot pump tot, , , 1 (5) where Nions is the total number of ions, Natoms the total number of atoms, Nwall tot, is the total number of particles in the wall (wall inventory), Qinj tot, the total external injection flux and Qpump tot, the total pumping flux. The time variation of the wall inventory is obtained by writing the global flux balance at the wall surface, =dN dt Q Q/ wall tot imp tot out tot, , , , where Qimp tot, is the total implanted particle flux and Qout tot, is the total outgassing flux. When >Q Q imp tot out tot, , , dN dt/ wall tot, is positive and is referred to as retention flux or rate ( Qret tot, ). The gas balance method uses Eq. (5) to calculate this retention flux [6]. In between pulses, plasma is shut down and no external injection is operating. The VV pressure P can be obtained using the ideal gas law =PV N k T0.5 atoms B , where V is the VV volume [ m3 ], kB the Boltzmann constant [ J.K 1 ] and T the gas temperature in the VV [K] . Deuterium in the VV are in the form of D 2 molecules which explains the 0.5 factor. Moreover, =Q PS k T2 /( ) pump tot VV eff B, , where SVV eff is the total effective D 2 pumping speed related to the VV pressure P . This pumping speed embeds all the contributions from the external active VV pumps (cryogenic, NBI box and turbo-molecular pumps). In this situation, Eq. (5) yields: = + V k T dP dt PS k T Q0.5 B VV eff B out tot, (6) Then the VV pressure can be obtained by considering the simulated Qout tot, and by integrating numerically Eq. (6). After the 4 th plasma, a 100 h resting time is simulated to study the short and long term outgassing behaviour of the JET wall in full-W configuration. The VV pressure is calculated using Eq. (6). SVV eff has been estimated through a calibrated gas injection (dry-run) performed at the beginning of the experimental session ( =S135 m .s VV eff 3 1 ). The JET VV volume is used ( =V185 m3 ) and the gas temperature in the VV is considered at the first-wall temperature ( = °T200 C ). In Fig. 5, the time evolution of the experimental post-pulse VV pressure and of the simulated one is displayed. The simulation clearly underestimates the VV pressure, by a factor of ∼ 7 right after the plasma shutdown. This factor decreases to a value ∼ 3 13 min later (end of pressure measurement). However, the simulated pressure follows a t0.74 trend during 100 h, close to the experimental trend ( t0.89 ). Such trends have also been observed experimentally in carbon machines [17], in JET with Be first-wall [18] as well as in other JET-ILW experiment [19]. The time variation of the wall trap inventories are also displayed in Fig. 5, with the traps differentiation defined in Section 3.2.2. A positive time variation indicates a filling of the traps while a negative one indicates that the traps are emptying, leading predominantly to outgassing. As explained in Section 3.3.1, after the plasma shutdown, the low temperature traps are emptying and a part of this detrapped particles are repopulating the high energy traps. This phenonemon lasts longer in the divertor region (up to 2000 s after the plasma shutdown) because the low energy trap retains more deuterium during plasma (due to the lower base temperature in this region). The short and long term outgassing are dominated by detrapping from trap 1 in the divertor, followed by detrapping from trap 3 in the first-wall. These traps are thus responsible for the dynamic retention of the wall. However, the quantity of desorbed D in between pulses is clearly underestimated in the simulation (1.3 × 10 22 at difference between the simulation and the experiment after 13 min). The inventory of the first-wall trap 3 at the Fig. 5. Time evolution of the simulated and experimental vacuum vessel pressures (a) and of the time variation of trap inventories ( dN dt/ wall , negative in (b) and positive in (c)). The first-wall and divertor trap inventories are arbitrarily differentiated according to their detrapping energies (cf. Section 3.2.2). J. Denis, et al. Nuclear Materials and Energy 19 (2019) 550–557 555 end of the 4 th plasma is high enough to entail such amount of outgassed particles between pulses (cf. Fig. 4b). This indicates that the dynamics of this trap is incorrect, probably due to a too high detrapping activation energy. 3.3.3. Wall outgassing and retention during pulses This section focuses on the wall dynamics during the pulses. In the upper part of Fig. 6, the time evolution of Qout tot, during the H-mode phase of the four consecutive plasmas is displayed. Qimp tot, during this phase is also plotted. As expected, Qout tot, tends to Qimp tot, , which is a sign of the wall surface saturation. Its time evolution is very similar plasma after plasma and shows from the second plasma an apparent quick saturation in less than a second. The lower part of Fig. 6 shows the time evolution of different retention flux: the total retention flux, =Q Q Q ret tot imp tot out tot, , , , and the retention flux related to the traps, defined as the positive time derivative of their respective inventory (similar to Section 3.3.2). If one omits the first plasma, which exhibits a transitory behaviour with a filling of high energy traps in both the divertor and the first-wall region, the following plasmas present a similar retention behaviour. At the beginning of these plasmas, the retention flux is clearly dominated by the first-wall high temperature trap and by the divertor low temperature trap. The latter one tends to see its contribution collapsing within a second while the retention flux in the divertor high temperature traps increases. Indeed, once the surface low temperature traps have reached saturation, D can diffuse in the bulk material and populate the high energy traps. This phenomenon declines plasma after plasma, as these traps get more and more saturated in the bulk, and will vanish after a series of plasmas. This saturation benefits to trap 1 which sees its dynamic retention slightly increasing plasma after plasma. However, even for the 4 th plasma, the dynamic retention in this trap still collapses and is mostly dominated by trapping in the first-wall trap 3. It probably results from a combination of a low implantation flux density in this region with the high density of trap 3 at the surface, which induces a slower trap filling. For the 4 th plasma, after two seconds of H-mode phase, the simulated total retention flux is ×1 10 D.s 21 1 and keeps on decreasing. The wall is therefore close to saturation. A direct confrontation to the experimental retention flux obtained from gas balance analysis is not conceivable as retention through Be-D co-deposition process is not simulated. Still, one can note that the order of magnitude of the obtained retention flux and its time evolution are consistent with what is experimentally observed [20]. The wall state during plasma can be extracted from such simulation, e.g. after two second of H-mode phase for the pulse number 4 (cf. Fig. 7). The divertor is completely saturated, with Rm ranging between 0.97 and 1 while the first-wall is partly saturated, with some points where Rm ranges between 0 and 0.8. These points match with locations where only atoms are striking (cf. Fig. 3) and therefore where the total particle flux density is smaller: the atom flux density is not sufficient to saturate trap 3 at the surface. 4. Conclusion A dynamic thermal desorption model has been developed for the SolEdge2D-EIRENE transport code package. A sequence of plasma discharges has been simulated in order to estimate an initial wall state for Fig. 6. Time evolution of total outgassing flux Qout tot, (a) and total retention flux Qret tot, (b) during the H-mode phase of each of the four pulses. The total implantation flux is also displayed in (a). The first-wall and divertor trap retention fluxes (b) are arbitrarily differentiated according to their detrapping energies (cf. Section 3.2.2). Fig. 7. Distribution along the JET wall of the molecular recycling coefficient Rm calculated by D-WEE after 2 s of H-mode. The positions of both pumping ducts, where Rm is not calculated by D-WEE but is directly forced in the SolEdge2D-EIRENE simulation, are also indicated. J. Denis, et al. Nuclear Materials and Energy 19 (2019) 550–557 556 coupled plasma-wall simulation in the case of a tungsten wall in the JET tokamak. Even though such wall is not representative of the one in JETILW, some conclusions can be drawn from this preliminary simulation: 1) The reaction-diffusion model is able to reproduce experimental trends (order of magnitude of the dynamic retention flux during pulses and time evolution of post-pulse outgassing), if the required trapping parameters and the overall wall condition (distribution of particle flux density, impact energy, etc.) are considered. Dynamic retention is set by traps with efficient detrapping around the PFCs base temperature between pulses. This traps cannot be detected in post-mortem analysis as they empty in the time period between the last plasma exposure and the analysis. Deuterium implantation at low-energy and subsequent in situ thermal desorption spectrometry [21] of post-mortem samples could be used to characterise these traps. The fitting of the experimental post-pulse pressure time evolution using the method described in this paper could also be envisaged. 2) In the simulation, the dynamic retention is sustained by a combination of low implanted particle flux density ( <10 m .s 20 2 1 ) region with high surface density traps, which lead to a progressive saturation of the wall surface. Such flux density are found in the firstwall region where no ions are impacting. JET-ILW first-wall is mainly made of Be (massive Be in limiters, Be coatings in recessed areas). In laboratory experiments, a saturation concentration of ∼ 0.3 D/Be in the implantation zone after D implantation at room temperature have been reported in [22,23]. Such saturation concentration remains constant up to a specimen exposure temperature of 400–500 °C [24]. Therefore, the dynamic retention could have been enhanced by considering a first-wall in Be which can house traps with a surface density four times higher than the W which was considered in the present simulation. 3) Eventually, the results question the assumption of =R1 m in the first-wall, which can lead to an overestimation of the main SOL ionization source in edge-plasma simulation. In our simulation, after two seconds of H-mode plasma, Rm can range between 0 and 0.8 depending on the local ion and atom conditions. Acknowledgements This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. This work was granted access to the HPC resources of Aix-Marseille Université financed by the project Equip@Meso (ANR10-10-EQPX-29-01) of the program ‹‹ Investissements d’Avenir ›› supervised by the Agence Nationale de la Recherche. Supplementary material Supplementary material associated with this article can be found, in the online version, at 10.1016/j.nme.2019.03.019 References [1] A.S. Kukushkin, et al., Physics requirements on fuel throughput in iter, J. Nucl. Mater. 415 (1, Supplement) (2011) S497–S500, https://doi.org/10.1016/j.jnucmat. 2010.08.050. Proceedings of the 19th International Conference on Plasma-Surface Interactions in Controlled Fusion [2] C. Grisolia, Plasma wall interaction during long pulse operation in tore supra, J. Nucl. Mater. 266–269 (1999) 146–152, https://doi.org/10.1016/S0022-3115(98) 00815-0. [3] S. Wiesen, et al., Impact of the jet iter-like wall on h-mode plasma fueling, Nucl. Fusion 57 (6) (2017) 066024. [4] H. Bufferand, et al., Numerical modelling for divertor design of the west device with a focus on plasma-wall interactions, Nucl. Fusion 55 (5) (2015) 053025. [5] V. Philipps, et al., Dynamic fuel retention and release under iter like wall conditions in jet, J. Nucl. Mater. 438 (2013) S1067–S1071, https://doi.org/10.1016/j. jnucmat.2013.01.234. Proceedings of the 20th International Conference on PlasmaSurface Interactions in Controlled Fusion Devices [6] T. Loarer, et al., Gas balance and fuel retention in fusion devices, Nucl. Fusion 47 (9) (2007) 1112. [7] E.A. Hodille, et al., Macroscopic rate equation modeling of trapping/detrapping of hydrogen isotopes in tungsten materials, J. Nucl. Mater. 467 (2015) 424–431, https://doi.org/10.1016/j.jnucmat.2015.06.041. [8] A.H.M. Krom, et al., Hydrogen trapping models in steel, Metall. Mater. Trans. B 31 (6) (2000) 1475–1482, https://doi.org/10.1007/s11663-000-0032-0. [9] J.F. Ziegler, et al., Srim - the stopping and range of ions in matter (2010), Nucl. Instrum. Methods Phys. Res. Sect. B 268 (11) (2010) 1818–1823, https://doi.org/ 10.1016/j.nimb.2010.02.091. 19th International Conference on Ion Beam Analysis [10] N. Ozisik, Heat Conduction - Second Edition, (1980). [11] D. Maillet, et al., Thermal Quadrupoles: Solving the Heat Equation through Integral Transforms, (2000). [12] J. Denis, et al., Wall surface temperature calculation in the soledge2d-eirene transport code, Phys. Scr. 2016 (T167) (2016) 014073. [13] A. Widdowson, et al., Overview of the jet iter-like wall divertor, Nucl. Mater. Energy 12 (2017) 499–505, https://doi.org/10.1016/j.nme.2016.12.008. Proceedings of the 22nd International Conference on Plasma Surface Interactions 2016, 22nd PSI [14] A. Widdowson, et al., Overview of fuel inventory in jet with the iter-like wall, Nucl. Fusion 57 (8) (2017) 086045. [15] E.A. Hodille, Study and modeling of the deuterium trapping in ITER relevant materials, UNIVERSITE D’AIX-MARSEILLE, 2016 Theses. [16] N. Fernandez, et al., Hydrogen diffusion and vacancies formation in tungsten: density functional theory calculations and statistical models, Acta Mater. 94 (2015) 307–318, https://doi.org/10.1016/j.actamat.2015.04.052. [17] B. Pégourié, et al., Deuterium inventory in tore supra: coupled carbon-deuterium balance, J. Nucl. Mater. 438 (2013) S120–S125, https://doi.org/10.1016/j. jnucmat.2013.01.019. Proceedings of the 20th International Conference on PlasmaSurface Interactions in Controlled Fusion Devices [18] V. Philipps, et al., Analysis of outgassing after joint european torus discharges under beryllium first wall conditions, J. Vacuum Sci. Technol. A 11 (2) (1993) 437–445, https://doi.org/10.1116/1.578750. [19] S. Brezinsek, et al., Fuel retention studies with the iter-like wall in jet, Nucl. Fusion 53 (8) (2013) 083023. [20] T. Loarer, et al., Comparison of long term fuel retention in jet between carbon and the iter-like wall, J. Nucl. Mater. 438 (2013) S108–S113, https://doi.org/10.1016/ j.jnucmat.2013.01.017. Proceedings of the 20th International Conference on Plasma-Surface Interactions in Controlled Fusion Devices [21] R. Bisson, et al., Dynamic fuel retention in tokamak wall materials: an in situ laboratory study of deuterium release from polycrystalline tungsten at room temperature, J. Nucl. Mater. 467 (2015) 432–438, https://doi.org/10.1016/j.jnucmat. 2015.07.028. [22] W.R. Wampler, Retention and thermal release of deuterium implanted in beryllium, J. Nucl. Mater. 123 (1) (1984) 1598–1602, https://doi.org/10.1016/00223115(84)90310-6. [23] H. Kawamura, et al., Retention of deuterium implanted in hot-pressed beryllium, J. Nucl. Mater. 176–177 (1990) 661–665, https://doi.org/10.1016/0022-3115(90) 90123-5. [24] R.A. Anderl, et al., Hydrogen isotope retention in beryllium for tokamak plasmafacing applications, J. Nucl. Mater. 273 (1) (1999) 1–26, https://doi.org/10.1016/ S0022-3115(99)00022-7. J. Denis, et al. Nuclear Materials and Energy 19 (2019) 550–557 557