# 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.

## 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

### Why Classical Plasma Simulations?

Classical simulations serve multiple purposes in AI4Plasma:

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

### 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)

---

## Arc Discharge Simulations

### 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)

### Governing Equations

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

#### 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)$

#### 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

#### 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

### 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

---

## Examples

### 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)

#### Physical Parameters

```python
# 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)

#### Solution Method

**Iterative Scheme**:

1. **Initialize**: Linear or parabolic temperature profile
2. **Update Properties**: Interpolate $\kappa(T), \sigma(T), \varepsilon(T)$ at current temperatures
3. **Compute Conductance**: $G = \int_0^R \sigma(T) r \, dr$ (trapezoidal integration)
4. **Calculate Electric Field**: $E = I / (2\pi G)$
5. **Solve Energy Equation**: Finite volume method on 1D radial mesh
6. **Under-Relaxation**: $T^{n+1} = T^n + \text{relax} \times (T^* - T^n)$
7. **Check Convergence**: $\text{RMS}(T^{n+1} - T^n) < \text{tol}$
8. **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)

#### 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³]
```

#### Run Example

```bash
python app/plasma/arc/solve_1d_arc_steady.py
```

#### 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

#### 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)

---

### 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.

#### Physical Parameters

```python
# 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

#### 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**:

1. **Initialize**: Load steady-state solution as $T^0$
2. **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

#### Initial Condition

Load from steady-state solution:

```python
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):

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

#### Run Example

```bash
python app/plasma/arc/solve_1d_arc_transient_explicit.py
```

#### 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)

#### Physical Insights

**Decay Phases**:
1. **Early (0-10 μs)**: Radiation-dominated cooling, velocity buildup
2. **Middle (10-50 μs)**: Convection enhances cooling, velocity peaks
3. **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. 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

#### Physical Parameters

```python
# 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

#### 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:

1. **Guess**: $T^{n+1,0} = T^n$
2. **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**:

1. **Initialize**: Load steady-state or analytical profile
2. **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
3. **Output**: Temperature evolution data

#### Run Example

```bash
python app/plasma/arc/solve_1d_arc_transient_noV.py
```

#### 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

#### 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 |

---

## Best Practices

### 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

### 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. 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

### 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₂

### 5. Output and Visualization

**Essential Plots**:
1. **Temperature Profile**: T vs. r at multiple times
2. **Velocity Profile**: V vs. r (if applicable)
3. **Property Evolution**: σ, κ, NEC vs. r
4. **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

---

## 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

---

## 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 |

---

## Advanced Topics

### Arc Quenching Ability

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

**Metrics**:
1. **Critical Current**: Minimum current for arc sustenance
2. **Decay Time Constant**: Time for temperature to drop to threshold
3. **Arc Voltage**: Higher voltage = better quenching

**Calculation Workflow**:
1. Run steady-state for various currents (10-500 A)
2. Compute arc voltage: $V = E \times L$ where $L$ is arc length
3. Plot V-I curve (arc characteristic)
4. Run transient decay from steady-state
5. 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

### 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**:
1. Load property tables for each gas component
2. Compute weighted averages at each temperature
3. Create composite property table
4. Run simulation with mixed properties

### 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))

---

## Applications

### Circuit Breaker Design

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

**Workflow**:
1. Define operating conditions (current, pressure)
2. Run steady-state for arc characteristics
3. Simulate decay after current zero
4. Evaluate quenching time
5. Compare gases/mixtures

**Deliverables**:
- Arc temperature distribution
- Cooling time constants
- Arc voltage characteristics
- Gas selection recommendations

### 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

### 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

---

## 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.

---

## Quick Start Guide

**New Users** — Start with these examples:

1. `solve_1d_arc_steady.py` — Learn basic steady-state arc physics
2. `solve_1d_arc_transient_noV.py` — Understand arc decay
3. `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