T Stindl
email: stindl@irs.uni-stuttgart.de

J Neudorfer

A Stock

M Auweter-Kurtz

C.-D Munz

S Roller

Rainer Schneider

Comparison of coupling techniques in a

published or not. The documents may come

Introduction

Numerical simulations of highly rarefied, fully ionized non-neutral plasma flows in complexshaped large-scale devices are playing a central role in quite different areas of applied science and technology. The required kinetic description of such flows is given by the time-dependent Vlasov equation in six-dimensional phase space and the full system of Maxwell equations. A powerful numerical method to tackle the non-linear Maxwell-Vlasov problem is the Particle-In-Cell (PIC) approach [START_REF] Hockney | Computer Simulation using Particles[END_REF][START_REF] Birdsall | Plasma Physics via Computer Simulation[END_REF], a technique with a long history over more than five decades [START_REF] Birdsall | [END_REF] which provides a statistical representation of the particle distribution functions in phase space. This numerical technique allows to retain most of the basic phenomena because no physical approximations to the fundamental equations are necessary, and collective effects are included in a self-consistent manner by coupling the charged particle distributions to the field equations via source terms. Within the PIC framework, the charged particle distributions are advanced in grid-free phase space according to the Lorentz equations -the characteristic theory solution of the Vlasov equation. For this purpose a multitude of well-understood schemes are available (see, e.g. [4] and references given therein). Furthermore, the Maxwell equations account for the evolution of the electromagnetic fields and are usually solved on a computational mesh in space. The majority of available PIC codes rely on the finite difference scheme proposed by Yee [5] for the Maxwell equations using a uniform staggered grid in space and time. However, relevant devices in practical applications are often complex shaped -like, for instance, the resonators of advanced gyrotrons -and thus regular computational grids are not well suited. Other approaches are based on finite elements as well as finite volume techniques which allow the use of unstructured grids (see, e.g. [6,7]), but are mostly second order accurate in both space and time and, hence, require very high resolution for practical problems. However, on fine computational meshes the dispersion properties of the mentioned techniques are not suitable to cover the electromagnetic wave propagation for long times correctly. A promising way out of these difficulties are field solvers based on discontinuous Galerkin and conforming finite element methods which have been proposed recently [8,9,[START_REF] Jund | IRMA Strasbourg[END_REF]. These new high-order approaches for PIC codes are able to approximate wave propagation over long distances and times without significant dispersion and dissipation errors, and, in contrast to finite difference schemes [START_REF] Greenwood | [END_REF], do not suffer from artificial numerical Cherenkov radiation [START_REF] Jacobs | Simulations of the Weibel instability with a high-order discontinuous Galerkin Particle[END_REF]. All introduced field solvers rely on explicit approximations to the Maxwell equations and thus underlie a (local) stability restriction known as the CFL condition, which links the time step size ∆t with the typical discretization ∆x according to ∆t < ∆x/c, where c denotes the speed of light. Beyond that Maxwell solver dependent restriction, explicit PIC procedures are limited by two additional accuracy conditions (see, e.g. [START_REF] Turner | [END_REF] and references given therein). The first constraint is closely related to the fastest electron response time and reads as ∆t < 2/ω pe , where ω pe is the electron plasma frequency. The second limitation states that the electron Debye length λ De should be resolved by the discretization: ∆x/λ De 1. A violation of this restriction may lead to the "aliasing instability" resulting from interpolation between the particles and the grid [START_REF] Jacobs | Simulations of the Weibel instability with a high-order discontinuous Galerkin Particle[END_REF]. Furthermore, it is advisable to check whether the number of simulation particles per Debye sphere is sufficient to obtain convergent simulation results. Especially, this aspect seems to be very important in the context of collisional particle simulations [START_REF] Turner | [END_REF]. Problems and difficulties arising with the computational grid for the Maxwell equations as well as the coupling of the fields with the charged particles can be circumvented in a fully grid-free framework. Very recently, Christlieb and coworkers proposed such a grid-free plasma simulation technique which is based on the boundary integral / treecode (BIT) method [14]. Here, the homogeneous solution of the Poisson equation is approximated by a boundary integral method while the particular solution of the elliptic equation is computed by a treecode. In that part of the BIT method the particles are subdivided in clusters and the direct particle-particle interaction is replaced by a particle-cluster interaction, which is determined by a multipole expansion of the free-space Green's function. Although the BIT technique is an attractive alternative to grid-based PIC simulations, it is not applicable in our cases, where the transient solution of the Maxwell equations is necessary. In the present paper we focus our attention to particle-mesh coupling techniques which are necessary for the exchange of information between mesh-based and mesh-free numerical models of the PIC code. As we will see, suitable coupling strategies enormously increase the effectiveness and efficiency of the PIC code which is a prerequisite for the simulation of complex applications, e.g. future high power micro-and millimeter wave sources for fusion applications. The structure of this article is as follows: The statement of the governing equations and the Particle-In-Cell frame work is given in Section 2. Afterwards, in Section 3, possible particle-mesh coupling strategies are discussed in more detail. Section 4 deals with simulation results of the Weibel instability obtained with different particle-mesh coupling techniques for deposition and interpolation. In Section 4 we also present results of a gyrotron resonator simulation. Finally, a conclusion and short outlook on future activities is given in Section 5.

The Physical-Mathematical and Numerical Model

The Governing Equations

In the context of simulating the self-consistent movement of charged particles in electromagnetic fields with Particle-in-Cell (PIC) methods it has been recognized [START_REF] Hockney | Computer Simulation using Particles[END_REF][START_REF] Birdsall | Plasma Physics via Computer Simulation[END_REF] that charge conservation is violated at the discrete level by the approximation procedure itself. For this reason it is expedient to model mathematically the electromagnetic wave propagation by the "Generalized Lagrange Multiplier" (GLM) Maxwell system [15,8] instead of the pure Maxwell equations. The basic idea of this approach is to couple the hyperbolic Maxwell equations (Maxwell's and Faraday's law) to the elliptic part (Gauss law and the absence of magnetic monopole statement) by two additional degrees of freedom Φ(x, t) and Ψ(x, t). For given charge and current density, ρ(x, t) and j(x, t), the GLM Maxwell equations in conservation form read as [16] ∂u ∂t

+ 3 β=1 K β ∂u ∂x β = g , (1)
where the state vector u(x, t)

= (E 1 , E 2 , E 3 , Ψ, B 1 , B 2 , B 3 , Φ)
T is given by the components of the electric field E(x, t) and magnetic induction B(x, t) as well as by the two additional variables Ψ and Φ. Furthermore, K β ∈ R 8×8 are block-structured matrices with constant entries reflecting the structure of the GLM system and g = -1 ǫ 0 (j 1 , j 2 , j 3 , 0, 0, 0, 0, -χ ρ) T denotes the source term vector with a constant χ associated with the divergence cleaning speed and the electric permittivity ǫ 0 . Clearly, the GLM Maxwell system has to be supplemented with appropriate initial and boundary conditions (for more details see [16]). The charged particle flow in the computational domain is mathematically fully described by the distribution function f (x, v, t) whose evolution in phase space is given by the Vlasov equation (see, e.g., [START_REF] Balecu | Transport Processes in Plasmas[END_REF]). Note that, for the sake of simplicity, we consider only one particle species of the plasma with charge q and mass m. The solution of the Vlasov equation obtained from the theory of characteristics fulfills the usual laws of classical mechanics (also called Lorentz equations)

dx ν (t) dt = v ν (t) and d dt [γ ν (t) v ν (t)] = q m [E(x ν , t) + v ν (t) × B(x ν , t)] , (2)
where the relativistic factor is given by γ

-2 ν = 1 -(v ν /c) 2 .
Here, the phase space coordinates of the "ν"-th particle are denoted by (v ν , x ν ) and the index ν runs over all simulation particles N P : ν = 1, . . . , N P . Note that in order to obtain the particles' Lorentz equations [START_REF] Birdsall | Plasma Physics via Computer Simulation[END_REF] we have already replaced the continuous particle distribution function f by its weak approximation f (for details, see [START_REF] Munz | [END_REF]). The self-consistent interaction of the simulation particles with the electromagnetic field is established by calculating the first two moments of f yielding the charge and current density,

ρ(x, t) = q N P ν=1 W ν S (x -x ν (t)) and j(x, t) = q N P ν=1 W ν v ν (t) S (x -x ν (t)) , (3)
where the individual particle weights are abbreviated by W ν , and S (x) denotes the B-spline based shape function [START_REF] Hockney | Computer Simulation using Particles[END_REF], which is introduced for renormalization reasons. Note that Jacobs and Hesthaven [8] proposed recently a polynomial shape function for charge and current density assignment which could be an attractive alternative to the B-spline approach.

General Numerical Framework

An attractive numerical method to solve the non-linear Maxwell-Lorentz (ML) model, ( 1)-( 3), is the PIC approach. The peculiarities of the PIC method are the powerful particle-mesh techniques which couple grid-based and mesh-free numerical models (see next section). The basic idea of this approach can be summarized as follows (see Fig. In each time step the electromagnetic fields, obtained by the numerical solution of the GLM Maxwell equations, are interpolated to the actual locations of these particles. According to the Lorentz force, the charges are redistributed and the new phase space coordinates are determined by solving the Lorentz equations numerically. To close the chain of self-consistent interplay, the particles have to be located with respect to the computational grid in order to determine the contribution of each charged particle to the changed charge and current density, which are the sources for the Maxwell equations in the subsequent time step.

Maxwell Solver: For the current ML simulation code a high-order discontinuous Galerkin (DG) method [START_REF]Discontinuous Galerkin Methods: Theory, Computation and Application[END_REF][START_REF] Hesthaven | Nodal Discontinuous Galerkin Methods[END_REF] was implemented for the numerical solution of the GLM Maxwell system (1).

The DG method has the advantage of maintaining a local scheme even for simulations of high convergence order. While for a Finite Volume method [START_REF] Fertig | Hybrid code development for the numerical simulation of instationary magnetoplasmadynamic thrusters High Performance Computing in Science and Engineering[END_REF] the reconstruction of a high-order solution using the cell mean values of neighboring elements includes information from the direct neighbors as well as more distant elements, the discontinuous Galerkin method uses local basis functions, eliminating the need for information from any other than the direct neighboring elements. This highly local approach even for high order simulations results in a scheme with an efficient parallelization [9]. The discontinuous Galerkin method has been investigated before in the two-dimensional PIC context by Jacobs and Hesthaven [8], who used a purely nodal DG method [START_REF] Hesthaven | Nodal Discontinuous Galerkin Methods[END_REF] with low storage explicit Runge-Kutta (LSERK) time stepping by Kennedy et al. [START_REF] Kennedy | [END_REF]. In contrast the present scheme uses a mixed nodal and modal DG method in space presented in [23] and a LSERK method in time. Originally, the current DG scheme was designed for a fully unstructured discretization of the three-dimensional domain using all types of elements (hexahedrons, tetrahedrons, pyramids and prisms) in the same mesh. However, for PIC we mainly use tetrahedrons and hexahedrons. In each element k the approximate element-local solution u k h (x, t) of the analytical one u k (x, t) takes the form

u k (x, t) ≈ u k h (x, t) = M i=1 ũk (x i , t)L i (x) = N n=1 ûk n (t)φ n (x) (4)
with the multivariate Lagrange basis {L i } i=1,...,M and the solution ũ(x i , t) evaluated at M nodal interpolation points (IPs) x i and the modal basis {φ n } n=1,...,N with the modal coefficients û(t) depending only on the time. The connection between ũ and û is given by

V û = ũ , where V i,j = φ j (x i ) (5)
is the Vandermonde matrix. In three-dimensions each element with the polynomial degree p needs N = (p + 1)(p + 2)(p + 3)/6 modal degrees of freedom (DOF). The number of IPs M does only equal N in case of linear tetrahedrons. For all other types of elements M > N . The choice of M and x i is discussed in [23]. Figures 2 and3 show the position of the nodal IPs for p = 2 elements. Obviously an approximation of Eq.( 1) by Eq.( 4) leads to a residual. Thus the

k ∂u k h ∂t φ j dx + 3 β=1 K β k ∂u k h ∂x β φ j dx ! = k gφ j dx, (6)
requiring the residual to be projected orthogonally (and thus optimally) onto the function space P p of the test functions φ j . A spatial integration by parts gives the weak form,

k ∂u k h ∂t φ j dx + ∂k F h • nφ j ds - 3 β=1 K β k u k h ∂φ j ∂x β dx = k gφ j dx , (7)
where F h denotes the intercell flux and n is an outward pointing unit vector defined at the boundary of the element. With the matrix operators

M ij = k φ i φ j dx = I, S β,ij = k L i ∂φ j ∂x β dx, E ij = ∂k L i φ j ds,
we can write the fully explicit scheme

I ∂ ûk h ∂t + E Fh • n - 3 β=1 K β S β ũk h = V g. ( 8
)
The mixed nodal and modal approach allows to use nodal values for the sources (ρ, j) and the intercell fluxes F h . The integration is done by Gauss quadrature in an initial phase of the computation, yielding a quadrature free approach.

Lorentz Solver: Two different time discretization methods for particle pushing in phase space have been implemented. The first one is the explicit time-centered leapfrog method introduced by Buneman [24], which is implemented according to the formulation proposed by Boris [START_REF] Boris | Proc. 4th Conf. on Num. Sim. of Plasmas[END_REF]. This scheme is second-order accurate and guarantees long-term phase space stability. The other time integrator is a low-storage fourth-order Runge-Kutta method [START_REF] Kennedy | [END_REF] enabling highorder consistency with the DG solver for the GLM Maxwell equations. Note that the memory usage of this approach is reduced by half compared to the classic fourth-order Runge-Kutta scheme. In general, the Boris leapfrog scheme is more stable, but the Runge-Kutta method gives a higher accuracy. The presented examples are computed only with the leapfrog method due to the lower computational effort compared to the RK scheme.

Strategies for Particle-Grid Coupling

According to the PIC methodology, the charged particles move freely in space while the electromagnetic fields are calculated on a computational grid (hereafter called DG grid; abbreviated by DGG). For the DG grid, structured and unstructured meshes can be used that decompose the computational domain, in general, in elements of different size. Depending on the order of the DG scheme, the interpolation points (IPs) are located at the edges, faces and the interior of each element (see Figs. 2 and3). Consequently, for the electromagnetic field calculations, IP-based charge and current densities have to be provided. For this purpose, the charge and current information of each charged simulation particle needs to be deposited onto the IPs of the elements of the computational grid. Vice versa, after the electromagnetic fields are determined on the IPs, this information has to be interpolated to the actual position of the particles. Note that the DG solver for the GLM Maxwell equations is especially designed to allow for high order accurate approximation on structured as well as general unstructured grids with polymorphic grid cells, i.e. tetrahedra, hexahedra, pyramids or prisms cells are treated in a generalized fashion, where only the IPs are important. Furthermore, the choice of suitable computational grids (structured, unstructured) as well as discretization elements depends strongly on the problem at hand, which is well-known, for instance, from boundary layer simulations.

The Single-Grid Approach

This strategy relies only on the DG grid and assumes that the particles are located within this, in general, unstructured grid [6]. The most basic, zero-order technique for deposition of charges and currents is the so-called cell mean value (CMV) approach. It consists of summing up the charge and current of all particles within a grid cell and applying the corresponding densities uniformly to all grid nodes (IPs) within a cell (see Fig. 4 left, for the two-dimensional situation). While this method is fast, it requires a very large number of simulation particles to avoid strong oscillations in the source terms. The inverse procedure of the CMV approach, namely to interpolate the information in the IPs to the particle position, consists in evaluating the given DG polynomials at the barycenter of each grid cell and assigning these electromagnetic fields to all particles within the tetrahedrons. However, depending on the spatial discretization, this approach can be highly inaccurate. In order to reduce the number of simulation particles but obtain accurate results, higher order methods for particle-grid coupling are necessary.

The Auxiliary Mesh Coupling Technique

A different technique for the deposition of charges and the interpolation of electromagnetic fields is obtained by superimposing a Cartesian background mesh (CBM) to the computational DGG (see right plot of Fig. 4). This enables high-order deposition and interpolation, which lead to increased accuracy and thus a reduction of the number of necessary simulation particles. The main idea is to split the intrinsic "particle-DGG coupling" into a strictly "particle-CBM coupling" and into an information exchange between the CBM and the DGG. One big advantage of this strategy is that the simulation particles can be very efficiently localized with respect to the Cartesian background mesh [START_REF] Westermann | [END_REF]. In order to deposit the charge and current density generated by the charged particles to the nodes of the CBM and, conversely, to interpolate the electromagnetic fields evaluated at the nodes of the CBM cells to the particles in the interior of these cells, we apply in the present context B-splines of the orders zero, one and three. These Left: Deposition via cell mean value method using the DG computational grid directly. Charge current densities of the particles are deposited uniformly within a cell. Center: Deposition from a particle onto the CBM and, vice versa, interpolation from the nodes of the CBM to the particle via first-order B-splines (area or rather volume weighting), where the weights according to the areas (volumes) between the particles and the mesh nodes are used. Right: Coupling of the CBM with the computational DG grid using area (volume) weighting. are constructed recursively and could easily be extended to two-and three-dimensional requirements [START_REF] De Boor | A Practical Guide to Splines[END_REF]. Note that zero-order B-splines in three dimensions are similar to the above mentioned CMV approach for unstructured grids and the first-order B-splines in two and three dimensions are known in literature as the area (see center plot of Fig. 4) and volume weighting method, respectively. The charge and current density available on the nodes of the CBM have to be provided on the IPs of the DGG. To interpolate the source terms from the nodes of a CBM grid zone to the IPs inside this cell, the B-spline approach is once again applied. Note that for the numerical examples presented below we used three-dimensional first-order B-splines, i.e. volume weighting. Since the location of the IPs is fixed at all times, the weights of the approximation points are computed only once during initialization. To obtain the electromagnetic fields at the particle positions after their computation by the Maxwell solver, it is possible to evaluate the fields at all actual particle positions. Thus, using the CBM it is possible to drastically lower the computational costs by computing the fields from the DG polynomials at the fixed nodes of the CBM and then interpolating them to the particle positions using the same B-spline approach. Moreover, the parameters for the evaluation of the DG polynomials at the nodes of the Cartesian background mesh can be retained for all time steps, leading to an additional decrease in computational costs. In classic coupling approaches the field solver grid is used to resolve field specific phenomena and particle driven phenomena, i.e. the length of electromagnetic waves and the Debye length. Usually both scales are not equal to each other and the resulting grid is a compromise. The size of the CBM can also be chosen w.r.t. the resolution of the particle driven phenomena, decoupling the above mentioned scales. In this context the CBM size has an impact on the resolution of certain particle driven phenomena in a plasma and thereby on the overall accuracy of the scheme.

Numerical Experiments

In the following the proposed particle-grid coupling strategies are numerically investigated with the fully relativistic, electromagnetic PIC solver, which is based on the time-dependent Maxwell-Lorentz system in six-dimensional phase space. For this, two examples are considered: At first, simulation results for the Weibel instability are discussed. Secondly, numerical computations of high power microwave generation in a gyrotron are presented.

Weibel Instability

In 1959 Weibel [28] has shown that an anisotropic electron velocity distribution

f (0) e (v) = 1 π 3/2 v v2 ⊥ exp - v 2 1 v2 - v 2 2 + v 2 3 v2 ⊥ , v2 = 2 k B T e m e , v2 ⊥ = 2 k B T ⊥e m e , (9)
with longitudinal and transversal thermal velocities, v and v⊥ , in a homogeneous, collisionless plasma supports unstable growth of electromagnetic waves. Since that time, the Weibel instability has been comprehensively with analytical [29,30] and numerical methods [31,8] and, very recently, applied for astrophysical scenarios of gamma-ray bursts and pulsar winds [32]. The basic characteristic of the kinetic Weibel instability is that "small scale" B-fields and current inhomogeneities are amplified and free energy the electron stream is converted to a "quasi-static" magnetic field (see Figure 5). In order to simulate the Weibel instability we use a computational domain of size 15 c ωpe ×15 c ωpe × c ωpe , where c and ω pe denote the speed of light and the electron plasma frequency, respectively. The simulation time is expressed in multiples of ω -1 pe and all energies shown normalized by the total initial energy. For all computations, a second-order accurate Maxwell solver, i.e. DG basis functions of degree p = 1, is chosen and periodic boundary conditions are applied. Initially, electron-ion pairs are randomly distributed in the computational domain, where the velocities of the electrons are drawn from the bi-Maxwellian distribution function (9) with v = 0.25 c and v⊥ = 0.05 c , while the ions are at rest (and mostly stay at rest during the simulation due to their comparatively high mass). The first numerical experiment is performed on a very fine computational DG grid with the CMV coupling technique and serves as a reference in the following. For this, the computational domain is discretized by 300 × 300 × 20 hexahedral grid cells, where the cell size is roughly the Debye sphere. Furthermore, since the CMV approach is used, 3.33 • 10 7 electron-ion pairs are distributed in the computational domain (which amounts to 37 particle per grid cell) to ensure smooth charge and current densities. The results of this simulation are depicted in Figure 5 (red full lines), where the time dependent evolution of the magnetic (a) and electric (c) energy stored in the electromagnetic fields and the kinetic (b) and total (d) energy is plotted. For the next numerical experiments, the number of simulation particles is drastically reduced by a factor of 37 to 9 • 10 5 particles per species without changing the overall charge and current density. Keeping the number of particles per grid cell constant at 37 particles per cell, the spatial resolution of the computational DGG is reduced to 90×90×6 hexahedral cells without changing the domain size. First, we repeat for this setup the simulation of the Weibel instability with the CMV coupling technique. The results of this computation are also plotted in Figure 5 as blue dashed lines. Afterwards, we perform simulations with the auxiliary mesh coupling approach by applying three-dimensional B-splines for assignment and interpolation and use a Cartesian background mesh of the same size as the DGG. For comparison purposes, the histories of the electromagnetic fields, the kinetic and total energy obtained for the first-(blue dashed-dotted) and third-(black dashed-dotted-dotted) order splines are also shown in Figure 5. The Figures 5a and5b and t ≈ 70 ω -1 pe is accompanied by a strong fall of the kinetic energy. While the simulation results obtained with the auxiliary mesh coupling technique and the CMV approach with low resolution show a pronounced peak in the B-field energy, the reference solution remains on the high energy level without a significant decrease after its maximum. It should be noted, though, that the peak in the magnetic energy curves of the simulations using a CBM better reflect the results of Morse and Nielson [31] and Jacobs et al. [START_REF] Jacobs | Simulations of the Weibel instability with a high-order discontinuous Galerkin Particle[END_REF]. The electric energy level of the simulations should be negligible since the plasma is quasi-neutral. As shown in Figure 5c, the CMV simulation for the low resolution grid and small number of simulation particles reveals a distinct electric energy level. In contrast, the first and thirdorder B-spline coupling method with CBM as well as the reference computation show an almost negligible electric energy level. Furthermore, the reference computation which uses the CMV coupling method shows a constant energy dissipation (Figure 5d) while the loss in total energy for the first or third-order spline coupling methods via the CBM is considerably smaller and occurs mainly in the process of the energy conversion around t ≈ 65 ω -1 pe . Jacobs et al. [START_REF] Jacobs | Simulations of the Weibel instability with a high-order discontinuous Galerkin Particle[END_REF] found that the remaining decrease of the total energy is due to the damping of the divergence error and can be eased by increasing the divergence cleaning speed. The effects of the size of the CBM have been studied, confirming the phenomena described by Jacobs et al. [START_REF] Jacobs | Simulations of the Weibel instability with a high-order discontinuous Galerkin Particle[END_REF] now for the use of the CBM. The number of particles per CBM cell generally increases the amount of energy converted from kinetic to magnetic energy. In contrast, very small numbers of particles per CBM cell generally introduce artificial electric energy through sampling errors. This results in poor conservation of the total energy for fine CBMs. For very fine CBMs with less than 5 particles per CBM cell, simulations with the volume weighting method became unstable while the use of third-order B-spline coupling showed stable operation even in this regime. This is due to the larger stencil of the higher order B-spline coupling. For the simulations shown here, the CBM discretized the domain in 100 × 100 × 6 grid cells. With this setup, the grid cells of CBM and DGG are about the same size with few corners of CBM and DGG being co-located. The CBM coupling technique using higher order B-splines not only yields results that compare favorably to the high resolution reference but also leads to a decrease in computational costs due to the considerable reduction of the number of simulation particles and the number of grid cells. Compared to the reference simulation, the low resolution CBM-based computations using B-Splines require only about 1.5% of the simulation time.

Resonator of a Gyrotron

The application of the high-order particle-grid coupling methods to the simulation of a gyrotron resonator is currently being examined. A schematic drawing of the gyrotron is shown in Figure 6. A gyrotron is a device for high power microwave generation. The emitted radiation is generated in the resonator by bunching electrons with cyclotron motion in a strong external magnetic field. Physically, the reason for the excitation of microwaves in the resonator is due to the relativistic electron cyclotron resonance maser instability. The geometry of the resonator used for our test case is shown in Figure 7. The setup is an azimuthal symmetric standard test case for gyrotron simulation presented in [START_REF] Illy | Untersuchung von Strahlinstabilitäten in der Kompressionszone von Gyrotron-Oszillatoren mit Hilfe der kinetischen Theorie und Zeitabhängiger Particle[END_REF]. A current of 15 A is simulated with 10 5 macro particles generating a stationary, azimuthal symmetric T E 0,3,1 mode. The grid consists of ∼ 4.5 million tetrahedral elements using a second-order DG space approximation, yielding ∼ 18 million DOFs. Simulation results, recorded at t = 1.6 ns, are depicted in Figure 8 and are obtained from a highly parallel computation on 128 cores of the "Nehalem" cluster (Intel Xeon X5560, 2.80 GHz, 8MB Cache, FSB 1333 MHz) at the High Performance Computing Center Stuttgart. Qualitatively the CMV method (left plot) shows a poor evolution of the modes whereas the volume weighting approach (second plot) indicates a slightly better shaping in the modes. The best result is obtained by the third-order B-splines (third plot) with a very good agreement with the comparative simulation (right plot). As a numerical artefact of the low order CMV deposition we can see crinkled fields where the particle beam is located. For the first-order splines the fields are getting smoother in Figure 8: Contour plots of a slice (y = 0) through the gyrotron resonator at 1.6 ns simulation time. Results from simulations with the CMV approach (left) and different CBM-based coupling methods (first-and third-order B-splines; second and third plot, respectively) are compared to a two-dimensional PIC simulation with rotational symmetry (plot on the right by courtesy of S. Illy [START_REF] Illy | Untersuchung von Strahlinstabilitäten in der Kompressionszone von Gyrotron-Oszillatoren mit Hilfe der kinetischen Theorie und Zeitabhängiger Particle[END_REF], note that the coordinates start at -0.02).

the beam region and for the third-order method we have nearly perfectly smooth shaped fields. Additionally we examine the energy and momentum conservation properties of our code by computing the total energy and momentum for the gyrotron example. Figure 9 exemplarily shows the normalized development in the total energy and total impulse in z-direction for the third order B-splines deposition and interpolation method. In the initial startup phase, the electromagnetic fields introduced by the particles are propagated within the resonator and interact with the particles until the working frequency of the gyrotron is reached for particles and fields. After this initial transient oscillation phase, a stationary level with good energy and momentum conservation is reached.

In conclusion, we have shown that demanding three-dimensional PIC simulations of complex large-scale engineering applications with the proposed high-order coupling techniques are possible and yield good results.

Conclusion and Outlook

In the presented work, two concepts for the coupling of particles with the computational DG grid are discussed and results from PIC simulations are compared. The simulation results demonstrate the ability of the high-order coupling methods, allowing the use of a smaller number of simulation particles and coarser computational grids resulting in a drastical reduction of the simulation time. High-order coupling methods are highly recommended to be used for the simulation of complex engineering applications, such as the gyrotron, to achieve qualitatively good results. The main future research activities will concern the optimization of the high-order discontinuous Galerkin PIC solver with respect to dynamical load balancing and memory utilization issues. Furthermore, new adaptivity concepts like local time stepping techniques are required for the great challenges of 'start-to-end' simulation of advanced electromagnetic wave sources and electrical propulsion systems.

Figure 1 :

 1 Figure 1: Standard PIC cycle

Figure 2 :

 2 Figure 2: p = 2 tetrahedron with 10 IPs.Figure 3: p = 2 hexahedron with 26 IPs.

Figure 3 :

 3 Figure 2: p = 2 tetrahedron with 10 IPs.Figure 3: p = 2 hexahedron with 26 IPs.

Figure 4 :

 4 Figure 4: Examples of different deposition and interpolation methods (simplified 2D view):Left: Deposition via cell mean value method using the DG computational grid directly. Charge current densities of the particles are deposited uniformly within a cell. Center: Deposition from a particle onto the CBM and, vice versa, interpolation from the nodes of the CBM to the particle via first-order B-splines (area or rather volume weighting), where the weights according to the areas (volumes) between the particles and the mesh nodes are used. Right: Coupling of the CBM with the computational DG grid using area (volume) weighting.

  reveal that all simulation experiments reproduce the typical "footprint" of the Weibel instability: The strong rise of the energy of magnetic induction between t ≈ 60 ω -1 pe (a) Magnetic energy (b) Kinetic energy (c) Electric energy (d) Total energy

Figure 5 :

 5 Figure 5: Comparison of energy histories obtained from the simulation of the Weibel instability with different coupling strategies. Reference solution using the CMV coupling method (red full) and the low resolution setup for CMV (blue dashed), first-(blue dashed-dotted) and third-order (black dashed-dotted-dotted lines) B-splines coupling method.

Figure 6 :

 6 Figure 6: Schematic drawing of a cut through the rotational symmetric gyrotron showing the resonator in the center.

Figure 7 :

 7 Figure 7: Geometry of the T E 0,3,1 mode gyrotron resonator, including the electron beam. (Figure by courtesy of S. Illy [33]).

Figure 9 :

 9 Figure 9: Development of the total energy (sum of kinetic and potential energies; red full line) and total momentum in z-direction (black dashed-dotted line) normalized to the initial values.

Acknowledgements

We gratefully acknowledge the Deutsche Forschungsgemeinschaft (DFG) for funding this research within the project "Numerical Modeling and Simulation of Highly Rarefied Plasma Flows". T. Stindl wishes to thank the Landesgraduiertenförderung Baden-Württemberg and the Erich-Becker-Stiftung, Germany, for their financial support. Computational resources have been provided by the Bundes-Höchstleistungsrechenzentrum Stuttgart.


K S C Peerenboom
email: k.s.c.peerenboom@tue.nl

J Van Dijk

W J Goedheer

G Degrez

M Van Der Mullen

J J A M Van Der Mullen

A finite volume model for multi-component diffusion in magnetically confined plasmas

In partially ionized, magnetically confined plasmas, the diffusive fluxes of the different species are coupled. Additionally, the fluxes are directionally coupled as well due to the Lorentz force. The challenge in modelling of multicomponent, magnetized plasmas is to take care of this coupling in the numerical method. In this paper, a complex form of the Stefan-Maxwell equations is used, to account for the coupling between the flow directions. To handle the coupling between the species fluxes in the finite volume method, a generalised, coupled form of the exponential scheme is used. The presented numerical method is applied to a magnetically confined hydrogen jet. The results show that the numerical method is capable of describing typical characteristics of magnetized plasmas, such as anisotropic diffusion and the presence of a pressure gradient sustained by the Lorentz force.

Introduction

Magnetically confined, partially ionized plasmas have several interesting applications, such as the plasma source Magnum-PSI, which is used for material research under ITER relevant conditions [START_REF] De Groot | Magnum-PSI, a plasma generator for plasma-surface interaction research in ITER-like conditions[END_REF]. To model such plasmas, a multi-component approach is needed, where the effect of the magnetic field on the diffusive fluxes of the different species is taken into account. The diffusive flux of a species in a partially ionized, magnetically confined plasma is not simply a linear function of the gradient of the density or mass fraction; the diffusive fluxes of all the species are coupled. In addition, the flow directions across and around the magnetic field lines are coupled via the Lorentz force.

The coupling between the flow directions across and around the field lines makes it more difficult to calculate the diffusion velocity from the Stefan-Maxwell equations. To overcome this difficulty, we will use a complex from of the Stefan-Maxwell equations, based on the approach presented in [START_REF] Giovangigli | Transport algorithms for partially ionized mixtures[END_REF][START_REF] Giovangigli | Multicomponent transport in weakly ionized mixtures[END_REF]. We will extend the approach of [START_REF] Giovangigli | Transport algorithms for partially ionized mixtures[END_REF][START_REF] Giovangigli | Multicomponent transport in weakly ionized mixtures[END_REF] to ambipolar diffusion in a magnetic field. To deal with the coupling between the diffusive fluxes of the species in the finite volume approach, we will use the discretization method of [START_REF] Peerenboom | A coupled discretization method for the continuity equations in multi-component mixtures[END_REF] where the coupling between the species fluxes is preserved.

The numerical method is tested on an illustrative example of Magnum PSI: an expanding hydrogen plasma. It is shown that the numerical method is capable of describing typical characteristics of Magnum PSI (and other magnetized plasmas), such as anisotropic diffusion and the presence of a pressure gradient sustained by the Lorentz force.

This paper is organised as follows. Section 2 explains how the diffusive fluxes in a magnetized plasma are obtained. In Section 3 these fluxes are discretized with the finite volume method in order to solve the species mass balances. Section 4 presents the model description; the results of the model are given in Section 5. The paper ends with a short conclusions section.

Stefan-Maxwell equations and constraints

In the presence of a magnetic field B, the Stefan-Maxwell equations for the diffusive velocity v k of species k can be written as [START_REF] Giovangigli | Transport algorithms for partially ionized mixtures[END_REF]:

l f kl ( v k -v l ) - n k q k p v k × B + y k p j × B = d k , (1)
with f kl the mutual friction coefficient between species k and l, d k the diffusion driving force, n k the number density, q k the charge, p the pressure, y k the mass fraction and j the current density. In the absence of pressure diffusion and thermophoresis, but in the presence of the ambipolar electric field E amb , the diffusion driving force is given by:

d k = -∇x k + n k q k p E amb + n k q k p v × B, (2)
with x k the mole fraction of species k and v the mass averaged velocity of the plasma. The friction coefficients f kl are nonlinear functions of the composition and can be written as:

f kl = x k x l D kl , f kk ≡ 0, (3)
with D kl the usual binary diffusion coefficient [START_REF] Chapman | The Mathematical Theory of Non-uniform Gases[END_REF][START_REF] Ferziger | Mathematical Theory of Transport Processes in Gases[END_REF]. Note that the friction coefficients are symmetric: f kl = f lk , and positive: f kl > 0. Introducing the friction matrix F = (F kl ), with the elements F kl given by:

F kl =    s f ks if k = l, -f kl if k = l, (4)
the Stefan-Maxwell equations can be rewritten as:

l F kl v l - n k q k p v k × B + y k p j × B = d k . (5)
Note that l F kl = 0, which demonstrates that F is singular and hence non-invertible.

Along, across and around the field lines

In the rest of the paper, we will assume that the plasma is cylindrically symmetric and that the magnetic field points in the axial direction. The Stefan-Maxwell equations can then be separated into the directions along, across and around the field lines (Figure 1). In the direction along the field lines ( ), the Stefan-Maxwell equations are not influenced by the magnetic field and can be written as:

l F kl v l = d k . (6)
The flow directions across (⊥) and around ( ) the field lines are coupled:

l F kl v ⊥ l - n k q k p v k B + y k j B p = d ⊥ k , (7)
l

F kl v l + n k q k p v ⊥ k B - y k j ⊥ B p = d k . (8)
The Stefan-Maxwell equations along the field lines can be solved independently from the other directions. The directions across and around the field lines need to be solved together, which can be done by putting together the ⊥ and direction in a 2N × 2N system, with N the number of species. However, from a computational viewpoint, it is favorable to use complex arithmetic [START_REF] Ehrlich | Complex matrix inversion versus real[END_REF].

To facilitate coupled solving of the directions across and around the field lines, the Stefan-Maxwell equations are made complex by adding i times equation [START_REF] Giovangigli | Mass conservation and singular multicomponent diffusion algorithms[END_REF] to equation [START_REF] Ehrlich | Complex matrix inversion versus real[END_REF]:

l F kl v ⊥ l + iv l - n k q k B p v k -iv ⊥ k + y k B p j -ij ⊥ = d ⊥ k + id k . ( 9
)
Substitution of j ⊥ = k n k q k v ⊥ k and j = k n k q k v k gives: l F kl v ⊥ l + iv l - n k q k B p v k -iv ⊥ k + y k B p l n l q l v k -iv ⊥ k = d ⊥ k + id k . ( 10
)
When we introduce the complex diffusion velocity v k = v ⊥ k + iv k , this can be written as:

l F kl v l + i n k q k p v k B -i y k p l n l q l v l B = d ⊥ k + id k . (11)
By adding ‡:

-

n k q k p B l y l v l + y k l n l q l p B l y l v l = 0, (12)
to each equation, the Stefan-Maxwell equations can be written in matrix vector form as:

(F + iF ) v = d ⊥ + id , (13)
where the matrix F is given by [START_REF] Giovangigli | Transport algorithms for partially ionized mixtures[END_REF][START_REF] Giovangigli | Multicomponent transport in weakly ionized mixtures[END_REF]:

F = (I -y ⊗ u) B (I -u ⊗ y) , B = diag n i q i B p , (14)
where u = (1, . . . , 1)

T , y = (y 1 , . . . , y N ) T , and ⊗ denotes the dyadic product of two vectors. The matrix F + iF is singular and has nullspace u.

The singularity in F + iF can be removed by applying the mass flux constraint k y k v k = 0 to the set of Stefan Maxwell equations [START_REF] Giovangigli | Mass conservation and singular multicomponent diffusion algorithms[END_REF]:

(F + iF + αy ⊗ y) v = d ⊥ + id , (15)
where α is a positive parameter. Inversion of (15) now gives:

v = (F + iF + αy ⊗ y) -1 d ⊥ + id = (G + iG ) d ⊥ + id . (16)
Note that instead of direct inversion of the Stefan-Maxwell system, it is possible to use iterative methods as discussed for example in [START_REF] Giovangigli | Transport algorithms for partially ionized mixtures[END_REF]. These iterative methods are computationally cheaper than direct inversion. We will not discuss this method here, since in our case the largest computational expense lies in solving the discretized species mass balances (Section 3), and not in the inversion of the Stefan-Maxwell system on each grid point. The real part of the expression ( 16) gives the diffusion velocity perpendicular to the magnetic field:

v ⊥ = Gd ⊥ -G d . ( 17
)
‡ This expression is added to bring the matrix F to a symmetric form. Since the expression added equals zero, it does not change the diffusion velocities.

The imaginary part gives the diffusion velocity around the field lines:

v = Gd + G d ⊥ . ( 18
)
The driving forces still contain the ambipolar field. In the next section, we will explain how the ambipolar field can be eliminated from the equation.

Zero-current constraint

For the direction across and along the field lines, the ambipolar (zero current) constraint is applied. We assume that there is cylindrical symmetry (with the magnetic field pointing in the z-direction). As a consequence, current can flow around the field lines without causing charge separation. The driving forces along, across and around the field lines can thus be written as:

d = d + z E amb p d ⊥ = d ⊥ + z E ⊥ amb p d = d (19)
The zero-current constraint for the and ⊥ direction can be written as:

z, v = 0, z, v ⊥ = 0, (20)
where z = (n 1 q 1 , . . . , n N q N ) T and denotes the inner product in species space. Substitution of the expression for v and v ⊥ gives:

z, G d + z E amb p = 0, z, G d ⊥ + z E ⊥ amb p -G d = 0, (21)
from which we can determine the expression for the ambipolar field:

E amb p = - z, G d z, Gz , E ⊥ amb p = - z, G d ⊥ -G d z, Gz . (22)
Backsubstitution of the ambipolar field in the expressions of the velocities across and around the field lines gives:

v = G - Gz ⊗ Gz z, Gz d , v ⊥ = G - Gz ⊗ Gz z, Gz d ⊥ -G - Gz ⊗ G z z, Gz d , v = G - G z ⊗ Gz z, Gz d ⊥ + G + G z ⊗ G z z, Gz d . ( 23
)
Note that the matrices relating the diffusion velocities to the driving forces are singular and have nullspace z. This problem was also identified by Giovangigli [START_REF] Giovangigli | Asymptotic stability of equilibrium states for ambipolar plasmas[END_REF] for the case without magnetic field, but in his paper no concrete matrix modifications were proposed. In Appendix A.2, we show appropriate matrix modifications; for a more detailed explanation of these modifications we refer to a future publication [START_REF] Peerenboom | On the ambipolar constraint in multi-component diffusion problems[END_REF]. Appendix A explains how the species mass balances are obtained from (23). In the next section, we discuss how the species mass balances are discretized.

Discretization

As explained in Appendix A, the stationary mass balance for species k in a twodimensional, axisymmetric symmetry can be written as:

l ∂ U kl + ∂ ⊥ U ⊥ kl y l -∂ Γ kl ∂ + ∂ ⊥ Γ ⊥ kl ∂ ⊥ y l = S k , (24)
with U = U kl , U ⊥ = U ⊥ kl , Γ = Γ kl and Γ ⊥ = Γ ⊥ kl the convection and diffusion matrices along and across the field lines as derived in Appendix A. It can be seen that the mass balances for all species form a coupled set of convection-diffusion equations. As shown in [START_REF] Peerenboom | A coupled discretization method for the continuity equations in multi-component mixtures[END_REF], these can not be discretized using the traditional (scalar) finite volume schemes.

To overcome this problem, a coupled form of the exponential scheme has been introduced [START_REF] Ten Thije | The finite volume complete flux scheme for one-dimensional advection-diffusion-reaction systems[END_REF][START_REF] Peerenboom | A coupled discretization method for the continuity equations in multi-component mixtures[END_REF]. In this coupled scheme, the solution vectors on the grid points are coupled via discretization matrices instead of discretization coefficients:

A C y C = A E y E + A W y W + A N y N + A S y S + s C ∆x, (25)
with s C the source on the central point C and ∆x the grid spacing. The discretization coefficients on the central(C), east(E), west(W ), north(N ) and south(S) nodal points are given by:

A E = 1 ∆x ΓB(P e ), A N = 1 ∆x ΓB(P n ) A W = 1 ∆x ΓB(-P w ), A S = 1 ∆x ΓB(-P s ) A C = A E + A W + A N + A S , (26)
with P the Péclet matrix, which is defined as P = ∆xΓ -1 U; the subscripts e, n, w and s denote the interface location where the Péclet matrix is calculated. The Bernoulli function B() is given by:

B(z) = z e z -1 . (27)
Calculation of the of the Bernoulli function of the Péclet matrix is done by diagonalization of the Péclet matrix:

P = VΛV -1 = V    λ 1 . . . λ N    V -1 , (28)
B(P) = VB(Λ)V -1 = V    B(λ 1 ) . . . B(λ N )    V -1 , (29)
with V the matrix with the eigenvectors of P as its columns and Λ a diagonal matrix with the eigenvalues λ i of the Péclet matrix. It can thus be seen that the coupled scheme boils down to: 1) Decoupling the system by changing to the basis of eigenvectors, 2) Discretization of the decoupled system with the scalar exponential scheme, 3) Transforming back to the original basis.

