Abstract

The impact of a flat-ended cylindrical rod onto a rigid stationary anvil, often known as the Taylor impact test, is studied. An axisymmetric model is developed to capture the deformation behaviour of the rod after impact. The most distinctive feature of the proposed model is that it takes into account the spatial and temporal variation of both longitudinal and radial deformation and consequently the strains and strain rates. The final deformed shapes and time histories of different field variables, as obtained from the model, are found to be in good agreement with corresponding experimental and numerical results reported in the literature. The proposed model is then used to formulate an inverse framework to estimate the Johnson–Cook constitutive parameters. In the inverse formulation, the objective function is constructed using the final deformed length and diameter at the impact end of the retrieved rod. Finally, the potential of the proposed model in estimating material parameters is illustrated through some examples.

1. Introduction

The impact of a cylindrical deformable rod on a rigid stationary surface is a classical problem in impact mechanics and has been the subject of interest for many years to both experimentalists and numerical analysts. This problem is frequently referred to as the Taylor impact problem in view of the pioneering work by Taylor [1]. His findings with those of Whiffen [2], Carrington & Gayler [3], which curved an elegant approach to estimate the plastic flow stress, still bears significance as a base line for understanding physical principles behind ductile metal behaviour under high strain rate. Their analysis uses data from recovered samples after several experiments on a flat-ended cylindrical rod impacting on a rigid massive anvil. The simplistic mathematical foundation made the test a reference point thereafter to estimate dynamic yield strength, as a function of density, impact velocity, original and final length, and length of the un-deformed region of the retrieved specimen and also as a validation test case for constitutive model development [4,5].

Since then, several attempts have been made to develop correlation between rate-dependent material constants at a high strain rate and the deformation characteristic of the Taylor specimen. Material constants were determined by feeding the experimentally measured deformation parameters (final deformed length, length of plastic zone, etc.) into the developed relation. Hawkyard et al. [6] equated kinetic energy in the impact event with the plastic work and obtained the mean dynamic yield strength of copper and mild steel. Therein a direction to determine static yield stress was also proposed. In his further study [7], Hawkyard simplified the idea of Taylor's theory by replacing the momentum equilibrium by an energy equilibrium equation across the plastic front. The strain hardening effect was also considered. Therein it was shown that the derived model exhibits a closer proximity to experimental observations when compared with Taylor's theory. Jones et al. [8] modified Taylor's equation of motion for the un-deformed portion of the rod by introducing a term relating to mass loss from the un-deformed part to the plastic zone. In [9], a method was presented to accommodate a more complex material description including the strain-hardening effect. It was assumed that the material on the impact side of the rigid-plastic interface is brought instantaneously to rest for low-velocity impact. After logical investigation and careful insight into the plastic front propagation, Jones et al. [10] tried to avert one very uncertain measurement of final un-deformed section length in the Taylor test by replacing it with a position measurement corresponding to a particular value of strain. In their subsequent work [11], an elementary theory of the Taylor impact test was developed wherein the state of stress in the cylindrical rod at high strain rate was obtained by establishing a linear relationship between the length of the rigid part and the final deformed length.

Being one-dimensional in nature, the above-mentioned models are unable to accurately capture the three-dimensional deformation behaviour of the impacting rod, especially at the impact end. Moreover, while determining the material constants, these models rely on a final un-deformed section length, which is one of the most unreliable measurements in the Taylor test [11,12]. Though some of the uncertainty associated with the un-deformed section length was removed in [10], it is desirable to have a correlation between the material constants and the deformation parameters that can be measured accurately in the Taylor test. Two such measurements are the final deformed length of the rod and deformed diameter at the impact end. This motivated the present study.

In this paper, an axisymmetric model for Taylor impact is derived. The plastic zone at the impact end is approximated by the frustum of a cone. Longitudinal elastic wave propagation and the follow-up by a more severe plastic front are considered. The condition of incompressible plastic deformation is made use of in order to derive expression for longitudinal and transverse velocity in the plastic zone. The procedure to incorporate strain-hardening and thermal softening effect is discussed. Strain rates are computed over each time-step to obtain better accuracy, unlike the commonly taken [1,11] average value over the entire process. The accuracy of the proposed model in capturing the deformation behaviour of the Taylor rod is demonstrated by comparing the prediction with some experimental and numerical results reported in the literature.

Next the potential of the proposed model in estimating Johnson–Cook (J–C) material parameters [13] is explored. The distinctive effects of strain hardening, strain-rate and thermal softening in the J-C model were estimated by earlier researchers using a combination of pseudo-static, split-Hopkinson bar and torsion tests at ambient and elevated temperatures. But noting the unique strain and strain rate variations along the length of the Taylor specimen over time, Johnson & Holmquist [14] devised an empirical method under bounded strain and strain rate. Batra & Kim [15] tried to characterize parameters for metals plastic flow by recursive trials on reconstruction of an experimentally obtained stress–strain curve. In the present study, J-C parameters are determined through an inverse formulation of the proposed model based on box-constrained nonlinear optimization. The final deformed length and diameter at the impact end of the rod are used to construct the objective function. This is unlike the work by Ozel & Karpat [16], wherein the objective function is constructed through a stress function which is experimentally expensive to obtain. Estimated J-C parameters for 4340 steel via the proposed model are found to be in good accordance with their reference values.

The paper is organized as follows. The model is derived and its salient features are discussed in §2. In §3, the developed model is validated through some experimental and numerical test cases. The inverse analysis based on the proposed model is shown in §4. Finally, conclusions are drawn in §5.

2. Theory

