Boundary conditions at the limiter surface obtained in the modelling of plasma wall interaction with a penalization technique A.Paredes, E. Serre, L. Isoardi, G. Chiavassa, G. Ciraolo, F. Schwander, Ph.Ghendrih*, Y. Sarazin*, P. Tamain* M2P2 UMR 6181, CNRS - Aix-Marseille Université, * CEA IRFM, 13108 Saint Paul-Lez-Durance, France.

Overview Summary The ANR project ESPOIR is concerned by a simulation effort of the plasma wall interaction based on first principle physics. In this frame, Isoardi et al. (JCP 2010) proposed a penalization technique to model solid plasma facing components that treats solid obstacle as a sink region corresponding to the strong plasma recombination in the solid state material. A major advantage of such a description is that the system can be solved in an obstacle free domain that allows the use of powerful numerical algorithms. Such a technique implemented in a minimal transport model for density and momentum appeared to exhibit a Mach-1 transition between the presheath and the limiter region. By analysing physics of detached plasmas that are governed both by strong recombination and plasma pressure decrease, as imposed by the penalization technique within the limiter region. we provide a unique control parameter A=Γ Γ csmi/Π Π that allows one to understand the results of the penalization technique for the Mach-1 transition

Introduction In numerical tools for solving three-dimensional fluid-like equations, surface boundary conditions on the plasma facing components are demanding, both to implement the various shapes and regarding the required CPU time. Bohm boundary conditions are usually used that constrain the particle flux to a sonic flow and determine the plasma energy flux lost by the plasma. To overcome this difficulty, we have proposed in ref. [1] an original penalization technique in the frame of a reduced model of transport equations for the ion density and the particle flux that allows one to treat the limiter as a region of the plasma governed by a very strong particle sink. In such system, the plasma-limiter interaction is no longer described by demanding to density and velocity to satisfy Bohm conditions. Penalized terms are added to density and momentum conservation equations such that the density and the parallel momentum are penalized toward zero inside the obstacle in order to mimic a sink for particles corresponding to the strong plasma recombination in the solid state material. The new system is handled now in an obstacle free domain for which fast solvers are available. Moreover, this technique allows an easy modification of the obstacle geometry and the associated boundary conditions. A very attractive feature showed in ref. [1] was that solutions provided a plasma velocity which is almost sonic at the boundaries obstacles as expected from the sheath conditions through the Bohm criterion. This result is surprising at first glance since the model does not introduce any departure from quasineutrality, the latter being the key ingredient of the sheath region that governs the Mach-1 criterion as plasma surface boundary condition. Analysing the physics of detached plasmas that are governed both by strong recombination and plasma pressure decrease, as imposed by the penalization technique within the limiter region allows one to understand the results of the penalization technique.

Physical modelling The plasma is addressed as a fluid. Let us consider the 1-D physics of the asymptotic matching of the presheath. The key physics are governed by the balance equations for particles and total plasma momentum, namely:

∂ t n + ∂ sΓ = S n

Sn and SΓ 1 (1) ∂ sΠ = SΓ source terms. mi This set of equations are governed by the variation of the particle flux and of the total plasma pressure:

∂tΓ +

TEM PLATE DESIGN © 2 0 0 8

www.PosterPresentations.co m

Γ = nMcs

(

Π = n(Te + Ti ) 1 + M 2

(2)

)

mi cs2 = (Te + Ti )

s is the curvilinear coordinate aligned along the plasma flow (close to the magnetic field direction since the plasma velocity transverse to the field lines is small in most situations of interest). Also, the momentum exchange between transverse and parallel directions is not explicitly taken into account in Eq.(1), although one could consider that this physics is included in the source term SΓ.

Penalization technique Introduced here in 1D for simplicity but 3D extension straigthforward. The limiter is modeled as a particles sink using a penalized system for Eq. (1) with Sn = S a localized density source and SΓ = 0. The limiter (or obstacle) Ω is represented by a mask function satisfying χ(s) = 1 if s in Ω and χ(s) = 0 elsewhere. Periodic bc are considered at s = 0 and s = 1. (n, Γ) depend now on the penalization parameter which must be small (η<< 1) in the following system :

∂ t nη + ∂ sΓη = (1 − χ )S − ∂ t Γη + (1 − χ )

1

η

