3. Plasma Simulation

This document provides comprehensive examples of direct plasma physics simulations using classical computational methods (finite volume, finite difference) for arc discharge and other plasma phenomena.

3.1. Overview

Plasma Physics Simulations use traditional numerical methods to solve the governing equations of plasma dynamics. Unlike Physics-Informed Neural Networks (PINNs), these approaches discretize the physical domain and equations directly, providing:

High accuracy: Well-established numerical methods with rigorous error bounds
Physical validation: Direct comparison with experimental measurements
Computational benchmarks: Reference solutions for validating ML-based approaches
Industrial relevance: Proven methods used in circuit breaker and plasma device design

3.1.1. Why Classical Plasma Simulations?

Classical simulations serve multiple purposes in AI4Plasma:

Ground Truth Generation: Provide reference data for training neural operators
Method Validation: Benchmark PINN and operator learning results
Physical Understanding: Reveal multi-scale phenomena and parameter sensitivities
Engineering Design: Established tools for industrial applications

3.1.2. Available Implementations

AI4Plasma provides classical simulation tools for:

Arc Discharge: Thermal plasma arcs in circuit breakers and switching devices
Corona Discharge: Non-equilibrium plasmas near sharp electrodes (future)
Streamer Physics: Ionization front propagation (future)

3.2. Arc Discharge Simulations

3.2.1. Physical Background

Arc discharge is a high-temperature, conducting plasma column occurring when electric current flows through ionized gas. Applications include:

Circuit Breakers: Arc quenching after current interruption
Welding: High-temperature plasma for metal joining
Lighting: Discharge lamps and plasma light sources
Space Re-entry: Plasma sheath on spacecraft

Key Physics:

Temperature range: 2,000 - 30,000 K
Pressure range: 0.1 - 100 bar (circuit breakers: 1-20 bar)
Current range: 1 - 100,000 A (industrial: 100-10,000 A)
Arc radius: 1 mm - 100 mm (typical: 5-20 mm)

3.2.2. Governing Equations

The arc plasma is modeled using energy and momentum balance in cylindrical coordinates (axisymmetric assumption).

3.2.2.1. Steady-State Energy Equation (Elenbaas-Heller)

\[ \frac{1}{r}\frac{d}{dr}\left(r \kappa(T) \frac{dT}{dr}\right) + \sigma(T) E^2 - S_{\text{rad}}(T) = 0 \]

Physical Terms:

Conduction: \(\frac{1}{r}\frac{d}{dr}(r\kappa\frac{dT}{dr})\) — Heat transport by thermal conductivity
Joule Heating: \(\sigma E^2\) — Ohmic dissipation from electric current
Radiation Loss: \(S_{\text{rad}} = 4\pi \varepsilon_{\text{nec}}(T)\) — Electromagnetic radiation emission

Electric Field Calculation:

Arc current constraint determines electric field:

\[ I = 2\pi \int_0^R \sigma(T) E \, r \, dr = 2\pi E G \]

where arc conductance \(G\) is:

\[ G = \int_0^R \sigma(T) r \, dr \]

Thus: \(E = I / (2\pi G)\)

3.2.2.2. Transient Energy Equation (Without Velocity)

\[ \rho C_p \frac{\partial T}{\partial t} = \frac{1}{r}\frac{\partial}{\partial r}\left(r\kappa\frac{\partial T}{\partial r}\right) + \sigma E^2 - S_{\text{rad}} \]

Additional Terms:

\(\rho(T)\): Mass density [kg/m³]
\(C_p(T)\): Specific heat at constant pressure [J/(kg·K)]
\(\frac{\partial T}{\partial t}\): Temporal temperature change

3.2.2.3. Transient Energy Equation (With Radial Velocity)

\[ \rho C_p \frac{\partial T}{\partial t} + \rho C_p V \frac{\partial T}{\partial r} = \frac{1}{r}\frac{\partial}{\partial r}\left(r\kappa\frac{\partial T}{\partial r}\right) + \sigma E^2 - S_{\text{rad}} \]