When a rod-shaped specimen impacts a rigid stationary anvil an elastic compressive wave generates at the impact interface and travels through the rod in longitudinal direction with a speed equal to the sound speed in the medium. When this compressive wave reaches the free trailing end, it is reflected as a tensile wave causing deceleration of the rod. For a sufficiently severe shock, i.e. once the magnitude of the compressive wave reaches the flow stress, the impact end undergoes plastic deformation. The plastic front with maximum stress magnitude equal to the flow stress starts propagating from the impact interface at a much lower speed. The travelling plastic front interacts with the reflected precursor elastic wave at some intermediate point along the length of the rod. Elastic wave then get reflected at the elastic–plastic interface and travels towards the rare end of the rod. This back and forth movement of the elastic wave causes deceleration of the rod.

A typical deformed shape of the Taylor specimen obtained through experiments [17] is shown in figure 1a. The entire rod may be divided into two parts separated by an elastic–plastic interface: the mushroom zone near the impact end where plastic deformation takes place and the un-deformed zone at the rear end, which is characterized by back and forth movement of the elastic wave. In the present model, the mushroom zone is approximated by a frustum of a cone which is attached to a cylindrical part representing the un-deformed zone as shown in figure 1a. Density, Young's modulus, Poisson's ratio and dynamic yield strength of the rod are denoted as ρ, E, ν and Yd, respectively. The initial velocity of the rod is v0.

Owing to symmetry with respect to the longitudinal axis, the three-dimensional model may be reduced to an axisymmetric model with l0 and D0 as the initial length and diameter of the rod, respectively, as shown in figure 1b. Let at any instant of time, h be the axial shortening, l be the extent of the un-deformed part, s be the extent of the plastic part and D be the diameter at the impact end of the rod, all being time-dependent entities. So overall length L at any instant of time is given by L=l+s=l0−h. The initial (t=0) and final (t=tf) conditions of (L,l,h,s,D) are (l0,l0,0,0,D0) and (Lf,lf,hf,sf,Df), respectively.

The present development is based on some assumptions, required to preserve the conciseness of the model, as mentioned below.

— Deformation due to elastic wave propagation is small when compared with its plastic counterpart and therefore neglected. Consequently, the trailing end remains un-deformed and the mushrooming zone deforms plastically.

— No unloading takes place in the mushrooming zone and the impact end of the rod always remain in contact with the anvil.

— Carrington & Gayler [3] showed through changes in micro-structure in the stressed zone that stress wave travels with a non-planar front. However, within the present mathematical framework in order to obtain integrated volume of the temporally varying shape and spatial extent of the plastic zone, the plastic front is assumed to be moving in a plane surface to the longitudinal direction. This assumption implicitly implies that the integrated volume would be equivalent to the volume of the plastic zone with a non-planar front.

— When the reflected precursor elastic wave interacts with the plastic front, ideally it gets reflected partially from the elastic–plastic interface and transmitted partially to reach the impacting interface [18]. However, in the present formulation it is assumed that the elastic wave completely reflects back from the plastic front, while the plastic wave continues propagating towards the rear end.

— Impact velocity is less than speed of the plastic front and also in a range such that the rod does not experience any fracture or material separation during the deformation process (limiting our focus to capturing of dynamic plastic behaviour only).

In the subsequent subsections, a mathematical model which relates the material constitutive (such as, ρ,E,ν,Yd along with strain hardening and thermal softening parameters) with the geometrical parameters (l0, D0, h, l, s, D), especially the mushrooming at the impact end, is derived. First, a kinematically admissible velocity field is assumed and an incompressibility condition is invoked in order to arrive at the longitudinal and radial velocity distribution in the mushrooming zone (§2a,b). The length of the un-deformed part of the rod is estimated by considering axial shortening (h) and plastic front movement in §2c. Finally in §2d, an expression for overall deceleration is derived by computing the change in velocity in the un-deformed part of the rod after a round trip of stress wave. To incorporate the important notion of strain hardening, the current state of strain is obtained via the spatial gradients (averaging spatially over plastic zone but preserving the temporal variation) of the velocity components in §2e. A critical aspect of implementing the derived model is discussed in §2f.

(a) Velocity profiles

Consider the rod is axisymmetric with the axis in x-direction (figure 1b). Let vx and vr be the velocity in the longitudinal and radial directions, respectively. The longitudinal velocity is assumed to remain uniform spatially in the un-deformed part of the rod. It gradually decreases from the plastic front and becomes zero at the impact interface. This spatial distribution of the longitudinal velocity has been found in good accordance with the numerical simulations through ABAQUS [19] and simulations in [18] (figure 2).

Assumed longitudinal velocity distribution. (Exact functional form of vx is derived in equation (2.14) of §2a.)

Since there is no transverse flow occurring, the radial velocity is zero in the un-deformed part of the rod. In the plastic region, the radial velocity, which causes the transverse deformation (bulging), is symmetric with respect to the longitudinal axis (axisymmetric assumption). Therefore, the boundary conditions on velocity profiles may be summarized (with known l, s and velocity of the rod v at a particular time instant) as follows. Longitudinal velocity is v all over the trailing un-deformed part, i.e.
vx(x,r)=v,∀x∈[−l,0],r∈[−D02,D02]2.1
and is zero at the impact surface, i.e.
vx(x,r)=0,atx=s,r∈[−D2,D2].2.2
The radial velocity is zero all over the trailing un-deformed part, i.e.
vr(x,r)=0,∀x∈[−l,0],r∈[−D02,D02]2.3
and is also zero along the centroidal axis in the plastic zone, i.e.
vr(x,r)=0,∀x∈[0,s],r=0.2.4

