Material characterization at high strain rate by Hopkinson bar tests and finite element optimization

Materials Science and Engineering A 487 (2008) 289–300

Material characterization at high strain rate by Hopkinson bar tests and finite element optimization M. Sasso a,∗ , G. Newaz b , D. Amodio a a

b

Mechanics Department, Polytechnic University of Marche, via Brecce Bianche, 60131 Ancona, Italy Mechanical Engineering Department, Wayne State University, Engineering Bldg., 48202 Detroit, MI, USA Received 16 May 2007; received in revised form 4 October 2007; accepted 10 October 2007

Abstract In the present work, dynamic tests have been performed on AISI 1018 CR steel specimens by means of a split Hopkinson pressure bar (SHPB). The standard SHPB arrangement has been modified in order to allow running tensile tests avoiding spurious and misleading effects due to wave dispersion, specimen inertia and mechanical impedance mismatch in the clamping region. However, engineering stress–strain curves obtained from experimental tests are far from representing true material properties because of several phenomena that must be taken into account: the strain rate is not constant during the test, the specimen undergoes remarkable necking, so stress and strain distributions are largely non-uniform, and the temperature increases because of plastic work. Experimental data have been post-processed using a finite element-based optimization procedure where the specimen dynamic deformation is reproduced. Optimal sets of material constants for different constitutive models (Johnson–Cook, Zerilli–Armstrong and others) have been computed by fitting, in a least mean square sense, the numerical and experimental load–displacement curves. © 2007 Elsevier B.V. All rights reserved. Keywords: Hopkinson bar; High strain rate; Thermoplastic instability; Constitutive modeling; Finite element optimization

1. Introduction The interest in the mechanical behavior of materials at high strain rates has increased in recent years, and by now it is well known that mechanical properties can be strongly influenced by the speed of load application. The split Hopkinson pressure bar (SHPB) has become a widely used experimental method to determine, for several classes of materials, the dynamic stress–strain curves at strain rates that range from 102 to 104 s−1 . Historically, the first who used long thin bars for measuring impact generated pressure waves, was Hopkinson [1]; later, this method was further developed by Kolsky [2]. More recently, the Hopkinson bar technique, initially focused on compression tests, has been modified and successfully applied to tensile [3,4] and torsion tests [5]. Moreover, in addition to metallic materials that are now routinely studied, its use has been extended to many other kind of materials, such as polymers, ceramics, composites and foams. In order to study this wide variety of materials, including brittle materials [6], a lot of studies have been conducted with the aim of

∗

Corresponding author. Tel.: +39 071 220 4520; fax: +39 071 220 4801. E-mail address: [email protected] (M. Sasso).

0921-5093/$ – see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.msea.2007.10.042

improving the experimental methodology and the response quality in the range of small deformations that typically is the weak spot of the SHPB technique. Anyway, also dealing with the large plastic deformations, distinctive of ductile metallic materials, a lot of work is still needed to completely understand and modeling the phenomena that are experimentally observed. In fact, as it will be shown in the next sections, the engineering stress–strain curves obtained from experimental tests are far from representing true material properties and, up to now, there is not any unified and universally accepted methodology to compute the actual material behavior from experimental data. In this work, an experimental SHPB device has been used to perform tension tests on steel specimens. The collected data has been post-processed using an optimization procedure based on a finite element model in order to estimate the material constants for some constitutive models that are generally adopted in impact applications. 2. Theoretical background The compression Kolsky bar, or split Hopkinson pressure bar, consists of two long metal bars. These bars sandwich a short cylindrical specimen, as schematically reproduced in Fig. 1.

290

M. Sasso et al. / Materials Science and Engineering A 487 (2008) 289–300

Fig. 1. Scheme of the SHPB device.

The free end of the first (input) bar is impacted by a projectile that usually consists in a third bar (striker bar); the compressive pulse generated by the impact propagates down the input bar into the specimen. Several reverberations take place within the specimen, causing the rise of its stress; in the meanwhile, a compressive pulse is transmitted to the second bar (output bar) and a tensile pulse is reflected back to the input bar. Conventionally, these pulses are measured by strain gages placed on the input and output bars. The bars must be designed to remain within their elastic limit throughout the test while the specimen undergoes large plastic deformation. Let the incident, reflected and transmitted pulses be denoted respectively by εI (t), εR (t) and εT (t). All of them are function of time and can be regarded also as displacement, strain and stress perturbations travelling along the bars with the sound speed of the material C0 and determining a particle velocity equal to C0 ε. By applying the elementary one-dimensional elastic wave theory and assuming that the specimen is in dynamic equilibrium, the average nominal strain δ, strain rate δ˙ and stress S along the gauge length of the specimen can be determined by applying the well known equations: ˙ = − 2C0 εR (t), δ(t) LS  2C0 t δ(t) = − εR (t) dt, LS 0 S(t) =

