Mixed Flow Ionospheric Aerodynamics

The background material laid out above motivates an investigation of the nature of the assumptions employed in Capon et al ([1]). This work therefore aims to address the research question
How does the presence of multiple species of plasma affect the ability to predict ionospheric aerodynamic forces?
This offers a number of avenues of study, including analysis of known satellites orbits, scaling parameter analysis, grain based plasma physics, and many others. For simplicity, the specific goals of this analysis are confined to:

1. 1.

   Extend set of parameters in Capon et al [22] to account for multi-species interactions in the formulation of a control surface.

2. 2.

   Simulation of multi-species flow using particle in cell code developed at University of New South Wales (UNSW) Canberra, and subsequent drag analysis.

In this way, the extensive parametric study of single species flows performed in [20] can be directly compared to mixed plasma flows, thereby allowing mixed species analysis to be an extension, rather than replacement, of the existing response surface model.

In studying mixed species plasma flows, the difference in interactions of mixed species compared to single species flows can be largely attributed to either differing charge-mass ratios[28], differing gas-surface interactions [10], or collisional chemical interactions[29]. In the plasmas being modelled (mesothermal LEO flows), the mean free path is on the order of kilometres [30], and therefore the analysis of chemical reactions is neglected in this work. Numerous gas-surface interaction models have been employed in previous works, and found to contribute significant variation between different simulated drag coefficients [10]. While this provides an exciting avenue for future research, this work holds a constant model of diffuse, neutralised reflection (elaborated on in Chapter 2), in order to allow direct comparison of results to those in [1]. The primary point of analysis will therefore be the mixing of ions with varying charge:mass ratios. For simplicity, analysis will be limited to mixes of only two species, chosen to be singly ionised H+ and O+. The choice of these species has 3 distinct motivators:

• •

  O+ and H+ represent the two extremes of charge to mass ratios in singly ionised particles in LEO.

• •

  O+ and H+ are the two most populous plasma species in the LEO regions under consideration (see Fig.1.2).

• •

  Use of O+ and H+ allow direct comparisons to the flows analysed in Capon et al [27], as these are all single species H+ or O+ flows.

Using these species, the parameters analysed were the ratio of oxygen to hydrogen ions ($\rho_{H+}:\rho_{O+}$), and the sensitivity of these results to variations in satellite charge. Specifically, analysis began with a representative O+ single species flow, and varied the ratio $\rho_{H+}:\rho_{O+}$ in steps, providing a series of simulated flows with a range of H+ and O+ concentrations. The flows were varied by either number density or mass density, investigating the effects of both parameters on drag distributions. These two series of simulations were then repeated using a lower body potential $\phi_{0}$, testing sensitivity of results to variations in body charge.

The DSMC method is a technique for simulating high Knudsen number flows, where the mean free path length of particles is sufficiently high for particle motion and collisions to be decoupled over sufficiently small time step $\Delta t$. The fundamental mechanic of DSMC is the superparticle method, which simulates the motion of particles directly, by grouping a given number of particles into a ’macroparticle’, and then tracking the motion and collisions of these macroparticles. As the ratio of real particles to simulated particles (a parameter referred to as $N_{Equiv}$ in this work) approaches 1 (each macroparticle representing one real particle), the simulation approaches a direct calculation of the position of each of the particles present, providing a computationally intense, but mathematically intuitive, approach to simulations. The general DSMC algorithm used in the dsmcFoam code implemented in OpenFOAM is as follows

1. 1.

   Initialisation - The mesh is generated using blockmesh utility from parameters given in blockmeshDict file, and a number of macroparticles is assigned to each cell in the mesh. This number is given by $\frac{V_{c}\cdot\rho_{k}}{N_{Equiv}}$ for cell volume $V_{c}$ and species density $\rho_{k}$ . These macroparticles are assigned representative temperatures and velocities from the initial flow conditions supplied, modified by the Maxwell-Boltzmann distribution, producing a collection of thermalised particles with the desired average macroscopic quantities. Each Cell is also assigned the list of macroparticles it contains.

2. 2.

   Particle Push - Iterating across all macroparticles, the macroparticles positions are updated ballistically: $\bf{x}_{k}\rightarrow\bf{x}_{k}+\bf{v}_{k}\Delta t$

3. 3.

   Update Cell Occupancy - Each Cells list of occupying particles is updated to reflect the changed positions in step 2.

4. 4.

   Collision partner selection - Within each cell, the average number of expected collisions is determined by iterating across the particles present, and (Using the No Time Counter (NTC) numerical method laid out in [29]), the average number of collisions expected for the present particles over a timestep $\delta t$ is determined, and an appropriate number of pairs of particles is selected to collide.

5. 5.

   Collisions - the particles in a collision pair have their internal and kinetic energy redistributed according to kinetic theory described in [29], with new velocity determined by relative kinetic energies, and new directions of travel randomised (while conserving total collision momentum).

