# Numerical Procedures

%

There are various numerical approaches to solving linear elastic fracture problems. For our purposes, the best classification is based on commercial availability. Special purpose computer programs, which make use of special properties of the fields in plane elasticity (see Chapter 4) and are usually usually not available commercially, may be very accurate, but are generally restricted to research by specialists. This is the case of the so-called boundary collocation, in which special power expansions of the unknowns are used for a particular problem and the coefficients of the expansion are determined so that the boundary conditions be satisfied only in some average sense.

Restricting attention to commercial programs, we may distinguish between special purpose programs and general purpose programs. Special purpose programs are specifically designed to deal with cracks and determine stress intensity factors, so that the user may access to post-process routines that will readily
compute the stress intensity factor from a basic numerical solution. General purpose programs are those available to solve general elasticity problems, which can be used, with special strategies, to solve fracture problems. With the general purpose programs, two basic strategies may be used: (1) incremental stiffness method, and (2) near-tip field fitting.

The incremental stiffness method essentially consists of determining the compliance for two different, but close, crack lengths a — Aa and a 4- Да. Then the energy release rate G is estimated as P2_ C(a 4- Да) – C(a – Aa) 2b 2A a

Any finite element, finite difference, or boundary element code may be used to produce the two compliance values for two close crack lengths. Numerical resolution and mesh refinement limit the accuracy of this procedure. In very general terms, there are two possible approaches: (1) use only one mesh and simulate crack extension by freeing one node, so that the crack extends by one element, or (2) modify the mesh for the second calculation in which the node at the crack tip is displaced by Aa. The first method is easy to use and does not require modification of the global stiffness matrix, but requires a fine mesh so that the numerical differentiation in (3.3.1) gives accurate results. The second method decouples the crack extension from the mesh size, but requires partial recalculation of the stiffness matrix.

Example 3.3.1 A commercial finite element code was used to analyze a single-edge cracked beam in three-point bending, with a span-to-depth ratio S/D = 4 (Guinea 1990). Half the beam was discretized so that 100 equally sized elements were placed along the crack plane. The crack length was varied by changing the boundary conditions along the nodes in the crack plane, from opening displacement prevented (no crack at this node) to load free (crack at this node). Computations were performed in plane stress with D — 1,6= 1, and E — 1, for a load P = 1 (in arbitrary, but consistent, units). Although, the purpose of the computations was other than determining Ki, the results can be used to examine the accuracy of the differential stiffness method.

Consider, for example, the case а = 0.5 D. The displacements computed lor crack lengths a i = 0.49D аг = 0.5ID were, respectively, 57.261 and 61.663, numerically identical to the compliance values (because P — 1). The energy release rate is then evaluated from (3.3.1) as Q re (12/2)(61.663 — 57.261)/0.02 = 110.1 in appropriate units. Now, since we always write Q = (P/bD)g(a), it turns out that in our calculation the numerical value of Q (with its arbitrary units) coincides with the dimensionless value of з(0.5). Therefore, g{0.5) r; 1 10. The stress intensity factor follows front Irwin’s equation as К і — J~E’Q к, x/TTO = 10.5. Now, since we write Kj = a tts/~Dk(a), we may easily find k{ 0.5) upon noting that for a span-to-depth ratio of 4, a/v = 6P/bD. The result is k(0.5) ss 10.5/6 = 1.75. This value is to be compared with that given by equations (3.1.1) and (3.1.2) (or Fig. 3.1.1) which give k(0.5) = 1.77. Thus, the numerical estimate turns out to be about 1.2% lower than the more accurate value. D