AB EB εT (t), AS

(1) (2) (3)

where AB and EB are the bar cross-sectional area and Young modulus, and where LS and AS represent the specimen original gauge length and cross-sectional area, respectively. Thus the strain rate and the strain in the specimen can be determined from the reflected pulse only, and the stress can be determined from the transmitted pulse only. Eliminating time t in Eqs. (1)–(3) by synchronizing the acquired signals, it is possible to detect the stress–strain curve “followed” by the specimen material. It must be noted that unfortunately the strain rate is not constant during the test and also that the computed stress and strain values are the engineering (or nominal) ones, because they are evaluated considering the initial specimen length and cross-sectional area. This treatment is based on some basic assumptions: • The waves propagating in the bars can be described by onedimensional elastic waves propagation theory. • The specimen inertia effect is negligible and the specimen is in dynamic equilibrium. • The stress and strain fields in the specimen are uniform. • The friction effect in the compression test is negligible.

These assumptions have been extensively studied in past decades by several authors [7–13], with the result that numerous experimental and numerical methods have been successfully established to satisfy these hypotheses or to overcome their imperfect realization in practical application. Theory and analysis of the tensile SHPB are basically identical to the compressive one [14], except for the change in shape and sign of the pulses. The principal drawback of tensile tests lies in the non-uniaxial stress and non-uniform strain distributions; in fact, the short dimension of the specimen accentuates the necking instability much more than in standard tensile tests performed on comparatively much longer specimens. On the other hand, tensile tests permit to reach a higher effective strain rate than in the compression tests carried out with the same equipment and under the same conditions of impacting mass and velocity. 3. Experimental device To perform tensile tests, the indirect tension method has been used [15]; in Fig. 2, a picture of the experimental apparatus is given together with two details of the clamping region. The threaded-ends of the impact tension specimen can be screwed directly into the Hopkinson bars or by means of two threaded collars. The first solution has been preferred because it reduces the unavoidable spurious reflections (due to mechanical impedance variation) of the first wave when it transits across the collars. Twisting the bars by hand, a small pre-tension load is applied on the screwed joints and full contact is achieved between the bar ends. In such a condition, the compressive strain pulse, initially generated in the input bar by the impact with the striker bar, travels through the contact surfaces from the input to the output bar with no plastic deformation of the specimen [16]. As the wave reaches the free end of the output bar, it is reflected back as a tensile wave; from now on, the same conditions of the compressive test are reproduced except for the sign of the waves: when the tensile pulse reaches the specimen, it is partially transmitted to the input bar and partially reflected as a compressive wave to the output bar. The input and output bars, made in AISI 4140 steel, can slide into special supports with low friction; the striker bar is fired by an air gun at speed of approximately 10 m/s generating intense stress waves up to 200 MPa. The principal characteristics of the equipment are summarized in Table 1. The bars have been chosen to have a small diameter and a high L/D ratio in order to reduce the dispersion of the waves due to the so-called Pochhammer–Chree effect. With the same aim, some discs of soft material (e.g. sheet paper) have been interposed at the striker bar–incident bar interface, acting as “pulse shaper” to suppress high frequency waves. Strain gauge rosettes, placed at

M. Sasso et al. / Materials Science and Engineering A 487 (2008) 289–300

291

Fig. 2. Overview of the experimental equipment. Table 1 Experimental rig data Length (mm) Striker bar Input bar Output bar

1219 3073 3073

Specimens

2, 4 and 8

Diameter (mm) 12.7 4

Material AISI 4140 steel, E = 207 GPa, density = 7800 kg/m3 , sound speed = 5150 m/s AISI 1018 CR steel

half-length of the bars, convert the strain into voltage signals that are acquired by a 12-bit digital oscilloscope; a couple of optical gates allows to trigger the signal acquisition and to measure the striker bar speed at the same time, that is very useful to verify the strain gages calibration. The adoption of different specimen lengths permits to perform tests at different strain rate with the same experimental set-up; the small specimen diameter makes the lateral and axial inertial stresses negligible (about 1 MPa only, computed by Davies and Hunter formula [11]) even if the optimal slenderness ratio is not fully respected.