6. 6.

   Reactions - If particle reactions are being used, these are resolved according to rules dictated by the type of reaction. In the cases studied here, the only reaction was a neutralisation reaction upon contact with charged surface, so the particles were reflected diffusely, and their charge was changed from +1 to 0.

7. 7.

   Update Cell occupancy - Each cell’s list of occupying particles is once again updated, this time to include changes in particle species from reaction interactions. The algorithm then returns to step 2.

Most of this algorithm is very intuitive, with the exception of the collisions method (known as No Time counter (NTC)). In this method, the collisions occuring within a cell are not necessarily based on the actual relative trajectories of particles, but rather the relative scattering cross section $\sigma$ of two particles. In this way, macroscopic collision properties are preserved, without the necessity of individually calculating the relative trajectories of every particle with every other, a computationally imposing task.

To simulate a satellite in mesothermal flow, the approximation of small satellite geometry as a cylinder was used. The principle advantages of a cylinder are its relative simplicity, and the ability to directly compare with orbital motion limited theory (Outlined in Chapter 5) The simulation was performed on the grid shown in Fig.2.3, representing the top half of the region being simulated (symmetry boundary conditions were employed to lower computational load), with the $obstacle$ patch marked, and the respective boundary conditions (BC) as described here.

As the simulation is symmetrical along the $xz$ axis (through $(0,0,0)$), and along the $yz$ axis, the bottom patch was rendered with symmetry boundary conditions, as were the front and back patches. the simulation was performed with only one cell in the z direction, but particles velocities and positions retained their respective z component, allowing them to pass into and out of the symmetry regions on the $front$ and $back$ faces.
When a particles motion passes through a symmetry patch, such as the $bottom$ patch, it is returned to the simulation with the component of its velocity in the direction of the outward face reversed (as in Fig.2.4) as if a particle with symmetry along that face was passing in from the other side. The patch marked $flow$ is treated with the flow boundary condition, and $obstacle$ face with corresponding obstacle condition. At each timestep, particles that pass out of faces with flow BC are removed, while a number of new macroparticles given by

$$N_{k,new}=A_{\perp}\cdot v_{flow}dt\left(\frac{\rho_{N,k}}{N_{Equiv}}\right)$$

(2.18)

is introduced, where $A_{\perp}$ is the cross sectional area of the face to the flow direction, $\rho_{N,k}$ and $v_{flow}$ are the number density and magnitude of flow velocity for species $k$, and all other parameters are as defined in Appendix B. When particles interact with faces using the obstacle BC, they are reflected diffusely, and charged particles are converted to their respective neutral particle. the momentum change in particles reflecting from the $obstacle$ patch is recorded and summed for each face, to be used in determining direct drag force on the object.

The $flow$ face is held at a constant potential $\phi$ determined by initial conditions (0V in this case), while symmetry faces boundary conditions have potential gradient set to zero, and the $obstacle$ face has a fixed potential (denoted as $\phi_{0}$, usually -50V in this thesis) The mesh density was determined alongside other key initial parameters with a numerical convergence study laid out in Chapter 3.

In order to identify steady state solution to these simulations, it is necessary to collect a large range of data as a function of time. Large datasets cause it to be grossly inneficient for this data to be recorded for anything other than the final, steady state conditions. To fulfill the requirements of monitoring the simulation at each timestep, it is therefore necessary to implement secondary methods, to record representative values for each timestep, separate to the writeFiles stage (In which the final state of the simulation is recorded for all fields). In the typical simulation, the total simulated time is 0.001s, while the algorithm iterates in timesteps $\delta t$ of $10^{-7}$s, so it impractical to record much more than a handful of representative values. Being the primary interest of this investigation, those values were chosen to be the total direct force on the satellite, and the charged (indirect) force. The direct (neutral) force on the satellite can be obtained as a sum of force moments from diffusely reflected particles, a functionality which is built in to the pdFoam solver already, providing net direct force values for each timestep.Taking a discrete approach to charged force calculation would be impractical

To determine the charged drag on the cylinder, we begin by considering the general momentum balance of a mesothermal plasma-body interaction laid out by [36]. This will be demonstrated with the inclusion of a magnetic field, although this will later be neglected under electrostatic assumptions for these simulations. In a mesothermal flow, the thermal ion and electron pressures are assumed to be negligible compared to the streaming energy, letting us write the ion and electron momentum equations as

$$n_{i}m_{i}\left[\frac{\partial\mathbf{v}_{i}}{\partial t}+(\mathbf{v}\cdot\nabla)\mathbf{v}_{i}\right]=q_{i}n_{i}(\mathbf{E}+\mathbf{v}_{i}\times\mathbf{B})$$

(2.19)

$$n_{e}m_{e}\left[\frac{\partial\mathbf{v}_{e}}{\partial t}\right]=q_{e}n_{e}(\mathbf{E}+\mathbf{v}_{e}\times\mathbf{B})$$

(2.20)

