Verification, Sensitivity, and Operating Limits of a GPU-Accelerated Planar WCSPH Model for Hydrodynamic Ram

1. Introduction

• HRAM loading matters for tank survival, but numerical studies often validate a single resolution and boundary treatment against one dataset.Earlier Lagrangian finite-element approaches also encountered severe mesh distortion and time-step restrictions under large fluid deformation.Source: §1 · p. 4
• Prior HRAM-SPH work rarely verifies boundary conditions against analytical references or reports resolution convergence, leaving agreement with experiments open to compensating errors.The literature also has not isolated the fluid–wall treatment’s contribution to far-field attenuation.Source: §1 · p. 7
• The paper extends an Adami ghost-particle boundary condition to a moving, decelerating projectile in a GPU-accelerated planar WCSPH solver.Its novelty is the moving-body HRAM setting and exhaustive characterization enabled by the 2D reduction, not GPU acceleration with an Adami-type boundary alone.Source: §1 · pp. 6, 8
• The study combines closed-form boundary checks, three-point resolution convergence, coefficient sensitivity sweeps, conservation diagnostics, GPU throughput scaling, and comparison with experimental and semi-analytical evidence.The planar particle count is roughly two orders of magnitude below an equivalent 3D discretization at the same spacing, making extensive characterization affordable.Source: §1 · pp. 7–8
• The planar reduction changes wave decay from the real 3D r^-1 pattern to r^-1/2, so raw pressure amplitudes overpredict experiment by roughly two orders of magnitude.The experimental comparison is therefore intended to assess timing, ordering, and sensitivity structure rather than peak-pressure agreement.Source: §1 · p. 9
• The prescribed projectile trajectory does not respond to fluid forces, producing fluid energy growth and making global energy conservation unavailable as a check.The paper also reports that the equation-of-state stiffness and cavitation cutoff prevent closed-form verification at the production operating point.Source: §1 · p. 9

2.1. SPH Fundamentals

• Fluid quantities are carried by Lagrangian particles and reconstructed as kernel-weighted sums over neighbours within the support radius.The solver uses the classical cubic B-spline kernel in its two-dimensional normalization, with smoothing length h = 1.3 Δx and support radius 2h = 2.6 Δx.Source: §2.1 · pp. 9–10

2.2. Governing Equations

• XSPH regularisation advects particles at a locally smoothed velocity to suppress tensile pairing and short-wavelength trajectory crossing.The blending coefficient is ε_XSPH = 0.02, with no appreciable bulk dissipation reported.Source: §2.2 · p. 10
• The symmetric specific-momentum balance uses an antisymmetric pairwise pressure term, conserving total linear momentum to machine precision.The Monaghan–Gingold dissipative kernel activates only when particle pairs approach, leaving expanding regions such as the projectile’s rarefaction wake dissipation-free.Source: §2.2 · pp. 10–11
• The model retains linear artificial viscosity with α = 0.03 and omits the quadratic term to avoid over-damping the finite-amplitude pressure pulse.The physical water-viscosity contribution is also retained, though it is negligible relative to artificial viscosity at the stated Reynolds numbers.Source: §2.2 · p. 11

2.3. Density Summation

• Direct density summation avoids continuity-equation drift in the projectile wake, while a [0.5, 1.5]ρ0 clamp limits pressure spikes without δ-SPH diffusion.The trade-off is a bounded inconsistency where bulk summed density meets Adami-assigned ghost density, identified as the likely source of non-monotone near-probe convergence.Source: §2.3 · pp. 11–12

2.4. Closure: Mie–Grüneisen Equation of State

• The Mie–Grüneisen closure uses a cavitation cutoff of −0.1 MPa and full liquid stiffness, which obstruct static-column and added-mass verification at the production operating point.The stiffness makes the static-column check numerically undiscriminating, while the cutoff clips the analytical solution’s suction side.Source: §2.4 · pp. 12–13

2.5. Time Integration

• At production resolution and 900 m/s, the CFL condition sets the leapfrog timestep to approximately 0.21 µs.The timestep uses the usual CFL condition.Source: §2.5 · p. 13

2.6.1. The Projectile: Adami Ghost Particles on a Moving Disk

• The moving disk uses 100 uniformly spaced ghost particles whose interpolated pressure, density, and velocity enforce no-penetration and pressure continuity without a contact-force penalty.Ghost density is recovered from pressure through the linear acoustic branch; total ghost line-mass is independent of particle count.Source: §2.6.1 · pp. 13–14

2.6.2. The Channel Walls: Two Treatments, Compared Directly

