1 Introduction and experimental overview

The thermal protection system of spacecrafts is one of the most expensive parts for space missions. Large efforts are undertaken to develop and use special materials for the thermal protection system (Blosser Reference Blosser1996). One idea was to use magnetic fields for heat flux reduction. If it is possible to reduce the heat flux by externally applied electromagnetic fields, one would be able to use simpler and cheaper materials. In 2002 the European Space Agency (ESA) started an investigation of heat flux mitigation by means of MHD (magnetohydrodynamics) in partially ionised argon (Gülhan et al. Reference Gülhan, Esser, Koch, Siebe, Riehmer, Giordano and Konigorski2009). Argon was chosen as test gas for this campaign due to its chemical inertness. This simplifies the analysis and makes the inclusion of complicated plasma chemistry unnecessary.

In the experiment an arc jet was used to create a flow field, with a typical Mach number between 5 and 10, which hits a target. The heat flux on the target surface was measured with and without a magnetic field. Applying a dipole-like external magnetic field a strong heat flux mitigation was observed (Gülhan et al. Reference Gülhan, Esser, Koch, Siebe, Riehmer, Giordano and Konigorski2009). However, even up to now, the experiments could not be reproduced by models based on Navier–Stokes fluid codes for neutrals with fluid models for the plasma components. One was not able to explain the main effects leading to the heat flux reduction.

One major problem for the modelling is the lack of information about the plasma component in the experiment. There exist only average values for density ( $n_{e}\approx 10^{11}\,\text{cm}^{-3}$ ) and temperature ( $T_{e}\approx 0.7~\text{eV}$ ), whereas the neutral component is well characterised using different diagnostic methods in terms of temperature and velocities, including axial and radial profiles at specific locations. The average neutral density $n_{n}$ is ${\approx}10^{15}\,\text{cm}^{-3}$ .

In addition, spectroscopic data and a video from the experiment are available (Detlev Konigorski, personal communications, January 2014). During the experiment a clear change in the emission was visible in photos looking into the vacuum chamber onto the blunt body. The flow illumination changed to a red colour when the magnetic field was applied, confirmed in observations by Kranc, Cambel & Yuen (Reference Kranc, Cambel and Yuen1969) in a similar experiment. The large bright area between the arc jet exit and the target shrank by activating the magnetic field and a small bright spot directly in front of the centre of the target appeared (see figure 6 of Gülhan et al. (Reference Gülhan, Esser, Koch, Siebe, Riehmer, Giordano and Konigorski2009)).

There are several groups that have tried to reproduce the experiment in simulations. They all applied Navier–Stokes codes or direct simulation Monte Carlo codes for the neutral gas description. In contrast to the complex and rather complete models of neutral transport, only fluid models for the description of the plasma were used, mostly with estimated transport coefficients (Katsurayama et al. Reference Katsurayama, Kawamura, Matsuda and Abe2008; Otsu et al. Reference Otsu, Katsurayama, Konigorski and Abe2012).

Using a self-consistent kinetic particle method instead of fluid models a more realistic description of neutral and plasma dynamics and their interaction is possible. Specifically, an existing particle-in-cell (PIC) code with Monte-Carlo collisions (Tskhakaya et al. Reference Tskhakaya, Matyash, Schneider and Taccogna2007) is used for this study.

Perturbed plasma source distribution functions are typical for plasma sources, e.g. those shown by Godyak & Piejak (Reference Godyak and Piejak1990) for RF (radio frequency) discharges or by Behringer & Fantz (Reference Behringer and Fantz1994) for DC (direct current) discharges. Such non-Maxwellian characteristics of the plasma components are also expected from the arc jet.

2 Simulation method

The well-established PIC method is combined with Monte Carlo collisions (PIC-MCC). In plasma physics, the PIC method is a widely accepted method to obtain a better understanding of the basic physical mechanisms of various systems because the PIC-MCC approach provides full insight into all (microscopic and macroscopic) parameters. A more detailed description of the PIC-MCC method itself can be found in several reviews (Birdsall & Langdon Reference Birdsall and Langdon2004; Tskhakaya et al. Reference Tskhakaya, Matyash, Schneider and Taccogna2007)

