Starfish Tutorial Part 5: Collisions and Chemical Reactions

Posted on December 12th, 2012
Previous Article :: Next Article

Introduction

Welcome to the final part of the Starfish tutorial. In the past lessons, we set up the computational domain, specified the potential solver parameters, added ion source, and included surface recombination. However, we have not yet considered gas interactions. The individual ion and neutral particles were able to pass right through without “seeing” each other. We will correct this omission in this part.

Starfish supports three types of material interactions: chemical reactions, MCC, and DSMC. The easiest way to differentiate between these is to think of them as fluid-fluid, particle-fluid, and particle-particle events, respectively. Chemical reactions operate solely with the density and temperature fields and are applicable to models described by a rate equation. They are useful for modeling production or destruction of material in processes such as ionization or recombination.

MCC, or Monte Carlo Collisions, are kinetic-fluid interactions. The source material collides with a target cloud. The number density of the target at the particle position is used to determine the collision probability. If the collision occurs, only the source particle is modified. The target is not affected by the collision. As such, momentum is not conserved. MCC is suitable for cases when the target material is sufficiently more dense than the source, such as when a rarefied ion beam interacts with a dense neutral cloud via the Charge Exchange (CEX) collision.

Finally, DSMC (Direct Simulation Monte Carlo) is a kinetic-kinetic collision process. This method collides particles with other particles in the simulation cell. Both energy and momentum are conserved. This method is suitable for modeling collisions in like gases, such as to model momentum exchange (MEX) collisions in the ion gas. Starfish implements the No Time Counter (NTC) method of Bird[1]. DSMC is the most computationally demanding of these three methods.

Specifying Interactions

Material interactions are specified in the interactions file which you have already seen before. We used this file previously to add surface recombination. We modify the interactions file as follows:

<!-- material interactions file -->
<material_interactions>
 
<!-- surface interaction -->
<surface_hit source="O+" target="SS">
<product>O</product>
<model>cosine</model>
<prob>0.9</prob>
<c_accom>0.5</c_accom>
<c_rest>0.9</c_rest>
</surface_hit>
 
<!-- chemical reaction -->
<chemistry model="ionization">
<sources>O,e-</sources>
<products>iO+,2*e-</products>
</chemistry> 
 
<!-- MCC charge exchange -->
<mcc model="cex">
<source>O+</source>
<target>O</target>
<sigma>const</sigma>
<sigma_coeffs>1e-18</sigma_coeffs>
</mcc> 
 
<!-- DSMC momentum exchange -->
<dsmc model="vhs">
<source>O+</source>
<target>O+</target>
<sigma>const</sigma>
<sigma_coeffs>1e-18</sigma_coeffs>
</dsmc> 
 
</material_interactions>

We have added three new interactions types to the list: chemistry, mcc, and dsmc. Of course, you could repeat these definitions to specify multiple collision pairs. We will now go through each interaction type in detail.

Chemical Reactions

The chemistry interaction type models interactions that can be written in terms of density change. As an example, consider the ionization process:
$$O+e^- \rightarrow O^+ +2e^-$$

We can write this process as follows:
$$\begin{array}{rcl}
\text{d}n_{O^+} &=& +k n_{O} n_{e^-} \text{d}t \\
\text{d}n_{O} &=& -k n_{O^+} n_{e^-} \text{d}t
\end{array}$$
ignoring electron density change. Here “k” is the ionization rate (which is typically function of the electron temperature) and the two “n” correspond to the densities of the atoms and electrons. The rate “k” is model-specific, and here we use an ionization rate from the dissertation of Fife[2]. Let’s take a look at what happens when we add

<chemistry model="ionization">
<sources>O,e-</sources>
<products>iO+,2*e-</products>
</chemistry>

to the interactions file. Since this reaction requires neutrals and electrons, we would expect the ionization to occur only in the region where neutrals are present – i.e. where the ions have recombined with the wall. This is indeed the case as we can see from Figures 1 and 2 below. Figure 1 shows the densities of the two source materials, while Figure 2 shows the density of the ions created by this reaction. As you have probably noticed above, we specified the ionization product to be “iO+”. This is simply a new material that was created with properties identical to the regular ions, and is being solely for diagnostics, to differentiate these newly born ions from the ones that are injected by the source. This is even more important in this case, since ionization is simply not an important process at the low densities encountered in this tutorial. We can see that the density of these new ions is 8 orders of magnitude smaller than the density of the primary ions!

neutral density electron density
Figure 1. Number density of neutrals and electrons
ionization density
Figure 2. Number density of ions created by the ionization reaction

You may have some questions at this point. Such as, where did the electrons come from? If you remember, at no point in this tutorial did we add electrons to the simulation. They were not needed, since we were using the hybrid model, in which electrons are assumed to follow the Boltzmann relationship and the potential is solved with the non-linear Poisson solver. Now, you could have specified electrons as actual particles and simulated them with the kinetic particle in cell method. But in case like ours, Starfish automatically creates a fluid electron species and adds them to the simulation. They can be then used as an actor in chemical reactions, or the target in MCC. The density of the electrons follows the model used by the potential solver, in this case, we are using the isothermal Boltzmann relationship, \(n_e = n_0 \exp [(\phi-\phi_0)/kT_{e,0}]\). The plot in Figure 1 shows exactly this expression.