that better approximate the average of the corresponding series; instead, the dashed curves represent the band of ±1S.D. from the average curves. It must be said that the striker bar speed was actually slightly different among each test, varying from 9.7 to 10.3 m/s (this is mainly due to the non-perfectly constant friction force between the internal surface of the gun and the striker bar), so the above reported deviations actually refer to tests carried out with slightly different conditions and strain rates. However, the general rule providing for flow stress increase at higher strain rates has been seen to remain valid also within the small variations occurred in each series, so that to slightly higher

4. Experimental results With the apparatus just described, dynamic tension tests have been performed on AISI 1018 CR steel specimens; the experimental results, in terms of engineering stress–strain curves computed by formulas (2) and (3), are shown in Fig. 3a, while the engineering strain rates are reported in Fig. 3b. As outlined, different average strain rates have been achieved by using different specimen lengths, precisely 2, 4 and 8 mm with 4 mm constant diameter: accordingly to (6), the shortest specimen (referred D4L2) showed the highest average strain rate of about 1500 s−1 , the median length specimen (D4L4) had 1000 s−1 , and the longest (D4L8) gave the lower strain rate of about 500 s−1 . Several tests for each specimen geometry have been carried out revealing an overall good repeatability of the tests themselves: considering five replications of each test, the stress standard deviations have been the 2.4, 2.0 and 5.1% of the average values for the D4L2, D4L4 and D4L8, respectively, while the strain rate standard deviations were 4.6, 3.5 and 3.1%. The continuous curves in Fig. 3 represent the experimental curves

Fig. 3. (a) Tensile tests on AISI 1018 CR steel: stress–strain curves. (b) Tensile tests on AISI 1018 CR steel: strain rate–strain curves.

292

M. Sasso et al. / Materials Science and Engineering A 487 (2008) 289–300

lies in the specimen geometry that is characterized by an abrupt change in diameter passing from the clamping to the gauge zone, and by a very small length to diameter ratio; this determines a high stress and strain non-uniformity along the entire specimen gauge length, which is so short that there is not enough space to recover a certain homogeneity. When a diffused necking occurs, as shown in Fig. 5, the stress and strain distributions become suddenly triaxial, and the stress–strain curves are characterized by an absolute maximum that is reached at the very beginning of the test at quite low strain level (Fig. 3a), after which the stress has a large reduction (geometrical softening). Under these conditions, the engineering measured data (S − δ) cannot be easily converted into true stress–logarithmic strain (σ–ε) relationship by means of the usual formulas: Fig. 4. Flow stress at 0.04 of strain vs. strain rate (engineering values).

strain rate corresponded a slightly higher stress–strain curve; this means that the standard deviations here computed represent a sort of upper limit, or worst case, of the system repeatability, and if one could perform a perfect control on the impact speed, he would probably measure even smaller deviations. Stress–strain curves are seen to move up as the strain rate increases, making the strain rate effect quite evident for the steel tested. Fig. 4 shows the flow stress at 0.04 strain at different strain rates, also including data coming from two quasi-static tests at 10−3 and 100 s−1 . From these dynamic results, a formidable strain rate sensitivity of the material (for the reason that strain rate is only varying from 500 to 1500 s−1 ) may be deduced, but it is worth to anticipate that the true strain rate levels experienced by the specimens, particularly in the necking region, are actually much higher than the approximate time averaged engineering strain rate reported above. Moreover, the increase of flow stress is not only due to the strain rate effect, but also to the triaxiality of the stress distribution: the nominal stress is proportional to the cross-section average of axial stress component but, being the radial and hoop components not negligible, especially for the shorter specimens, the effective stress increment can be assumed to be less than that observed in nominal curves of Fig. 3a. This fact will be further discussed in Section 7. The main objective of this kind of tests is to describe the material behavior by means of constitutive models, whose constants are determined by fitting the experimental data. The difficulties of applying this approach moving from dynamic tensile tests,

A0 , A   A0 ε = log . A

σ=S

(4a) (4b)

In fact, while (4b) can be considered applicable even if the strain distribution is not constant across the section, Eq. (4a) provides only the average axial true stress; to compute the average effective stress the Bridgman formulation should be applied: σ = σaxial

    2R a  −1 1+ . log 1 + a 2R

(5)

