Add component database layer and validated physics core

- 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>
This commit is contained in:
jze9
2026-07-06 19:45:32 +05:00
parent 04233aefd7
commit 2046dcba10
17 changed files with 698 additions and 0 deletions

1
.gitignore vendored
View File

@@ -11,3 +11,4 @@ dist/
*.sqlite3-shm
out/
.DS_Store
.claude/settings.local.json

View File

@@ -31,4 +31,6 @@ ChipDip, Cable.ru и др.).
- [ ] **Этап 5 — Хранилище результатов**: `storage/schema.py` + `database.py` — SQLite (WAL), таблица `runs` со всеми прогонами (успех/провал + честная причина), однопроцессный писатель поверх многопроцессной очереди.
- [ ] **Этап 6 — Поиск и оптимизация**: `optim/search_space.py` (геном переменной длины), `objective.py` (КПД), дешёвый квазистатический предфильтр, `optim/sweep.py` (Monte Carlo/LHS, миллионы прогонов), `optim/evolutionary.py` (ГА + coordinate-descent полировка) — всё пишется в общую таблицу `runs`.
- [ ] **Этап 7 — Отчётность**: `report/plots.py`, `summary.py`, `bom.py`, раздел "ограничения модели", `cli.py` (`gausse sweep/evolve/simulate/report`).
- [ ] `report/animate.py` — анимация одного прогона: положение снаряда в трубе по времени + визуализация поля/тока каждой катушки (свечение/интенсивность цвета ~ ток), сохранение в GIF/MP4 (matplotlib `FuncAnimation`). Нужна по запросу пользователя — "графика где будет показана симуляция пролёта цилиндра по трубе и электромагнитные поля в каждый момент времени".
- [ ] **Этап 8 — Сквозная проверка**: резюмируемый прогон на реальной базе компонентов, проверка честной записи в SQLite, финальный отчёт (КПД, скорость, стоимость, сравнение датчиков) с разделом ограничений.
- [ ] **Этап 9 — GPU-ускорение массового sweep (сервер с GTX 1070)**: после того как CPU/`scipy.solve_ivp`-модель провалидирована тестами (Этап 2-4) — батч-версия интегратора с ФИКСИРОВАННЫМ шагом (RK4/полу-неявная схема), считающая сразу N траекторий параллельно как один тензор (`cupy`, если доступна CUDA, иначе векторизованный `numpy`/`numba`), для прогона по-настоящему миллионов конфигураций на сервере. Важно: сначала корректность на CPU, потом скорость на GPU — численные результаты GPU-пути должны сверяться с CPU-эталоном на контрольной выборке, чтобы ускорение не подменило точность честными числами "для галочки".

View File

@@ -0,0 +1,61 @@
"""Загрузка базы реальных компонентов из JSON и фильтры для поиска."""
import json
from dataclasses import dataclass, field
from pathlib import Path
from gausse.components.schema import (
CapacitorSpec,
ProjectileMaterialSpec,
SensorSpec,
SwitchSpec,
WireSpec,
)
DEFAULT_DATA_DIR = Path(__file__).parent / "data"
def _load_json(path: Path) -> list[dict]:
if not path.exists():
return []
return json.loads(path.read_text(encoding="utf-8"))
@dataclass
class ComponentDatabase:
wires: list[WireSpec] = field(default_factory=list)
capacitors: list[CapacitorSpec] = field(default_factory=list)
switches: list[SwitchSpec] = field(default_factory=list)
sensors: list[SensorSpec] = field(default_factory=list)
projectile_materials: list[ProjectileMaterialSpec] = field(default_factory=list)
@classmethod
def load(cls, data_dir: Path = DEFAULT_DATA_DIR) -> "ComponentDatabase":
return cls(
wires=[WireSpec(**row) for row in _load_json(data_dir / "wires.json")],
capacitors=[
CapacitorSpec(**row) for row in _load_json(data_dir / "capacitors.json")
],
switches=[SwitchSpec(**row) for row in _load_json(data_dir / "switches.json")],
sensors=[SensorSpec(**row) for row in _load_json(data_dir / "sensors.json")],
projectile_materials=[
ProjectileMaterialSpec(**row)
for row in _load_json(data_dir / "projectile_materials.json")
],
)
def wires_by_material(self, material: str) -> list[WireSpec]:
return [w for w in self.wires if w.material == material]
def capacitors_rated_at_least(self, voltage_v: float) -> list[CapacitorSpec]:
return [c for c in self.capacitors if c.voltage_v >= voltage_v]
def switches_rated_for(self, current_a: float, voltage_v: float) -> list[SwitchSpec]:
return [
s
for s in self.switches
if s.max_current_a >= current_a and s.max_voltage_v >= voltage_v
]
def sensors_by_kind(self, kind: str) -> list[SensorSpec]:
return [s for s in self.sensors if s.kind == kind]

