Numerical Simulation of Ductile Fracture in Cylindrical Tube Impacted Against a Rigid Surface

Abstract

Impact of tubes has an extensive application in crashworthiness. Analytical and numerical studies on axial impact of tubes have been carried out mostly at low-to-moderate impact velocities at which only the buckling of tubes is observed. However, experimental observations show that tubes fracture when the impact velocity is high. The objective of this article is to study the fracture of cylindrical tubes impacted against a rigid surface by numerical simulation. First, through-the-thickness fracture at the impacted end is simulated at the impact velocity of 350 m/s. Then, a detailed parametric study is carried out with respect to various parameters like the friction coefficient, tube thickness, tube length, and impact velocity. With an increase in the friction coefficient, the spread of fracture is more toward the outer edge. Further, the spread of fracture toward the inner edge takes place in an inclined manner. At higher tube thickness, the spread of fracture ceases after some time without reaching the inner edge. The tube length is found to have no effect on the fracture. It is observed that the fracture spreads faster with the impact velocity.

Keywords

tube ductile fracture continuum damage mechanics

INTRODUCTION

The response of tubes subjected to axial impact loading has been studied mostly at low-to-moderate impact velocities. When the impact velocity is not sufficiently large for the tube to fracture, the tube fails by buckling, that is called the axial crushing of the tube. The book by Jones (1989) describes, in detail, the response of tubes subjected to compressive axial impact. Depending on the geometry of the tube, the loading conditions, and the properties of the tube material, there are three types of bucking. If the length of the tube is too large compared to its diameter, the tube buckles in Euler mode, that is called the global bending. At lower impact velocities, the buckling is formed progressively from either of the tube end. This buckling is called the dynamic progressive buckling. At higher impact velocities, usually the entire length of the tube wrinkles before a large radial displacement develops in the tube. This is called the dynamic plastic buckling.

The problem of buckling of tubes has been analyzed numerically using commercial finite element codes by several researchers (Karagiozova et al., 2000; Karagiozova and Jones, 2001; Galib and Limam, 2004; Tarigopula et al., 2006). Most of the analyses have considered tubes of cylindrical cross-section. Further, the loading condition used is a stationary tube impacted by a moving rigid mass. The tube material is considered to be elasto-plastic strain hardening. Most of the studies have included the strain rate effects but not the thermal effects. Karagiozova et al. (2000) have used the finite element code ABAQUS/Explicit to simulate buckling of high-strength aluminum alloy tube at the impact velocities of 4–180 m/s. Karagiozova and Jones (2001) also have used ABAQUS/Explicit code to simulate buckling of steel and aluminum alloy tubes at the impact velocities of 5–100 m/s. Galib and Limam (2004) have used the finite element code RADIOSS to simulate buckling of Al–Mg–Si alloy tube at the maximum impact velocity of 60 km/h (i.e., 50/3 m/s). Tarigopula et al. (2006) have used the finite element code LS-DYNA to simulate buckling of square tube made of high strength steel alloy DP800 at the impact velocities of 5–15 m/s.

Only few studies have been conducted to study the fracture of tubes impacting against a rigid surface. Wang and Lu (2002) studied the impact of aluminum and mild steel tubes against a rigid surface at velocities of 137–385 m/s both experimentally as well as numerically using ABAQUS/Explicit code. The tube material is considered to be elasto-plastic strain hardening. The strain rate effect on the yield stress is incorporated using the Cowper–Symonds relation. However, the thermal effect is not included. Their experimental results show that, at lower impact velocities, thinner tubes buckle axisymmetrically, while thicker tubes mushroom at the impacted end before buckling. The important experimental result is that, at higher impact velocities, the tubes fracture, developing only a single through-the-thickness crack initially. Eventually, many such cracks develop leading to petaling. Their numerical results fall into three categories: (1) development of folds at the impacted end for thinner tubes at low impact velocities, (2) mushrooming and folds at the impacted end for all tubes at medium impact velocities, and (3) mushrooming at the impact end and wrinkling for thicker tubes at higher impact velocities. Wang and Lu (2002) have mentioned that they could not simulate the dynamic fracture using the ABAQUS/Explicit code. Galib et al. (2006) studied the impact of square aluminum tubes at velocities of 7–15.2 m/s both experimentally as well as numerically using RADIOSS code. The loading condition used is a stationary tube impacted by a moving rigid mass. The tube material is considered to be elasto-plastic strain hardening. The strain rate or thermal effects on the yield stress are not considered. Lemaitre's CDM model (Lemaitre, 1985) is incorporated in the analysis in an uncoupled manner. That is, the effect of damage on the material behavior is not included. However, the damage is used for predicting fracture initiation through the Lemaitre's damage growth law and the critical damage fracture initiation criterion. The experimentally observed fracture occurs at the corners of the impacted end. The predicted fracture location using the critical damage criterion matches with this result.