Continuity Equation (mass conservation):

\[ \frac{\partial (\rho r)}{\partial t} + \frac{\partial (\rho V r)}{\partial r} = 0 \]

Momentum-Derived Velocity:

\[ \frac{\partial V}{\partial r} = -\frac{1}{\rho^2 C_p}\frac{d\rho}{dT}\left[r(\sigma E^2 - S_{\text{rad}}) + \frac{\partial}{\partial r}\left(r\kappa\frac{\partial T}{\partial r}\right)\right] \]

Physical Interpretation:

Radial velocity driven by radial pressure gradients
Important in high-pressure systems (P > 10 bar)
Convective heat transport affects arc decay rate

3.2.3. Material Properties

Temperature-dependent plasma properties are essential for accurate modeling:

Transport Coefficients:

Electrical Conductivity \(\sigma(T)\): Increases sharply above ~5000 K
Thermal Conductivity \(\kappa(T)\): Complex behavior with reaction/ionization contributions
Viscosity \(\mu(T)\): Temperature-dependent (not always included)

Thermodynamic Properties:

Density \(\rho(T)\): Decreases with temperature (ideal gas approximation)
Specific Heat \(C_p(T)\): Peaks near dissociation/ionization temperatures
Enthalpy \(h(T)\): Energy content including chemical reactions

Radiation Properties:

Net Emission Coefficient (NEC) \(\varepsilon_{\text{nec}}(T, R)\): Depends on temperature AND arc radius
Accounts for self-absorption (optically thick vs. thin plasmas)
Critical for energy balance in arc core

Data Sources:

Experimental measurements (arc tunnels, shock tubes)
Chemical equilibrium calculations (NASA CEA, NIST databases)
Boltzmann equation analysis for electron properties

3.3. Examples

3.3.1. 1. Steady-State Arc Discharge

File: app/plasma/arc/solve_1d_arc_steady.py

Purpose: Solve the Elenbaas-Heller equation for steady-state thermal arc temperature distribution.

Physical Problem:

An electric arc at steady state where Joule heating balances conductive and radiative losses. The temperature profile establishes a quasi-equilibrium between energy input and dissipation mechanisms.

Typical Applications:

Characterizing arc plasma properties
Initial condition for transient simulations
Arc quenching ability assessment
Parametric studies (current, pressure, gas composition)

3.3.1.1. Physical Parameters

# Gas selection
gas = 'SF6'          # Options: 'SF6', 'Air', 'CO2', 'N2', 'Ar'

# Arc parameters
I = 200              # Arc current [A]
R = 10e-3           # Arc radius [m] (10 mm)
Tb = 2000           # Boundary temperature [K]

# Numerical parameters
mesh_num = 500      # Spatial mesh cells
relax = 0.1         # Relaxation factor (0 < relax ≤ 1)
converge_tol = 1e-6 # Convergence tolerance [K]
max_ite = 6000      # Maximum iterations

Parameter Guidelines:

Current: 10-10,000 A (typical circuit breakers: 100-1000 A)
Radius: 1-50 mm (determined by electrode geometry)
Boundary Temperature: 300-3000 K (depends on cooling)
Relaxation Factor: 0.05-0.3 (smaller = more stable, slower convergence)

3.3.1.2. Solution Method

Iterative Scheme:

Initialize: Linear or parabolic temperature profile
Update Properties: Interpolate \(\kappa(T), \sigma(T), \varepsilon(T)\) at current temperatures
Compute Conductance: \(G = \int_0^R \sigma(T) r \, dr\) (trapezoidal integration)
Calculate Electric Field: \(E = I / (2\pi G)\)
Solve Energy Equation: Finite volume method on 1D radial mesh
Under-Relaxation: \(T^{n+1} = T^n + \text{relax} \times (T^* - T^n)\)
Check Convergence: \(\text{RMS}(T^{n+1} - T^n) < \text{tol}\)
Repeat: Until convergence or max iterations

Finite Volume Discretization:

\[ \frac{r_{i+1/2}\kappa_{i+1/2}(T_{i+1} - T_i)}{\Delta r} - \frac{r_{i-1/2}\kappa_{i-1/2}(T_i - T_{i-1})}{\Delta r} + r_i \Delta r (\sigma_i E^2 - S_{\text{rad},i}) = 0 \]

Boundary Conditions:

At \(r = 0\) (axis): \(\frac{dT}{dr} = 0\) (symmetry)
At \(r = R\) (boundary): \(T = T_b\) (Dirichlet)

3.3.1.3. Data Files Required

Thermodynamic Properties (gas_p1.dat):

Format: Temperature, Density, Enthalpy, Cp, Electrical Conductivity, Thermal Conductivity

# T(K)    rho(kg/m³)  h(J/kg)      Cp(J/kg/K)  sigma(S/m)   kappa(W/m/K)
300       5.963       0.0          520.0       0.0          0.0134
500       3.578       104000       610.0       0.0          0.0220
1000      1.789       415000       890.0       0.001        0.0450
...
30000     0.060       45600000     8200.0      18000.0      6.5000

Net Emission Coefficient (gas_p1_nec.dat):

Format: Temperature × Radius matrix of NEC values

# First row: Radius values [m]
# First column: Temperature values [K]
# Matrix: NEC(T, R) [W/m³]

3.3.1.4. Run Example

python app/plasma/arc/solve_1d_arc_steady.py

3.3.1.5. Expected Output

Console Output:

======================================================================
1D Stationary Arc Model - SF6 Plasma
======================================================================

Loading plasma property data for SF6...
  - Thermodynamic data loaded: 150 temperature points
  - Temperature range: 300 - 30000 K
  - NEC data loaded: 150 temperatures × 50 radii
  - Radius range: 0.500 - 50.000 mm

Interpolating NEC for arc radius R = 10.0 mm...
  - NEC interpolated for 150 temperature points

Initializing arc model...
  - Arc current: 200 A
  - Arc radius: 10.0 mm
  - Mesh cells: 500
  - Boundary temperature: 2000 K

Solving 1D stationary arc equation...
  - Relaxation factor: 0.1
  - Convergence tolerance: 1.0e-06 K
  - Maximum iterations: 6000
----------------------------------------------------------------------
Iteration    0: RMS error = 1.234e+03 K, Max T = 12000.0 K
Iteration  100: RMS error = 5.678e+01 K, Max T = 21543.2 K
Iteration  200: RMS error = 2.345e+00 K, Max T = 22987.6 K
...
Iteration  850: RMS error = 8.765e-07 K, Max T = 23456.8 K
----------------------------------------------------------------------

Solution obtained successfully!

Temperature profile statistics:
  - Maximum temperature: 23456.78 K (at axis)
  - Minimum temperature: 2000.00 K
  - Boundary temperature: 2000.00 K

Electrical properties:
  - Arc conductance: 0.2341 S
  - Electric field: 136.4 V/m
  - Voltage drop: 1.364 V (over 1 cm)
  - Power dissipation: 272.8 W

Output Files:

results/sta_arc1d_SF6.csv: Temperature and property distributions
results/sta_arc1d_SF6.png: Temperature profile plot
results/sta_arc1d_SF6_properties.png: Conductivity, thermal conductivity plots

Typical Results:

Peak temperature: 20,000-25,000 K (for SF₆ at 200 A)
Temperature drops sharply near boundary (boundary layer ~0.5-2 mm)
Iteration count: 500-2000 (depends on relaxation factor)
Computation time: 5-30 seconds

3.3.1.6. Physical Insights

Temperature Distribution:

Maximum at arc center (r = 0)
Steep gradient near boundary (high heat loss zone)
Nearly flat in core (energy generation dominates)

Energy Balance:

Core region: Joule heating ≈ Radiation loss
Boundary region: Conduction dominates
Total power: \(P = E \times I\) (must match integrated source terms)

Scaling Laws:

Higher current → Higher peak temperature
Larger radius → Higher peak temperature (reduced heat loss/volume)
Higher pressure → Modified properties (NEC increases, σ varies)

3.3.2. 2. Transient Arc Discharge with Radial Velocity (Explicit Method)

File: app/plasma/arc/solve_1d_arc_transient_explicit.py

Purpose: Simulate time-dependent arc behavior including convective heat transport from radial gas flow.

Physical Scenario:

After current interruption in a circuit breaker, the arc decays as thermal energy dissipates through conduction, radiation, and convection. Radial velocity develops due to pressure gradients from non-uniform heating.

When to Use This Method:

High-pressure arcs (P > 10 bar)
Convection-dominated decay
Accurate blast flow modeling
Arc-chamber interaction studies

Key Feature: Explicit time integration is RECOMMENDED for better stability compared to implicit methods.

3.3.2.1. Physical Parameters

# Gas and geometry
gas = 'SF6'
I = 0                # Current = 0 for pure decay
R = 10e-3           # Arc radius [m]

# Time integration
dt = 1e-9           # Time step [s] (1 nanosecond)
step_num = 100000   # Number of steps (100 μs total)
save_freq = 100     # Save every 100 steps

# Numerical
mesh_num = 500      # Spatial mesh
Tb = 2000           # Boundary temperature [K]
enable_joule = False # Joule heating (False for decay)

Time Step Selection:

Stability Criterion: \(\Delta t < \frac{(\Delta r)^2}{2\alpha}\) where \(\alpha = \kappa/(\rho C_p)\) is thermal diffusivity
Typical: \(\Delta t = 10^{-9}\) to \(10^{-7}\) s
Smaller steps for steep gradients or high conductivity

3.3.2.2. Solution Method

Explicit Euler Time Stepping:

\[ T_i^{n+1} = T_i^n + \Delta t \left[\frac{1}{\rho C_p}\left(\nabla \cdot (\kappa \nabla T)\right)_i + \frac{\sigma E^2 - S_{\text{rad}}}{\rho C_p}\right]^n \]

Velocity Update:

Integrate velocity equation from continuity and momentum balance:

\[ V_{i+1/2}^{n+1} = V_{i-1/2}^{n+1} + \Delta r \left(\frac{\partial V}{\partial r}\right)_i^{n+1} \]

Algorithm:

Initialize: Load steady-state solution as \(T^0\)
Time Loop: For each time step \(n = 0, 1, ..., N-1\):
- Update properties: \(\kappa^n, \sigma^n, \rho^n, C_p^n, \varepsilon^n\) from \(T^n\)
- Compute conductance: \(G^n\)
- Calculate electric field: \(E^n = I / (2\pi G^n)\)
- Compute temperature flux: \(q_i^n = -\kappa_i \frac{T_{i+1} - T_i}{\Delta r}\)
- Compute velocity gradient: \(\frac{\partial V}{\partial r}\) from momentum equation
- Integrate velocity: \(V^{n+1}\) from axis outward
- Update temperature: \(T^{n+1}\) with convection term \(V \frac{\partial T}{\partial r}\)
- Apply boundary conditions
- Save if \(n \mod \text{save\_freq} = 0\)

Advantages of Explicit Method:

Simple implementation
Guaranteed stability with proper time step
No matrix inversion required
Easy to parallelize

3.3.2.3. Initial Condition

Load from steady-state solution:

initial_file = './app/plasma/arc/results/sta_arc1d_SF6.csv'
dat = pd.read_csv(initial_file)
x_list = dat['R(m)'].values
T_list = dat['T(K)'].values
Tfunc_init = interpolate.interp1d(x_list, T_list, kind='cubic')

Or define analytically (parabolic profile often reasonable):

T_center = 15000  # K
Tfunc_init = lambda r: (1 - (r/R)**2) * (T_center - Tb) + Tb

3.3.2.4. Run Example

python app/plasma/arc/solve_1d_arc_transient_explicit.py

3.3.2.5. Expected Output

Console Output:

======================================================================
1D Transient Arc Model WITH Radial Velocity (Explicit Method) - SF6
======================================================================

Loading plasma property data for SF6...
  - Thermodynamic data loaded: 150 temperature points
  - NEC data loaded: 150 temperatures × 50 radii

Loading initial temperature profile...
  - Reading from file: ./app/plasma/arc/results/sta_arc1d_SF6.csv
  - Initial profile loaded: T_max = 23456.78 K

Initializing transient arc model with radial velocity...
  - Arc current: 0 A (Decay mode - no Joule heating)
  - Arc radius: 10.0 mm
  - Solution method: EXPLICIT (recommended for stability)

Time integration parameters:
  - Time step: 1.00e-09 s
  - Number of steps: 100000
  - Total simulation time: 1.00e-04 s (100.00 μs)
  - Save frequency: every 100 step(s)
----------------------------------------------------------------------
Time step      0: t = 0.000e+00 s, T_max = 23456.8 K, V_max = 0.0 m/s
Time step   1000: t = 1.000e-06 s, T_max = 22873.5 K, V_max = 45.3 m/s
Time step   5000: t = 5.000e-06 s, T_max = 18234.7 K, V_max = 123.7 m/s
Time step  10000: t = 1.000e-05 s, T_max = 13456.2 K, V_max = 201.5 m/s
...
Time step 100000: t = 1.000e-04 s, T_max = 4532.8 K, V_max = 87.3 m/s
----------------------------------------------------------------------

Simulation completed successfully!
  - Final maximum temperature: 4532.8 K
  - Final maximum velocity: 87.3 m/s
  - Total computation time: 3.45 minutes

Output Files:

results/tra_arc1d_SF6_explicit.mat: MATLAB format with T(r,t), V(r,t)
results/tra_arc1d_SF6_explicit_*.png: Temperature snapshots
results/tra_arc1d_SF6_explicit.gif: Animated evolution

Typical Results:

Temperature decay: 23,000 K → 4,000 K in 100 μs
Peak velocity: 50-250 m/s (subsonic flow)
Velocity develops quickly (first ~10 μs)
Computation time: 2-10 minutes (depends on step count)

3.3.2.6. Physical Insights

Decay Phases:

Early (0-10 μs): Radiation-dominated cooling, velocity buildup
Middle (10-50 μs): Convection enhances cooling, velocity peaks
Late (>50 μs): Diffusion-dominated, velocity decreases

Velocity Distribution:

Maximum near arc boundary (strongest gradients)
Near-zero at axis (symmetry)
Outward flow (positive radial velocity)

Convection Effects:

Accelerates cooling by ~20-50% compared to no-velocity case
More important at high pressure
Reduces arc lifetime (beneficial for circuit breakers)

3.3.3. 3. Transient Arc Discharge without Radial Velocity

File: app/plasma/arc/solve_1d_arc_transient_noV.py

Purpose: Simulate arc decay considering only thermal diffusion (no convection).

Physical Justification:

Valid for low-pressure arcs (P < 5 bar)
Early decay phases where flow hasn’t developed
Simplified analysis and faster computation
Conservative estimate (slower cooling than with convection)

When to Use:

Low-pressure systems
Initial arc development (before flow)
Rapid parameter studies
Validation and benchmarking

3.3.3.1. Physical Parameters

# Gas and geometry
gas = 'SF6'
I = 0                # Current = 0 for decay
R = 10e-3           # Arc radius [m]

# Time integration
dt = 1e-6           # Time step [s] (1 microsecond - larger than with velocity)
step_num = 1000     # Number of steps (1 ms total)
save_freq = 1       # Save every step

# Numerical
mesh_num = 500               # Spatial mesh
sweep_max_num = 10           # Sweep iterations per step
sweep_res_tol = 1e-6        # Sweep convergence [K]
Tb = 2000                    # Boundary temperature [K]
enable_joule = False         # Joule heating

Time Step Selection:

Larger steps possible without velocity stiffness
Typical: \(\Delta t = 10^{-6}\) to \(10^{-5}\) s
Implicit method allows larger steps than explicit

