2.1 Flow

Slender impactors have a cylindrical shaft (shank) capped by a curved nose. The flow of the material being penetrated by such impactors, relative to the impactor, can be approximated as a superposition of point sources upon a uniform background flow. In the simplest case (Tate Reference Tate1978; Yarin et al. Reference Yarin, Rubin and Roisman1995), for an impactor penetrating at speed $u(t)$ at time  $t$ , the relative velocity field of the impacted material may be obtained from the stream function

(2.4) $$\begin{eqnarray}\unicode[STIX]{x1D6F9}(r,z,t)=(R^{2}\cos \unicode[STIX]{x1D719}-2r^{2})u(t)/4,\end{eqnarray}$$

where $r$ and $z$ are as in figure 2(a), $\unicode[STIX]{x1D719}=\arctan (r/z)$ and $R$ is the radius of the impactor’s shank. The stream function $\unicode[STIX]{x1D6F9}$ represents the superposition of the flow field due to a growing spherical cavity of volumetric strength $\unicode[STIX]{x03C0}u(t)R^{2}$ per unit time and a background flow of rate $u(t)$ . Figure 2(b) shows the streamlines associated with (2.4) in the frame of the impactor. Note that, in a fixed frame, it is the impactor that penetrates at $u(t)$ which, in turn, is controlled by (2.1).

From (2.4), with $\tilde{r}=\sqrt{r^{2}+z^{2}}$ , we obtain velocity and strain-rate fields in the impacted material relative to the impactor, respectively,

(2.5) $$\begin{eqnarray}\boldsymbol{u}(r,z,t)=u(t)(\sin \unicode[STIX]{x1D719}\,\hat{\boldsymbol{e}}_{r}+\cos \unicode[STIX]{x1D719}\,\hat{\boldsymbol{k}})\frac{R^{2}}{4\tilde{r}^{2}}+u(t)\hat{\boldsymbol{k}},\end{eqnarray}$$

and

(2.6) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D63F}(r,z,t) & = & \displaystyle \frac{u(t)R^{2}}{4\tilde{r}^{3}}(\cos ^{2}\unicode[STIX]{x1D719}-2\sin ^{2}\unicode[STIX]{x1D719})\hat{\boldsymbol{e}}_{r}\otimes \hat{\boldsymbol{e}}_{r}+\frac{u(t)R^{2}}{4\tilde{r}^{3}}\hat{\boldsymbol{e}}_{\unicode[STIX]{x1D703}}\otimes \hat{\boldsymbol{e}}_{\unicode[STIX]{x1D703}}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{u(t)R^{2}}{4\tilde{r}^{3}}(\sin ^{2}\unicode[STIX]{x1D719}-2\cos ^{2}\unicode[STIX]{x1D719})\hat{\boldsymbol{k}}\otimes \hat{\boldsymbol{k}}-\frac{3u(t)R^{2}}{4\tilde{r}^{3}}\sin \unicode[STIX]{x1D719}\cos \unicode[STIX]{x1D719}(\hat{\boldsymbol{e}}_{r}\otimes \hat{\boldsymbol{k}}+\hat{\boldsymbol{k}}\otimes \hat{\boldsymbol{e}}_{r}),\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where $\otimes$ denotes the outer product. The relative acceleration field $\boldsymbol{a}(r,z,t)=a_{r}\hat{\boldsymbol{e}}_{r}+a_{\unicode[STIX]{x1D703}}\hat{\boldsymbol{e}}_{\unicode[STIX]{x1D703}}$ is obtained by computing $\unicode[STIX]{x2202}\boldsymbol{u}/\unicode[STIX]{x2202}t+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}$ . In particular,

(2.7) $$\begin{eqnarray}a_{r}(r,z,t)=\frac{R^{2}r}{4\tilde{r}^{3}}\left\{\dot{u}(t)-\frac{3z}{\tilde{r}^{2}}u(t)^{2}-\frac{R^{2}}{2\tilde{r}^{3}}u(t)^{2}\right\}.\end{eqnarray}$$

We note that (i) strain rates change with the impactor’s speed $u(t)$ , (ii) far from the impactor the flow field becomes uniform, i.e. $\boldsymbol{u}\rightarrow u(t)\hat{\boldsymbol{k}}$ when $\tilde{r}\gg R$ , and (iii) the strain-rate tensor’s magnitude $|\unicode[STIX]{x1D63F}|=\sqrt{\unicode[STIX]{x1D60B}_{ij}\unicode[STIX]{x1D60B}_{ij}}=\sqrt{6}uR^{2}/4\tilde{r}^{3}$ so that, at any $z$ , $|\unicode[STIX]{x1D63F}|\sim r^{-3}$ as $r\rightarrow \infty$ .