Where $\mathbf{E}$ and $\mathbf{B}$ are the electric and magnetic field vectors, $n_{i}$ and $n_{e}$ the ion and electron number density, and $\mathbf{v}$ the flow velocity [37]. Taking the sum of these two equations, and introducing charge density as $\rho_{c}=(\sum_{k}^{K}Z_{k}q_{e}n_{k})-q_{e}n_{e}$ and current density $\mathbf{J}=\rho_{c}\mathbf{v}$ gives

$$(m_{i}n_{i}+m_{e}n_{e})\frac{\partial\mathbf{v}}{\partial t}+m_{i}n_{i}(\mathbf{v}\cdot\nabla)\mathbf{v}=\rho_{c}\mathbf{E}+\mathbf{J}\times\mathbf{B}$$

(2.21)

Next we introduce the divergence form of Gauss’ law and the Maxwell-Ampere equation to express $\rho_{c}$ and $\mathbf{J}$ in terms of $\mathbf{E},\mathbf{B}$

$$\rho_{c}=\epsilon_{0}\nabla\cdot\mathbf{E}$$

(2.22)

$$\mathbf{J}=\frac{1}{\mu_{0}}(\nabla\times\mathbf{B})-\epsilon_{0}\frac{\partial\mathbf{E}}{\partial t}$$

(2.23)

Equations 2.22 and 2.23 can be substituted into Eqn.2.21 to give

$$(m_{i}n_{i}+m_{e}n_{e})\frac{\partial\mathbf{v}}{\partial t}+m_{i}n_{i}(\mathbf{v}\cdot\nabla)\mathbf{v}=\epsilon_{0}(\nabla\cdot\mathbf{E})\mathbf{E}+\frac{1}{\mu_{0}}(\nabla\times\mathbf{B})\times\mathbf{B}-\epsilon_{0}\frac{\partial\mathbf{E}}{\partial t}\times\mathbf{B}$$

(2.24)

using the product rule and Faradays law, the time derivative of the electric field can be written in terms of $\partial/\partial t(\mathbf{E}\times\mathbf{B})$[1]

$$\frac{\partial}{\partial t}(\mathbf{E}\times\mathbf{B})=\frac{\partial\mathbf{E}}{\partial t}\times\mathbf{B}+\mathbf{E}\times\frac{\partial\mathbf{B}}{\partial t}=\frac{\partial\mathbf{E}}{\partial t}\times\mathbf{B}-\mathbf{E}\times(\nabla\cdot\mathbf{E})$$

(2.25)

substituting Eqn.2.25 into Eqn.2.24 gives

$$\begin{split}(m_{i}n_{i}+m_{e}n_{e})\frac{\partial\mathbf{v}}{\partial t}+m_{i}n_{i}(\mathbf{v}\cdot\nabla)\mathbf{v}&=\epsilon_{0}[(\nabla\cdot\mathbf{E})\mathbf{E}-\mathbf{E}\times(\nabla\times\mathbf{E})]\\ +\frac{1}{\mu_{0}}[(\nabla\cdot\mathbf{B})\mathbf{B}-&\mathbf{B}\times(\nabla\times\mathbf{B})]-\epsilon_{0}\frac{\partial}{\partial t}(\mathbf{E}\times\mathbf{B})\end{split}$$

(2.26)

The curls in Eqn.2.26 can be eliminated by using

$$\mathbf{A}\times(\nabla\times\mathbf{A})=\frac{1}{2}\nabla(\mathbf{A}\cdot\mathbf{A})-(\mathbf{A}\cdot\nabla)\mathbf{A}$$

(2.27)

Eqn.2.26 then becomes

$$\begin{split}(m_{i}n_{i}+m_{e}n_{e})\frac{\partial\mathbf{v}}{\partial t}+m_{i}n_{i}(\mathbf{v}\cdot\nabla)\mathbf{v}&=\epsilon_{0}[(\nabla\cdot\mathbf{E})\mathbf{E}+(\mathbf{E}\cdot\nabla)\mathbf{E})]\\ +\frac{1}{\mu_{0}}[(\nabla\cdot\mathbf{B})\mathbf{B}+&(\mathbf{B}\cdot\nabla)\mathbf{B}]-\epsilon_{0}\frac{\partial}{\partial t}(\mathbf{E}\times\mathbf{B})\end{split}$$

(2.28)

This can be simplified considerably with the use of the Maxwell stress tensor $\mathbf{\overline{T}}$

$$\mathbf{\overline{T}_{i,j}}=\epsilon_{0}\left(E_{i}E_{j}-\frac{1}{2}\delta_{i,j}E^{2}\right)+\frac{1}{\mu_{0}}\left(B_{i}B_{j}-\frac{1}{2}\delta_{i,j}B^{2}\right)$$

(2.29)

where $\delta_{i,j}$ is a kronecker delta. taking the divergence of $\mathbf{\overline{T}}$ and introducing the Poynting vector $\mathbf{S}=\mu_{0}^{-1}(\mathbf{E}\times\mathbf{B})$, Eqn.2.28 can now be written as

