Automated dynamic fracture procedure for explicit time integration brick finite elements

An elastodynamic fracture model has been implemented in the explicit finite element software DYNA3D that predicts energy release rates of stationary, through-thickness, 3-D cracks in linearly elastic materials. This work is part of an on-going effort to investigate implementation of automated fracture models in DYNA3D. It follows the implementation of a linear elastic fracture model that is capable of simulating automatic crack propagation without user intervention. The current model uses a pathindependent volume integral expression obtained by modifying an expression developed earlier for 2-D crack problems to compute the dynamic energy release rate. It is implemented for 3-D solid (brick) elements. Domain integral method is used to develop the volume integral expression. Domain integral form of the expression is particularly wellsuited for applications with the finite element method as it overcomes the difficulty associated with defining contours around the crack tip. Also, it does not involve elements around the crack front, thereby leading to better accuracy when using the finite element method for crack analysis. The implementation of the model has two basic steps search for elements in a chosen integration volume, and numerical evaluation of the integral expression. The integration volume to be used is input by means of two values one for number of rings of elements around the crack front to be ignored, and the other for outer limit of the volume. Some mechanical field quantities in the integral expression are not available in DYNA3D’s brick element implementation. These values are determined for the integration volume elements and stored as additional history variables, if needed, during the numerical evaluation phase. Numerical examples to verify the accuracy of the current model are presented. INTRODUCTION Advancements in non-destructive examination techniques have made it increasingly obvious that all structures contain cracks or crack-like flaws. Hence, fracture mechanics has become a significant part of the structural design process and thereby an active subject of research. Int. Jnl. of Multiphysics Volume 1 · Number 1 · 2007 33 * Associate Professor, author to whom correspondence should be addressed (ala.tabiei@uc.edu) † Graduate Research Assistant Linear elastic fracture mechanics (LEFM) is the branch of fracture mechanics that deals with behavior of bodies of nominally elastic materials containing cracks subjected to quasistatic loading conditions. LEFM ceases to be valid when significant plastic deformation precedes failure or time becomes an important variable or both. For example, inertia effects are important when the load changes abruptly or the crack grows rapidly and LEFM has been found to be inadequate for such problems. As a result, significant advances have been made in incorporating other types of material behavior and dynamic analysis in fracture mechanics. Many experimental approaches have been devised to study the process of dynamic fracture. However, due to short spans of the fracture events, it is difficult to directly measure higher-order physical quantities like energy distribution, instantaneous energy release rate, and dynamic stress field close to the crack tip. This drawback can be overcome by using computational approaches. Therefore, the advancement of dynamic fracture mechanics relies heavily on parallel advances in computational methods for dynamic fracture simulation [1]. A detailed review of the impact of computational technology on furthering the understanding of fundamental fracture phenomena can be found in [2]. Among the widely used computational methods, the finite element method is well established as a powerful and versatile numerical technique for solving solid mechanics problems. Together with the latest computing technology, it is very well suited for advanced fracture mechanics problems. However, in the analysis of dynamic crack phenomena by means of computational methods, such as the finite element method, a fundamental difficulty is encountered in efforts to compute values of crack tip energy flux versus time or amount of crack growth. The difficulty arises from the fact that, on the one hand, the crack tip energy flux is defined in terms of values of field quantities for points arbitrarily close to the crack tip while, on the other hand, it is precisely for points near the crack tip for which the accurate calculation of field quantities is most difficult [3]. A detailed review of the major advances in finite element theory and implementation which addresses the problems of fracture mechanics can be found in [4]. The first crack tip contour integral expression for elastodynamic energy release rate was proposed by Atkinson and Eshelby [5], who argued that the form for dynamic growth should be the same as the quasi-static growth with the elastic energy density replaced by the total mechanical energy density, that is, the elastic energy plus the kinetic energy [3]. Kostrov and Nikitin [6] and Freund [7] subsequently derived the equivalent integral expression for dynamic energy release rate in terms of crack tip stress and deformation fields directly from the field equations of elastodynamics. Nakamura et al. [3] developed an integral expression for crack tip energy flux in terms of near tip mechanical fields which is applicable to a broad range of material response. More significantly, they showed that many seemingly different path-independent integrals that were developed earlier can be extracted from their expression by invoking appropriate restrictions on material response and crack tip motion. They also gave a couple of modified forms of their expression that are more suitable for analysis of dynamic crack phenomena by means of computational methods. One such expression is a domain integral which circumvents the difficulty associated with the accurate calculation of field quantities for points arbitrarily close to the crack tip. Li et al. [8], and Moran and Shih [9] also developed domain integrals by converting line integrals to surface integrals in 2-D problems and surface integrals to volume integrals in 3-D problems, which then enabled the use of stresses and strains at the element integration points. Charoenphan et al. [10] applied a similar procedure to the 2-D expression developed by Nakamura et al. [3] and developed a domain integral expression for the dynamic energy release rate of through-thickness 3-D 34 Automated dynamic fracture procedure for explicit time integration brick finite elements


