A finite deformation continuum\discrete model for the description of fragmentation and damage in brittle materials

https://doi.org/10.1016/S0022-5096(98)00027-1Get rights and content

Abstract

A dynamic finite element analysis of large displacements, high strain rate deformation behavior of brittle materials is presented in total Lagrangian coordinates. A continuum\discrete damage model capable of capturing fragmentation at two size scales is derived by combining a continuum damage model and a discrete damage model for brittle failure. It is assumed that size and distribution of potential fragments are known a priori, through either experimental findings or materials properties, and that macrocracks can nucleate and propagate along the boundaries of these potential fragments. The finite deformation continuum multiple-plane microcracking damage model accounts for microcracks within fragments. Interface elements, with cohesive strength and reversible unloading before debonding, between potential fragments describe the initiation of macrocracks, their propagation, and coalescence leading to the formation of discrete fragments. A surface-defined multibody contact algorithm with velocity dependent friction is used to describe the interaction between fragments and large relative sliding between them. The finite element equations of motion are integrated explicitly using a variable time step. Outputs are taken at discrete time intervals to study material failure in detail.

The continuum\discrete damage model and the discrete fragmentation model, employing interface elements alone, are used to simulate a ceramic rod on rod impact. Stress wave attenuation, fragmentation pattern, and overall failure behavior, obtained from the analyses using the two models, are compared with the experimental results and photographs of the failing rod. The results show that the continuum\discrete model captures the stress attenuation and rod pulverization in agreement with the experimental observations while the pure discrete model underpredicts stress attenuation when the same potential fragment size is utilized. Further analyses are carried out to study the effect of potential fragment size and friction between sliding fragments. It is found that compared with the continuum\discrete damage model, the discrete fragmentation model is more sensitive to the multi-body discretization.

Introduction

The study of ceramics under high velocity impact and penetration has gained prominence over the years owing to their high compressive strength and low density which make them more suitable for armor applications when compared to other conventional materials. There have been attempts to determine the fracture characteristics and threshold conditions for dynamic fragmentation of ceramics through particle impact experiments (Evans, 1978; Field, 1988; Shockey et al., 1990a; Andrews and Kim, 1997). Also, there have been numerous experimental studies to evaluate and characterize the performance of ceramics under various impact and boundary conditions (Shockey et al., 1990b; Espinosa et al., 1992; Hauver et al., 1993, Hauver et al., 1994; Woodward et al., 1994; Espinosa, 1995; Orphal and Franzen, 1996, Orphal and Franzen, 1997a; Orphal et al., 1997b, Xu and Espinosa, 1997; Brar et al., 1997). The outcome of these experiments has been that ceramics can be extremely efficient in defeating long rod penetrators as long as they are properly confined. In this case, high compressive strength and friction between fragmented pieces, in the comminuted zone, is achieved. The more efficient the confinement, the better the resistance of ceramics to penetration. However, many phenomena observed during deformation of the material under such high strain rate still remain unanswered, e.g., the effect of material processing route and grain size on penetration resistance, coupling between structural design and ceramic response, etc. These issues are expected to play a significant role in deriving useful design criteria and scaling experimental results to field application. Achieving these objectives by experiments alone is a difficult if not impossible task. Hence, recourse is sought to a combined experimental and numerical study.

The use of finite element analyses to determine the performance of ceramic targets under high velocity impact and penetration is relatively recent. The majority of the efforts up to now have been devoted to modeling the complicated response of ceramics in the post-failure regime, i.e., in the presence of cracks. Brittle failure in ceramics is the result of initiation and propagation of cracks, marked by fragmentation and stress attenuation. The various models available until now to describe this behavior can be categorized into two classes ; those based on a continuum approach, and those based on a discrete approach. The continuum approach (Johnson and Cook, 1985; Addessio and Johnson, 1989; Curran et al., 1990; Steinberg, 1992; Johnson and Holmquist, 1992; Rajendran and Grove, 1992a; Rajendran, 1992b, etc.) is based on homogenizing the cracked solid and finding its response by degrading the elasticity of the material. The fundamental assumption in these models is that the inelastic strains are caused by microcracks whose evaluation during loading degrades the strength of the material. This degradation is defined in terms of damage moduli whose evaluation under compressive, as well as tensile loading, is formulated using the generalized Griffith criterion. In addition, some of these models account for the initiation of cracks, coalescence, friction between fragments in the comminuted zone, etc. However, some of these phenomenological models have shortcomings in that they cannot describe damage induced anisotropy and also, that their parameters are difficult to measure experimentally. To overcome these limitations, (Espinosa, 1992, Espinosa, 1995) derived a multiple-plane microcracking model based on the micromechanics of solids. In this model, the constitutive response of the material is obtained from fundamental quantities that can be determined experimentally such as grain size and fracture toughness. In addition, since the dynamic growth of microcracks is described independently on each orientation, damage induced anisotropy and rate effects are naturally incorporated in the model. The model has been used extensively (Espinosa et al., 1996, Espinosa et al., 1997, Espinosa et al., 1998a) to study the high velocity impact response of ceramics under various conditions. In spite of these developments, continuum models have been criticized because they require assumptions on the size and distribution of microcracks to start with, and because they cannot describe the growth of dominant cracks leading to failure, which are not suitable for homogenization. Computationally, continuum models possess the problem of very large deformation under high strain rates, specially for higher order elements. As a result, elements distort and integration time step reduces below an acceptable level. One of the possible means to circumvent this problem is to use adaptive meshing (Camacho and Ortiz, 1996; Espinosa et al., 1996,Espinosa et al., 1998c) by which the elements in the highly deformed zone are remeshed into regular shaped elements.