$$m_{i}n_{i}(\mathbf{v}\cdot\nabla)\mathbf{v}-\nabla\mathbf{\overline{T}}=-(m_{i}n_{i}+m_{e}n_{e})\frac{\partial\mathbf{v}}{\partial t}-\epsilon_{0}\mu_{0}\frac{\partial\mathbf{S}}{\partial t}$$

(2.30)

For a steady state solution, we allow the time derivatives terms on the RHS to tend to zero. To obtain a control surface formulation of the resulting equation, we then consider surface $S$ containing volume $V$, defined by outward normal vector $\mathbf{\hat{n}}$,, such that equation 2.30 becomes

$$\int_{S}n_{i}m_{i}(\mathbf{v}\cdot\mathbf{\hat{n}})\mathbf{v}dS-\int_{S}\mathbf{\overline{T}}\cdot\mathbf{\hat{n}}dS=0$$

(2.31)

Where the two integrals can be taken to represent ion momentum, and the electromagnetic stress (from left to right) respectively. This equation implies that the variation in mechanical momentum between flows into and out of the volume are caused by net electromagnetic (and mechanical) interactions within the volume (conservation of momentum in an electromagnetic field must account for the momentum stored by the fields). The influence of the plasma on a body can then be found by defining a second surface $S_{2}$ to represent the charged body, such that 2.30 becomes

$$\int_{S}(n_{i}m_{i}(\mathbf{v}\cdot\mathbf{\hat{n}})\mathbf{v}-\mathbf{\overline{T}}\cdot\mathbf{\hat{n}})dS+\int_{S_{2}}(n_{i}m_{i}(\mathbf{v}\cdot\mathbf{\hat{n}})\mathbf{v}-\mathbf{\overline{T}}\cdot\mathbf{\hat{n}})dS_{2}=0$$

(2.32)

Where the second integral is the force exerted by a plasma on the body defined by $S_{2}$. This gives the force on a surface S as

$$\mathbf{F_{c}}=-\int_{S}n_{i}m_{i}(\mathbf{v}\cdot\mathbf{\hat{n}})\mathbf{v}dS+\int_{S}\mathbf{\overline{T}}\cdot\mathbf{\hat{n}}dS$$

(2.33)

The first term in this is analogous to the discrete summation method employed in calculating direct force, replacing a summation of individual momentum exchanges with an integral over the average mechanical momentum exchange. The second term provides a convenient way to calculate indirect force for any face on the body surface, requiring only the electric field and surface normal at that point.

In order to implement the indirect force calculations, the integral in Eqn.2.33 is converted to a summation, for the program to perform at each step of the simulation.

$$\textbf{f}_{c}=\sum_{i}^{S}\textbf{T}_{i}\cdot\mathbf{n}_{i}$$

(2.34)

Where the index i represents a face in the surface S, and $\mathbf{n}_{i}$ is the normal vector for this face. In order to calculate each entry in this integral, the electric field is extracted (in our case we used the electric field calculated in the potential assignment step of the algorithm in Fig.2.2), over the surface of interest (the $obstacle$ patch in this case). The value for $\hat{n}$ can be obtained with a loop over the surface of interest, taking each face element, and for a face defined by points (a, b, c, d), taking the normal area vector as the cross product $N=\vec{ba}\times\vec{bc}$ shown in Fig.2.5. This can be done with a loop over the face set for the surface of interest and the ’SArea’ method in openFoam.

Two approaches were considered in the method of calculating the product $\textbf{T}_{i}\cdot\hat{n}_{i}$, with the first being the piecewise definition, in which the $x,y,z$ elements of $\hat{n}$ and E are extracted, and the Indirect force for each face is calculated as the vector with entries $(F_{E,x},F_{E,y},F_{E,z})$ (with subscript $C$ used instead of $i$ to clearly distinguish from index i in Eqn.2.34):

	$\displaystyle F_{C,x}=(E_{x}^{2}-\frac{\textbf{E}^{2}}{2})n_{x}+E_{x}E_{y}n_{y}+E_{x}E_{z}n_{z}$		(2.35)
	$\displaystyle F_{C,y}=E_{x}E_{y}n_{x}+(E_{y}^{2}-\frac{\textbf{E}^{2}}{2})n_{y}+E_{y}E_{z}n_{z}$
	$\displaystyle F_{C,z}=E_{x}E_{z}n_{z}+E_{y}E_{z}n_{z}+(E_{z}^{2}-\frac{\textbf{E}^{2}}{2})n_{z}$

This method suffers from inefficiencies associated with converting the vector fields into 3 scalar fields each, then converting them back after performing the additions. The second approach to this calculation solves these efficiency issues by keeping E and $\hat{n}$ as vector fields, and performing the calculation through openFoam vector manipulation methods instead. This gives the expression:

$$\textbf{F}_{C}=(\textbf{E}\textbf{E}^{T}-\frac{1}{2}{\bf{\mathbf{1}}}|\textbf{E}|^{2})\cdot\mathbf{n}$$

