Numerical approach to dimple fracture under very short pulse loading •

Dimple fracture of aluminum alloy 7075-T651 subjected to various stress intensity pulses with durations of 20, 40 and 80 μs is investigated by discrete dislocation dynamics and the Gurson-Tvergaard-Needleman constitutive model (GTN-model). A dominant void ahead of macro crack tip is modeled by a micro crack. The evolution of principal stress component along the inner ligament is computed by the sum of the crack tip stress field and the interaction stress field of dislocations emitted from both the micro crack and macro crack. The yielding flow of GTN-model is used to estimate void volume fraction, which is used for a fracture initiation criterion. The numerical results of the dimple fracture are obtained in a reasonable agreement with the experimental ones.


INTRODUCTION
Dimple fracture usually results from nucleation, growth, and coalescence of microscopic voids that initiate at inclusions and second-phase particles.In the review of ductile fracture [1], the initiation criterion of dimple fracture could be considered in terms of the critical crack tip opening displacement (CTOD) exceeding the mean void initiating particle spacing.A certain critical condition must be achieved not only at the crack tip, but also it must be achieved over a characteristic distance.The presence of smaller-scale voids in the ligament between the crack and the larger void was taken into account by using the GTN-model to describe void growth and coalescence process [2].The fracture initiation was assumed to occur once the failure condition was satisfied everywhere in the inner ligament.
The investigation of dimple fracture mechanisms of aluminum alloy 7075-T651 subjected to short pulse loading has been carried out successfully by Homma and Rizal [3][4][5].
Three stress pulse durations of 20, 40 and 80 μs are generated for three kinds of specimen configurations with the same cross section of 40 mm height and 10 mm thickness, and three different lengths (spans) of 50, 100 and 180 mm respectively.The fracture surface is examined by a scanning electron microscope (SEM) and the observation result shows that the fracture surface is entirely covered by dimples, and that a large void was initiated ahead of the crack tip prior to the small voids.The large void is considered to play a dominant role on the crack initiation and was fined as a dominant void.The crack tip is blunted further with the load increment and the nucleated dominant void grows.It is considered that the fracture initiation takes place when a Trang 46 dominant void and a crack tip coalesce with each other by fine voids nucleated as a void sheet.
The microscopic morphological features are identical or the same fracture mechanism takes place under all the stress pulse durations of 20, 40 and 80 μs.The dominant void is nucleated at the same distance from the crack tip for three pulse durations.Moreover, the dynamic fracture toughness remarkably increased for the pulse duration less than 40 μs and the corresponding CTOD also indicates the similar tendency.
The effect of impulsive loading duration on the plasticity ahead of crack tip has been successfully investigated by discrete dislocation dynamics (DDD) [6].This DDD model has been also extended to cleavage fracture under short pulse loading [7].This work aims to obtain the insight into time effect on the initiation of dimple fracture subjected to short pulse loading.The previous DDD model for plasticity is extended to a two-dimensional DDD coupling crack model, in which a micro crack is placed at the nucleation site of the dominant void.The void volume fraction in the inner ligament between a main crack tip and the micro crack tip is calculated by the GRN-model and postulated as a fracture initiation criterion.The numerical analysis is carried out to simulate the experiment done by one of the authors [5] , and the numerical and the experimental results are compared to examine the model validity.

Dislocation model
Because it is extremely complicated to describe a real dislocation movement in a polycrystalline material, a simplified DDD model, in which discrete dislocations are confined to the movement along a straight slip line in a homogeneous and isotropic material, was developed and it could successfully represent the loading time effect on the CTOD [6].The previous DDD model is extended to a twodimensional DDD coupling crack model to investigate the dimple fracture initiation.
The model consists of a main crack and an embedded finite micro crack located ahead of the main crack and along the crack plane.A micro crack is modeled to represent a dominant void ahead of a main crack.Therefore, the center of the micro crack is 100 μm distant from the main crack tip, which is equal to the experimental dominant void position [3][4][5].Both the cracks are assumed to be in an infinite elastic plane and in infinite thickness in the z-direction so that the plane strain condition prevails.The dislocation model is illustrated in Fig. 1.

