Arc Density Constraint Model

This example demonstrates the ArcDensityConstraintModel which enforces energy density limits in the plasma arc to prevent excessive thermal loading.

==================================================
Arc Density Constraint Model
==================================================
Energy Density: 3.18e+05 J/m³
Maximum Allowed: 1.00e+09 J/m³
Margin: 1.00e+09 J/m³
Constraint Satisfied: Yes
==================================================
[OK] Energy density is within safe limits

Generating plot: Energy Density vs Energy per Pulse...

import matplotlib.pyplot as plt
import numpy as np
import openmdao.api as om

from paroto.models.arc_density import ArcDensityConstraintModel

# Create OpenMDAO problem with max energy density parameter
prob = om.Problem()
prob.model.add_subsystem(
    "density",
    ArcDensityConstraintModel(model_parameters={"max_energy_density": 1e9}),  # 1 GJ/m³
)
prob.setup()

# Set parameters for typical operating conditions
prob.set_val("density.energy_per_pulse", 10e-3)  # 10 mJ
prob.set_val("density.thermal_diameter", 0.002)  # 2 mm arc diameter
prob.set_val("density.arc_length", 0.01)  # 10 mm arc length

# Run model
prob.run_model()

# Extract results
rho_E = prob.get_val("density.energy_density")[0]
margin = prob.get_val("density.density_margin")[0]
satisfied = prob.get_val("density.density_constraint_satisfied")[0]

print("=" * 50)
print("Arc Density Constraint Model")
print("=" * 50)
print(f"Energy Density: {rho_E:.2e} J/m³")
print("Maximum Allowed: 1.00e+09 J/m³")
print(f"Margin: {margin:.2e} J/m³")
print(f"Constraint Satisfied: {'Yes' if satisfied > 0.5 else 'No'}")
print("=" * 50)

if satisfied > 0.5:
    print("[OK] Energy density is within safe limits")
else:
    print("[FAIL] Energy density exceeds maximum threshold")
    print("  Reduce energy per pulse or increase arc volume")

# Plot energy density vs energy per pulse
print("\nGenerating plot: Energy Density vs Energy per Pulse...")
energy_range = np.linspace(1e-3, 50e-3, 50)  # 1 mJ to 50 mJ
density_values = []
max_density = 1e9  # J/m³

for energy in energy_range:
    prob.set_val("density.energy_per_pulse", energy)
    prob.run_model()
    density_values.append(prob.get_val("density.energy_density")[0])

plt.figure(figsize=(8, 5))
plt.plot(
    energy_range * 1e3, np.array(density_values) / 1e9, "b-", linewidth=2, label="Energy Density"
)
plt.axhline(y=1.0, color="r", linestyle="--", linewidth=2, label="Max Allowed (1 GJ/m³)")
plt.axvline(x=10, color="g", linestyle=":", alpha=0.7, label="Nominal (10 mJ)")
plt.xlabel("Energy per Pulse (mJ)")
plt.ylabel("Energy Density (GJ/m³)")
plt.title("Arc Energy Density vs Pulse Energy")
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()