(2.36)

This vector based solution structure allows for both more compact code, and more efficient and scalable solutions, making it the ideal choice to meet the constraints of the indirect force monitoring code.
This method was implemented as an optional method in pdFoam, allowing it to be removed in cases where the additional data was not worth the increased computational cost (Although the computational cost of this method is $<1\%$ of the cost of the charge assignment method alone in the meshes used here, so this is not significant enough to ever justify disabling this method). A sample of output data is shown in Fig.2.6, in which the average of charged and direct force in the x direction for each timestep is plotted over the raw values. Notably, the charged force has considerably less noise than that of the direct force, as the direct force is determined by particles impacting the cylinder in one time interval, while the charged force is determined by the local potential, which shows much less time variation than a measurement of particle density at a point. The sharp rise at the start of the plot is a feature of the formation of the sheath, with the quasi-neutral plasma (which provides no charged force, and allows little room for ions to accelerate, leaving a reduced direct force) in the initial conditions separating to form a plasma sheath (steady state condition).

A numerical model is convergent if a sequence of model solutions with increasingly refined conditions approaches a fixed value [38]. Likewise, a numerical model is consistent if and only if the sequence converges to the solution of the physical equations governing this particular phenomena. In the case of PIC simulations, the key approximations made are the discretization of the timestep, the mesh, and the macroparticle approximation, splitting continous space and time into discrete timesteps and cells for the algorithm to iterate through, and replacing a large number of discrete particles with a smaller number of macroparticles (all of which are outlined in Chapter 2.)
This means that if this model is both convergent and consistent, the solution given by our simulations should approach the real value of the systems modelled as the cell size and timestep are lowered, and the number of macroparticles used is raised. As a finer mesh and smaller timestep is used, the computational cost increases accordingly, and it is the goal of the Numerical convergence study undertaken here to determine a set of parameters for an optimal balance of computational cost to simulation fidelity to be used in further simulations.

For this study, the mesh density was scaled by a parameter SX, where an SX of 1 corresponds to the mesh shown in Fig.2.3, and the number of cells in each of the $x$ and $y$ directions scales roughly linearly with an increasing SX, giving a relationship between total number of cells and SX parameter shown in Fig.3.1. The time interval for each step of the algorithm was denoted $\Delta t$, while the ratio of real particles to simulated particles was denoted $N_{Equiv}$.

The main bounds imposed on these parameters were those imposed by computational cost at the fine end, and the Courant Condition at the coarse end. The Courant condition is a necessary condition for the convergence of the solutions of the PDEs laid out in Chapter 4. It states that for a given flow velocity $u$, spatial dimension $\Delta x$ and time interval $\Delta t$,

$$C=\frac{u\Delta t}{\Delta x}\leq C_{max}$$

(3.1)

Where $C_{max}=1$ is usually used for a time stepping solver. Physically, this means that any propagating wave with velocity $u$ will spend at least one time step in each consecutive cell on its path, allowing the behaviour of these signals to be resolved.

A range of cases were prepared, with conditions laid out in table below, varying both the SX and $N_{Equiv}$ parameters, in the range 4 - 56 and $10^{3}$ - $10^{6}$ respectively.

Constant	Value	Unit
$v_{\infty}$	$7500$	$ms^{-1}$
$\rho_{H+}$	$2\times 10^{10}$	$m^{-3}$
$\rho_{O+}$	$2\times 10^{10}$	$m^{-3}$
$\phi_{B}$	$-50$	$V$
$r_{0}$	0.3	$m$

These cases were run to a steady state, and the charged and direct force in the x direction were measured using the methods in section 2. In the case of mesh parameters and Macroparticle convergence, the average force values varied significantly enough between similar cases that the standard deviation ($\sigma$) of force values for the last 50$\%$ of the simulation was used as a surrogate value to quantify convergence. $\sigma$ was defined as in Eqn. 3.2, where $f_{t}$ is the integrated force at timestep $t$, N is the total number of timesteps sampled, and $\mu$ is the mean . As the number of timesteps taken after each case reaches steady state is constant between the cases, the $\sigma$ values for that case become an effective measure of quantifying variation within the force values.

$$\sigma_{f}=\sqrt{\frac{\sum_{t}(f_{t}-\mu)^{2}}{N}}$$

(3.2)

As the set of parameters grows arbitrarily fine, the value of $\sigma$ approaches a converged value $\sigma_{0}$. The criteria chosen to qualify a set of parameters as adequately converged was the value of $\sigma$ corresponding to those parameters being within 10$\%$ of $\sigma_{0}$. Fig.3.3 shows variation of $\sigma_{0}$ with SX, for a number of different $N_{Equiv}$ values.

