EXTRACTION OF FRACTURE MECHANICS PARAMETERS FROM FEM ANALYSIS: ALGORITHMS AND PROCEDURES

Full Text

Introduction

Defects and cracks play a decisive role in characterizing the strength and failure of these materials and structures [1–5]. Therefore, gaining insight on cracking processes is of crucial importance [1]. The first step in analyzing any fracturing process is to determine the crack tip asymptotic fields in order to characterize the stress, deformation and displacement near the crack tip, which requires the coefficients (unknowns) of the crack tip asymptotic field to be determined via reliable methods [1–10]. The first terms of the crack tip stress series expansion in isotropic linear elastic materials are singular, and, hence, dominant, in the proximate vicinity of the crack tip. Therefore, in the singular dominant zone the first terms are sufficient to characterize the crack tip fields. However, at further distances from the crack tip, the importance of the higher order terms become evident [6–11]. Thus, precise and simple algorithms are needed to reliably calculate coefficients in the multi-parameter crack-tip fields. Numerical methods and in particular finite element method (FEM) [6–11] allow us to extract the crack tip parameters. Moreover, it is worth noting that even determining the stress in-tensity factor is still the subject of investigations. For instance, in [12] direct extraction of stress in-tensity factors by a high-order numerical manifold method is realised. The proposed in [12] stress intensity factor (SIF) extraction method is shown to yield highly accurate results even without mesh refinement. Formulas extracting SIFs of the biharmonic equations on cracked domains with clamped (or simply supported or free) boundary conditions along the crack faces are derived in [13]. In [13] it is shown the iteration methods quickly converge and the proposed enrichment method yields highly accurate stress intensity factors. It is also demonstrated that for a known true solution, the extraction formulas yield exact stress intensity factor. Thus, the determination of SIF still raises questions. The determination of higher order coefficients requires more accurate approaches [14–24].

In this paper, a method for calculating the parameters of fracture mechanics based on finite-element analysis is proposed and tested. The efficiency of the method for extracting parameters of fracture mechanics is shown.

 

1. The multi-parameter crack tip stress field expansion

   The main objective of this paper is the numerical determination of higher-order coefficients of WE. The polar coordinate system r, θ is introduced and centered at the crack tip. In polar coordinates the Williams series solution for the near crack – tip stress field has the form [18; 23; 24]

   2

   σij (r, θ) = ∑

   ∞

   ∑ amf (k) (θ)rk/2−1, (1)

   k

   m=1 k=−∞

   m,ij

   k

    

   m,ij

    

   m,ij

    

   where index m is associated to the fracture mode; am are coefficients related to the geometric configuration, loads and fracture modes; f (k) (θ) are angular functions depending on stress components and mode. Analytical expressions for angular eigenfunctions f (k) (θ) are available [16; 17]:

   f (k)

    

   1,11(θ) = k [(2 + k/2 + (−1)k ) cos(k/2 − 1)θ − (k/2 − 1) cos(k/2 − 3)θ] /2,

   f (k)

    

   1,22(θ) = k [(2 − k/2 − (−1)k ) cos(k/2 − 1)θ + (k/2 − 1) cos(k/2 − 3)θ] /2,

   f (k)

    

   1,12(θ) = k [− (k/2 + (−1)k ) sin(k/2 − 1)θ + (k/2 − 1) sin(k/2 − 3)θ] /2,

   f (k)

    

   2,11(θ) = −k [(2 + k/2 − (−1)k ) sin(k/2 − 1)θ − (k/2 − 1) sin(k/2 − 3)θ] /2,

   f (k)

    

   2,22(θ) = −k [(2 − k/2 + (−1)k ) sin(k/2 − 1)θ + (k/2 − 1) sin(k/2 − 3)θ] /2,

   f (k)

    

   2,12(θ) = k [− (k/2 − (−1)k ) cos(k/2 − 1)θ + (k/2 − 1) cos(k/2 − 3)θ] /2.

    

   (2)

   The multi-parameter fracture mechanics concept consists in the idea that the crack-tip stress field is described by means of WE (1). In this work the central crack in an infinite plane medium is considered. Analytical determination of coefficients in crack-tip expansion for a finite crack in an infinite plane medium is given in [23; 24]:

   a1

    

   2n+1 = (−1)

   22

    

   n+1 (2n)!σ∞

   , a1 = −σ∞/4, a1

    

   = 0 (3)

    

   for Mode I crack loading,

   a1

   

   23n+1/2(n!)2(2n − 1)an−1/2

    

   n (2n)!σ∞

   2 22 2k

   

   2n+1 = (−1)

   , a

    

   12 1

   23n+1/2(n!)2(2n − 1)an−1/2 2k

   = 0 (4)

   for Mode II crack loading. The analytical solution (3), (4) allows us to validate the proposed method since one can compare the numerical results with the analytical ones. The crack length is less than the width

   Вестник Самарского университета. Естественнонаучная серия. 2020. Том 26, № 1. С. 69–77

   Vestnik of Samara University. Natural Science Series. 2020, vol. 26, no. 1, pp. 69–77 71

    

   and height of the plate. It is shown that the higher order terms in WE can play significant role in the description of the crack tip fields. Nowadays, various techniques are used to determine the parameters that characterize the crack-tip stress field. Now one can enumerate analytical [15; 16; 23; 24], experimental [7–11; 25; 26] and numerical [19] methods. One of the promising methods is FEM. One of the numerical examples discussed below is the large plate with the central crack. The finite element solution will be obtained and the results will be compared with the analytical formulae (3) and (4).

    