Let Dx be the diameter of the rod at a distance x from the plastic front (figure 1b) towards the impact surface. To begin with, for the radial velocity at any point {x,r}∈{[0,s]×[−Dx/2,Dx/2]} a complete quadratic polynomial as given in equation (2.5) is assumed. The dimensionless entities x/s and r/Dx are chosen as the independent arguments of the polynomial. The unknown constants, αi may be determined by satisfying equations (2.1)–(2.4)
vrv=α0+α1xs+α2rDx+α3xsrDx+α4x2s2+α5r2Dx2.2.5
Substituting equations (2.3) and (2.4) into equation (2.5), one may have α0=α1=α2=α4=α5=0. Therefore, the expression for vr may be reduced to vr=α3(x/s)(r/Dx)v. As for the assumed profile, radial velocity vr, which causes the bulging at the impact end varies linearly with x, Dx may be written as
Dx=D0+D−D0sx=D0+βx(where,β=D−D0s).2.6
Finally,
vr=α3xsrD0+βxv(∵Dx=D0+βx).2.7

As mentioned earlier, the plastic deformation in the mushrooming zone is assumed to be incompressible. In axisymmetric formulation, the condition for incompressibility (zero divergence of the velocity field) may be approximated as
∂vr∂r+vrr+∂vx∂x=0.2.8
Using equation (2.7) in equation (2.8), one may get
∂vx∂x=−(∂vr∂r+vrr)=−2α3xsvD0+βx2.9
and
⇒vx=−2α3vs[D0+βxβ2−D0β2ln⁡(D0+βx)]+f(r),2.10
where, f(r) is the integration constant independent of x.

Satisfying equation (2.1) (i.e. vx=v at x=0), f(r) may be obtained as
f(r)=v+2α3vβ2sD0[1−ln⁡D0].2.11
Then, substituting f(r), equation (2.10) may be re-written as
vx=v−2α3vβ2s[βx−D0ln(D0+βxD0)].2.12
Next satisfying equation (2.2) (i.e. vx=0 @ x=s), α3 may be obtained as
α3=β2s2[(D−D0)−D0ln⁡(D/D0)].2.13
Finally, substituting α3 from equation (2.13) into equations (2.12) and (2.7), the longitudinal and radial velocity profile in the plastic zone (with, γ=D/D0) may be, respectively, obtained as
vx=v−vD0(γ−1−ln⁡γ)[βx−D0ln(D0+βxD0)]2.14
and
vr=β2vxr2D0(D0+βx)(γ−1−ln⁡γ).2.15
Here, a point to note is that the local profiles of longitudinal and radial velocity components (equations (2.14) and (2.15)) in the mushrooming zone are clearly functions of the velocity (v) of the un-deformed part. Though not explicit in the above equations, the velocity distributions will also depend on the other obvious influencing material constitutive parameters, because the current diameter, D (and, hence, β and γ explicitly appearing in the equations), is a direct outcome governed by those material properties as will be seen in the subsequent subsections.

(b) Mushrooming at impact end

The radial velocity at (x=s, r=0.5D) represents the rate of change of the radius at the impact end. Therefore,
d(radius)dt=β2vs(0.5D)2D0(D0+βs)(γ−1−ln⁡γ)=14βv(γ−1)(γ−1−ln⁡γ),2.16
where γ=D/D0, β=(D−D0)/s. Finally, the rate of change of the diameter at the impact end, currently known as the bulging rate, may be obtained as
dDdt=2d(radius)dt=12βv(γ−1)(γ−1−ln⁡γ).2.17

(c) Length of un-deformed part and axial shortening

The length of the un-deformed part of the rod reduces due to the movement of the plastic front towards the free end. Therefore, the rate of change of length of the un-deformed part may be related to the plastic front speed as
dldt=−Cp.2.18
Now, as shown in figure 2, the longitudinal velocity remains uniform (=v) in the un-deformed part and hence no axial shortening takes place there. In the plastic zone, the longitudinal velocity varies according to equation (2.14). Therefore, the rate of axial shortening which takes place only in the plastic zone may be denoted by the velocity difference at either end of the plastic zone as
dhdt=vx|x=0−vx|x=s=v.2.19
The length of the plastic zone increases with movement of the plastic front and decreases due to axial shortening in the plastic zone. Therefore, the rate of change of the length of the plastic zone may be written as
dsdt=Cp−ddt(h)=Cp−v.2.20
To have meaningful positive evolution of plastic zone s, the model is only considered to be valid when v<Cp.

(d) Deceleration model

Deceleration of the un-deformed part is caused by the repeated back and forth travelling of the longitudinal elastic wave with front speed Cl=(E/ρ). From the one-dimensional wave propagation theory, the relation between the magnitude of stress wave and the longitudinal velocity at any material point yields v=σCl/E. As the maximum stress magnitude in the plastic medium can at most reach the flow stress value, σ is replaced by Yd, safely ignoring the initial magnitude build-up in a very quick time-span. Then from [18], the change in particle velocity over one round trip of the shock-wave may be written as
Δv=(−ClYdE)−(ClYdE)=−2ClYdE.2.21
Now, let Δt1 be the time taken by the elastic compressive wave to travel from the plastic front to the free end. Therefore,
ClΔt1=l.2.22
The movement of the plastic front towards the free end results in reduction of the un-deformed zone, and therefore when returning the distance travelled by the reflected tensile wave is less. If Δt2 be the time taken by the reflected tensile wave to travel from the free end to the plastic front one may write
ClΔt2=l+∫0Δt1+Δt2(dldt)dt.2.23
Combining equations (2.22) and (2.23), one may have
ClΔt=l+[l+∫0Δt(dldt)dt],2.24
where Δt=Δt1+Δt2 is the total time for the round trip of the wave. Substituting equations (2.18), equation (2.24) may be recast as
ClΔt=2l−∫0ΔtCpdt=2l−CpΔt.2.25
Here, Cp is assumed to be invariant within the small time-frame of Δt. So, Δt may be finally expressed as
Δt=2lCl+Cp.2.26
Now from equations (2.21) and (2.26), with Δt→0 in the un-deformed part of the rod
dvdt≃limΔt→0ΔvΔt=limΔt→0−2ClYd/E2l/(Cl+Cp)=−Cl(Cl+Cp)Yd.lE2.27
It is shown in appendix A that this final form conserves the linear momentum.

