78af686e6d
Разрывный sign(v) в силе трения заставлял адаптивный solve_ivp бесконечно дробить шаг у v≈0 (chattering): одна конфигурация считалась 35+ минут, скорость эволюции упала с ~430/с до 3-5/с. Стандартная регуляризация sign(v)->tanh(v/0.01) во всех путях (CPU-разряд, подлёт, GPU numpy, fused cupy-ядро) + регрессионный тест на гладкость силы в нуле. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
167 lines
7.0 KiB
Python
167 lines
7.0 KiB
Python
"""Физреализм v2: скин/близость (Доуэлл), трение+воздух, импульсные рейтинги ключей."""
|
||
|
||
import math
|
||
|
||
import numpy as np
|
||
|
||
from gausse.components.database import ComponentDatabase
|
||
from gausse.components.schema import (
|
||
CapacitorSpec,
|
||
ProjectileMaterialSpec,
|
||
SensorSpec,
|
||
SwitchSpec,
|
||
WireSpec,
|
||
)
|
||
from gausse.physics.ac_resistance import dowell_ac_factor, skin_depth_m
|
||
from gausse.physics.constants import switch_pulse_limit_a
|
||
from gausse.sim.stage import ProjectileConfig, StageConfig, build_stage_physics, run_stage
|
||
|
||
DB = ComponentDatabase.load()
|
||
|
||
|
||
def make_projectile() -> ProjectileConfig:
|
||
steel = ProjectileMaterialSpec(
|
||
name="steel", density_kg_m3=7850.0, mu_r=200.0, b_sat_tesla=1.8, price_per_kg=100.0, source="test"
|
||
)
|
||
return ProjectileConfig(material=steel, diameter_m=0.008, length_m=0.02)
|
||
|
||
|
||
def make_stage() -> StageConfig:
|
||
wire = WireSpec(
|
||
part_id="w1", material="copper", gauge_mm=0.8, insulation_od_mm=0.85,
|
||
resistivity_ohm_m=1.68e-8, max_current_a=10.0, price_per_m=3.0, source="test",
|
||
)
|
||
capacitor = CapacitorSpec(
|
||
part_number="c1", capacitance_uf=100.0, voltage_v=400.0, esr_ohm=0.05, max_current_a=300.0,
|
||
price=300.0, source="test",
|
||
)
|
||
switch = SwitchSpec(
|
||
part_number="sw1", kind="MOSFET", max_current_a=500.0, max_voltage_v=500.0,
|
||
on_resistance_ohm=0.02, on_voltage_drop_v=None, turn_on_time_ns=50.0, price=50.0, source="test",
|
||
)
|
||
sensor = SensorSpec(part_number="s1", kind="optical", propagation_delay_ns=500.0, price=20.0, source="test")
|
||
return StageConfig(
|
||
wire=wire, capacitor=capacitor, switch=switch, sensor=sensor,
|
||
tube_od_m=0.01, turns_per_layer=20, layers=4,
|
||
sensor_to_coil_distance_m=0.02, charge_voltage_v=350.0,
|
||
)
|
||
|
||
|
||
# --- Доуэлл ---
|
||
|
||
|
||
def test_skin_depth_copper_50hz():
|
||
# медь на 50 Гц: δ ≈ 9.3 мм (классическое справочное значение)
|
||
delta = skin_depth_m(1.68e-8, 2 * math.pi * 50)
|
||
assert abs(delta - 9.2e-3) < 0.5e-3
|
||
|
||
|
||
def test_dowell_limits():
|
||
# низкая частота -> фактор стремится к 1
|
||
assert dowell_ac_factor(0.001, 0.00108, 8, 1.68e-8, 2 * math.pi * 1.0) < 1.001
|
||
# больше слоёв -> сильнее эффект близости (монотонность)
|
||
omega = 2 * math.pi * 2000
|
||
f = [dowell_ac_factor(0.002, 0.00208, m, 1.68e-8, omega) for m in (1, 4, 10)]
|
||
assert f[0] < f[1] < f[2]
|
||
# всегда >= 1
|
||
assert all(x >= 1.0 for x in f)
|
||
|
||
|
||
def test_dowell_series_matches_formula_at_boundary():
|
||
# непрерывность на стыке ряда и полной формулы (Δ≈0.25)
|
||
omega_lo, omega_hi = None, None
|
||
# подберём частоты чуть ниже/выше границы Δ=0.25 для d=1мм
|
||
d, rho = 0.001, 1.68e-8
|
||
for omega in np.logspace(2, 6, 4000):
|
||
delta = skin_depth_m(rho, omega)
|
||
big = (math.pi / 4) ** 0.75 * d / delta # porosity~1
|
||
if big < 0.25:
|
||
omega_lo = omega
|
||
elif omega_hi is None:
|
||
omega_hi = omega
|
||
break
|
||
f_lo = dowell_ac_factor(d, d, 6, rho, omega_lo)
|
||
f_hi = dowell_ac_factor(d, d, 6, rho, omega_hi)
|
||
assert abs(f_lo - f_hi) < 0.01
|
||
|
||
|
||
def test_thick_wire_many_layers_significant_ac_factor():
|
||
# толстый провод (2мм) в 10 слоёв на кГц — эффект близости должен быть
|
||
# существенным (>1.3), иначе канал потерь не работает
|
||
f = dowell_ac_factor(0.002, 0.00208, 10, 1.68e-8, 2 * math.pi * 1500)
|
||
assert f > 1.3
|
||
|
||
|
||
# --- трение и воздух ---
|
||
|
||
|
||
def test_stage_with_friction_slower_than_ideal():
|
||
"""Выходная скорость с трением/воздухом ниже, чем была бы без них."""
|
||
stage = make_stage()
|
||
projectile = make_projectile()
|
||
res = run_stage(-0.06, 3.0, stage, projectile)
|
||
assert res.feasible, res.reason
|
||
# подлёт 3 м/с на ~3 см: только трение уже отъедает заметную скорость
|
||
phys = build_stage_physics(stage, projectile)
|
||
assert phys.circuit_params.retard_const_n > 0
|
||
assert phys.circuit_params.drag_coeff_n_per_mps2 > 0
|
||
|
||
|
||
def test_slow_projectile_stalls_before_sensor():
|
||
"""Медленный снаряд далеко от датчика останавливается трением — честный отказ."""
|
||
stage = make_stage()
|
||
projectile = make_projectile()
|
||
res = run_stage(-1.0, 0.15, stage, projectile) # 15 см/с за метр до катушки
|
||
assert not res.feasible
|
||
assert "трением" in res.reason
|
||
|
||
|
||
def test_energy_balance_with_friction():
|
||
"""Баланс: E_in = E_ост.конденсатора + E_потерь(включая трение) + ΔE_кин."""
|
||
stage = make_stage()
|
||
projectile = make_projectile()
|
||
res = run_stage(-0.06, 3.0, stage, projectile)
|
||
assert res.feasible, res.reason
|
||
balance = res.energy_remaining_cap_j + res.energy_dissipated_j + res.kinetic_energy_delta_j
|
||
assert abs(balance - res.energy_in_j) < 0.02 * res.energy_in_j
|
||
|
||
|
||
# --- импульсные рейтинги ключей ---
|
||
|
||
|
||
def test_pulse_ratings_loaded_from_datasheets():
|
||
by_pn = {s.part_number: s for s in DB.switches}
|
||
assert by_pn["IRFP254PBF"].pulse_current_a == 92.0 # IDM из даташита
|
||
assert by_pn["BT151-650R"].pulse_current_a == 100.0 # ITSM
|
||
assert by_pn["IRG4PC50F"].pulse_current_a == 280.0 # ICM
|
||
for s in DB.switches:
|
||
limit = switch_pulse_limit_a(s)
|
||
assert limit >= s.max_current_a # импульсный >= продолжительного
|
||
|
||
|
||
def test_pulse_limit_fallback_without_datasheet():
|
||
from dataclasses import replace
|
||
|
||
sw = replace(DB.switches[0], pulse_current_a=None)
|
||
# без даташита — консервативный множитель по типу ключа
|
||
assert switch_pulse_limit_a(sw) == sw.max_current_a * {"SCR": 10.0, "MOSFET": 4.0, "IGBT": 3.0}[sw.kind]
|
||
|
||
|
||
def test_friction_force_smooth_at_zero_velocity():
|
||
"""Регрессия на завис решателя: сила трения обязана быть ГЛАДКОЙ в v=0.
|
||
|
||
Разрывный sign(v) заставлял solve_ivp бесконечно дробить шаг у v≈0
|
||
(одна конфигурация считалась 35+ минут на сервере).
|
||
"""
|
||
from gausse.physics.circuit import StageCircuitParams, retarding_force_n
|
||
|
||
p = StageCircuitParams(
|
||
capacitance_f=1e-4, r_total_ohm=0.1, mass_kg=0.01,
|
||
retard_const_n=0.034, drag_coeff_n_per_mps2=2e-5,
|
||
)
|
||
# непрерывность: около нуля сила ~ 0, нечётная, без скачка
|
||
assert abs(retarding_force_n(p, 1e-6)) < 1e-5
|
||
assert abs(retarding_force_n(p, 1e-6) + retarding_force_n(p, -1e-6)) < 1e-12
|
||
# вдали от нуля выходит на полное трение
|
||
assert retarding_force_n(p, 1.0) > 0.9 * p.retard_const_n
|