On the contrary, models based on a discrete approach (Needleman, 1988; Biner, 1994; Lawn, 1994; Gong, 1995; Xu and Needleman, 1995; Camacho and Ortiz, 1996, etc.) nucleate cracks, and follow their propagation and coalescence during the deformation process. In this respect, the Camacho and Ortiz, 1996model is one of the most elaborate and powerful models currently available. In their model, conical and longitudinal cracks can nucleate at any node in a finite element mesh when the resolved normal or shear stress at the node exceeds an effective fracture stress. Cracks are nucleated by duplicating nodes and propagated, along the element boundary, by opening further nodes. Adaptive remeshing is used to provide a reach enough set of possible fracture paths around the crack tip. The forces at the cracked surfaces are brought to zero in accordance with the Griffith criterion accounting for unloading, if any, before reaching the critical fracture opening. Following this procedure, fragments are generated as cracks coalesce in a closed path. Thereafter, the model accounts for contact and friction between fragments. Realizing the computer resources required for carrying out three dimensional analyses, they carry out two dimensional axisymmetric analyses accounting for radial cracks based on a continuum approach. The major disadvantages of the discrete models are that they are extremely computer intensive and become infeasible as the number of cracks increases. At the limit, meshes with element size of the order of the grain size are required to capture all possible crack nucleating sites.

Besides all these developments, a fact remains that the experimental and numerical work have been carried out more or less independently of each other. There have been very few attempts (Espinosa et al., 1998b) to use meaningful experimental results to motivate the numerical modeling and vice versa. The experimental evidence (Simha et al., 1995; Espinosa and Brar, 1995; Espinosa et al., 1998b; Brar et al., 1997) is used herein in an attempt to simplify the modeling procedure while preserving all the physics of the problem. It is known that numerous cracks get generated inside grains and at grain boundaries in ceramics during the early stage of impact and penetration. Only a fraction of them grow and finally propagate, mostly along the grain boundaries, to form major flaws (henceforth referred to as macrocracks), while the rest of them get arrested during the deformation. The coalescence of these macrocracks in a closed path gives rise to fragments. This phenomenon is evident in Fig. 1Fig. 2Fig. 3Fig. 4

Fig. 1 shows the schematic arrangement of ceramic rod on rod impact together with an optical photograph of the damaged end of the rod. The figure shows fragments, bounded by cracks, whose average size has been measured to be more than 100 μm. In some cases, the major dimension of the fragments exceeds 400 μm. An SEM micrograph of the ceramic microstructure within a fragment is shown in Fig. 2. The micrograph confirms that each fragment contains numerous microcracks of length about 10–20 μm. Fig. 3 shows a schematic of yet another experiment carried out to investigate interface defeat in ballistic experiments. Fig. 4 shows a photograph of the cross-section obtained after cutting the penetrated TiB2 ceramic sample. The fragment size measured from the center line is found to vary in the range of 200–500 μm.

Microcracks that are not able to propagate and transform into a major crack can be modeled using a continuum approach so that the resulting degradation in the strength of the material is accounted for. A discrete approach can be used simultaneously to model crack coalescence and formation of fragments. To this effect, one simplification is made in the present analyses, based on the above experimental findings, by assuming that the macrocracks initiate only at potential fragment boundaries, i.e., the microcracks inside fragments do not propagate to the fragment boundaries and transform into macrocracks. This assumption makes it possible to model individual ceramic fragments with a continuum approach, and their boundaries by a discrete approach. The contact-interface element of Espinosa et al., 1997inserted along fragment boundaries can further eliminate the computation overheads in creating the microcracks and the remeshing to provide the crack path. The above assumption of macrocracks only along fragment boundaries and the existence of fragments containing a large number of microcracks, allows the analyst to consider a cluster of grains as a single fragment. Hence, analyses can be carried out with few hundreds to few thousands of fragments instead of millions of grains. The size and distribution of fragments can be determined from the initial grain structure of ceramics, distribution of macrocracks in recovered samples after shock and validation of the analysis results with the experimental results. This further concretizes the concept of a combined experimental–numerical study for such complex problems.

The following sections briefly describe the multiple-plane microcracking model of Espinosa, 1995extended to a large deformation formulation and the contact-interface model of Espinosa et al. (Espinosa et al., 1998c) including rate and temperature effects followed by the analyses of the ceramic rod on rod impact carried out based on discrete as well as continuum-discrete models. The comparison of the analyses results with the experimental results clearly marks the advantages and disadvantages of the two approaches and the effect of size of fragments and friction between fragments.

In our formulation, the field equations describing the material response of a body use a total Lagrangian approach. Considering a solid with volume B0 in the reference configuration, and a deformation process characterized by the mapping x (X, t), a material point initially at X will be located at x = X+u after deformation, in which u is the displacement vector, as shown in Fig. 5. A displacement based finite element formulation is obtained from the weak form of the momentum balance or dynamic principle of virtual work. At time t, the weak form is given byB0[∇0T00(b0a)]·ηdB0=0B0T0:∇0ηdB0−∫B0ρ0(b0a)·ηdB0−∫S0σt·ηdS0=0where T0 is the first Piola–Kirchhoff stress tensor at time t ; b0, a, and t are the body force vector, acceleration vector, and boundary traction vector on volume B0 and boundary S0σ, respectively. Virtual displacement field η is assumed to be admissible, and ρ0 represents the material density per unit volume in the reference configuration. The symbol ∇0 denotes the material gradient with respect to the reference configuration, and : is used to denote the inner product between second order tensors, e.g., A : BAi jBji, where the summation convention on repeated indices is implied. Another form of the weak form of the momentum balance, in terms of spatial quantities, is given byB0τ:∇sηdB0−∫B0ρ0(b0a)·ηdB0−∫S0σt·ηdS0=0in which superscript s stands for the symmetric part of the tensor, τ is the Kirchhoff stress and ∇ is the spatial deformation tensor. The Kirchhoffs stress is related to first Piola–Kirchhoff stress T0 and second Piola–Kirchhoff stress tensor S by τ = FT0 = FSFT, where F is the deformation gradient at time t.