Analysis of the $\sigma_{F}$ values for the charged and direct force showed a strong convergence with changing $N_{Equiv}$, (Fig.3.3) (given by use the coarsest set of parameters), while an increasing SX value only accounted for a convergence of 10$\%$ or less of the maximum value (Fig.3.2), with more than half of this variation occuring in the range SX$<16$. As the normalised $\sigma$ values taken above this mesh parameter seem to not be significantly affected by changing SX, the mesh was considered to be converged at an SX of 30 (and a value of SX=40 was chosen to be used in final simulations, allowing better resolution of flow features in the fore-wake).

The standard deviation values showed a strong convergence as $N_{Equiv}$ decreased, in both the direct and indirect forces (Fig.3.3). To identify the point at which $\sigma$ had converged to within the $10\%$ margin selected, the curves for $\sigma_{D}$ and $\sigma_{C}$ were normalised (setting the maximum value of each curve to 1), and fit to an empirical function for ease of interpolation. The function chosen was given by $\sigma=A/x^{n}+\sigma_{0}$, for arbitrary scaling parameter A, final converged standard deviation $\sigma_{0}$, and variable of interest $x$.

The data and fit were presented in terms of total number of macroparticles, rather than $N_{Equiv}$, as the simulations of different densities aim to keep the same number of macroparticles per cell, by adjusting the $N_{Equiv}$ parameter accordingly.

In Fig.3.4, the final value $\sigma_{0}$ towards which $\sigma$ converges was found to be $10^{-17}N\approx 0$ in the direct case, and $0.004N$ in the indirect case, implying the standard deviation approaches zero as the number of macroparticles rises to infinity (Realistically, this is not the case, as the number of macroparticles used is bounded by the upper limit where each macroparticle represents one real particle). This fit allows the identification of the $10\%$ margin as being $\approx 4.5\times 10^{5}$ macroparticles.

A selection of cases was prepared with these new values for SX and $N_{Equiv}$, with the time step for one iteration of the algorithm varying from $5\times 10^{-6}$s to $5\times 10^{-9}$s between cases. Each case was run until a time of 0.001s as before, upon which convergence was complete When measuring convergence of the time step $\Delta t$, the standard deviation becomes an inneffective tool, as the number of timesteps being sampled increases as $\Delta t$ decreases, resulting in the standard deviation growing with a finer timestep, rather than converging. As a result, the direct and charged force averages must be used directly, resulting in Fig.3.5. It is immediately obvious that the dominant source of variation for timesteps of $10^{-6}$s or smaller is random variation within the simulation, while the timestep of $5\times 10^{-6}s$ shows significant variation from these, in both charged and direct force values. This is to be expected, as the average cell size is around 1mm, so the courant number for the main flow of this simulation would be $C\approx 10$, well above the convergence criteria $C\leq 1$. While the force values appeared to converge for $\Delta t=10^{-6}s$, the time step used in further simulations was chosen to be $10^{-7}$s, allowing the courant condition to still be satisfied for particle speeds of up to 70$\%$ greater magnitude than the bulk flow, such as those in the close fore-wake region.

This derivation functions on the assumptions of electrostatic interactions of a collisionless, unmagnetised multi-species plasma with multiply charged ions with a body at a fixed potential with respect to a quasi-neutral freestream plasma. Electron distribution is modelled as a Boltzmann electron fluid [24], an isothermal, inertia-less electron fluid. Quantities referencing the satellite body are denoted with the subscript $0$, while quantities referring to the main flow (without perturbations) are denoted with $\infty$.

As in section 2, Equations 2.2 and 2.23 describe the evolution of a system in terms phase space distribution function $f$, potential $\phi$, and space-charge density $\rho$ given by

$$\rho=\left(\sum_{k}^{K}q_{k}n_{k}\right)-q_{e}n_{e}\quad n_{k}=\int f_{k}d\mathbf{c}$$

(4.1)

where the electron density $n_{e}$ in the Boltzmann electron fluid model is described by

$$n_{e}=n_{e,\infty}exp\left[\frac{q_{e}\phi(x)}{k_{B}T_{e}}\right]$$

(4.2)

Where $n_{e,\infty}$ is given as the sum of the respective ionic charges in a quasi-neutral freestream (as net charge density in freestream=0).

$$n_{e,\infty}=\frac{1}{q_{e}}\sum^{K}_{k}q_{k}n_{k,\infty}$$

(4.3)

substituting the expressions for charge densities from Eqns 4.1-4.3 into equation 2.23 gives a pair of expressions describing distribution of a $K$ species plasma denoted by subscript $k$

$$\frac{\partial f_{k}}{\partial t}+\mathbf{c}_{k}\cdot\nabla_{x}f_{k}-\frac{q_{k}}{m_{k}}\nabla_{x}\phi\cdot\nabla f_{k}=0$$

(4.4)

$$\nabla^{2}\phi=-\frac{q_{e}}{\epsilon_{0}}\sum^{K}_{k}\frac{q_{k}}{q_{e}}\left(n_{k}-n_{k,\infty}exp\left[\frac{q_{e}\phi}{k_{B}T_{e}}\right]\right)$$

(4.5)