Another question you may have is, how is the temperature calculated? In this case the electrons were isothermal, but other models are available as well. As an aside, in the case of general fluid materials, the temperature is computed from the energy solver. In the case of kinetic particles, velocity is computed from the root mean squared of the mean cell velocity variation, \(\sqrt{(1/n)(g_1^2+g_1^2+\ldots g_n^2} = \sqrt{3kT/m}\), where \(g_i=|v_i-\bar{v}|\). You can see an example of this below in Figure 3. The flood shows the computed temperature of the ions created by the ionization process. In our definition, temperature basically indicates how much each particle deviates from the average bulk flow. Well we can see that there is very little variation (the white region below lower limit cut off) until the ions get to the wake region. Here we end up with mixing of the two flows from both sides of the sphere, which then correspondingly gives us large variation in particle velocities, and thus increased temperature. The colored contour lines correspond to the potential and the solid black lines are the velocity streamlines. You can see that now a small potential hill forms that reflects some of the incoming ions.

ion temperature
Figure 3. Temperature variation of the ions created by the ionization reaction. Contour lines correspond to the potential and the black lines are velocity streamlines.

Finally, you may also be wondering how the change in density works with kinetic materials. The “dn” is converted into a volumetric source. After the density changes are computed for all chemical reactions, Starfish loops through the cells and creates in each cell as many particles as needed to satisfy the new density increase. The particles are created randomly in the cell and have initial velocity corresponding the Maxwellian distribution at the target gas temperature. This process will create particles with specific weights less than the default one prescribed in the materials.xml file. Similarly, negative densities are used to destroy particles. The density is converted to a number loss, and random particles from the cell are destroyed until enough are removed. If the remaining mass loss is less than a mass of a single macroparticle, the specific weight of a particle is reduced accordingly.

MCC

Let’s now move onto MCC, which was described in more detail in a previous article. To specify MCC interactions, we add the following to the list of interactions

<mcc model="cex">
<source>O+</source>
<target>O</target>
<sigma>const</sigma>
<sigma_coeffs>1e-18</sigma_coeffs>
</mcc>

Just as in the case of the chemical reaction, we specify the source and target, as well as the collision model. In this case, we’ll be using the CEX handler. This handler models the electron exchange between a fast ion and a slow neutral resulting in a fast neutral and a slow ion,
$$O^+_{fast} + O_{slow} \rightarrow O_{fast} + O^+_{slow}$$

But since this is MCC, we don’t actually modify the neutrals. Instead, the neutrals are used to compute the collision probability following \(P=1-\exp(-n_ng\sigma\Delta t)\), where \(n_n\) is the neutral density, \(g\) is the relative velocity between the ion and the neutral (with the neutral assumed to be stationary), and \(\sigma\) is the collision cross-section. The sigma model is typically a function of relative velocity, with a number of models existing to describe different collision events. A classic model for modeling CEX is the model of Rapp and Francis. However here, for simplicity, we use a constant cross-section. But even with this relatively large value, collisions are still going to be a very rare event due to the low gas densities. So just to demonstrate this effect, let’s go ahead and specify a dense background neutral environment. In the materials file, let’s modify the oxygen atom by adding an init tag,

<material name="O" type="kinetic">
<molwt>16</molwt>
<charge>0</charge>
<spwt>2e5</spwt>
<init>nd_back=2e18</init>
</material>

This background density will be added to any density from the actual kinetic particles. We can also turn off the potential solver. One way is to replace the Poisson solver with a constant electric field model with zero components,

<!-- set zero electric field -->
<solver type="constant-ef">
<comps>0,0</comps>
</solver>

This will allow us to see the effect of collisions. The source loads a cold ion beam, and hence, in the absence of forces and collisions, the ions should continue moving in a straight line. Collisions will scatter the motion. You can see this comparison for yourself below in Figure 4. You can see that once the collisions are enabled, we both start seeing diffusion of ions into the wake behind the cylinder, and also the overall density of ions increases. This is due to the presence of many slow ions that are taking a long time to leave the simulation.

ion density, no collisions ion density, CEX collisions
Figure 4. Ion density without (left) and with (right) CEX collisions enabled

DSMC

DSMC is demonstrated in more detail in a separate tutorial.

References

[1] Bird, G.A. Molecular Gas Dynamics and the Direct Simulation of Gas Flows, Oxford Science Publications, 2003
[2] Fife, J., Hybrid-PIC Modeling and Electrostatic Probe Survey of Hall Thrusters, Ph.D. Dissertation, Massachusetts Institute of Technology, 1998.

Subscribe to the newsletter and follow us on Twitter. Send us an email if you have any questions.

Leave a Reply

You can use these tags: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong> <pre lang="" line="" escaped="" cssfile=""> In addition, you can use \( ...\) to include equations.