3.3.3.2. Solution Method

Implicit Time Integration with Sweeping:

Crank-Nicolson or backward Euler scheme:

\[ \rho C_p \frac{T_i^{n+1} - T_i^n}{\Delta t} = \frac{1}{r}\frac{\partial}{\partial r}\left(r\kappa^{n+\theta}\frac{\partial T^{n+\theta}}{\partial r}\right) + \sigma^{n+\theta} E^2 - S_{\text{rad}}^{n+\theta} \]

where \(\theta \in [0, 1]\) controls implicitness (0 = explicit, 1 = fully implicit, 0.5 = Crank-Nicolson).

Sweep Iterations:

Nonlinear properties require iteration within each time step:

Guess: \(T^{n+1,0} = T^n\)
Sweep Loop: For \(k = 0, 1, ..., k_{\max}\):
- Update properties at \(T^{n+1,k}\)
- Solve linear system for \(T^{n+1,k+1}\)
- Check convergence: \(\text{RMS}(T^{n+1,k+1} - T^{n+1,k}) < \text{tol}\)
- If converged, proceed to next time step

Algorithm:

Initialize: Load steady-state or analytical profile
Time Loop: For \(n = 0, 1, ..., N-1\):
- Sweep Loop: Until convergence:
  - Update properties from current temperature guess
  - Assemble coefficient matrix (using FiPy)
  - Solve linear system: \(A \mathbf{T}^{n+1} = \mathbf{b}\)
  - Check sweep convergence
- Save if needed
Output: Temperature evolution data

3.3.3.3. Run Example

python app/plasma/arc/solve_1d_arc_transient_noV.py

3.3.3.4. Expected Output

Console Output:

======================================================================
1D Transient Arc Model (No Radial Velocity) - SF6 Plasma
======================================================================

Loading plasma property data for SF6...
  - Thermodynamic data loaded: 150 temperature points
  - NEC data loaded: 150 temperatures × 50 radii

Loading initial temperature profile...
  - Reading from file: ./app/plasma/arc/results/sta_arc1d_SF6.csv
  - Initial profile loaded: T_max = 23456.78 K

Initializing transient arc model...
  - Arc current: 0 A (Decay mode - no Joule heating)
  - Arc radius: 10.0 mm
  - Mesh cells: 500

Time integration parameters:
  - Time step: 1.00e-06 s
  - Number of steps: 1000
  - Total simulation time: 1.00e-03 s (1000.00 μs)
  - Joule heating: Disabled
----------------------------------------------------------------------
Step    0: t = 0.000e+00 s, T_max = 23456.8 K, Sweeps = 0
Step   10: t = 1.000e-05 s, T_max = 21234.5 K, Sweeps = 3
Step   50: t = 5.000e-05 s, T_max = 16789.3 K, Sweeps = 4
Step  100: t = 1.000e-04 s, T_max = 13567.2 K, Sweeps = 4
...
Step 1000: t = 1.000e-03 s, T_max = 3245.7 K, Sweeps = 2
----------------------------------------------------------------------

Simulation completed successfully!
  - Final maximum temperature: 3245.7 K
  - Average sweeps per step: 3.5
  - Total computation time: 45.2 seconds

Output Files:

results/tra_arc1d_noV_SF6.mat: Temperature evolution T(r,t)
results/tra_arc1d_noV_SF6_*.png: Snapshots at selected times
results/tra_arc1d_noV_SF6_comparison.png: With/without Joule heating

Typical Results:

Slower decay than with velocity
Temperature: 23,000 K → 3,000 K in 1 ms (vs. ~200 μs with velocity)
Sweep iterations: 2-5 per time step
Computation time: 30-60 seconds

3.3.3.5. Comparison: With vs. Without Velocity

Feature	With Velocity	Without Velocity
Cooling Rate	Faster (~5×)	Slower
Time Step	Small (ns)	Larger (μs)
Computational Cost	Higher	Lower
Physical Realism	High (P > 10 bar)	Moderate
Use Case	High-pressure, accurate	Low-pressure, fast