(e) Thermo-mechanical behaviour

In many practical cases, the material may not behave as perfectly plastic and the flow stress may significantly depend on plastic strain and strain rate. Furthermore, in a highly dynamic event such as this Taylor impact case, the material may even undergo thermal melting. The need of incorporating thermo-mechanical strain-rate sensitivity and subsequent material hardening/softening in Taylor's impact problem was presented by Woodword et al. [20]. The procedure to incorporate such a constitutive law in the proposed model is discussed in this section. Towards this, two commonly used material models, viz. J-C model [13] and the dislocation-based Zerilli–Armstrong (Z-A) model for body-centred cubic (BCC) slip [21] as given, respectively, in equations (2.28) and (2.29), are considered.
Yd=(A+Bϵ¯pln)[1+Cln(ϵ¯˙plϵ˙0)][1−(T−TrTm−Tr)m]2.28
and
Yd=C0+C1exp⁡(−C3T+C4Tlog⁡ϵ¯˙pl)+C5(ϵ¯pl)n,2.29
where (ABnCm) and (C0C1C3C4C5n) are material parameters for the J-C and Z-A models, respectively. In the J-C model (equation (2.28)), ϵ˙0 is the reference strain at which A, B, n are characterized, Tr is the ambient temperature and Tm is the melting temperature. Now the temporal variation of plastic strain ϵ¯pl, plastic strain rate ϵ¯˙pl and temperature T are computed as follows.

Knowing the spatial distribution of the longitudinal and the radial velocity (equations (2.14) and (2.15), respectively) in the plastic zone, components of strain rate may be obtained as
ϵ˙xx=∂vx∂x=−β2vxD0(γ−1−lnγ)(D0+βx),2.30ϵ˙rr=∂vr∂r=vrr=ϵ˙θθ=−β2vx2D0(γ−1−lnγ)(D0+βx)2.31andϵ˙xr=12(vx∂r+∂vr∂x)=−β2vr2(γ−1−lnγ)(D0+βx)2.2.32
Ignoring the initial rapid elastic build-up, these strain rates actually denote the plastic components of the strain tensor effective within the plastic zone. Though the components of the strain rate are spatial function of (x,r), an average strain rate value computed at the centre of the plastic zone is used for calculating effective plastic strain. Thus substituting (x=0.5s, r=0) in equations (2.30)–(2.32), strain rate values at central point of the plastic zone are obtained as
ϵ˙xx=−βv(γ−1)D0(γ+1)(γ−1−lnγ)=−2D+D0dDdt,2.33ϵ˙rr=ϵ˙θθ=−βv(γ−1)2D0(γ+1)(γ−1−lnγ)=12ϵ˙xx2.34andϵ˙xr=0.2.35
Now the effective plastic strain rate may be computed as
ϵ¯˙pl=23ϵ˙αβϵ˙αβ=23(ϵ˙xx2+ϵ˙rr2+ϵ˙θθ2+2ϵ˙xr2)=∣ϵ˙xx∣=2D+D0dDdt.2.36
For the given spatial distribution of strains (equations (2.30)–(2.32)), the expression of effective strain rate at (x=0.5s, r=0) is found to be coincident with the volume averaged (over the whole mushrooming zone, the size varying with time) effective rate value. Then the effective plastic strain ϵ¯pl may be obtained by integrating this above rate over time. Finally, the computed effective plastic strain and strain rate are used in equations (2.28) and (2.29).

Next the remaining required input variable is current temperature T. Assuming an adiabatic condition, the temperature evolution from ambient value is computed from the rate of heat generation as
dTdt=χσ:ϵ˙plρCp=χYd.ϵ¯˙plρCp,2.37
where, Taylor-Quinney [22] empirical constant χ denotes the fraction of plastic work converted into heat (normally taken as 0.9) and ρ, Cp denote current density and specific heat density at constant pressure, respectively. The numerator indicates the current state of stress σ bounded by a closed convex set of admissible stresses Eσ=int(Eσ)⋃∂Eσ={(Yd)∈S∣f:=(3J2−Yd)≤0}, where J2 is the second invariant of deviatoric stress components. So at yield surface ∂Eσ, during the plastic flow, the plastic work is computed as (Yd.ϵ¯˙pl).

(f) Corrections for bulging rate