The principle of PIC is to follow so-called super-particles. Each of these represents the same number of physical particles following the same trajectory as single particles due to their identical mass-to-charge ratio, resulting in the same Lorentz force. The electric field for the Lorentz force is calculated self-consistently on a spatially equidistant grid from the Poisson equation. In the electrostatic approximation, magnetic field effects are included in the Lorentz force, but are determined just by their external sources (e.g. magnets). Corrections from plasma currents are neglected. Particles experience collisions according to Monte Carlo collision algorithms. This allows one to implement all relevant collision types. For this work electron–neutral elastic, ionisation, excitation, Coulomb collision, neutral–neutral and charge exchange collisions are implemented. Effective ionisation and radiation energy losses are included in the model. Neutrals are treated as plasma particles with zero charge. Technically, this is identical to a DSMC (direct simulation Monte Carlo) model for the neutrals.

In the experiment, the neutral density is approximately four orders of magnitude larger than the electron density (Gülhan et al. Reference Gülhan, Esser, Koch, Siebe, Riehmer, Giordano and Konigorski2009). If one uses the same super-particle factor for plasma and for neutrals, one runs into computational memory limits due to the large number of neutral particles.

Due to the relatively low electron energy compared with ionisation energies above $15~\text{eV}$ , there are only a very small number of ionisation collisions. Due to the very low ionisation degree, the neutrals can be assumed to be practically decoupled in density from the plasma component.

Therefore, different weights for neutrals and plasma particles are taken care of by scaling the cross-section of collisions with neutrals linearly. The cross-section for neutral–neutral collisions is increased quadratically by this factor. These scalings are technically equivalent to modifying the super-particle factor for neutrals. This allows one to reduce the number of pseudo-particles for neutrals to a moderate number. In this work a scaling factor of $10^{3}$ is used.

The spatial domain is two-dimensional (axial and radial). Due to collisions three components in velocity space need to be resolved. The domain consists of 400 axial and 140 radial cells. With a cell size of $\unicode[STIX]{x0394}r\approx 1.05~\text{mm}$ this leads to a domain of $42~\text{cm}\times 14~\text{cm}$ . The cell size and the time step $\unicode[STIX]{x0394}t\approx 10^{-12}~\text{s}$ are chosen in order to resolve the Debye lengths and the electron cyclotron frequency. For this system the cyclotron frequency ( $\unicode[STIX]{x1D708}\approx 7\times 10^{10}~\text{Hz}$ ) in the region with strong magnetic fields is higher than the plasma frequency ( $\unicode[STIX]{x1D714}_{p,e}\approx 3\times 10^{10}~\text{Hz}$ ).

The target is placed directly at the symmetry axis at the right side of the domain. It has a radius of $3.6~\text{cm}$ and its surface is located at $26~\text{cm}$ according to the experimental conditions. The target has a dielectric surface made out of quartz ( $\unicode[STIX]{x1D716}=4.27$ ) and the inner part is metallic. All boundaries except the left one have a fixed potential of zero. The left domain boundary is set to $0.5~\text{V}$ , to represent the exit of the arc jet and ensure that particles move towards the target. At outer domain boundaries all particles are absorbed. Plasma particles that hit the target are also absorbed, only neutrals are reflected at the surface of the target. This represents the inertness against chemical reactions of the noble gas argon.

The size of the domain is scaled down by a factor of 100 to overcome spatial resolution limits. One applies a self-similarity scaling keeping the relevant non-dimensional parameters constant, namely the ratio of system length to gyroradii and to mean free path. The first represents the influence of the magnetisation by the external magnetic field, the second the collisional effects. By increasing all densities and the magnetic field by the same factor, the ratio of system length to gyroradii and to mean free path is conserved. This ensures the physical correctness of the simulation as proven by Taccogna et al. (Reference Taccogna, Longo, Capitelli and Schneider2005).

The external magnetic field was calculated using the finite-element magnetic method solver (FEMM) (Meeker Reference Meeker2010). The topology of the dipole-like magnetic field is shown in figure 1.

Figure 1. Magnetic field in the simulation domain.