The zero-velocity streamline, which emanates from the stagnation point $z=z_{s}=-R/2$ , is an ovoid of Rankine:

(2.8) $$\begin{eqnarray}r=F(z)=\left\{\left(R^{2}-z^{2}+z\sqrt{z^{2}+2R^{2}}\right)/2\right\}^{1/2};\end{eqnarray}$$

this is shown in figure 2(b). We observe that the zero-velocity streamline imitates the shape of an impactor with a near-cylindrical shank and a curved nose very well. In fact, if the impactor’s shape was given by (2.8), then the flow would satisfy our expectation that the material velocity normal to the impactor’s surface vanishes. Of course, general torpedo-shaped impactors (for example, figure 1) will not strictly satisfy (2.8). Subsequently, we will ignore these differences in shape, and take the impactor’s shape to be the profile generated by (2.8). This assumes that modelling the impactor’s exact profile may not greatly affect its dynamics. This simplifying assumption is also a practical one, considering that the complex rheology of geomaterials is itself not well understood and mathematical descriptions incorporate many idealizations. Finally, we note that, when crucial, it is possible to systematically improve our approximation of the impactor’s shape by introducing point vortex sources (Tate Reference Tate1986a ,Reference Tate b ).

The kinematic description (2.5) works well given the impactor’s speed and slenderness, and because material flow around the impactor is severely confined. However, (2.5) will be unable to capture flow detachment or changes in the flow introduced by fluid friction. Nevertheless, we surmise that the kinematic description (2.5) will yield useful answers, at least at leading order, for any rheology, provided the impactor is slender and penetrates at a fast rate; for example, this held true for projectile impact into metals (Yarin et al. Reference Yarin, Rubin and Roisman1995).

2.2 Rheology

The next step is to select appropriate constitutive models for geomaterials such as soil and clay. Soils and clays may both be modelled as pressure-dependent Bingham fluids (Oldroyd Reference Oldroyd1947; Prager Reference Prager1961). Such fluids have a yield criterion that separates solid- and fluid-like behaviour. Yielding in soils depends on the local pressure, but this is not so in clays, wherein yielding is regulated by a maximum shear strength  $K$ . In either case, the rheology postyield may be represented as

(2.9) $$\begin{eqnarray}\unicode[STIX]{x1D748}=-p\mathbf{1}+\sqrt{2}K\unicode[STIX]{x1D63F}/|\unicode[STIX]{x1D63F}|,\end{eqnarray}$$

where $\unicode[STIX]{x1D748}$ is the postyield stress tensor, $p=-(\text{tr}\,\unicode[STIX]{x1D748})/3$ is the pressure and $\mathbf{1}$ is the identity tensor. The constitutive law above is rate-independent, and should be compared with the usual Amontons–Coulomb friction law, with $K$ playing the role of a ‘friction coefficient’.

The shear strength in water-saturated clays increases with vertical depth. For simplicity we set

(2.10) $$\begin{eqnarray}K=R_{f}k(h-z),\end{eqnarray}$$

where $R_{f}$ accounts for the rate dependence of clay (Casagrande & Wilson Reference Casagrande and Wilson1951) and the shear strength gradient $k$ models the variation with depth of the maximum shear stress. For example, in seabed clay $k\approx 1.5~\text{kPa}~\text{m}^{-1}$ (Freeman, Murray & Schüttenhelm Reference Freeman, Murray and Schüttenhelm1988). Following Mitchell & Soga (Reference Mitchell and Soga2005) we take

(2.11) $$\begin{eqnarray}R_{f}=1+\unicode[STIX]{x1D706}\max \{\ln |\dot{\unicode[STIX]{x1D6FE}}/\dot{\unicode[STIX]{x1D6FE}}_{ref}|,0\},\end{eqnarray}$$

where the rate parameter $\unicode[STIX]{x1D706}$ is a constant that estimates the increase in the clay’s strength per log cycle, and $\dot{\unicode[STIX]{x1D6FE}}$ and $\dot{\unicode[STIX]{x1D6FE}}_{ref}$ are magnitudes of the actual and reference strain rates, respectively. We will estimate $\dot{\unicode[STIX]{x1D6FE}}$ as $u/(2R)$ , where $u$ is the impactor’s penetration speed. Laboratory tests (O’Loughlin et al. Reference O’Loughlin, Richardson, Randolph and Gaudin2013) on clay estimate $\unicode[STIX]{x1D706}$ to be in the range 0.2–1 for $\dot{\unicode[STIX]{x1D6FE}}_{ref}=0.17$ , but suggest that in the field $\unicode[STIX]{x1D706}$ may lie within 0.15–0.2. Henceforth, we will set $\dot{\unicode[STIX]{x1D6FE}}_{ref}=0.17$ in all our computations. Finally, the ‘max’ in (2.11) is introduced to ensure that rate effects are included only if the strain rate $\dot{\unicode[STIX]{x1D6FE}}$ is greater than the reference strain rate $\dot{\unicode[STIX]{x1D6FE}}_{ref}$ .