The discretization of Eq. (3)defines a system of nonlinear ordinary differential equations which can be solved for the updated deformation xn+1. A displacement-based finite element formulation is obtained by expressing field variables at any point in an element as a function of nodal quantities and the element shape functions in the reference configuration. For the two dimensional analyses considered herein, they are given as :uia=a=1NENNa(ξ,η)uaivia=a=1NENNa(ξ,η)vaiaia=a=1NENNa(ξ,η)aaiwhere uia, via, and aia are the displacement, velocity, and acceleration at the point of interest. Na (ξ, η) are the shape functions relating nodal quantities uai, vai, and aai to the point of interest and the summation is taken over the number of nodes in the element, i.e., NEN. The vector ξ, η contains the natural coordinates of the point of interest in the isoparametric master element. For the triangular elements, these natural coordinates are related to the area coordinates as L1 = ξ, L2 = η and L3 = 1−ξη.

Substitution of the discretized variables into Eq. (3)leads to the following system of ordinary differential equations,Miajbajb=fextia−fintiawhere,fintia=eBeona,jτijdB0=eBeoBTτdB0fextia=eBeoρ0boiNeadB0+∫∂Seot̄iNadS0Miajb=eBeoδijρ0NebNeadB0in which fintia, fextia, and Miajb are the internal nodal forces, the external nodal forces, and the lumped mass matrix, respectively. In the above equations, summation is taken over all the elements in the mesh according to their nodal connectivities.

The large deformation continuum response of ceramics in the presence of cracks is described through a microcracking multiple-plane model based on a dilute approximation (Taylor model ). The formulation is an extension of the small deformation multiple-plane microcracking model given in Espinosa, 1995, assuming that the displacements and rotations are large while the strains are moderate. This assumption allows us to adopt the rate form of the work conjugate second Piola–Kirchhoff stress tensor and Green–Lagrangian strain tensor with an additive decomposition of the later into an elastic and inelastic part, i.e.,Ėij=Ėeij+Ėcijwhere Ėe is the elastic part of the Lagrangian strain rate, and Ėc is the inelastic part arising from the presence of microcracks within the solid. Based on this decomposition, the elasticity law can be formulated as ;Ṡij=LijklĖeklin which Li jkl is the fourth order material stiffness tensor.

The basic assumption in the microcracking multiple-plane model is that microcracking and\or slip can occur on a discrete number of orientations (see Fig. 5) Slip plane properties ( friction, initial size, density, etc.) and their evolution are independently computed on each plane. In contrast to scalar representations of damage, the present formulation is broad enough to allow the examination of damage induced anisotropy and damage localization in the interpretation of impact experiments.

For a representative volume B0, in the underformed configuration, of an elastic solid containing penny-shaped microcracks with a density N (k) , the average Lagrangian inelastic strain rate tensor can be written as ;Ėcij=k=113N(k)π(a(k))212(b̄ḃ(k)in(k)j+n(k)ib̄ḃ(k)j)+k=113N(k)πa(k)ȧ(k)(b̄(k)in(k)j+n(k)ib̄(k)j)where the subindex k is used to label the orientations, a (k) denotes the radius of a microcrack on orientation k, n (k) the corresponding unit normal, and b̄ (k) the average displacement jump vector across A (k) .

If the resolved normal traction acting on the microcracks on orientation k is tensile, the average displacement jump vector resulting from an applied second Piola–Kirchhoff stress field S is given byb̄(k)i=1A(k)A(k)b(k)idA=16(1−ν2)3E(2−ν)a(k)(2Sijn(k)j−νSjln(k)jn(k)ln(k)i),in which E and ν are the Youngs modulus and Poissons ratio of the uncracked solid, and ak is the radius of the penny-shaped microcracks on orientation k. By contrast, if the normal traction is compressive, the microcracks are closed and the average displacement jump is given byb̄ki=32(1−ν2)3πE(2−ν)akfkiwhere fk is the effective shear traction vector on orientation k given byfki=(Skτ+μSkn)(nτ)ki.In the above equation, μ is the friction coefficient of the microcrack faces, Skτ and Skn are the resolved shear stress and the normal stress acting on microcracks with orientation k, respectively, and nkτ is the unit vector in the direction of the resolved shear traction. Embodied in Eq. (12)is the notion that fk provides the effective driving force for the sliding of the microcracks. If τ(k)μσ(k)n a sticking condition occurs. Hence, irreversible displacement jumps develop making the deformation process strongly nonlinear and history dependent.