To express Eqns 4.4-4.5 in dimensionless form, for each parameter $\psi$ with units $U_{n}$, a substitution is made such that $\Psi=\frac{\psi}{\psi_{0}}$, where $\psi_{0}$ is an characteristic value for the corresponding parameter, also using units $U_{n}$, and $\Psi$ is a dimensionless variable corresponding to that unit. As an example, the velocity $\mathbf{c}_{k}$ becomes $\mathbf{c}_{k}=v_{k,0}\mathbf{u}_{k}$, where $v_{k,0}$ is a representative velocity (freestream velocity in this case), and $\mathbf{u}_{k}$ is a dimensionless representation of velocity. This allows measurement with respect to the characteristic values, so $\mathbf{u}_{k}$ is 1 in the freestream, and greater than 1 in regions in which it has experienced an acceleration, respectively. The full list of substitutions is given here:

	$\displaystyle\mathbf{x}=$	$\displaystyle r_{0}\mathbf{X}$	$\displaystyle t=$	$\displaystyle\tau/\omega_{0}$
	$\displaystyle\mathbf{c}_{k}=$	$\displaystyle v_{k,0}\mathbf{u}_{k}$	$\displaystyle f_{k}$	$\displaystyle=n_{k,\infty}F_{k},$
	$\displaystyle\phi$	$\displaystyle=\phi_{0}\mathbf{\Phi},$	$\displaystyle q_{k}$	$\displaystyle=q_{e}Z_{k}$
	$\displaystyle n_{k}=$	$\displaystyle n_{k,\infty}N_{k},$	$\displaystyle\nabla_{x}=$	$\displaystyle\hat{\nabla}_{X}/r_{0}$
	$\displaystyle\nabla_{c}=$	$\displaystyle\hat{\nabla}_{u}/v_{k,0}$

For $r_{0}$, $\omega_{0}$, $v_{0}$, $\phi_{0}$ as characteristic length, frequency, velocity, and potential respectively. $\mathbf{X},\tau,\mathbf{u},F,\Phi,Z,$ and $N_{k}$ are dimensionless parameters, and $\nabla_{x}$, $\nabla_{u}$ are the dimensionless gradient operators in physical and velocity space, respectively. Substituting these transformations into Eqn.4.4, Eqn4.5 gives

$$\left(\frac{\omega_{0}r_{0}}{v_{k,0}}\right)\frac{\partial F_{k}}{\partial t}+\mathbf{u}_{k}\cdot\hat{\nabla}_{X}F_{k}-Z_{k}\left(\frac{q_{e}\phi_{0}}{m_{k}v^{2}_{k,0}}\right)\hat{\nabla}_{X}\Phi\cdot\hat{\nabla}_{u}F_{k}=0$$

(4.6)

$$\hat{\nabla}^{2}_{X}\Phi=\sum^{K}_{k}-Z_{k}\frac{r^{2}_{0}q_{e}n_{k,\infty}}{\epsilon_{0}\phi_{0}}\left(N_{k}-\exp{\left[\left(\frac{q_{e}\phi_{0}}{k_{B}T_{e}}\right)\Phi\right]}\right)$$

(4.7)

The ion drift ratio is introduced as $S_{k}=v_{k,0}/v_{k,t}$, where $v_{k,0}$ is the ion drift velocity (the relative velocity between the charged body and the ion freestream), and $v_{k,t}$ is the ion thermal velocity, given by [45]

$$v_{k,t}=\sqrt{\frac{2k_{B}T_{k,\infty}}{m_{k}}}$$

(4.8)

allowing us to write Eqn.4.6 as

$$\left(\frac{\omega_{0}r_{0}}{v_{k,B}}\right)\frac{\partial F_{k}}{\partial t}+(1+S_{k}^{-1})\mathbf{u}_{k}\cdot\hat{\nabla}_{X}F_{k}-Z_{k}\left(\frac{q_{e}\phi_{0}}{m_{k}v^{2}_{k,0}}\right)\frac{1}{(1+S_{k}^{-1})}\hat{\nabla}_{X}\Phi\cdot\hat{\nabla}_{u}F_{k}=0$$

(4.9)

From inspection, there are 5 dimensionless parameters that govern the behaviour of Eqns 4.9-4.7,

	$\displaystyle\alpha_{k}=$	$\displaystyle-Z_{k}\left(\frac{q_{e}\phi_{0}}{m_{k}v_{k,B}^{2}}\right),$	$\displaystyle\mu_{e}=$	$\displaystyle\left(\frac{q_{e}\phi_{0}}{k_{B}T_{e}}\right)$
	$\displaystyle\xi_{k}=$	$\displaystyle-Z_{k}\left(\frac{r^{2}_{0}n_{k,\infty}q_{e}}{\epsilon_{0}\phi_{0}}\right),$	$\displaystyle\Omega_{k}=$	$\displaystyle\left(\frac{\omega_{0}r_{0}}{v_{k,B}}\right)$		(4.10)
	$\displaystyle S_{k}=$	$\displaystyle v_{k,B}\sqrt{\frac{m_{k}}{2k_{B}T_{k,\infty}}}$