However, its theoretical hypothesis are not fully satisfied and Eqs. (4b) and (5) require the knowledge of the instantaneous cross-section radius a and profile curvature radius R at the neck section that are very hard to be measured in such dynamic experiments. Furthermore, two other aspects must be taken into consideration in order to assess the actual material properties: • The temperature rise is not negligible in this kind of tests, because the rapidity of the deformation makes the thermal diffusion distance comparable with the specimen size; so high strain rate processes can be considered adiabatic, and the large plastic deformation work is transformed into heat with almost instantaneous temperature increment [17]. • The strain rate is not constant during the tests, and not even uniformly distributed into the specimen. Because of these difficulties, fitting of experimental data has been carried out by a finite element optimization procedure,

Fig. 5. Two tensile specimens after the test.

M. Sasso et al. / Materials Science and Engineering A 487 (2008) 289–300

that will be described in Section 6. In the next section, a brief description of the constitutive models here adopted is presented. 5. Constitutive models for the material A variety of constitutive models exists for the representation of engineering materials behavior when the strain rate and the temperature effects cannot be neglected; their objective is to fit the experimental data such as shown in Fig. 4 with only one equation. A review of the most used models is given in [18]; in this paper some of them have been implemented. At low and constant strain rate, and at room temperature, metals are known to work harden along a path usually described with the so-called “parabolic hardening” equations: σ = (σ0 + kεn ),   ε n σ = σ0 1 + , εy

(6a) (6b)

where σ 0 is the yield stress, k and n are the work hardening coefficient and exponential factors, respectively; often this relationship is further simplified into the Hollomon equation: n

σ = kε .

(7)

For mild steel it has been found that the effects of temperature and strain rate are basically separated; Johnson and Cook [19] proposed the phenomenological constitutive relation:   ε˙ n (1 − T ∗m ), (8) σ = (A + Bε ) 1 + C log ε˙ 0 which has been frequently used in impact and dynamic analysis due to its simplicity. In Eq. (8), A, B and n are the three coefficients describing the quasi-static behavior of the material accordingly to (6a), while C and m account for the strain rate hardening and thermal softening effects; ε˙ 0 is a user-defined reference strain rate, generally taken as unity, and T* is the homologous temperature defined by T∗ =

T − Troom , Tmelt − Troom

(9)

where T is the absolute temperature, and Troom and Tmelt represent the room and melting absolute temperatures, respectively. The main limitations with this model lie in the fact the strain and strain rate hardening and the temperature softening phenomena are uncoupled, and may be this is the principal reason that makes the Johnson–Cook model to be considered to work well at moderate strain rate level, not much over 102 s−1 ; anyway, it is considered a “milestone” of constitutive modeling. Many authors proposed additional empirical equation, among them Klopp–Clifton [20] and Litonski [21], whose models are respectively expressed by  n  m  −τ εp ε˙ T σ = σ0 , (10) ε˙ 0 Troom εy   ε n σ = σ0 1 + (1 − aT )(1 + b˙ε)m , (11) εy

293

taking into account all the phenomena previously mentioned in a manner slightly different from the Johnson–Cook relationship. While the Johnson–Cook constitutive relation is purely empirical, Zerilli and Armstrong [22] applied dislocation dynamic concepts to develop a constitutive model able to account for strain hardening, strain rate and temperature effects in a coupled way; they expressed the flow stress of a material as σ = σG + σ ∗ (T, ε˙ ),

(12)

where σ G represents the athermal barriers, determined by the structure of the material, encountered by dislocations on their course, and σ* is due to the thermally activated barriers, that is the barriers that can be overcome by thermal energy. For bcc metals, the σ* is essentially determined by the force opposing the dislocation movement at atomic level and the consequent overall lattice potential; so the thermally activated motion is not dependent on the strain, and the strain hardening becomes independent of strain rate and temperature, leading to the well known formula: σ = C0 + C1 exp(−C3 T + C4 T log ε˙ ) + C5 εn (for bcc metals).

(13)

In literature, this model has been found to work well up to very high strain rate levels, even those reached in the Taylor cylinder impact test. These constitutive models have been implemented into the finite element optimization procedure that will be described in the next section. All of them account for the thermal softening so they are able to describe, even if in a macroscopic way, the strain localisation in shear bands that could lead to the so-called “thermoplastic instability” if the softening rate turns out to be higher than the hardening rate. For comparison, besides the aforementioned material models, finite element simulations have been performed also using the quasi-static material power law curve (6a) and the Perzyna model, expressed by   m  ε˙ σ = σ0 (ε) 1 + , (14) C that accounts for the strain rate, but not for the temperature effect. In Eq. (14) C and m are the coefficients governing the strain rate effect, and σ 0 (ε) represents the quasi-static material properties, modeled by (6a). 6. Finite element optimization Within the optimization module of the Ansys® finite element commercial code [23], the experimental tests described in previous sections have been numerically reproduced. Preliminary finite element simulations showed that the lateral extremities of the specimen remain within the elastic limit and no relevant error is committed if they are considered infinitely rigid. Accordingly to this simplifying assumption and taking advantage of the symmetrical geometry, only the “useful” part of a half-length specimen has been modeled in order to reduce the computation time.