In order to compute the inelastic strain tensor at all times, it becomes necessary to follow the evolution of the microcrack radius, ak, in the selected orientations. Following Freund, 1990, an equation of evolution for a in the case of mixed mode loading can be derived, Espinosa et al., 1992, Espinosa, 1995, viz.,ȧk=m±cR[1−(KIC\Kkeff)]≥0in which n± and m± are phenomenological material constants which may have different values in tension and compression, cR is the Rayleigh wave speed, KIC is the material toughness, and Kkeff(S,ak) is an effective stress intensity factor which is a function of the stress state S and the microcrack size a (k) . For mixed mode conditions, Kkeff is derived by considering an average energy release rate associated with an increase in radius of the microcracks, namely,Gk=101−ν2E[K2I+K2II+K2III\(1−ν)]dθfrom which the following expression for Kkeff is obtained,Kkeff=GkE1−ν2.

The general structure of these constitutive equations corresponds to that of a solid with a damage-induced anisotropic stress–strain relation with elastic degradation. In particular, the effective behavior of the solid is predicted to be rate dependent due to crack kinetics effects. From a computational standpoint, this insures numerical reliability and mesh independence, according to Needleman, 1988, and Espinosa, 1989. This is in contrast to quasi-static formulations of damage for which the governing equations become ill-posed in the softening regime, as in Sandler and Wright, 1984. Details about the stress update algorithm, assuming the above additive decomposition of strain rates into an elastic and cracking part, can be found in Espinosa, 1995.

It should be pointed out that this inelastic model is a continuum model in which material damage results from microcracking. If the material is subjected to a predominantly tensile stress state, microcracks along orientations perpendicular to the direction of maximum tensile stresses will grow according to , . In this case, significant dilation is expected due to mode I crack opening. If a predominantly compressive state of stress with shear is imposed, then crack opening is inhibited but inelasticity is manifested by the growth of penny-shaped cracks in modes II and III (shear modes).

The continuum model is integrated with a dynamic contact\interface model for predicting the discrete response of the material with regard to the cracks initiating at and propagating along the fragment boundary. A versatile multi-body contact model for explicit dynamic analysis has been developed. The contact algorithm is based on predicting accelerations assuming no contact and later correcting the accelerations of the surface nodes so that the surfaces do not interpenetrate. In addition, a velocity dependent friction model is included. This surface-based contact algorithm allows contact between bodies that undergo large relative displacements as they move. It also allows the easy incorporation of a velocity-dependent friction model where the friction coefficient is made a function of pressure and temperature. A detailed description of the contact algorithm can be found in Espinosa et al., 1996, Espinosa et al., 1998c).

The contact algorithm is augmented with interface elements to simulate the initiation of cracks and subsequent large sliding, opening and closing of the cracked surfaces. The model is based on the interface model proposed by Tvergaard, 1990for quasi-static calculations. It assumes that a perfect interface between two surfaces carries forces that oppose separation and shear between them until debonding. After debonding, the two surfaces will behave as distinct identities. The propagation of a crack can thus be simulated as the consecutive failure of interface elements between potential fragments. The magnitude of the opposing forces before debonding is a function of the relative normal and shear displacement jumps between the two surfaces. The normal and tangential tractions are given as,Tn=unδnF(λc),TtutδtF(λc)F(λc)=274Tmax(1−2λc2c),for0≤λc≤1where, Tmax is the maximum interface traction and λc is a non-dimensional parameter representing an effective normalized displacement jump given by,λc=unδn2+utδt2.In the above equation, un and ut are the normal and tangential displacement jumps at the interface. δn and δt are critical displacement jump values at which the interface breaks, i.e., crack initiation. Tmax, δn and δt can be easily determined from energy arguments based on the mode I critical strain energy release rate of the interface. It is evident that the value of λc varies from 0 to 1 with λc = 1 defining interface failure. Four node quadrilateral elements are embedded as the interface elements between fragments. The normal and tangential forces are computed depending upon the state of stress at the interface, as shown in Fig. 6 and Table 1.

As long as the value of λc remains less than unity, the interface normal traction is computed from the contact model in the case of compression-shear and from the interface model in the case of tension-shear. The shear traction is computed from the interface model. These normal and tangential tractions contribute to the internal nodal forces so that the fintia of Eq. (6)gets modified as ;fintia=eBeoBTτdB0+finteia=eBeoBTτdB0+inteSinteeNTTdSwhere finteia is the force contributed by interface element to node i which is calculated by integrating the interface tractions over the midsurface of the element. For values of λc larger than unity, signifying initiation of crack, only the contact model is used to compute the interface tractions in both states. The dependence of the interface model on the materials properties and switching from interface to contact and vice versa make the approach realistic to simulate crack initiation, unloading if any before opening of cracks and subsequent large sliding at the newly created surfaces.

Two possibilities exist for the description of interface friction unloading. The first, based on the existence of a potential, describes unloading following the loading path. The second, based on an irreversible process, may describe unloading connecting the maximum displacement jump and the origin. The irreversible interface law (Pandolfi and Ortiz, 1997) could be used in problems where interface damage occurs. In this law, the maximum displacement jump is used to represent history effects. Both possible interface laws are schematically shown in Fig. 7. The reversible model is used in the calculations reported later.

Rate and temperature effects in the interface description can be easily incorporated in terms of the λc parameters as given by,Tmax=T0max1+βlnλ̇cλ̇c,01−[T−T0][Tm−T0]xwithT0max=48Gc27δn.In the above expression, Tmax is the maximum interface traction at the current displacement jump rate λc and current temperature T, T0max is the maximum interface traction at reference displacement jump rate λc,0 and reference temperature T0, and Tm is the melting point of the material. The parameter β and x can be identified through specially designed experiments. Since there are three parameters and only one energy equation, a proper characterization requires impact experiments, i.e., spallation and compression shear, in which the cohesive law can be properly characterized. Similarly, the functional dependence of the interface traction on displacement jumps and temperature needs intensive experimental characterization.