INTRODUCTION
Advancements in non-destructive examination techniques have made it increasingly obvious that all structures contain cracks or crack-like flaws.Hence, fracture mechanics has become a significant part of the structural design process and thereby an active subject of research.
Linear elastic fracture mechanics (LEFM) is the branch of fracture mechanics that deals with behavior of bodies of nominally elastic materials containing cracks subjected to quasistatic loading conditions.LEFM ceases to be valid when significant plastic deformation precedes failure or time becomes an important variable or both.For example, inertia effects are important when the load changes abruptly or the crack grows rapidly and LEFM has been found to be inadequate for such problems.As a result, significant advances have been made in incorporating other types of material behavior and dynamic analysis in fracture mechanics.
Many experimental approaches have been devised to study the process of dynamic fracture.However, due to short spans of the fracture events, it is difficult to directly measure higher-order physical quantities like energy distribution, instantaneous energy release rate, and dynamic stress field close to the crack tip.This drawback can be overcome by using computational approaches.Therefore, the advancement of dynamic fracture mechanics relies heavily on parallel advances in computational methods for dynamic fracture simulation [1].A detailed review of the impact of computational technology on furthering the understanding of fundamental fracture phenomena can be found in [2].
Among the widely used computational methods, the finite element method is well established as a powerful and versatile numerical technique for solving solid mechanics problems.Together with the latest computing technology, it is very well suited for advanced fracture mechanics problems.However, in the analysis of dynamic crack phenomena by means of computational methods, such as the finite element method, a fundamental difficulty is encountered in efforts to compute values of crack tip energy flux versus time or amount of crack growth.The difficulty arises from the fact that, on the one hand, the crack tip energy flux is defined in terms of values of field quantities for points arbitrarily close to the crack tip while, on the other hand, it is precisely for points near the crack tip for which the accurate calculation of field quantities is most difficult [3].A detailed review of the major advances in finite element theory and implementation which addresses the problems of fracture mechanics can be found in [4].
The first crack tip contour integral expression for elastodynamic energy release rate was proposed by Atkinson and Eshelby [5], who argued that the form for dynamic growth should be the same as the quasi-static growth with the elastic energy density replaced by the total mechanical energy density, that is, the elastic energy plus the kinetic energy [3].Kostrov and Nikitin [6] and Freund [7] subsequently derived the equivalent integral expression for dynamic energy release rate in terms of crack tip stress and deformation fields directly from the field equations of elastodynamics.Nakamura et al. [3] developed an integral expression for crack tip energy flux in terms of near tip mechanical fields which is applicable to a broad range of material response.More significantly, they showed that many seemingly different path-independent integrals that were developed earlier can be extracted from their expression by invoking appropriate restrictions on material response and crack tip motion.They also gave a couple of modified forms of their expression that are more suitable for analysis of dynamic crack phenomena by means of computational methods.One such expression is a domain integral which circumvents the difficulty associated with the accurate calculation of field quantities for points arbitrarily close to the crack tip.Li et al. [8], and Moran and Shih [9] also developed domain integrals by converting line integrals to surface integrals in 2-D problems and surface integrals to volume integrals in 3-D problems, which then enabled the use of stresses and strains at the element integration points.Charoenphan et al. [10] applied a similar procedure to the 2-D expression developed by Nakamura et al. [3] and developed a domain integral expression for the dynamic energy release rate of through-thickness 3-D cracks which they then implemented for shell element in the explicit finite element software DYNA3D for modeling crack growth in composite material shell structures.
The primary objective of the present work is to implement the domain integral expression developed by Charoenphan et al. [10] for 3-D solid (brick) elements in DYNA3D.This work is part of an on-going effort to investigate the implementation of fracture models in DYNA3D to simulate automatic crack propagation using brick elements.DYNA3D is an explicit finite element software used for analyzing the dynamic response of 3-D solids and structures.It uses element removal techniques to simulate the failure process.The crack opening profile therefore cannot be modeled.Furthermore, the stress based failure methodology is not able to describe failure accurately.Therefore, it is important to implement fracture theories for solving crack propagation problems [11].The present work follows the implementation of LEFM models by Tabiei and Wu [11].Their models have the capabilities to simulate automatic crack propagation without user intervention.Details of the implementation of the current model and verification examples to validate it are given in this report.