2. Finite element over-deterministic method

   As it is noted in [1] the basic principle of the finite element over-deterministic method is the use of a large number of FE data points in order to calculate the crack tip parameters. This is done by forming an algebraic system of equations where the number of equations is more than the number of unknowns. In this case the over-deterministic system of equations is encountered. In the framework of using the over-deterministic method to determine the coefficients of (1) nodal stresses can provide the necessary set of equations. The over deterministic technique assumes more equations than unknowns in order to obtain more accurate values. This implies that one can form an over deterministic system. Taking data from different points at different distances from the crack tip is allowed as higher order terms are included in the stress equations. The algorithm is implemented using in the mathematical software Maple. One can use the approach described in [12] and one can present eqn. (1) in the matrix form as

   σ = CA (5)

   The closed form solution of (5) for the unknown vector of fracture mechanics parameters can be written as

   A = (CT C)−1 CT σ (6)

   where (CT C)−1CT is the pseudo-inverse of C. The coefficients are estimated by minimizing the objective function which is of quadratic form for stress expression in terms of unknown parameters:

   T

   J (A) = (σ − CA)

   (σ − CA) /2. (7)

   Table 1

   Coefficients of multi-parameter Williams series expansion for a plate with a central crack of small length

    

   Fracture mechanics parameters	error
   a1 1/2 1 = 4.909MPam	0.01%
   a1 2 = −2.449MPa	0.09%
   a1 −1/2 3 = 2.484MPam	0.13%
   a1 −3/2 5 = −0.6236MPam	0.22%
   a1 −5/2 7 = 0.3112MPam	0.31%
   a1 −7/2 9 = −0.1951MPam	0.35%
   a1 −9/2 11 = 0.1361MPam	0.44%
   a1 −11/2 13 = −0.1056MPam	0.54%
   a1 −13/2 15 = 0.0786MPam	0.68%

    