Fig. 1 Schematic illustration of the dislocation model
Following the current specimen configuration, the main crack length 2a is 40 mm in x-direction.In the model, a slip system is considered at only one crack tip, because effects of dislocations on a slip system at the other crack tip could be negligibly small from estimation from simple calculation of stresses at a position 40 mm distant from a dislocation in an infinite plane as mentioned in the previous work [7].Therefore, the main crack tip plasticity is represented by only two arrays of discrete edge dislocations.The slip lines, on which the edge dislocations will be emitted, are symmetric and inclined against the crack plane at angles ±α.The Burgers vectors are parallel to the slip lines and correspondent to the unit vectors as: (cosα, sinα), (-cosα, sinα) for the inclined slip lines of the number 1 and 2 respectively.
The micro crack is used to model the presence of dominant void ahead of the crack tip.The length of the micro crack is taken as 20 μm, which is nearly equal to a half of the dominant void size observed in the experiment [3].It is natural to consider that the micro crack is located in the plastic zone ahead of the main crack tip.Therefore, four symmetric arrays of discrete edge dislocations are activated at both tips of the micrio crack.It is assumed that two slip lines exist at each crack tip and are inclined against the crack plane at angles ±α' and ±(π-α').The unit Burgers vector is discomposed as (cosα', sinα'), (-cosα', sinα'), (-cosα', -sinα') and (cosα', -sinα') for the inclined slip lines of the number 1 to 4, respectively.

2.2.
STRESS ANALYSIS ON A DISLOCATION AND DISLOCATION MOVEMENT

