Physical tube geometry, 200g projectiles, and single-pulse discharge

Three fixes from user feedback:

1. Single-pulse discharge (user spotted 3 current humps for one stage):
   a thyristor fires ONCE per shot -- you can't recharge the cap in
   microseconds. The discharge now terminates at the first current zero
   OR first local minimum (where the slug starts pumping current back),
   whichever comes first. Residual coil energy at cutoff is accounted as
   freewheel-diode dissipation using the exact saturating magnetic energy
   integral, so energy still balances. Verified: 3 humps -> 1.

2. Real tube geometry: the genome now carries tube INNER diameter (bore,
   the projectile flies through) and wall thickness; outer diameter =
   inner + 2*wall = the coil's inner diameter (which drives the field).
   The projectile must fit the bore (diameter < inner - clearance).

3. Projectiles up to 200 g: diameter to 28mm, length to 150mm, with mass
   capped at 200g (length clamped by density).

Detail/BOM now report inner/outer/wall tube diameters, projectile mass in
grams, and whether it fits the bore. Same single-pulse + eddy physics
mirrored into the GPU batch integrator. Evolutionary polish updated for
the new tube params. All test groups pass.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
This commit is contained in:
jze9
2026-07-07 15:33:47 +05:00
parent 4d16825f7d
commit 35cbbb04c6
8 changed files with 180 additions and 51 deletions

View File

@@ -70,13 +70,25 @@ def _ln_cosh(xp, z):
return az + xp.log1p(xp.exp(-2 * az)) - xp.log(xp.asarray(2.0)) return az + xp.log1p(xp.exp(-2 * az)) - xp.log(xp.asarray(2.0))
def _derivatives(xp, q, i, x, v, p: BatchDischargeParams): def _overlap(xp, x, p: BatchDischargeParams):
half_span = (p.coil_length_m + p.slug_length_m) / 2 half_span = (p.coil_length_m + p.slug_length_m) / 2
w = p.smoothing_width_m w = p.smoothing_width_m
s1 = 0.5 * (1 + xp.tanh((x + half_span) / w / 2)) s1 = 0.5 * (1 + xp.tanh((x + half_span) / w / 2))
s2 = 0.5 * (1 + xp.tanh((half_span - x) / w / 2)) s2 = 0.5 * (1 + xp.tanh((half_span - x) / w / 2))
overlap = s1 * s2 return s1 * s2, (s1 * s2 / w) * (s2 - s1)
d_overlap = (s1 * s2 / w) * (s2 - s1)
def _magnetic_energy(xp, x, i, p: BatchDischargeParams):
"""W = ½L_air·I² + L_iron·overlap·(I·g G) — точная (с насыщением) энергия поля."""
overlap, _ = _overlap(xp, x, p)
z = i / p.i_sat_a
g = p.i_sat_a * xp.tanh(z)
g_integral = p.i_sat_a**2 * _ln_cosh(xp, z)
return 0.5 * p.l_air_h * i**2 + p.l_iron_coeff * overlap * (i * g - g_integral)
def _derivatives(xp, q, i, x, v, p: BatchDischargeParams):
overlap, d_overlap = _overlap(xp, x, p)
z = i / p.i_sat_a z = i / p.i_sat_a
tanhz = xp.tanh(z) tanhz = xp.tanh(z)
@@ -118,6 +130,7 @@ def integrate_batch_discharge(
v = xp.array(v0, dtype=xp.float64) v = xp.array(v0, dtype=xp.float64)
done = xp.zeros(n, dtype=bool) done = xp.zeros(n, dtype=bool)
past_peak = xp.zeros(n, dtype=bool) # ток уже прошёл пик и начал спадать
exit_v = xp.array(v0, dtype=xp.float64) exit_v = xp.array(v0, dtype=xp.float64)
exit_x = xp.array(x0, dtype=xp.float64) exit_x = xp.array(x0, dtype=xp.float64)
exit_q = xp.array(q0, dtype=xp.float64) exit_q = xp.array(q0, dtype=xp.float64)
@@ -142,19 +155,33 @@ def integrate_batch_discharge(
x_new = x_old + dt / 6 * (k1[2] + 2 * k2[2] + 2 * k3[2] + k4[2]) x_new = x_old + dt / 6 * (k1[2] + 2 * k2[2] + 2 * k3[2] + k4[2])
v_new = v_old + dt / 6 * (k1[3] + 2 * k2[3] + 2 * k3[3] + k4[3]) v_new = v_old + dt / 6 * (k1[3] + 2 * k2[3] + 2 * k3[3] + k4[3])
# потери R dt (трапеция) и пиковый ток — только для активных # потери I²·R_eff (с вихревыми × overlap) и пиковый ток — только для активных
step_diss = 0.5 * (i_old**2 + i_new**2) * params.r_total_ohm * dt overlap_old, _ = _overlap(xp, x_old, params)
overlap_new, _ = _overlap(xp, x_new, params)
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
energy_diss = energy_diss + xp.where(active, step_diss, 0.0) 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)) 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))
# нуль тока: ток был положительным и стал <=0 -> разряд закончился. # ОДИН импульс: обрыв на нуле тока ИЛИ на первом локальном минимуме
# exit-значения берём интерполяцией СТАРОГО и НОВОГО состояния в точке I=0. # (снаряд начал подкачивать ток обратно). Что раньше.
crossed = active & (i_old > 0) & (i_new <= 0) crossed = active & (i_old > 0) & (i_new <= 0)
local_min = active & past_peak & (i_new > i_old) & ~crossed
cut = crossed | local_min
# zero-crossing: интерполяция к I=0; local_min: берём состояние минимума (old)
frac = xp.where(crossed, i_old / (i_old - i_new + 1e-30), 0.0) frac = xp.where(crossed, i_old / (i_old - i_new + 1e-30), 0.0)
exit_v = xp.where(crossed, _interp(v_old, v_new, frac), exit_v) exit_v = xp.where(crossed, _interp(v_old, v_new, frac), exit_v)
exit_x = xp.where(crossed, _interp(x_old, x_new, frac), exit_x) exit_x = xp.where(crossed, _interp(x_old, x_new, frac), exit_x)
exit_q = xp.where(crossed, _interp(q_old, q_new, frac), exit_q) exit_q = xp.where(crossed, _interp(q_old, q_new, frac), exit_q)
done = done | crossed exit_v = xp.where(local_min, v_old, exit_v)
exit_x = xp.where(local_min, x_old, exit_x)
exit_q = xp.where(local_min, q_old, exit_q)
# остаточная энергия катушки при обрыве на минимуме -> в потери (freewheel)
residual = _magnetic_energy(xp, x_old, i_old, params)
energy_diss = energy_diss + xp.where(local_min, residual, 0.0)
done = done | cut
# продвигаем только ещё активные конфигурации # продвигаем только ещё активные конфигурации
q = xp.where(active, q_new, q_old) q = xp.where(active, q_new, q_old)

