A Three ‐ Dimensional Finite Volume Arbitrary Lagrangian ‐ Eulerian code for plasma simulations

A Three-Dimensional Finite Volume Arbitrary Lagrangian-Eulerian simulation code was developed to study different plasma physics problems in 3D+t. The code is based on a complex multi-component species program with transport and radiation terms written and applied to plasma and fusion physics problems. Three different examples are shown: double-base chemical propellant combustion, ignition and propagation of a thermonuclear detonation wave, and, the development of the Kelvin-Helmholtz (KH) instability in local plane slab models of the magnetopause, showing the response of a background equilibrium to the excitation by finite amplitude perturbations generated upstream.


INTRODUCTION
A Three-Dimensional Finite Volume Arbitrary Lagrangian-Eulerian simulation code was developed, that includes ion viscosity, thermal conduction (electrons and ions), magnetic diffusion, thermonuclear production or chemical reactions (including local and non-local capabilities for the deposition of the charged thermonuclear products), Bremsstrahlung radiation, EOS (from the ideal gas to the degenerate electron gas).
Nowadays there exist several efficient numerical codes in the subject.They are the result of many years of research, test and development as they are able to deal with very complex physics and geometry and most important, they are supported with empirical validation.However they are not always accessible mainly because of the high costs of their licenses.On the other hand the construction of own computational codes brings the following important advantages: a) to get a deeper knowledge of the physical processes involved and the numerical methods used to simulate them, and b) more flexibility to adapt the code to particular situations in a more efficient way than a closed general code would.
These advantages have motivated the present work, which is intended to set a starting point in the subject.The present work will show in brief the distinct aspects that have been worked out in order to cover the simulation objective.