Model description

As an illustrative example for Magnum PSI, we will present a two-dimensional convection-diffusion problem for a hydrogen plasma in a magnetic field. Note that the aim of this example is not to give a complete, self-consistent model of Magnum PSI, but to test the ability of the method to describe its characteristic properties. These properties include magnetic confinement of the charged particles and the presence of a radial pressure gradient. Simulations are performed with the plasma modelling platform Plasimo for different values of the magnetic field. More information on Plasimo can be found in [START_REF] Van Dijk | The plasma modelling toolkit plasimo[END_REF]. In the simulations, we solve mass and momentum balances for the plasma bulk and each modelled species individually. The energy balance is not solved, instead we assume a constant temperature of 10000 K for both the electrons and the heavy particles.

Flow field

The flow field is calculated by solving the mass and momentum balance for the bulk plasma. The momentum balance is given by:

∇ • (ρ v v) = -∇p + ∇ • τ + j × B, (30)
with τ the viscosity tensor and ρ the mass density which is related to the species densities via:

ρ = k n k m k . ( 31
)
The current in the Lorentz force term j × B on the right hand side of the momentum balance needs to be calculated from the diffusion velocities. Since we assumed that the plasma is ambipolar and axisymmetric, there are no radial and axial currents, but there is a current in azimuthal direction:

j = k n k q k v k . ( 32
)
This azimuthal current with an axial magnetic field gives a radial Lorentz force on the flow. Overall mass conservation is given by:

∇ • (ρ v) = 0. ( 33
)
Plasimo uses the SIMPLE algorithm [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF] in combination with Karki corrections [START_REF] Karki | Pressure based calculation procedure for viscous flows at all speeds in arbitrary configurations[END_REF] to solve equations ( 30) and (33) for compressible flows. Velocities are calculated on a co-located grid. Majumdar interpolation [START_REF] Majumdar | Role of underrelaxation in momentum interpolation for calculation of flow with nonstaggered grids[END_REF] is used to avoid oscillations.

Composition

In the composition calculation, it is assumed that the plasma consists of H, H + and electrons. Chemical reactions are neglected. The mass and momentum (Stefan-Maxwell eq.) balances are solved as described in Sections 2 and 3.

Geometry and boundary conditions

The geometry and grid of the model are depicted in figure 2. The geometry is cylindrically symmetric and the magnetic field points in the axial direction. The small cylinder on the left represents the plasma source. The right part of the figure gives the vessel and a gas exhaust. Note that grid stretching is applied in the arc region and the gas exhaust.

At the walls, no-slip boundary conditions are applied to both (radial and axial) velocity components. At the inlet of the arc, and at the gas exhaust, the pressure is prescribed to 15000 Pa and 1000 Pa respectively; homogeneous Neumann conditions are applied for the velocity components. At the symmetry axis homogeneous Neumann boundary conditions are applied to both the velocity and the pressure.