3.4. Best Practices

3.4.1. 1. Material Property Data

Data Quality:

Use validated property tables from literature (NIST, NASA, experiments)
Ensure smooth interpolation (avoid oscillations)
Check physical consistency (positive σ, κ, ρ, Cp)
Verify units (SI standard)

Temperature Range:

Extend beyond expected arc temperatures
Include low-temperature region (300-2000 K near boundary)
High-temperature region (up to 30,000 K for arc core)

Pressure Effects:

NEC strongly depends on pressure (via arc radius)
Higher pressure → Higher NEC → More radiation
Use pressure-appropriate property tables

3.4.2. 2. Mesh and Discretization

Spatial Resolution:

Minimum: 200 cells (coarse, fast)
Recommended: 500-1000 cells (good accuracy)
High precision: 2000+ cells (research)

Mesh Refinement:

Concentrate points near boundary (steep gradients)
Uniform mesh acceptable for cylindrical geometry
Check grid independence (repeat with 2× cells)

Boundary Layer:

Ensure ~10-20 cells in temperature drop region
Critical for accurate heat flux calculation

3.4.3. 3. Time Integration

Stability:

Explicit: \(\Delta t < \frac{(\Delta r)^2}{2\alpha_{\max}}\) where \(\alpha = \kappa/(\rho C_p)\)
Implicit: Larger steps allowed, but accuracy may suffer
Use explicit for velocity, implicit for pure diffusion

Accuracy:

Start with small steps, gradually increase if stable
Monitor solution smoothness and energy conservation
Typical time scales:
- Arc ignition: 1-100 ns
- Current zero: 1-10 μs
- Post-arc decay: 10-1000 μs

3.4.4. 4. Convergence and Validation

Steady-State Convergence:

Monitor RMS error (should decrease exponentially)
Check maximum temperature (should stabilize)
Verify energy balance: \(\int (σE^2 - S_{\text{rad}}) dV = 0\)

Transient Validation:

Compare with experiments (if available)
Energy conservation: Track total energy over time
Physical consistency: Temperatures should be monotonic

Parameter Studies:

Sweep current: 50 A to 5000 A
Vary radius: 5 mm to 50 mm
Different gases: SF₆, Air, CO₂, N₂

3.4.5. 5. Output and Visualization

Essential Plots:

Temperature Profile: T vs. r at multiple times
Velocity Profile: V vs. r (if applicable)
Property Evolution: σ, κ, NEC vs. r
Time Series: T_max vs. t, Energy vs. t

Animation:

Create GIF showing temperature evolution
Useful for presentations and debugging
Export frames at regular intervals

Data Export:

CSV for plotting and analysis
MATLAB .mat for post-processing
HDF5 for large datasets

3.5. Performance Benchmarks

Typical computational performance (standard PC):

Simulation Type	Mesh Size	Time Steps	Time/Step	Total Time	Memory
Steady-State	500	~1000 iter	0.02 s	20 s	<100 MB
Transient (No V)	500	1000	0.05 s	50 s	<200 MB
Transient (With V, Explicit)	500	100,000	0.002 s	200 s	<500 MB
Transient (With V, Implicit)	500	10,000	0.05 s	500 s	<300 MB

Hardware: Intel i7 CPU, 16 GB RAM (no GPU acceleration)

Optimization Tips:

Use compiled property interpolation (NumPy vectorization)
Reduce save frequency for long simulations
Pre-compute property tables on finer grid
Consider Numba or Cython for hotspots

3.6. Troubleshooting Guide

Problem	Possible Causes	Solutions
Steady-state not converging	Poor initial guess, large relaxation	Reduce relax to 0.05-0.1, better T_init
Negative temperatures	Time step too large, instability	Reduce dt, check boundary conditions
Oscillatory solution	Insufficient mesh resolution	Increase mesh_num, check property smoothness
Very slow convergence	Temperature-dependent stiffness	Use implicit method, adaptive time stepping
Unphysical velocities	Pressure gradient singularities	Smooth temperature profile, check ρ(T)
High memory usage	Saving too frequently	Increase save_freq, reduce step_num