2.3 Water-saturated clays

The computation of the impactor’s dynamics will now be demonstrated in the specific case of the penetration into seabeds of deep-penetrating anchors of the kind shown in figure 1. For this, we will employ the rheology of seabed clay as modelled through (2.9)–(2.11). Analysis of impacts into soils would follow the same procedure, but with a different rheology, and this will be taken up in the future.

The stress tensor during seabed penetration is found by combining (2.6) and (2.9):

(2.12) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D748} & = & \displaystyle -p\mathbf{1}+\frac{K}{\sqrt{3}}\{(\cos ^{2}\unicode[STIX]{x1D719}-2\sin ^{2}\unicode[STIX]{x1D719})\hat{\boldsymbol{e}}_{r}\otimes \hat{\boldsymbol{e}}_{r}+\boldsymbol{e}_{\unicode[STIX]{x1D703}}\otimes \hat{\boldsymbol{e}}_{\unicode[STIX]{x1D703}}\nonumber\\ \displaystyle & & \displaystyle +\,(\sin ^{2}\unicode[STIX]{x1D719}-2\cos ^{2}\unicode[STIX]{x1D719})\hat{\boldsymbol{k}}\otimes \hat{\boldsymbol{k}}-3\sin \unicode[STIX]{x1D719}\cos \unicode[STIX]{x1D719}(\hat{\boldsymbol{e}}_{r}\otimes \hat{\boldsymbol{k}}+\hat{\boldsymbol{k}}\otimes \hat{\boldsymbol{e}}_{r})\!\},\end{eqnarray}$$

where $\hat{\boldsymbol{e}}_{r}$ and $\hat{\boldsymbol{k}}$ are shown in figure 2(a) and $\hat{\boldsymbol{e}}_{\unicode[STIX]{x1D703}}=\hat{\boldsymbol{k}}\times \hat{\boldsymbol{e}}_{r}$ . The normal traction $\unicode[STIX]{x1D70E}_{nn}$ to be utilized in (2.3) is given by $\hat{\boldsymbol{n}}\boldsymbol{\cdot }\unicode[STIX]{x1D748}\boldsymbol{\cdot }\hat{\boldsymbol{n}}$ , where $\hat{\boldsymbol{n}}$ is the normal to the anchor’s surface defined by (2.8); see figure 2(a). However, the tangential traction $\unicode[STIX]{x1D70F}$ in (2.3) is not estimated well by $\hat{\boldsymbol{n}}\boldsymbol{\cdot }\unicode[STIX]{x1D748}\boldsymbol{\cdot }\hat{\boldsymbol{t}}$ , where $\hat{\boldsymbol{t}}$ is the unit tangent in the $\hat{\boldsymbol{e}}_{r}-\hat{\boldsymbol{k}}$ plane, as the velocity field (2.5) corresponds to an inviscid flow. Hence, we postulate that $\unicode[STIX]{x1D70F}$ is proportional to the maximum shear strength  $K$ , i.e.

(2.13) $$\begin{eqnarray}\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D6FC}K,\end{eqnarray}$$

where the adhesion factor $\unicode[STIX]{x1D6FC}$ is obtained from experiments – for example, Kaolin clay has $\unicode[STIX]{x1D6FC}=0.4$ (O’Loughlin et al. Reference O’Loughlin, Richardson, Randolph and Gaudin2013). Combining (2.3) with (2.8), (2.12) and (2.13), we obtain

(2.14) $$\begin{eqnarray}T(z,t)=\frac{F^{\prime }(2\unicode[STIX]{x1D70E}_{rz}F^{\prime }-\unicode[STIX]{x1D70E}_{rr}-\unicode[STIX]{x1D70E}_{zz}F^{\prime 2})}{(1+F^{\prime 2})^{3/2}}+\frac{\unicode[STIX]{x1D6FC}K}{\sqrt{1+F^{\prime 2}}},\end{eqnarray}$$