View File

@@ -0,0 +1,61 @@
"""Типизированные описания реальных, продающихся в розницу компонентов."""
from dataclasses import dataclass
from typing import Literal
@dataclass(frozen=True)
class WireSpec:
part_id: str
material: Literal["copper", "aluminum"]
gauge_mm: float
insulation_od_mm: float
resistivity_ohm_m: float
max_current_a: float
price_per_m: float
source: str
@dataclass(frozen=True)
class CapacitorSpec:
part_number: str
capacitance_uf: float
voltage_v: float
esr_ohm: float
max_current_a: float
price: float
source: str
@dataclass(frozen=True)
class SwitchSpec:
part_number: str
kind: Literal["SCR", "MOSFET", "IGBT"]
max_current_a: float
max_voltage_v: float
on_resistance_ohm: float | None
on_voltage_drop_v: float | None
turn_on_time_ns: float
price: float
source: str
@dataclass(frozen=True)
class SensorSpec:
part_number: str
kind: Literal["inductive", "optical", "hall"]
propagation_delay_ns: float
price: float
source: str
sensitivity_v_per_mps: float | None = None
threshold_v: float | None = None
@dataclass(frozen=True)
class ProjectileMaterialSpec:
name: str
density_kg_m3: float
mu_r: float
b_sat_tesla: float
price_per_kg: float
source: str

View File

@@ -0,0 +1,66 @@
"""ОДУ разряда конденсатора через ключ и катушку, связанное с механикой снаряда.
Состояние: [Q, I, x, v] — заряд конденсатора, ток в катушке, положение и
скорость снаряда. Индуктивность зависит от x (проход снаряда), поэтому
уравнение на dI/dt выведено из закона сохранения потокосцепления
Ψ = L(x)*I:
V_c = Q / C
dQ/dt = -I
dI/dt = (V_c - I*R_total - I*(dL/dx)*v) / L(x)
dv/dt = F(x, I) / m
dx/dt = v
"""
from dataclasses import dataclass
import numpy as np
from gausse.physics.force import force_on_slug_newtons
from gausse.physics.inductance import CoilInductanceModel
Q, I, X, V = range(4)
@dataclass(frozen=True)
class StageCircuitParams:
capacitance_f: float
r_total_ohm: float
mass_kg: float
mu_eff: float
total_turns: int
coil_length_m: float
b_sat_tesla: float
def derivatives(
t: float,
state: np.ndarray,
inductance_model: CoilInductanceModel,
params: StageCircuitParams,
) -> list[float]:
q, current, x, v = state
l_x = inductance_model.l_of_x(x)
dl_dx = inductance_model.dl_dx(x)
v_c = q / params.capacitance_f
d_q = -current
d_i = (v_c - current * params.r_total_ohm - current * dl_dx * v) / l_x
force = force_on_slug_newtons(
current_a=current,
dl_dx_unsaturated=dl_dx,
mu_eff=params.mu_eff,
total_turns=params.total_turns,
coil_length_m=params.coil_length_m,
b_sat_tesla=params.b_sat_tesla,
)
d_v = force / params.mass_kg
return [d_q, d_i, v, d_v]
def zero_current_crossing_event(t: float, state: np.ndarray, *_args) -> float:
return state[I]
zero_current_crossing_event.terminal = True
zero_current_crossing_event.direction = -1.0