294

M. Sasso et al. / Materials Science and Engineering A 487 (2008) 289–300 Table 2 Constrained and unconstrained optimization parameters

Johnson–Cook Zerilli–Armstrong Klopp–Clifton Litonski Perzyna

Total parameters

Constrained parameters

Design variables

5 6 5 6 5

3 3 3 3 3

2 3 2 3 2

experimental one:

 2 ε [Snum (δ) − Sexp (δ)] , S.D. = N

Fig. 6. Virtual specimens geometry with three significant nodes highlighted.

Transient thermo-mechanical analysis have been performed on the three specimen geometries here considered; in Fig. 6, the final deformed profiles of the specimens are compared with the undeformed ones (dotted lines). The virtual specimens, meshed with plane axial-symmetric elements, have been subjected to the experimentally measured elongations imposing on the upper bound nodes the vertical nodal displacements v(t) computed by v(t) = δ(t)LS ,

(15)

where δ(t) is the measured engineering strain (function of time) and LS is the initial length of the specimen. The lower bound nodes are constrained to slide along the horizontal x direction only, representing the half-length symmetry plane, while the left bound nodes are constrained to slide along the vertical y direction only, representing the axial-symmetry axis. The temperature increment during the deformation has been computed considering that the entire plastic work is transformed into heat with attendant temperature increase; this leads to the following equation: ε σ dε PLWK T = = 0 , ρCp ρCp

(16)

where PLWK denote the plastic work done (per volume unit), ρ and Cp are the material density and specific heat, respectively. Clearly, this is a case of inverse procedure, where the actual unknown is not the structure behavior, but the material properties, that still have to be computed. In this paper, such a problem has been solved with an optimization approach. Generally speaking, the optimum design is the design that minimizes an objective function; in the specific case, the optimum design has been defined as the set of material coefficients that permits to minimize the standard deviation of the numerically obtained engineering stress–strain curve from the

(17)

where Snum (δ) and Sexp (δ) are the numerical and experimental engineering stresses, δ is the engineering strain and N is the number of simulation sub-steps that has been forced to be equal to the number of experimental sample points (1 ␮s time spaced). The numerical engineering stress is obtained by dividing the total reaction force in the lower bound nodes of the virtual specimen by the initial cross-sectional area. This is an unconstrained minimization problem, where the design variable space is constituted by the material coefficients, and the objective function is expressed by the average of the standard deviation computed for the three different geometry (D4L2, D4L4, D4L8): OBJ =

S.D.D4L2 + S.D.D4L4 + S.D.D4L8 . 3

(18)

The only minimization constraints are represented by the upper and lower limit of the variability range of the material constants. Running the analysis for different material coefficients, different OBJ values are obtained; the objective function is then approximated by a polynomial function of the design variables, as follows: OBJ(x) ∼ = a0 +

n i=1

ai xi +

n

bij xi xj ,

(19)

i,j=1

where n is the number of design variables x and the coefficient a0 , ai and bij are determined by a weighted least mean square fitting of the already computed OBJs. The polynomial function in (19) is the function to be actually minimized. As already outlined, the constitutive models illustrated in the previous section have been adopted, and the number of nondependent parameters has been reduced imposing the constraint that the equation should approximate the experimental data obtained from quasi-static (˙ε = 1 s−1 , T = Troom = 273 K) test; thus, only the temperature and strain rate related coefficients represent the design variables xi of the optimization procedure. Table 2 summarizes the situation about the optimization variables, while the constrained parameter values are reported in Table 3. A brief note for the Zerilli–Armstrong model is required: the last term in Eq. (13) represents the strain hardening, while the first two terms have been considered to correspond to the

M. Sasso et al. / Materials Science and Engineering A 487 (2008) 289–300

295

Table 3 Optimization results Quasi-static constraint

Variable range

Best set

Minimum OBJ (MPa)

Perzyna Clifton Johnson–Cook Litonski

