Dynamic enlargement of a hole in a sheet: crater formation and propagation of cylindrical shock waves

TL;DR

This paper provides a complete theoretical solution for the dynamic expansion of a cylindrical hole in a strain hardening elastoplastic sheet, revealing shock wave formation at high velocities and implications for protective material design.

Contribution

It introduces the first comprehensive theoretical model of dynamic crater formation in plane-stress conditions, including shock wave analysis and constitutive sensitivity.

Findings

Shock waves occur at lower expansion velocities in plane-stress deformation.

Material compressibility has negligible effect on shock wave formation.

The model reduces to classical solutions in the quasistatic limit.

Abstract

Predicting the shape of a crater formed by high velocity impact is of interest in several fields. It can aid in design of more efficient protective structures, in forensic analysis of bullet holes, and in understanding the effects of meteorite impact in both space systems and in extreme geological events. In this paper we present, for the first time, a complete theoretical solution of the dynamic plane-stress problem. We consider the steady-state expansion of a cylindrical hole in a strain hardening elastoplastic sheet and find that a self-similar field emerges if the `specific cavitation energy' is constant. It is shown that at the quasistatic limit this solution reduces to available classical solutions, while at high expansion velocities shock waves can appear. Investigation of the constitutive sensitivities of the expansion field is conducted and compared with available results for the spherical field which is commonly applied to predict resistance to high velocity penetration. It is shown that shock waves appear at significantly lower expansion velocities, in the plane-stress deformation pattern, for which material compressibility is found to have a negligible effect. This insensitivity can be taken advantage of in the future for design of light weight protective layers by incorporating porosity.

Tables

Table 1. Table 1: Material properties of investigated representative materials. Material 1 serves as a reference material. The critical (dimensionless) velocity at which a shock wave first appears is denoted by m c subscript 𝑚 𝑐 m_{c} . Recall that σ y subscript 𝜎 𝑦 \sigma_{y} is dimensionless with respect to the elastic modulus.

Material	$\nu$	$n$	${\sigma }_y$	$m_c$
1	0.3	0.3	0.003	0.15
2	0.3	0.1	0.003	0.08
3	0.1	0.3	0.003	0.15
4	0.5	0.3	0.003	0.15
5	0.3	0.3	0.001	0.11

Equations

$$\xi =\frac{r}{a}.$$

$$V=\frac{v}{\dot{a}}.$$

$$\frac{d()}{dr}=({)}^{\prime }\frac{1}{a},\dot{()}=({)}^{\prime }(V-\xi )\frac{\dot{a}}{a}\text{where}({)}^{\prime }=\frac{d()}{d\xi },$$

$$\frac{\partial (h{\bar{\sigma }}_r)}{\partial r}+\frac{h}{r}({\bar{\sigma }}_r-{\bar{\sigma }}_{\theta })=\frac{\partial (\rho hv)}{\partial t}+\frac{1}{r}\frac{\partial (r\rho h{v}^2)}{\partial r}.$$

$$\frac{\partial (\rho h)}{\partial t}+\frac{1}{r}\frac{\partial (r\rho hv)}{\partial r}=0.$$

$$\frac{\partial (h{\sigma }_r)}{\partial r}+\frac{h}{r}({\sigma }_r-{\sigma }_{\theta })=\frac{\rho }{E}h\dot{v}.$$

$$(h{\sigma }_r{)}^{\prime }+\frac{h}{\xi }({\sigma }_r-{\sigma }_{\theta })=m^2\frac{\rho }{{\rho }_0}h(V-\xi )V^{\prime },$$

$$V^{\prime }+\frac{V}{\xi }+(V-\xi )\left(\frac{h^{\prime }}{h}+\frac{{\rho }^{\prime }}{\rho }\right)=0,$$

$$m=\frac{\dot{a}}{C_E}.$$

$${\sigma }_e^2={\sigma }_r^2-{\sigma }_r{\sigma }_{\theta }+{\sigma }_{\theta }^2.$$

$$\{\begin{array}{ll}{\dot{\epsilon }}_r& ={\dot{\sigma }}_r-\nu {\dot{\sigma }}_{\theta }+{\eta }_r{\dot{\epsilon }}_p,\\ {\dot{\epsilon }}_{\theta }& ={\dot{\sigma }}_{\theta }-\nu {\dot{\sigma }}_r+{\eta }_{\theta }{\dot{\epsilon }}_p,\\ {\dot{\epsilon }}_z& =-\nu ({\dot{\sigma }}_r+{\dot{\sigma }}_{\theta })-({\eta }_r+{\eta }_{\theta }){\dot{\epsilon }}_p,\end{array}$$

$${\eta }_r=\frac{2{\sigma }_r-{\sigma }_{\theta }}{2{\sigma }_e},{\eta }_{\theta }=\frac{2{\sigma }_{\theta }-{\sigma }_r}{2{\sigma }_e}.$$

$${\epsilon }_r=\ln \left(\frac{\partial r}{\partial R}\right),{\epsilon }_{\theta }=\ln \left(\frac{r}{R}\right),{\epsilon }_z=\ln \left(\frac{h}{h_0}\right).$$

$${\epsilon }_r+{\epsilon }_{\theta }+{\epsilon }_z=\ln \left(\frac{{\rho }_0}{\rho }\right).$$

$${\epsilon }_r+{\epsilon }_{\theta }+{\epsilon }_z=3\beta {\sigma }_h,$$

$$\frac{\rho }{{\rho }_0}=e^{-3\beta {\sigma }_h}.$$

$${\dot{\epsilon }}_r=\frac{dv}{dr},{\dot{\epsilon }}_{\theta }=\frac{v}{r},{\dot{\epsilon }}_z=\frac{\dot{h}}{h},$$

$$V^{\prime }=(V-\xi )\left({\sigma }_r^{\prime }-\nu {\sigma }_{\theta }^{\prime }+\alpha {\eta }_r{\sigma }_e^{\prime }\right),$$

$$\frac{V}{\xi }=(V-\xi )\left({\sigma }_{\theta }^{\prime }-\nu {\sigma }_r^{\prime }+\alpha {\eta }_{\theta }{\sigma }_e^{\prime }\right),$$

$$\frac{h^{\prime }}{h}={\epsilon }_z^{\prime }=-\nu ({\sigma }_r^{\prime }+{\sigma }_{\theta }^{\prime })-\alpha ({\eta }_r+{\eta }_{\theta }){\sigma }_e^{\prime }.$$

$$\alpha =\frac{d{\epsilon }_p}{d{\sigma }_e}.$$

$${\sigma }_r^{\prime }+{\sigma }_r{\epsilon }_z^{\prime }+\frac{1}{\xi }({\sigma }_r-{\sigma }_{\theta })=m^2{e}^{-3\beta {\sigma }_h}(V-\xi )V^{\prime }.$$

$${\sigma }_e^{\prime }={\eta }_r{\sigma }_r^{\prime }+{\eta }_{\theta }{\sigma }_{\theta }^{\prime },$$

$$V^{\prime }=(V-\xi )\left((1+\alpha {\eta }_r^2){\sigma }_r^{\prime }-(\nu -\alpha {\eta }_r{\eta }_{\theta }){\sigma }_{\theta }^{\prime }\right),$$

$$\frac{V}{\xi }=(V-\xi )\left((1+\alpha {\eta }_{\theta }^2){\sigma }_{\theta }^{\prime }-(\nu -\alpha {\eta }_r{\eta }_{\theta }){\sigma }_r^{\prime }\right),$$

$${\epsilon }_z^{\prime }=-(\nu +\alpha {\eta }_r({\eta }_r+{\eta }_{\theta })){\sigma }_r^{\prime }-(\nu +\alpha {\eta }_{\theta }({\eta }_r+{\eta }_{\theta })){\sigma }_{\theta }^{\prime }.$$

$$(1-{\sigma }_r(\nu +\alpha {\eta }_r({\eta }_r+{\eta }_{\theta }))){\sigma }_r^{\prime }-{\sigma }_r(\nu +\alpha {\eta }_{\theta }({\eta }_r+{\eta }_{\theta })){\sigma }_{\theta }^{\prime }+\frac{1}{\xi }({\sigma }_r-{\sigma }_{\theta })=m^2{e}^{-3\beta {\sigma }_h}(V-\xi )V^{\prime }.$$

$${\epsilon }_p=\{\begin{array}{ll}0& {\sigma }_e<{\sigma }_y,\\ {\sigma }_y{\left(\frac{{\sigma }_e}{{\sigma }_y}\right)}^{\frac{1}{n}}-{\sigma }_e& {\sigma }_e\le {\sigma }_y,\end{array}$$

$$\begin{array}{rl}\Delta =(V-\xi )& [1-\nu (1+\nu ){\sigma }_r+\alpha ((1-2\nu {\sigma }_r){\eta }_{\theta }^2-{\sigma }_r{\eta }_r({\eta }_r+{\eta }_{\theta }))\\ & -m^2{e}^{-3\beta {\sigma }_h}(V-\xi {)}^2(1-{\nu }^2+\alpha ({\eta }_r^2+{\eta }_{\theta }^2-2\nu {\eta }_r{\eta }_{\theta }))]=0.\end{array}$$

$$m^2{e}^{-3\beta {\sigma }_h}(V-\xi {)}^2=\frac{1-\nu (1+\nu ){\sigma }_r+\alpha ((1-2\nu {\sigma }_r){\eta }_{\theta }^2-{\sigma }_r{\eta }_r({\eta }_r+{\eta }_{\theta }))}{1-{\nu }^2+\alpha ({\eta }_r^2+{\eta }_{\theta }^2-2\nu {\eta }_r{\eta }_{\theta })}.$$

Full text

1 Introduction

In 1947 G. I.  Taylor reports on an investigation on the enlargement of a circular hole in a thin plastic sheet (Taylor, 1948). Considering a perfectly plastic material, Taylor combines the Tresca yield condition with the Mises flow rule to arrive at what we would refer to today as a non-associated plasticity model, and studies the quasistatic expansion process. He compares his results with those from an unpublished communication with Hance A. Bethe111This communication is apparently an extension of earlier work by Bethe (1941). and concludes that Bethe’s model is “inconsistent with any theory of plasticity”. Later, R.  Hill (1950) considers the same problem and subjects Taylor’s work to similar critisism. He claims that “Taylor’s method of integration was insufficiently accurate” and that “the almost exact agreement… found experimentally by Taylor is certainly fortuitous”. All three authors considered situations in which the hole expansion process occurs symmetrically under plane stress conditions, such that “the metal near the hole piles up into a thickened crater”. To examine if the symmetry assumption holds, Taylor conducted experiments in lead sheets. The expansion was achieved by a series of tapering sections of steel cones that penetrated the sheet while being rotated. Images of two representative samples where shown (Taylor, 1948), and have been included in Appendix A of this paper.

The controversy between Bethe, Taylor, and Hill is revisited in more detail by Masri et al. (2010), and resolved by employing modern constitutive models and integration methods. It revolves around the values of the radial pressure applied at the hole $(p)$, and the deformed thickness of the ‘crater’ $(h)$. These two values combine to predict the work invested in the quasistatic expansion process, which is proportional to their product $(ph)$. More recently, Cohen et al. have shown that, in fact, both $p$ and $h$ can become singular at the cavity wall, while their multiplication arrives at an asymptotic value (Cohen and Durban, 2013a, 2010; Cohen et al., 2009). Accordingly, they define the ‘specific cavitation energy’ - $s_{c}=ph/h_{0}$, with $h_{0}$ representing the initial thickness of the sheet. This asymptotic quantity can be interpreted as the incremental work needed to create a new unit of hole volume, during self-similar expansion. It provides a measure for the resistance of a solid to penetration.

Motivated by the wartime need for prediction of penetration resistance of protective metal targets, the studies by Bethe, Taylor, and Hill, provided a meaningful estimation of the energy invested in the perforation process. However, their results did not account for the dynamics of the process and, in particular, the role of inertia in resisting high velocity perforation. In fact, it is during the same years that the first theories of dynamic propagation of plastic deformation were being formulated. Perhaps most notable is the work of von Karman and Duwez (1950) that considered tensile impact of a bar. Quite surprisingly, to this day, the fundamental problem and the quasistatic solutions presented by Bethe, Taylor, and Hill, have not been extended to dynamic expansion. This is despite the fact that its solution can extend to explain several additional phenomena including extreme geophysical events, such as meteorite impact (Melosh and Collins, 2005; Goodier, 1965), volcanic eruptions and earthquakes (Gudmundsson, 2016, 2014), and the growing threat of orbital debris on space operations (Johnson, 2010; Smirnov and Kondratyev, 2009; Hopkins and Swift, 1970; Hosseini and Abbas, 2006).

In absence of solutions for the dynamic expansion problem in plane-stress conditions, Goodier (1965), followed the notion of Bishop et al. (1945), to determine the indentation resistance of metal plates by employing the spherically symmetric cavity expansion field. To account for inertia, he used solutions by Hill (1950) and Hopkins (1960). This methodology has been shown to provide good prediction for deep penetration problems in various materials (Forrestal et al., 1995; Forrestal and Tzou, 1997; Masri, 2010; Chen and Li, 2002; Ben-Dor et al., 2005; Warren and Forrestal, 1998; Johnsen et al., 2018; Vorobiev et al., 2007; Masri and Durban, 2009; Gabi et al., 2013), and over the years has developed into an active research discipline. For thin sheets, other studies continue to employ the quasistatic plane-stress expansion field pioneered by Bethe, Taylor, and Hill, as recently reviewed by Ryan et al. (2018). In particular, models based on the ‘specific cavitation energy’ (Cohen et al., 2010a; Masri, 2014, 2015; Cohen and Durban, 2013a, 2010; Cohen et al., 2009), have been identified as providing the best performance.

In this paper we present, for the first time, a complete solution of the dynamic plane-stress problem, which reduces to the solutions by Bethe, Taylor, and Hill at the quasistatic limit. This solution is achieved by identifying that when a constant ‘specific cavitation energy’ applied to the sample, is larger than the quasistatic value $(s_{d}>s_{c})$, the cylindrical hole will expand in a self-similar manner at a constant velocity. Additionally, it is shown that if the expansion velocity exceeds a critical value, a cylindrical shock wave will form. Solutions of the entire field, including the shock discontinuity are obtained and sensitivity to the constitutive properties (i.e. yield stress, elastic compressibility, and strain hardening) are examined. In the next section we begin by presenting the problem setting. Then, in section 3, the governing equations are derived. In section 4, the solution procedure is presented, including derivation of an analytical solution for the elastic zone. In section 5 we present results and provide a discussion in comparison with available results for the quasistatic field, and for the dynamic spherical cavity expansion field. Finally, we conclude in section 6, and discuss directions for future work.

2 Problem setting and the self-similar field

Consider a cylindrical cavity embedded in a plane sheet being expanded radially at constant expansion velocity $\dot{a}$ by application of internal (dimensionless) pressure $p$, as illustrated in Fig. 1. The sheet, of undeformed thickness $h_{0}$, spans indefinitely in the plane, hence disturbances are not reflected back towards the cavity. It will be shown that in this setting, assuming plane-stress conditions and permitting changes in thickness, a self-similar field emerges in which the ‘specific cavitation energy’ (Cohen and Durban, 2013a, 2010; Cohen et al., 2009) invested in expanding the cavity remains constant in time, although the applied pressure and deformed thickness at the cavity wall become singular.

Assuming symmetry about the axis of the cylindrical cavity, the material velocity $v$ is in the radial direction and we denote the current and reference radial coordinates by $r$ and $R$, respectively. A self-similar field must be independent of the current radius of the cavity $a$, and thus, if it exists, all field variables can be written as functions of a single independent variable, the dimensionless radial coordinate

[TABLE]

As was previously suggested for the spherical field in Durban and Fleck (1997) we may define the dimensionless velocity along the radial coordinate by

[TABLE]

Hence, at the cavity wall $(r=a)$, we have both $\xi=1$ and $V=1$. We can now write the useful transformations

[TABLE]

for the radial and time differentiation of a field variable in the self-similar state of cavity expansion, respectively.

As shown for the spherical cavity expansion field (Cohen et al., 2010b; Cohen and Durban, 2013b; Durban and Masri, 2004), for materials with a definite yield point, in absence of remote loading, we expect that at the steady-state an elastic precursor, at $\xi_{e}$, propagates into the undeformed stationary regime, followed by an elastoplastic interface, at $\xi_{i}$, that separates the elastic region from the elastoplastic region, in which the material has yielded. It was shown in (Cohen et al., 2010b; Cohen and Durban, 2013b) that in the spherical setting, there exists a critical expansion velocity $\dot{a}$, at which a plastic shock wave may appear $(\xi_{p}$). In the present study we investigate the analogous field in a plane-stress deformation pattern, hence we have an additional field variable; the current thickness of the sheet, denoted by $h$. The propagating singularity surface, identified as a shock wave, implies a finite jump in field variables and as such induces a localized jump in the thickness of the sheet. A similar geometrical discontinuity was obtained by Knowles (2002), in the investigation of shock wave propagation due to tensile impact of rubber.

In the next sections we will show that a self-similar cavity expansion field exists in the plane-stressed deformation pattern. It is associated with a constant level of energy input invested in hole expansion (i.e. specific cavitation energy), while the applied radial pressure at the cavity wall is not finite, in contrast to the extensively studied spherical pattern.

3 Governing equations

Continuity. We begin by writing the continuity equations in the most general form. This will be useful later in deriving the jump conditions across a shock discontinuity.

In the plane-stress deformation pattern, though changes in thickness of the sheet are permitted, motion is assumed to occur only in the plane. Conservation of momentum in the axially symmetric field is thus written for radial motion

[TABLE]

which employs the underlying assumption, that the field is uniform through the thickness of the sheet. Here $\rho$ is the current density and $(\bar{\sigma}_{r},\bar{\sigma}_{\theta})$ are, respectively, the radial and circumferential principal stress components. In the following formulation we will remove the bar notation to denote the dimensionless values, henceforth all stress components in the formulation have been nondimensionalized with respect to the elastic modulus $E.$

The second continuity equation requires conservation of mass

[TABLE]

By combining equations (4) and (5) we may rewrite the equation of motion in the more compact Navier-Stokes form

[TABLE]

We now insert relations (1) - (3) into (5) and (6) to arrive at the self-similar versions of the equation of motion

[TABLE]

and conservation of mass

[TABLE]

respectively. Here we have used the dimensionless expansion velocity

[TABLE]

where $C_{E}=\sqrt{E/\rho_{0}}$ is the wave velocity in a long elastic rod and $\rho_{0}$ is the density of the material in its undeformed state.

Constitutive response. Following the standard principles of Mises flow theory plasticity, we write the dimensionless Mises effective stress in the present plane-stress field as

[TABLE]

Employing an associated flow rule, we take the total logarithmic strain rate as the sum of a Hookean hypoelastic part and a plastic part, and invoke the principle of plastic power equivalence, which for the present deformation pattern, leads to a set of three scalar constitutive equations

[TABLE]

where, for compactness, we have defined

[TABLE]

Here $(\varepsilon_{r},\varepsilon_{\theta},\varepsilon_{z})$ are the radial, circumferential and out of plane logarithmic strain components, respectively, $\nu$ is the Poisson’s ratio, and $\varepsilon_{p}$ is the effective (logarithmic) plastic strain; a known function of the effective stress. We limit our analysis to rate independent material response and neglect thermal effects. Hence, at this point in the formulation it is most instructive to proceed with an arbitrary hardening relation $\varepsilon_{p}=\varepsilon_{p}(\sigma_{e})$, which can be determined by the standard tension test for a material of interest.

Strain kinematics. The logarithmic strain components in the present deformation pattern can be written as

[TABLE]

and thus the dilatation is

[TABLE]

which is an integrated form of the continuity equation (8). On the other hand, by adding relations (11) and performing straight forward integration over time we find that

[TABLE]

where the plastic branch does not participate due to plastic incompressibility. Here $\sigma_{h}=(\sigma_{r}+\sigma_{\theta})/3$ is the hydrostatic stress, and $\beta=1-2\nu$ is the compressibility parameter. Combining the above two relations yields the useful result which is independent of the loading path

[TABLE]

By differentiating the strain components (13) we have the strain rates

[TABLE]

which, in combination with the self-similar field relations (1)-(3), can be inserted into (11) to arrive at the self-similar form of the constitutive relations

[TABLE]

[TABLE]

[TABLE]

Note that since $\varepsilon_{p}=\varepsilon_{p}(\sigma_{e})$ we have replaced $\varepsilon^{\prime}_{p}=\alpha\sigma^{\prime}_{e}$ where

[TABLE]

We can now rewrite the equation of motion (7) with the aid of (16) and (20) as

[TABLE]

Summary of equations. At this point we have a system of four first-order nonlinear differential equations (18)-(20) and (22) with the five unknown field variables $(\sigma_{r},\sigma_{\theta},\sigma_{e},\varepsilon_{z},V)$. An additional relation is provided by definition of the effective Mises stress (10), which implies also a relation between the plastic strain and the principal stress components. It is possible to reduce this system by removing the dependence on $\sigma_{e}$ and $\varepsilon_{z}$. First, by differentiating equation (10) we can write

[TABLE]

which is substituted into equations (18)-(20) to arrive at the forms

[TABLE]

[TABLE]

[TABLE]

Now, substituting (26) in (22) reads

[TABLE]

In summary, we have a set of three differential equations (24), (25) and (27) with three unknown field variables $(\sigma_{r},\sigma_{\theta},V)$ that depend on the single independent variable $\xi$. Recall that $\eta_{r},\eta_{\theta}$ are known functions of $\sigma_{r},\sigma_{\theta}$ defined in (12) and with the effective stress given in (10). Once a solution is obtained, the density ratio $\rho/\rho_{0}$ is given from (16), and the thickness ratio $h/h_{0}$ can be obtained by integration of (26) with the definition in (13)3. By (21), $\alpha$ is assumed to be a known function of $\sigma_{e}$. To represent a broad range of plastic response, in the present investigation we will employ the power law

[TABLE]

where $\sigma_{y}$ is the dimensionless yield stress and $n$ is the hardening index.

Shock Condition. The governing system of equations (22) and (24)-(26) may become singular if the determinant of its coefficients vanishes, which after some algebra reads

[TABLE]

From the above result, it is first noticed that singularity necessarily appears at the cavity wall where $V=\xi$. This singularity is associated with the unbounded levels of stress that exist even at the quasistatic limit (Cohen and Durban, 2013a, 2010; Cohen et al., 2009). Nonetheless, it appears only at the boundary of the body and thus does not involve a discontinuity. On the other hand, the second root of the above equation can occur within the material region. Accordingly, we can write a shock condition in the form

[TABLE]

If this singularity appears, the local representation of the governing equations is no longer sufficient and jump conditions are required.

Jump Conditions. If a shock wave appears, in accordance with the shock condition (30), field variables may vary discontinuously across the propagating front where the governing equations are replaced with jump conditions that assure conservation of mass and momentum. If the field is self similar, a shock wave propagates at a constant dimensionless radial location denoted here by $\xi_{w}=r_{w}/a$ and with constant velocity,

[TABLE]

according to (3). Considering a frame of reference that moves with the discontinuity, in the steady-state field, we may write the transformation $\delta\xi=\xi-\xi_{w}$, such that the discontinuity is located at $\delta\xi=0$, and any field variable $g$ can be rewritten in the form $g(r,t)=\hat{g}(\xi-\xi_{w})$. Then, by substituting this transformation into the continuity equations (4) and (5) and performing integration between $\delta\xi=0^{-}$ and $\delta\xi=0^{+}$ we arrive at the jump condition

[TABLE]

for continuity of momentum, and

[TABLE]

for continuity of mass. The above two relations are plane-stress forms of the Rankine-Hugoniot jump conditions. Both are readily written in terms of the dimensionless radial coordinate and velocity of the self-similar field, and the square brackets denote the jump in the of the enclosed quantity across the discontinuity (\left\llbracket g\right\rrbracket=g^{+}-g^{-}).

Since present system consists of four differential equations (24)-(27) with four unknowns $(\sigma_{r},\sigma_{\theta},\varepsilon_{z},V)$ , two additional jump conditions are needed to fully define the jump across the singularity. First we identify that to permit discontinuity, singularity of the field equations must hold on both sides of the shock. According to the shock condition (30), this translates to the jump condition

[TABLE]

An additional requirement that must be accounted for is compatibility, which for the present problem translates to continuity of the circumferential strain. Hence, we return to the constitutive relation (11)2, which we now write in incremental form, and require $\rm d\varepsilon_{\theta}=0$ across the shock to write

[TABLE]

With the aid of (10), the above result can be rewritten in the form

[TABLE]

Now, by rearranging, we find that variations of stress throughout plastic deformation at constant circumferential strain, obey the first order nonlinear differential equation

[TABLE]

Hence, for any given initial state (in front of the shock) represented by the pair $(\sigma^{+}_{r},\sigma^{+}_{\theta})$, isothermal compression, by variation of $\varepsilon_{r}$ at a constant $\varepsilon_{\theta}$, follows a unique path up to the final state (behind the shock), represented by the pair $(\sigma^{-}_{r},\sigma^{-}_{\theta})$, which can be obtained by straight forward integration. Although departures from the isothermal compression curve, defined by (37), are permitted within the shock, and lead to dissipation, in absence of rate dependence and thermo-mechanical coupling, the end states are assumed to exist on the curve. This provides us with a correspondence condition between the two states and thus dictates permissible jumps across the shock discontinuity. In connection with available shock wave theories (Davison and Graham, 1979), the above relation essentially provides a differential analog to the Hugoniot curve.

Boundary Conditions. Since the remote field is undisturbed and stationary, we may write the remote boundary conditions as

[TABLE]

At the cavity wall we prescribe the expansion velocity, which by definition implies

[TABLE]

which serves as a compatibility requirement.

For a given expansion velocity $m$, a physical self-similar solution exists if it is associated with a finite level of applied load, which in the present field, is identified with the dynamic specific cavitation energy

[TABLE]

Though the boundary conditions (38) and (39) do not seem to depend on the expansion velocity, it is the singularity of the equations at the remote field $(\xi=\xi_{e})$ that inserts an additional unknown, which is directly related to $m$. In the next sections we present the solution procedure, beginning with an analytical solution of the remote elastic zone, which appears if the material has a definite yield point.

4 Solution procedure

4.1 The elastic zone

In the elastic zone the plastic strain identically vanishes, thus $\alpha\equiv 0$ and equations (24)-(26) reduce to

[TABLE]

[TABLE]

[TABLE]

Upon subtracting equation (42) from (41) it is possible to perform integration and after some algebra to arrive at the relation for the material velocity

[TABLE]

this is inserted back into (42) to write

[TABLE]

Relations (41), (43) and (44) can now be inserted into the equation of motion (22) to write

[TABLE]

A further simplification is achieved by recognizing that in the elastic zone the dimensionless stress components are small compared with unity $|\sigma_{r}|,|\sigma_{\theta}|\ll 1$, and as previously suggested by Durban and Masri (2004) for the spherical field, we may write the linearized form of the above two relations

[TABLE]

[TABLE]

The determinant of coefficients of the above system is

[TABLE]

Notice that, according to boundary condition (38), in the remote field the stress components vanish and the system of equations is homogenous. Hence, a nontrivial solution at $\xi=\xi_{e}$ can be obtained only if the determinant vanishes $(\Delta=0)$, thus further implying that the radial location of the elastic precursor is given by

[TABLE]

and since in the dynamic field $\xi<\xi_{e}$ we consistently have $\Delta^{2}\geqslant 0$. It is worth noting that, according to relation (3)2 in combination with (9), a wave at $\xi=\xi_{w}$ propagates with the velocity

[TABLE]

which by (50) implies that cylindrical waves in a thin linearly elastic plate propagate faster than longitudinal waves in a long elastic bar $(v_{e}/C_{E}=1/\sqrt{1-\nu^{2}}>1)$. Additionally, we find that the elastic wave propagates at finite velocity, even if the material is incompressible $(\nu=1/2)$. This is due to the effective compressibility induced by the out of plane deformation. In the next section it will be shown that this effective compressibilty dominates the response, such that the effect of the elastic compressiblity becomes negligible.

We now proceed to solve the linear system of equations (47) and (48) to arrive at the two relations for the radial and circumferential stress components

[TABLE]

[TABLE]

where we have readily applied the stress free boundary condition (38) at $\xi=\xi_{e}$ (with $\Delta=0$). The remaining integration constant $c$ must be obtained by solving the entire field, to comply with the compatibility condition at the cavity wall (39) for a given expansion velocity $m$.

Notice that the principal stress difference is

[TABLE]

Hence, in combination with (44), we find that the cavity boundary conditions $(V=\xi=1)$ cannot be satisfied assuming a purely linearly elastic response. In other words, a linearly elastic solution to the self-similar expansion problem does not exist.

Now, we may obtain the out of plane component of strain by integration of (43), which upon substitution of the above stress components along with the boundary condition (38) reads

[TABLE]

Similarly, at the linear limit, the velocity field (44) reduces to

[TABLE]

Yield. The elastic zone terminates once the effective stress arrives at the yielding limit, $\sigma_{e}=\sigma_{y}$. By assigning $\sigma_{e}=\sigma_{y}$ in (10) and substituting the principal stress components with the relations (52) and (53), we can find the location of the elastoplastic interface $\xi_{i}$

[TABLE]

where $\Delta_{i}=\Delta(\xi_{i})$.

In the range $1\leqslant\xi<\xi_{i}$, the material flows plastically. Due to the nonlinearity in this regime, we shall solve the governing system of equations in the elastoplastic zone by numerical integration, as will be discussed in the following section. It will be shown that under certain conditions shock waves may appear in this zone and applying jump conditions across the discontinuity is imperative to obtain a solution of the entire field.

4.2 The elastoplastic zone

Once plasticity sets in, the solution proceeds by straight forward numerical integration222Integration is performed using a forth-order Runge-Kutta method. of the system of four first order nonlinear differential equations (24)-(27), to find the four unknown field variables $(\sigma_{r},\sigma_{\theta},\varepsilon_{z},V)$. This is done by applying a shooting method such that for a given expansion velocity $(m)$ the integration constant $c$, in the elastic zone equations (52)-(57), is chosen so that the compatibility condition at the cavity wall (39) is satisfied. If the cavity expansion rate is sufficiently high, then the governing system of equations may become singular at a given radial location $\xi=\xi_{p}$ within the elastoplastic zone. That cylindrical singularity surface is identified as a shock wave. If a shock wave appears, then jump conditions (32)-(34) are applied, along with the numerical integration of (37), to obtain the isothermal compression curve, to resolve the entire field.

5 Results and Discussion

Following the solution procedure described in the previous section, results are obtained for self-similar expansion in sheets of different constitutive properties and at various expansion velocities. Table 1 summarizes the set of representative material properties used to investigate the constitutive sensitivities of the expansion dynamics and the critical velocity at which a shock wave first appears $(m_{c})$.

The self-similar fields that develop in sheets of Materials 1, 2 and 5, subjected to hole expansion at different rates (as represented by $m$) are shown in Figure 2. For all materials, it is immediately observed that as expansion velocities increase, the field concentrates near the cavity wall $(\xi=1)$ and the gradients steepen. Beyond a critical expansion velocity $(m_{c}$) a discontinuity is observed. By comparing fields for Materials 1 and 2, (with $n=0.3$ and $n=0.1$, respectively) it is found that material hardening also contributes to the localization of the response near the cavity wall. Hence, the critical velocity for appearance of a shock is highly sensitive to the hardening index (Table 1).

As the quasistatic limit is approached $(m\to 0)$, we recover the results obtained in earlier studies (Cohen et al., 2009; Masri et al., 2010). Therein, it was shown that for perfectly plastic materials (with $n\to 0$), the solutions obtained by Taylor (1948) and Hill (1950) are recovered. Those results are comparable to the black curves for Material 2 in Figure 2, and exhibit non-intuitive trends. Both the radial and the circumferential stresses are shown to be non-monotonic. The radial stress decreases in magnitude as the cavity wall is approached, while the circumferential stress changes signs throughout the field. Similar non-monotonic behaviors are observed for all materials and can be explained by the particular deformation pattern in which circumferential expansion and radial compression compete in promoting hole expansion. Nonetheless, in all cases, the effective stress (derived from Eq. (10)) varies monotonically, as shown in Figure A.2 of the Appendix section. Additionally, it is observed that the principal stress components and the effective stress, all approach a singularity at the cavity wall.

Next, we examine the variations in thickness for the different fields. For the material with lower hardening (Material 2) we find that the edge of the ‘crater’ is a thin knife edge, as described by Taylor (1948). That sharp increase in thickness at the cavity wall appears also for the hardening materials (Material 1 and 5), however only at high expansion velocities. For moderate expansion velocities a peak value is observed (as shown also in (Cohen et al., 2009; Masri et al., 2010) for the quasistatic limit). This peak value is accompanied by change in slope of the radial stress, as shown in the enlarged region for Material 1.

Variations in density are shown in Figure A.2 of the Appendix section. Quite surprisingly, it is found that over a broad range of parameters, elastic compressiblity leads to negligible changes in density of the order of $1\%$. To further examine this nearly isochoric deformation, we compare results obtained for different values of the Poisson’s ratio $\nu$ (i.e. Materials 1, 3 and 4) and we find a negligible effect on the resulting field. This insensitivity is shown in Figure 2 by comparing the thickness profiles for Material 1 (with $\nu=0.3$) and Material 4 (with $\nu=0.5$). Discrepancies between the responses of the two materials are observed only in near the edge of the crater. It is found that the peak value of thickness does not appear for the incompressible response. Overall, it is apparent that the thickening of the plate is a preferable mechanism for cavity expansion, in comparison with volume change, and the effect of material compressibility is localized to the vicinity of the crater333Note that similar insensitivity to Poisson’s ratio was obtained in all field variables and for $0<\nu<0.5$. The comparison in Figure 2 (shown by the yellow dashed lines) serves as a representative example. .

After examining the expansion fields for different representative materials, and identifying the singularity that arises at $\xi=1$, we now examine in Figure 3 the radial variation of the nominal pressure

[TABLE]

It is shown that in all cases the nominal pressure arrives at a finite value when the edge of the crater is approached. According to equation (40), this value is the ‘specific cavitation energy’ - $s_{d}=s(\xi\to 1)$.

To further investigate the dependence of the specific cavitation energy on the expansion velocity, we plot curves of $s_{d}=s_{d}(m)$ in Figure 4. Here we further confirm the insensitivity of the expansion field to the elastic compressibility, by observing that the curves for Materials 1, 3 and 4 are indistinguishable. Curves for Materials 2 and 3, that have lower hardening index and yield stress, respectively, expand at lower values of $s_{d}$, with difference more pronounced at the quasistatic limit (an inset with an enlarged view of this limit is shown). The curve for Material 2 is limited to values of $m\leq 0.15$. This limitation is due to numerical difficulties that arise in the high gradient zone along with high sensitivity to the value of integration constant $c$, at high expansion velocities444Recall that $c$ arises in the solution of the elastic field (52)-(57).. This increased sensitivity is observed in Figure 4, by the formation of an asymptotic limit for the integration constant $c$. Essentially, at these extreme expansion velocities, small perturbations in the far field response, as dictated by $c$, are associated with large changes in the near field response.

To compare the results obtained for plane-stress expansion with the commonly applied spherical field results, we have plotted in Figure 4 the ‘specific cavitation energy’ for the spherical field (from Cohen et al. (2010b))555Note that for the spherical field the specific cavitation energy is identical to the cavitation pressure (i.e. $s_{d}=p$). Find comparison between quasistatic response in Cohen and Durban (2013a). with $n=0.3$, $\sigma_{y}=0.003$ and with two extreme values of the Poisson’s ratio $\nu=0$ and $\nu=0.5$. It is found that while the incompressible limit provides a strict upper bound in the entire range, spherical expansion with lower values of $\nu$ can provide a good approximation in an intermediate range of expansion velocities. Nonetheless, in the spherical field the appearance of a shock is delayed to higher expansion velocities in the range $\ 0.4<m_{c}<0.7$ in comparison to $0.08\leq m_{c}\leq 0.15$ the present plane-stress field, for representative Materials in Table 1.

6 Conclusions

High velocity expansion of a cylindrical hole in a sheet is investigated in the context of large strain $J_{2}$ plasticity, accounting for strain hardening and elastic compressibility in a plane-stress deformation pattern. While earlier studies have been limited to quasistatic expansion, or adopt dynamic results obtained for spherically symmetric deformation, this study reveals that self-similar expansion does exist in the plane-stress field and captures the shape of the crater that is formed and its dependence on both the expansion velocity and the constitutive properties of the penetrated material. It is found that a steady-state is reached when the energy invested in creating new hole volume - the ‘specific cavitation energy’, is constant. Solutions of the entire field are obtained, including jump conditions that are applied across a shock discontinuity, if it appears. Study of the constitutive sensitivities of the resulting field reveals the influence of material hardening and yield stress on the dynamic expansion, and shows that the effective compressibility, induced by out of plane deformation, dominates over the effect of elastic compressibility of the material. At the quasistatic limit, results are shown to agree with available solutions. In the dynamic range, comparison of the penetration resistance (the specific cavitation energy) obtained for the present field, with available results for the spherical field, shows significant geometric sensitivity; the spherical field is highly sensitive to elastic compressibility and exhibits shock waves at expansion velocities that are of order of magnitude higher.

There are a number of limitations in the present analysis, which give rise to a multitude of additional open questions that should be subject for future work. The present formulation is restricted to in plane motion and does not account for out of plane bending or plug formation. Thermo-mechanical coupling is not accounted for, nor is strain rate dependence. Incorporating the later effect can potentially be achieved by applying the methodology used in (Santos et al., 2019), for the spherical field. Thermal effects can be accounted for by extending the work in Masri (2014) to include inertia. The insensitivity to material compressiblity, identified in this work, suggests an avenue for exploration of the dynamic performance of light weight protective layers, by incorporating porosity, as in Cohen and Durban (2012, 2013a, 2013b); Czarnota et al. (2017). Additional constitutive models can also be explored, to study the response of reinforced concrete or geological materials (Vorobiev et al., 2007; Durban and Papanastasiou, 1997), or the effect of mechanically induce phase transformations (Rodríguez-Martínez et al., 2014).

By providing prediction of the geometry of the crater and its dependence on the penetration velocity, the results obtained in the present study can be used for forensic analysis of projectile penetration in ductile metals. The dependence of penetration resistance on expansion velocity and the constitutive properties of the target material can inform the design of more efficient protective layers.

Acknowledgement

The author wishes to acknowledge the support of the Army Research Office and Dr Ralph A. Anthenien, Program Manager, under award number W911NF-19-1-0275.

Appendix