View File

@@ -0,0 +1,6 @@
"""Физические константы и справочные удельные сопротивления."""
MU_0 = 4e-7 * 3.141592653589793 # Гн/м, магнитная проницаемость вакуума
RESISTIVITY_COPPER_OHM_M = 1.68e-8 # при ~20°C
RESISTIVITY_ALUMINUM_OHM_M = 2.82e-8 # при ~20°C

View File

@@ -0,0 +1,36 @@
"""Сила, действующая на ферромагнитный снаряд: F = 0.5*I^2*dL/dx с насыщением.
Насыщение — грубая, явно приближённая поправка (v1): оцениваем поле внутри
катушки простой соленоидной формулой и, если оно превышает B_sat материала
снаряда, пропорционально ослабляем вклад сердечника в dL/dx. Реальное
насыщение нелинейно и зависит от геометрии — здесь используется линейный
клэмп как консервативная эвристика, а не точная кривая намагничивания.
"""
from gausse.physics.constants import MU_0
def solenoid_field_estimate_tesla(
mu_eff: float, total_turns: int, coil_length_m: float, current_a: float
) -> float:
turns_per_m = total_turns / coil_length_m
return MU_0 * mu_eff * turns_per_m * abs(current_a)
def saturation_scale(b_estimate_tesla: float, b_sat_tesla: float) -> float:
if b_estimate_tesla <= 0:
return 1.0
return min(1.0, b_sat_tesla / b_estimate_tesla)
def force_on_slug_newtons(
current_a: float,
dl_dx_unsaturated: float,
mu_eff: float,
total_turns: int,
coil_length_m: float,
b_sat_tesla: float,
) -> float:
b_estimate = solenoid_field_estimate_tesla(mu_eff, total_turns, coil_length_m, current_a)
scale = saturation_scale(b_estimate, b_sat_tesla)
return 0.5 * current_a**2 * dl_dx_unsaturated * scale

View File