Literature review shows that there are many numerical simulations of the buckling of tubes under impact loading at low or moderate impact velocities. However, the simulation of fracture of tubes at higher impact velocities does not seem to have been attempted so far except by Galib et al. (2006). Hence, the objective of this study is to study the fracture of cylindrical tubes in impact problems by numerical simulation. First, through-the-thickness fracture at the impacted end is simulated at the impact velocity of 350 m/s. The predicted fracture pattern is found to be consistent with experimental results available in the literature. Then, the detailed parametric study is carried out with respect to various parameters like the friction coefficient (at the impact surface), tube thickness, tube length, and impact velocity.

Lemaitre's (1985) continuum damage mechanics (CDM) model is used to incorporate the effects of void nucleation and growth on the material behavior. The fracture initiation is assumed to take place when the damage reaches the critical value. The damage growth law, as well as the critical damage value, is obtained from the experimental results of LeRoy et al. (1981). The Johnson and Cook (1985) model is used to incorporate the strain rate and thermal effects on the yield stress of the material. The temperature rise, needed to estimate the thermal effects, is estimated by the transient heat conduction equation with the incremental plastic work (per unit volume) as the heat source and the frictional heat as the boundary condition at the impact surface.

DYNAMIC THERMO-ELASTO-PLASTIC FORMULATION

Incremental Strain–Displacement Relation

The incremental logarithmic strain measure, used in this formulation, is defined by (Bathe, 1996)

(1)

where δ_ij is the Kronecker's delta and

^{t + Δ t}_{t} λ_{i}

the principal incremental stretches.

Incremental Thermo-elasto-plastic Constitutive Relation for Damaged Material

While developing the incremental thermo-elasto-plastic constitutive relation for damaged material, it is assumed that the material is isotropic and the elastic behavior is linear. The von-Mises yield function (^tF₁) is used as the plastic potential in obtaining the plastic stress–strain relationship. The dependence of yield stress (^tσ_Y) on equivalent plastic strain ( $^{t} ɛ_{eq}^{pL}$ ), equivalent plastic strain rate ( $^{t} {\overset{\cdot}{ɛ}}_{eq}^{pL}$ ), and temperature (^tT) is approximated by a Johnson and Cook (1985) model. Further, it is assumed that the elastic ( $_{t} Δ ɛ_{ij}^{eL}$ ), plastic ( $_{t} Δ ɛ_{ij}^{pL}$ ), and thermal strains ( $_{t} Δ ɛ_{ij}^{TL}$ ) are additive. Thermal properties and some mechanical properties (Young's modulus E, Poissons's ratio ν, and density ρ) are assumed to be independent of the temperature. Then, the incremental thermal-elasto-plastic stress–strain relationship becomes

(2)

where ^tR_ij is given by

(3)

and the elastic-plastic constitutive tensor

^{t} C_{ijkl}^{EP}

is given by

(4)

Here,

C_{ijkl}^{E}

is the elastic constitutive matrix, ^tH ′ the derivative of ^tσ_Y with respect to

^{t} ɛ_{eq}^{pL}

, and ^tD the damage. In this study, the damage at a point at a time, assumed isotropic, is defined as the area void fraction in a plane at that point at that time. The tensor _tΔσ_ij on the left side of Equation (2) must be an objective incremental stress tensor. To make it objective, the stress-updating procedure discussed by Varadhan (1997) is followed. Radial backward return algorithm (Gautam and Dixit, 2010), which in this study has been modified to take into account the effects of damage, temperature and strain rate, is used for evaluating the stress in Equation (2).

Incremental Damage Growth Law

In this study, the damage growth law based on the experimental results of LeRoy et al. (1981) is used (Gautam and Dixit, 2010). The law, written in an incremental form, is given as

(5)

where c_d, a₁, and a₂ are the material constants obtained by fitting the above equation through the experimental results on damage versus equivalent plastic strain in tension tests by LeRoy et al. (1981). The quantity −^tY is the dissipative part of the thermodynamic force at time t corresponding to the rate of the damage variable ^tD. The expression for −^tY is given as (Lemaitre, 1985):

(6)

where E is the Young's modulus, ν the Poisson's ratio, ^tσ_eq the equivalent stress, and ^tσ_m the mean stress. The ratio

(\frac{^{t} σ_{m}}{^{t} σ_{e q}})

is defined as the triaxiality. The values of the constants are shown in Table 1.

Table 1.

Values of the material constants in the damage growth law (Equation (5)).

