Theory

PIBS Theoretical Model

This document derives the physics and numerics actually implemented in pibs/ballistics/, so you can read the code (or trust its output) with a clear picture of what it's computing. It follows the "Serebryakov system" of reduced-variable interior ballistics, as used in Soviet/PRC ballistics literature (see References).

Everything here is a direct read of base_gun.py, gun.py, recoilless.py and prop/prop.py — not a black-box description.

1. Scope: what "interior ballistics" means here

Interior ballistics covers the phase from primer ignition to shot exit: propellant burns, gas pressure builds, the shot accelerates down the bore. PIBS solves the forward problem (given a design, compute its pressure/velocity trace, §1-§10) and the inverse/constrained problem (given performance targets, solve for the grain web size and barrel length, §11) — both available in either GUI.

The model is 0-dimensional in the combustion chamber (space-mean pressure) with a 1-dimensional correction for the pressure gradient along the bore (§6). It is not a CFD simulation — there's no radial structure, no turbulence, no ignition transient. That's the right trade-off for design-level accuracy in seconds rather than hours.

2. Reduced (dimensionless) variables

Rather than integrate pressure, travel, velocity and time in SI units directly, the solver non-dimensionalizes everything against charge/gun-specific scales. This is what "reduced form" means in the README, and it's why gun.py variable names look like p_bar, l_bar, v_bar, t_bar — _bar denotes the dimensionless version.

Quantity	Scale	Reduced variable
Pressure	p_scale = f·Δ	p_bar = p / p_scale
Travel	l_0 = V_0 / S (chamber volume ÷ bore area — "equivalent chamber length")	l_bar = l / l_0
Velocity	v_j = √(2 f w / (θ φ m))	v_bar = v / v_j
Time	t_scale = l_0 / v_j	t_bar = t / t_scale
Burnt fraction	web-normalized	Z = e / e_1 ∈ [0, Z_b]

where:

• f — propellant force (specific energy of combustion, J/kg)
• Δ = w / V_0 — loading density (charge mass / chamber volume)
• w — charge mass, m — shot mass
• θ — reduced adiabatic index of propellant gas (θ = γ − 1)
• φ — the fictitious mass factor (§5): shot mass alone under-accounts for how much mass the pressure has to accelerate, because the propellant gas itself has inertia.
• e_1 — grain half-web (half-thickness of the burning layer); Z=1 is when the primary grain shape burns through, Z=Z_b ≥ 1 is complete combustion (see §3 for multi-perforated grains).

v_j is the theoretical maximum velocity the charge could impart to the shot in the idealized limit of complete, instantaneous burning with no losses — it's the natural velocity scale for the problem.

3. Propellant charge model

3.1 Form functions

As a grain burns, its burning surface area changes with the fraction of web consumed — this is the "form function". PIBS uses the standard cubic form function:

ψ(Z) = χ·Z·(1 + λ·Z + μ·Z²)          for 0 ≤ Z ≤ 1
ψ(Z) = χₛ·Z·(1 + λₛ·Z)               for 1 < Z ≤ Z_b   (secondary/sliver phase)
σ(Z) = dψ/dZ                          (relative burning surface)

ψ(Z) is the burnt mass fraction; σ(Z) its derivative, i.e. the burning-surface area relative to the initial grain surface. The coefficients (χ, λ, μ) — and (χₛ, λₛ, Z_b) for grains that keep burning past Z=1 — are not free parameters; they're derived from grain geometry:

• Simple geometries (SimpleGeometry in prop/prop.py): sphere, cylinder, tube, strip. Two shape ratios α = e₁/b, β = e₁/c (web vs. the grain's other two characteristic dimensions) fully determine (χ, λ, μ); these grains burn out at Z_b = 1 (no secondary phase — a sphere just gets smaller and vanishes).
• Multi-perforated geometries (MultPerfGeometry): 7/14/19-perf cylinders/rosettes/hexagons. These start out progressive (burning surface increases as the perforations enlarge) until the web between perforations is consumed (Z=1), after which the grain fragments into slivers that burn out degressively for 1 < Z ≤ Z_b. This progressive-then-degressive shape is why multi-perf propellant is standard in high-performance guns — it front-loads surface area growth to partially offset the falling pressure as chamber volume increases with shot travel.

A Propellant can mix two grain populations (main + auxiliary, with independently set geometry, size and web ratio) — used for igniter/ auxiliary charges in some of the bundled examples — plus an optional combustible-cartridge fraction (case material that itself burns and contributes gas).

3.2 Burn rate law

Grains regress at a rate governed by Vieli's/Muraour's law — linear in pressure raised to a (propellant-specific) exponent:

dZ/dt = u₁·p^n / e₁

u₁ (burn rate coefficient) and n (pressure exponent) come from the propellant composition (propellants.csv — 40 bundled compositions sourced from published references), separate from geometry (grain shape) — the two combine multiplicatively in the Propellant class (composition × geometry = grain).

4. Equation of state and space-mean pressure

Propellant gas is modeled with the Nobel–Abel (covolume) equation of state — ideal gas with a correction α for the finite volume the gas molecules themselves occupy at high density. Combining this with conservation of energy (chemical energy released by burning = kinetic energy of shot + gas + work done against bore friction and ambient counter-pressure) yields the space-mean pressure as an algebraic function of the current state, without needing to integrate a separate energy equation:

p_bar = (ψ − v_bar²) / (l_bar + l_ψ_bar)

l_ψ_bar = 1 − Δ·[(1−ψ)/ρ_p + α·ψ]

l_ψ_bar is the "unburnt-propellant and covolume" correction to the available expansion length — it accounts for the volume still occupied by solid, unburnt propellant ((1−ψ)/ρ_p) and the covolume of the gas already produced (α·ψ). This is Gun.f_p_bar() in the code, and it's evaluated fresh at every integration step — it is not itself an ODE state variable.

5. Equation of motion

With p_bar available algebraically, the reduced-variable ODE system (time domain) is:

dZ/dt_bar     = √(θ/2b) · p_bar^n              (burn progress)
dl_bar/dt_bar = v_bar                           (kinematics)
dv_bar/dt_bar = (θ/2)·(p_bar − p_ad_bar(v_bar)) (Newton's 2nd law, reduced)

b is a dimensionless grouping of the burn-rate/geometry/energy terms:

b = S²e₁² / (f·φ·w·m·u₁²) · (f·Δ)^(2(1−n))

(S = bore area). p_ad_bar(v_bar) is the ambient aerodynamic counter-pressure (§7).

gun.py actually solves this system in three different independent variables depending on which is most numerically convenient at each stage: ode_t (time), ode_l (travel — needed once v_bar is no longer small, since 1/v_bar blows up at t=0), and ode_z (burnt fraction — used to robustly locate burnout/exit even when dZ/dt→0 near full burn). This is why the source has three near-identical ode_* methods; they're the same physics, re-parameterized for numerical robustness at different phases of the shot.

The φ (fictitious mass) factor

Newton's law needs some mass to divide by, but only the shot's mass is truly rigid — the propellant gas behind it also accelerates and carries kinetic energy. Rather than track the gas velocity field explicitly, PIBS lumps this into a correction factor:

φ = φ₁ + λ₂·cc·(w/m)

φ₁ = 1/(1 − drag_coefficient) is the separate correction for bore friction/engraving resistance (the "Resistance %" input). λ₂ and the chambrage correction cc come from the pressure-gradient solution, §6 — this is where the charge-to-shot mass ratio and gun geometry feed back into the acceleration itself, not just the pressure distribution.

6. Pressure gradient along the bore (the Lagrange problem)

Space-mean pressure (§4) is a single number, but at any instant the local pressure is highest at the breech (accelerating the gas column behind it) and lowest at the shot base — this is the classical "Lagrange problem" of interior ballistics: gas of nonzero mass filling a tube, accelerated by pressure with a fixed shot at one end.

PIBS solves this via two dimensionless distribution coefficients λ₁, λ₂, computed once per gun from one of three closed-form/ semi-closed-form models (the "Gradient" setting, Solutions in ballistics/__init__.py):

• SOL_LAGRANGE: the classical simplifying assumption — gas density uniform along the tube at every instant. Gives the closed form λ₁ = 1/2, λ₂ = 1/3 directly, no iteration needed.

• SOL_PIDDUCK: Pidduck's self-similar solution for a gas with adiabatic index θ+1, solved via the integral equation

  ∫₀¹ (1 − Ω·ξ²)^[1/(k−1)] dξ  =  [w/(2φ₁m)]·[(k−1)/k]·(1−Ω)^[k/(k−1)] / Ω
  

  for Ω (root-found with the Dekker method) — w/(φ₁m) is the charge mass relative to the friction-corrected ("fictitious") shot mass, not the raw shot mass. λ₁, λ₂ follow from Ω by further quadrature. This is the more physically complete model and is the default.

• SOL_MAMONTOV: the k→1 (isothermal-gas) limit of the same Pidduck family, evaluated as a separate closed form to avoid the numerical 0/0 at k=1.

(pidduck() in gun.py implements the shared machinery for the latter two; see its docstring for the exact source citation.)

Given λ₁, λ₂ and a chambrage correction cc (which accounts for the chamber being wider than the bore, χ_k = chamber area / bore area), the breech and shot-base pressures are recovered from the space-mean pressure:

P_s / P  =  1 / (1 + λ₂'·w/(φ₁·m))
P_b / P  =  (φ₁·m + λ₂'·w) / (φ₁·m + λ₁'·w)

with λ₁', λ₂' further adjusted for how far the shot has travelled (the chamber's relative contribution to the total gas column shrinks as travel increases). Pressure at an arbitrary axial probe point between breech and shot is then interpolated assuming a quadratic (Lagrangian) profile in the gas velocity field — this is what powers the "Barrel Trace" plot in the desktop app and the tube-strength calculation.

This is why "Peak Avg Pressure", "Peak Breech Pressure" and "Peak Shot Pressure" are three different numbers occurring at three different times — each is found by an independent Golden Section Search over the relevant pressure trace (§8), not read off the same peak.

7. Ambient/aerodynamic drag

Once the propellant is exhausted (or even before, for high-velocity shots), the atmosphere in front of the shot pushes back. This is modeled as a counter-pressure term derived from the shot's velocity relative to the local speed of sound, using the ambient density, pressure and adiabatic index inputs (func_p_ad_bar() in base_gun.py) — a standard supersonic-projectile drag correction. It can be switched off (setting ambient density to zero), in which case this term vanishes and the gun is solved as if firing in vacuum ahead of the shot.

8. Characteristic points and numerical methods

The solver doesn't just uniformly sample the trajectory — it explicitly locates named events, each via the appropriate numerical technique:

Event	Meaning	Method
SHOT_START	ignition, Z=Z₀ (see below)	initial condition
FRACTURE	Z=1, primary grain shape consumed (multi-perf only)	direct integration to Z=1
BURNOUT	Z=Z_b, propellant fully consumed	direct integration to Z=Z_b
SHOT_EXIT	l = l_g (barrel length reached)	direct integration to l_bar=l_g_bar
PEAK_AVG_P / PEAK_BREECH_P / PEAK_SHOT_P	maxima of the respective pressure trace vs. time	Golden Section Search (gss()), independently per pressure definition

Z₀ (the burnt fraction at ignition) is itself root-found: given a user-specified start pressure (the pressure needed to shear the shot's rotating band / begin motion), PIBS solves ψ(Z₀) = ψ₀ where ψ₀ comes from inverting the Nobel-Abel EOS at that pressure and zero travel — via the Dekker method (dekker(), a bisection/secant/ inverse-quadratic hybrid similar to Brent's method).

The ODE system itself is integrated with an adaptive embedded Runge–Kutta–Fehlberg scheme (rkf()), which self-adjusts step size to a user-specified relative tolerance (-log10(ε) in the UI — e.g. a value of 4 means tolerance 1e-4) and supports early-abort conditions (used to detect pressure exceeding the 600 MPa validity ceiling of the Nobel-Abel EOS, or a squib condition where the shot stalls in the bore).

9. Recoilless guns

A recoilless gun vents propellant gas backward through a nozzle, countering the breech reaction so the weapon needs no recoil mechanism. Recoilless (in recoilless.py) extends the conventional model with two additional coupled state variables tracked through every integration step:

• η (eta) — cumulative mass fraction of propellant gas that has flowed out through the nozzle (as opposed to still filling the chamber/bore behind the shot)
• τ (tau) — a reduced temperature/energy ratio, since gas leaving through the nozzle carries away energy, which the simple closed-form pressure equation (§4) no longer captures on its own

The space-mean pressure becomes p_bar = τ/(l_bar+l_ψ_bar)·(ψ−η) (compare to §4 — note the (ψ−η) in place of ψ: only the gas mass retained in the gun contributes to pressure) and dη, dτ are integrated alongside Z, l_bar, v_bar.

The recoilless condition — zero net reaction force — is enforced by sizing the throat area, not by an extra force-balance ODE. At construction time, Recoilless.__init__ computes the required dimensionless throat area s_j_bar from the nozzle expansion ratio (a_bar) and nozzle efficiency (χ₀, an empirical loss factor for non-ideal nozzle flow) via the standard isentropic converging- diverging nozzle thrust-coefficient relation, then checks it physically fits within the breech face (s_j_bar ≤ chambrage ratio) — raising a clear error if it doesn't (an infeasible design, not a numerical failure).

Because a meaningful fraction of the propellant's chemical energy is deliberately vented rather than used to accelerate the shot, recoilless guns have markedly lower thermal/ballistic efficiency than closed-breech guns for the same propellant — confirmed by the bundled examples: the 105mm M40A1 recoilless design solves to ~7% thermal efficiency vs. ~33% for the 76mm ZiS-3 conventional gun (see the User Guide for the full comparison). This is the expected physical trade-off for eliminating recoil, not a modeling artifact.

10. Efficiencies

Three efficiency figures are derived from the solved trace (GenericResult.get_eff()):

thermal efficiency    te = (v_muzzle / v_j)²
ballistic efficiency  be = te / φ
piezometric efficiency pe = ½·φ·m·v_muzzle² / (P_max·S·l_g)

te measures how much of the theoretical maximum velocity (§2) was actually realized. be further corrects for the fictitious-mass factor φ, isolating losses attributable to charge-to-shot mass ratio and pressure gradient rather than incomplete combustion. pe measures how efficiently peak pressure was "used" — a gun whose pressure trace is flat and high for the whole travel (ideal) has pe→1; one with a sharp early spike that decays quickly wastes barrel length at sub-peak pressure and has low pe.

11. Constrained design, optimization, and structural sizing

Everything above solves the forward problem: given a design, compute its performance. The pieces below build on that same core solver to solve the inverse problem (given performance targets, find a design) and to size the resulting barrel. All four are available in both the desktop app and the web GUI.

11.1 Constrained design (the inverse problem)

ConstrainedGun/ConstrainedRecoilless (constrained_gun.py, constrained_recoilless.py) solve for grain web size and barrel length given a target velocity v_d and pressure p_d, via nested root-finding:

1. For a chosen charge-to-shot mass ratio w/m and chamber loading density Δ/ρ_p ("load fraction"), the charge mass and chamber volume — and hence l_0, Δ — are fixed.
2. The web (via half-web e_1) is root-found with the Dekker method so that the chosen pressure constraint point (space-mean, breech, or shot-base — the "Pressure Constraint Point") exactly peaks at the design pressure p_d: the search probes outward from a small starting web, doubling/halving until the sign of p_peak(e_1) − p_d brackets a root, then converges on it.
3. With e_1 fixed, the gun is integrated forward from that peak point until velocity reaches v_d, giving the required barrel length l_g.
4. Because the chambrage correction cc (§6) depends on l_g, and the pressure-gradient coefficients used in step 2 depend on cc, steps 2-3 repeat with the newly-converged cc until l_g itself stops changing (bounded by max_iteration).

An error here ("Solution requires excessive tube length...", "Design velocity exceeded before peak pressure point...") means the chosen charge ratio and loading density can't physically reach the target within the barrel-length search ceiling — not a numerical failure.

11.2 Optimization: minimum bore volume / minimum barrel length

Constrained.find_min_v() treats load fraction as a free variable instead of a fixed input: it searches it with Golden Section Search between the minimum load fraction (below which the design pressure is unreachable even in the closed-bomb limit — psi_0 would exceed 1, §8) and the maximum (found by maximum_load_fraction(), which probes upward until the constrained solve in §11.1 stops converging), minimizing either total chamber+barrel volume (MIN_BARR_VOLUME) or barrel length alone (MIN_PROJ_TRAVEL).

11.3 Tube strength & autofrettage

Given the solved pressure trace, Gun.structure() sizes barrel wall thickness at each axial probe point using thick-walled cylinder theory under the Tresca (maximum shear stress) yield criterion, against the highest pressure ever reached at that axial position across the entire trace (not just when the shot happens to pass it — the pressure envelope over time, via to_px_u() from §6), scaled by a safety factor.

Monoblock (single-piece barrel) — for inner:outer radius ratio k:

k = (1 − 2p/σ_y)^(−1/2)

Thickening the wall (raising k) lets it contain a higher p without the bore exceeding yield strength σ_y — but this construction is mathematically impossible above p = σ_y/2, regardless of wall thickness (the bore itself yields under a purely elastic stress distribution before the wall can help).

Autofrettage — the barrel is first pressurized beyond yield so it retains a residual compressive pre-stress at the bore once depressurized; the working pressure must first overcome that pre-stress before the bore sees any net tension. The optimal autofrettage radius ratio:

m = exp(p / σ_y),   k = m  (full autofrettage)

has no p = σ_y/2 ceiling — which is why high-pressure guns are routinely autofrettaged rather than built monoblock. Barrel volume (and hence estimated tube mass) is integrated from these radius ratios along the bore.

11.4 Guidance diagrams

A guidance diagram maps every design — across a grid of charge-to-shot mass ratios and chamber loading densities — that satisfies a given velocity/pressure target, showing how barrel length and chamber volume trade off across that space. It's §11.1's constrained solve repeated at every grid point: for each charge ratio, the maximum feasible load fraction is found (maximum_load_fraction()), then every load fraction from the minimum up to that maximum (in steps) is solved individually, keeping only the combinations that converge to a physically valid design.

The desktop app parallelizes this sweep across CPU cores (multiprocessing.Pool, since a single diagram can mean solving hundreds of individual designs). The web GUI instead runs it sequentially (guide_graph(sequential=True, ...)) — spawning worker processes from inside a Streamlit rerun is unsafe on Windows' spawn start method, since each worker would try to re-import and re-run the Streamlit script itself — so keep the sweep grid modest there.

References

As cited in the source (docstrings, pidduck() in gun.py, and the bundled example descriptions):

• 金志明 (2014). 《枪炮内弹道学》(Interior Ballistics of Guns).
• 鲍廷钰，邱文坚 (1995). 《内弹道学》(Interior Ballistics).
• 王连荣，张佩勤 (1987). 《火炮内弹道计算手册》(Gun Interior Ballistics Calculation Handbook), National Defense Industry Press.
• 兵器工业部第二管理局 (1982). 《国产火炮手册》(Domestic Gun Handbook).
• Multiple Russian-language ordnance service manuals (cited per-example in pibs/examples/*.json).