Физреализм v2: скин/близость (Доуэлл), трение+воздух, паспортные импульсные токи ключей

- physics/ac_resistance.py: R_AC обмотки по Доуэллу на частоте импульса
  1/sqrt(LC) — для толстого провода во многих слоях потери в разы выше DC
- трение о трубку (0.35·m·g) + аэродинамика (0.5·rho·Cd·A·v^2) во всей
  динамике: CPU ОДУ разряда, подлёт к датчику (с событием остановки),
  GPU numpy-путь, fused cupy-ядро, аналитический coast GPU-sweep'а
- SwitchSpec.pulse_current_a: паспортные ITSM/IDM/ICM из даташитов вместо
  generic-множителей; отчёт теперь различает превышение продолжительного
  рейтинга (норма для импульса) и импульсного предела (отбраковка) — фикс
  вводившего в заблуждение флага switch_current_over_limit
- КПД теперь может быть слегка отрицательным (трение съело больше, чем
  добавила слабая катушка) — это честно
- MODEL_VERSION -> gausse-physics-v2

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
This commit is contained in:
jze9
2026-07-08 00:57:07 +05:00
parent 83a4a7b805
commit d4affa9574
15 changed files with 417 additions and 59 deletions

View File

@@ -8,8 +8,9 @@
"on_voltage_drop_v": 1.5,
"turn_on_time_ns": 10000.0,
"price": 140.0,
"source": "REAL: КУ202Н, chipdip.ru (https://www.chipdip.ru/product/ku202n) 140₽ / procontact74.ru (https://procontact74.ru/44-tiristory/ku202n-tiristor-10a-400v-44011/) 189₽; Von=1.5В с даташита; время включения ОЦЕНКА ~10мкс (типично для тиристора этого класса, точное значение в листинге не приведено)",
"url": "https://www.chipdip.ru/product/ku202n"
"source": "REAL: КУ202Н, chipdip.ru (https://www.chipdip.ru/product/ku202n) 140₽ / procontact74.ru (https://procontact74.ru/44-tiristory/ku202n-tiristor-10a-400v-44011/) 189₽; Von=1.5В с даташита; время включения ОЦЕНКА ~10мкс (типично для тиристора этого класса, точное значение в листинге не приведено); ITSM=100А (ударный неповторяющийся, 10мс) по справочнику КУ202Н",
"url": "https://www.chipdip.ru/product/ku202n",
"pulse_current_a": 100.0
},
{
"part_number": "BT151-650R",
@@ -20,8 +21,9 @@
"on_voltage_drop_v": 1.7,
"turn_on_time_ns": 1500.0,
"price": 190.0,
"source": "REAL: chipdip.ru (https://www.chipdip.ru/product/bt151-650r), ~180-200₽ (BT151-600R вариант подтверждён на 180₽); Von и время включения ОЦЕНКА (типично для тиристора этого класса, не указаны продавцом)",
"url": "https://www.chipdip.ru/product/bt151-650r"
"source": "REAL: chipdip.ru (https://www.chipdip.ru/product/bt151-650r), ~180-200₽ (BT151-600R вариант подтверждён на 180₽); Von и время включения ОЦЕНКА (типично для тиристора этого класса, не указаны продавцом); ITSM=100А (даташит NXP, полусинус 10мс)",
"url": "https://www.chipdip.ru/product/bt151-650r",
"pulse_current_a": 100.0
},
{
"part_number": "IRFP254PBF",
@@ -32,8 +34,9 @@
"on_voltage_drop_v": null,
"turn_on_time_ns": 30.0,
"price": 39.9,
"source": "REAL: procontact74.ru (https://procontact74.ru/01-elektronnyie-komponentyi-/43-tranzistory/polevyie--mosfet-/58406-tranzistorirfp254pbfmosfetnch250v23ato247acoriginalbydemontajnogi10mm/); ВНИМАНИЕ: демонтированная (Б/У) деталь, не новая партия; время включения ОЦЕНКА ~30нс (типично для этого класса MOSFET)",
"url": "https://procontact74.ru/01-elektronnyie-komponentyi-/43-tranzistory/polevyie--mosfet-/58406-tranzistorirfp254pbfmosfetnch250v23ato247acoriginalbydemontajnogi10mm/"
"source": "REAL: procontact74.ru (https://procontact74.ru/01-elektronnyie-komponentyi-/43-tranzistory/polevyie--mosfet-/58406-tranzistorirfp254pbfmosfetnch250v23ato247acoriginalbydemontajnogi10mm/); ВНИМАНИЕ: демонтированная (Б/У) деталь, не новая партия; время включения ОЦЕНКА ~30нс (типично для этого класса MOSFET); IDM=92А (импульсный ток стока, даташит)",
"url": "https://procontact74.ru/01-elektronnyie-komponentyi-/43-tranzistory/polevyie--mosfet-/58406-tranzistorirfp254pbfmosfetnch250v23ato247acoriginalbydemontajnogi10mm/",
"pulse_current_a": 92.0
},
{
"part_number": "IRFP250PBF",
@@ -44,8 +47,9 @@
"on_voltage_drop_v": null,
"turn_on_time_ns": 35.0,
"price": 320.0,
"source": "REAL: chipdip.ru (https://www.chipdip.ru/product/irfp250), новая деталь, 320₽ (276₽ от 15шт, 245₽ от 150шт); Rds(on)=0.085Ом при 18А/10В с даташита; время включения ОЦЕНКА",
"url": "https://www.chipdip.ru/product/irfp250"
"source": "REAL: chipdip.ru (https://www.chipdip.ru/product/irfp250), новая деталь, 320₽ (276₽ от 15шт, 245₽ от 150шт); Rds(on)=0.085Ом при 18А/10В с даташита; время включения ОЦЕНКА; IDM=120А (импульсный ток стока, даташит)",
"url": "https://www.chipdip.ru/product/irfp250",
"pulse_current_a": 120.0
},
{
"part_number": "IRFP4468PBF",
@@ -56,8 +60,9 @@
"on_voltage_drop_v": null,
"turn_on_time_ns": 40.0,
"price": 610.0,
"source": "REAL: chipdip.ru (https://www.chipdip.ru/product/irfp4468pbf), 610₽ (534₽ от 15шт); паспортный максимум 290А при 25°C, взято консервативное практическое значение 180А; Rds(on)=0.0026Ом при 180А/10В с даташита; время включения ОЦЕНКА",
"url": "https://www.chipdip.ru/product/irfp4468pbf"
"source": "REAL: chipdip.ru (https://www.chipdip.ru/product/irfp4468pbf), 610₽ (534₽ от 15шт); паспортный максимум 290А при 25°C, взято консервативное практическое значение 180А; Rds(on)=0.0026Ом при 180А/10В с даташита; время включения ОЦЕНКА; IDM=1080А (импульсный ток стока, даташит)",
"url": "https://www.chipdip.ru/product/irfp4468pbf",
"pulse_current_a": 1080.0
},
{
"part_number": "IRG4PC50F",
@@ -68,7 +73,8 @@
"on_voltage_drop_v": 1.8,
"turn_on_time_ns": 300.0,
"price": 504.0,
"source": "REAL: chipdip.ru (https://www.chipdip.ru/product/irg4pc50f-8048990896), ~504₽; Vce(sat) и время включения ОЦЕНКА (типично для этого семейства IGBT, не подтверждено со страницы товара)",
"url": "https://www.chipdip.ru/product/irg4pc50f-8048990896"
"source": "REAL: chipdip.ru (https://www.chipdip.ru/product/irg4pc50f-8048990896), ~504₽; Vce(sat) и время включения ОЦЕНКА (типично для этого семейства IGBT, не подтверждено со страницы товара); ICM=280А (импульсный ток коллектора, даташит)",
"url": "https://www.chipdip.ru/product/irg4pc50f-8048990896",
"pulse_current_a": 280.0
}
]
]