@@ -0,0 +1,113 @@
"""Индуктивность многослойного соленоида с ферромагнитным сердечником.
Приближённая инженерная модель (v1): воздушная индуктивность по формуле
Уилера для многослойных катушек, скорректированная на присутствие
ферромагнитного снаряда через коэффициент размагничивания конечного
цилиндра (аппроксимация вытянутым сфероидом). Точная картина требует МКЭ
(FEA) — это явно вынесено как ограничение модели, см. PLAN.md.
"""
import math
from dataclasses import dataclass
import numpy as np
INCH_M = 0.0254
@dataclass(frozen=True)
class WindingGeometry:
coil_length_m: float
radial_depth_m: float
mean_radius_m: float
total_turns: int
total_wire_length_m: float
def winding_geometry(
tube_od_m: float, wire_od_m: float, turns_per_layer: int, layers: int
) -> WindingGeometry:
"""Геометрия обмотки из параметров намотки, а не свободной "длины катушки"."""
coil_length_m = turns_per_layer * wire_od_m
radial_depth_m = layers * wire_od_m
inner_radius_m = tube_od_m / 2
total_wire_length_m = 0.0
for i in range(layers):
layer_radius_m = inner_radius_m + (i + 0.5) * wire_od_m
total_wire_length_m += 2 * math.pi * layer_radius_m * turns_per_layer
mean_radius_m = inner_radius_m + radial_depth_m / 2
return WindingGeometry(
coil_length_m=coil_length_m,
radial_depth_m=radial_depth_m,
mean_radius_m=mean_radius_m,
total_turns=turns_per_layer * layers,
total_wire_length_m=total_wire_length_m,
)
def air_core_inductance_wheeler(mean_radius_m: float, length_m: float, depth_m: float, turns: int) -> float:
"""Формула Уилера для многослойного воздушного соленоида (Гн).
Точность ~1% при a, l, c одного порядка величины (Wheeler, 1928).
"""
a_in = mean_radius_m / INCH_M
l_in = length_m / INCH_M
c_in = depth_m / INCH_M
l_uh = (0.8 * a_in**2 * turns**2) / (6 * a_in + 9 * l_in + 10 * c_in)
return l_uh * 1e-6
def demagnetizing_factor_prolate(aspect_ratio: float) -> float:
"""Коэффициент размагничивания вдоль длинной оси вытянутого сфероида.
aspect_ratio = длина / диаметр снаряда. При aspect_ratio -> 1 (сфера)
возвращает 1/3; при aspect_ratio -> inf стремится к 0.
"""
m = max(aspect_ratio, 1.0 + 1e-6)
root = math.sqrt(m**2 - 1)
return (1 / (m**2 - 1)) * ((m / (2 * root)) * math.log((m + root) / (m - root)) - 1)
def effective_permeability(mu_r: float, demagnetizing_factor: float) -> float:
return mu_r / (1 + demagnetizing_factor * (mu_r - 1))
def _sigmoid(z: np.ndarray) -> np.ndarray:
return 0.5 * (1 + np.tanh(z / 2))
@dataclass(frozen=True)
class CoilInductanceModel:
"""L(x) и аналитическая dL/dx для катушки с проходящим ферромагнитным снарядом.
x — положение центра снаряда относительно центра катушки (м), катушка
считается центрированной в x=0 в своей локальной системе координат.
"""
l_air_h: float
coil_length_m: float
slug_length_m: float
mu_eff: float
smoothing_width_m: float
def _edges(self) -> tuple[float, float]:
half_span = (self.coil_length_m + self.slug_length_m) / 2
return -half_span, half_span
def overlap_fraction(self, x: np.ndarray) -> np.ndarray:
e1, e2 = self._edges()
w = self.smoothing_width_m
return _sigmoid((x - e1) / w) * _sigmoid((e2 - x) / w)
def d_overlap_dx(self, x: np.ndarray) -> np.ndarray:
e1, e2 = self._edges()
w = self.smoothing_width_m
s1 = _sigmoid((x - e1) / w)
s2 = _sigmoid((e2 - x) / w)
return (s1 * s2 / w) * (s2 - s1)
def l_of_x(self, x: np.ndarray) -> np.ndarray:
return self.l_air_h * (1 + (self.mu_eff - 1) * self.overlap_fraction(x))
def dl_dx(self, x: np.ndarray) -> np.ndarray:
return self.l_air_h * (self.mu_eff - 1) * self.d_overlap_dx(x)

11
tests/fixtures/data/capacitors.json vendored Normal file
View File

@@ -0,0 +1,11 @@
[
{
"part_number": "test-cap-450v-220uf",
"capacitance_uf": 220.0,
"voltage_v": 450.0,
"esr_ohm": 0.05,
"max_current_a": 200.0,
"price": 150.0,
"source": "fixture"
}
]

View File

@@ -0,0 +1,10 @@
[
{
"name": "test-steel",
"density_kg_m3": 7850.0,
"mu_r": 200.0,
"b_sat_tesla": 1.8,
"price_per_kg": 100.0,
"source": "fixture"
}
]

18
tests/fixtures/data/sensors.json vendored Normal file
View File

@@ -0,0 +1,18 @@
[
{
"part_number": "test-optical",
"kind": "optical",
"propagation_delay_ns": 500.0,
"price": 20.0,
"source": "fixture"
},
{
"part_number": "test-inductive-pickup",
"kind": "inductive",
"propagation_delay_ns": 0.0,
"price": 5.0,
"source": "fixture",
"sensitivity_v_per_mps": 0.02,
"threshold_v": 0.1
}
]

13
tests/fixtures/data/switches.json vendored Normal file
View File