• Channel-wall ghost particles reuse the disk interpolation kernel with zero wall velocity and acceleration, while the hard free-slip baseline lets near-wall density and pressure collapse.The baseline particle density within approximately 2h of the wall falls to 0.5–0.7 ρ0.Source: §2.6.2 · p. 14

2.6.3. Velocity Ramp

• The projectile velocity ramps linearly from zero over 5 µs to avoid an artificial impulsive-entry shock.The ramp is short relative to the 50–90 µs acoustic transit time to either probe, and its influence on extracted pressure is confined to t < tramp.Source: §2.6.3 · p. 15

2.7. Projectile Dynamics

• The quadratic drag law uses the physical 3D projectile parameters, including mass 8.01 g, radius 6.25 mm, and Cd = 0.05.The mass and drag coefficient are not rescaled to a per-unit-depth 2D projectile.Source: §2.7 · p. 15
• The projectile trajectory is integrated as a 3D rigid-body ODE decoupled from the planar fluid, which receives no momentum feedback from the fluid.The matched trajectory supplies ghost-particle kinematics but does not independently predict projectile motion.Source: §2.7 · pp. 15–16
• The drag-law trajectory reproduces Varas et al. ’s velocity decay and penetration history at both impact speeds, with Figure 2 showing the 900 m/s comparison.The drag coefficient is calibrated against the experiment’s implied deceleration, so the trajectory is a matched boundary condition.Source: §2.7 · pp. 15–16, Figure 2 · p. 16

2.8. GPU Implementation

• Numba compiles the compute-heavy GPU steps, including neighbor search, density summation, ghost interpolation, momentum evaluation, drift, integration, and wall enforcement.The implementation uses a CUDA-executable path on an NVIDIA TITAN RTX and provides a multicore CPU fallback.Source: §2.8 · p. 16
• Figure 3 compares projectile velocity decay and position over time at 600 m/s.The caption identifies the two plotted comparisons but does not report their outcomes.Source: Figure 3 · p. 17
• The workflow separates verification against closed-form solutions from the later comparison with Varas et al., with Figure 4 depicting the solver pipeline and evidence tracks.The flow proceeds from governing-equation discretization through time integration to verification, convergence/sensitivity, and validation.Source: §2.8 · p. 17, Figure 4 · p. 18
• At each timestep, the CFL-limited step is recomputed from the maximum velocity across fluid and ghost particles before the leapfrog half-step is committed.Figure 4 makes this timestep check an explicit decision node.Source: §2.8 · p. 17

2.9. Verification Protocols

• A padded-domain disk test compares surface pressure and net force against a derived potential-flow solution for an impulsively started circular cylinder.The planar added-mass coefficient is the full displaced-fluid mass per unit depth, unlike the one-half factor for a 3D sphere.Source: §2.9 · p. 19
• A stationary-disk acoustic-pulse test compares reflected-to-incident amplitude with the normal-incidence steel-water reflection coefficient, approximately 0.92.The stated steel properties give an impedance ratio Zs/Zw of approximately 25; this test, like the added-mass case, isolates the boundary operator from channel geometry and drag.Source: §2.9 · p. 19
• At production settings, the −3.75 MPa added-mass suction peak is clipped by the −0.1 MPa cavitation cutoff, so comparison is dominated by the clamp rather than the boundary.The boundary operator can still be checked in reduced-c0 and reduced-U regimes where the obstructions lift.Source: §2.9 · p. 20

• Resolution convergence is assessed at Δx = 2.0, 1.5, and 1.0 mm with h = 1.3Δx and v0 = 900 m/s, using standard Richardson/GCI construction.The apparent order is estimated from the three resolutions.Source: §2.10 · p. 20
• Ghost-particle count and artificial viscosity are swept independently at both velocities while all other settings are held fixed.The tested values are Ng = 50, 100, 200 and α = 0.01, 0.03, 0.06, 0.10.Source: §2.10 · p. 20
• Total mass and mechanical energy are logged at every production-run output step, and GPU throughput is measured across approximately 10^4–10^6 particles.The throughput range is chosen to check whether wall-clock scaling is sensible across particle counts.Source: §2.10 · p. 20

• Normalized RMSE over a common window supplements peak pressure and arrival-time comparisons against Varas et al..The additional metric checks pressure-trace decay shape and cumulative impulse, which a peak and arrival time alone do not establish.Source: §2.11 · pp. 20–21
• Varas et al. pressure traces are digitized with WebPlotDigitizer v4.6, calibrated to printed axis ticks, and sampled at ≥40 points per trace.Three independent-operator redigitizations estimate uncertainty at approximately ±4% of local pressure and ±1 µs in time; timing uncertainty is carried as an error bar for near-probe timing.Source: §2.11 · p. 21