PLASMA EQUATIONS
We use a 3D, time dependent, two fluids (ions and electrons), two-temperature model.The equations are (MKS units except where noted) Continuity 0 v U w wU t (1) Electric field transformation where U is the mass density, n e (n i ) the electron (ion) number density, Y k the mass concentration of species, Z k the overall reaction rate for the species, q D k flux vector diffusion for species, v the fluid (ion) velocity, p the total pressure (ions plus electrons), H e (H i ) the specific electron (ion) internal energy, T e (T i ) the electron (ion) temperature (in eV), q e (q i ) the electron (ion) heat flux, 3 the viscosity tensor (including artificial viscosity for numerical purposes), E (E') the electric field in the laboratory (plasma) frame, B the magnetic field, j the current density, R is the momentum transfer between ions and electrons, W eq the equipartition time, e is the elementary charge, m e the electron mass, m p the proton mass, and k the Boltzman's constant.
The above plasma equations are completed with the equation of state.

NUMERICAL RESOLUTION
We have used a Lagrangian Finite Volume method with hexaedric cells that move at arbitrary velocity.The notation and definition of the geometry is described Ref. 1.
The main advantage of this method is that the geometry is imposed to the cells and not to the operators.This means, for example, that a general problem with cylindrical geometry is calculated in Cartesian coordinates.
The numerical scheme is based on the Finite Volumes method.The integration in time is sequential.This means that each process is integrated with a different uncoupled method during a time step.
Diffusion processes (of species concentration and energy) are integrated with ADI methods (Alternating Direction Implicit), full implicit or explicit methods according to the case.
Hydrodynamics is integrated with a new variant of the ICE method (Implicit Continuous Eulerian) treated in Ref. 2. For the numerical integration we assume that a variable is constant over the volume or surface where it is defined.
Source terms from chemical or nuclear reactions are explicitly integrated in time.
Each integration is performed in a Lagrangian way, thus temporally avoiding the convective terms.Then, after each calculation cycle, mesh vertices are moved arbitrary over the fluid.This is done in order to dynamically adapt the mesh to the solution.As a consequence, convective terms are restored.This reorganization technique allows the calculation to automatically increase spatial resolution where steep gradients of any flux variable are produced.The numerical scheme is implemented according to the following general procedure: 1. Read constants and calculation parameters.

FINITE VOLUME FORMULATION
This method is very similar to finite difference and consists of integrating conserved quantities over a cell.
For a general fluid property I, the following identity where integration is performed over a cell of volume V with boundary S moving at an arbitrary speed w.
For example, setting v U I , and using the momentum equation The integration domain is represented with a structured irregular mesh, with fixed connectivity made of hexaedric cells.Each cell is surrounded by 6 faces and 8 vertices.Fluid variables are assigned to staggered locations in the mesh.Pressure, internal energy, density, species concentration, cell volume and mass are all assigned to cell centers.Coordinates (x,y,z) and velocities (u x ,u y ,u z ) are assigned to cell vertices.
Mesh is constructed with GID pre and post processing software, then recalculated with the Near Orthogonal Grid (NOG) method.

TIME INTEGRATION
Step 1. Integration of hydrodynamic terms.By "hydrodynamic terms" we refer only to the Euler equations included into the full set of conservation equations.So in this part of the calculation fluid variables are updated in time assuming that the other processes (diffusion, chemical or nuclear reaction, convection) are temporally frozen during one time step.
The calculations necessary to advance a solution one step in time 't are separated in two distinct phases.Phase One consists of an explicit Lagrangian calculation, except mesh vertices are not moved.Phase Two consists of an iteration that adjusts the pressure gradient forces to the advanced time level.This phase, which is optional, eliminates the usual Courant-like numerical stability condition that limits sound waves to travel no further than one cell per time step.In an explicit method pressure forces can be transmitted only one cell each time step, that is, cells exert pressure forces only on neighboring cells.When the time step is chosen so large that sound waves should travel more than one cell, the one cell limitation is inaccurate and a catastrophic instability develops.To overcome this instability time-advanced pressure gradients may be used.Unfortunately the time advanced pressures depend on the accelerations and velocities computed from those pressures, so an iterative calculation of the set of equations is required.Physically, iteration offers a means by which pressure signals can traverse across more than one cell in a unit time step.
For regions where a low Mach number is expected, a good rule of thumb obtained from discussions kept in www.CFD-online.com is that the relative error in the pressure field, must be less than the square of the Mach number.So for M | 10 -3 , the pressure should be iterated to one part in a million.If we include a safety factor of 0.1, just to ensure convergence, it would be necessary to run the calculation with 8 bytes in a floating point format number (double precision), since 4 bytes offer only about 6 significant figures.
Step 2. In the second sequential step we consider the time advance of the internal energy and species concentration due to thermal conduction through Fourier's law, species diffusion through Fick's law and viscosity through Newton's law.In this step hydrodynamics, convection and chemical or nuclear reactions are kept frozen.
In what follows we give the basis of the numerical implementation only for the thermal conduction, however, the other diffusive processes are treated almost exactly in the same way.
The total internal energy (mH into cell i will change due to the contribution of diffusion flux (q) through cell boundaries according to: All time derivatives are forward in time, transport coefficients are evaluated at cell faces.Because T and Y are assigned to cell centers, a weighted average between adjacent cells is performed in order to obtain values at the face.
In order to express the diffusive flux, the temperature gradient must be considered at cell boundaries.On the other hand, it is performed a local linearization between temperature and internal energy (with fixed concentration in the thermal conduction case).
All the above mentioned assumptions allow us to implement explicit, or implicit or a regular ADI (Alternating Implicit Direction) 3,4 methods.During each ADI stage, several tridiagonal systems must be solved, one for each row of cells lying in a column.This is done with the Thomas algorithm.
Step 3. In the final sequential step chemical or nuclear reactions are considered.They represent the source terms in the energy and species concentration equations.In the present work they are explicitly integrated in time.This is tentatively done because it is possible that the explicit treatment would not be the principal cause of the time step restriction.
Note that due to the time splitting procedure the overall truncation error is O(¨r,¨z,¨t).

MESH REZONING
The following considerations led to implement an adaptive technique to correct the nodal concentration dynamically.
Lagrangian cell methods are not for describing flows undergoing large distortions, because cells may undergo severe deformations.
The physical solution may present steep gradient of temperature, or species concentration, or heat release in different regions.Moreover, the localization of these regions is unknown at the beginning of the calculation, and they may be originated or not depending on chemical or nuclear reactions, shock fronts or any other condition.
The adaptive method consists of shifting mesh vertices over the fluid in order to keep a reasonable mesh structure and increase the spatial resolution where the physical solution demands.
Each vertex is shifted to a new location therefore exchanging mass, species, energy and momentum among the surrounding cells.The following procedure is applied.
Step 1. Calculate the new coordinates and then calculate the volume exchange between cells Step 2. Associated with the volume exchanged there will also be a mass, species concentration and total energy exchange.The mass, concentration or energy per unit mass assigned to this volume can be computed in various ways.The use of a simple average of the quantities on either side of the lines leads to a computational instability, but a stable calculation can be obtained by weighting the average in favor of the value in the cell from which the quantity is subtracted.This is the upstream or donor cell convective flux approximation.
Step 3. A shift in a vertex is accompanied by a momentum exchange between neighborhood vertices.A similar donor cell procedure is used for momentum exchange among vertices.
Step 4. Once updated the conserved cell and vertex properties, the following quantities must be retrieved: a) cell volumes from the new vertex coordinates, b) new cell densities, c) new species concentration, d) vertex masses from cell masses by averaging among adjacent cells, and, e) vertex velocities from new momentum.