At the symmetry axis and gas exhaust, homogeneous Neumann conditions are applied for all species. At the other boundaries the composition is prescribed. In the arc, the mass fraction of H + is set to 0.2, while in the vessel the H + mass fraction is set to 1e-4. The mass fractions of the electrons and H at the wall are determined by constraints k q k y k /m k = 0 and k y k = 1, respectively.

Results

The model was run with magnetic strengths of 0, 1 and 3 Tesla. The increasing magnetic field strength has the most pronounced effect on the H + density, as shown in figure 3. Comparing the simulations with and without magnetic field, it can clearly be seen that radial diffusion at the exit of the plasma source is strongly reduced in the simulations with magnetic field. The simulation at 3 Tesla shows more pronounced magnetic confinement (narrower beam) than the simulation at 1 Tesla. Qualitatively, this is comparable to the experimental results of [START_REF] Van Rooij | Thomson scattering at Pilot-PSI and Magnum-PSI[END_REF][START_REF] Shumack | Rotation of a strongly magnetized hydrogen plasma column determined from an asymmetric balmer-beta spectral line with two radiating distributions[END_REF]. Additionally, at 3 Tesla, a hollow profile may be observed in the centre of the beam.

The degree of magnetization can be expressed in the Hall parameter, which is defined as the ratio between the Larmor frequency and the collision frequency of the electrons:

H = Ω e ν e = q e B m e ν e , (34)
with Ω e the Larmor frequency and ν e the collision frequency. In figure 4, a radial profile of the Hall parameter, pressure and ambipolar field can be seen. As can be seen, the Lorentz force on the bulk flow due to the azimuthal current gives rise to a radial pressure gradient. The increased pressure in the centre of the beam increases the collision frequency and thereby reduces the Hall parameter, leading to greater radial diffusion of the ions. We believe this is the mechanism which causes the hollow profile in the ion density, and the flattened profile of the pressure and the Hall parameter.

Whereas in an unmagnetized plasma the ambipolar field points outward, since the ions are the less mobile species, it is negative for magnetized plasmas, where the electrons are the less mobile species. This effect is clearly visible in figure 4. Due to the hollow profile in the ion density of the 3 Tesla simulation, there is a small region where the ambipolar field is positive.

In figure 5, the axial plots of the pressure and the Hall parameter can be seen. When the plasma is magnetized, the pressure on the axis is larger in the expansion region, due to the earlier mentioned radial Lorentz force. At the inlet of the plasma source where the pressure is high the Hall parameter is low. As the plasma expands the Hall parameter increases, because of the decreasing pressure.

Conclusions

We discussed a method for the modelling of multi-component, magnetized plasmas with a finite volume scheme. In this method, both the coupling between the flow directions and the species is conserved. The coupling between the flow directions is taken into account by using the Stefan Maxwell equations in complex form. In the discretization, the coupling between the species was taken into account by using a coupled form of the exponential scheme. As a test case, we used the method to model an expanding hydrogen jet under increasing magnetic field strengths. The results showed expected characteristics of magnetized jets, such as magnetic confinement of the ions, radial pressure gradients and reversal of the ambipolar field.