• The planar model represents a 12.5 mm steel sphere as a disk entering a water-filled channel 75 mm in half-width and 750 mm long at 900 or 600 m/s.At Δx = 1.5 mm, the production lattice contains approximately 50,000 particles before disk-interior particles are removed, roughly two orders of magnitude fewer than a full 3D lattice at the same spacing.Source: §3 · p. 21
• PTn and PTf are placed at x = 30 mm and 105 mm, respectively, matching the projected transducer locations reported by Varas et al..Both lie at the tank mid-height, y = 75 mm, clear of the near-wall kernel-truncation band.Source: §3 · p. 21
• The near probe samples the entry-shock near field, while the far probe samples the reflection-dominated far field.Together, the stations bracket the two loading regimes distinguished by the paper’s four-phase decomposition.Source: §3 · p. 22

4.1. Verification Status: What the Closed-Form Checks Can and Cannot Establish

• Two of three designed closed-form checks cannot be passed at the production operating point, so the boundary is characterized but not fully verified there.The hydrostatic signal is overwhelmed by density-error amplification under the stiff equation of state; the disk suction solution is clamped by cavitation, and the steel–water reflection coefficient was not measured. Reduced-closure checks assess the same boundary operator but leave its interaction with the production closure unresolved.Source: §4.1 · pp. 22–23

4.2. How Much Does Resolution Actually Matter?

• The near-probe peak is non-monotone across meshes, while the far-probe peak refines monotonically and supports only a coarse apparent order estimate.PTn changes 785.8 → 850.3 → 830.5 MPa, so Richardson/GCI is not defensible; the production value carries a practical ±32 MPa (∼4%) mesh uncertainty. PTf changes 658 → 714.9 → 737.5 MPa, with apparent order porder ≈2.7 and an extrapolated limit ≈752 MPa, about 2% above production. The bounded ghost–fluid density inconsistency likely explains the near-boundary behavior.Source: §4.2 · pp. 24–25, Figure 5 · p. 25, §4.6 · pp. 30–31

4.3. Does the Choice of Timestep Quietly Change the Answer?

• Changing CFL across a 4× range changes results by 3.5% or less, indicating the production timestep is not the limiting factor.The fixed-resolution production case was repeated at CCFL = 0.125, 0.25, and 0.50 to isolate temporal from spatial error.Source: §4.3 · p. 26, Table 2 · p. 26

4.4. Is the Particle Field Actually Smooth, or Just Reported That Way?

• Across the production run, the disorder index remains bounded and high-wavenumber pressure power stays a small fraction of the total.The diagnostics track wake-region nearest-neighbour disorder and pressure-field power above k > π/2h along the probe axis, providing a direct view of particle-field smoothness.Source: §4.4 · pp. 26–27, Figure 6 · p. 27

4.5. Sensitivity: Which Knobs Actually Matter

• Ghost-particle count has little effect once boundary sampling is adequate, whereas artificial viscosity affects PTn peak pressure more strongly at 600 m/s than at 900 m/s.Ghost mass normalization makes total ghost line-mass independent of Ng, and Ng = 100 adequately overlaps within kernel support. The retained α = 0.03 is the smallest swept value that suppresses spurious oscillations at both speeds; the 600 m/s overprediction relative to experimental attenuation moves in the expected direction with α, while raw planar amplitudes remain distinct from geometry-driven experimental peaks.Source: §4.5.1 · pp. 27–28, Figure 7 · p. 28, §4.5.2 · pp. 28–29

4.6. Numerical Provenance of the Implementation Pipeline

• The implementation reprojects ghost data before each spatial update, while a bounded ghost–fluid density inconsistency at the interface likely drives non-monotone near-probe convergence.The operator sequence is G, then Lh, then T∆t, with ghost data tied to current projectile kinematics. Because the reference sound speed dominates throughout and Ma never exceeds approximately 0.6, viscosity-driven numerical dissipation is proportionally larger at 600 m/s than at 900 m/s.Source: §4.6 · pp. 30–31

4.7. Conservation and Repeatability