View File

@@ -83,6 +83,9 @@ class SwitchSpec:
price: float
source: str
url: str | None = None
# паспортный импульсный ток из даташита (ITSM тиристора / IDM MOSFET /
# ICM IGBT). Если None — применяется типовой множитель SWITCH_SURGE_FACTOR.
pulse_current_a: float | None = None
@dataclass(frozen=True)

View File

@@ -35,6 +35,8 @@ class BatchDischargeParams:
smoothing_width_m: "any"
i_sat_a: "any" # ток насыщения; NO_SATURATION_I_SAT где железа нет
r_eddy_coeff_ohm: "any" # вихревое сопротивление снаряда (× overlap)
retard_const_n: "any" # сухое трение о трубку (Н, против движения)
drag_coeff_n: "any" # аэродинамика: сила = drag_coeff·v² (против движения)
def params_from_models(xp, models, circuit_params) -> BatchDischargeParams:
@@ -62,6 +64,8 @@ def params_from_models(xp, models, circuit_params) -> BatchDischargeParams:
smoothing_width_m=col([m.smoothing_width_m for m in models]),
i_sat_a=col(i_sat),
r_eddy_coeff_ohm=col([getattr(c, "r_eddy_coeff_ohm", 0.0) for c in circuit_params]),
retard_const_n=col([getattr(c, "retard_const_n", 0.0) for c in circuit_params]),
drag_coeff_n=col([getattr(c, "drag_coeff_n_per_mps2", 0.0) for c in circuit_params]),
)
@@ -104,7 +108,8 @@ def _derivatives(xp, q, i, x, v, p: BatchDischargeParams):
v_c = q / p.capacitance_f
d_q = -i
d_i = (v_c - i * r_eff - dlambda_dx * v) / dlambda_di
d_v = force / p.mass_kg
retard = xp.sign(v) * (p.retard_const_n + p.drag_coeff_n * v * v) # трение + воздух
d_v = (force - retard) / p.mass_kg
d_x = v
return d_q, d_i, d_x, d_v
@@ -128,7 +133,7 @@ def _get_fused_step(cp):
s2 = 0.5 * (1.0 + cp.tanh((hs - x) / w * 0.5))
return s1 * s2, (s1 * s2 / w) * (s2 - s1)
def _der(q, i, x, v, hs, w, isat, lair, liron, C, rt, red, m):
def _der(q, i, x, v, hs, w, isat, lair, liron, C, rt, red, fc, dc, m):
ov, dov = _ovl(x, hs, w)
z = i / isat
tz = cp.tanh(z)
@@ -140,15 +145,15 @@ def _get_fused_step(cp):
reff = rt + red * ov
dq = -i
di = (q / C - i * reff - liron * dov * g * v) / dl_di
dv = (liron * dov * G) / m
dv = (liron * dov * G - cp.sign(v) * (fc + dc * v * v)) / m
return dq, di, v, dv
@cp.fuse()
def step(q, i, x, v, dt, hs, w, isat, lair, liron, C, rt, red, m):
a = _der(q, i, x, v, hs, w, isat, lair, liron, C, rt, red, m)
b = _der(q + dt * 0.5 * a[0], i + dt * 0.5 * a[1], x + dt * 0.5 * a[2], v + dt * 0.5 * a[3], hs, w, isat, lair, liron, C, rt, red, m)
c = _der(q + dt * 0.5 * b[0], i + dt * 0.5 * b[1], x + dt * 0.5 * b[2], v + dt * 0.5 * b[3], hs, w, isat, lair, liron, C, rt, red, m)
d = _der(q + dt * c[0], i + dt * c[1], x + dt * c[2], v + dt * c[3], hs, w, isat, lair, liron, C, rt, red, m)
def step(q, i, x, v, dt, hs, w, isat, lair, liron, C, rt, red, fc, dc, m):
a = _der(q, i, x, v, hs, w, isat, lair, liron, C, rt, red, fc, dc, m)
b = _der(q + dt * 0.5 * a[0], i + dt * 0.5 * a[1], x + dt * 0.5 * a[2], v + dt * 0.5 * a[3], hs, w, isat, lair, liron, C, rt, red, fc, dc, m)
c = _der(q + dt * 0.5 * b[0], i + dt * 0.5 * b[1], x + dt * 0.5 * b[2], v + dt * 0.5 * b[3], hs, w, isat, lair, liron, C, rt, red, fc, dc, m)
d = _der(q + dt * c[0], i + dt * c[1], x + dt * c[2], v + dt * c[3], hs, w, isat, lair, liron, C, rt, red, fc, dc, m)
qn = q + dt / 6.0 * (a[0] + 2 * b[0] + 2 * c[0] + d[0])
in_ = i + dt / 6.0 * (a[1] + 2 * b[1] + 2 * c[1] + d[1])
xn = x + dt / 6.0 * (a[2] + 2 * b[2] + 2 * c[2] + d[2])
@@ -210,7 +215,8 @@ def integrate_batch_discharge(
q_new, i_new, x_new, v_new, overlap_old, overlap_new = fused(
q_old, i_old, x_old, v_old, dt, hs, params.smoothing_width_m,
params.i_sat_a, params.l_air_h, params.l_iron_coeff, params.capacitance_f,
params.r_total_ohm, params.r_eddy_coeff_ohm, params.mass_kg,
params.r_total_ohm, params.r_eddy_coeff_ohm,
params.retard_const_n, params.drag_coeff_n, params.mass_kg,
)
else:
k1 = _derivatives(xp, q_old, i_old, x_old, v_old, params)
@@ -228,6 +234,10 @@ def integrate_batch_discharge(
r_eff_old = params.r_total_ohm + params.r_eddy_coeff_ohm * overlap_old
r_eff_new = params.r_total_ohm + params.r_eddy_coeff_ohm * overlap_new
step_diss = 0.5 * (i_old**2 * r_eff_old + i_new**2 * r_eff_new) * dt
# работа трения/воздуха: P = (F_тр + drag·v²)·|v|
fric_old = (params.retard_const_n + params.drag_coeff_n * v_old**2) * xp.abs(v_old)
fric_new = (params.retard_const_n + params.drag_coeff_n * v_new**2) * xp.abs(v_new)
step_diss = step_diss + 0.5 * (fric_old + fric_new) * dt
energy_diss = energy_diss + xp.where(active, step_diss, 0.0)
peak_current = xp.maximum(peak_current, xp.where(active, xp.abs(i_new), 0.0))
past_peak = past_peak | (active & (i_old > 1.0) & (i_new < i_old))

View File

@@ -27,7 +27,7 @@ from gausse.gpu.batch_integrator import integrate_batch_discharge, params_from_m
from gausse.optim.objective import MODEL_VERSION, build_detail
from gausse.optim.progress_log import ProgressLogger, default_log_path
from gausse.optim.search_space import SearchBounds, decode, genome_to_dict, sample_genome
from gausse.physics.constants import SWITCH_SURGE_FACTOR
from gausse.physics.constants import switch_pulse_limit_a
from gausse.sim.coilgun import CoilgunResult
from gausse.sim.stage import build_stage_physics
from gausse.storage.database import insert_runs, open_connection
@@ -51,6 +51,23 @@ class _State:
reason: str | None = None
def _coast_velocity(v0: float, distance_m: float, retard_const_n: float, drag_coeff: float, mass_kg: float) -> float:
"""Скорость после подлёта на distance_m с трением и воздухом (точное решение).
m·v·dv/dx = (a + b·v²) ⇒ v²(d) = (v0² + a/b)·e^{2bd/m} a/b (b>0),
либо v² = v0² 2ad/m при b=0. Возвращает 0.0, если снаряд останавливается.
Та же физика, что численный подлёт в sim/stage.py (там ОДУ, тут квадратура).
"""
if distance_m <= 0:
return v0
a, b, m = retard_const_n, drag_coeff, mass_kg
if b <= 0:
u = v0 * v0 - 2.0 * a * distance_m / m
else:
u = (v0 * v0 + a / b) * math.exp(-2.0 * b * distance_m / m) - a / b
return math.sqrt(u) if u > 0 else 0.0
def _prepare_stage(stage, projectile, cur_x_global, cur_v, coil_center_global):
"""Готовит разряд ступени: физика + аналитический fire-state из текущего (x,v)."""
if cur_v <= 1e-6:
@@ -59,15 +76,22 @@ def _prepare_stage(stage, projectile, cur_x_global, cur_v, coil_center_global):
sensor_global = coil_center_global - stage.sensor_to_coil_distance_m
if sensor_global < cur_x_global - 1e-9:
return None, "датчик позади снаряда (не сработает)"
# подлёт до датчика — с трением о трубку и воздухом (та же физика, что на CPU)
cp_ = phys.circuit_params
v_at_sensor = _coast_velocity(
cur_v, sensor_global - cur_x_global, cp_.retard_const_n, cp_.drag_coeff_n_per_mps2, projectile.mass_kg
)
if v_at_sensor <= 1e-3:
return None, "снаряд остановлен трением о трубку, не долетев до датчика"
if stage.sensor.kind == "inductive":
if stage.sensor.sensitivity_v_per_mps * cur_v < (stage.sensor.threshold_v or 0.0):
if stage.sensor.sensitivity_v_per_mps * v_at_sensor < (stage.sensor.threshold_v or 0.0):
return None, "инд. датчик: сигнал ниже порога"
fire_delay = (stage.sensor.propagation_delay_ns + stage.switch.turn_on_time_ns) * 1e-9
x_fire_global = sensor_global + cur_v * fire_delay
x_fire_global = sensor_global + v_at_sensor * fire_delay
x_fire_local = x_fire_global - coil_center_global
q0 = phys.capacitance_f * stage.charge_voltage_v
return {
"phys": phys, "q0": q0, "x_fire": x_fire_local, "v_fire": cur_v,
"phys": phys, "q0": q0, "x_fire": x_fire_local, "v_fire": v_at_sensor,
"energy_in": 0.5 * phys.capacitance_f * stage.charge_voltage_v**2,
"coil_center": coil_center_global, "stage": stage, "mass": projectile.mass_kg,
}, None
@@ -136,10 +160,10 @@ def run_gpu_sweep(
stage = p["stage"]
if not bool(feas[k]):
st.alive = False; st.reason = f"ступень {s_idx}: разряд не скоммутировался"; continue
surge = stage.switch.max_current_a * SWITCH_SURGE_FACTOR.get(stage.switch.kind, 4.0)
surge = switch_pulse_limit_a(stage.switch)
if float(peak_i[k]) > surge:
st.alive = False
st.reason = f"ступень {s_idx}: пиковый ток {float(peak_i[k]):.0f}А > предела ключа ({surge:.0f}А)"
st.reason = f"ступень {s_idx}: пиковый ток {float(peak_i[k]):.0f}А > импульсного предела ключа ({surge:.0f}А)"
continue
ev = float(exit_v[k])
st.kinetic_delta_j += 0.5 * p["mass"] * (ev**2 - p["v_fire"] ** 2)

View File

@@ -15,11 +15,13 @@ import numpy as np
from gausse.components.database import ComponentDatabase
from gausse.optim.search_space import Genome, SearchBounds, decode, genome_to_dict
from gausse.physics.ac_resistance import dowell_ac_factor
from gausse.physics.constants import switch_pulse_limit_a
from gausse.physics.inductance import air_core_inductance_wheeler, winding_geometry
from gausse.sim.coilgun import CoilgunResult, run_coilgun
from gausse.storage.schema import RunRecord
MODEL_VERSION = "gausse-physics-v1"
MODEL_VERSION = "gausse-physics-v2" # v2: скин/близость (Доуэлл), трение+воздух, паспортные импульсные токи ключей
@dataclass
@@ -115,6 +117,10 @@ def build_detail(config, coilgun_result: CoilgunResult, db: ComponentDatabase, i
l_air_h = air_core_inductance_wheeler(
geometry.mean_radius_m, geometry.coil_length_m, geometry.radial_depth_m, geometry.total_turns
)
_char_omega = 1.0 / math.sqrt(l_air_h * stage.capacitor.capacitance_uf * 1e-6)
r_ac_factor = dowell_ac_factor(
stage.wire.gauge_mm / 1000, wire_od_m, stage.layers, stage.wire.resistivity_ohm_m, _char_omega
)
cap = stage.capacitor
# cap может быть батареей (CapacitorBank) или одиночным конденсатором (CapacitorSpec)
@@ -170,7 +176,10 @@ def build_detail(config, coilgun_result: CoilgunResult, db: ComponentDatabase, i
},
"electrical": {
"charge_voltage_v": stage.charge_voltage_v,
"wire_resistance_ohm": r_wire,
"wire_resistance_dc_ohm": r_wire,
# скин + эффект близости на частоте импульса (Доуэлл)
"wire_ac_resistance_factor": r_ac_factor,
"wire_resistance_ac_ohm": r_wire * r_ac_factor,
"air_inductance_h": l_air_h,
},
}
@@ -181,13 +190,16 @@ def build_detail(config, coilgun_result: CoilgunResult, db: ComponentDatabase, i
peak_current_a = (
float(np.max(np.abs(r.discharge_i))) if r.discharge_i is not None and len(r.discharge_i) else None
)
# диагностика по току: хватает ли батареи и ключа на пиковый ток разряда.
# Это предупреждение, а не жёсткий отказ: в базе — паспортный (не импульсный)
# максимум тока, а одиночный выстрел кратковременный, поэтому честнее
# флажок, чем ложная отбраковка. Параллельные конденсаторы уже дают
# реальный выигрыш через сниженный ESR (меньше потери), см. capacitor_bank.
# Диагностика по току. Жёсткая граница для ключа — ИМПУЛЬСНЫЙ рейтинг
# (ITSM/IDM/ICM из даташита): её превышение уже даёт feasible=False в
# run_stage. Пик выше ПРОДОЛЖИТЕЛЬНОГО рейтинга при одиночном коротком
# импульсе — норма (так работают все coilgun'ы), это ИНФОРМАЦИОННЫЙ
# флаг, а не противоречие с feasible. Для батареи конденсаторов — тоже
# предупреждение (в базе паспортный ripple, не импульсный максимум).
switch_pulse_a = switch_pulse_limit_a(stage.switch)
bank_over = peak_current_a is not None and peak_current_a > cap.max_current_a
switch_over = peak_current_a is not None and peak_current_a > stage.switch.max_current_a
switch_over_cont = peak_current_a is not None and peak_current_a > stage.switch.max_current_a
switch_over_pulse = peak_current_a is not None and peak_current_a > switch_pulse_a
entry["outcome"] = {
"reached": True,
"feasible": r.feasible,
@@ -197,8 +209,12 @@ def build_detail(config, coilgun_result: CoilgunResult, db: ComponentDatabase, i
"t_sensor_s": r.t_sensor_s,
"t_fire_s": r.t_fire_s,
"peak_current_a": peak_current_a,
"switch_pulse_limit_a": switch_pulse_a,
"bank_current_over_limit": bank_over,
"switch_current_over_limit": switch_over,
# пик > продолжительного рейтинга: НОРМА для одиночного импульса
"switch_current_over_continuous_rating": switch_over_cont,
# пик > импульсного (даташит): такой конфиг отбраковывается
"switch_current_over_pulse_limit": switch_over_pulse,
"peak_field_estimate_tesla": r.peak_field_estimate_tesla,
"saturation_warning": r.saturation_warning,
"energy_in_j": r.energy_in_j,

View File

@@ -0,0 +1,61 @@
"""AC-сопротивление обмотки: скин-эффект + эффект близости (метод Доуэлла).
Разряд конденсатора через катушку — импульс с характерной частотой
ω ≈ 1/√(L·C) (сотни Гц — единицы кГц). На этой частоте ток вытесняется
к поверхности провода (скин-эффект), а в многослойной намотке соседние
слои дополнительно вытесняют ток друг у друга (эффект близости). Для
толстого провода и многих слоёв эффективное сопротивление в 1.53 раза
выше DC — этот канал потерь раньше не учитывался и завышал КПД.
Модель — классический Доуэлл (Dowell, 1966) для m-слойной обмотки,
круглый провод сведён к эквивалентной фольге:
Δ = (π/4)^{3/4} · (d/δ) · √η, δ = √(2ρ/(ω·μ₀)) — глубина скина,
η = d_голый/d_изолир — коэффициент заполнения слоя,
F_R = Δ·[φ₁(Δ) + (2(m²1)/3)·φ₂(Δ)],
φ₁ = (sh2Δ+sin2Δ)/(ch2Δcos2Δ), φ₂ = (shΔsinΔ)/(chΔ+cosΔ).
ЧЕСТНО об ограничениях: ток разряда — затухающий импульс, а не синусоида;
берём одну характерную частоту ω=1/√(LC) (первая гармоника импульса), без
разложения в спектр. Доуэлл предполагает плотную рядовую намотку и
пренебрегает кривизной. Точность — десятки процентов, но это несравнимо
честнее, чем DC-сопротивление, занижавшее потери в разы для толстого
провода во многих слоях.
"""
import math
from gausse.physics.constants import MU_0
def skin_depth_m(resistivity_ohm_m: float, omega_rad_s: float) -> float:
"""Глубина скин-слоя δ = √(2ρ/(ω·μ₀)) для немагнитного проводника (Cu/Al)."""
return math.sqrt(2.0 * resistivity_ohm_m / (omega_rad_s * MU_0))
def dowell_ac_factor(
wire_bare_d_m: float,
wire_insulated_d_m: float,
n_layers: int,
resistivity_ohm_m: float,
omega_rad_s: float,
) -> float:
"""F_R = R_AC/R_DC ≥ 1 для m-слойной обмотки на частоте ω (Доуэлл)."""
if omega_rad_s <= 0.0 or n_layers < 1 or resistivity_ohm_m <= 0.0:
return 1.0 # ρ=0 (идеальный провод): R_DC=0, фактор не имеет смысла
delta = skin_depth_m(resistivity_ohm_m, omega_rad_s)
porosity = min(wire_bare_d_m / wire_insulated_d_m, 1.0)
d_eff = (math.pi / 4.0) ** 0.75 * wire_bare_d_m * math.sqrt(porosity)
big_delta = d_eff / delta
m = float(n_layers)
if big_delta < 0.25:
# ряд малых Δ (φ-формулы теряют точность из-за сокращения знаков):
# F_R ≈ 1 + (5m²1)/45 · Δ⁴
return 1.0 + (5.0 * m * m - 1.0) / 45.0 * big_delta**4
two_d = 2.0 * big_delta
phi1 = (math.sinh(two_d) + math.sin(two_d)) / (math.cosh(two_d) - math.cos(two_d))
phi2 = (math.sinh(big_delta) - math.sin(big_delta)) / (math.cosh(big_delta) + math.cos(big_delta))
f_r = big_delta * (phi1 + (2.0 * (m * m - 1.0) / 3.0) * phi2)
return max(f_r, 1.0)

View File

@@ -29,6 +29,15 @@ class StageCircuitParams:
# отражённое сопротивление вихревых токов снаряда; действует × overlap(x)
# (только пока снаряд в катушке). См. physics/losses.py.
r_eddy_coeff_ohm: float = 0.0
# тормозящие силы на снаряд: сухое трение о трубку (константа, Н) и
# аэродинамическое сопротивление (× v², Н·с²/м²). Обе против движения.
retard_const_n: float = 0.0
drag_coeff_n_per_mps2: float = 0.0
def retarding_force_n(params: StageCircuitParams, v: float) -> float:
"""Тормозящая сила (трение + воздух) со знаком ПРОТИВ скорости."""
return float(np.sign(v)) * (params.retard_const_n + params.drag_coeff_n_per_mps2 * v * v)
def effective_resistance(inductance_model: CoilInductanceModel, params: StageCircuitParams, x) -> float:
@@ -51,7 +60,7 @@ def derivatives(
d_q = -current
d_i = (v_c - current * r_eff - dlambda_dx * v) / dlambda_di
force = inductance_model.force_newtons(x, current)
d_v = force / params.mass_kg
d_v = (force - retarding_force_n(params, v)) / params.mass_kg
return [d_q, d_i, v, d_v]

View File

@@ -15,3 +15,29 @@ RESISTIVITY_ALUMINUM_OHM_M = 2.82e-8 # при ~20°C
# Это оценки, а не паспортные surge-рейтинги конкретных деталей.
SWITCH_SURGE_FACTOR = {"SCR": 10.0, "MOSFET": 4.0, "IGBT": 3.0}
CAPACITOR_SURGE_FACTOR = 5.0
def switch_pulse_limit_a(switch) -> float:
"""Импульсный предел тока ключа для одиночного выстрела.
Приоритет — паспортный импульсный рейтинг из даташита (ITSM тиристора,
IDM MOSFET, ICM IGBT), занесённый в базу компонентов. Только если его
нет — консервативный типовой множитель SWITCH_SURGE_FACTOR к
продолжительному току.
"""
pulse = getattr(switch, "pulse_current_a", None)
if pulse:
return float(pulse)
return switch.max_current_a * SWITCH_SURGE_FACTOR.get(switch.kind, 4.0)
G_ACCEL_M_S2 = 9.81
# Сухое трение скольжения сталь↔пластик (ПВХ/оргстекло): типично 0.30.4.
# Труба горизонтальна: нормальная сила = m·g. Оценка, вынесена для уточнения
# по реальному замеру (наклонная плоскость с реальной трубкой и снарядом).
FRICTION_COEFF_SLUG_TUBE = 0.35
AIR_DENSITY_KG_M3 = 1.2 # при ~20°C, уровень моря
# Cd плоского торца цилиндра в свободном потоке ~0.81.2. ЧЕСТНО: в узкой
# трубе есть ещё поршневой эффект (снаряд толкает столб воздуха) — он НЕ
# учтён и занижает потери на высоких скоростях в длинной трубе.
DRAG_COEFF_CYLINDER = 0.85

View File

@@ -50,6 +50,12 @@ def eddy_reflected_resistance_ohm(
def hysteresis_energy_j(slug_radius_m: float, slug_length_m: float, n_cycles: float = 1.0) -> float:
"""Энергия гистерезиса за выстрел (Дж) — вторичный канал, обычно << вихревых."""
"""Энергия гистерезиса за выстрел (Дж) — вторичный канал, обычно << вихревых.
СОЗНАТЕЛЬНО не подключена в динамику: для типичных снарядов это
~0.0050.05 Дж на выстрел (<0.1% энергии банки) — на порядки меньше
вихревых потерь, уже сидящих в контуре через r_eddy. Вкручивать её в
ОДУ значило бы изображать точность, которой у модели нет.
"""
volume = math.pi * slug_radius_m**2 * slug_length_m
return HYSTERESIS_LOSS_J_PER_M3 * volume * n_cycles

View File

@@ -29,7 +29,7 @@ def _sech2(z: np.ndarray) -> np.ndarray:
def make_optical_sensor_event(x_sensor_m: float) -> SensorEvent:
def event(t: float, state: np.ndarray) -> float:
def event(t: float, state: np.ndarray, *_args) -> float:
return state[0] - x_sensor_m
event.terminal = True
@@ -43,7 +43,7 @@ def make_inductive_sensor_event(
threshold_v: float,
width_m: float,
) -> SensorEvent:
def event(t: float, state: np.ndarray) -> float:
def event(t: float, state: np.ndarray, *_args) -> float:
x, v = state
bump = _sech2((x - x_sensor_m) / width_m)
return sensitivity_v_per_mps * v * bump - threshold_v

View File

@@ -28,7 +28,14 @@ from gausse.physics.inductance import (
effective_permeability,
winding_geometry,
)
from gausse.physics.constants import SWITCH_SURGE_FACTOR
from gausse.physics.ac_resistance import dowell_ac_factor
from gausse.physics.constants import (
AIR_DENSITY_KG_M3,
DRAG_COEFF_CYLINDER,
FRICTION_COEFF_SLUG_TUBE,
G_ACCEL_M_S2,
switch_pulse_limit_a,
)
from gausse.physics.force import saturation_scale, solenoid_field_estimate_tesla
from gausse.physics.losses import eddy_reflected_resistance_ohm
from gausse.physics.sensors import (
@@ -114,8 +121,19 @@ def _switch_equivalent_resistance_ohm(
return drop_v / max(i_ref, 1e-6)
def _ballistic_derivatives(t: float, state: np.ndarray) -> list[float]:
return [state[1], 0.0]
def _coast_derivatives(t: float, state: np.ndarray, retard_const_n: float, drag_coeff: float, mass_kg: float) -> list[float]:
"""Подлёт к датчику: не баллистика в вакууме, а с трением о трубку и воздухом."""
v = state[1]
retard = retard_const_n + drag_coeff * v * v
return [v, -math.copysign(retard, v) / mass_kg if v != 0.0 else 0.0]
def _stall_event(t: float, state: np.ndarray, *_args) -> float:
return state[1] - 1e-3 # снаряд практически остановился (трение съело скорость)
_stall_event.terminal = True
_stall_event.direction = -1.0
@dataclass(frozen=True)
@@ -133,23 +151,31 @@ class StagePhysics:
r_eddy_coeff_ohm: float
capacitance_f: float
mu_eff: float
r_wire_dc_ohm: float = 0.0
r_ac_factor: float = 1.0 # R_AC/R_DC обмотки (скин + близость, Доуэлл)
def build_stage_physics(stage: StageConfig, projectile: ProjectileConfig) -> StagePhysics:
wire_od_m = stage.wire.insulation_od_mm / 1000
geometry = winding_geometry(stage.tube_od_m, wire_od_m, stage.turns_per_layer, stage.layers)
r_wire = _wire_resistance_ohm(stage.wire, geometry.total_wire_length_m)
r_wire_dc = _wire_resistance_ohm(stage.wire, geometry.total_wire_length_m)
capacitance_f = stage.capacitor.capacitance_uf * 1e-6
l_air_h = air_core_inductance_wheeler(
geometry.mean_radius_m, geometry.coil_length_m, geometry.radial_depth_m, geometry.total_turns
)
char_omega = 1.0 / math.sqrt(l_air_h * capacitance_f)
# скин + эффект близости: на частоте импульса обмотка сопротивляется
# сильнее, чем по DC (для толстого провода во многих слоях — в разы)
r_ac_factor = dowell_ac_factor(
stage.wire.gauge_mm / 1000, wire_od_m, stage.layers, stage.wire.resistivity_ohm_m, char_omega
)
r_wire = r_wire_dc * r_ac_factor
r_switch = _switch_equivalent_resistance_ohm(
stage.switch, capacitance_f, l_air_h, stage.charge_voltage_v
)
r_total_ohm = r_wire + r_switch + stage.capacitor.esr_ohm
demag = demagnetizing_factor_prolate(projectile.aspect_ratio)
mu_eff = effective_permeability(projectile.material.mu_r, demag)
char_omega = 1.0 / math.sqrt(l_air_h * capacitance_f)
r_eddy_coeff = eddy_reflected_resistance_ohm(
slug_radius_m=projectile.diameter_m / 2,
slug_length_m=projectile.length_m,
@@ -168,11 +194,15 @@ def build_stage_physics(stage: StageConfig, projectile: ProjectileConfig) -> Sta
total_turns=geometry.total_turns,
b_sat_tesla=projectile.material.b_sat_tesla,
)
mass_kg = projectile.mass_kg
frontal_area_m2 = math.pi * (projectile.diameter_m / 2) ** 2
circuit_params = StageCircuitParams(
capacitance_f=capacitance_f,
r_total_ohm=r_total_ohm,
mass_kg=projectile.mass_kg,
mass_kg=mass_kg,
r_eddy_coeff_ohm=r_eddy_coeff,
retard_const_n=FRICTION_COEFF_SLUG_TUBE * mass_kg * G_ACCEL_M_S2,
drag_coeff_n_per_mps2=0.5 * AIR_DENSITY_KG_M3 * DRAG_COEFF_CYLINDER * frontal_area_m2,
)
return StagePhysics(
inductance_model=inductance_model,
@@ -182,6 +212,8 @@ def build_stage_physics(stage: StageConfig, projectile: ProjectileConfig) -> Sta
r_eddy_coeff_ohm=r_eddy_coeff,
capacitance_f=capacitance_f,
mu_eff=mu_eff,
r_wire_dc_ohm=r_wire_dc,
r_ac_factor=r_ac_factor,
)
@@ -223,14 +255,21 @@ def run_stage(
else:
sensor_event = make_optical_sensor_event(x_sensor_m)
cp = phys.circuit_params
flight_sol = solve_ivp(
_ballistic_derivatives,
_coast_derivatives,
(0.0, stage.sensor_time_budget_s),
[entry_x_m, entry_v_mps],
events=sensor_event,
args=(cp.retard_const_n, cp.drag_coeff_n_per_mps2, mass_kg),
events=(sensor_event, _stall_event),
rtol=1e-8,
atol=1e-10,
)
if len(flight_sol.t_events[1]) > 0 and len(flight_sol.t_events[0]) == 0:
return StageResult(
feasible=False,
reason="снаряд остановлен трением о трубку, не долетев до датчика",
)
if len(flight_sol.t_events[0]) == 0:
return StageResult(
feasible=False,
@@ -284,7 +323,12 @@ def run_stage(
# её гасит обратный диод (freewheel), учитываем как потери, чтобы энергия
# сходилась. Считаем ТОЧНО (с насыщением), а не ½LI².
residual_inductor_j = float(inductance_model.magnetic_energy(x_final, i_final))
energy_dissipated_j = resistive_diss + residual_inductor_j
# работа трения о трубку и воздуха ЗА РАЗРЯД: P = (F_тр + F_возд)·|v|
v_during = discharge_sol.y[3]
friction_diss = float(
np.trapezoid((cp.retard_const_n + cp.drag_coeff_n_per_mps2 * v_during**2) * np.abs(v_during), discharge_sol.t)
)
energy_dissipated_j = resistive_diss + residual_inductor_j + friction_diss
kinetic_before_j = 0.5 * mass_kg * v_fire_mps**2
kinetic_after_j = 0.5 * mass_kg * v_final**2
@@ -310,13 +354,13 @@ def run_stage(
# Батарея конденсаторов — только предупреждение (см. build_detail):
# электролиты переносят высокий импульсный ток, а в базе — паспортный
# непрерывный/ripple рейтинг, поэтому жёстко браковать по нему нечестно.
switch_surge_a = stage.switch.max_current_a * SWITCH_SURGE_FACTOR.get(stage.switch.kind, 4.0)
switch_surge_a = switch_pulse_limit_a(stage.switch)
if peak_current_a > switch_surge_a:
return StageResult(
feasible=False,
reason=(
f"пиковый ток {peak_current_a:.0f} А превышает импульсный предел ключа "
f"{stage.switch.part_number} ({switch_surge_a:.0f} А = {stage.switch.max_current_a:.0f} А × surge)"
f"{stage.switch.part_number} ({switch_surge_a:.0f} А, даташит/surge)"
),
t_sensor_s=t_sensor_s,
t_fire_s=t_sensor_s + fire_delay_s,

View File

@@ -17,7 +17,9 @@ def test_evaluate_feasible_genome_has_fitness_equal_efficiency():
if result.feasible:
found_feasible = True
assert result.fitness == result.efficiency
assert 0.0 <= result.efficiency <= 1.0
# КПД <= 1 всегда; может быть СЛЕГКА отрицательным: трение о трубку
# за время разряда может съесть больше, чем слабая катушка добавила
assert -0.01 <= result.efficiency <= 1.0
assert result.cost_rub > 0
break
assert found_feasible, "ни один из 50 случайных геномов не оказался реализуемым"

147
tests/test_realism_v2.py Normal file
View File

@@ -0,0 +1,147 @@
"""Физреализм 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]

View File

@@ -61,13 +61,14 @@ def test_stage_is_feasible_with_realistic_components():
assert result.feasible, result.reason
def test_energy_conserved_exactly_when_lossless():
def test_energy_conserved_exactly_when_electrically_lossless():
lossless_capacitor = CapacitorSpec(
part_number="c-lossless", capacitance_uf=100.0, voltage_v=400.0, esr_ohm=0.0,
max_current_a=300.0, price=300.0, source="test",
)
# снаряд без вихревых потерь: огромное удельное сопротивление -> R_eddy≈0,
# чтобы «без потерь» действительно означало отсутствие всех каналов диссипации
# снаряд без вихревых потерь: огромное удельное сопротивление -> R_eddy≈0.
# Электрических потерь нет, но МЕХАНИЧЕСКИЕ (трение о трубку + воздух)
# есть всегда — они и составляют весь energy_dissipated_j.
no_eddy_steel = ProjectileMaterialSpec(
name="steel-no-eddy", density_kg_m3=7850.0, mu_r=200.0, b_sat_tesla=1.8,
price_per_kg=100.0, source="test", resistivity_ohm_m=1e12,
@@ -78,8 +79,9 @@ def test_energy_conserved_exactly_when_lossless():
capacitor=lossless_capacitor, projectile=no_eddy_projectile,
)
assert result.feasible, result.reason
assert result.energy_dissipated_j == pytest.approx(0.0, abs=1e-9)
balance = result.energy_remaining_cap_j + result.kinetic_energy_delta_j
# только трение/воздух: заметно меньше 1% энергии банки за миллисекунды разряда
assert 0.0 <= result.energy_dissipated_j < 0.01 * result.energy_in_j
balance = result.energy_remaining_cap_j + result.energy_dissipated_j + result.kinetic_energy_delta_j
assert balance == pytest.approx(result.energy_in_j, rel=1e-4)

View File

@@ -39,7 +39,9 @@ def test_run_detail_has_full_physics_and_distances(tmp_path):
assert "inter_stage_gaps_m" in detail
assert "coil_center_positions_m" in detail
stage0 = detail["stages"][0]
assert "electrical" in stage0 and "wire_resistance_ohm" in stage0["electrical"]
assert "electrical" in stage0 and "wire_resistance_ac_ohm" in stage0["electrical"]
# v2: скин/близость — AC-сопротивление обмотки выше DC
assert stage0["electrical"]["wire_ac_resistance_factor"] >= 1.0
assert "winding" in stage0 and "total_turns" in stage0["winding"]
assert "sensor_to_coil_distance_m" in stage0