2046dcba10
- components/schema.py + database.py: typed loader for real retail parts (wire, capacitors, switches, sensors, projectile materials) - physics/inductance.py: multilayer solenoid inductance (Wheeler) with ferromagnetic-slug coupling via smooth overlap model and analytic dL/dx - physics/force.py: F = 0.5*I^2*dL/dx with a saturation clamp - physics/circuit.py: coupled [Q,I,x,v] discharge ODE - Validated against textbook RLC analytical solutions (under/over/critically damped) and energy conservation, not just spot-checked by eye - PLAN.md updated with animation (per-run field visualization) and GPU-accelerated batch sweep as explicit later stages Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
117 lines
4.2 KiB
Python
117 lines
4.2 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):
|
||
model = _make_frozen_model()
|
||
params = StageCircuitParams(
|
||
capacitance_f=C_F,
|
||
r_total_ohm=r_ohm,
|
||
mass_kg=1.0,
|
||
mu_eff=1.0,
|
||
total_turns=100,
|
||
coil_length_m=0.05,
|
||
b_sat_tesla=1e9, # эффективно отключает клэмп насыщения для этого теста
|
||
)
|
||
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)
|