An explicit central-difference integration algorithm is being used to integrate the system of ordinary differential equations in time. The algorithm, accounting for acceleration corrections due to contact and equivalent nodal forces arising from interface elements, is summarized in Table 2. As in any initial boundary value problem, initial displacements and velocities u0 and v0 are required. Initial accelerations a0 are calculated from initial applied forces fext0, and initial internal forces fint0 from the equation of motion at time t.

At each time step n, the nodal accelerations must first be corrected for any time-dependent changes in the traction boundary conditions. Then, a displacement predictor at time n+1, is computed using the corrected acceleration, displacements and velocities at time step n. Modified accelerations at time n are computed based on the corrected acceleration and changes in accelerations resulting from surface contact\interface determined from the displacement predictor at n+1. Updated displacements at n+1 are used in the update of stresses and the computation of internal forces. Lastly, accelerations and velocities at time n+1 are obtained completing the time integration scheme.

This explicit integration method is very useful for studies in which high rates of loading are expected. The time steps used by these explicit calculations are limited by stability, so care must be taken in finite deformation dynamic calculations to ensure that waves do not propagate through the mesh faster than the material wave speeds. To this end, the time step is calculated dynamically from the maximum element frequency in the mesh ωmax, such that Δt ≤ 2\ωmax. Flanagan and Belytschko, 1984, derived the following estimate of ωmax for an N-noded isoparametric element,ω2max≤Nλ̂+2μ̂ρBiIBiIA2in which BiIBiI is the trace of [B] [B] T, and the area A is found as CIJXIYJ, whereCIJ=∫∫A∂NI∂ξ∂NJ∂η∂NI∂η∂NJ∂ξdηdξ,and NI are the element shape functions.

After developing a finite element model for large deformation brittle failure, validation by simulating specially designed experiments is required. Moreover, accessing its ability to capture the observed material inelasticity and failure mechanisms is desirable. Espinosa and Brar, 1995and Simha et al., 1995studied damage and fragmentation of AD-998 and AD-94 sintered alumina bars manufactured by Coors Porcelain Company, Golden, CO. The experiments consisted of the co-axial impact of a ceramic rod against another rod made of identical materials with diameter 12.7 mm. The length of the impactor rod was 80 mm, while that of the target rod was 170 mm. Moreover, in-material stress measurements were performed with manganin gauges (Type C-880113, Micro-Measurements, Raleigh, NC) embedded at 10 rod diameters from the impact surface. The gauges were backed by 25.4 or 50 mm rods of the same material as the front piece. The back pieces of the assembled bar targets were set in a lexan disk using epoxy in order to align the target for a planar impact. A coaxial trigger pin was also set through a hole in the lexan ring to trigger a manganin gauge bridge circuit and a high speed Imacon camera operated at 105 frames\s. The bars were painted black so that cracks and faults could be distinguished during the failure event. Measured stress profiles revealed the material remains elastic when the impact velocity is below about 100 m\s. By contrast, when the impact velocity reaches 304 m\s, the maximum compressive stress is sustained for only 200 μs with a main pulse duration of approximately 1.8 μs followed by a tail. Damage evolution was observed in a sequence of photographs taken every 10 μs with the Imacon camera. Well defined patterns of longitudinal cracks are also observed in a region approximately two bar diameters in size, symmetrically formed from the impact surface. A violent radial expansion with a debris of fine particles results from accumulated damage and pulverization of the material.

The results obtained from the rod on rod experiments are used to validate the numerical model described earlier. The model is further used to carry out parametric analyses and establish the adequacy of the proposed model. The dimensions of the two ceramic rods used in the analyses are the same as mentioned above. All the analyses are carried out at an impact velocity of 304 m\s.

As stated earlier, it is important to measure the initial grain size of the ceramic material to estimate the smallest possible fragment size (cluster of grains) based on experimental observations. Based on this estimate, fragments are generated in the ceramic bars with random morphology. The application of the model to actual grain sizes and fragments is the subject of current research work. For discretization, it is assumed that potential fragments are of quasi-hexagonal shape with random size and their major axes are aligned parallel to the axis-of-symmetry. In order to minimize the computation time, potential fragments are generated only in a portion of the length of the two bars where severe fragmentation was observed. The rest of the bars length is discretized continuously without potential fragments.

A pre-processor has been developed and integrated in our version of FEAP (Zienkiewicz and Taylor, 1991a, Zienkiewicz and Taylor, 1991b) for generating potential fragments and the interface elements between them. The bar radius, the axial length and the number of fragments to be generated along the two directions, i.e., X- and Y-, are given as input to the pre-processor. The pre-processor calculates the average fragment size based on the input and generates the hexagons. Subsequently, the vertices of the hexagon are displaced randomly maintaining the outer boundaries. The direction and the amount of displacement for each node are generated using a standard random number generator. Each hexagon is then discretized into six 6-noded quadratic triangular elements. On account of symmetry at the center and to maintain the continuity of the material, partial hexagons are generated at the boundaries which are discretized into two or three triangles. Four noded interface elements are embedded along the edges of hexagons by picking two edge nodes from the hexagon on the right side of an edge and remaining two nodes from the hexagon on the left side of the edge.