Material	c _d	a₁ (MPa⁻¹)	a₂ (MPa⁻¹)	D _cr
AISI1045	5.83 × 10⁻⁰³	4.04 × 10⁻⁰⁴	2.60	0.50

The procedure to find the values of c_d, a₁, and a₂ and the critical damage from the experimental results of LeRoy et al. (1981) is described in Gautam and Dixit (2010).

The above incremental damage growth law (Equations (5) and (6)) is used both for the positive triaxiality as well as for the negative triaxiality between (−1/3 ≤ (^tσ_m/^tσ_eq) ≤ 0). If the triaxiality (^tσ_m/^tσ_eq) is less than −1/3, the damage increment is made zero. Thus, the damage increment (at a Gauss point) is computed only when the triaxiality is greater than −1/3 as suggested by Bao and Wierzbicki (2005). However, no closure of the initially developed cracks is incorporated for the negative triaxiality. After calculating the incremental damage, the damage, at a Gauss point, is updated as ^t+ΔtD = ^tD + _tΔD, where ^t+ΔtD, is the damage at time t + Δt.

Integral Form of Equation of Motion

An integral form of the equation of motion at time t + Δt is used for the finite element formulation. It is given by the virtual work expression as discussed in Bathe (1996). Since, configuration at time t + Δt is unknown, following Bathe (1996), the integral is transformed to domain at time t. This is given, after some simplifications, by the following expression

(7)

where ^tρ is density at time t, ^t+Δtü_i the acceleration at time t, _tΔu_i the incremental displacement vector, _tΔε_ij and _tΔη_ij the linear and non-linear parts of the incremental Green–Lagrange strain tensor, respectively, δ the variation, and ^t+ΔtF the virtual work of the specified traction, contact traction, and body force.

Finite Element-finite Difference Scheme

The domain is discretized into a number of 8-noded brick elements and the incremental displacement field is approximated over each element by trilinear shape functions. Then Equation (7) leads to the following system of second-order ordinary differential equations:

(8)

where ^t[M ] is the global mass matrix, ^t[K] the global stiffness matrix, ^t{f} the global internal force vector, ^t+Δt{F} the global external force vector at t + Δt, ^t{F^T} the global external force vector due to thermal strain, and ^t{F^M} the global external force vector due to change in plastic properties. Further, _t{ΔU} is the global incremental displacement vector and ^t+Δt{Ü} the global acceleration vector at t + Δt.

To convert Equation (8) into a system of algebraic equations, the Newmark's implicit finite difference scheme as described in Bathe (1996) is used. Using this scheme and some algebraic manipulations, Equation (8) reduces to

(9)

where the effective stiffness matrix ^t[K_d] and the effective incremental force vector _t{ΔF_d} are given by

(10)

(11)

where a₀, a₁, and a₂ are the Newmark's constants and

^{t} {\overset{\cdot}{U}}

the global velocity vector. This algebraic system, (Equation (9)), is solved to obtain the incremental displacement vector _t{ΔU} which, however, becomes only an approximate solution of the governing equations. To minimize the error of the approximate solution, the modified Newton–Raphson algorithm (Bathe, 1996) is used. A scheme to remove spurious fluctuations in triaxiality within a finite element is used (Gautam and Dixit, 2010; Gautam et al., 2010).

THERMAL FORMULATION

The three-dimensional, transient, heat conduction equation is given as:

(12)

where it is assumed that the conductivity k, density ρ, and specific heat c do not depend on the temperature. The heat generated

^{t} \overset{\cdot}{q}

is considered as the sum of

{\overset{\cdot}{q}}_{p}

(heat generated per unit volume per unit time due to plastic deformation of interval Δt) and

{\overset{\cdot}{q}}_{f}

(heat generated per unit area per unit time at the contact surface due to friction):

(13)

The expression for

^{t} {\overset{\cdot}{q}}_{p}

is given by

(14)

where β is the fraction of plastic work converted into heat. Usually, β lies in the range 0.85 ≤ β ≤ 0.95 (Farren and Taylor, 1925; Gao and Wagoner, 1991). The expression for

^{t} {\overset{\cdot}{q}}_{f}

is given by

(15)

where ^tτ_s is the frictional or shear stress at the contact surface and ^tV_s the tangential relative velocity at the contact surface. The factor γ is taken as a function of the effusivities of the two materials in contact (Lin and Lin, 1990; Xing and Makinouchi, 2002).

The temperature field over the 8-noded brick element is also approximated by the same trilinear shape functions. This leads to the system of first-order ordinary differential equations. Using the direct integration method, the system of first-order ordinary differential equations is converted into a system of algebraic equations as

(16)

where

