7f040c3cf5
1. Iron saturation now in the dynamics (was the "13 Tesla" bug). Flux linkage lambda(x,I) = L_air*I + L_iron*overlap(x)*g(I) with g(I)=I_sat*tanh(I/I_sat) saturating at the current where iron reaches B_sat. Both the circuit back-EMF (dlambda/dx) and the force (coenergy, dW'/dx) derive from the SAME lambda, so energy stays conserved (0.025% error) AND the field caps at B_sat instead of running to 13 T. Reduces to the old 0.5*I^2*dL/dx in the low-current limit. On the config that reported 60% efficiency / 11.6 T, it now gives 40% / 1.90 T. 2. Switch surge-current limit is now a hard feasibility constraint: a config whose peak discharge current exceeds the switch's pulse rating (continuous * surge factor per device kind) is infeasible -- otherwise the optimizer "wins" with configs that vaporize their own thyristor (e.g. 119 A through a 12 A BT151). Capacitor current stays a warning (electrolytics tolerate pulses; DB has continuous not pulse ratings). Regression tests added for both (field cap near saturation, energy conservation under strong discharge). 87 tests pass. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
114 lines
4.1 KiB
Python
114 lines
4.1 KiB
Python
"""Проверка ОДУ разряда против аналитического решения простого RLC-контура.
|
||
|
||
dL/dx заморожена в ноль (mu_eff=1.0), поэтому механика отвязана от контура
|
||
и [Q, I] эволюционируют как классический последовательный RLC-разряд
|
||
конденсатора. Сравниваем численное решение solve_ivp с учебным аналитическим
|
||
решением во всех трёх режимах затухания.
|
||
"""
|
||
|
||
import math
|
||
|
||
import numpy as np
|
||
import pytest
|
||
from scipy.integrate import solve_ivp
|
||
|
||
from gausse.physics.circuit import Q, StageCircuitParams, derivatives
|
||
from gausse.physics.inductance import CoilInductanceModel
|
||
|
||
L_H = 100e-6
|
||
C_F = 200e-6
|
||
V0 = 400.0
|
||
Q0 = C_F * V0
|
||
|
||
|
||
def _make_frozen_model() -> CoilInductanceModel:
|
||
return CoilInductanceModel(
|
||
l_air_h=L_H, coil_length_m=0.05, slug_length_m=0.03, mu_eff=1.0, smoothing_width_m=0.001
|
||
)
|
||
|
||
|
||
def _integrate(r_ohm: float, t_span: tuple[float, float], t_eval: np.ndarray):
|
||
# mu_eff=1.0 -> вклад железа нулевой, λ = L_air·I, чистый RLC-контур
|
||
model = _make_frozen_model()
|
||
params = StageCircuitParams(
|
||
capacitance_f=C_F,
|
||
r_total_ohm=r_ohm,
|
||
mass_kg=1.0,
|
||
)
|
||
state0 = [Q0, 0.0, 0.0, 0.0]
|
||
return solve_ivp(
|
||
derivatives,
|
||
t_span,
|
||
state0,
|
||
args=(model, params),
|
||
t_eval=t_eval,
|
||
method="RK45",
|
||
rtol=1e-10,
|
||
atol=1e-14,
|
||
)
|
||
|
||
|
||
def test_underdamped_matches_analytical_solution():
|
||
r_ohm = 0.05
|
||
alpha = r_ohm / (2 * L_H)
|
||
omega0 = 1 / math.sqrt(L_H * C_F)
|
||
omega_d = math.sqrt(omega0**2 - alpha**2)
|
||
|
||
t_eval = np.linspace(0, 5 / alpha, 200)
|
||
sol = _integrate(r_ohm, (t_eval[0], t_eval[-1]), t_eval)
|
||
|
||
expected_q = Q0 * np.exp(-alpha * t_eval) * (
|
||
np.cos(omega_d * t_eval) + (alpha / omega_d) * np.sin(omega_d * t_eval)
|
||
)
|
||
expected_i = Q0 * (omega0**2 / omega_d) * np.exp(-alpha * t_eval) * np.sin(omega_d * t_eval)
|
||
|
||
assert np.allclose(sol.y[Q], expected_q, rtol=1e-3, atol=Q0 * 1e-4)
|
||
assert np.allclose(sol.y[1], expected_i, rtol=1e-3, atol=abs(expected_i).max() * 1e-3)
|
||
|
||
|
||
def test_overdamped_matches_analytical_solution():
|
||
r_ohm = 5.0
|
||
alpha = r_ohm / (2 * L_H)
|
||
omega0 = 1 / math.sqrt(L_H * C_F)
|
||
root = math.sqrt(alpha**2 - omega0**2)
|
||
s1, s2 = -alpha + root, -alpha - root
|
||
|
||
t_eval = np.linspace(0, 5 / abs(s1), 200)
|
||
sol = _integrate(r_ohm, (t_eval[0], t_eval[-1]), t_eval)
|
||
|
||
a_coef = -Q0 * s2 / (s1 - s2)
|
||
b_coef = Q0 * s1 / (s1 - s2)
|
||
expected_q = a_coef * np.exp(s1 * t_eval) + b_coef * np.exp(s2 * t_eval)
|
||
expected_i = -(a_coef * s1 * np.exp(s1 * t_eval) + b_coef * s2 * np.exp(s2 * t_eval))
|
||
|
||
assert np.allclose(sol.y[Q], expected_q, rtol=1e-3, atol=Q0 * 1e-4)
|
||
assert np.allclose(sol.y[1], expected_i, rtol=1e-3, atol=max(abs(expected_i).max(), 1e-9) * 1e-3)
|
||
|
||
|
||
def test_critically_damped_matches_analytical_solution():
|
||
alpha = 1 / math.sqrt(L_H * C_F) # R = 2*sqrt(L/C) => alpha = omega0
|
||
r_ohm = 2 * alpha * L_H
|
||
|
||
t_eval = np.linspace(0, 5 / alpha, 200)
|
||
sol = _integrate(r_ohm, (t_eval[0], t_eval[-1]), t_eval)
|
||
|
||
expected_q = Q0 * (1 + alpha * t_eval) * np.exp(-alpha * t_eval)
|
||
expected_i = alpha**2 * Q0 * t_eval * np.exp(-alpha * t_eval)
|
||
|
||
assert np.allclose(sol.y[Q], expected_q, rtol=1e-3, atol=Q0 * 1e-4)
|
||
assert np.allclose(sol.y[1], expected_i, rtol=1e-3, atol=max(abs(expected_i).max(), 1e-9) * 1e-3)
|
||
|
||
|
||
def test_energy_conserved_when_lossless():
|
||
"""С R=0 энергия конденсатора должна переходить в энергию катушки без потерь."""
|
||
r_ohm = 0.0
|
||
t_eval = np.linspace(0, 1e-3, 500)
|
||
sol = _integrate(r_ohm, (t_eval[0], t_eval[-1]), t_eval)
|
||
|
||
cap_energy = sol.y[Q] ** 2 / (2 * C_F)
|
||
inductor_energy = 0.5 * L_H * sol.y[1] ** 2
|
||
total_energy = cap_energy + inductor_energy
|
||
initial_energy = Q0**2 / (2 * C_F)
|
||
|
||
assert np.allclose(total_energy, initial_energy, rtol=1e-4)
|