For neutrals a uniform background density (representing background pressure effects) is combined with the main particle source at the left side of the domain near the axis with a radial extension up to $1~\text{cm}$ . The location of the source is shown in green in figure 1. Neutrals are injected with a Maxwellian-distributed velocity with a drift of $2070~\text{m}~\text{s}^{-1}$ in axial direction according to the mean velocity measured by microwave interferometry. For every sampled velocity it is ensured that particles leave the source with an initial angle smaller than $12\,^{\circ }$ to reproduce the angular distribution of the particle source in the experiment. The characterisation of the plasma species coming from the arc jet is not known experimentally. Therefore, variation of the source parameters was used to represent the experiment as well as possible. Finally, the following plasma source distributions were used. Ions have similar source distributions as the neutrals.

Electrons are injected with a two-Maxwellian distribution. The first one has a shift in energy of approximately $1~\text{eV}$ to represent the measured averaged temperature from experiment and a width of $1~\text{eV}$ . A second Maxwellian with the same shift and variable width and source strength is added. The parameters are fitted to get as close as possible to the experimental results without a magnetic field. The best result are obtained for a source strength of $1/3$ of the bulk source and a width of $20~\text{eV}$ . The second contribution is necessary to provide some electrons with higher energies of above $13~\text{eV}$ , otherwise there is no line emission in the visible range as observed in the photos and in the spectroscopic analysis of the experiment. According to Kranc et al. (Reference Kranc, Cambel and Yuen1969) the red emission is due to excitation of metastable argon atoms with energies of approximately $11~\text{eV}$ . The resulting electron energy distribution function is shown in figure 2 for different axial positions.

The number of injected particles is chosen to approximate the average density for electrons ( $n_{e}=n_{i}\approx 2.8\times 10^{11}~\text{cm}^{-3}$ ) and for neutrals ( $n_{n}\approx 1.6\times 10^{15}~\text{cm}^{-3}$ ). For a typical run approximately $4\times 10^{6}$ super-particles for each plasma species and $6\times 10^{6}$ super-particles for neutrals are used.

Figure 2. Electron axial velocity distribution functions for different axial positions averaged over $10^{5}$ time steps ((a): without magnetic field, (b) with magnetic field).

During the experiment strong fluctuations were observed when the magnetic field was applied. To reflect this observation in the model in the case of an applied magnetic field, radial anomalous diffusion is implemented. In addition to the Bohm diffusion (Bohm et al. Reference Bohm, Burhop, Massey and Williams1949), which scales like

(2.1) $$\begin{eqnarray}D\propto \frac{T}{B},\end{eqnarray}$$

with the temperature $T$ and the magnetic field $B$ , a constant transport is taken into account. This contribution is motivated by the fact that the arc jet itself delivers a turbulent flow due to angular gas insertion. Movies of the experiment give the impression that this turbulence is strongly increased within the source by the activation of the magnetic field. For cylindrical geometry the radial diffusion has to be implemented correctly as derived by Bormann, Brosens & De Schutter (Reference Bormann, Brosens and De Schutter2001):

(2.2) $$\begin{eqnarray}p_{\pm }^{r}=\frac{1}{2}\left(1\pm \frac{1}{2k}\right),\end{eqnarray}$$

with $k=r/\sqrt{2D\times \unicode[STIX]{x0394}t}$ for a diffusion coefficient $D$ . The diffusion coefficient is set to $10~\text{m}^{2}~\text{s}^{-1}$ for plasma particles and $1~\text{m}^{2}~\text{s}^{-1}$ for neutrals. The estimation for the coefficients is based on diffusion coefficients from magnetic fusion (Perkins et al. Reference Perkins, Barabaschi, Boucher, Cordey, Costley, Deboo, Diamond, Fujisawa, Greenfield and Hogan1994). Bohm et al. (Reference Bohm, Burhop, Massey and Williams1949) showed that turbulence in a plasma can be created or at least be enhanced by magnetic fields. The non-local effect is also observed in systems like the cylindrical Hall thruster from PPPL (Princeton Plasma Physics Laboratory) (Ellison, Raitses & Fisch Reference Ellison, Raitses and Fisch2012).