(17)

(18)

where ^t{T} is the global temperature vector, [C] the global specific heat matrix, [K]^Thermal the global conductivity matrix, and {Q} the global heat flux vector. Here, η is the parameter of the finite difference scheme which, in this study, is chosen as 0.5. Equation (16), which is a set of linear algebraic equations, is solved by the Gauss elimination method. Equation (16) is solved on the deformed configuration at time t + Δt. This is done after solving Equation (9) iteratively to determine _t{ΔU} and then updating the domain. The boundary condition on the free surface is the convective heat loss from the tube to the surrounding atmosphere.

CONTACT FORMULATION

The dynamic, large deformation, thermo-elasto-plastic, updated Lagrangian finite element formulation developed above can be used for contact analysis with appropriate modifications. The most important modification is the development of an additional set of finite element equations (involving the contact stiffness matrix) relating the unknown contact forces and displacements (Zhong, 1993). The contact force formulation in this study employs a node-to-surface interface model (Zhong, 1993), since this study deals with a large deformation problem including slipping at the contact surface. Lagrange multiplier method is used to enforce the contact constraints (Zhong, 1993).

RESULTS AND DISCUSSION

The developed formulation is first validated for static and dynamic cases by simulating the fracture in tension test of AISI1045 steel (Gautam and Dixit, 2010) and two fracture patterns (viz., tensile splitting and confined fracture) in the Taylor rod impact test (Gautam et al., 2010). Then, this formulation is used to study the through-the-thickness fracture in cylindrical tubes impacted against a rigid surface at the impact velocity of 350 m/s. Finally, the formulation is employed to carry out the parametric study with respect to various parameters like the friction coefficient (μ_f), tube length (L₀), impact velocity (V), and tube thickness (t₀).

Geometric and Material Properties of the Tube

The geometry of the tube is taken from Wang and Lu (2002) and is shown in Figure 1. The outer diameter of the tube (D₀) is fixed at 12.55 mm. The values of the tube length (L₀), friction coefficient (μ_f), impact velocity (V), and tube thickness (t₀) for various cases are mentioned in the respective sections. The number of element in the tube is 19008 and the number of nodes 22428. The numbers of element chosen along the length and thickness of the tube are 88 and 6, respectively. The material of the tube is AISI1045 steel.

As stated earlier, the dependence of the yield stress (^tσ_Y) on the equivalent plastic strain (

^{t} ɛ_{eq}^{pL}

), the equivalent plastic strain rate (

^{t} {\overset{\cdot}{ɛ}}_{eq}^{pL}

), and the temperature (^tT) is assumed to be governed by the following Johnson and Cook (1985) model:

(19)

Here, ^oσ_Y is the initial yield stress, K the strain-hardening coefficient, n the strain-hardening exponent, T_m the melting temperature, m^p the thermal-softening exponent, B^p the strain rate-hardening coefficient,

{\overset{\cdot}{ɛ}}_{ref}^{pL}

the reference strain rate, and T_ref the reference temperature. Table 2 shows the mechanical and thermal properties of AISI1045 steel (LeRoy et al., 1981; Jaspers and Dautzenberg, 2002) as well as parameters like T_ref, β, and γ.

Table 2.

Material properties of AISI1045 steel and other parameters.

E(GPa)	ν	ρ(kg/m³)	^oσ_Y(MPa)	K(MPa)	n	m ^p	B ^p
210	0.30	7800	302	796	0.59	1.0	1.0
ϵ^pL_ref	T _m (°C)	k	c (W/m°C)	α (J/kg°C)	T_ref (°C)	β	γ
1.0	1460	52	432.60	1.10 × 10⁻⁵	25	0.90	0.10

The values of the coefficients c_d, a₁, and a₂ in the damage growth law Equation (5) are given in Table 1. The critical values of the damage variable (D_cr) are estimated from the experimental results of LeRoy et al. (1981) on tension tests. The value of the damage at which the graph of area void fraction versus strain becomes almost vertical is taken as the critical value.

As stated earlier, the incremental damage growth law (Equations (5) and (6)) is used both for the positive triaxiality as well as for the negative triaxiality between (−1/3 ≤ (^tσ_m/^tσ_eq) ≤ 0). If the triaxiality (^tσ_m/^tσ_eq) is less than −1/3, the damage increment is made zero. Thus, the damage increment (at a Gauss point) is computed only when the triaxiality is greater than −1/3 (Bao and Wierzbicki, 2005). However, no closure of the initially developed cracks is incorporated for the negative triaxiality. Further, the damage increment is also made zero, if the damage has already reached the critical value (D_cr). After calculating the incremental damage, the damage, at each Gauss point, is updated as:

(20)

Then, the damage is averaged over each element. Next, the critical damage criterion is used to identify the failed elements. That is, if the average damage over any element reaches the critical value (D_cr), that element is considered as failed. Finally, the spread of fracture is simulated by tracing the path of the failed elements. While doing so, the failed elements are assumed to have zero stiffness.

Since the problems are axisymmetric, the damage distribution also becomes axisymmetric. However, as seen from experimental results of Wang and Lu (2002), fracture in the impact of tubes is not axisymmetric. This non-axisymmetric fracture occurs due to either material imperfection or imperfect tube geometry or non-perpendicular impact, etc. To simulate a non-axisymmetric fracture mode, the material imperfection, in the form of some small initial damage and equivalent plastic stain, is introduced at some locations. The material imperfection is introduced along some radial lines selected randomly. It is introduced in all the elements on both sides of the radial lines. An initial damage value of 0.03 and equivalent plastic strain value of 0.1 are assumed for these elements.

Ductile Fracture Simulation in Tube Impacted at Velocity of V = 350 m/s

In this section, ductile fracture at V = 350 m/s is studied. The tube dimensions are chosen as L₀ = 62.75 mm and t₀ = 0.78 mm. The coefficient of friction μ_f is assumed to be 0.05. Figure 2 shows the simulated tube with fracture at the impacted end. The failed elements (i.e., the elements where the damage has reached the critical value) have been removed.

Figure 3(a)–(f) shows the deformed configurations of the tube at various time steps. For the sake of clarity, only the mesh of 32 elements along the length is shown. As can be seen from the figure, till t = 4 µs, the deformation consists of only the mushrooming of the impacted end. At t = 8 µs, it is observed that both the inner and outer edges of the impact face lift. This is also clear from Figure 4 which shows the z displacement at the outer edge (i.e., point C of Figure 1) and at the inner edge (i.e., point A of Figure 1) at the impact face. It is clear that the outer edge of the tube lifts as soon as the tubes impacts the rigid surface. However, the inner edge lifts little later. This is also clear from the deformed configuration of Figure 3(b) which shows that the lifting is more at the outer edge. As the deformation proceeds, the material at the impacted face comes to rest in z-direction. However, the material above the impacted face still keeps moving down with much higher velocity and piles up on the material below (t = 8 µs). Hence, there is a bulging of the tube just above the impacted face. As the deformation progresses, the thickness of the tube wall increases and grows along the length of the tube (Figure 3(c)–(f)). As the damage at the impact face reaches the critical value, a crack initiates and grows as shown in Figure 3(e) and (f). In these figures, the failed elements (i.e., the elements where the damage has reached the critical value) have been darkened.

Figure 5(a)–(f) shows the progress of fracture along the thickness direction in r–z plane where r-direction is along the initial imperfection. In these figures, the left side is the inner surface, the right side the outer surface, and the bottom the impact face. Only six layers of the elements along the length are shown. Initially, the fracture starts at point B of the impact face (i.e., at radial distance = r₀ + t₀/2) (Figure 5(a) and (b)). As the deformation progresses (Figure 5(c)), the crack moves toward the outer edge. Later, it also moves toward the inner edge but in an inclined manner signifying a slant crack (Figure 5(d)). Finally, all the elements along the thickness direction upto certain height fail resulting in a through-the-thickness crack. The reason for the manner of propagation is explained in the next paragraph.

Figure 6 shows the damage growth at three locations on the impact face along a radial line at which the initial imperfection is provided. It is seen that the damage reaches the critical value first at point B (i.e., at radial distance = r₀ + t₀/2), next at point C (i.e., at the outer edge), and finally at point A (i.e., at the inner edge). Thus, the fracture initiates first at point B. The reason why the damage reaches the critical value first at B is that the growth of equivalent plastic strain is the fastest at point B (Figure 7). This happens, because this point experiences a lot of deformation due to the movement of the outer and the inner walls of the tube in opposite directions (Figure 8). Also, during the initial phase of the impact, the triaxiality at point B is much higher than that at point A or point C (Figure 9). This is because of the effect of stress waves propagating in axial and radial directions. Next, we explain why the crack propagates toward the outer edge. Initially, the growth of equivalent plastic strain is same at the outer and inner edges. However, when unloading takes place at the inner edge (Figure 10), the plastic strain stops to grow but continues to grow at the outer edge. Further, the triaxiality at the outer edge remains above the cut-off value of −1/3 for most of the time (Figure 9). This combination of the sustained equivalent plastic strain growth and the triaxiality being above the cut-off value leads to the propagation of crack toward the outer edge. At the inner edge, the equivalent plastic strain grows only after t = 10 µs (Figure 7) when the movement of the inner edge changes the direction from radially outward to radially inward (Figure 8). This leads to the propagation of the crack toward the inner surface after t = 20 µs.