substituting the dimensionless parameters in Eqn.4.10 into Eqn.4.7 and 4.9 gives

$$\Omega_{k}\frac{\partial F_{k}}{\partial\tau}+(1+S_{k}^{-1})\mathbf{u}_{k}\cdot\hat{\nabla}_{X}F_{k}+\frac{\alpha_{k}}{(1+S_{k}^{-1})}\hat{\nabla}_{X}\Phi\cdot\hat{\nabla}_{u}F_{k}=0$$

(4.11)

$$\hat{\nabla}_{X}^{2}\Phi=\sum^{K}_{k}\xi_{k}(N_{k}-\exp[\mu_{e}\Phi])$$

(4.12)

Reducing the 4+5K quantities ($\epsilon_{0}$, $q_{e}$, $T_{e}$, $k_{0}$, as well as separate $m_{k}$, $q_{k}$, $T_{k}$ and $n_{k}$ for each ion species k) to 1+4K dimensionless parameters ($\mu_{e}$, plus separate $\alpha_{k}$, $\xi_{k}$, $S_{k}$, and $\Omega_{k}$ for each ion species k. The relative effect of each species on the potential $\Phi$ can be made more explicit by introducing the dimensionless parameters $\chi$ and $B_{k}$

$$\chi=\sqrt{\sum^{K}_{k}\xi_{k}}$$

(4.13)

$$\beta_{k}=\frac{\xi_{k}}{\chi^{2}}$$

(4.14)

Allowing Eqn.4.12 to be expressed as

$$\hat{\nabla}_{X}^{2}\Phi=\chi^{2}\left(\sum^{K}_{k}(\beta_{k}N_{k})-\exp[\mu_{e}\Phi]\right)$$

(4.15)

The physical meaning, and relative significance of these various parameters is laid out below

Mathematically, the independent dimensionless parameters in (4.10) form a complete set of $\Pi$ groups to describe the K-species system of Vlasov-Maxwell equations [22] (Eqn. 4.4 and Eqn.4.5). A physical intuition for each of these parameters is laid out here, with special attention paid to $\alpha_{k}$,$\beta_{k}$,$\chi$, the parameters of interest for this study.

In Capon et al [27], a response surface is developed to find drag as a function of $\alpha_{k}$ and $\chi$. While this response surface is only designed for single species flows, it provides a physics based function to equate wake features to drag coefficients, and therefore an important basis for study of mixed flows. Findings by Capon et al[20] also indicated that the bulk of the variation in drag forces is accounted for by variation in $\alpha_{k}$ and $\chi$, and proposed a response surface to capture this variation as a function of the form

$$C_{D,C}=f(\alpha_{k})+f(\chi)+f(\alpha_{k},\chi)$$

(4.21)

Where $C_{D,C}$ was a drag coefficient representing both charged and neutral drag forces. $f(\alpha_{k})$ and $f(\chi)$ are intended to capture variations in $C_{D,C}$ from $\alpha$ and $\chi$ respectively, with the third term accounting for coupled effects. The form of the functions for $\alpha_{k}$ and $\chi$ are given as

$$f(\alpha_{k})=\mathscr{A}(1+2\alpha_{k})^{0.5-a},\quad f(\chi)=\mathscr{B}\chi^{-1}$$

(4.22)

The term $f(\alpha_{k})$ in this has an analogous form to the orbital motion limited impact parameter [47]. In this case, $\mathscr{A}$ is a geometry dependent constant, while $a$ is to account for the reduction in direct drag caused by ion accelerations. The second term can be made more clear by substituting $\chi=\frac{r_{0}}{\lambda_{\phi}}$, where $\lambda_{\phi}$ is given by Eqn.4.19 into the Child Langmuir law for sheath thickness ($d_{sh}$) [21] around a high voltage cathode, given by

$$d_{sh}=\frac{2^{5/4}}{3}\sqrt{\frac{\epsilon_{0}\phi_{0}}{q_{i}n_{i}}}=\frac{2^{5/4}}{3}\lambda_{\phi}=\frac{2^{5/4}r_{0}}{3}\chi^{-1}$$

(4.23)

thus the term $\chi^{-1}$ is proportional to the ratio of sheath thickness to body size, $d_{sh}/r_{0}$. In the case of a drag coefficient, this is effectively a measurement of the collection surface of the charged body, roughly analogous to the relationship of incident surface area to $C_{D}$ in a neutral flow. The coupling term in Eqn.4.22 is given by

$$f(\alpha_{k},\chi)=\mathscr{C}\frac{2\alpha_{k}^{c}+\mathscr{D}}{1+\chi}$$

(4.24)

where $\mathscr{B},\mathscr{C},\mathscr{D},$ and $c$ are all fitted constants. The response surface was formulated by fitting the constants in Eqns.4.22-(4.24) to the drag coefficients of a range of single species, mesothermal plasma flows, giving an empirical response surface (Fig.4.1)