The quadratic triangular elements in the continuous region of the two bars are generated using a pre-processor already in-built in this version of FEAP. One additional module has been developed to connect the hexagons and continuous regions by embedding interface elements between them.

Two types of meshes have been used in the analyses. As shown in Fig. 8 with zooms centered at the impact plane, one is a coarse mesh having 100 hexagons, 2713 nodes, 923 triangular elements and 510 interface elements. A fine mesh shown in the same figure having 304 hexagons, 7464 nodes, 2564 triangular elements and 1656 interface elements was also used.

The properties of the alumina ceramic and the parameters for the contact and interface elements used in the various analyses are given in Table 3Table 4Table 5. The ceramic multiple-plane continuum damage model parameters are selected such that a damage threshold is properly captured and, at the same time, they are in agreement with values reported in the literature. As pointed out in Espinosa and Brar, 1995, it should be noted that microstructural differences have an effect on model parameters. For instance, the existence of a second phase at the grain boundaries controls parameters μ, n±, m± and KIC. The grain size controls the values of a0 and maximum crack density of each orientation. The choice of different crack densities as a function of orientation is motivated by experimental data (Espinosa et al., 1992) indicating that the density of active microcracks is a function of the mechanism responsible for crack nucleation. In the case of predominantly compressive stress states, glass at the grain boundaries determines the early stages of inelasticity and acts as a precursor for the development of triple point microcracks. In the case of predominantly tensile stress states, cracks are mainly nucleated due to grain boundary decohesion at a much lower stress level. In the present axisymmetric calculations, plane 1 is a plane on which positive normal tractions result from unloading waves emanating from the bar surface. By contrast, planes 6 and 11 are mainly shear cracks.

The interface parameters were selected consistently with the assumed fracture toughness of the material, i.e., KIC = 1.7 MPa m or Gc = 78.75 N\m. Assuming a δc = 1 μm, a quasi-static Tmax = 140 MPa is obtained. In the present analyses, a reversible law is used and no rate or temperature effects are included due to the limited experimental data available.

Two dimensional axisymmetric analyses for the rod on rod impact have been carried out to validate the model with the experimental results and to study the effect of various parameters. The continuum-discrete approach as mentioned earlier, as well as the discrete approach using the interface elements, are employed to analyze the hexagon region of the two bars, while only continuum damage is used to analyze the continuous region. The effect of fragment size is analyzed by carrying out analyses with two different mesh sizes as given earlier. Analyses have also been carried out with three different values of the friction coefficients between fragments to study their effect on the response during deformation. The various analyses carried out are summarized in Table 6.

As mentioned earlier, the continuum damage model possesses the problem of very large distortion of elements which has been circumvented before by adaptive meshing of the highly deformed zone. In the present case too, the fragments deform excessively when the continuum\discrete model is used. The remeshing of grains and reinsertion of interface elements between them due to mesh changes are under development. The problem has been circumvented for the present analyses by using erosion (Espinosa et al., 1998a). Fragments in which all elements attain a critical value of erosion parameter at their integration points are eroded. The value of 2.0 for the effective inelastic strain considered as erosion parameter has been used in all the analyses with the continuum\discrete model.

As given below, the results obtained from various analyses are discussed in terms of the continuum-discrete vs discrete damage models, effect of fragment size and friction and the fracture patterns predicted with the two models. The results obtained from the analyses are compared with the experimental results in terms of the measured in-material stress histories and the recorded fragmentation pattern with high speed photography.

As mentioned in the previous section, failure in ceramics is the result of initiation and propagation of cracks. As a result of which, stress attenuation and fragmentation take place. In general, stress attenuation can be described reasonably well by means of continuum models. However, due to the inherent assumptions, they cannot describe crack coalescence leading to discrete fragmentation. Hence, a complete description of brittle failure can be obtained either by means of discrete fragmentation models or continuum\discrete models.

In order to examine the features of continuum\discrete and discrete models in predicting the response during impact, analyses are carried out with a coarse as well as a fine mesh. The analyses with the coarse mesh have been carried out without friction between fragments, while a value of friction coefficient of 0.2 is used in the case of the fine mesh. The axial stress at 120 mm from the impact plane is retrieved from analyses results. This stress is compared with the experimentally measured stress and is plotted in Fig. 9Fig. 10 for the two types of meshes, respectively. An analysis has also been carried out with the continuum damage model alone for comparison purposes. Its stress history is also given in Fig. 9.

It is found that the discrete model reasonably predicts the stress increase and its peak value, but overpredicts the pulse duration and fails to capture the strong stress attenuation observed experimentally. By contrast, the pure continuum model provides a good overall fit for the set of parameters previously reported. However, the continuum\discrete model predicts the pulse duration and shock attenuation with reasonable accuracy but underestimates the peak stress. It is not surprising that the continuum\discrete model results in a lower peak stress than the one obtained with a pure continuum damage model. In fact, since the same crack density is used in both cases, the continuum\discrete model possesses a larger crack surface area per unit volume in view of the point that, in addition to the microcrack density used in the multiple-plane model, macrocracks are allowed in the formation of fragments. It can be inferred from the plot that the continuum\discrete model can capture not only the measured stress profile, but also the pulverization observed experimentally in the impact region (see later discussion on fragmentation pattern and ceramic pulverization near the impact region).

An important aspect of the failure of ceramics is the role of friction between fragments. It is known that the shear strength of comminuted ceramics is controlled by the friction properties between fragments. For this reason, confinement pressure also play a dominant role in the material shear resistance. In the case of rod on rod impact, the confinement pressure is generally low. However, friction between fragments is important at the early stages of damage.