Some idea of the fracture locus corresponding to the incremental damage growth law of Equations (5) and (6) can be obtained by plotting the graph of equivalent plastic strain versus triaxiality upto fracture at points A, B, and C of the impact face (Figure 1). This graph is obtained by combining Figures 7 and 9 and is shown in Figure 11. Fracture states of the three points are marked by open circles in the figure. For positive triaxiality, the fracture strain (i.e., the equivalent plastic strain at fracture) seems to decrease with triaxiality. However, none of the points has zero triaxiality at fracture. Therefore, it is difficult to say anything about the trend of variation of the fracture strain for negative triaxiality in the range of −1/3 ≤ (^tσ_m/^tσ_eq) ≤ 0.

Effect of Various Parameters on Tube Fracture

A detailed parametric study is carried out to study the effect of various parameters on fracture. The parameters selected are (1) tube length (L₀), (2) coefficient of friction (μ_f) between the tube and the rigid surface, (3) impact velocity (V), and (4) tube wall thickness (t₀). First, the parametric study with L₀ is carried out at three different values of L₀ namely, 62.750, 31.375 and 20.917 mm. In this case, the values of other parameters are chosen as: μ_f = 0.05, V = 350 m/s, and t₀ = 0.78 mm. Then, the parametric study with μ_f is carried out at four different values of μ_f namely, 0.0, 0.025, 0.05, and 0.10. In this case, the values of other parameters are chosen as: V = 350 m/s, t₀ = 0.78 mm, and L₀ = 62.75 mm. Then, the parametric study with V is carried out at three different values of V namely, 300, 350, and 400 m/s. In this case, the values of other parameters are chosen as: t₀ = 0.78 mm, L₀ = 62.75 mm, and μ_f = 0.05. Finally, the parametric study with t₀ is carried out at two different values of t₀ namely, 0.78 and 1.3 mm. In this case, the values of other parameters are chosen as: L₀ = 62.75 mm, μ_f = 0.05, and V = 350 m/s.

Effect of Length

It is seen that the damage growth at various locations on the impact face hardly changes with change in length. This shows that the fracture is a local phenomenon and is not affected by the length of the tube.

Effect of Friction

Figure 12 shows the damage growth for various values of friction coefficient at four locations on the impact face along a radial line at which the initial imperfection is provided. It is seen that, at the outer edge (i.e., at point C), the friction coefficient does not have much effect on the damage growth (Figure 12(a)). At the right of point B also, similar trend is observed except that, for μ_f = 0, the damage values are less than those at other values of μ_f (Figure 12(a)). When one looks at the damage growth at the left of point B (Figure 12(c)), it is found that the damage decreases significantly with an increase in the friction coefficient. The same trend is observed at the inner edge, i.e., at point A (Figure 12(d)). Next, we explain the reason for such a behavior by studying the growth of equivalent plastic strain and triaxiality.

The damage growth at point C does not depend much on the friction coefficient because both the equivalent plastic strain and triaxiality at this point (Figures 13(a) and 14(a)) do not change much with the friction coefficient. For the same reason, the damage growth at the right of point B also does not get affected much due to the friction coefficient except for μ_f = 0. For μ_f = 0, the damage values are less than those at other values of μ_f, because the triaxiality is less than the cut-off value of −1/3 for most of the initial periods (Figure 14(b)). At the left of point B, the equivalent plastic strain decreases with the friction coefficient (Figure 13(c)) while the triaxiality is less than the cut-off value for μ_f = 0.025, 0.05, and 0.1 (Figure 14(c)). This explains why the damage at this point decreases with the friction coefficient. At point A, the triaxiality is less than the cut-off value for the initial period (except for μ_f = 0) while the equivalent plastic strain decreases with the friction coefficient at higher values of t. Because of this, the damage at point A decreases with the friction coefficient.

Figures 15–17 show the progress of fracture along the thickness direction of initial imperfection for various friction cases, i.e., μ_f = 0, 0.025, and 0.1, respectively. The result for the case of μ_f = 0.05 is presented earlier in Section (Figure 5). It is seen that till t = 4 µs, no fracture initiation takes place for μ_f = 0 case. At t = 8 µs, only a single layer of elements fails for the case of μ_f = 0, while three layers of elements fail for the other three cases. As the time progresses, the fracture spreads toward the outer edge, the spread being more at higher values of the friction coefficient (see Figures 15(c), 16(c), 5(c), and 17(c) at t = 12 µs). At t = 16 µs, the fracture for the case of μ_f = 0 progresses in a symmetric manner while for the case of μ_f = 0.025 and 0.05 it progresses in symmetric manner only after t > 16 µs. However, for μ_f = 0.1, the fracture never spreads to the inner edge of the impact face and progresses in an inclined manner.