A = 520, B = 269, n = 0.282 σ 0 = 525, εy = 0.0026, n = 0.0646 A = 520, B = 269, n = 0.282 σ 0 = 520, εy = 0.0027, n = 0.0666 C0 = 26, C5 = 269, n = 0.282

C = 14,030 ± 279, m = 0.292 ± 0.021 m = 0.0392 ± 0.0017, τ = 0.9124 ± 0.0599 C = 0.0476 ± 0.0021, m = 0.5530 ± 0.0250 a = 0.0021 ± 0.0003, b = 0.0687 ± 0.0116, m = 0.0676 ± 0.0037 C1 = 1703 ± 149, C3 = 4.53 × 10−3 ± 0.0011, C4 = 1.83 × 10−4 ± 0.00004

47 36 36 39

Zerilli–Armstrong

C: 6000–16,000, m: 0.1–1 m: 0.01–0.08, τ: 0.2–1.5 C: 0.09–0.23, m: 0.3–0.9 a: 0.0006–0.0025, b: 0.001–0.5, m: 0.001–0.5 C1 : 0–2000, C3 : 10−5 to 10−2 , C4 : 10−6 to 0−3

yield σ 0 of Eq. (11); so from the static test it is not found the value of C0 , but only its dependence on the other terms: C0 = 520 − C1 exp(−273C3 ). Of course, for a first evaluation of the coefficients in (19), a preliminarily investigation of the design variable space is needed. To have a good estimation of the shape of the object function, some simulations have been performed with material constant sets representing equally spaced points inside the design variable range; thus each variable range has been divided

± ± ± ±

1.5 2.0 2.0 1.5

35 ± 2.0

into 5 intervals for the two-parameter models and into 3 for the three-parameter models, requiring 52 and 33 simulations, respectively. Afterwards, the local minimum of the objective approximation function can be computed and the corresponding design variable set is reintroduced into the simulation code obtaining a new OBJ value, that is used to update the coefficients of function (19), and to compute a new local minimum. This is repeated iteratively until the new solutions are not able to reduce, with a relatively small tolerance (1%), the final OBJ value; gen-

Fig. 7. Objective functions visualization for the applied models.

296

M. Sasso et al. / Materials Science and Engineering A 487 (2008) 289–300

erally only two or three additional simulations are required to perform this iterative procedure. 7. Comparing results An overview of the principal results of the FE optimization is given in the next pictures and in Table 3. In Fig. 7, the objective functions are plotted as 3D meshed surfaces with some contour lines for the two-parameter models Johnson–Cook, Klopp–Clifton and Perzyna; the local minimums at the optimal design variable values, specified in Table 3, are clearly visible. The contour lines refer to isohypses of the 3D objective function in geometrical progression starting from the minimum + 4 MPa to minimum + 80 MPa. The shaded surfaces instead, refer to the three parameters Litonski and Zerilli–Armstrong models, and represent

Fig. 8. Numerical results for the optimum design set.

isohypses of the 4D objective functions. The surfaces are concentric encapsulated ellipsoids, grouping all the material coefficient combinations that give the same OBJ value, with the centre corresponding to best design set whose values are reported in Table 3. The black inner surface represents the minimum + 1 MPa isohypse of the 4D OBJ approximation function, while the others are 15 MPa equally spaced in increasing order. The optimization procedure seems to be stable with respect to the experimental data scattering, without significant variation of the iterations number. Further, the dispersion of the experimental data leads to a variability range of the optimal parameters (also reported in Table 3) that is seen to be quite small, except for the Zerilli–Armstrong model; this is probably due to the exponential nature of the law, that makes its parameters somewhat sensitive to the scattering of the experimental data they are requested to fit, but this does

Fig. 9. Stress state along the neck cross-section at 20% engineering strain.

M. Sasso et al. / Materials Science and Engineering A 487 (2008) 289–300

not jeopardy the overall matching of the numerical model with the real data, since the OBJ variation is observed to be in line with other models. For the best sets of material constants, the numerical engineering curves have been compared with the experimental ones (Fig. 8). It can be seen the good matching among the experimental stress–strain curve and all the simulated curves, except those obtained with the Perzyna and the quasi-static (denoted as “static” in the legend) parabolic hardening models; in particular the “static” curves remained constantly below the experimental curves for all the three geometries here considered, with a comprehensive standard deviation of 121 MPa, and showed a negative slope much lower than the experimental evidence. This is confirmed by the Perzyna curves that move up for increasing strain rate but, as for the “static” ones, they cannot accurately reproduce the stress decay because they does not account for any material softening. On the contrary, the models that account for both strain rate and thermal effects are able to provide a better matching of the experimental data, with the standard deviation indicated in Table 3. It means that the main phenomena occurred in the specimen have been correctly modeled; strong reduction of engineering stress after the initial peak cannot be explained only by the geometric softening due to necking, but a certain effective softening, probably due to temperature rising, must be admitted. Furthermore, the curves of Fig. 8 give a numerical “a posteriori” confirmation of what has been mentioned in Section 4: actually the increase of the experimental stress–strain curves is not only due to the strain rate effect. In fact, even in the simu-