3. Numerical examples

   The first example is the plate with the small central crack. In this work, 2D finite element analysis (FEA) of cracked specimens is carried out using Abaqus software to estimate SIF, T-stresses and coefficients of higher-order terms of WE. The analysis is done with 8-noded plane strain elements. The quarter point element is used to capture square root singularity at the crack tip. The center crack model is of dimension 400 mm × 400 mm having a crack of 10 mm length. The mesh pattern around the crack tip is kept very fine to capture the high-stress gradient. The mesh convergence is achieved with 72 elements along circumferential and 60 along the radial direction. In total, there are 13 344 elements. The typical finite element mesh is shown in fig. 1. To determine the higher order coefficients the stress tensor components from the nodes belonging to concentric circles are used. One can use different number of concentric circles. A class of numerical experiments with different numbers of concentric circles has been realized. The minimum number of stress tensor components was 219 since one can use the only circle with the following stress tensor components σ11, σ12, σ22. Increasing the number of considered concentric circles surrounding the crack tip one can enhance the dimension of the

   Bachareva Yu.N., Mironov A.V., Petrova D.M. Extraction of fracture mechanics parameters from fem analysis...

   72Бахарева Ю.Н., Миронов А.В., Петрова Д.М. Извлечение параметров механики разрушения...

    

   system (5). The maximum number of equations in (5) in the numerical experiments performed was 3492 from which the first fifteen coefficients of WE have been obtained.

    

   

   Fig. 1. Typical mesh containing singular elements near the crack tips

    

   The results of extraction of the coefficients of WE in the vicinity of the crack tip are given in Table 1 where the first column shows the coefficients of WE obtained from FEA whereas the second one shows the error in FEM comparatively with the analytical results given by formulae (4) and (5) for an infinite plate with the central crack.

    

   Crack tip fracture parameters for the SCB specimen

   Table 2

    

   Fracture mechanics parameters	n = 2	n = 4	n = 8
   KI (MPam1/2)	23.90	23.99	24.00
   KII (MPam1/2)	0.45	0.41	0.40
   a1 2(MPa)	-0.44	-0.456	-0.457
   a1 −1/2 3(MPam )		0.145	0.146
   a1 −1 4(MPam )		0.001	0.000
   a1 −3/2 5(MPam )			0.021
   a1 −2) 6(MPam			0.006
   a1 −5/2 7(MPam )			0.0004
   a1 −3 8(MPam )			0.0002

    

4. Extraction of the coefficients of the Williams series expansion for the semicircular bend specimen from the FEM analysis

   In this part of the paper the semicircular bend (SCB) specimen with an inclined crack shown in fig. 2 is studied. The following notations are adopted. P is the applied load, S is loading span in the SCB specimen, a is crack length, α is crack inclination angle. The semi-circular bend specimen subjected to three-point bending has received much attention in recent years for measuring the mixed mode I/II fracture resistance [16–19]. In this work, 2D FEA of semidisks with vertical crack and inclined notches is carried out using Abaqus software. To estimate SIF, T-stress and higher-order terms and verify the experimental results obtained FEM calculations have been employed. The analysis is done with 8-noded strain elements. The results of FEM analysis are shown in fig. 3,4. Fig. 3 (left) shows the distribution of the von Mises stress intensity. Fig. 3(right) shows the distribution of the stress component σ11. Figure 4 shows the distribution of stress component σ22.

    

5. Extraction the coefficients of the WE near the crack tip by digital photoelasticity

Photoelasticity is a whole field experimental technique to obtain stress fields in both 2-D and 3-D elasticity problems [18; 20; 25; 26]. Digital photoelasticity method has rapidly progressed in the last few years and has matured into an industry-friendly technique. Recently there has been a lot of works devoted to various aspects of the method and its applications [18; 20; 27]. The experimental setup is shown in fig. 5 (left). The experimental isochromatic fringe patterns in the plate with the central crack are shown in fig. 5. The

Вестник Самарского университета. Естественнонаучная серия. 2020. Том 26, № 1. С. 69–77

Vestnik of Samara University. Natural Science Series. 2020, vol. 26, no. 1, pp. 69–77 73

 

P

 



a

 

S

 

Fig. 2. Geometry of the semi-circular bend specimens

 

 

Fig. 3. Distributions of the von Mises equivalent stress (left) and the stress component σ11 (right)

 

over deterministic method has been applied to the experimental data obtained from the photoelasticity observations. The stress optic law relates the fringe order N and the in-plane principal stresses σ1, σ2 as Nfσ /t = σ1 − σ2, where fσ is the material stress fringe and t is the thickness of the specimen. The results of calculations are given in Table 3. +++

 

Conclusions

In this paper, we propose and describe an algorithm for constructing the stress field expansion coefficients at the crack tip from finite element calculation data. The algorithm is tested on several examples and the results are compared with the results of the photoelastic experiments. The comparison showed good agreement between the values of the coefficients of the multi-parametric asymptotic expansion. It is shown that higher approximations in the asymptotic expansion are especially significant when processing the entire set of experimental information. The example problems emphasise that the use of multi-parameter stress field is a practical necessity to apply concepts of Fracture Mechanics to solve real life engineering problems. It is

 

 

Fig. 4. The distribution of the stress component σ22

Bachareva Yu.N., Mironov A.V., Petrova D.M. Extraction of fracture mechanics parameters from fem analysis...

74Бахарева Ю.Н., Миронов А.В., Петрова Д.М. Извлечение параметров механики разрушения...

 

Fig. 5. Experimental setup of transmission photoelasticity

 

 

 

Fig. 6. Isochromatic images for 85 kg, 125 kg and 135 kg

 

Table 3

Coefficients of the Williams series expansion for the plate with the central crack with the geometric parameters as in the experimental photoelasticity method

the photoelasticity method FEM analysis

a1 1/2 1

1/2

1 = 7.2528MPam

a1

a1 = 7.2527MPam

1

2 = −2.7516MPa a2 = −2.7516MPa

a1 −1/2 1

−1/2

3 = 2.1406MPam

a1 −1

a3 = 2.0163MPam

1 −1

4 = −0.3370MPam

a4 = −0.3021MPam

a1 −3/2 1

−3/2

5 = −0.2844MPam

a5 = −0.2757MPam

a1 −2 1 −2

6 = −0.0919MPam

a6 = −0.0985MPam

a1 −5/2 1

−5/2

7 = 0.0765MPam

a7 = 0.0712MPam

a1 −3 1 −3

8 = 0.0255MPam

a1

−7/2

a8 = 0.0019MPam

1

−7/2

9 = −0.0340MPam

a9 = −0.0015MPam

a1 −4 1 −4

10 = 0.0255MPam

a10 = 0.0019MPam

a1 −9/2 1

−9/2

11 = 0.0098MPam

a11 = 0.0077MPam

a1 −5 1 −5

12 = 0.0019MPam

a12 = 0.0012MPam

a1 −11/2 1

−11/2

13 = 0.0056MPam

a13 = 0.00509MPam

a1 −6 1 −6

14 = 0.0008MPam

a1

−13/2

a14 = 0.0007MPam

1

−13/2

15 = 0.00181MPam

a15 = 0.00147MPam

shown that the multi-parameter ansatz allows us to collect data from a larger zone which helps to obtain accurate values of fracture mechanics parameters.

About the authors

Y. N. Bakhareva

Samara National Research University

Author for correspondence.
Email: bakharevayu@yandex.ru
ORCID iD: 0000-0002-6482-504X

PhD, associate professor of the Department of Mathematical Modeling in Mechanics

Russian Federation

A. V. Mironov

Samara National Research University

Email: Mironov.AV2020@yandex.ru
ORCID iD: 0000-0002-0666-9878

Master’s Degree student

Russian Federation

D. M. Petrova

Email: petrova.DA.2020@yandex.ru
ORCID iD: 0000-0002-8264-3426

Master’s Degree student

References

Supplementary files