Tutorial 6: Calculating the density and heat capacity of an O2 plasma

The most common use of minplascalc is the calculation of thermophysical properties of plasmas in LTE. In this example we’ll look at the thermodynamic properties \(\rho\) and \(C_P\).

The relatively simple case of a pure oxygen plasma is useful for demonstration and validation purposes. As in Tutorial 5, a JSON mixture file must be created by the user to specify the plasma species present and the relative proportions of elements. We’ll use the example mixture file shipped with minplascalc, located at minplascalc/data/mixtures/O2_plasma.json.

We start by loading up the matplotlib graphing module, numpy for array operations, and minplascalc.

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import minplascalc as mpc

Next, we create a minplascalc Mixture object using the path to the mixture file as before.

mixture = mpc.Mixture(str(mpc.MIXTUREPATH / "O2_plasma.json"), 1000, 101325)

Next, set a range of temperatures to calculate the thermophysical properties at - in this case we’re going from 1000 to 25000 K. Also initialise lists to store the property values at each temperature

temperatures = np.linspace(1000, 25000, 100)
density = []
cp = []

Now we can perform the property calculations. We loop over all the temperatures setting the mixture object’s temperature attribute to the appropriate value, initialising the species numbers using the Mixture object’s initialiseNi(Ni) function (which takes a list of the initial guesses of the \(N_i\) value for each species), and calculating the LTE composition using the object’s solveGfe() function. Once we have the LTE composition, we can calculate the plasma density by calling the Mixture object’s calculateDensity() function.

The Mixture object also has a calculate_heat_capacity() function. Internally, this makes two additional function calls to solveGfe() at temperatures slightly above and slightly below the target temperature. A numerical derivative then gives the LTE value of \(C_P\) at that point.

Note that execution of these calculations is fairly compute intensive and as a result the following code snippet may take several seconds to complete.

for T in temperatures:
    mixture.T = T
    mixture.initialise_ni([1e20 for n in range(len(mixture.species))])
    mixture.solve_gfe()
    density.append(mixture.calculate_density())
    cp.append(mixture.calculate_heat_capacity())

Now we can visualise the properties by plotting them against temperature, to see how they vary.

fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(10, 5))
ax0.set_title("$\mathregular{O_2}$ plasma density at 1 atm")
ax0.set_xlabel("T, K")
ax0.set_ylabel("$\mathregular{\\rho, kg/m^3}$")
ax0.semilogy(temperatures, density);
ax1.set_title("$\mathregular{O_2}$ plasma heat capacity at 1 atm")
ax1.set_xlabel("T, K")
ax1.set_ylabel("$\mathregular{C_P, J/kg \cdot K}$")
ax1.plot(temperatures, cp);
plt.tight_layout()

The results obtained using minplascalc are comparable to other data for oxygen plasmas in literature, for example Boulos et al 1994 (see README for full reference). In particular the position and size of the peaks in \(C_P\), which are caused by the highly nonlinear dissociation and first ionisation reactions of O2 and O respectively, are accurately captured.