View File

@@ -39,7 +39,7 @@ def _tournament_select(
def _continuous_vector(genome: Genome) -> list[float]: def _continuous_vector(genome: Genome) -> list[float]:
vec = [genome.tube_od_m] vec = [genome.tube_inner_d_m, genome.tube_wall_m]
for stage in genome.stages: for stage in genome.stages:
vec.append(stage.sensor_to_coil_distance_m) vec.append(stage.sensor_to_coil_distance_m)
vec.append(stage.charge_voltage_fraction) vec.append(stage.charge_voltage_fraction)
@@ -52,7 +52,9 @@ def _continuous_vector(genome: Genome) -> list[float]:
def _apply_continuous_vector(template: Genome, vec: list[float]) -> Genome: def _apply_continuous_vector(template: Genome, vec: list[float]) -> Genome:
genome = copy.deepcopy(template) genome = copy.deepcopy(template)
idx = 0 idx = 0
genome.tube_od_m = vec[idx] genome.tube_inner_d_m = vec[idx]
idx += 1
genome.tube_wall_m = vec[idx]
idx += 1 idx += 1
for stage in genome.stages: for stage in genome.stages:
stage.sensor_to_coil_distance_m = vec[idx] stage.sensor_to_coil_distance_m = vec[idx]

View File

@@ -88,7 +88,7 @@ def decoded_summary(genome: Genome, db: ComponentDatabase) -> dict:
} }
def build_detail(config, coilgun_result: CoilgunResult, db: ComponentDatabase, initial_x_m: float, initial_v_mps: float) -> dict: def build_detail(config, coilgun_result: CoilgunResult, db: ComponentDatabase, initial_x_m: float, initial_v_mps: float, tube_inner_d_m: float | None = None, tube_wall_m: float | None = None) -> dict:
"""Максимально подробная запись: параметры + вычисленная физика + результат каждой ступени. """Максимально подробная запись: параметры + вычисленная физика + результат каждой ступени.
Именно это пишется в SQLite (decoded_summary_json), чтобы по базе можно Именно это пишется в SQLite (decoded_summary_json), чтобы по базе можно
@@ -197,16 +197,24 @@ def build_detail(config, coilgun_result: CoilgunResult, db: ComponentDatabase, i
entry["outcome"] = {"reached": False} entry["outcome"] = {"reached": False}
stages_detail.append(entry) stages_detail.append(entry)
tube_od_m = config.stages[0].tube_od_m if config.stages else None
return { return {
"tube_od_m": config.stages[0].tube_od_m if config.stages else None, "tube": {
"inner_diameter_m": tube_inner_d_m, # бор — снаряд летит внутри
"wall_thickness_m": tube_wall_m,
"outer_diameter_m": tube_od_m, # = внутренний диаметр катушки
},
"tube_od_m": tube_od_m, # для обратной совместимости
"tube_length_m": (coil_centers_m[-1] + 0.05) if coil_centers_m else None, "tube_length_m": (coil_centers_m[-1] + 0.05) if coil_centers_m else None,
"projectile": { "projectile": {
"material": projectile.material.name, "material": projectile.material.name,
"diameter_m": projectile.diameter_m, "diameter_m": projectile.diameter_m,
"length_m": projectile.length_m, "length_m": projectile.length_m,
"mass_kg": projectile.mass_kg, "mass_kg": projectile.mass_kg,
"mass_g": projectile.mass_kg * 1000,
"mu_r": projectile.material.mu_r, "mu_r": projectile.material.mu_r,
"b_sat_tesla": projectile.material.b_sat_tesla, "b_sat_tesla": projectile.material.b_sat_tesla,
"fits_in_bore": (tube_inner_d_m is None) or (projectile.diameter_m < tube_inner_d_m),
}, },
"launch": {"initial_x_m": initial_x_m, "initial_v_mps": initial_v_mps}, "launch": {"initial_x_m": initial_x_m, "initial_v_mps": initial_v_mps},
"n_stages": len(config.stages), "n_stages": len(config.stages),
@@ -220,7 +228,7 @@ def evaluate(genome: Genome, db: ComponentDatabase, bounds: SearchBounds) -> Eva
config, initial_x_m, initial_v_mps = decode(genome, db, bounds) config, initial_x_m, initial_v_mps = decode(genome, db, bounds)
cost_rub = compute_cost_rub(config, db) cost_rub = compute_cost_rub(config, db)
result = run_coilgun(config) result = run_coilgun(config)
detail = build_detail(config, result, db, initial_x_m, initial_v_mps) detail = build_detail(config, result, db, initial_x_m, initial_v_mps, genome.tube_inner_d_m, genome.tube_wall_m)
if not result.feasible: if not result.feasible:
n_stages = len(config.stages) n_stages = len(config.stages)

View File

@@ -16,7 +16,8 @@ from gausse.components.schema import CapacitorBank
from gausse.sim.coilgun import CoilgunConfig from gausse.sim.coilgun import CoilgunConfig
from gausse.sim.stage import ProjectileConfig, StageConfig from gausse.sim.stage import ProjectileConfig, StageConfig
TUBE_WALL_CLEARANCE_M = 0.001 # зазор между снарядом и внутренней стенкой трубы (бором), чтобы снаряд летел свободно
PROJECTILE_BORE_CLEARANCE_M = 0.0005
@dataclass(frozen=True) @dataclass(frozen=True)
@@ -31,12 +32,18 @@ class SearchBounds:
sensor_to_coil_distance_m_max: float = 0.05 sensor_to_coil_distance_m_max: float = 0.05
inter_stage_gap_m_min: float = 0.01 inter_stage_gap_m_min: float = 0.01
inter_stage_gap_m_max: float = 0.10 inter_stage_gap_m_max: float = 0.10
tube_od_m_min: float = 0.009 # ТРУБА: внутренний диаметр (бор — снаряд летит внутри) и толщина стенки.
tube_od_m_max: float = 0.02 # Внешний диаметр = внутренний + 2×стенка = внутренний диаметр КАТУШКИ.
tube_inner_d_m_min: float = 0.005
tube_inner_d_m_max: float = 0.030
tube_wall_m_min: float = 0.001
tube_wall_m_max: float = 0.003
# СНАРЯД: диаметр обязан быть < внутр. диаметра трубы; масса до 200 г.
projectile_diameter_m_min: float = 0.004 projectile_diameter_m_min: float = 0.004
projectile_diameter_m_max: float = 0.008 projectile_diameter_m_max: float = 0.028
projectile_length_m_min: float = 0.01 projectile_length_m_min: float = 0.01
projectile_length_m_max: float = 0.03 projectile_length_m_max: float = 0.15
projectile_mass_max_kg: float = 0.2
charge_voltage_fraction_min: float = 0.5 charge_voltage_fraction_min: float = 0.5
charge_voltage_fraction_max: float = 1.0 charge_voltage_fraction_max: float = 1.0
initial_launch_velocity_mps: float = 3.0 initial_launch_velocity_mps: float = 3.0
@@ -48,6 +55,12 @@ class SearchBounds:
cap_parallel_max: int = 8 cap_parallel_max: int = 8
def _max_length_for_mass(diameter_m: float, density_kg_m3: float, max_mass_kg: float) -> float:
"""Максимальная длина снаряда, при которой масса не превышает max_mass_kg."""
area = math.pi * (diameter_m / 2) ** 2
return max_mass_kg / (density_kg_m3 * area)
@dataclass @dataclass
class StageGene: class StageGene:
wire_idx: int wire_idx: int
@@ -71,11 +84,17 @@ class ProjectileGene:
@dataclass @dataclass
class Genome: class Genome:
tube_od_m: float tube_inner_d_m: float # бор трубы (снаряд летит внутри)
tube_wall_m: float # толщина стенки трубы
stages: list[StageGene] stages: list[StageGene]
inter_stage_gaps_m: list[float] inter_stage_gaps_m: list[float]
projectile: ProjectileGene projectile: ProjectileGene
@property
def tube_od_m(self) -> float:
"""Внешний диаметр трубы = внутренний диаметр катушки."""
return self.tube_inner_d_m + 2 * self.tube_wall_m
def _clip(value: float, lo: float, hi: float) -> float: def _clip(value: float, lo: float, hi: float) -> float:
return max(lo, min(hi, value)) return max(lo, min(hi, value))
@@ -101,36 +120,47 @@ def sample_stage_gene(db: ComponentDatabase, bounds: SearchBounds, rng: random.R
def sample_genome(db: ComponentDatabase, bounds: SearchBounds, rng: random.Random) -> Genome: def sample_genome(db: ComponentDatabase, bounds: SearchBounds, rng: random.Random) -> Genome:
tube_od_m = rng.uniform(bounds.tube_od_m_min, bounds.tube_od_m_max) tube_inner_d_m = rng.uniform(bounds.tube_inner_d_m_min, bounds.tube_inner_d_m_max)
max_diameter = min(bounds.projectile_diameter_m_max, tube_od_m - TUBE_WALL_CLEARANCE_M) tube_wall_m = rng.uniform(bounds.tube_wall_m_min, bounds.tube_wall_m_max)
# диаметр снаряда < внутр. диаметра трубы минус зазор
max_diameter = min(bounds.projectile_diameter_m_max, tube_inner_d_m - PROJECTILE_BORE_CLEARANCE_M)
diameter_m = rng.uniform(bounds.projectile_diameter_m_min, max(max_diameter, bounds.projectile_diameter_m_min)) diameter_m = rng.uniform(bounds.projectile_diameter_m_min, max(max_diameter, bounds.projectile_diameter_m_min))
material_idx = rng.randrange(len(db.projectile_materials))
density = db.projectile_materials[material_idx].density_kg_m3
# длина ограничена и границами, и максимальной массой 200 г
len_cap = min(bounds.projectile_length_m_max, _max_length_for_mass(diameter_m, density, bounds.projectile_mass_max_kg))
length_m = rng.uniform(bounds.projectile_length_m_min, max(len_cap, bounds.projectile_length_m_min))
n_stages = rng.randint(bounds.min_stages, bounds.max_stages) n_stages = rng.randint(bounds.min_stages, bounds.max_stages)
stages = [sample_stage_gene(db, bounds, rng) for _ in range(n_stages)] stages = [sample_stage_gene(db, bounds, rng) for _ in range(n_stages)]
gaps = [ gaps = [
rng.uniform(bounds.inter_stage_gap_m_min, bounds.inter_stage_gap_m_max) rng.uniform(bounds.inter_stage_gap_m_min, bounds.inter_stage_gap_m_max)
for _ in range(n_stages - 1) for _ in range(n_stages - 1)
] ]
projectile = ProjectileGene( projectile = ProjectileGene(material_idx=material_idx, diameter_m=diameter_m, length_m=length_m)
material_idx=rng.randrange(len(db.projectile_materials)), return Genome(
diameter_m=diameter_m, tube_inner_d_m=tube_inner_d_m, tube_wall_m=tube_wall_m,
length_m=rng.uniform(bounds.projectile_length_m_min, bounds.projectile_length_m_max), stages=stages, inter_stage_gaps_m=gaps, projectile=projectile,
) )
return Genome(tube_od_m=tube_od_m, stages=stages, inter_stage_gaps_m=gaps, projectile=projectile)
def repair(genome: Genome, db: ComponentDatabase, bounds: SearchBounds) -> Genome: def repair(genome: Genome, db: ComponentDatabase, bounds: SearchBounds) -> Genome:
"""Приводит геном в границы после мутации/скрещивания (клэмп, не отбраковка).""" """Приводит геном в границы после мутации/скрещивания (клэмп, не отбраковка)."""
genome.tube_od_m = _clip(genome.tube_od_m, bounds.tube_od_m_min, bounds.tube_od_m_max) genome.tube_inner_d_m = _clip(genome.tube_inner_d_m, bounds.tube_inner_d_m_min, bounds.tube_inner_d_m_max)
max_diameter = max(genome.tube_od_m - TUBE_WALL_CLEARANCE_M, bounds.projectile_diameter_m_min) genome.tube_wall_m = _clip(genome.tube_wall_m, bounds.tube_wall_m_min, bounds.tube_wall_m_max)
genome.projectile.diameter_m = _clip(
genome.projectile.diameter_m, bounds.projectile_diameter_m_min, max_diameter
)
genome.projectile.length_m = _clip(
genome.projectile.length_m, bounds.projectile_length_m_min, bounds.projectile_length_m_max
)
genome.projectile.material_idx = genome.projectile.material_idx % len(db.projectile_materials) genome.projectile.material_idx = genome.projectile.material_idx % len(db.projectile_materials)
# снаряд обязан влезать в бор трубы
max_diameter = max(genome.tube_inner_d_m - PROJECTILE_BORE_CLEARANCE_M, bounds.projectile_diameter_m_min)
max_diameter = min(max_diameter, bounds.projectile_diameter_m_max)
genome.projectile.diameter_m = _clip(genome.projectile.diameter_m, bounds.projectile_diameter_m_min, max_diameter)
# длина в границах И под массой ≤ 200 г
density = db.projectile_materials[genome.projectile.material_idx].density_kg_m3
len_cap = min(bounds.projectile_length_m_max, _max_length_for_mass(genome.projectile.diameter_m, density, bounds.projectile_mass_max_kg))
genome.projectile.length_m = _clip(genome.projectile.length_m, bounds.projectile_length_m_min, max(len_cap, bounds.projectile_length_m_min))
for stage in genome.stages: for stage in genome.stages:
stage.wire_idx %= len(db.wires) stage.wire_idx %= len(db.wires)
stage.capacitor_idx %= len(db.capacitors) stage.capacitor_idx %= len(db.capacitors)
@@ -187,7 +217,7 @@ def decode(genome: Genome, db: ComponentDatabase, bounds: SearchBounds) -> tuple
capacitor=capacitor, capacitor=capacitor,
switch=switch, switch=switch,
sensor=sensor, sensor=sensor,
tube_od_m=genome.tube_od_m, tube_od_m=genome.tube_od_m, # внешний диаметр трубы = внутр. диаметр катушки
turns_per_layer=gene.turns_per_layer, turns_per_layer=gene.turns_per_layer,
layers=gene.layers, layers=gene.layers,
sensor_to_coil_distance_m=gene.sensor_to_coil_distance_m, sensor_to_coil_distance_m=gene.sensor_to_coil_distance_m,
@@ -209,7 +239,9 @@ def mutate(genome: Genome, db: ComponentDatabase, bounds: SearchBounds, rng: ran
child = copy.deepcopy(genome) child = copy.deepcopy(genome)
if rng.random() < rate: if rng.random() < rate:
child.tube_od_m = rng.uniform(bounds.tube_od_m_min, bounds.tube_od_m_max) child.tube_inner_d_m = rng.uniform(bounds.tube_inner_d_m_min, bounds.tube_inner_d_m_max)
if rng.random() < rate:
child.tube_wall_m = rng.uniform(bounds.tube_wall_m_min, bounds.tube_wall_m_max)
for stage in child.stages: for stage in child.stages:
if rng.random() < rate: if rng.random() < rate:
@@ -295,7 +327,8 @@ def crossover(parent_a: Genome, parent_b: Genome, rng: random.Random) -> Genome:
tube_source = parent_a if rng.random() < 0.5 else parent_b tube_source = parent_a if rng.random() < 0.5 else parent_b
return Genome( return Genome(
tube_od_m=tube_source.tube_od_m, tube_inner_d_m=tube_source.tube_inner_d_m,
tube_wall_m=tube_source.tube_wall_m,
stages=stages, stages=stages,
inter_stage_gaps_m=gaps, inter_stage_gaps_m=gaps,
projectile=copy.deepcopy(projectile_source.projectile), projectile=copy.deepcopy(projectile_source.projectile),
@@ -304,7 +337,8 @@ def crossover(parent_a: Genome, parent_b: Genome, rng: random.Random) -> Genome:
def genome_to_dict(genome: Genome) -> dict: def genome_to_dict(genome: Genome) -> dict:
return { return {
"tube_od_m": genome.tube_od_m, "tube_inner_d_m": genome.tube_inner_d_m,
"tube_wall_m": genome.tube_wall_m,
"stages": [vars(s) for s in genome.stages], "stages": [vars(s) for s in genome.stages],
"inter_stage_gaps_m": genome.inter_stage_gaps_m, "inter_stage_gaps_m": genome.inter_stage_gaps_m,
"projectile": vars(genome.projectile), "projectile": vars(genome.projectile),
@@ -312,8 +346,16 @@ def genome_to_dict(genome: Genome) -> dict:
def genome_from_dict(data: dict) -> Genome: def genome_from_dict(data: dict) -> Genome:
# обратная совместимость: старые геномы в БД несли tube_od_m без стенки
if "tube_inner_d_m" in data:
tube_inner_d_m = data["tube_inner_d_m"]
tube_wall_m = data["tube_wall_m"]
else:
tube_wall_m = 0.0015
tube_inner_d_m = max(data.get("tube_od_m", 0.012) - 2 * tube_wall_m, 0.005)
return Genome( return Genome(
tube_od_m=data["tube_od_m"], tube_inner_d_m=tube_inner_d_m,
tube_wall_m=tube_wall_m,
stages=[StageGene(**s) for s in data["stages"]], stages=[StageGene(**s) for s in data["stages"]],
inter_stage_gaps_m=list(data["inter_stage_gaps_m"]), inter_stage_gaps_m=list(data["inter_stage_gaps_m"]),
projectile=ProjectileGene(**data["projectile"]), projectile=ProjectileGene(**data["projectile"]),

View File

@@ -61,3 +61,24 @@ def zero_current_crossing_event(t: float, state: np.ndarray, *_args) -> float:
zero_current_crossing_event.terminal = True zero_current_crossing_event.terminal = True
zero_current_crossing_event.direction = -1.0 zero_current_crossing_event.direction = -1.0
def current_local_min_event(t: float, state: np.ndarray, inductance_model, params) -> float:
"""Первый локальный МИНИМУМ тока после пика — там снаряд начинает
«подкачивать» ток обратно (генераторный режим за центром катушки).
Тиристор проводит ОДИН импульс: как только ток, спадая, перестаёт падать
и пошёл бы вверх, ключ запирается (в реальности так и есть — за один
выстрел катушка стреляет один раз, конденсатор не перезаряжается за мкс).
Событие = dI/dt, ловим переход снизу вверх (direction=+1).
"""
q, current, x, v = state
dlambda_di = inductance_model.dlambda_di(x, current)
dlambda_dx = inductance_model.dlambda_dx(x, current)
r_eff = effective_resistance(inductance_model, params, x)
v_c = q / params.capacitance_f
return (v_c - current * r_eff - dlambda_dx * v) / dlambda_di
current_local_min_event.terminal = True
current_local_min_event.direction = 1.0

View File

@@ -175,6 +175,16 @@ class CoilInductanceModel:
"""Сила из coenergy: F = ∂/∂x ∫₀ᴵ λ dI' = L_iron·overlap'(x)·G(I).""" """Сила из coenergy: F = ∂/∂x ∫₀ᴵ λ dI' = L_iron·overlap'(x)·G(I)."""
return self.l_iron_coeff * self.d_overlap_dx(x) * self._g_integral(i) return self.l_iron_coeff * self.d_overlap_dx(x) * self._g_integral(i)
def magnetic_energy(self, x: np.ndarray, i: np.ndarray) -> np.ndarray:
"""Энергия магнитного поля W = ∫₀ᴵ i·(∂λ/∂i) di (с учётом насыщения).
= ½·L_air·I² + L_iron·overlap(x)·(I·g(I) G(I)). Нужна для точного
учёта остаточной энергии катушки при обрыве разряда (одиночный импульс).
"""
return 0.5 * self.l_air_h * i**2 + self.l_iron_coeff * self.overlap_fraction(x) * (
i * self._g(i) - self._g_integral(i)
)
def effective_field_current_a(self, i: float) -> float: def effective_field_current_a(self, i: float) -> float:
"""'Насыщенный' эффективный ток для оценки поля (поле ~ g(I), кап B_sat).""" """'Насыщенный' эффективный ток для оценки поля (поле ~ g(I), кап B_sat)."""
return float(self._g(np.asarray(float(i)))) return float(self._g(np.asarray(float(i))))

View File

@@ -214,33 +214,45 @@ def run_stage(
r_eddy_coeff_ohm=r_eddy_coeff, r_eddy_coeff_ohm=r_eddy_coeff,
) )
q0 = capacitance_f * stage.charge_voltage_v q0 = capacitance_f * stage.charge_voltage_v
# Тиристор проводит ОДИН импульс: разряд обрывается на первом возврате тока
# к нулю ИЛИ на первом локальном минимуме тока (снаряд начал подкачивать ток
# обратно) — что раньше. Иначе модель гоняла бы катушку несколько раз за
# выстрел, что физически невозможно (конденсатор не перезарядить за мкс).
discharge_sol = solve_ivp( discharge_sol = solve_ivp(
circuit.derivatives, circuit.derivatives,
(0.0, stage.discharge_time_budget_s), (0.0, stage.discharge_time_budget_s),
[q0, 0.0, x_fire_m, v_fire_mps], [q0, 0.0, x_fire_m, v_fire_mps],
args=(inductance_model, circuit_params), args=(inductance_model, circuit_params),
events=circuit.zero_current_crossing_event, events=(circuit.zero_current_crossing_event, circuit.current_local_min_event),
rtol=1e-8, rtol=1e-8,
atol=1e-11, atol=1e-11,
) )
if len(discharge_sol.t_events[0]) == 0: # какое из событий оборвало разряд первым
term_candidates = []
for ev_idx in range(2):
if len(discharge_sol.t_events[ev_idx]) > 0:
term_candidates.append((discharge_sol.t_events[ev_idx][0], discharge_sol.y_events[ev_idx][0]))
if not term_candidates:
return StageResult( return StageResult(
feasible=False, feasible=False,
reason="разряд не скоммутировался (ток не вернулся к нулю) в пределах временного бюджета", reason="разряд не скоммутировался (ток не вернулся к нулю) в пределах временного бюджета",
t_sensor_s=t_sensor_s, t_sensor_s=t_sensor_s,
t_fire_s=t_sensor_s + fire_delay_s, t_fire_s=t_sensor_s + fire_delay_s,
) )
term_candidates.sort(key=lambda p: p[0])
q_final, i_final, x_final, v_final = discharge_sol.y_events[0][0] q_final, i_final, x_final, v_final = term_candidates[0][1]
energy_in_j = 0.5 * capacitance_f * stage.charge_voltage_v**2 energy_in_j = 0.5 * capacitance_f * stage.charge_voltage_v**2
energy_remaining_cap_j = q_final**2 / (2 * capacitance_f) energy_remaining_cap_j = q_final**2 / (2 * capacitance_f)
# потери = I²·R_eff(x), где R_eff включает вихревые потери снаряда (× overlap) # потери = I²·R_eff(x), где R_eff включает вихревые потери снаряда (× overlap)
overlap_during = inductance_model.overlap_fraction(discharge_sol.y[2]) overlap_during = inductance_model.overlap_fraction(discharge_sol.y[2])
r_eff_during = r_total_ohm + r_eddy_coeff * overlap_during r_eff_during = r_total_ohm + r_eddy_coeff * overlap_during
energy_dissipated_j = float( resistive_diss = float(np.trapezoid(discharge_sol.y[1] ** 2 * r_eff_during, discharge_sol.t))
np.trapezoid(discharge_sol.y[1] ** 2 * r_eff_during, discharge_sol.t) # при обрыве на минимуме тока в катушке остаётся энергия магнитного поля —
) # её гасит обратный диод (freewheel), учитываем как потери, чтобы энергия
# сходилась. Считаем ТОЧНО (с насыщением), а не ½LI².
residual_inductor_j = float(inductance_model.magnetic_energy(x_final, i_final))
energy_dissipated_j = resistive_diss + residual_inductor_j
kinetic_before_j = 0.5 * mass_kg * v_fire_mps**2 kinetic_before_j = 0.5 * mass_kg * v_fire_mps**2
kinetic_after_j = 0.5 * mass_kg * v_final**2 kinetic_after_j = 0.5 * mass_kg * v_final**2

View File

@@ -4,7 +4,7 @@ import pytest
from gausse.components.database import ComponentDatabase from gausse.components.database import ComponentDatabase
from gausse.optim.search_space import ( from gausse.optim.search_space import (
TUBE_WALL_CLEARANCE_M, PROJECTILE_BORE_CLEARANCE_M,
SearchBounds, SearchBounds,
crossover, crossover,
decode, decode,
@@ -25,8 +25,15 @@ def test_sample_genome_respects_bounds():
genome = sample_genome(DB, BOUNDS, rng) genome = sample_genome(DB, BOUNDS, rng)
assert BOUNDS.min_stages <= len(genome.stages) <= BOUNDS.max_stages assert BOUNDS.min_stages <= len(genome.stages) <= BOUNDS.max_stages
assert len(genome.inter_stage_gaps_m) == len(genome.stages) - 1 assert len(genome.inter_stage_gaps_m) == len(genome.stages) - 1
assert BOUNDS.tube_od_m_min <= genome.tube_od_m <= BOUNDS.tube_od_m_max assert BOUNDS.tube_inner_d_m_min <= genome.tube_inner_d_m <= BOUNDS.tube_inner_d_m_max
assert genome.projectile.diameter_m <= genome.tube_od_m - TUBE_WALL_CLEARANCE_M assert genome.tube_od_m == pytest.approx(genome.tube_inner_d_m + 2 * genome.tube_wall_m)
# снаряд влезает в бор трубы
assert genome.projectile.diameter_m <= genome.tube_inner_d_m - PROJECTILE_BORE_CLEARANCE_M + 1e-9
# масса снаряда не больше 200 г
import math as _m
mat = DB.projectile_materials[genome.projectile.material_idx]
mass = _m.pi * (genome.projectile.diameter_m / 2) ** 2 * genome.projectile.length_m * mat.density_kg_m3
assert mass <= BOUNDS.projectile_mass_max_kg + 1e-6
for stage in genome.stages: for stage in genome.stages:
assert 0 <= stage.wire_idx < len(DB.wires) assert 0 <= stage.wire_idx < len(DB.wires)
assert BOUNDS.turns_per_layer_min <= stage.turns_per_layer <= BOUNDS.turns_per_layer_max assert BOUNDS.turns_per_layer_min <= stage.turns_per_layer <= BOUNDS.turns_per_layer_max
@@ -53,8 +60,8 @@ def test_mutate_keeps_genome_within_bounds():
genome = mutate(genome, DB, BOUNDS, rng, rate=0.5) genome = mutate(genome, DB, BOUNDS, rng, rate=0.5)
assert BOUNDS.min_stages <= len(genome.stages) <= BOUNDS.max_stages assert BOUNDS.min_stages <= len(genome.stages) <= BOUNDS.max_stages
assert len(genome.inter_stage_gaps_m) == len(genome.stages) - 1 assert len(genome.inter_stage_gaps_m) == len(genome.stages) - 1
assert BOUNDS.tube_od_m_min <= genome.tube_od_m <= BOUNDS.tube_od_m_max assert BOUNDS.tube_inner_d_m_min <= genome.tube_inner_d_m <= BOUNDS.tube_inner_d_m_max
assert genome.projectile.diameter_m <= genome.tube_od_m - TUBE_WALL_CLEARANCE_M + 1e-9 assert genome.projectile.diameter_m <= genome.tube_inner_d_m - PROJECTILE_BORE_CLEARANCE_M + 1e-9
# decode должен всегда успевать без исключений после repair # decode должен всегда успевать без исключений после repair
decode(genome, DB, BOUNDS) decode(genome, DB, BOUNDS)