297

lations performed with the parabolic hardening constitutive law where the strain rate effect is not really considered (Eq. (6a)), the nominal curves moves up as the specimen length decreases, with a maximum stress of 600, 640 and 730 MPa for the D4L8, D4L4 and D4L2, respectively; this means that also the stress distribution in the specimen contributes to the raising of the engineering curves. The situation is summarized in Fig. 9, where the stress components and the triaxiality (computed as the hydrostatic to Von Mises equivalent stress ratio) are plotted along the node A–node B necking section, as computed by the Zerilli–Armstrong model at 0.2 of nominal elongation. While Von Mises equivalent stresses are quite constant along the section path, with average values of 815, 791 and 763 MPa for the D4L2, D4L4 and D4L8 specimens respectively, the axial, radial and circumferential components are not uniformly distributed, especially in the shortest specimen which is characterized by a much higher triaxiality (>1.3 in the core). It is to note that in general the equivalent stress is lower than the average axial component, which is particularly high in the D4L2 specimen; this is another reason, beside the strain rate effect, why the shortest specimen follows the highest nominal flow curve. In Fig. 10 the temperature increase is shown as computed by the Zerilli–Armstrong model (the other models give analogous results) for the three specimen geometries at their maximum elongation; the D4L8 specimen, reaching lowest strain and strain rate levels, undergoes the lowest temperature increase, while in the other ones the temperature rises to relatively high values, about 200 K for the D4L4 and more than 230 K for the D4L2

Fig. 10. Temperature increase (K).

298

M. Sasso et al. / Materials Science and Engineering A 487 (2008) 289–300

specimen. These temperatures are enough to establish a sort of thermoplastic instability into the material, that happens when the thermal softening rate exceed the strain hardening rate; for a constant strain rate, this is mathematically expressed by       dσ ∂T ∂σ ∂σ + = ≤ 0, (20) dε ∂ε T ∂T ε ∂ε where dσ/dε represents the resulting slope of the true stress–strain curve, (∂σ/∂ε)T is the strain hardening rate at constant temperature, (∂σ/∂T)ε is the thermal softening rate and (∂T/∂ε) represents the temperature increase relevant to a ∂ε strain increment. The last term, as well summarized in [14], can be computed in closed form if the Johnson–Cook model is adopted with thermal softening exponent m set to unity; but generally the dependence of temperature on the plastic strain must be evaluated by numerically integrating Eq. (16). Fig. 11 shows the isothermal quasi-static true stress–strain curve compared with the adiabatic curves computed with the optimal constitutive models coefficients, and considering a 103 s−1 constant strain rate. It can be noticed how the dynamic

Fig. 11. True stress–logarithmic strain curves at 103 constant logarithmic strain rate.

curves, at relatively low strain, are much above the isothermal curve because of strain rate hardening, but after a critical strain level is reached, the material starts to soften; at higher strain, the dynamic curves go even below the static one. It must be said that

Fig. 12. Logarithmic strain and strain rate histories from simulation.

M. Sasso et al. / Materials Science and Engineering A 487 (2008) 289–300

the critical strain, at which the thermoplastic instability occurs, is differently modeled by the constitutive laws used: it varies from 0.075 for the Zerilli–Armstrong to 0.12 for the Litonski model. In Fig. 12 the strain and strain rate histories of the A, B and C nodes are reported, confirming that the true strain rate inside the specimen is much higher than the overall average value deduced from experimental engineering data (e.g. in the core of the D4L2 specimen the strain rate is almost permanently between 4000 and 5000 s−1 ). The curves also show how the strain rate distribution is widely irregular, highlighting very different histories among the A, B, and C nodes of the specimens individually considered. The strain rate profiles vary from geometry to geometry as well; in particular, while the intermediate and the longest specimens show similar behaviors, in the D4L2 the corresponding nodes have very different histories with the node B, that lies on the

299