Effect of Impact Velocity

Next, we study the effect of impact velocity on the damage growth at the same four locations. It is observed that the damage growth increases with impact velocity. Hence, the plots of damage growth, equivalent plastic strain, and triaxiality are not presented.

Figures 18 and 19 show the progress of fracture along the thickness direction of initial imperfection for V = 300 and 400 m/s, respectively. The result for the case of V = 350 m/s is presented earlier in section (Figure 5). It is seen that for the impact velocity of 300 m/s also the fracture starts at point B. Then, it grows upward (Figure 18(b)) and later, it spreads to the outer edge (Figure 18(e)). No fracture is observed at the inner edge. For the impact velocity of V = 400 m/s, no fracture is observed upto t < 6 µs because the triaxiality is less than the cut-off value of −1/3. However, at t = 6 µs, the fracture initiates at a point between B and C (i.e., near the outer edge) and spreads toward the outer edge (Figure 19(b) and (c)). At t = 12 µs, the fracture reaches the outer edge and another fracture initiates near the inner edge at a point between A and B (Figure 19(d)). Finally, all the elements along the thickness direction upto a certain height fail resulting in a through-the-thickness crack (similar to the one observed for the impact velocity of V = 350 m/s). It is worth noticing that for V = 400 m/s, initially, the fracture growth is slow but later it grows at a much faster rate.

Effect of Tube Thickness

Figure 20 shows the damage growth for two different tube thicknesses at the same four locations. It is seen that, the damage decreases with the thickness at point C. This is because the equivalent plastic strain (Figure 21(a)) decreases with the thickness while the triaxiality does not change much with the thickness (Figure 22(a)). On the other hand, at two points around B, the damage increases with the thickness after t = 2 µs (Figure 20(b) and (c)). This happens in spite of the equivalent plastic strain decreasing with the thickness for a long time in the initial period at these two points (Figure 21(b) and (c)). This increase in the damage with the thickness is due to a very high value of the maximum triaxiality (of the order of 2.5) at these two points (Figure 22(b) and (c)). At point A, the damage growth for t₀ = 1.3 mm stops after t ≥ 15 µs (Figure 20(d)). This happens in spite of the triaxiality for t₀ = 1.3 mm exceeding the cut-off value of −1/3 after t = 7 µs (Figure 22(d)). The reason for this zero growth is that the growth of the equivalent plastic strain for t₀ = 1.3 mm also stops at around the same time (Figure 21(d)).

Figure 23 shows the progress of fracture along the thickness direction of initial imperfection for t₀ = 1.3 mm. The result for t₀ = 0.78 mm is presented earlier in section (Figure 5). As shown earlier for t₀ = 0.78 mm, the fracture starts at t = 4 µs at point B while for t₀ = 1.3 mm, the fracture initiates both at point B as well as at a point between A and B. As the time progresses (from t = 8 µs to 12 µs), the number of failed elements increases with the thickness. However, after t ≥ 20 µs, the spread of fracture for t₀ = 1.3 mm ceases without reaching the inner edge while for t₀ = 0.78 mm, the fracture spreads to the inner edge.

CONCLUSIONS

A damage coupled thermo-elasto-plastic formulation is employed for the simulation of ductile fracture in tubes of AISI 1045 steel impacted against a rigid surface. The Lemaitre's CDM model, along with a damage growth law based on the experimental results of LeRoy et al. (1981), is used for this purpose. The Johnson and Cook (1985) model is used to incorporate the strain rate and thermal effects on the yield stress. The specific conclusions of this study are as follows:

At the impact velocity of V = 350 m/s, fracture is observed at the impacted end. The impacted end mushrooms but no buckling is observed. Damage reaches the critical value first at the mid-point of the tube wall at the impacted end. Therefore, the fracture first starts at the mid-point of the tube wall, then spreads to the outer edge and then to the inner edge forming a through-the-thickness crack. The observed fracture pattern is found to be consistent with the experimental results on steel tubes reported in the literature (Wang and Lu, 2002).

Parametric study of various parameters:

Effect of tube length (L₀): It is found that the tube length has no effect on the damage growth and fracture at the impacted end.

Effect of friction coefficient (μ_f):

– At the outer edge, the friction does not have much effect on the damage growth. However, at the inner edge, the damage decreases with the friction.

– With increasing friction, the spread of fracture is more toward the outer edge of the impact face. Further, the spread of fracture toward the inner edge takes place in an inclined manner.