• The total energy rises 3.18% by 150 µs, because the prescribed, one-way-coupled projectile boundary supplies work without feeling the fluid reaction.Energy is closed within 0.40% for the first ≈75 µs; over the full run projectile kinetic energy falls by 8.53 × 10^5 J/m while fluid energy rises by 1.14 × 10^6 J/m. The fluid gain is 1.33 times the projectile loss. Global energy conservation is therefore not an available check; a boundary-work audit is the appropriate diagnostic. The reported pressure histories lie within the early near-closed budget window, but this does not prove local pressure is free of coupling artefacts; the acceleration contribution is argued to be sub-percent during the shock phase. Rotating or offsetting the lattice changes PTn peak pressure by no more than 1.8%, while PTf is more sensitive.Source: §4.7 · pp. 31–34, Figure 9 · p. 34

4.8. Speed: Where Does the Time Actually Go

• Production-mesh throughput is near-linear in particle count, and the 150 µs event runs in about 14 s on one consumer-class GPU.The three meshes span N = 2.8 × 10^4 to 1.1 × 10^5, with times of 0.132, 0.230, and 0.576 min and near-constant per-particle cost of 4.6–5.1 × 10^-6 min/particle. The planar run is reported as two to three orders of magnitude faster than the equivalent 3D discretization, enabling the section’s parameter sweeps.Source: §4.8 · pp. 34–35, Figure 10 · p. 35

4.9. Now, Finally, the Comparison Against Experiment

• The near-probe wave arrives ∼8 µs ahead of the experimental trace at both velocities, while near/far and velocity orderings follow expected patterns; peak-pressure differences are not validation metrics.The planar model overpredicts pressure at both probes by two orders of magnitude, so the comparison characterizes timing and qualitative ordering rather than 3D amplitude agreement.Source: §4.9 · p. 37, Table 6 · p. 38

4.10. Where This Sits Relative to a Different Kind of Solver

• The planar model trades out-of-plane relief and direct 3D amplitude agreement for a rigid, closed-form-verifiable boundary and roughly two orders fewer particles.Planar wave amplitude decays as r^-1/2 rather than the 3D r^-1; the planar-to-spherical amplitude ratio is ≈2.2 at the near probe and ≈4.1 at the far probe, and the actual overestimate can be larger for the extended, sustained projectile source.Source: §4.10 · pp. 39–40

4.11. Positioning Against the Broader HRAM Numerical Literature

• The literature comparison frames the contribution as a combination of verification depth, convergence characterization, and speed, while acknowledging that deformable-wall methods capture wall deformation this rigid-wall solver omits.Closed-form boundary-condition verification was attempted but did not close, a caveat made explicit in the comparison.Source: §4.11 · p. 41

4.12. What This Solver Is, and Is Not Yet, Ready For

• The solver is suited to early-time (≲75 µs) relative screening studies, but not yet to absolute 3D peak-pressure prediction, late-time drag- and cavitation-phase impulse, or permanent wall-deformation and failure prediction.The stated limits follow from planar wave decay, one-way-coupling energy drift beyond ∼75 µs, and the rigid-wall idealization, respectively.Source: §4.12 · pp. 41–42

5. Conclusions

• The planar solver was examined against analytical solutions before experimental comparison, with joint convergence, coefficient-sensitivity, conservation, and lattice-orientation reporting enabled by its modest particle count.The experimental comparison is explicitly framed with the model’s dimensional simplification stated.Source: §5 · p. 42
• Hydrostatic and added-mass checks encounter limits shaped by equation-of-state stiffness, hydrostatic head, and the cavitation cutoff; impedance verification needs a second material.The verification programme characterizes rather than closes the boundary condition’s operating envelope.Source: §5 · p. 42
• Prescribed projectile motion without backreaction accompanies modest total energy growth, with the fluid absorbing more energy than the projectile relinquishes.The diagnostics identify global energy conservation as an inappropriate check for this one-way-coupled formulation and point to two-way coupling as an extension.Source: §5 · pp. 42–43
• Planar geometry, more than the boundary condition, governs near-probe peak agreement and wave-arrival timing; the model overpredicts far-probe PTf and raw amplitudes diverge from 3D experiments.The far-probe result does not establish the wall-ghost treatment as a recovery mechanism.Source: §5 · p. 43
• The full campaign, including reported sweeps, completes in a small fraction of the roughly eighteen minutes per equivalent 3D case, with planar particle counts about two orders of magnitude smaller.The throughput advantage is attributed to the smaller planar particle count.Source: §5 · p. 43
• The proposed roadmap prioritizes axisymmetric geometry and two-way projectile coupling; lower c0, a relaxed cavitation cutoff, and tests above 900 m/s would broaden the validated envelope.Quantitative reliability of full-field pressure, velocity, and density data depends on the axisymmetric and two-way coupling extensions.Source: §5 · pp. 43–44