@@ -0,0 +1,13 @@
[
{
"part_number": "test-scr-25a",
"kind": "SCR",
"max_current_a": 25.0,
"max_voltage_v": 600.0,
"on_resistance_ohm": null,
"on_voltage_drop_v": 1.5,
"turn_on_time_ns": 2000.0,
"price": 80.0,
"source": "fixture"
}
]

22
tests/fixtures/data/wires.json vendored Normal file
View File

@@ -0,0 +1,22 @@
[
{
"part_id": "test-cu-0.5",
"material": "copper",
"gauge_mm": 0.5,
"insulation_od_mm": 0.55,
"resistivity_ohm_m": 1.72e-8,
"max_current_a": 3.0,
"price_per_m": 5.0,
"source": "fixture"
},
{
"part_id": "test-al-0.5",
"material": "aluminum",
"gauge_mm": 0.5,
"insulation_od_mm": 0.55,
"resistivity_ohm_m": 2.82e-8,
"max_current_a": 2.2,
"price_per_m": 2.0,
"source": "fixture"
}
]

View File

@@ -0,0 +1,116 @@
"""Проверка ОДУ разряда против аналитического решения простого 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)

View File

@@ -0,0 +1,45 @@
from pathlib import Path
from gausse.components.database import ComponentDatabase
FIXTURES = Path(__file__).parent / "fixtures" / "data"
def test_loads_all_component_types():
db = ComponentDatabase.load(FIXTURES)
assert len(db.wires) == 2
assert len(db.capacitors) == 1
assert len(db.switches) == 1
assert len(db.sensors) == 2
assert len(db.projectile_materials) == 1
def test_wires_by_material_filters_correctly():
db = ComponentDatabase.load(FIXTURES)
copper = db.wires_by_material("copper")
assert len(copper) == 1
assert copper[0].part_id == "test-cu-0.5"
def test_capacitors_rated_at_least_filters_by_voltage():
db = ComponentDatabase.load(FIXTURES)
assert len(db.capacitors_rated_at_least(400.0)) == 1
assert len(db.capacitors_rated_at_least(500.0)) == 0
def test_switches_rated_for_current_and_voltage():
db = ComponentDatabase.load(FIXTURES)
assert len(db.switches_rated_for(20.0, 500.0)) == 1
assert len(db.switches_rated_for(30.0, 500.0)) == 0
def test_sensors_by_kind():
db = ComponentDatabase.load(FIXTURES)
assert len(db.sensors_by_kind("optical")) == 1
assert len(db.sensors_by_kind("inductive")) == 1
def test_missing_data_dir_returns_empty_database(tmp_path):
db = ComponentDatabase.load(tmp_path)
assert db.wires == []
assert db.capacitors == []

35
tests/test_force_model.py Normal file
View File

@@ -0,0 +1,35 @@
import pytest
from gausse.physics.force import (
force_on_slug_newtons,
saturation_scale,
solenoid_field_estimate_tesla,
)
def test_field_estimate_scales_with_current_and_turns():
b_low = solenoid_field_estimate_tesla(mu_eff=100, total_turns=100, coil_length_m=0.05, current_a=10)
b_high = solenoid_field_estimate_tesla(mu_eff=100, total_turns=100, coil_length_m=0.05, current_a=100)
assert b_high > b_low
def test_saturation_scale_is_one_below_bsat():
assert saturation_scale(b_estimate_tesla=0.5, b_sat_tesla=1.8) == pytest.approx(1.0)
def test_saturation_scale_clamps_above_bsat():
scale = saturation_scale(b_estimate_tesla=3.6, b_sat_tesla=1.8)
assert scale == pytest.approx(0.5)
def test_force_is_reduced_once_saturated():
common = dict(dl_dx_unsaturated=1e-3, mu_eff=200, total_turns=300, coil_length_m=0.05, b_sat_tesla=1.8)
force_unsaturated = force_on_slug_newtons(current_a=5, **common)
force_saturated = force_on_slug_newtons(current_a=500, **common)
# без клэмпа сила росла бы как I^2 (в 10000 раз); с насыщением рост должен быть намного меньше
assert force_saturated / force_unsaturated < 5000
def test_force_zero_at_zero_current():
common = dict(dl_dx_unsaturated=1e-3, mu_eff=200, total_turns=300, coil_length_m=0.05, b_sat_tesla=1.8)
assert force_on_slug_newtons(current_a=0, **common) == pytest.approx(0.0)

82
tests/test_inductance.py Normal file
View File

@@ -0,0 +1,82 @@
import math
import numpy as np
import pytest
from gausse.physics.inductance import (
CoilInductanceModel,
air_core_inductance_wheeler,
demagnetizing_factor_prolate,
effective_permeability,
winding_geometry,
)
def test_winding_geometry_scales_with_turns_and_layers():
geo = winding_geometry(tube_od_m=0.02, wire_od_m=0.001, turns_per_layer=50, layers=3)
assert geo.total_turns == 150
assert geo.coil_length_m == pytest.approx(0.05)
assert geo.radial_depth_m == pytest.approx(0.003)
assert geo.total_wire_length_m > 0
# больше слоёв -> больше провода при тех же витках на слой
geo_more_layers = winding_geometry(0.02, 0.001, 50, 5)
assert geo_more_layers.total_wire_length_m > geo.total_wire_length_m
def test_wheeler_inductance_is_positive_and_grows_with_turns():
l_small = air_core_inductance_wheeler(mean_radius_m=0.015, length_m=0.05, depth_m=0.003, turns=50)
l_large = air_core_inductance_wheeler(mean_radius_m=0.015, length_m=0.05, depth_m=0.003, turns=200)
assert l_small > 0
assert l_large > l_small
def test_demagnetizing_factor_bounds():
# почти сфера (aspect_ratio -> 1): классическое значение 1/3
assert demagnetizing_factor_prolate(1.001) == pytest.approx(1 / 3, abs=0.02)
# очень вытянутый стержень: коэффициент размагничивания стремится к 0
assert demagnetizing_factor_prolate(200.0) < 0.01
def test_effective_permeability_limits():
# без размагничивания (Nd=0) эффективная проницаемость равна мю материала
assert effective_permeability(mu_r=500.0, demagnetizing_factor=0.0) == pytest.approx(500.0)
# полное размагничивание (Nd=1) гасит эффект проницаемости до 1
assert effective_permeability(mu_r=500.0, demagnetizing_factor=1.0) == pytest.approx(1.0)
@pytest.fixture
def model() -> CoilInductanceModel:
return CoilInductanceModel(
l_air_h=100e-6,
coil_length_m=0.05,
slug_length_m=0.03,
mu_eff=50.0,
smoothing_width_m=0.001,
)
def test_inductance_peaks_when_slug_centered(model):
l_center = model.l_of_x(np.array(0.0))
l_far = model.l_of_x(np.array(0.2))
assert l_center > l_far
assert l_far == pytest.approx(model.l_air_h, rel=1e-3)
def test_dl_dx_vanishes_at_symmetric_center(model):
assert model.dl_dx(np.array(0.0)) == pytest.approx(0.0, abs=1e-9)
def test_dl_dx_sign_flips_across_center(model):
# слева от центра индуктивность растёт (dL/dx > 0), справа убывает (< 0)
assert model.dl_dx(np.array(-0.01)) > 0
assert model.dl_dx(np.array(0.01)) < 0
def test_analytic_dl_dx_matches_finite_difference(model):
xs = np.linspace(-0.06, 0.06, 25)
h = 1e-7
analytic = np.array([model.dl_dx(x) for x in xs])
finite_diff = np.array(
[(model.l_of_x(x + h) - model.l_of_x(x - h)) / (2 * h) for x in xs]
)
assert np.allclose(analytic, finite_diff, rtol=1e-3, atol=1e-9)