The changing diameter along with changing plastic front speed (Cp) as a result of strain hardening make it essential to integrate the governing ODEs (equation (a)–(f) in table 1), numerically with respect to time. Now the changing diameter influences β[=(D−D0)/s] and γ[=D/D0] directly. As shown in figure 1b, β is actually the measure of the angle created due to bulging from the original rectangular projected shape of the plastic zone. Therefore, β ideally increases from zero at the time of impact and reaches a final stable value βf[=(Df−D0)/sf] when impact is ceased. However, the presence of s in the denominator of β may result in instability in computation, especially during the initial phase of the impact, when s is very small but increasing very rapidly. In order to avoid such instability, an average value βav=0.5βf is reached in the following way. First equations (2.17)–(2.19), (2.27), (2.36) and (2.37) are integrated in time with assumed βav=0.5 [or βf=1] and the final un-deformed length lf and the final length of the plastic zone sf are obtained. Next incompressibility condition is imposed to estimate βav equating initial and final volume as
π4D02l0=π4D02lf+π12sf(D02+Df2+D0Df).2.38
Solving equation (2.38), the value of Df and that of βf are obtained as
Df=D02[{12sf(l0−lf)−3}−1]2.39
and
⇒βav=12Df−D0sf=D04sf[{12sf(l0−lf)−3}−3].2.40

Computed βav as per equation (2.40) after the first run is then used and another set of time-integrations is performed to obtain the actual response time-history and final deformed shape.

For the better understanding of the proposed model, the governing equations are summarized and implementation steps are given in table 1. In all steps, time integration is to be performed till the velocity of the rear end becomes zero and the entire rod comes to rest. For a given material, the speed of plastic front Cp may be computed from the square root of the density averaged slope of the stress–strain curve ((∂σ/∂ϵ)/ρ). However, as in the present case, the plastic front is compressive in nature and moves as a plane in the longitudinal direction; the simplified expression Cp=Yd/ρ as given in [23] is taken.

3. Validation of the proposed model

The model developed in the previous section is validated through some experimental and numerical results reported in the literature. A series of experiments on a cylindrical rod made of 4340 steel impacting a rigid stationary anvil with different impact velocities was conducted in [12]. Later hydro-code (MSPH) simulations of similar impact processes were performed in [24]. The proposed model is first validated with these test cases (table 2a). To start with, the J-C constitutive behaviour (equation (2.28)) is incorporated, as was used in MSPH [24]. Corresponding J-C parameters are given in table 2b.

Computed final length of the rod and diameter of the impact face are compared with their corresponding experimental values [12] in table 3. Numerically computed results through MSPH [24] and ABAQUS simulations are also shown in the same table. Next the Z-A constitutive behaviour (equation (2.29)) is also incorporated in the present model and corresponding deformations are shown in table 3. The Z-A material parameters as taken from [16] are given in 2c.

It can readily be seen that the predictions via the proposed model are in very good agreement with their experimental and numerical counterparts. It is observed that the present set of ODEs yields closer agreement with the experimental data when Z-A constitutive behaviour is incorporated. The marginal deviations in all the cases are in comparable scales with the small measurement error in experiments and theoretical assumptions made to derive the present mathematical formulation. This strengthens the claim of resolving the final deformed shape subjected to certain material parameters and initial conditions but with more significant ease than that of computationally exhaustive simulations.

Now, other than the final geometry, the actual time-evolution of deformation in the Taylor specimen with varying strain and strain rates must be validated before venturing confidently into predicting material constitutive behaviour [25]. So, to examine the time-resolved behaviour, the temporal variation of the desired parameters (diameter, length and velocity) are investigated. To this end, time histories for the diameter of the impact end, the deformed length and the spatially averaged velocity of un-deformed part for 37.97 mm long rod with initial striking velocity 181 m s−1 are computed with the proposed model (with J-C and Z-A material models) and compared with that of obtained via hydro-code simulation [24] and with that of FEM simulations through ABAQUS [19], in figures 3 and 4.

It may be seen that the predicted time histories are also in very good agreement with their numerical counterparts. The small discrepancies in the velocity time history (figure 4) are due the fact that the MSPH result, as digitized from Batra & Zhang [24], indicates the velocity time history of a point at the rear end of the rod—so it is stepwise, while the velocity computed from the proposed model is the average value over the un-deformed part of the rod and hence a continuous smooth deceleration is obtained. Whereas the velocity time history obtained from ABAQUS follows the MSPH trend, but with oscillations due to several reflections of stress-waves as the velocity profile is extracted on a single node at the trailing edge and on the axis of symmetry.

In the case of the Z-A model, the time-histories depict somewhat hardened behaviour as compared to simulations with the J-C model, as was also observed by Forde et al. [26]. But, investigating table 3, the Z-A model is seen to be closer to experimental observations, though with a slight softening; i.e. the Z-A model-based deformation is between J-C model-based deformation and experimental test data. In the absence of grain-size definition in the present macro-level model, the Hall–Petch (H-P) relation of deformation twinning, as suggested by Armstrong [27] to nullify this difference, could not be incorporated. However, for the test cases considered in this study, it is observed that the deviation of plastic deformation with respect to their experimental counterpart is in a negligible range. Thus, in the present context of a simplified approach with macro-level plastic behaviour, these micro-mechanics-based fine details yielding nonlinearity in a bulged profile were ignored and the conoid shape of a plastic zone with straight edges was assumed, as depicted earlier in figure 1b.