INTEGRAL EXPRESSION FOR DYNAMIC ENERGY RELEASE RATE IN LINEARLY ELASTIC MATERIALS
The energy approach to fracture analysis states that crack extension occurs when the energy available for crack growth is sufficient to overcome the resistance of the material.Griffith [12] was the first to propose the energy criterion for fracture, but Irwin [13] is primarily responsible for developing the present version of this approach: the energy release rate, G , which is defined as the rate of change in potential energy with crack area for a linear elastic material.At the moment of fracture, G = G C , the critical energy release rate, which is a measure of fracture toughness. ( where ͟ is potential energy and A is crack area [14]. When time becomes an important variable in fracture problems, inertia effects become significant and dynamic fracture analysis is required for practical and reasonably accurate solutions.For the application of energy approach to dynamic fracture problems, Nakamura et al. [3] developed the following 2-D integral expression for dynamic energy release rate of a plane crack in a linear elastic body advancing in the x 1 (1) direction (normal to the crack tip): (2 where U is the strain energy density, T is the kinetic energy density, σ ij are the components of stress, u i are the components of displacements, ρ is the mass density of the material, n j are the components of the unit normal vector to a contour Γ 0 around the crack tip which begins on one crack face and ends on the opposite face, and A 0 is the area enclosed by Γ 0 , another contour Γ around the crack tip but closer, and a portion of crack faces between the end points of the contours.The contours and area are shown in figure (1).This expression is valid for remote contour paths in the absence of body forces in linearly elastic materials.Note that a comma in the sub-script of a term, in the above expression as well as others in this report, denotes a partial derivative with respect to the spatial direction represented by the component that follows, i.e., .A unique feature of this expression is the fact that it does not involve mechanical field quantities close to the crack tip which makes it particularly well suited for implementation in the finite element method.However, it is recommended only for problems in which the ratio v/c 2 is less than 0.5 where v is the instantaneous crack tip speed and c 2 is shear wave speed in the material.Nakamura et al. [3] attributed this condition to the second and fifth terms in the expression.These terms are typically comparable in magnitude but of opposite signs leading to their net contribution to the energy release rate being small and a potential source of computational error.At lower speeds their contribution to the energy release rate is insignificant relative to the other terms, however, the contribution becomes significant at higher crack tip speeds.
In 3-D crack problems, an expression similar to (2) for the dynamic energy release rate can be written as: (3) where contour surface S 1 is a through-thickness cross-sectional area encircling the crack front, B is thickness of the structure, δ ij is kronecker delta, m j are the components of the unit normal vector of S 1 , V 1 is the volume enclosed by S 1 excluding the area in the vicinity of the crack front, and all other terms are same as before.It has to be noted that the above expression corresponds to a through-thickness crack in a uniform thickness structure.
There is one major issue though, with the above expression related to applications with the finite element method.Ideally, for a finite element model using this expression, the surface S 1 should be defined such that it passes through the element integration points as stresses and strains are determined most accurately at these points in finite element analyses.Though this is certainly possible, it leads to other inconveniences.A second option is to pass it through element boundaries and determine the corresponding values by interpolation and/or extrapolation from nodal values.This method however, leads to a loss in accuracy.Another option is to convert the surface integral into a volume integral using the divergence theorem but this is also not viable as the volume thus obtained would involve elements close to the crack front which would again induce inaccuracies.
This difficulty is circumvented by using finite domain integrals which were introduced by Kishimoto, Aoki and Sata [15,16], Atluri [17], and Nishioka and Atluri [18,19] that uses the concept of virtual crack extension methods.The procedure that converts equation ( 3) into a domain integral form more suitable for finite element analysis of 3-D crack problems is summarized in the next section.

INTEGRAL EXPRESSION FOR DYNAMIC ENERGY RELEASE RATE IN LINEARLY ELASTIC MATERIALS USING DOMAIN INTEGRAL METHOD
The dynamic energy release rate expression in equation ( 3) is converted to a domain integral form as it is more convenient for use with the finite element method.The domain integral form is obtained by converting the surface integral in (3) into a volume integral involving elements on the outer boundary of volume V 1 only as follows.
Consider a uniform thickness 3-D plate with a through-thickness crack, as shown in figure (2).Surfaces S 1 , S 2 , S 3 , and S 4 , volume V 1 , direction x 1 , assumed as the direction in which the crack is advancing, and unit normal vectors, m, pointing outward from the surfaces, are also shown.

Figure 2 Integration path for the domain integral evaluation
Consider now the surface formed by S 1 U S 2 U S 3 U S 4 .Obviously, the surface integral of equation ( 3) evaluated over S 1 alone will differ from the same evaluated over S 1 U S 2 U S 3 U S 4 by the contributions of the additional surfaces S 2 , S 3 , and S 4 .If these contributions are eliminated, then the integral value along S 1 and the combined surface will be the same.This can be achieved as summarized below.Though it appears that the contributions of three surfaces have to be eliminated, actually, it is only left to eliminate the contribution of S 3 as the integrals on surfaces S 2 and S 4 vanish since on these surfaces σ ij m j = 0 and m 1 = 0. Defining components m j = q 1 n j , equation ( 3) can be written as: (4) Now, if q 1 is chosen as a scalar function that is unit value on S 1 and zero on S 3 , the surface integral on S 3 will be zero and its contribution is eliminated.Therefore, the expression for dynamic energy release rate is transformed to: (5) where S is now the surface formed by S 1 U S 2 U S 3 U S 4 and enclosing the volume V, q 1 is a scalar function defined as: and all the other terms are same as before.
Finally, applying divergence theorem to equation ( 5), the energy release rate is written in its domain integral form as:

6)
where V B is the volume of elements enclosed by S 1 minus the volume of elements enclosed by S 3 , and V A is the volume of elements on the outer boundary of V B with S 1 as their outer surface.Figure (3) shows a simple example of integration volumes V A and V B that can be used to determine the dynamic energy release rates for a semi-infinite crack using equation (6).

IMPLEMENTATION OF THE MODEL IN DYNA3D
As mentioned earlier, the current model is implemented for the eight-node brick element in DYNA3D.This element uses one-point reduced integration by default.For this element, equation ( 6) can be discretized as: (7) where J is the element Jacobian relating the element volume in local coordinates to the volume in global coordinates and all other terms are as given before.The factor of 8 multiplying the terms originally in the integral expression is the weight for the single integration point, (0,0,0) in local element coordinates.
For a linear elastic material, strain energy density U and kinetic energy density T can be written as: Substituting equations ( 8) and ( 9) in equation ( 7) and expanding the index notation leads to the equation: (10) Note that the vector and tensor quantities in the above expression are in a local coordinate system and not in the global coordinate system.The local coordinate system is a Cartesian system with the axes along the directions normal (1-direction), bi-normal (2-direction) and tangent (3-direction) to the crack front as shown in figure (4).It would seem like the evaluation of the energy release rate expression in equation ( 10) is a straight forward process.However, this is not the case as only two of the quantities needed to evaluate the expression are readily available for the brick element as implemented in DYNA3D, namely the strain energy density and the stresses.All the other quantities are determined at each time step and stored as history variables, if needed, when using the current model and then used.A term-by-term description of the determination of the quantities is given later in this section.
The current model uses the same methodology to model actual 3-D crack shapes in finite element models as the earlier LEFM model, [11].The crack shapes are modeled with edges and surfaces of brick elements.The cracks are decomposed into a series of "sub-cracks" at the crack front each of which is defined by three nodes with the identification numbers 1, 2 and 3, with node 2 on the crack front, as shown in figure (5).The input phase subroutines of DYNA3D have been modified to read the crack input parameters and allocate memory for the fracture model arrays.Implementation of the current model has two major steps: determination of elements in the integration volumes V A and V B , and evaluation of the integral expression after determining the required quantities of the integration volume elements.
The first stage of the implementation is the determination of elements in V A and V B .They are determined based on two input values, the number of rings of elements around the crack front that are excluded from the evaluation (elements enclosed by surface S 3 ), henceforth denoted by nofrings1, and the ring number of the outer most layer of elements in the chosen volume, henceforth denoted by nofrings2.The input phase subroutines are modified to read nofrings1 and nofrings2 along with the thickness of the structure as additional inputs from the input file.New subroutines have been created to find the elements in the integration volumes, and nodes on the outer surface S 1 integration path based on nofrings1, nofrings2, and nodes on the crack front.At the first time step, these elements and nodes are found and stored in separate arrays.
Starting from the second time step, for each crack in the structure, the required field quantities are determined for the elements in the integration volumes at their integration points, (0,0,0) in local element coordinates, and dynamic energy release rate is computed using equation (10).As mentioned earlier, some of the required field quantities are readily available in DYNA3D while others are determined in the implementation.A brief description of the procedure used to determine the various quantities follows.
Density ρ : Density of the material is available from the finite element model.
Strain energy density U: Strain energy density of all elements are stored in the main data array of DYNA3D.Values for the integration volume elements are retrieved from this array for the computation.
Kinetic energy density T: In each element, components of velocity of the integration point, v ip α , in global coordinates, are first determined by interpolating nodal velocities using the shape functions as follows: (11) where v i α are the components of nodal velocities, and φ j are the shape functions at the integration point.The nodal velocities in global coordinates are also not available at the current time step, only the values lagging by half a time step are available.Therefore, before using the above relation, the current values are obtained using equation (14) given later in this section.
Components of the integration point velocity in local crack coordinates, , are obtained using the relation: (12) where T jα are elements of the transformation tensor that relates unit vectors of the global and local crack coordinate systems.
Finally, kinetic energy density of the element is determined using equation ( 9) in the previous section.
Stresses σ ij : Stresses in the elements, in global coordinates, are available in the main array.Each integration volume element's values are retrieved and transformed to local crack coordinates using the following coordinate transformation relation for tensors: (13) where , are elements of the stress tensor in local crack and global coordinates respectively, and T iα , T iβ are again elements of the transformation tensor.

Velocity gradients
. u i,1 : Velocity gradients at the integration points of elements are not readily available.In each element, they are determined from nodal velocities in global coordinates lagging by half a time step, which are the only related values available in the main array.First, the nodal velocities at the current time step are obtained using the relation: (14) where , are components of nodal velocities at current time step and half a step behind respectively, a n α are components of nodal accelerations at current time step, and ⌬t is the current time step.Then, elements of the velocity gradient tensor at the integration points, in global coordinates, v ip α,β , are obtained using the relation: (15) where v j α are the components of nodal velocities in global coordinates, and φ j,β are the elements of the strain displacement matrix at the integration point.Finally, the required velocity gradient at the integration point, in local crack coordinates, .u i,1 , is obtained using the relation: (16) where T iα and T 1β are again components of the transformation tensor.
Displacement gradients u i,1 : Knowing the velocity gradients, increments of the displacement gradients, in local crack coordinates, ∆u ip i,1 , are obtained using the relation: Then, the total displacement gradient at the current time step is obtained using the relation: (18) where n u ip i,1 , n-1 u ip i,1 , are the displacement gradients at the current time step and previous time step respectively.
As obvious from equation ( 18), these values are stored as additional history variables for the elements in the integration volumes.
Accelerations ü : Accelerations at the integration points of elements are determined from nodal accelerations which are in-turn obtained from resultant nodal forces.The nodal accelerations at time step n are determined using the relation: (19) where j ü n α are components of acceleration of node j, m j is mass of node j, ( j F n α ) ext , ( j F n α ) hg , ( j F n α ) int are components of external, hourglass, and internal nodal forces at node j, C is a damping co-efficient, and j v n-1/2 α are the components of velocity of node j lagging by half a time step.
Then, components of acceleration of the integration point of an element in global coordinates, ü ip α , is obtained using the relation: (20) where ü j α are components of nodal accelerations, and φ j are the shape functions at the integration point.
Finally, acceleration components in local crack coordinates, ü ip j , are then obtained using the relation: (21) where T jα are again elements of the transformation tensor.
Jacobian J: Jacobian for the brick elements is equal to one-eight their volume per their definition.
Gradients of q: Gradients of the q function with respect to the local crack axes, q 1,1 , q 1,2 , q 1,3 , are obtained by introducing a function in the global coordinates as: where ê 1 is the unit vector in the local crack co-ordinate system's x 1 (1) direction.It is assumed to be a nodal property, thereby an 8×3 array for each element, in the implementation.The value of at the integration point of an element is then obtained by interpolation, using the element shape functions as: (22) where are the components of at the integration point, are the corresponding components at the nodes, and φ j are the element shape functions.
Gradients of the components at the integration point, , are then obtained using the relation: (23) where φ j,β are the elements of the strain displacement matrix.
Finally, the required gradient values of in crack coordinates are obtained using transformation relations for tensors as follows: (24 where q 1,i are the required gradients of the components of at the integration point in crack coordinates, and T 1α , T iβ are again elements of the transformation tensor.

MODEL VERIFICATION
Four stationary crack problems are used to validate the current model.Though the domain integral expression used in this model is capable of predicting energy release rates for cracks propagating at speeds up to one-half the shear wave speed in the material, currently, it is not suited for simulations of propagating cracks.The automatic remeshing strategy implemented by Tabiei and Wu [11] to model the crack growth explicitly uses the strategy of delete-andfill process: first, a group of elements in a region around crack front is deleted, then the crack is extended and finally, the local domain is refilled with new elements.Since this process involves the instantaneous release of some nodal constraints, spurious high-frequency oscillations are observed in the finite element solutions making the model unsuitable for simulations of crack propagation problems.
All geometries in the following examples are meshed with brick elements as the current model has been implemented for use with brick elements only.A couple of quasi-static problems used to study the effect of mesh density and prove the path-independence of the finite domain integrals are presented first followed by a couple of dynamic problems.

ENERGY RELEASE RATE FOR A SEMI-INFINITE CRACK UNDER QUASI-STATIC CONDITIONS
Long duration simulations in explicit dynamic codes are an acceptable way of simulating quasi-static problems.Hence, for long duration simulations, the dynamic energy release rate predicted by the current model should converge to the quasi-static solutions.This condition is used in this example and the next for validation of the current model and to study the effect of mesh density.
In this example, a mode-I crack in a semi-infinite plate, as shown in figure (6), whose edges at finite length are displacing at a quasi-static rate under plane-strain conditions, is considered.Freund [20] determined the energy release rate in this case to be: where G is energy release rate, u 0 is displacement of the finite edges, H is half-width of the plate, E is Young's modulus of the material, and v is Poisson's ratio of the material.
Three finite element models of a plate of length 2540 mm, width 508 mm, thickness 6.35 mm, Young's modulus 68.9 GPa, and Poisson's ratio 0.3, shown in figure (6), are created with 500, 2000, and 8000 elements.Displacement at each gripped edge is prescribed as a linearly increasing function of time.For a displacement of 0.635 mm, the analytical solution is found from equation (25) to be 120.3KJ/m 2 .The simulation time for the finite element analysis is chosen to be 4 milliseconds which is long enough for the response of this problem to be quasi-static.The energy release rates predicted by the current model are given in Table (1).
As observed from values in the table, the energy release rates predicted by the current dynamic model are very close to the analytical solution, with the maximum error being 1.91% for the 200×40×1 (finest) mesh with nofrings1 = 5 and nofrings2 = 15 and 18, the two largest domains.This is a reasonable error considering the fact that there is an accumulation of errors associated with the various field quantities in the integral expression with the increase in the number of elements in the integration volumes.Also, since it is a fine mesh and nofrings1 = 5 , more elements close to the crack front are used in the computation, which as mentioned is a major source of error in fracture analysis using the finite element method.Another point to be noted from the values is that there is very little variation in the results predicted by the current model with change in integration volumes.The standard deviations are found to be 0.501076, 0.272261, and 0.666184 for the meshes with 500, 2000, and 8000 elements respectively.This minimal variation proves the path independence or in this case domain independence of the current model, which is a significant feature of integral expressions for the energy release rate.

ENERGY RELEASE RATE FOR A CENTRAL HORIZONTAL CRACK IN A SQUARE PLATE UNDER QUASI-STATIC CONDITIONS
The model for this example is a square plate with a horizontal central crack subjected to remotely applied uniform traction σ = 100MPa .The geometry, material properties and boundary conditions are shown in figure (7).The exact value of stress intensity factor K 1 for this case is 4.72 MPa m 1/2 (Tada [21]) and the corresponding energy release rate is 101.366J/m 2 for plane strain.The traction is applied as a ramp load first and then kept constant after 1 sec.for a total of 2.0 sec.simulation, as shown in figure (8).Four different meshes with 1280, 2880, 5120 and 8000 elements are used to study the mesh sensitivity.Figure (9) shows the energy release rate curves obtained for these meshes.As expected of finite element analyses, the values converge close to the exact solution with mesh refinement, the 8000 element mesh converging to a value of approximately 99.26 J/m 2 , an error of 2.07%.This shows that mesh refinement will provide convergence of the energy release rate.Figure (10) shows results of the test for path-independence, the curves obtained for different sets of values for nofrings1 and nofrings2, that is, different integration volumes, for the 8000 elements model.The different curves are denoted by two numbers with the first one denoting nofrings1 and the second one nofrings2.Better than the previous example, no difference is seen in the energy release rate values obtained for all the different cases considered proving the path-independence of the model.

ENERGY RELEASE RATE FOR A SEMI-INFINITE CRACK SUBJECTED TO A SINGLE STEP PULSE
In this example again, the case of a semi-infinite crack with an exact solution is used for validation.Freund [22] obtained the following analytical solution for a semi-infinite crack in an infinite plane subjected to a tensile stress wave normal to the crack plane.At t = 0, the stress wave arrives at the crack plane and at a later time t = τ , the crack starts to propagate.The analytical solution for the dynamic stress intensity factor in the interval 0 < t < τ (i.e. the crack is stationary) is: where c 1 is the longitudinal wave speed in the medium, σ a is the magnitude of the tensile pulse with step function time dependence, and B is a function of Poisson's ratio given in [22].
For plane strain, the dynamic energy release rate is then related to the dynamic stress intensity factor as follows: (27 where v is the crack tip speed, and A(v) is a universal function of crack tip speed.For a stationary crack (v = 0), A is unity.
To simulate the infinite plane problem for which equation ( 26) is applicable, the half-width of the model, H , must be large enough so that within the period of interest, reflected waves from the boundaries do not reach the vicinity of the crack tip.In addition, the length of the model (in x 1 direction) must be sufficiently large compared to H so that within the interval of interest, the crack tip region is unaffected by the unloading waves emanating from the remote edges.As the geometry of the model used in example (1), shown in figure (6), is large enough to satisfy these geometric constraints, it is again used for this validation example also with the same mesh densities as given earlier.
A step pulse is applied at t = -H/c 1 at the edge Simulation time in DYNA3D input file = 0.085 ms The variation of dynamic energy release rate with time for this problem, obtained using the three finite element models with 8000, 2000, and 500 elements are shown in figure (11) along with the exact solution.Two major observations can be made from the curves.First, overall best results for this stress wave problem are obtained for the model with the finest mesh i.e. 8000 elements.Second, considering this finite element model's results only, it can be observed that errors are very low, less than 5%, after t > 0.4 H/c 1 and substantially larger at shorter times.
Both these variations are not unexpected from the current model for this stress wave problem.Since the energy release rate is evaluated in this model using an integral expression that involves mechanical field quantities of a subset of elements of the finite element model, the predicted results depend to a great extent on the number of elements of the integration domain through which the refracted wave from the crack-tip has traveled -more the number of elements better are the results.Comparing coarse meshes to finer ones, clearly the number of elements that the refracted wave has reached is lesser in the case of the former due to the relatively larger element sizes.Therefore, the finest mesh yields the best results for this problem.Considering a particular mesh only, at lower times the wave is in elements closer to the crack-tip and some of these elements, as discussed earlier, are not included in the computation of the energy release rate.This leads to the inclusion of only few elements through which the wave has passed in the computation of the energy release rate during these times leading to higher error in the predicted values.

ENERGY RELEASE RATE OF A CENTRAL CRACK IN A FINITE PLATE SUBJECTED TO DYNAMIC LOADING
This example is another test for the current model under dynamic conditions.A rectangular plate with a central crack under plane strain conditions is examined.The plate, shown in figure (12), is loaded dynamically in the axial direction by a uniform tension P(t) with Heaviside-function time dependence as shown in figure (13).The material is linearly elastic with shear modulus G = 76.923GPa, Poisson's ratio v, and mass density ρ = 5000 Kg/m 3 .The thickness of the plate is assumed to be 0.4 mm and the crack length is 2a = 4.8 mm.A finite element model comprising of 5000 (100×50×1 ) elements and 10324 nodes is used for the simulation.This problem was solved numerically by Chen [23] using finite differences and has been frequently used as a reference to validate other methods.
Variation of a normalized energy release rate G/G ref where predicted by the current model and those obtained from Chen's [23] results are shown in figure (14).Since both are numerical results, no comments can be made about the accuracy of the current model based on these reference values.The energy release rates predicted by the current model are found to be lower than the reference values in general, especially close to the peak value.Overall, good agreement is observed between the values validating the current implementation.

CONCLUSIONS
A dynamic fracture model has been successfully implemented in the explicit finite element software DYNA3D.This work is part of an on-going effort to investigate the implementation of automated fracture models in DYNA3D for simulation of crack propagation problems using 8-node solid (brick) elements.The current model computes the energy release rates of through-thickness 3-D cracks in linear elastic bodies using brick elements.It uses a domain integral expression that is well suited for use with the finite element method.The implementation comprises of two basic parts -search of elements in an user-input integration domain, and numerical evaluation of the energy release rate expression.Details of the implementation are provided in this report.Numerical examples of stationary cracks under different loading conditions are presented to validate the model and prove its accuracy.Only stationary crack problems are presented since this model is not suitable for analysis of crack propagation problems at the moment even though the energy release rate expression used in the model is applicable to crack propagation problems with crack-tip speeds up to one-half the shear wave speed in the material.When crack propagation problems are simulated, the energy release rate begins to oscillate at a very high frequency from the first instant when a crack's length is incremented.This oscillatory behavior is attributed to instantaneous release of nodal constraints during explicit modeling of the crack growth and can be overcome by using a special scheme to release the nodes gradually.Implementation of such a gradual nodal release technique is currently in progress and details of this work will be presented in a following report.The work presented in this report is only the preliminary stage in the author's effort to implement a dynamic fracture model in DYNA3D that can be used in a wide range of applications.In addition to the gradual nodal release technique, further capabilities like the ability to simulate automatic crack propagation in different types of materials, non-self-similar propagation, etc. will be added to the current model and reported in the future.

Figure 3
Figure 3 Example of a finite domain for the domain integral expression for dynamic energy release rate

Figure 4
Figure 4 Local coordinates at a 3-D crack front

Figure 5
Figure 5 Decomposition of a 3-D crack tip

Figure 6
Figure 6 Mode -I crack in a semi-infinite plate

Figure 7 Figure 8
Figure 7 Geometry, material properties, and boundary conditions of central horizontal crack problem

Figure 9
Figure 9 Central horizontal crack in a square plate: G Vs t for different meshes

Figure 10
Figure 10 Central horizontal crack in a square plate: test for path independence

x 2 =
H of the model, and the pulse arrives at the crack plane ( x 2 = 0 plane) at t = 0 .Other input data are as follows: • Magnitude of the tensile pulse σ a = 1GPa • Longitudinal wave speed in material (Aluminum alloy) c 1 = 5787 mm/ms , • If pulse is applied at edge x 2 = H = 254 mm , then time taken for longitudinal wave to reach the other edge (x 2 = -H ) and reflect = 2H/c 1 = 0.08778 ms • Time taken for stress wave to reach crack tip = H/c 1 = 0.04389 ms •

Figure 11
Figure 11 Variation of normalized energy release rate with normalized time for a semi-infinite crack in an infinite plane subject to a single step pulse

Figure 12
Figure 12 Rectangular plate with a central crack under dynamic loading; a = 0.24 cm

Figure 14
Figure 14 Variation of normalized energy release rate with time for a rectangular plate with a central crack subjected to a dynamic load

Table 1 :
Energy release rates for a mode-I crack in a semi-infinite plate