The near-tip field fitting consists of making use of the known near-tip behavior of the stress and displacement or crack opening fields to make an estimate for Kj. It can make use, for example, of the stress distribution ahead of the crack tip, which is known to behave as 022 = Ki/sfhtr, where 022 is the stress normal to the crack plane, and r is the distance to the crack tip. This means that a plot of OijsJlitr vs. r should tend to К/ as r approaches zero. It is also possible to use the displacement field, particularly the crack opening, which is known to behave as w = %K/ s/r/E’fbt. Therefore, the limit of wE’sfbt j4fr as r —> 0 is also Kj.   Example 3.3.2 The results of the nodal reactions along the uncracked ligament, or the crack opening distribution, may be used to make a near-tip field fit. We use Kj — limr_,o aiis/lttv. and agree to write Kj — aNs/T)k(a) and for a given value of a, we define Figure 3.3.1 Plot of normal stresses times fr vs. r (r =distance to crack tip). Extrapolation to zero gives the dimensionless stress intensity factor.

then k(a) = lim,._o tc(r). Plotting the nodal values of ren vs. the distance to the crack tip, and extrapolating to zero, we get an estimate of k(a). This was used in the finite element computations described in the previous example (now for a/D = 0.5). The nodal normal stresses were obtained as ом,, ~ Rn/l>h, where Rn is the nodal reaction and h the width of the elements. From this, the nodal values of к were obtained and plotted as shown in Fig. 3.3.1. The extrapolated value gives k{ 0.5) яз 1.65. This value differs by 7% from the more accurate value k{0.5) = 1.77 obtained from equations (3.1.1) and (3.1.2) (or Fig. 3.1.1). D

The foregoing examples show two of the ways to determine Кj and Q from numerical results. The determination of Ki from the crack opening profile is left as one of the exercises. The general experience is that the differential stiffness method is more accurate for a given mesh size. This is probably due to the cancellation of constant errors in the differentiation process. However this method requires two computations, while the near-tip field fitting requires only one, although this is really not a problem with the kind of computers available today.

Getting good results (less than 5’% error) with near-tip analysis requires extremely fine meshes, because of the difficulty in representing the crack tip singularity with ordinary finite elements. Indeed, careful studies of convergence by Wilson (1971) and Oglesby and Lamackey (1972) showed that the near-lip approximate solution may not converge to the analytical solution whatever the mesh refinement. To solve this problem, one needs special elements whose shape functions include a r~singular term.

Various singular finite elements have been developed (see, e. g., Aliabadi and Rooke 1991), but most of them incorporate special shape function and require specially designed finite element codes. A remarkable exception is the sorcalled quarter-node isoparametric element (Barsoum 1975, 1976; Henshell and Shaw 1975). In this formulation, a standard 8-node isoparametric quadrilateral element is collapsed, as shown in Fig. 3.3.2a, to a triangular quarter-point element. The vertex corresponding to the collapsed nodes 1-7-8 becomes the singular point, and a r-1/2 singularity is achieved by placing nodes 2 and 6 at a quarter (from the singular vertex) of the radial sides of the triangle. These elements are placed in a rosette around the crack tip as shown in Fig. 3.3.2b.

The stress intensity factor may be evaluated from the displacement fields of any of the elements, but most usually K[ is obtained from the values of the crack opening evaluated at the two nodes along the crack faces. With this method, values of К j accurate within a few percent may be obtained without much mesh refinement. However, recent recommendations by ESIS Technical Committee 8 (1991) suggest, again, that best results for stress intensity factors are obtained if energetic approaches based on the determination of G are used instead of near-tip fields.

The differential stiffness method is not the-only way to determine G■ The J-integral and other path – independent integral expressions may be (and have been) used to determine the energy release rate. This has the advantage that the evaluation of J is made using values of the fields at points far from the crack tip, where the errors are expected to be smaller. It also avoids numerical differentiation, and a single  computation is enough. It requires, however, special postprocessing routines, both in finite elements codes and in boundary clement codes.

Although finite element codes dominate the market, commercial codes based on boundary elements have recently become available. They yield /</-values of much higher accuracy than finite element codes. The main advantage of this kind of formulation is that only the boundary of the elements must be modeled, so that the number of degrees of freedom is greatly reduced. This is, of course, achieved at the cost of larger complexity of the code, especially the postprocessing. In particular, handling cracks may require special formulations and special postprocessing which are outside the scope of this book (for details, see Aliabadi and Rooke 1991).