Although in this paper it was assumed that the plasma is cylindrically symmetric and the magnetic field points in the axial direction, it should, in principle, be possible to apply the method to more general magnetic configurations. Fluxes could be calculated in terms of the directions along, across and around the field lines; after that the fluxes could be transformed to the directions of the cell faces. However, in this transformation, special attention should be paid to the consistency of fluxes at the interface of two adjacent cells, which is of utmost importance in the finite volume method.

⊥BFigure 1 .

 1 Figure 1. The directions along ( ), across (⊥) and around ( ) the magnetic field lines.

Figure 2 .

 2 Figure 2. Geometry and grid of the model. 1) inlet of the plasma source, 2) wall, 3) gas exhaust, 4) wall/end plate, 5) symmetry axis.

Figure 3 .

 3 Figure3. H + densities for different values of the magnetic field. The stronger the magnetic field, the narrower (more confined) the plasma beam becomes. In the 3 Tesla simulation, a slightly hollow profile can be observed.

Figure 4 .Figure 5 .

 45 Figure 4. Radial plots at the axial position z=0.1.In the pressure, the ambipolar field as well as in the Hall parameter, the increasing confinement at increasing magnetic fields can be observed.

Appendix A. The species mass balances

The stationary mass balance for species k can be written as:

with ρ the mass density, y k the mass fraction of species k, v the mass averaged velocity of the plasma, v k the diffusion velocity and S k the mass source term. Rewriting (A.1) in components along and across the field lines gives:

In this appendix, we will rewrite the total mass flux to:

where U , U ⊥ , Γ and Γ ⊥ are the convection and diffusion matrix parallel and perpendicular to the field, respectively.

Appendix A.1. Conversion to mass fractions

To be able to solve the species mass balances for the mass fractions y k , the diffusion driving forces need to be expressed in the massfractions and the gradients of the mass fractions. Specification of the driving forces gives:

The gradients of the mole fractions can be converted to the gradients of the mass fractions via a matrix:

Special care has been taken to make M non-singular [START_REF] Giovangigli | Mass conservation and singular multicomponent diffusion algorithms[END_REF]. The vector z can be converted to y via:

with N = diag (ρq k /m k ). Note that the driving forces proportional to the gradient of the mass fractions will end up in the diffusion matrix. The forces proportional to the mass fractions will end up in the convection matrix together with the mass averaged velocity.

Appendix A.2. Regularization

In order to keep the diffusion matrix Γ non-singular, the matrix G -Gz⊗Gz z,Gz should be non-singular. It can, however, be seen that z is in the nullspace of this matrix. To remove this singularity, the matrices need to be regularized. This can be done by the following modification:

with w = M -T q and q = (q 1 /m 1 , . . . , q N /m N ) T . The coefficient β is a positive parameter.

Appendix A.3. Conversion to mass fluxes

The diffusive velocity vectors v and v ⊥ can be converted to the diffusive mass fluxes J and J ⊥ by multiplication with the species mass densities:

with R = diag (ρy k ).


Sophie ; Roux
email: sophie.roux@inrets.fr

Jimmy ; Armoogum
email: jimmy.armoogum@inrets.fr

Jean-Loup Madre
email: jean-loup.madre@inrets.fr

Roux Sophie

Armoogum Jimmy

Madre Jean-Loup

DYNAMIC OF CAR OWNERSHIP AND CAR USE IN FRANCE SINCE THE 1960S

Keywords: motorization, car ownership, car use, driving license, behaviour, dynamic, surveys

Transportation has experienced extraordinary growth during the last century and still evolves throughout the twenty-first century. To clarify the changes to come, studies of the history of past habits and behaviors are essential. The knowledge of mobility in France is partly based on National Travel Surveys (NTS). In the past these surveys have been conducted five times (1966-67, 1973-74, 1981-82, 1993-94 and 2007-08), which gives five photographies of mobility for households living in France. The definitions and principles of these surveys 12th WCTR, July 11-15, 2010 -Lisbon, Portugal Dynamic of car ownership and car use in France since the 1960s ROUX, Sophie; ARMOOGUM, Jimmy; MADRE, Jean-Loup haven't been modified since the beginning which makes the measurement of structural changes easier. The goal of this paper is to analyze the behavioral changes related to the motorization since the 1960s. The analysis of changes will be focused on the evolution of household motorization (non-motorized, with one car and multi motorized) and the number of driving license holders since the 1960s.

INTRODUCTION

Transportation has experienced extraordinary growth during the last century and still evolves throughout the twenty-first century. To clarify the changes to come, study of the history of past habits and behavior is essential. The knowledge of mobility in France is partly based on National Travel Surveys. In the past these surveys have been conducted five times (1966-67, 1973-74, 1981-82, 1993-94 and 2007-08), which gives five photographies of mobility for households living in France. NTS are rich in information, but they are generally run every ten to twelve years, it is difficult to measure the trends in behaviors and understand the story of a subtle history of mobility. The protocol of NTS has changed which requires the development of tools to analyze mobility in comparative static. However the definitions and principles of these NTS haven't evolved which facilitates the measurement of structural changes. The goal of this paper is to analyze behavioral changes related to the motorization since the 1960s. After a short description of the French National Travel Survey, we will analyze the evolution of the households' car fleet, and then a phenomenon directly linked to the evolution of the fleet, the numbers of individuals who hold a driving license. In the last 60 years we have known a major change as at the beginning of the period we have a rapid growth in term of household motorization and in the last 20 years car fleet is rising by multi-equipment. Thus we will see these dynamic throughout the analysis of these five cross-sections surveys 1 , it will provide an historical perspective for France.

THE FRENCH NATIONAL TRAVEL SURVEYS

Once per decade, the Ministry of Transport and the National Institute of Statistics use to conducting a National Travel Survey with the scientific support of INRETS. Each survey is fitted into a continuation of the previous to allow one to measure the structural changes. The aim of these surveys is the description of short and long distance trips made by households living in France, as well as their access to and use of public and private transport means. In France, The NTS is the only survey completed at the nation level describing all trips made by households living in France, whatever their purpose, length, mode, season or time of day. These surveys also emphasize the knowledge of the vehicle fleet, use and of (Armoogum & all, 2007).

The definitions and principles of these NTS haven't evolved [START_REF] Armoogum | Mise en perspective des enquêtes nationales transport 1981-82, 1993-93 et 2007-08[END_REF]. However, as one might envisage, some field of surveys and issues have changed over the 5 surveys requiring us to take a common field to all surveys. Thus, the results presented below will only state that cars without taking into account the light-utility vehicles. Similarly, the variables of geographic zoning in France have changed significantly since the first survey. The results presented below are in a constant zoning, the reference being the last zoning dating from the 1999 census. The definition of the zones area is as follow: city center, suburb and rural area. These zones area can be use only for the 3 last surveys.

EVOLUTION OF THE HOUSEHOLDS' MOTORISATION: MULTI-MOTORIZATION AND FLEET OF VEHICLES OLDER

Cars have penetrated our lifestyle in the last century. The household motorization is marked by two major evolutions (Fig. 1). The first is linked with the acquisition or not of a car by the household. In 1966, slightly more than half of households had a car. Until now, this proportion increases, more than 80% of households had a car in 2007. The motorization had a strong growth between 1966 and 1981. The share of non-motorized households decreases of less than 18 points between these 2 dates, that is an average annual decline of 3.0%. The second evolution is linked with the number of cars available in each household. Since mid 60's we note a continuous growth of the share of the households multi-motorized. In 1966, less than 6% of households had at least two cars this figure is about one third in 2007. This last trend pushes the average number of cars per household. What are the factors that may affect the household's motorization°? Henry [START_REF] Fraisse | L'équipement des ménages en voiture. Comparaison des résultats de l'enquête nationale transport et communication 1993-94 et du panel SOFRES[END_REF] showed that household car equipment is explained by the influence of economic characteristics (income), demographic (age, sex, etc.) and spatial (habitat type, housing environment etc.). For this initial exploratory work on the history of the motorization, we will only take an interest in household size and the zone of residence. Does the number of cars increases proportionally with the size of the household? Is it higher when households live in rural areas?