necking cross-section, undergoing a lower strain rate than the internal node C. Because of all these issues, the numerical and experimental engineering curves can be regarded as the “rough” response of the specimen, made by the overall summation of the different contributes given by different parts of the specimen. The fact that the sample nodes here considered reach strain values very different among each other, experiencing very different strain rate levels even among a single geometry, actually means that they are forced to follow different effective true stress–strain curves during their deformation process. Moreover, because of the strain rate variation, the effective true stress–strain relationship followed by a certain elementary part of the specimen is not even represented by the adiabatic curves of Fig. 11, but it must be considered as the envelope of the different adiabatic curves corresponding to the different strain rate levels experimented by that elementary part of the specimen. The effective true stress–strain curves occurred at nodes A, B and C for each specimen, as numerically obtained with the Zerilli–Armstrong model, are given in Fig. 13. 8. Conclusions In this work, high strain rate tension tests have been performed on AISI 1018 CR steel specimens, by means of a SHPB device; afterwards, the same tests have been numerically reproduced with a finite element commercial code simulating the large deformation process underwent by the specimens. Analysis of the numerical results showed that the specimens are interested by widely non-uniform stress and strain distributions; moreover, different internal points of the specimen have been demonstrated to experience largely different strain rate levels at any given time, meaning that the engineering stress–strain curves represent only the “boundary conditions” of the specimen behavior, and they are not quantitatively representative of the material properties. An optimization procedure permitted to obtain the material coefficients of some common constitutive models, providing good correspondence between numerical and experimental nominal curves. In this way it has been possible to take into account for the coupled effects of the material, characterized by strain and strain rate hardening and thermal softening, and of the particular specimen geometries, that lead to a “geometrical softening” because of large reduction in area and simultaneously to an apparent hardening because of high triaxiality. These results encourage for future works aiming to calibrate, with a similar approach, a damage evolution model of the material, giving the possibility not only to characterize the flow stress–strain behavior but also to assess for the failure mechanisms. References

Fig. 13. True stress–logarithmic strain curves from simulation.

[1] B. Hopkinson, Philos. Trans. R. Soc. London A 213 (1914) 437–456. [2] H. Kolsky, Proc. Phys. Soc. London B 62 (1949) 676–700. [3] J. Harding, E.O. Wood, J.D. Campbell, J. Mech. Eng. Sci. 2 (1960) 88–96.

300 [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16]

M. Sasso et al. / Materials Science and Engineering A 487 (2008) 289–300 G.H. Staab, A. Gilat, Exp. Mech. 31 (1991) 232–235. J.L. Lewis, J.D. Campbell, Exp. Mech. 12 (1972) 520–524. D.J. Frew, M.J. Forrestal, W. Chen, Exp. Mech. 42 (1) (2002) 93–106. H. Kolsky, Stress Waves in Solids, Dover Publications, 1952. P.S. Follansbee, C. Frantz, J. Eng. Mater. Technol. 106 (1983) 61–66. D.A. Gorham, J. Phys. E: Sci. Instrum. 16 (1983) 477–479. H. Meng, Q.M. Li, Int. J. Impact Eng. 28 (2003) 537–555. E.D.H. Davies, S.C. Hunter, J. Mech. Phys. Solids 11 (1963) 155–179. I.D. Bertholf, C.H. Karnes, J. Mech. Phys. Solids 231 (1975) 19. D.A. Gorham, J. Phys. D: Appl. Phys. 24 (1991) 1489–1492. U.S. Lindholm, L.M. Yeakley, Exp. Mech. 8 (1968) 1–9. T. Nicholas, Exp. Mech. 21 (5) (1981) 177–185. G.B. Broggiato, F. Campana, M. Sasso, Proceedings of the 12th ICEM, 2004.

[17] M.A. Meyers, Dynamic Behavior of Materials, J. Wiley & Sons Eds, 1994, pp. 362–379. [18] J.A. Zukas, High Velocity Impact Dynamics, J. Wiley & Sons Eds, 1990, pp. 127–296. [19] G. Johnson, W. Cook, Proceedings of the 7th International Symposium on Ballistics, 1983, pp. 541–547. [20] R. Klopp, R. Clifton, T. Shawki, Mech. Mater. 4 (1985) 375–385. [21] J. Litonski, Bull. Acad. Polon. Sci. 25 (1) (1977) 1–8. [22] F. Zerilli, R. Armstrong, J. Appl. Phys. 61 (5) (1987) 1816–1825. [23] Ansys® Release 9.0 Documentation, Advanced Analysis Techniques Guide.