3.7. Advanced Topics

3.7.1. Arc Quenching Ability

Definition: Measure of gas effectiveness at extinguishing arcs (crucial for circuit breakers).

Metrics:

Critical Current: Minimum current for arc sustenance
Decay Time Constant: Time for temperature to drop to threshold
Arc Voltage: Higher voltage = better quenching

Calculation Workflow:

Run steady-state for various currents (10-500 A)
Compute arc voltage: \(V = E \times L\) where \(L\) is arc length
Plot V-I curve (arc characteristic)
Run transient decay from steady-state
Extract time constant: \(\tau = -t / \ln(T/T_0)\)

References:

IEEE Std 242-2001 (Buff Book)
L. Zhong et al., IEEE Trans. Plasma Sci., 2019

3.7.2. Multi-Component Mixtures

Approach: Interpolate between pure gas properties

For gas mixture with fractions \(x_i\):

\[ \kappa_{\text{mix}} = \sum_i x_i \kappa_i(T) \]

(more sophisticated: Wilke’s formula for transport properties)

Implementation:

Load property tables for each gas component
Compute weighted averages at each temperature
Create composite property table
Run simulation with mixed properties

3.7.3. Radiation Models

Optically Thin: NEC independent of radius (valid for small arcs)

Optically Thick: NEC depends on R (self-absorption)

Improved Model: Solve radiative transfer equation

\[ \frac{dI_\nu}{ds} = \kappa_\nu (B_\nu - I_\nu) \]

(typically done offline, tabulated as NEC(T, R))

3.8. Applications

3.8.1. Circuit Breaker Design

Goal: Optimize gas composition and chamber geometry for fast arc quenching

Workflow:

Define operating conditions (current, pressure)
Run steady-state for arc characteristics
Simulate decay after current zero
Evaluate quenching time
Compare gases/mixtures

Deliverables:

Arc temperature distribution
Cooling time constants
Arc voltage characteristics
Gas selection recommendations

3.8.2. Welding Process Optimization

Goal: Maintain stable arc for consistent weld quality

Parameters:

Current: 50-500 A
Arc length: 2-10 mm
Shielding gas: Ar, Ar-CO₂, He

Simulation Outputs:

Arc temperature (affects melting)
Arc pressure (affects penetration)
Heat flux to workpiece

3.8.3. Plasma Torch Design

Goal: Achieve desired gas temperature and flow for material processing

Features:

High current: 100-1000 A
Confined arc (small radius)
Gas injection (modeled as velocity BC)

Outputs:

Exit gas temperature
Thermal efficiency
Electrode erosion estimates

3.9. References

[1] L. Zhong, Y. Cressault, and P. Teulet, “Evaluation of Arc Quenching Ability for a Gas by Combining 1-D Hydrokinetic Modeling and Boltzmann Equation Analysis,” IEEE Trans. Plasma Sci., vol. 47, no. 4, pp. 1835-1840, 2019.

[2] L. Zhong, Q. Gu, and S. Zheng, “An improved method for fast evaluating arc quenching performance of a gas based on 1D arc decaying model,” Physics of Plasmas, vol. 26, no. 10, p. 103507, 2019.

3.10. Quick Start Guide

New Users — Start with these examples:

solve_1d_arc_steady.py — Learn basic steady-state arc physics
solve_1d_arc_transient_noV.py — Understand arc decay
solve_1d_arc_transient_explicit.py — Full physics with convection

For Circuit Breaker Engineers:

Focus on transient decay simulations
Compare different SF₆ alternatives (eco-friendly gases)
Evaluate quenching metrics

For Plasma Physicists:

Use for validation of ML methods (PINN, DeepONet)
Generate training data for neural operators
Explore multi-scale phenomena

For Method Developers:

Benchmark against analytical solutions
Test numerical schemes
Develop adaptive algorithms