The average number of cars per household grows strongly with the size of household: it is 0.6 for one-person households and nearly 1.7 for households with three persons or more in 2007 (Fig. 2). Indeed, the motorization of a household is all more profitable when it concerns the largest number of people. Over the period 1981 and 2007, the average number of cars per household increases parallel for households of two persons and three persons or more (+0.6 point). However, the single-person households have the highest growth rates between these 2 dates: 90% on the overall of the period against 63% for two persons households and 31% for larger households size. We can think that, for multi-persons households of, a car is the family transport mode. However, this strong shift for one-person households can be explained partly by the use of the car that is more and more individualized and on the other hand, certainly by a cohort effect. Indeed, in 1981, one person households aged 65 and over had 0.13 vehicles in average against 0.46 in 2007. The elderly households keep their past habits and still use their cars today [START_REF] Bodier | La voiture : une habitude qui se prend jeune[END_REF]. When examining the share of the number of cars depending on the zone of residence (city center, suburb and rural area), we find the same trends as those described above: decrease of the share of non-motorised households and an increase of multi-cars households in the 3 zones. However, these two phenomena evolve in different ways depending on area of residence [START_REF] Bussière | Vers la saturation ? Une approche démographique de l'équipement des ménages en automobile dans trois régions urbaines[END_REF]. The share of households having no car is highest in the cities centers and is the lowest for multi-motorization regardless of the period (Fig. 3). We can think that the public transport more developed in city center competes with the automobile. Conversely, where access to public transport is more difficult, particularly in Dynamic of car ownership and car use in France since the 1960s ROUX, Sophie; ARMOOGUM, Jimmy; MADRE, Jean-Loup rural areas, the share of households with at least one car is more important, as well as multimotorization (Fig. 4). In 2007, car ownership is almost generalized in rural areas. The motorization of households living in suburbs is between these two extremes, where the transport network is more developed than in rural areas but less than in cities centers (Fig. 5). Our previous work on the capital region (Madre, Armoogum 1996), show that the motorisation increases with the distance from the city center. We found that the level of motorisation for the city center of Paris was the lowest in France but this figures reach the highest level for France in the outer suburbs of Paris.

CONDITION OF ACQUISITION AND AGE OF CARS: AN AGEING FLEET

In 1966, about 50% of the cars were acquired new by households, they are less than 28% in 2007 (Fig. 6). Second-hand cars take up an increasingly important place in the equipment of Dynamic of car ownership and car use in France since the 1960s ROUX, Sophie; ARMOOGUM, Jimmy; MADRE, Jean-Loup French households. Between 1973 and 2007, the share of cars purchased second-hand has increased by 34%. We can thus conclude that the vehicle fleet of households tends to get older over time (Bodier & all, 1995) [START_REF] Hivert | Observation en continu du comportement automobile des ménages et de la structure du parc », INRETS, p 4. Dynamic of car ownership and car use in France since the 1960s ROUX[END_REF]. The sharing out of age of cars and the age of cars at the time of purchase tends to increase the average age of cars (Fig. 7 and8). Actually, during the year 1966, the average age of cars was 5.4 years against 8.8 years in 2007 and the average age of the cars at the time of purchase has increased from 2.7 years to 4.2 years between 1981 and 2007. (Fig. 9). Technological advances have greatly contributed to the increase in the average age of cars (cars are more robust). Households now keep their vehicles longer than before (Jeger, 2001). But one can also wonder whether this increase is not the combined result of the multi-motorization of households and a different use of the second or third vehicle, less recent. Indeed, older car had lower annual mileage (Fig. 10). This phenomenon is identical over time in all surveys. We can think that households would tend to use older vehicles to make short trips, thereby lengthening the time of life of the car, and to use the main car more frequently. The household would tend to share their mileage between several vehicles (Join-Lambert, 1981). The evolution of household car ownership is a consequence of the democratization of driver's licence. Indeed, the more people hold a driving licence, the more the population are likely to drive and the more they are likely to increase the car fleet by buying cars. The study of driver's license holding therefore goes hand in hand with motorization.

DRIVING LICENCE HOLDERS: MORE AND MORE, MORE AND MORE WOMEN, AND MORE AND MORE SENIORS

Since 1966, the rate of driving licenses holders is constantly increasing. It is nearly 42% in 1966 and has doubled in 2007, nearly 84% among people aged 18 and over (Fig. 11). A cohort effect can be suspected. According to the biographical section of the 2007 NTS, the analysis of obtaining the driver's license by cohort sheds a twofold light. The first concerns the proportion of respondents having a driver's license: more the cohort is young and more the proportion of persons in possession of the driver's license is high. The second evolution is that of the age when obtaining the driver's license. It decreases with the youth of the cohort (Papon & all, 2009). When comparing the proportion of men and women who have driver's license, we observe that for different period men have more often driver's license than women (Fig. 12). Men are more than 70% in 1973 to hold it and more than 91% in 2007. But the major factor for the difference in proportion of driver's licences between both genders: it is decreasing over time from nearly 37 points in 1973 to 15 points in 2007. Women were 33% to have their license in 1973 and over 76% in 2007. Women in 2007 have slightly exceeded the level reached by men in 1973. The driver, male "privilege" in the 60s [START_REF] Fagnani | Activités féminines et transports urbains[END_REF], has faded over the surveys. Before 1997, young people could get their driving license for free when in France we used to have a military service (which was mandatory for young male). With the end of this military service, there was a decrease in the number of driving licenses issued to young male especially those living in rural areas (Avrillier & all, 2008). We can probably also attribute the increased number of women who have a driver's license during Dynamic of car ownership and car use in France since the 1960s ROUX, Sophie; ARMOOGUM, Jimmy; MADRE, Jean-Loup the 2nd half of the 20th century to the movement in pursuit of emancipation and empowerment of women, particularly professional. Arguably, the commute pushes the massive car use, and, in the recent decades this correlation is very evident among women. This car use is probably due to the poor access to public transport. At the same age, the proportion of people with a driving license is higher and higher over time (with one exception, 18-24 years living in suburbs and rural areas) and for people over 55 years, the differences between surveys are highest (Fig. 13, 14 and 15). People aged over 65 have experienced the development of the automobile during their middle age, and thus have practiced driving late or not (Gallez & all, 1997) at all especially women. The democratization of the driver's license, its non-renewable phenomenon in France, and a lower average age of obtaining a driver's license led to the fact that people who had the driver's license in the early 80s have always a driver's license today. Overall, the proportion of people with driving licenses is higher for people living in a rural area, and then follows those living in suburbs and finally those living in the city center. We see logically the same Dynamic of car ownership and car use in France since the 1960s ROUX, Sophie; ARMOOGUM, Jimmy; MADRE, Jean-Loup phenomenon already described for the motorization of households. More public transport is developed, the less people have the driver's license.

CONCLUSION

Since the '60s, the motorization of households has rapidly increased resulting by a decrease in non-motorization of households accompanying by an increase in the multi-motorization and an aging of the car fleet in household. We note a decrease in the driver's license age and the increasing share of women who hold it is link to the evolution of the motorization.

Figure 2 :

 2 Figure 2: Average number of cars per household by household size at different period

  Figure 3: Share of number of cars ownership in the city center

Figure 5 :

 5 Figure 5: Share of number of cars ownership in rural areas

Figure 9 :Figure 10 :

 910 Figure 9: Average age of car and average age of car at the time of purchase at different period

Figure 11 :Figure 12 :

 1112 Figure 11: Proportion of individuals aged 18 and over having a driver's license at different period

Figure 13 :Figure 14 :Figure 15 :

 131415 Figure 13: Proportion of individuals to hold a driving license depending on age in the city center at different period

12th WCTR, July 11-15, 2010 -Lisbon, Portugal

Dynamic of car ownership and car use in France since the 1960s

ROUX, Sophie; ARMOOGUM, Jimmy; MADRE, Jean-Loup Indeed, the multi-motorisation could allow the second partner (who was most often a woman [START_REF] Fagnani | Activités féminines et transports urbains[END_REF]) to drive a car, especially to reach their workplace. The growth of the motorization and it's diffusion to the women allows the urban sprawls. Moreover, the possession of a car is now closely linked to the type of territory where we live. However, these surveys do not allow on the one hand to follow individual changes of behavior as data are collected once per decade on a new sample, and on the other hand to describe the history of mobility continuously. The method that is traditionally used to catch life stories is the biographic survey. In France, the 2007-08 issue of the NTS innovates as compared with the previous ones, among other additions, by the collection of biographic retrospective data giving way to the analysis of individuals' past behavior over their lifetime. By providing longitudinal data on the interviewees' lives, these data will allow richer analyses on mobility behavior and habits, and present a historical perspective on travel means in France in the past.