COMPUTATIONAL FRACTURE MECHANICS WHY? - to compute fracture mechanics parameters (SIF, G) in 2D and 3D configurations; - to compute J integral and CTOD in elastic-plastic analyses; - to simulate crack growth (under general mixed-mode conditions); - to solve special problems: dynamic fracture, ductile fracture, cohesive fracture, fracture at interfaces, …. NUMERICAL METHODS: To determine the distribution of stresses and strains in a cracked body subject to external loads or displacements (when close-form analytical solutions are not available) - The Finite Difference Method; - The Finite Element Method; - The Boundary Element Methods. Computational FM 1 TECHNIQUES FOR COMPUTING FRACTURE MECHANICS PARAMETERS POINT MATCHING APPROACHES: 1) STRESS MATCHING 2) DISPLACEMENT MATCHING (The Displacement Correlation Technique) 3) STRESS FUNCTION MATCHING ENERGY APPROACHES: 4) THE GLOBAL ENERGY RELEASE METHOD 5) THE STIFFNESS DERIVATIVE TECHNIQUE (for FEM) 6) CONTOUR INTEGRATION 7) THE ENERGY DOMAIN INTEGRAL METHOD 8) THE MODIFIED CRACK-CLOSURE INTEGRAL TECHNIQUE Computational FM 2 POINT MATCHING APPROACHES: Idea: to correlate numerical solutions for stresses or displacements at specific locations with analytic solutions depending on the stress intensity factors. Hp: homogeneity, isotropy linearly elastic material 1) STRESS MATCHING hp: mode I loading. KI = lim ( σ 22 2π r r →0 ) (θ = 0) 2) DISPLACEMENT MATCHING hp: mode I loading. KI = κ = 3−4ν κ = (3−ν)/(1+ν) 2G lim 2π u 2 κ +1 r → 0 r (θ = π ) (plane strain) (plane stress) better than 1): higher precision in the numerical calculation of displacements (continuity between elements) Computational FM 3 3) MATCHING STRESS FUNCTIONS (The Boundary Collocation Method) - The Airy stress function is expressed in terms of complex polynomials. The coefficients of the polynomials are inferred from nodal quantities. SIF's are inferred from the stress functions. (highly cumbersome) Drawbacks: - High degree of mesh refinement is required for engineering accuracy even for simple geometry, loading and a single crack. - Use of special elements at the crack tip that exhibit the 1/√r singularity is necessary (enriched elements, quarter point elements). Computational FM 4 2) THE DISPLACEMENT CORRELATION TECHNIQUE (Ref. Chan, Tuba and Wilson, Eng. Fract. Mech., 2(1),1970) Hp: homogeneity, isotropy linearly elastic material General form of displacement function near crack tip: KI θ θ r sin (κ + 1) - 2cos 2 + .... 2G 2π 2 2 K r θ θ cos (κ - 1) - 2 sin2 + .... - II 2 2G 2π 2 u2 = KI r θ θ cos (κ - 1) + 2sin2 + .... 2 2 2G 2π K r θ θ sin (κ + 1) + 2 cos 2 + .... + II 2G 2π 2 2 u1 = κ = 3−4ν κ = (3−ν)/(1+ν) (plane strain) (plane stress) Consider a FE mesh at crack tip: Computational FM 5 Pure Mode I and plane stress: u2 = u1 = KI r θ 2 θ sin - cos 2 + .... G 2π 2 1+ ν 2 KI r θ 1 -ν θ cos + sin2 + .... G 2π 2 1 + ν 2 taking θ = ± π and r = rAB: K *I = 2π G(1+ ν ) (u2B - u 2C ) rAB 4 Pure Mode II and plane stress: u2 = K II G r θ ν -1 θ cos + sin2 + .... 2π 2 1 + ν 2 u1 = K II G θ 2 θ r sin + cos 2 + ... 2 2π 2 1+ ν taking θ = ± π and r = rAB: K *II = 2π G(1+ ν ) (u1B - u1C ) rAB 4 Mixed Mode and plane stress: Mode I and mode II problems uncouple at θ = ± π. Consequently previous equations are still applicable. Computational FM 6 For Quarter Points Singular Elements: (Shih, de Lorenzi and German, 1976, IJF, 12, 647-651) - If FE mesh is composed of quarter point elements at crack tip: K *I = 2π G [4(u2B - u2C ) - (u2D - u 2E )] rAD κ + 1 K *II = 2π G [4(u1B - u1C ) - (u1D - u1E )] rAD κ + 1 (from shape functions of element and analytical expressions for COD and CSD). Computational FM 7 FINITE ELEMENT MESH DESIGN - In homogeneous, isotropic, elastic material, the stress field at crack tip exhibits a 1/√r singularity. - The FE mesh must be such that the stress field characteristics are reproduced in the numerical solutions. Solutions: - Mesh refinement at crack tip (slow convergence of local parameters) - Use of “enriched” elements at crack tip (Wilson element; Tracey element, …): stress / displacement variations around crack tip are embedded in the shape functions of element (Tracey, 1971, Eng. Fract. Mech., 3, 255-266) (edge displacement incompatibility with surrounding elements) - Use of Quarter Point Singular Elements: in 8-nodes isoparametric quadrilateral elements, placing the mid-side nodes on or outside the ¼ of the side causes the Jacobian of transformation to become non positive definite. The 1/√r singularity is obtained. (Barsoum,1976, IJNME, 10, 25; Henshell and Shaw, 1975, IJNME, 9, 495507) (only nodal coordinate input data are altered. Element satisfies essential converge criteria) Computational FM 8 - If the 8-nodes element is first degenerated to an isoparametric triangle with 3 nodes at crack tip and then distorted by moving the side nodes, improved results are obtained: 1/√r singularity also within element not only along edges. - In 3D, analogous results obtained by distorting three dimensional 20-nodes isoparametric elements (Ingraffea and Manu,1980, IJNME, 15, 1427-1445) Computational FM 9 QUARTER POINT ELEMENTS: PROPERTIES Computational FM 10 Computational FM 11 ENERGY APPROACHES: 1) THE GLOBAL ENERGY RELEASE METHOD: Hp: 2D linear elastic body. From energy balance: G =- dW da W = total potential energy per unit width. - Perform two analyses, one for crack length a, the second for a +∆a, for constant load. Compute G from: G =- ∆W ∆L ∆U ∆U = = ∆a ∆a ∆a ∆a L = potential of the applied loads per unit width in a 2D body; U = elastic strain energy stored in the body per unit width - For pure mode I, compute SIF from: K I = EG (plane stress) EG KI = (plane strain) 1-ν 2 (hp: isotropy, homogeneity) Computational FM 13 Drawbacks: - two analyses required for KI computation. At least three analyses for KI and KII computations. Four analyses for KI, KII and KIII. - SIF's not directly computed. Postprocessing of results is required. - For mixed-mode fracture it is very difficult to separate G into its mode I and mode II components. Advantages: - no special crack-tip elements are necessary. Relatively coarse meshes can be used. - When used in conjunction with the Stiffness Derivative Technique, only partial reanalysis is necessary. Computational FM 14 2) THE STIFFNESS DERIVATIVE TECHNIQUE (Parks, IJF, 10, 1974, 487-502) - Formulated in terms of FE stiffness matrix. Not compatible with boundary element analysis. - Total potential energy of FE solution: W = 1/2 uT[K]u - uTP [K] = structure stiffness matrix; u = vector of nodal displacements; P = vector of applied nodal forces. - Differentiate to obtain G: ∂uT ([K]u − P) - 1 uT ∂[K] u + uT ∂P G =2 ∂a ∂a ∂a where first term on right hand side must be zero. - If applied load is independent of crack length and pure mode I: K I2 1 ∂[K] u G= = - uT E' 2 ∂a Computational FM 15 -If only elements surrounding the crack tip are modified by crack advance ∆a (see figure), then [K] remains unchanged in the region outside: K I2 1 G= = - uT E' 2 ∂[K ic ] 1 T u ≅ u ∑ ∂a 2 i=1 nc nc ( ) 1 ∑ ∆a [Kic ]a+ ∆a − [K ic ]a u i=1 [Kci] = stiffness matrix of the ith element surrounding the crack tip (within Γ1 in the figure); nc= number of such elements; [Kci]a element stiffness matrix for original crack configuration; [Kci]a+∆a element stiffness matrix for new crack configuration; Computational FM 16 6) CONTOUR INTEGRATION The J integral can be evaluated numerically along a contour (in 2D) or a surface (in 3D) surrounding the crack tip. Then in an elastic body: G =J Advantages: - Path independence enables the user to calculate J at a remote contour (surface) where numerical accuracy is greater. - No mesh refinement is required. - J can be separated into its mode I and mode II components to deduce KI and KII. Drawbacks: - Contour/surface integrals are difficult to implement numerically and the accuracy is not high. Solution: - The Energy Domain Integral: the J integral is formulated in terms of area (for 2D) or volume (for 3D) integrals. Better accuracy. Easier to implement. Computational FM 17 7) THE ENERGY DOMAIN INTEGRAL METHOD (Shih, Moran and Nakamura, 1986, IJF, 30, 79-102) The J integral for a 2D body is: ∂u ∂u J = ∫ Uddx 2 - T ds = ∫ Udn1 - σ ij i n jds ∂x1 ∂x1 Γ Γ - Consider the close path Γn = Γ0 + ΓS+ + ΓS- - Γ . - Introduce a weight function q(x1,x2) that is equal to unity on Γ and zero on Γ0 and Γs. The J integral along the new path is: ∂u J = - ∫ Udδ 1j - σ ij i q n jds ∂x1 Γn Computational FM 18 - Apply the divergence theorem to transform the integral along the close contour into a domain integral: ∂Ud ∂ ∂ui ∂ui ∂q ∂q σ q dA J = - ∫ Ud -σ dA - ∫ ∂x1 ij ∂x1 ∂x j ∂x1 ∂x j ij ∂x1 A A where the second term vanishes for elastic problems. - Implement the domain integral into a finite element code (e.g. Franc). Mode separation: (Bui, J.M.P.S., 31(6), 439-448, 1983) - The displacement field in the crack-tip region is decomposed into symmetric and antisymmetric components: u = usymm. + uantisymm. - The J integral is evaluated twice: once with the symmetric displacements to find J1 and a second time with the antisymmetric displacements to find J2. - For elastic problems: GI = J1 GII = J2 SIF’s are found from GI and GII. Computational FM 19 8) THE MODIFIED CRACK-CLOSURE INTEGRAL TECHNIQUE (Rybicki and Kanninen, 1977) - The strain energy release rate G is estimated in terms of the work done by the stresses ahead of the crack tip over the displacements produced by the introduction of a virtual crack extension (Irwin concept of crack-closure integral): lim 2 ∆ 1 G = σ 22 (x1,0) u2 ( ∆a - x1,0)dx 1 ∫ I ∆a → 0 ∆a 0 2 a lim 2 ∆ 1 G = σ 12 (x1,0) u1( ∆a - x1 ,0)dx 1 II ∆a → 0 ∆a 0∫ 2 a - Two numerical analyses are required to obtain the stress field ahead of the crack before propagation and to compute the displacement field after a virtual crack extension is introduced. Computational FM 20 - Simplification proposed by Rybicki and Kanninen: the displacement field ahead of the crack is approximated by the the displacement field behind the crack tip. The problem is then solved with one single analysis step. - In the finite element approach the crack closure integrals are rewritten in terms of equivalent nodal forces and relative nodal displacements (different expressions for each crack tip element). Advantages: - no assumption of isotropy or homogeneity around the crack tip is necessary. Computational FM 21 MIXED-MODE CRACK GROWTH - The direction of crack growth is defined using the Maximum Circumferential Stress Criterion. (Erdogan and Sih, 1963, J. Basic Engrg., 85, 519-527) - The crack is propagated of ∆a along the calculated direction. - Remeshing: a new mesh is generated around the new crack tip. Different techniques. Franc3D code: Solid Modeling Techniques to represent the desired geometrical aspects of the structure explicitly. The discretization is mapped onto this original solid model. Remeshing: the solid model remains valid and the new mesh is remapped into the solid model and inherits the necessary attributes. (see Franc3D manual for more details and references) - New fracture parameters at a + ∆a are calculated. Computational FM 22
© Copyright 2024