The nonlinearity in the above time-histories is due to the treatment of variable strain rate and corresponding strain hardening in place of assuming a constant strain rate as explained in §2e, generally in coherence with that of high-fidelity MSPH and ABAQUS simulations. These temporal variation of strains and strain rates are important as they do not remain the same throughout the impact event. When impact occurs, the impact end of the rod instantaneously comes to rest and experiences a steep velocity gradient which results in high strain rates at the impact end. Moreover, during the initial phase, since the plastic front has not yet moved sufficiently towards the rear end of the rod, the thickness of the plastic zone (s) is small, and therefore the plastic deformation is localized over a smaller region. However, as time progresses, owing to the increase in thickness of the plastic zone and deceleration of the rear end of the rod, the velocity gradient over the plastic zone decreases, which results in a decrease in strain rate. In order to illustrate whether the proposed model is able to capture this phenomenon, spatial distribution of longitudinal strain along the rod axis at different time instants is shown in figure 5. It is evident from the figure that the information regarding spatial and temporal variation of strain rates—unique attributes related to the Taylor test—can be extracted from the present model. Furthermore, their nonlinear distribution is found to be consistent with the physical phenomena experimentally monitored using high-speed photography and VISAR [26].

Now for the considered tests case (i.e. 37.97 mm×dia:7.595 mm rod impacting at 181 m s−1) the sensitivities of the parameters φ=(ABnCm) on Df and Lf are investigated and compared with that of an axisymmetric model simulated through ABAQUS [19]. As evident from figure 6, the sensitivity of different material parameters as obtained from the proposed model is found to be coherent with that of ABAQUS. In the case of parameter ‘m’, the sensitivities on Df and Lf are found to be much less both in ABAQUS and the present model and are thus not shown here.