In order to investigate the effect of friction between fragments, analyses have been carried out with three values of the friction coefficient as 0.0, 0.2 and 0.5 using both the continuum\discrete and discrete models and the coarse mesh. Axial stress histories predicted by the two models with the 0.2 and 0.5 friction coefficients together with that obtained with the zero coefficient discussed above are plotted in Fig. 11Fig. 12Fig. 12 respectively. The fragmented bars at 10, 20, 30 and 40 μs are shown in Fig. 13 for the three coefficients obtained by the discrete model. A comparison with high speed photographs, at these times, taken during the experiment is also shown in the figure. The experimental photographs clearly show longitudinal cracks and an increasing radial expansion, near the impact surface, leading to the pulverization of the ceramic.

As observed in Fig. 11, Fig. 12 the friction coefficient has a lesser effect in the case of continuum\discrete model compared to the discrete model. The difference in the peak stress predicted by continuum\discrete model is about 0.8 GPa when the friction coefficient is increased from 0–0.5. But the pulse duration and stress attenuation are not altered significantly. In the case of the discrete model, the peak stress increases by about 2.2 GPa as the friction coefficient is increased to 0.5. Also, the pulse duration increases and the stress attenuation is reduced significantly. As the value of friction coefficient is definitely larger than zero, the response predicted by the discrete model is in sharp contrast to that obtained experimentally. It should be noted that in the case of the continuum\discrete damage model, a significant fraction of the total dissipated energy occurs within the fragments due to the growth of a large density of penny-shaped cracks. In the calculations, the internal friction coefficient used in the description of the multiple-plane model was not modified.

It is observed in Fig. 13 that the discrete model predicts overall fragmentation consistently with the experimental results only at low values of μ. While the discrete model is more sensitive to the friction coefficient between fragments, it does not require fragment erosion to advance the calculation. As a result, a better understanding of the evolution of fragment initiation and evolution is attained. It is observed that the amount of radial expansion and rotation in the fragmented zone is very much reduced with the increase in friction coefficient. The response predicted by the model when the friction coefficient is set to 0.2, matches reasonably well with the fragmentation progression as recorded by high speed photographs. Similar results have been obtained with the continuum\discrete model (not plotted here) with the difference that it requires a small amount of fragments erosion to advance the solution.

It is also important for the model to predict fracture patterns leading to failure in agreement with the experimental findings. For this purpose in this subsection, the fracture pattern obtained from the discrete and continuum\discrete models are given in Fig. 14Fig. 15Fig. 16Fig. 17Fig. 14, Fig. 15, Fig. 16, Fig. 17 respectively, obtained with the coarse and fine meshes at the early stages of fragmentation, first 12 μs. It is observed, in both models, that fracture of ceramic rods starts from the outer surface and propagates towards the core and away from the impact surface. This is consistent with the wave propagation theory in inelastic media. Together with the compressive shock wave generated at the impact surface, tensile rarefaction waves are generated at the rod free surface which tail the shock wave with the sound speed of the ceramic material. These rarefaction waves are primarily responsible for the initiation of cracks near the surface. So long as the rarefaction wave from the other end of the bars has not arrived to the hexagon region, failure of both bars, impactor and target, remains more or less similar in pattern. However, due to arrival of the unloading pulse, the core of the shorter impactor bar produces additional fragments in the impactor. These phenomena can be seen in the case of the coarse, as well as fine mesh, in Fig. 14, Fig. 15, Fig. 16, Fig. 17h;17.

It is found that in both models, there is not excessive deformation and rotation of grains, in the failure zone, during the first 12 μs. Fragments displace radially outward giving the step pattern of the failure surface at the boundary which is consistent with the experimental observation as evident from the high speed photographs. The photographs show that the crack density gradually diminishes away from the impact surface. This behavior is predicted to some approximation by both models, though the deformation and rotation of grains in the simulations are appreciable. A larger gradient in crack density is observed in the continuum\discrete damage model.

An important observation to be made from these plots is that there are distinct planes at which radial bulging is observed. The length interval at which this takes place is the same in the case of the coarse and fine meshes indicating the phenomena to be related to a material behavior. One possible reason for this may be the attenuation of the tensile rarefaction radial waves due to release of energy in crack initiation and enlargement. Additional failure can take place only by the new rarefaction wave generated down the length of the bars. This observation is important from the experiment at viewpoint in deciding the location of the gauge for stress and radial velocity measurement. The constant length interval in the case of coarse and fine meshes indicates that the pulse duration and the attenuation of waves appear independent of the mesh size.

As mentioned earlier, one of the objectives of the present research is to develop a model which minimizes the run-time and yields results consistent with experimental observations. Fig. 18Fig. 19 compare the axial stress obtained from coarse and fine meshes with the discrete and continuum\discrete models, respectively. The friction coefficient in the case of continuum\discrete model is 0.2 while the discrete model uses zero friction coefficient. It is observed that there is not appreciable difference in the peak stress, pulse duration and stress attenuation predicted by both models for a given value of friction coefficient, between fragments, as the mesh is refined. Hence, the two models are predicted to be almost mesh independent for the two fragment sizes considered in the present analyses. This shows that the assumption of fragments containing a large number of grains is not far from realism in computational analyses.