EXAMPLES Double-base propellant combustion wave
The double-base propellant is a solid energetic material commonly used for propulsion systems and gas generators.Because its relatively large stability under room conditions, an external heat source (igniter) must be provided to initiate the propellant burning.During the ignition process, heat is transported by conduction inside the solid.As consequence of the corresponding temperature increase, rate reactions of the degradation reactions begins to accelerate.This endothermic process of decomposition generates reactive species.
As degradation evolves, the solid concentration drops down from one to zero, then the igniter is no longer required, and combustion could reach selfsustained conditions.As a consequence of convection and mass conservation through the burning surface, a chemically reacting flow emerges from the interface toward the gas phase.
The emerging molecules react each other in the gas phase, developing great amounts of heat in the flame zone.The heat release profile is not evenly distributed because of the difference between characteristic times of the reactions involved.This is the principal cause of the flame structure and the steep gradients in the temperature profile.
Fig. 1 shows the flame structure of a double-base propellant.It can be qualitatively observed that the luminous flame approaches the burning surface and burning rate increases as pressure increases. 5n order to increase spatial resolution where steep temperature gradients and high chemical heat release are produced the rezone velocity U with which vertices are moved over the fluid is calculated using a distribution function for vertex concentration according to 6 1 where c is a normalization factor such that Then the velocity U for each vertex becomes: where W, is a time relaxation factor established to keep vertex shifts controlled without diminishing the time step.This calculation of U is 1D, where the coordinate x represents the path length following one computational direction.
The flame fronts are efficiently captured by the mesh as soon as they appear or they move through the integration domain.Fig. 2 shows the steady state temperature profiles for different bulk pressures in the gas phase near the burning surface.This behavior is in agreement with the experimental results of Ref. 5. The hydrodynamics is best tested in a 2D-plane channel, where the propellant lies at the inlet section.The calculation is able to deal with transonic flows as well as subsonic ones, as can be seen in Fig. 3.

Thermonuclear detonation wave
Fast implosions of annular current sheaths (fast Z pinches) have made important progress in reaching high energy densities.New possibilities for a broad spectrum of experiments, from X-ray generation to controlled thermonuclear fusion are devised.At the same time experiments have shown that a cylindrical liner can successfully confine a magnetized, hot plasma column.The technical parameters necessary for ICF applications being inaccessible so far, thus, numerical simulations provide an insight into the structure of both the liner and the plasma, giving the range of parameters required for ignition and detonation in such a liner-pinch system 7 .Spark ignition appears to be a more promising approach to ICF than volume ignition.The idea of spark ignition by a hot spot in a Z-pinch is not new, coming out of a long and varied experience with Z-pinch necking (plasma focus, fiber pinch, gas puff, etc.).It seems that the logical way of using ignition in DT is that of amplifying the fusion energy in a conical channel to a level at which a detonation can be triggered in advanced fuels 8 .Fig. 4 shows a typical result, corresponding to physically feasible situations.Even the mesh suffers from large distortion, a fully Lagrangian calculation was used.

Development of the KH instability
As an example of a pure Eulerian calculation, a 3D numerical simulation of the large amplitude perturbations and waves at the duskside low latitude boundary layer of the magnetopause generated by an interplanetary tangential discontinuity is presented in this Conference 9 .In Fig. 5 temperature and magnetic lines plot show that the KH waves appear already after a few theoretical e-folding time, due to the great oscillations generated by the tangential discontinuity stress that has reduced the amplification time.See Ref. 9 for more details.

CONCLUSIONS
The code was successfully applied to different problems.The sequential implementation resulted of practical use, because each physical process was treated with a different numerical method.
The time step was restricted by efficiency considerations rather than stability conditions.When time step becomes larger, the calculation progress faster because fluid variables, transport properties, etc., are updated with less frequency.But the implicit stage (phase two) of the hydrodynamic advance requires more iterations to perform.An optimum time step for which computational effort is minimum has been used.

2 areFIGURE 1 .
FIGURE 1. Flame structures of a double-base propellant.Then, new coordinates (x*,y*) for each vertex are calculated in such a way that the Lagrangian coordinate K, defined as

FIGURE 2 .
FIGURE 2. Calculated temperature profiles in the gas phase under different bulk pressures in steady conditions.

FIGURE 3 .
FIGURE 3. Density plot of the Mach number in steady state, showing the capability of the code to deal subsonic and supersonic regimes in the same integration domain.

FIGURE 4 .
FIGURE 4. Time evolution of the temperature plot of a DT thermonuclear detonation wave propagating in a cylindrical channel.
, and repeat from step 5 until the end of the calculation is reached.Steps 1 to 4 are performed once at the beginning of the calculation.Steps 5 to 16 are performed once at each time step.Boundary conditions are implemented with phantom cells and vertices.
7. Update Y and H due to chemical reactions.8. Update Y, and H due to diffusion process (conductivity, diffusion, viscosity) 9. Update p, H, coordinates and velocities due to