In order to determine how the bullet’s kinetic energy changes from the data presented above, we use models describing the velocity decay in Newtonian, viscoelastic and shear thickening fluids.

Inertial model

For a Newtonian fluid, the force F D counteracting on the object depends on the Reynolds number, which is the ratio between the inertial and viscous forces of the system:

$$Re=\frac{{\rm{inertial}}\,{\rm{forces}}}{{\rm{viscous}}\,{\rm{forces}}}=\frac{\rho {r}_{b}v}{\eta },$$ (1)

where ρ is the density of the fluid, r b the radius of the object (and hence the characteristic length scale of our system), η the viscosity and v the impact velocity of the object. For small Re (\(Re\ll 1\)), the viscous forces dominate. In this regime, F D is given by the Stokes drag equation for a spherical object6:

$${F}_{D}=6\pi \eta {r}_{b}\frac{dz(t)}{dt},$$ (2)

where z(t) is the position of the object at time t (with z(0) = 0 describes the moment of impact). For large Re (\(Re\gg 1\)), the inertial forces dominate and F D can be approximated with the drag equation:

$${F}_{D}=-\,{C}_{d}\pi \rho {r}_{b}^{2}\,{(\frac{dz(t)}{dt})}^{2},$$ (3)

with C d defined as the hydrodynamic drag coefficient of the object. Using r b = 3 mm and v = 100 m/s, we obtain a Reynolds number during impact on the order of 105 for water, which is indeed in the large-Re regime. To determine the Re regime of cornstarch and PVA the complex viscosity η* is used20:

$${\eta }^{\ast }=\frac{1}{\omega }\sqrt{{(G^{\prime} )}^{2}+{(G^{\prime\prime} )}^{2}}$$ (4)

Where ω is the angular frequency used in the rheometer measurements. Assuming that the object impact generates a high strain on the liquid during impact, the Reynolds numbers of cornstarch and PVA can be calculated. Using the storage and loss modulus (Fig. 1) for strains larger than 500%, the Reynolds number for both cornstarch and PVA are on the order of 10. Although the Reynolds number of both fluids are lower compared to water, the inertial forces still dominate the viscous forces during impact. Therefore, the drag should be inertial and not viscous (Eq. (3)). Using Newton’s second law, one can then obtain a differential equation for the velocity of the object, using \(v\equiv \frac{dz(t)}{dt}\):

$${m}_{b}\frac{dv}{dt}=-\,{C}_{d}\pi \rho {r}_{b}^{2}{v}^{2},$$ (5)

with m b the mass of the object. The solution for the impact velocity follows as:

$$v=\frac{{v}_{0}}{1+t/{\tau }_{in}},$$ (6)

where v 0 is the impact velocity and τ in an inertial velocity decay time. The theoretical value of the inertial velocity decay time is defined as \({\tau }_{in}=\frac{{m}_{b}}{{C}_{d}\pi {r}_{b}^{2}{v}_{0}\rho }\). This model will be referred to as the inertial model.

Viscoelastic model

To model a viscoelastic fluid, it is important to note that the fluid can show both viscous and/or elastic behaviour depending on the timescale of the impact20. One of the simplest and best-known models that describes viscoelasticity is the Maxwell model28,29, in which the viscous and elastic responses of a viscoelastic fluid are approximated with a damper (dashpot) and spring connected in series. The dashpot, with viscosity η, describes the viscous behaviour while the spring, with elasticity E, describes the elastic behaviour of the fluid. In this model, it is assumed that the dashpot and spring behave according to Newton’s and Hooke’s law, respectively, meaning that both η and E are independent of the magnitude of the introduced stress. A further assumption is that both the elastic and viscous responses are uncorrelated. Using the Maxwell model, the material response on the impact can be described with the differential equation of a damped harmonic oscillator:

$$\frac{{d}^{2}z(t)}{d{t}^{2}}+2{\omega }_{0}\zeta \frac{dz(t)}{dt}+{\omega }_{0}^{2}z(t)=\mathrm{0,}$$ (7)

Where \({\omega }_{0}\equiv \sqrt{\frac{E}{{m}_{b}}}\) is defined as the undamped angular frequency of the oscillator and \(\zeta \equiv \eta \mathrm{/2}\sqrt{{m}_{b}E}\) the damping ratio. Depending on the elasticity and viscosity of the Maxwell model, the response to an object impact can either be underdamped (ζ < 1) or overdamped (ζ > 1). If the system is underdamped, the bullet would oscillate within the fluid with a decrease in amplitude over time. However, as the bullet is completely stopped in both cornstarch and PVA (Fig. 3b) without oscillating, we assume that the system is overdamped and the viscous dissipation of the dashpot dominates. With this assumption, a solution of Eq. 7 can be found where the velocity decays exponentially with time:

$$v(t)={v}_{0}{{\rm{e}}}^{-\gamma t}$$ (8)

Where \(\gamma \equiv \frac{\eta }{2{m}_{b}}\) and v 0 is (again) the impact velocity of the object.

Goff et al. studied object impact on foams30 by using an adaptation of the Kelvin-Voigt model (a model where the spring and dashpot are connected in parallel instead of in series)20 where the plasticity of the foam was taken into account31. In contrast with Goff et al., plasticity does not play a role in our experiments as the fluids are not irreversibly deformed during impact. However, whether the viscous and elastic responses in the model are correlated or not might have an influence on the velocity decay in viscoelastic fluids at low speeds. For example, Akers et al. showed that in low velocity impacts (~2.43 m/s) on viscoelastic micellar fluids, the object oscillates. This can be described as a damped oscillation, but how the oscillations change over time is dependent on whether the elastic and viscous responses are correlated (Kelvin-Voigt model) or not (Maxwell model). However, as both cases predict an exponential decay for the velocity at high impact velocities, this is beyond the scope of our paper.

Added Mass model

Finally, we used the added mass model given by Waitukaitis and Jaeger13 to describe the shear thickening response of a suspension on an impacting object. In this model, the response on the external stress is approximated as an inelastic collision in which, as mentioned in the introduction, the objects pushes the suspended particles together to form a solidified plug. This plug pulls the surrounding liquid downward like a ‘snowplough’13, adding mass over time that the object needs to displace. Using Newton’s first law, one can obtain the following second-order differential equation:

$$({m}_{b}+{m}_{a}(t))\,\frac{{d}^{2}z(t)}{d{t}^{2}}=-\,(\frac{d{m}_{a}(t)}{dt}\frac{dz(t)}{dt})-{F}_{ext},$$ (9)

where m a (t) equals the increasing mass over time and F ext an external force which, in our case, is the gravitational force. Then, following the same logic as Waitukaitis et al., the added mass over time can be approximated with a cone-like region around the impacting object, increasing with time:

$${m}_{a}(t)=\frac{1}{3}\pi \rho C{({r}_{b}+kz(t))}^{2}kz(t),$$ (10)

with r b the radius of the impacting object with a circular surface area. C and k are defined as the “added-mass coefficient” and “front coefficient”, respectively (for details, see the paper of Waitukaitis and Jaeger13). The former is introduced to compensate for the fact that the impacted medium is a liquid instead of a solid, and the latter as the distance between the solidification front of the added mass and the object’s position13.

It is important to note that Waitukaitis and Jaeger only tested this model for a limited range of low impact velocities (0.2 < v < 2 m/s), where the object only deforms the surface rather than penetrating the material13. At high impact velocities, the object penetrates the liquid and creates an air cavity (Fig. 2b). The added mass model does not take the detailed form of the interface into account, and quantitative discrepancies may therefore occur.