It can be argued that it is surprising that not much mesh sensitivity is observed in the case of the pure discrete fragmentation model. In fact, it is expected that as the number of potential fragments in the mesh is increased, the stress history would approach the experimentally measured stress history. A possibility for the lack of significant increase in stress attenuation behind the wave front, with fragment size reduction, (see Fig. 18), may be due to the fact that an extremely fine mesh, with potential fragment sizes of the order of the material grain size is required to capture the initiation of all possible microcracks. Certainly, in the present analyses, the reduction in the fragment size in the case of fine mesh is far from achieving the actual grain size of the material.

A novel way of predicting structural response of ceramics during impact capable of capturing brittle failure features has been presented by combining two classes of damage models, viz., continuum and discrete damage models. The model so arrived, and named continuum\discrete model, is able to predict the response of ceramics in agreement with experimental observations. The concept of arriving at meaningful finite element simulations of such highly nonlinear class of problems, in a reasonable time, may be well served as established by comparing simulation results with experimental results. The simulations presented here establish the advantages of the continuum\discrete damage model over the discrete damage model in terms of the structural response of ceramics, e.g., deformation, stress history, fracture pattern, etc. The arguments in favor of the two classes of damage models may be applied simultaneously to the continuum\discrete model. However, one argument against the model may be related to the assumption of the initial fragment size and interface elements between them providing pre-determined crack paths. The discrete damage model with the capability of initiating as many cracks as required with their propagation in the correctly predicted direction, or the continuum\discrete damage model with the fragment size of the order of actual grain size may seem to be scientifically more accurate. However, such analyses would be so computer intensive that they may prove to be far from realism for field design applications. Also, such analyses would be carried out omitting useful experimental findings. On the contrary, potential fragment shape and size can be reasonably chosen based on experimental findings as discussed in the introduction. This further concretizes the concept of combined experimental\numerical investigation for such highly nonlinear and complex class of problems. The concept of choosing a random pattern of potential fragments, a priori, and embedding interface elements between them, i.e., prefixing the direction for crack propagation, may find application readily in many other problems for which the position of crack initiation and propagation are more or less well defined, e.g., machining of ceramics, delamination in fiber reinforced composite materials (Lu, 1998), etc.

The interface law plays a significant role in the continuum\discrete model. In principle, any interface law can be used after experimental determination of its parameters through well designed experiments, viz., plate impact experiments and Kolsky bar experiments (see contributions by Clifton and Klopp, 1985; Duffy, 1980). The rate effect included in our interface formulation is expected to eliminate the need for regularization through nonlocal analyses, thus making the prediction of dynamic failure mesh independent.

One disadvantage of the continuum\discrete model is that with the accumulation of damage, element distortion within fragments may become excessive. In such case, remeshing of elements within each fragment is required if erosion is not desirable. Such remeshing with a consistent mesh transfer operator can be easily implemented. Hence, it may be concluded that the present methodology provides a powerful tool for predicting the structural response of brittle materials in agreement with experimental measurements and observationsGurson, 1977, Zukas, 1990.

Section snippets

Acknowledgements

This research was supported by the National Science Foundation through grants MSS-9309006, MSS 9311289, CMS-9523113, the Army Office of Scientific Research through Grant No DAAH04-96-1-0142, and ARO-MURI Grant No. DAAH04-96-1-0331.

References (49)

  • G.T Camacho et al.

    International Journal of Solids and Structures

    (1996)
  • H.D Espinosa

    International Journal of Solids and Structures

    (1995)
  • H.D Espinosa et al.

    Journal of the Mechanics and Physics of Solids

    (1995)
  • H.D Espinosa et al.

    International Journal of Solids and Structures

    (1998)
  • G.R Johnson et al.

    Eng. Frac. Mech

    (1985)
  • B.R Lawn

    Acta. Met. et Mat

    (1994)
  • A Needleman

    Comput. Meth. App. Mech. Engg

    (1988)
  • D.L Orphal et al.

    International Journal of Engineering

    (1996)
  • D.L Orphal et al.

    International Journal of Impact Engineering

    (1997)
  • D.L Orphal et al.

    International Journal of Impact Engineering

    (1997)
  • D.A Shockey et al.

    International Journal of Impact Engineering

    (1990)
  • V Tvergaard

    Materials Science and Engineering

    (1990)
  • R.L Woodward et al.

    International Journal of Impact Engineering

    (1994)
  • F.L Addessio et al.

    LA-UR-

    (1989)
  • Andrews, E. W. and Kim, K. S. (1997) Threshold conditions for dynamic fragmentation of ceramic particles. International...
  • S.B Biner

    Proc. Symp. Materials Research Society

    (1994)
  • Brar, N. S., Yuan, G., Espinosa, H. D. and Zavattieri, P. D. (1997) Experimental study of interface defeat in confined...
  • Clifton, R. J. and Klopp, R. W. (1985) Pressure-shear plate impact testing, In MetalsHandbook : Mechanical Testing 8,...
  • D Curran et al.

    International Journal Impact Engineering

    (1990)
  • Duffy, J. (1980) The Dynamic Plastic Deformation of Metals : A Review, AFML-TR-82-4024, Wright-Patterson Air Force...
  • Espinosa, H. D. (1989) Finite element analysis of stress induced damage in ceramics. MSc thesis, Brown University,...
  • Espinosa, H. D. (1992) Micromechanics of the dynamic response of ceramics and ceramic composites. Ph.D. thesis, Brown...
  • H.D Espinosa et al.

    J. Hard Mat.

    (1992)
  • Espinosa, H. D., Emore, G. L. and Zavattieri, P. D. (1996) Computational modeling of geometric and material...
  • Cited by (0)

    View full text