where $F^{\prime }=\text{d}F/\text{d}z$ and $T$ changes with time and location $z$ along the anchor. The first term on the right captures the contribution to $T$ from normal traction, while the second term is due to tangential traction. The total vertical force $N$ on the anchor at any time $t$ may now be obtained from (2.2):

(2.15) $$\begin{eqnarray}N(t)=2\unicode[STIX]{x03C0}\int _{z_{s}}^{L}T(z,t)F(z)\,\text{d}z,\end{eqnarray}$$

where $z=z_{s}$ locates the anchor’s tip; see figure 2(a).

To express $\unicode[STIX]{x1D748}$ , and hence $N$ , only in terms of the material parameters of the seabed and the anchor’s geometry and dynamics, we need to find the pressure  $p$ . For this, we consider linear momentum balance in the $r$ direction:

(2.16) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{rr}}{\unicode[STIX]{x2202}r}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{rz}}{\unicode[STIX]{x2202}z}+\frac{\unicode[STIX]{x1D70E}_{rr}-\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D703}\unicode[STIX]{x1D703}}}{r}=\unicode[STIX]{x1D70C}_{s}a_{r},\end{eqnarray}$$

where $\unicode[STIX]{x1D70C}_{s}$ is the density of seabed clay. We then utilize (2.7) and (2.12) to obtain

(2.17) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}r}=-\unicode[STIX]{x1D70C}_{s}\frac{R^{2}r}{4\tilde{r}^{3}}\left\{\dot{u}(t)-\frac{3z}{\tilde{r}^{2}}u(t)^{2}-\frac{R^{2}}{2\tilde{r}^{3}}u(t)^{2}\right\}-\frac{\sqrt{3}Kr}{\tilde{r}^{2}}\left(2+\frac{z}{K}\frac{\unicode[STIX]{x2202}K}{\unicode[STIX]{x2202}z}\right),\end{eqnarray}$$

where the last term is due to $K$ ’s variation with depth and we recall that $\tilde{r}=\sqrt{r^{2}+z^{2}}$ . At any time, with the anchor’s deceleration $\dot{u}$ provisionally known from (2.1), the above equation is integrated along constant $z$ -lines from the anchor’s surface at $r=F(z)$ to the plastic boundary at $r=r_{p}(z)$ – beyond which the clay remains solid and does not flow – to obtain the pressure $p_{a}:=p|_{r=F(z)}$ on the anchor’s surface. Utilizing (2.10), we find

(2.18) $$\begin{eqnarray}p_{a}=\unicode[STIX]{x1D70C}_{s}g(h-z)+\left[-\frac{\unicode[STIX]{x1D70C}_{s}R^{2}\dot{u}}{4\tilde{r}}+\unicode[STIX]{x1D70C}_{s}\left(\frac{zR^{2}}{4\tilde{r}^{3}}+\frac{R^{4}}{32\tilde{r}^{4}}\right)u^{2}+\sqrt{3}K\frac{2h-3z}{h-z}\ln \tilde{r}\right]_{r=F(z)}^{r=r_{p}(z)},\end{eqnarray}$$

where we have taken the far-field pressure to be lithostatic, i.e. due to the seabed’s weight, $z=h$ is the location of the sea floor with respect to the anchor (figure 2 a) and $[f(r)]_{r=F(z)}^{r=r_{p}(z)}=f(r_{p})-f(F)$ for any function $f(r)$ . The first term in the brackets in (2.18) contributes to the added mass arising from the clay’s inertia while the second leads to the coefficient $\unicode[STIX]{x1D6FD}_{1}$ in (1.1). Finally, the sea’s weight acts on both the seabed and the anchor and may, therefore, be ignored.

To utilize (2.18) the location $r_{p}(z)$ of the plastic boundary has to be estimated. By definition the seabed clay does not yield for $r>r_{p}$ . An exact estimate for the plastic boundary requires the solution for the equilibrium of the unyielded material. We do not pursue this complex plasticity calculation, given the relative simplicity of our overall description, and because material behaviour and system geometry are rather ill-constrained in actual scenarios. Instead, we note from (2.6) that the strain-rate’s magnitude drops off as $r^{-3}$ for large  $r$ , so that we expect the seabed to be relatively undisturbed not too far from the anchor. We then set $r_{p}/R\approx 4$ , a value which is consistent with the large-deformation, elasto-plastic computations of Sabetamal et al. (Reference Sabetamal, Carter, Nazem and Sloan2016), and with estimates derived from the analytical solution of the dynamical expansion of a spherical cavity in pressure-dependent materials (Durban & Masri Reference Durban and Masri2004).