Effect of impact velocity (V):

– The damage growth and spread of fracture take place at a faster rate at higher impact velocities.

– For the impact velocity of V = 300 m/s, the fracture spreads till the outer edge but does not spread till the inner edge. For the impact velocity of V = 400 m/s, the second fracture initiation takes place at a point closer to the inner edge. The two fractures eventually coalesce to form a through-the-thickness crack.

Effect of tube thickness (t₀):

– The damage at the outer edge decreases with the thickness. At the inner edge, the damage growth at higher thickness stops after some time.

– As a result, the spread of fracture at higher thickness ceases after some time without reaching the inner edge.

Footnotes

Figures 3 and 11 appear in color online:

ACKNOWLEDGMENTS

The authors thank Rajkumar, Master's student of Mechanical Engineering at IIT Kanpur, for his continued support toward putting the manuscript in order and all the corrections.

References

Bao

Wierzbicki

(2005). On the Cut-off Value of Negative Triaxiality for Fracture, Engineering Fracture Mechanics, 72: 1049–1069.

Bathe

(1996). Finite Element Procedures, New Delhi, Prentice Hall of India.

Farren

Taylor

(1925). The Heat Developed During Plastic Extension of Metal, Proceedings of Royal Society of London. Series A, 107: 422–451.

Galib

Limam

(2004). Experimental and Numerical Investigation of Static and Dynamic Axial Crushing of Circular Aluminum Tubes, Thin Walled Structures, 42: 1103–1137.

Galib

Limam

Combescure

(2006). Influence of Damage on the Prediction of Axial Crushing Behavior of Thin-walled Aluminum Extruded Tubes, International Journal of Crashworthiness, 11: 1–12.

Gao

Wagoner

(1991). A Simplified Model of Heat Generation During the Uniaxial Tensile Test, Metallurgical Transactions A, 18A: 1001–1009.

Gautam

Babu

Dixit

(2010). Ductile Fracture Simulation in the Taylor Rod Impact Test Using Continuum Damage Mechanics, International Journal of Damage Mechanics, (doi: 10.1177/1056789509357119).

Gautam

Dixit

(2010). Ductile Failure Simulation in Spherodized Steel Using Continuum Damage Mechanics Coupled Finite Element Formulation, International Journal of Computational Methods, 7: 319–348.

Jaspers

SPFC

Dautzenberg

(2002). Material Behavior in Conditions Similar to Metal Cutting: Flow Stress in the Primary Shear Zone, Journal of Materials Processing Technology, 122: 322–330.

10.

Johnson

Cook

(1985). Fracture Characteristics of Three Metals Subjected to Various Strains, Strain Rates, Temperatures and Pressures, Engineering Fracture Mechanics, 21: 31–48.

11.

Jones

(1989). Structural Impact, Cambridge, Cambridge University Press.

12.

Karagiozova

Alves

Jones

(2000). Inertia Effects in Axisymmetrically Deformed Cylindrical Shells Under Axial Impact, International Journal of Impact Engineering, 24: 1083–1115.

13.

Karagiozova

Jones

(2001). Dynamic Effects on Buckling and Energy Absorption of Cylindrical Shells Under Axial Impact, Thin Walled Structure, 39: 583–610.

14.

Lemaitre

(1985). A Continuous Damage Mechanics Model For Ductile Fracture, Transactions of ASME: Journal of Engineering Materials and Technology, 107: 83–89.

15.

LeRoy

Embury

Edward

Ashby

(1981). A Model of Ductile Fracture Based on The Nucleation and Growth of Voids, Acta Metallurgica, 29: 1509–1522.

16.

Lin

(1990). An Investigation of a Coupled Analysis of a Thermo-elasticplastic Model During Warm Upsetting, International Journal of Machine Tools and Manufacture, 30: 599–612.

17.

Tarigopula

Langseth

Hopperstad

Clausen

(2006). Axial Crushing of Thin-walled High-strength Steel Sections, International Journal of Impact Engineering, 32: 847–882.

18.

Varadhan, S.N. (1997) Dynamic Large Deformation Elasto-plastic Analysis of Continua, Master's Thesis, Department of Mechanical Engineering, IIT Kanpur.

19.

Wang

(2002). Mushrooming of Circular Tubes Under Dynamic Axial Loading, Thin Walled Structure, 40: 167–182.

20.

Xing

Makinouchi

(2002). FE Modeling of Thermo-elasto-plastic Finite Deformation and its Application in Sheet Warm Forming, Engineering Computations, 19: 392–410.

21.

Zhong

(1993). Finite Element Procedures for Contact-Impact Problems, Oxford, Oxford University Press.