- Plasma detachment with momentum losses (SΓ<0): Eq. (8) becomes for variations of the particle flux and plasmas pressure governed by sΓ=-νΓn, sΠ=-νΠmiΓ

(

(

)

FIG 1: One-dimensional profiles of the numerical (squares) and analytical (solid line) solutions for density (left) and parallel momentum (right) (b). Vertical dotted lines indicate the limiter location.

- The ratio M = Γ/n perfectly matches MΩ on Fig. 2 for two functions MΩ (see Figure caption). - The solution in the plasma is not significantly influenced by MΩ - In particular, M ± Bohm-like conditions find at the limiter with a relatively reduced sensitivity to the chosen function M.

( )

By introducing (4) in Eq.(3 ) and by identifying the terms: order η−1 order η0 ∂ t n + ∂ s Γ = (1 − χ )S − χn~ χn = 0 (6) (5) 1 ~ χ (Γ − nM Ω ) = 0 ∂ t Γ + (1 − χ ) ∂ s Π = χ Γ − nM Ω mi

(

)

- In the plasma domain, χ = 0 and Eq.(6) reduces to the initial Eq.(1) with Sn = S and S Γ = 0. - Inside the limiter, χ = 1, and Eq.(5) becomes: n = 0, Γ = nMΩ n = 0 and Γ= 0 in the limiter, i.e, the physical plasma conditions in the limiter. Numerical method Strong gradients in the vicinity of the obstacle requires the use of a robust shock-capturing finite difference method. - numerical fluxes : shock capturing scheme based on the Roe approximate Riemann solver [3], associated with a robust third order ENO spatial reconstruction [4]. -Time integration of Eq.(3) : second order TVD Runge-Kutta scheme [4] with an explicit treatment of the convective terms and an implicit treatement of the penalization terms

Analytical analysis Without loss of generality, time variations are included in generalised source terms so that steady state like equations are obtained from Eq. (1) with the new source terms Sn = Sn − tn and SΓ = SΓ − tΓ: (6)

[

FIG 2: Magnification of the solution profiles around the limiter showing that the Mach-1 condition is robust to MΩ.. M profiles (*) for (a) linear MΩ and (b) hyperbolic tangent MΩ (solid lines). Vertical dotted lines indicate the limiter location

3D simulations The tokamak is considered with a cylindrical poloidal cross section of radius r, poloïdal angle and toroïdal angle . A 3D limiter is defined by (r, , ) = 1 if (r, , ) in Ω and (r, , ) = 0 if (r, , ) elsewhere. Magnetic field lines intersecting or not the limiter are illustrated in the (θ, ϕ) plane in Fig. 3.

FIG 3. Schematic representation of the magnetic field lines together with the 3D limiter in the (θ, ϕ) plane. All the lines intersect the limiter except the pink one.

Eq. (1) is generalized to 3D with source terms containing radial diffusion mechanisms affecting both (n, Γ). A negative density gradient is fixed at the core boundary. Finite difference (ENO)-Fourier approximations are used

]

1 1 + ε (1 − A² )1/ 2 A Π 1 − ε (1 − A² )1 / 2 n= (Te + Ti ) M=

[

]

(7)

as a function of a unique control parameter A = 2Γ Γ mics/Π Π varying between −1 and 1 ε=1

supersonic branch

ε=-1 subsonic branch (plasma presheath)

- In the standard presheath with steady state conditions, one finds that: A is monotonic the plasma remains on the subsonic branch - Access to the supersonic branch requires a reversal of the variation of A: achieved on the Debye scale for a non-neutral plasma such that the sheath electric field governs an increase of the plasma pressure. variations of A govern the n-variations as well as the bifurcations between the subsonic and supersonic branches at A²= 1.

∂n = +∞ ∂A

∂M = +∞ ∂A

FIG 6: Subsonic and supersonic Mach and density profiles as a function of A.

1 1 1 1 1 ∂ s A = ∂ sΓ − ∂ sΠ + ∂ s (Te + Ti ) A Γ Π 2 (Te + Ti )

FIG 4: Density and parallel momentum profiles in the parallel direction along a magnetic field line intersecting the limiter (red line of Fig 3.)

A=1

transition to supersonic branch will occur upstream from the limiter Α monotonic transition to supersonic branch requires the sheath electric field

In the region with negative particle source ( sΓ<0), subsonic branch provided that :

νΠ

sin 2 (a ) ≥ β

sA>0

on the

Transition to a supersonic flow in the presheath must therefore be associated to a large temperature variation.

- Penalization In the presheath one considers Sn > 0 and SΓ 0 A is monotonic (from Eq. 8). One observes the well-known plasma density decrease towards the limiter (see Fig 1). Conversely, one obtains an increasing Mach number that is consistent with a constant plasma pressure. In the limiter region, Sn < 0 (particle sink). Particle sink enforces a decreasing n as in the presheath together with a reversal of M variation. On the subsonic branch, reversing the variation of M requires to reverse the variation of the control parameter A (M, A) exhibit a maximum. However, n must switch from being a decreasing function of A to an increasing function of A. As readily seen of Fig. 6 this transition can only occur at the point A = 1 and therefore at M = 1.

Conclusions

-Penalization technique recovers bc similar to that governed by the asymptotic matching of the presheath conditions to the Bohm condition within the sheath. -This asymptotic matching results from the bifurcation from the subsonic to the supersonic regime (n, M): only achieved on the Debye scale where the departure from quasineutrality reverse the variation of A. -In penalization, the different physics between the presheath and the limiter the reverse variation of A. Continuity of n and of its derivative in the limiter boundary layer can then only be achieved by stepping through the same bifurcation point for n. The plasma solution then remains subsonic (Fig. 6) and exhibits both a (n, M) decreasing in the limiter region as required by the penalization.

Numerical results show the expected behaviour for n and Γ (Fig. 4, 5), in particular a Mach number close to 1 at the limiter boundary.

For convenience, in the following we shall consider the case A hence Γ 0. According to the definition of A:

=0

1 ∂s A > 0 A

2ν Γ

Sn and SΓ source terms.

From Eq. (2), M and n are obtained as roots of second-order polynomials:

Numerical set up Initial conditions n(s, t = 0) = 1 and Γ(s, t = 0) = 0. Limiter size is set at s = 0.1, and particle source S=2 Grid points number m = 100. Results - Excellent agreement with the analytical solution of ref.[1] on Fig. 1. - Inside the limiter, (n, Γ) 0, following the prescribed behavior.

FIG 5: Iso-surface of zero density colored by the Mach number in a tokamak with a circular poloidal cross section. The grey surface shows a 3D limiter modelled using penalization.

1 ∂ s Π = SΓ mi

s ∈ [0,1], t ∈ [0, ∞]

The penalized terms in Eq.(3) enforce the initial boundary conditions on density and momentum: ~ ~ nη = n + ηn~, Mη = M + ηM Γη = Γ + ηΓ + O η 2 (4)

1 ∂s A A

∂sΓ = Sn

(3)

ΓΩ bounds the ratio Γ/n inside the limiter for numerical stability since ΓΩ / n=MΩ where MΩ is imposed inside the limiter (non identically zero) and is zero outside. For the same reason, the pressure force must be canceled within the limiter.

(9)

with β=1+0.5 sTΓ/T sΓ

nη

1 1 ∂ s Πη = − Γη − ΓΩ η mi

)

1/ 2 1 1 ν ∂ s A = ∂ s Γ β − Π 1 + ε 1 − A2 4ν Γ Γ A

0

(8)

- Plasma detachment with constant temperature and SΓ~0: Transition to supersonic branch: s A=0 sΓ=0 namely the ionisation front where there is the transition from a positive to a negative particle source governed by ionisation

- In detached plasmas, with strong recombination rates and subject to strong plasma pressure decrease, transition to a supersonic presheath can only occur at the ionisation front (Sn changes sign), provided the plasma momentum loss is relatively small with respect to the particle source on either sides of the ionisation front

References

[1] Isoardi L., Chiavassa G., Ciraolo G., Haldenwang P., Serre E., Ghendrih Ph., Sarazin Y., Schwander F. & Tamain P. J. Comp. Phys., J. Comp. Phys., 229(6), 2220-2235, 2010 [2] LeVeque, Finite, Cambridge University Press, (2002) pp 130-132 [3] Shu C.W. & Osher S.J. , J. Comput. Phys., 83 (1989), pp. 32-78.

Acknowledgements Financial support from ANR is acknowledged, project ANR-09-BLAN-0035-01, http://sites.google.com/site/anrespoir/