—

Use Stine-Watson model for a DC thermal plasma in H₂.

This example demonstrates how to use the Stine-Watson model to compute the thermal efficiency, plasma radius, and voltage of a thermal plasma torch.

The model is based on the work of [Stine1962], and is extended to include radiation effects. It can describes the temperature distribution in a DC plasma torch, with gas flow stabilized by water-cooled walls.

It is based on the following assumptions:

• Typical assumptions for a DC thermal plasma torch:

  • The arc is in local thermodynamic equilibrium.

  • The arc is in steady-state and is axially symmetrical.

  • Thermodiffusion effects are neglected.

  • Gravity and viscous dissipation effects are negligible.

  • Thermodynamic and transport properties depend only on temperature (negligible pressure gradients).

  • Negligeable magnetic field.

  • No effect of charged particles (\(q = 0\)) / neutral plasma.

  • Current is assumed to be constant along the arc.

• Strong assumptions ((@reader: to be discussed) !)

  • Flows are laminar (\(\vec{v} = v_z \hat{e}_z\)).

  • There are two zones, as represented on Figure 3:

    • a hot channel, where Joule Power is deposited, and heat is evacuated through radiation and conduction, as well as in the augmentation of enthalpy of the gas. It is assumed that this hot region is a cylinder of length \(L\) and of radius \(r_0\).

    • a cold channel, which surrounds the hot one, where there is only heat conduction and change of enthalpy (no radiation and no Joule power in this region). It is assumed that this cold region is a cylinder with a hole, of length \(L\), of inner radius \(r_0\) and outer radius \(R\).

  • Furthermore, it will be assumed that the arc length L is larger than the torch radius R (L≫R). It should also be noted that the arc length L is a input of the model.

• Assumptions to close the system:

  • On ratio of flows: The total mass flow is \(\dot{m}\) is known, what is not known is the ratio of the mass flow in the hot region \(\dot{m}_h\) over the total mass flow \(\dot{m}\). This ratio is denoted by \(\alpha\) and is a parameter of the model. It is assumed that the mass flow in the cold region is \(\dot{m}_c = \dot{m} - \dot{m}_h\).

  • On the radius of the hot channel \(r_0\): The radius of the hot channel is not known, but it is determined by continuity of thermal fluxes at the interface between the hot and cold regions.

Notes

The equation solved is the following:

\[\begin{split}\rho_m v_z \frac{\partial h}{\partial z} = \frac{1}{r} \frac{\partial}{\partial r} \left(r \kappa(T) \frac{\partial T}{\partial r}\right) +\frac{\partial}{\partial z}\left(\kappa(T) \frac{\partial T}{\partial z}\right) \\ +\frac{\sigma(r, z) I^2}{\left(\int_0^R 2 \pi r \sigma(r, z) d r\right)^2} -\vec{\nabla} \cdot \overline{q^{r a d}}\end{split}\]

where:

• \(\rho_m\) is the mass density of the gas [kg/m³],

• \(v_z\) is the axial velocity of the gas [m/s],

• \(h\) is the enthalpy of the gas [J/kg],

• \(T\) is the temperature of the gas [K],

• \(\kappa(T)\) is the thermal conductivity of the gas [W/m/K],

• \(\sigma(T)\) is the electrical conductivity of the gas [S/m],

• \(I\) is the current [A], assumed to be constant along the arc,

• \(\vec{\nabla} \cdot \overline{q^{r a d}}\) is the radiative power [W/m^3].

The radiative power is computed using the Net Emission Coefficient ([Gueye2017]), i.e. \(\vec{\nabla} \cdot \overline{q^{r a d}} = 4 \pi \varepsilon_N(\theta)\).

Introducing the integrated thermal conductivity \(\Theta(T) = \int_0^T \kappa(s) ds\), and using geometry simplification, the equation can be rewritten as:

\[\rho_m v_z \frac{\partial h(\theta)}{\partial z} = \frac{1}{r} \frac{\partial}{\partial r}\left(r \frac{\partial \theta}{\partial r}\right) +\frac{\sigma(\theta) I^2}{\left(\int_0^R 2 \pi r \sigma(r, z) d r\right)^2} -4 \pi \varepsilon_N(\theta)\]

The following approximations can then be made on thermal conductivity, enthalpy and electrical conductivity:

\[\begin{split}\begin{gathered} \left(\rho_m v_z\right)(r)=\left\{\begin{array}{c} \alpha \frac{\dot{m}^{\text {tot }}}{\pi r_0^2} \text { for } 0<r<r_0 \\ (1-\alpha) \frac{\dot{m}^{\text {tot }}}{\pi R^2-\pi r_0^2} \text { for } r_0<r<R \end{array}\right. \\ \sigma(\theta)= \begin{cases}0 & \text { for } \theta<\theta_\sigma \\ a_\sigma\left(\theta-\theta_\sigma\right) & \text { for } \theta>\theta_\sigma\end{cases} \\ h(\theta)=a_h \theta \end{gathered}\end{split}\]

Therefore,

• In the cold region (\(r_0<r<R\)):

  \[(1-\alpha) \frac{\dot{m}^{t o t}}{\pi R^2-\pi r_0^2} a_h \frac{\partial \theta}{\partial z} =\frac{1}{r} \frac{\partial}{\partial r}\left(r \frac{\partial \theta}{\partial r}\right)\]

• In the hot region (\(0<r<r_0\)):

  \[\alpha \frac{\dot{m}^{t o t}}{\pi r_0^2} a_h \frac{\partial \theta}{\partial z} =\frac{1}{r} \frac{\partial}{\partial r}\left(r \frac{\partial \theta}{\partial r}\right) +\frac{a_\sigma\left(\theta-\theta_\sigma\right) I^2} {\left(\int_0^{r_0} 2 \pi r a_\sigma\left(\theta-\theta_\sigma\right) d r\right)^2} -a_{\varepsilon}\left(\theta-\theta_\sigma\right)\]

This can be solved analytically with the following boundary conditions:

• At the radial boundary of the hot channel, the conductivity starts to increase, i.e. is no more zero. Said otherwise, the radius of the channel r_0 is such that \(\theta(r=r_0, z)=\theta_{\sigma}\).

• At the outer radius of the cold channel, the wall is cooled by water, and the temperature is assumed to be null and constant, i.e. \(\theta(r=R, z)=0\).

• The origin of the z axis is defined to be the boundary of the hot channel: \(\forall 0<r<r_0, \theta(r, z=0)=\theta_{\sigma}\)

• For the cold region, the temperature is assumed to be constant along the z axis: \(\forall r_0<r<R, \frac{\partial \theta}{\partial z}(r, z)=0\).

Tags: plasma Stine-Watson thermal efficiency voltage plasma radius H2 hydrogen CO2 CH4 thermal plasma DC plasma VAC flow rate

# This is an option for the online documentation, so that each image is displayed separately.
# sphinx_gallery_multi_image = "single"

Import the required libraries.

import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns

from rizer.thermal_plasma.fit_LTE_data import FitLTEData
from rizer.thermal_plasma.stine_watson import StineWatson

# Set the style of the plots.
sns.set_theme("poster")

Get the Stine-Watson coefficent for each gas, at different pressures.

# sw_data = FitLTEData(
#     gas_name_transport="H2",
#     gas_name_radiation="H2",
#     pressure_atm=1,
#     source_transport="Boulos2023",
#     source_radiation="Gueye2017",
#     emission_radius_mm=10,
#     max_temperature_fit=12_000,  # K,
# )

sw_data = FitLTEData(
    gas_name_transport="80CO2_20CH4",
    gas_name_radiation="CO2",
    pressure_atm=1,
    source_transport="Wu2016_Niu2016",
    source_radiation="Billoux2012",
    emission_radius_mm=5,
    max_temperature_fit=12_000,  # K,
)

a_sigma, theta_sigma, a_h, a_epsilon = sw_data.fit_all()
print(f"a_sigma = {a_sigma:.4e} (S/m)/(W.m^-1)")
print(f"theta_sigma = {theta_sigma:.4e} W.m^-1")
print(f"a_h = {a_h:.4e} (J/kg)/(W.m^-1)")
print(f"a_epsilon = {a_epsilon:.4e} (W.m^-3.sr^-1)/(W.m^-1)")


gas_nice = {
    "H2": r"$\mathregular{H_2}$",
    "80CO2_20CH4": r"$\mathregular{80\%CO_2-20\%CH_4}$",
}

• 
• 
• 

/home/runner/work/rizer/rizer/rizer/thermal_plasma/fit_LTE_data.py:321: UserWarning: In `FitLTEData.fit_nec`, using previously fitted theta_sigma.
  warn("In `FitLTEData.fit_nec`, using previously fitted theta_sigma.")
a_sigma = 4.1068e-01 (S/m)/(W.m^-1)
theta_sigma = 7.8882e+03 W.m^-1
a_h = 4.4774e+03 (J/kg)/(W.m^-1)
a_epsilon = 6.6216e+04 (W.m^-3.sr^-1)/(W.m^-1)

Define the parameters.

# Define the flow rate.
m_dot_kg_per_day = 50.4  # [kg/day], mass flow rate.
m_dots = [m_dot_kg_per_day / (24 * 3600)]  # [kg/s], mass flow rate.
if sw_data.thermo_transport_data._gas_name == "80CO2_20CH4":
    V_dot_m3_per_h = np.array([70, 100])  # [m^3/h], volume flow rate.
    rho = 1.56  # [kg/m^3], density of 80CO2-20CH4 at 300K, 1atm.
    m_dots = [
        V_dot
        / 3600  # [m^3/s], volume flow rate.
        * rho
        * (273 / 300)
        for V_dot in V_dot_m3_per_h
    ]  # [kg/s], mass flow rate.

# Define the torch and plasma geometry.
Ls = [0.2, 0.3, 0.4]  # [m], plasma length, value for BO300.
R = 0.025  # [m], torch radius, value for BO300.
alpha = 0.65  # [0-1] parameter to control the mass flow going in the hot region.
if alpha < 0 or alpha > 1:
    raise ValueError("Alpha must be between 0 and 1.")

# Define the current, supposed to be constant accross the plasma.
currents = np.arange(100, 310, 10, dtype=float)  # [A]

# Plot options
colors = ["blue", "red"]
linestyles = ["-.", "--", "-"]
assert len(colors) == len(m_dots)
assert len(linestyles) == len(Ls)

Solve the Stine-Watson model.

# Plot plasma radius as a function of current.
fig_r0, ax_r0 = plt.subplots(figsize=(12, 8), layout="constrained")
# Plot thermal efficiency as a function of current.
fig_eff, ax_eff = plt.subplots(figsize=(12, 8), layout="constrained")
# Plot total voltage as a function of current.
fig_v, ax_v = plt.subplots(figsize=(12, 8), layout="constrained")

for mass_flow_rate, color in zip(m_dots, colors):
    for L, linestyle in zip(Ls, linestyles):
        # Computing with radiation.
        plasma_radius = np.zeros_like(currents)
        thermal_efficencies = np.zeros_like(currents)
        voltages = np.zeros_like(currents)

        for i, current in enumerate(currents):
            # Define the Stine-Watson model.
            sw_model = StineWatson(
                torch_radius=R,
                arc_length=L,
                current=current,
                total_mass_flow=mass_flow_rate,
                ratio_of_hot_mass_flow_over_total_mass_flow=alpha,
                gas_data=sw_data,
            )

            # Get plasma radius, thermal efficiency and voltage.
            plasma_radius[i] = sw_model.r_0
            thermal_efficencies[i] = sw_model.thermal_efficiency()
            voltages[i] = sw_model.total_voltage()

        ax_r0.plot(
            currents,
            plasma_radius * 1e3,
            label=f"L={L} m, ṁ={mass_flow_rate * 3600 * 24:.1f} kg/day",
            color=color,
            linestyle=linestyle,
        )
        ax_eff.plot(
            currents,
            thermal_efficencies * 100,
            label=f"L={L} m, ṁ={mass_flow_rate * 3600 * 24:.1f} kg/day",
            color=color,
            linestyle=linestyle,
        )
        ax_v.plot(
            currents,
            voltages,
            label=f"L={L} m, ṁ={mass_flow_rate * 3600 * 24:.1f} kg/day",
            color=color,
            linestyle=linestyle,
        )


# Set the labels and title.
properties = f"R = {R} m\n" + r"$\alpha$ = m_hot/m_tot = " + f"{alpha}"
# .. for plasma radius.
fig_r0.suptitle(
    f"Plasma radius [mm] versus current [A], for {gas_nice[sw_data.thermo_transport_data._gas_name]} "
    f"at {sw_data.thermo_transport_data._pressure_atm} atm",
    fontsize=28,
)
ax_r0.set_title(properties, fontsize=14, loc="left")
ax_r0.set_xlabel("Current [A]")
ax_r0.set_ylabel("Plasma radius [mm]")
ax_r0.legend()


# .. for thermal efficiency.
if "80CO2_20CH4" == sw_data.thermo_transport_data._gas_name:
    europlasma_data = np.array(
        [
            [241.78, 82.6],  # [A, %]
            [291.72, 81.9],  # [A, %]
        ]
    )
    europlasma_currents = europlasma_data[:, 0]
    europlasma_thermal_efficiencies = europlasma_data[:, 1]
    ax_eff.plot(
        europlasma_currents,
        europlasma_thermal_efficiencies,
        "r*",
        ms=20,
        label="EuroPlasma (BO300, 100 Nm3/h)",
    )

    europlasma_data_2 = np.array(
        [
            [142.65, 87.6],  # [A, %]
            [192.58, 85.0],  # [A, %]
            [292.46, 81.0],  # [A, %]
        ]
    )
    europlasma_currents_2 = europlasma_data_2[:, 0]
    europlasma_thermal_efficiencies_2 = europlasma_data_2[:, 1]
    ax_eff.plot(
        europlasma_currents_2,
        europlasma_thermal_efficiencies_2,
        "bv",
        ms=20,
        label="EuroPlasma (BO300, 70 Nm3/h)",
    )
fig_eff.suptitle(
    f"Thermal efficiency versus current [A], for {gas_nice[sw_data.thermo_transport_data._gas_name]} "
    f"at {sw_data.thermo_transport_data._pressure_atm} atm",
    fontsize=28,
)
ax_eff.set_title(properties, fontsize=14, loc="left")
ax_eff.set_xlabel("Current [A]")
ax_eff.set_ylabel("Thermal efficiency [%]")
ax_eff.legend()

# .. for total voltage.
if "80CO2_20CH4" == sw_data.thermo_transport_data._gas_name:
    europlasma_data = np.array(
        [
            [171.90, 1914.5],  # [A, V]
            [192.02, 1787.4],
            [211.90, 1585.3],
            [212.15, 1702.6],
            [241.60, 1637.5],
            [251.66, 1569.0],
            [271.53, 1520.2],
            [291.66, 1494.1],
        ]
    )
    europlasma_currents = europlasma_data[:, 0]
    europlasma_voltages = europlasma_data[:, 1]
    ax_v.plot(
        europlasma_currents,
        europlasma_voltages,
        "r*",
        ms=20,
        label="EuroPlasma (BO300, 100 Nm3/h)",
    )

    europlasma_data_2 = np.array(
        [
            [142.45, 1683.1],  # [A, V]
            [153.01, 1611.4],
            [162.33, 1556.0],
            [171.90, 1513.6],
            [191.78, 1428.9],
            [212.39, 1347.5],
            [232.27, 1292.1],
            [242.58, 1259.5],
            [252.88, 1233.4],
            [272.52, 1194.3],
            [292.39, 1155.2],
        ]
    )
    europlasma_currents_2 = europlasma_data_2[:, 0]
    europlasma_voltages_2 = europlasma_data_2[:, 1]
    ax_v.plot(
        europlasma_currents_2,
        europlasma_voltages_2,
        "bv",
        ms=20,
        label="EuroPlasma (BO300, 70 Nm3/h)",
    )

fig_v.suptitle(
    f"Total voltage [V] versus current [A], for {gas_nice[sw_data.thermo_transport_data._gas_name]} "
    f"at {sw_data.thermo_transport_data._pressure_atm} atm",
    fontsize=28,
)
ax_v.set_title(properties, fontsize=14, loc="left")
ax_v.set_xlabel("Current [A]")
ax_v.set_ylabel("Total voltage [V]")
ax_v.legend()

plt.show()

• 
• 
• 

Plot the temperature map.

# Set the parameters for the model.
sw_model.I = 100.0  # [A]

# Plot the theta map.
sw_model.plot_theta_map()

# Plot the temperature map.
sw_model.plot_temperature_map()

• 
•