For the further validation of the proposed model, next the predictions are compared with some other test results given in [28]. The rod is composed of aluminium alloy 7A04-T6 (E=69.35 GPa, ν=0.31 and ρ=2850 kg m−3) with D0=12.62 mm. The material constitutive law to capture the plastic behaviour is taken to be a modified J-C constitutive law with Voce empirical formula [29] as
Yd=fpl[1+Cln(ϵ˙plϵ˙0)][1−(T−TrTm−Tr)m],3.1
where fpl is defined with A1, A2, t1, σu, ϵu and w as the material parameters, given in table 4, as follows:
fpl={A1+A2[1−exp−ϵpl/t1],ifϵpl<(ϵu−A1E)σu[w(1+ϵpl−ϵu+A1E)+(1−w)(ϵplϵu)ϵu],ifϵpl≥(ϵu−A1E).3.2

Comparisons of the computed final length of the rod and the diameter of the impact face with their corresponding experimental values are given in table 5. For the present focus of capturing the plastic response only in a simplistic form, the velocities considered are such that they will not initiate material damage or fracture. A very good agreement between the experimental values and the prediction via the proposed model is again observed.

4. Parameter estimation via an inverse analysis

If Taylor’s original idea [1] is to be followed, the average dynamic yield limit can be estimated by feeding the deformed length parameters (lf and Lf) in the following Taylor formula as
Yd∗=ρv022l0−lfl0−Lf1ln⁡(l0/lf).4.1
In addition to the problem in precise determination of lf the major limitation of the mathematical model that leads to equation (4.1) is that it fails to capture the thermo-mechanical behaviour of the material. Therefore, equation (4.1) is restricted only to material which exhibits rigid-perfect plasticity. Since most of the materials do not behave as perfectly plastic especially at high strain rate the Taylor test in its classical form is being used only for validation rather than a material test. However, the beauty of Taylor’s work perhaps lies in the fact that it showed how an estimate of dynamic plasticity can be made through a simple test, and therefore it still bears significance in the context of metal plasticity. Keeping this in mind, an improved mathematical model for the Taylor test was derived in the previous sections. It was shown that the proposed model not only predicts the final deformed shape accurately, it also captures their time evolution reasonably well. Now the main purpose of this section is to further explore the potential of the proposed model in order to have better understanding of material behaviour especially at high strain rate.

To start with in order to draw the similitude between the proposed model and the Taylor approach for perfectly plastic material, the differential equations ((a)–(h) in table 1) are solved with constant yield value. By parametric variation of the constant yield limit ((500 600 700 800 900 1000 1100) MPa), the deformed geometries (Lf and lf) are obtained from the present model and fed into equation (4.1). The obtained estimate ((485 573 676 780 885 990 1097) MPa) are indeed very close to the assumed values. This shows that Taylor’s formula is a special case of the present model if the yield value is assumed to be constant.

Now, as stated above, real materials do not behave as perfectly plastic and thus it is important to understand the actual evolution of yield strength as a function of different thermo-mechanical variables. Owing to nonlinearity in the equation, closed form expressions for different material parameters are not possible. Therefore in this study, material parameters estimation is performed through an inverse technique based on the proposed model. For demonstration purposes, estimations of J-C constitutive parameters (ABnCm) are considered here. However, the strategy may be extended to other material models such as the Z-A model or the MTS model.

The inverse problem is posed under an optimization framework with a standard least-squares based objective functional as follows:
Λ(φ)=12∑i[(Df,ic(φ)−Df,im)2(Df,im)2+(Lf,ic(φ)−Lf,im)2(Lf,im)2],4.2
where, D()f,i and L()f,i are, respectively, the final diameter and overall length (=l0−hf) of the rod for ith test situation. Superscripts ‘m’ and ‘c’ denote the measured data from experiments and computed data from the governing forward equations, respectively (table 1). Each test case here represents a particular set of initial diameter, length and impact velocity. φ=(ABnCm)T denotes the decision variable in the minimization problem. It is to be noted that the objective functional is constructed using final total length and diameter at the impact end and thus the uncertainties associated with measurement of the un-deformed part is avoided.

The inverse problem is cast as a bound constrained optimization problem as follows:
minφ{Λ(φ)|φL≤φ≤φU}(subject to equation (a)–(h) in table 1),4.3
where, φU and φL denotes the upper and lower bounds, respectively, of the decision variable. The optimization is performed via a standard gradient based quasi-Newton method [30]. The gradient of the cost functional (equation (4.2)) is computed via first-order optimality criteria as
g=∂Λ∂φ=∑i[(∂Df,ic∂φ)(Df,ic−Df,im)(Df,im)2+(∂Lf,ic∂φ)(Lf,ic−Lf,im)(Lf,im)2],4.4
where (∂Dcf,i/∂φ) and (∂Lcf,i/∂φ) represents parameter sensitivities at the final time for the ith test case. These sensitivities are shown in figure 6 and some patterns for each of the parameters were obtained. Following Taylor’s idea [1], these sensitivity patterns may be used to inversely calculate the dynamic yield limit. But a more rational approach would be to incorporate accurate mathematical sensitivity under an optimization framework. Thus to arrive at the final value of these sensitivities at the end of the time-evolution of D and L(=l0−h), a sensitivity based differential equation is formed and solved recursively. Considering the time discretized form of the proposed forward model with yk={lkvkhkDkϵ¯plkTk}T, the discrete form of the sensitivity differential equation is as follows:
∂yk+1∂φ=∂yk∂φ+Δtg(yk,φ,tk)4.5
where
g:=∂y˙k∂φ={−∂Cpk∂φ(1+CpkCl)∂Cpk∂lk[Cpklk∂lk∂φ−(2−CpkCl)∂Cpk∂φ]∂vk∂φβ2[∂Dk∂φ{vkD01−(D0/Dk)−ln⁡(Dk/D0)(ψk)2}+∂vk∂φ(D0kψk)]β(D¯0k)ψk[vk∂Dk∂φ{2D0D¯0k−(D0k)2Dkψk}+(D0k)∂vk∂φ]χρCp[ϵ¯˙plk∂Ydk∂φ+Ydkβ(D¯0k)ψk(vk∂Dk∂φ{2D0D¯0k−(D0k)2Dkψk}+(D0k)∂vk∂φ)]}4.6
In equation (4.6), the superscript ‘k’ denotes the kth time-step, ∂Cpk/∂φ=(1/2ρYdk)(∂Ydk/∂φ) and
∂Ydk∂φ={[1+Cln⁡ϵ˙plk][1−(Tk∗)m](ϵplk)n[1+Cln⁡ϵ˙plk][1−(Tk∗)m]B(ϵplk)nln⁡ϵplk[1+Cln⁡ϵ˙plk][1−(Tk∗)m][A+B(ϵplk)n][1−(Tk∗)m]ln⁡ϵ˙plk−[A+B(ϵplk)n][1+Cln⁡ϵ˙plk](Tk∗)mln⁡(Tk∗)}T4.7
with D0k=Dk−D0, D¯0k=Dk+D0, ψk=Dk−D0−D0ln⁡(Dk/D0) and Tk∗=(Tk−Tr)/(Tm−Tr). Starting with a zero initial value (i.e. ∂y0/∂φ=0), equation (4.5) is solved along with the governing equations (a)–(g) of table 1 and the final value of ∂y/∂φ is used to estimate ∂Dcf/∂φ and ∂Lcf/∂φ. After obtaining sensitivities, material parameter φ is updated and the procedure is repeated till the cost function is sufficiently reduced. The iterations are performed via the function ‘fmincon’ available in commercial software package Matlab ([31]), wherein the cost function (equation (4.2)) and gradient values (equation (4.4)) are supplied for obtaining material parameter update.

Next the potential of the model along with the inverse formulation is illustrated through J-C parameter estimation of 4340 steel. Density, shear and bulk modulus are assumed to be known apriori as given in table 2. Length and diameter of the un-deformed projectile are taken from l0∈(20,40) mm and D0∈(10,15) mm, respectively, with varying L/D ratio between 1.33 and 4. Two velocity ranges 75–115 m s−1 and 230–270 m s−1 are taken for denoting low strain rate and high strain rate, respectively, generating a total of 180 test cases altogether. Now in the absence of real-life experimental observations for all these test data, two methods are employed to generate Dmfs and Lmfs.

First, the proposed forward model itself (table 1) is used to generate the experimental data for a known material, 4340 steel. Now there are three measurements to be made in each of the test cases—final length, deformed diameter at the impact end and initial impact velocity. Noting the accuracy feasible in measuring geometric dimensions of the retrieved samples through Laser transducers [32], the tolerance of velocity measurement does pose scope for pollution in the experimental data. In order to handle such possibilities, a perturbation (i.e. measurement noise) is added to the initial velocities for all test cases while simulating the inverse problem. Thus, data pollution comes in terms of initial velocities which finally affect the computed geometry after impact. The measurement noise is added as follows:
vnoisy=vexact[1+α∗(r)],4.8
where α is noise intensity and r denotes a random variable uniformly distributed in [−1, 1].

Secondly, an axisymmetric finite-element model (FEM) is implemented in the ABAQUS [19] software package with J-C plasticity for 4340 steel to generate synthetic data. Computed Dmf and lmf through FE analysis are used as measurements in the inverse problem. It may be noted here that the mathematical model used to generate measurements is essentially a numerical approximation (via FEM) of the original problem and the proposed model used for material parameter estimation is entirely different. This is in sharp contrast to the previous experiment, wherein the same model is used for generating measurements and parameter estimation.

Thus, for the present inverse procedure, the model dependencies in estimating plasticity parameters are completely avoided. In both cases, bounds on the decision variable are taken as φL=[0.5×109 0.3×109 0.2 0.1×10−1 0.1×101]T and φU=[1.0×109 0.7×109 0.3 0.2×10−1 0.2×101]T. The mean of the upper and lower bounds is taken as the initial guess of the material parameters to start the inverse problem. However, in test cases considered in the present analysis, no significant sensitivity of the initial guess to the reconstructed parameters is observed.

Table 6 shows reconstructed parameters for two cases. For the first set of measurements 0%, 1% and 5% (Cases I, II and III) noise intensities are considered. It is evident from the reconstruction result that the inverse procedure can appropriately identify plasticity parameters with different orders of magnitude.

It is also noted that the choice of proper bounds is very important for identifying parameters of different magnitude which can be done apriori from a crude knowledge about the material model to be identified. As expected, with the increase in data pollution the solution also deviates from the true value. The order of deviation depends on the sensitivity of the decision variables on the constructed objective function and thus affects the reconstructed solution. In the case of measurement obtained from ABAQUS, the reconstruction results deviate from the reference value especially for parameter A. As shown previously in figure 3, the final length and diameter obtained through ABAQUS do deviate a bit from the proposed model. This results deviation (though not out of range) of the estimated values from the reference values when experimental values are simulated through ABAQUS. However, as noted from table 6, inverse formulation with the proposed axisymmetric model can reasonably identify the J-C parameters of 4340 steel. The parameter estimation may further be improved with an efficient inverse algorithm [33,34].

Finally, the estimated dynamic yield strengths of 4340 steel (rod of l0=37.97 mm, D0=7.595 mm impacting with a velocity of 181 m s−1) for different measurement data are compared with their reference behaviour in figure 7. Obtained stress–strain behaviour as apparent from figure 7 is adequately indicative of the robustness of the proposed algorithm in estimating J-C plasticity parameters. As expected, measurements generated from ABAQUS dynamic yield strengths do deviate from reference value. However, the characteristics remain the same with the increase in effective plastic strain.

5. Closure

An axisymmetric model for the impact of a flat-ended cylindrical rod on a rigid stationary surface is developed in this study. Plastic deformation causing bulging at the impact end, propagation of plastic front along the length of the rod, spatial and temporal variation of plastic strains and strain rates constitute the basis of the model. The proposed model gives an implicit relationship between material yield strength and deformation of the impacting rod in the form of a set of ODEs in time. The derived model is based on some assumptions which are directly or indirectly related to the impact velocity. Higher impact velocity may cause hydrodynamic effects (if, v>Cp) or shock propagation (if, v>K/ρ, K being bulk modulus) or even fracture, which violate the premise of the present model. On the other hand, if the impact velocity is very small such that elastic deformation and unloading play a dominant role, the prediction may be affected since these effects are neglected in the present analysis. It is important to note that the efficacy of the present model has to be assessed in the context of Taylor’s impact test, where impact velocities are always taken within a range which is different for different materials. In that range the present model gives a reasonably good prediction. This is demonstrated by comparing the obtained estimate with some experimental and numerical results reported in the literature. It is shown that the proposed model is efficient in not only predicting the final deformed shape of the rod but also in capturing the response development over time.

Next the proposed model is taken to identify J-C parameters via a constrained nonlinear optimization scheme. The relationship between material strength and bulged diameter at the impact end is explored in the inverse formulation. This is in contrast with the one-dimensional models which rely on the length of an un-deformed section whose experimental measurement in the Taylor test is uncertain. Moreover, the proposed model takes into account the temporal variation of strain, strain rates and temperature effect unlike the existing model [1,11], wherein a time-average strain rate value is used. The major advantage of such temporal variation is that, in place of only mean dynamic yield strength, the actual evolution of dynamic yield strength as a function of strain and strain rate may be obtained. The accuracy of the model with the inverse framework in estimating J-C parameters to capture that evolution is demonstrated. It is also shown that the model is capable of estimating true material constants with reasonable accuracy when a different forward model (here FEM-based ABAQUS simulation) is used to generate the measurement data. While the present inverse formulation effectively illustrates the potential of the proposed model, a more detailed exploration with different material models (including micro-mechanical effects) is yet to be addressed. The present inverse formulation relies on the deformation parameters obtained from the final recovered rod alone. As the proposed model can also predict deformation time history with reasonable accuracy, it is anticipated that an inverse formulation enriched with intermediate profiles obtained from high-speed photography [26] would provide more information about the material behaviour. This may constitute the content of a future work.

Appendix A. Momentum conservation

Taking the trailing rigid part (length, l) of the rod, with cross-sectional area (Ar=(π/4)D02) under consideration, the mass M(=ρArl) moves at v and subjected to yielding pressure [σ(=ρClv)=Yd] at the elastic–plastic interface. So the deceleration equation (2.27) is rearranged as
dvdt+YdCpρCll=−Ydρl[∵Cl=Eρ]A 1
Substituting, Yd=−ρClv in LHS
ldvdt−vCp=−YdρA 2⇒ldvdt+vdldt=−Ydρ[∵dldt=−Cp]A 3⇒ρArldvdt+vddt(ρArl)=−YdArA 4and⇒ddt(Mv)=−YdAr,A 5
which, invariably, is the statement of linear momentum conservation.

. 2005A history of the Taylor test and its present use in the study of lightweight materials. In Design and use of light-weight materials (eds Teixeira-DiasF, DoddB, LachE, SchultzP), pp. 12–24. Aveiro, Portugal: University of Aveiro.

. 1983A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures. In Seventh Int. Symp. on Ballistics, The Hague, The Netherlands 19–21 April, pp. 541–547. International Ballistics Society.