Main crack slip line
First, stresses applied to a dislocation on slip lines activated at the main crack tip will be considered in case where a micro crack does not exist.
The stress components are calculated by complex potential functions for the plane strain condition.Two complex potential functions φ(z) and ψ(z) of the variable z=x+iy are used as follows: where the superposed bar denotes the complex conjugate, , the prime (') and ('') denote the first and the second derivative of the complex potential functions.
A dislocation will contribute stress field surrounding it, namely dislocation stress field.Therefore, if a dislocation exists near a crack, it will contribute the traction stress and shear stress on the crack surface.Because of the zero traction stress and shear stress on the crack surface condition, V. Vitek et al. [8] composed the image complex stress potential functions from stress functions for a dislocation in an infinite plane to eliminate the traction stress and the shear stress on the crack surface.
To determine the stress field of the crack, including the applied stress σ a , the complex potentials are given as: where a is the half length of crack, and The stress field of an edge dislocation with Burgers vector (b x , b y ) located at (x d , y d ) in a smooth infinite plane is determined by: , G is shear modulus, and ν is Poisson's ratio.
In case of a cracked infinite plane, the image stress field can be determined by the complex stress functions as: The stress components for an edge dislocation on the inclined slip plane ahead of crack tip subjected to a uniform applied stress can be determined by summation of each stress component calculated from the equations (1) to (4).Then, the shear stress acting on an inclined slip plane rotated from the crack plane by an angle α can be found as: After n dislocations have been emitted, the resolved shear stress applied to the ith dislocation located at zi1 on an inclined slip line number 1, is obtained as: where z kj is the position of the k th dislocation on the inclined slip line of the number j.The equation consists of three parts: the shear stress caused by the crack subjected to the uniform applied stress σ a , the image shear stress caused by the dislocations at z i1 and z i2 on the slip lines of the number 1 and 2 respectively, and the summation of direct interaction shear stresses induced by all the other dislocations and their image stresses.

Dislocation emission and movement from the crack tip
The main crack is subjected to a semi-sinusoidal applied load function of time.To obtain converged and consistent results, the applied load-time history is divided into very short time intervals.
For the symmetry, simultaneous phenomenon occurs on two slip lines.At each time step, for example on the inclined slip line number 1, the resolved shear stress at the source position, τ res (z 01 ) is computed.A new dislocation will be emitted and allowed to move a away from the source when τ res (z 01 ) exceeds a critical value.The friction shear stress, τ f is assigned to be a critical shear stress for dislocation emission, following the same consideration as the other researchers [9,10].Moreover, the resolved shear stresses of all the other emitted dislocations are also computed and compared with τ f to determine the dislocation movement.Dislocation is feeling the resolved shear stress larger than τ f will move in the inclined slip direction with a velocity based on the following empirical expression [9,10]: The appropriate material parameters of A and m have been chosen to achieve a reasonable agreement between the numerical results and the experimental ones.
At the beginning of every time interval, resolved shear stresses are calculated and then the traveling velocities of all the emitted dislocations are calculated.Traveling distances of all the dislocations are evaluated to determine new positions of the dislocations at the end of every time interval.The new positions of dislocations are determined and updated for the next step through an interactive process by summation of the previous positions with the obtained traveling distances.
Concerning to the micro crack, it is approximately subjected to the corresponding applied stress, σ yy , which is assumed to be uniform around it and equal to sum of the stress fields from the tip of main crack and the emitted dislocations in the plastic zone [11].For each time step, the development of dislocation arrays around the main crack is obtained.Then the local tensile stress, σ yy at the position of x R =100 μm distant from the main crack tip to the center of micro crack is calculated by: The simulation of micro crack is definitely identical with the main crack.The resolved shear stress affects on an edge dislocation on an inclined slip line is also calculated by Eq. ( 6), but four slip lines (j=1÷4) are included instead of only two slip lines as in the main crack simulation case.Dislocation emission occurs simultaneously on four inclined slip lines at the sources, motion and positions of dislocations are also taken into account and updated as the above process.
The computing time is found to be significantly affected by the number of dislocations.The bundled technique is employed to reduce the computing time [12].A super dislocation of rank 2 is formed out or bundled of three dislocations of rank 1 after the ninth dislocation has been emitted.This process is repeated such whenever the ninth dislocation rank n∈N is created.

The GTN-model
Observation results of dimple fracture surfaces of aluminum alloy 7075-T651 [3,4,5] indicated that crack initiation is a result of coalescence between a dominant void and a crack tip.The coalescence is caused by nucleation of many fine voids in the inner ligament between the main crack and the dominant void, and their subsequent linkage.Therefore, to describe the dimple fracture, the nucleation, growth and coalescence process of fine voids in the ligament must be considered.
The modified yield condition form of GTNmodel for a void-containing or porous material can be expressed as [13]: is mean normal stress, σ M is an equivalent tensile flow stress representing the actual microscopic stress state in the matrix material, q 1 and q 2 are extra material parameters.
The void volume fraction will be used for an initiation criterion of dimple fracture in this model.The rate of void volume fraction in GTNmodel is given as: is the rate of mean plastic strain, A is the scaling coefficient, and M σ& is the rate of equivalent tensile flow stress.
The scaling coefficient A is given as: where f N , ε Ν   and  s N are a mean volume fraction of particles causing void nucleation, a mean value of strain distribution, and a standard deviation of the distributed strain, respectively.
To evaluate the evolution of void volume fraction f, the parameters, where the superscript dot denotes the material time derivative and E t is the tangential slope in the uni-axial true stress-natural plastic strain curve at the current stress level σ M .
The rate of equivalent plastic work in the matrix material is assumed to be equal to the macroscopic rate of plastic work and expressed as: where p ij ε& is the plastic part of the macroscopic strain rate indicated later.
Substituting the incremental work hardening relation ( 12) into ( 13) the rate of the equivalent flow stress in the matrix material is given as: Next, the rate of plastic volume strain p kk ε& must be calculated.According to the classical plasticity, the plastic flow rule is given as [2]: where Λ is the plastic flow factor, and determined from consistency condition.
The consistency condition requires that during plastic loading the yield condition must be satisfied as: The first term represents the change caused by the changes in macroscopic stress.The second represents the hardening effect in the matrix material, whereas the third term includes the contribution caused by changes in void volume fraction.
Substituting (10), ( 13) into ( 16) and solving for Λ in (15) finally we obtain: where The hardening rule of the matrix material in the uni-axial tensile test is assumed to be specified by the form [2]: where σ y is the tensile yield stress of the matrix material and N is the strain-hardening exponent.

Model parameters
In order to simulate, some numerical parameters are chosen by the previous experimental results of the aluminum alloy 7075-T651 [3][4][5].Young's modulus, E is 71 GPa and the Poisson ratio, ν is 0.33.The mechanical properties of 7075 aluminum alloy are substantially insensitive against the current strain rate [4].Therefore, the dynamic yield strength in the simulation is reasonably consistent and taken as 525 MPa for three pulse durations.
The friction shear stress τ f , which is provided by the material to resist movement of dislocations along the slip line, can be taken arbitrarily, but it is considered to be closely related to uni-axial yield strength.However, in the current simulation, the effect of time is considered and yield strength is insensitive against time.Therefore, τ f is taken consistent and equal to 50 MPa for three pulse durations.The magnitude of Burgers vector is chosen to be a common to DDD simulations and taken as b=0.25 nm.A dislocation is nucleated at a specific point near the crack tip on the slip line is a successful and realistic approach to modeling the dislocation nucleation problem [12].The position of dislocation emission source is determined as the nearest one that can provide stable numerical results and taken as 200b distant from the crack tip and is identical for both main crack and micro crack.
The inclined slip angles between the slip lines and the crack plane α and α' of both the main crack and micro crack are chosen to be identical and equal to 70.5 0 , which is the direction of the maximum shear stress of the crack tip field.The material parameters A and m in the dislocation velocity law in Eq. ( 7) are taken as 1.5x10 -4 and 2, respectively.Time interval ∆t is taken so small as 0.1 ns to prevent a dislocation from overshooting the others at each time step or to make sure that numerical stability is guaranteed.
The key input parameter to the GTN-model is the initial, f 0 .Certainly, depending on differences of microstructures, this value could vary in the range of zero to 0.06 [14].The initial void volume fraction, f 0 = 0.02 is used, so that the expected critical value f c is nearly equal to the average value of void volume fraction to use in the inner ligament for the failure criterion [2].
In order to use the GTN-model, the strain hardening behavior has to be determined.The strain-hardening exponent, N is taken as 0.03 for aluminum alloy 7075-T6 [15].The magnitude of σ M is assumed to be varied at a constant time step from σ y to the maximum value, so that the condition for the onset of necking is nearly satisfied [16].By using the relation between σ M and p M ε in Eq.( 22), the onset of necking condition is satisfied when σ M reaches the value of 1.045σ y .Following Aravas [2] a typical set of values for a material parameters in Eq. ( 9) are taken as q 1 =1.5, q 2 =1, The applied stresses to the main crack in a infinite plane are quasi-static uniform stresses that yield the identical time histories of the stress intensity obtained by the dynamic fracture toughness tests under three stress pulse durations [3][4][5].

RESULTS AND DISCUSSION
At first, the local tensile stress σ yy applied to the micro crack is calculated at the position 100 μm distant from the main crack tip, corresponding to the center of the micro crack, for three pulse durations by using Eq. ( 8).The applied uniform stresses σ a to the main crack are the quasi-static uniform stresses that yield sinusoidal stress intensity-time histories similar to ones obtained in the dynamic fracture toughness K Id tests for three stress pulse durations.The sinusoidal stress intensity-time histories have the amplitude of 24 MPam -1/2 , which is dynamic fracture toughness for 80 μs stress pulse duration.For each time step, dislocation arrays on the slip lines at the main crack tip are updated and then, σ yy applied to the micro crack is calculated as shown in Fig. 2.
Local tensile stress σ yy is composed of the tensile stress caused by the stress field of main crack, the interaction stress field of the emitted dislocations from it including the stress field of dislocations and the image stress field of dislocations.The micro crack is assumed to be subjected to the corresponding σ yy , which is approximated to be uniform around of it.Similar to the main crack simulation process, dislocations are emitted on the slip lines ahead of micro crack and the stress component ahead of it is also calculated by the distribution of crack tip stress field and the interaction of emitted dislocations stress field.
The principal stress component along the inner ligament between the main crack and micro crack is calculated by the contribution of the crack tip stress field and the interaction stress field of emitted dislocations from both of the cracks.For instances, the stress component as a function of time for three pulse durations is shown in Fig. 3 at the inner ligament position of x L =50 μm distant from the main crack tip.The maximum value of stress component along the inner ligament distance for three pulse durations is shown in Fig. 4. Stress component tends to increase at both crack tips caused by the stress field of cracks.For the identical magnitude of the applied load, under the longer pulse duration, the smaller value of stress component is obtained.This could be explained that the longer loading time is, the more dislocations are emitted so that the macroscopic stress component is affected by the surrounding emitted dislocations.
In generally, stress component is gradually increased as a function of time for the whole of pulse duration as in Fig. 3.However, for a very small step of time it can be varied or does not follow the trend of increasing caused by the computational calculation errors.To calculate the rate of stress component, it is necessary to be approximated.For simplification, the stress component is assumed to be approximated as a semi-sinusoidal function of time.The semisinusoidal stress component-time profile is divided into small intervals.The rate of stress component as a function of time for three pulse durations is shown in Fig. 5 at the inner ligament position of x L =50 μm.The shorter pulse duration is, the higher stress component rate is obtained.
10) must be determined.The equivalent tensile flow stress σ M is related to the effective plastic strain state p M ε by the incremental relation:

Fig. 2 Fig. 3
Fig. 2 Local tensile stress contributed by main crack as a function of time

Fig. 6
Fig. 6 Void volume fraction along the inner ligament distance