Next Article in Journal
Graphene Nanoplatelets Impact on Concrete in Improving Freeze-Thaw Resistance
Next Article in Special Issue
Understanding the Fracture Behaviors of Metallic Glasses—An Overview
Previous Article in Journal
An End-to-End Deep Learning Image Compression Framework Based on Semantic Analysis
Previous Article in Special Issue
An Elastic Interface Model for the Delamination of Bending-Extension Coupled Laminates
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Numerical Evaluation of SIFs of 2-D Functionally Graded Materials by Enriched Natural Element Method

Department of Naval Architecture and Ocean Engineering, Hongik University, Jochiwon, Sejong 30016, Korea
Appl. Sci. 2019, 9(17), 3581; https://doi.org/10.3390/app9173581
Submission received: 29 July 2019 / Revised: 23 August 2019 / Accepted: 25 August 2019 / Published: 1 September 2019
(This article belongs to the Special Issue Fatigue and Fracture of Non-metallic Materials and Structures)

Abstract

:

Featured Application

Prediction of crack propagation of functionally graded materials (FGMs).

Abstract

This paper presents the numerical prediction of stress intensity factors (SIFs) of 2-D inhomogeneous functionally graded materials (FGMs) by an enriched Petrov-Galerkin natural element method (PG-NEM). The overall trial displacement field was approximated in terms of Laplace interpolation functions, and the crack tip one was enhanced by the crack-tip singular displacement field. The overall stress and strain distributions, which were obtained by PG-NEM, were smoothened and improved by the stress recovery. The modified interaction integral M ˜ ( 1 , 2 ) was employed to evaluate the stress intensity factors of FGMs with spatially varying elastic moduli. The proposed method was validated through the representative numerical examples and the effectiveness was justified by comparing the numerical results with the reference solutions.

1. Introduction

In the late 1980s, a new material concept called functionally graded material (FGM) was proposed to resolve the inherent problem of traditional lamination-type composites [1]. The sharp material discontinuity across the layer interface causes the stress concentration, which may trigger the layer delamination. This crucial stress concentration can be significantly minimized by inserting a graded layer between two distinct homogeneous material layers [2,3]. This is because the material discontinuity completely disappears according to the material composition gradation, where the constituent particles of two base materials are mixed up in a random microstructure within a graded layer to maximize the desired performance [4,5,6]. Naturally, a functionally graded material is an inhomogeneous material, with spatially non-uniform material properties characterized by continuity and functionality. In addition to the suppression of stress concentration, the material concept of FGM rapidly and continuously spread throughout engineering and fields [7,8,9,10].
Early research efforts were concentrated on material characterization, fabrication, modeling, and analysis [1,11,12]. This was because the mechanical behaviors of FGMs are governed by the geometric dimensions and orientation, microstructure, and volume fractions of constituent particles. Recently, however, the concern toward the crack problems has increased because the structural failure of FGMs is dominated by micro-cracking [7,13,14]. In this regard, an accurate numerical prediction of stress intensity factors and the crack propagation has been an essential research subject [15,16]. For these subjects, an accurate reproduction of the 1 r singularity in the near-tip stress field in highly heterogeneous media becomes a key task [17,18,19,20].
To evaluate the stress intensity factors of FGMs with cracks, one can consider the use of the well-known J or M integral methods. However, these conventional indirect integral methods cannot reflect the spatially varying material properties of FGMs. The studies on the fracture mechanics of inhomogeneous bodies were initiated in the 1970~80s by assuming the spatially varying elastic modulus as an exponential function [17,21]. The standard integral methods were modified and/or refined to reflect the spatial non-uniformity of material properties by subsequent investigators. The most works were made by utilizing the finite element method [15,22,23,24,25]. However, since late 1990s, the employment of meshless methods to crack problems of inhomogeneous bodies has been actively advanced, particularly for the computation of SIFs by the modified M integral method [26,27]. Here, the extension of the natural element method (NEM) is worth noting, even though it was restricted to 2-D homogenous material [28].
The natural element method was introduced to overcome the demerits of conventional meshless methods [29], the difficulty in enforcing the displacement constraint and the numerical integration. The Laplace interpolation functions in NEM strictly obey the Kronecker delta property so that the imposition of displacement constraint becomes easy. In addition, Delaunay triangles, which were produced during the process for introducing Laplace interpolation functions, also serve as a background mesh for the numerical integration. In particular, PG-NEM can further improve the numerical integration accuracy by maintaining the consistency between the Delaunay triangle and the integration region [30]. Although Laplace interpolation functions provide the high smoothness of C 1 -continuity, there is still room for further improvement in capturing the high stress singularity at the crack tip.
In this context, this paper introduces an enriched PG-NEM to explore whether and how much the enrichment of interpolation function increases the prediction reliability of stress intensity factors for FGMs. The validity of enrichment was reported for homogeneous materials [31,32], but it was rarely reported for inhomogeneous materials. The trial function is enriched by the asymptotic displacement fields of mode I and II, and the approximated overall stress field is enhanced by the patch recovery technique. The proposed method was validated through the illustrative numerical examples and its effectiveness was quantitatively evaluated.

2. 2-D Inhomogeneous Cracked Bodies

2.1. Linear Elasticity of 2-D Cracked Bodies

Figure 1 represents a 2-D linearly elastic isotropic inhomogeneous material with an edge crack which is contained within a bounded domain Ω 2 with the boundary Ω = Γ D Γ N Γ c ¯ . Here, Γ D and Γ N indicate the displacement and force boundary regions, and Γ c = Γ c + Γ c ¯ denotes the traction-free crack surface. As a representative non-homogeneous material, FGMs are characterized by the spatially varying elastic modulus E and Poisson’s ratio ν over the bounded domain Ω . For the mathematical description purpose, two Cartesian co-ordinate systems are employed, { x , y } for the 2-D linear elasticity problem and { x , y } for the SIF evaluation of crack respectively. Assuming the crack surface is traction-free and neglecting the body force b , then the displacement field u ( x ) in the Cartesian coordinate system { x , y } is subjected to the equilibrium equations
σ = 0 i n Ω
with the displacement constraint
u = u ^ o n Γ D
and the force boundary condition
σ n = { t ^ o n Γ N 0 o n Γ c ±
Here, σ are the Cauchy stresses, n the outward unit vector normal to Ω , and t ^ the contour traction. When the displacement and strains are assumed to be small, the Cauchy strain ε is constituted to the Cauchy stress σ via the ( 3 × 2 ) gradient-like operator L such that
ε = ε ( u ) = L u
Letting D be the constitutive tensor, the stresses and strains are constituted by
σ = D : ε
Note that the displacement, strains, and stresses are calculated based on the co-ordinate system { x , y } and transformed into the co-ordinate system { x , y } .
For a homogeneous cracked body, the energy release rate per unit crack propagation along the x axis can be estimated by the path-independent J integral defined by
J = Γ ( W δ 1 j σ i j u i x 1 ) n j d s
using the indicial notation (i.e., x 1 = x and x 2 = y ), the Dirac delta function δ 1 j , and the strain-energy density W = σ ε / 2 = ε i j D i j k l ε k l / 2 . Here, Γ indicates an arbitrary closed path, which surrounds the crack tip counterclockwise. As shown in Figure 2, it is expanded to C = Γ + Γ c + Γ c + + Γ o in order to generate a grayed donut-type region A, in which a smooth function q ( x ) ( 0 q 1 ) is introduced to recast the integral into an equivalent domain form [33]. The function q , called by weighting function, becomes unity on Γ , zero on Γ o , and arbitrary value between 0 and 1 within the grayed donut region A. Then, the above Equation (6) can be expanded as following
J = A ( σ i j u i x 1 W δ 1 j ) q x j d A + A x j ( σ i j u i x 1 W δ 1 j ) q d A
according to the divergence theorem, together with n = n i on Γ in C . By further expanding the second term on the right hand side, Equation (7) becomes
J = A ( σ i j u i x 1 W δ 1 j ) q x j d A + A ( σ i j x j u i x 1 + σ i j 2 u i x j x 1 σ i j ε i j x 1 1 2 ε i j D i j k l x 1 ε k l ) q d A
But, the second integral on the right-hand side of Equation (8) becomes zero according to the equilibrium (1), compatibility (4), and the material uniformity. Therefore, the J integral for homogeneous materials becomes the area integral defined by
J = A ( σ i j u i x 1 W δ 1 j ) q x j d A
But, for non-homogeneous materials, the last material derivation term within the second integrand of Equation (8) does not become zero. Therefore, Equation (8) ends up with a more general J ˜ integral [34], which is given by
J ˜ = A ( σ i j u i x 1 W δ 1 j ) q x j d A A 1 2 ε i j D i j k l x 1 ε k l q d A
The last term in Equation (10) becomes extremely small as a contour Γ 0 shrinks to the crack tip, but its contribution is not negligible when the domain of integration A is reasonably large.

2.2. Modified Interaction Integral M ˜ ( 1 , 2 )

In order to extract K I and K I I from J integral, one needs to employ the interaction integral, in which two equilibrium states of a cracked body are considered. State 1 is the real equilibrium state of a body that is subjected to the prescribed boundary conditions, while state 2 stands for an auxiliary equilibrium state which would be the asymptotic near-tip fields for modes I or II. Another equilibrium state, called state S, could be established by combining these two states, for which the J ˜ integral in Equation (10) is rewritten as [21]
J ˜ ( S ) = A ( ( σ i j ( 1 ) + σ i j ( 2 ) ) ( u i ( 1 ) + u i ( 2 ) ) x 1 W ( S ) δ 1 j ) q x j d A + A x j ( ( σ i j ( 1 ) + σ i j ( 2 ) ) ( u i ( 1 ) + u i ( 2 ) ) x 1 W ( S ) δ 1 j ) q d A
with three different strain energy densities, W ( 1 ) , W ( 2 ) and W ( 1 , 2 ) , defined by
W ( S ) = 1 2 ( σ i j ( 1 ) + σ i j ( 2 ) ) ( ε i j ( 1 ) + ε i j ( 2 ) ) = 1 2 σ i j ( 1 ) ε i j ( 1 ) + 1 2 σ i j ( 2 ) ε i j ( 2 ) + 1 2 ( σ i j ( 1 ) ε i j ( 2 ) + σ i j ( 2 ) ε i j ( 1 ) ) = W ( 1 ) + W ( 2 ) + W ( 1 , 2 )
By utilizing the equilibrium (1) (i.e., σ i j ( ) / x j = 0 ) and the compatibility (4) (i.e., ε i j ( ) = u i ( ) / x j ), Equation (11) can be further simplified as
J ˜ ( S ) = A ( ( σ i j ( 1 ) + σ i j ( 2 ) ) ( u i ( 1 ) + u i ( 2 ) ) x 1 ( W ( 1 ) + W ( 2 ) + W ( 1 , 2 ) ) δ 1 j ) q x j d A   + A 1 2 ( ε i j ( 1 ) D i j k l x 1 ε k l ( 1 ) + σ i j ( 1 ) ε i j ( 2 ) x 1 σ i j ( 2 ) x 1 ε i j ( 1 ) + σ i j ( 2 ) ε i j ( 1 ) x 1 σ i j ( 1 ) x 1 ε i j ( 2 ) ) q d A = J ˜ ( 1 ) + J ˜ ( 2 ) + M ˜ ( 1 , 2 )
Here,
J ˜ ( 1 ) = A ( σ i j ( 1 ) u i ( 1 ) x 1 W ( 1 ) δ 1 j ) q x j d A 1 2 A ε i j ( 1 ) D i j k l x 1 ε k l ( 1 ) q d A
J ˜ ( 2 ) = A ( σ i j ( 2 ) u i ( 2 ) x 1 W ( 2 ) δ 1 j ) q x j d A
denote the J ˜ integrals for state 1 and state 2 respectively, and
M ˜ ( 1 , 2 ) = A ( σ i j ( 1 ) u i ( 2 ) x 1 + σ i j ( 2 ) u i ( 1 ) x 1 W ( 1 , 2 ) δ 1 j ) q x j d A + A 1 2 ( σ i j ( 1 ) ε i j ( 2 ) x 1 σ i j ( 2 ) x 1 ε i j ( 1 ) + σ i j ( 2 ) ε i j ( 1 ) x 1 σ i j ( 1 ) x 1 ε i j ( 2 ) ) q d A
indicates the modified interaction integral. All the quantities in Equation (16) are evaluated based on a coordinate system aligned to the crack tip, and the identification of domain of integration A and the weighting function q ( x ) will be described in details in Section 3.2.
Since two J ˜ integrals for inhomogeneous cracked bodies which are subject to mixed-mode loading also represent the energy release rates, one can obtain
J ˜ ( 1 ) = 1 E ¯ t i p ( K I ( 1 ) 2 + K I I ( 1 ) 2 )
J ˜ ( 2 ) = 1 E ¯ t i p ( K I ( 2 ) 2 + K I I ( 2 ) 2 )
And
J ˜ ( S ) = J ˜ ( 1 ) + J ˜ ( 2 ) + 2 E ¯ t i p ( K I ( 1 ) K I ( 2 ) + K I I ( 1 ) K I I ( 2 ) )
where E ¯ t i p at the crack tip is E t i p for plane stress, while it is E t i p / ( 1 ν 2 ) for plane strain. By equating Equation (13) with Equation (19), one can obtain
M ˜ ( 1 , 2 ) = 2 E ¯ t i p ( K I ( 1 ) K I ( 2 ) + K I I ( 1 ) K I I ( 2 ) )
In 2-D linear fracture mechanics, the closed-form near-tip displacement fields are available for mode I and mode II [35]. The stress intensity factor K I ( 1 ) for mode I can be obtained by letting state 2 be the pure mode-I asymptotic field (i.e., K I ( 2 ) = 1 and K I I ( 2 ) = 0 ):
M ( 1 ,   ModeI ) = 2 E ¯ t i p K I ( 1 )
In a similar manner, the stress intensity factor K I I for mode-II can be also determined.

3. Enriched Petrov-Galerkin Natural Element Method

3.1. Enriched NEM Approximation

The boundary value problem expressed by Equations (1)–(3) is converted to the weak form according to the virtual work principle: Find u ( x ) such that
Ω ε ( v ) : σ ( u ) d Ω = Γ N t ^ v d s
for all the test displacements v ( x ) in the Cartesian coordinate system { x , y } . Next, a non-convex NEM grid N E M , shown in Figure 3a, is constructed using N nodes and a number of Delaunay triangles. Then, for a given NEM grid, trial and test displacement fields u ( x ) and v ( x ) for PG-NEM approximation are interpolated as
u h ( x ) = J = 1 N u ¯ J ϕ J ( x ) + k 1 [ Q 11 ( x ) Q 12 ( x ) ] + k 2 [ Q 21 ( x ) Q 22 ( x ) ]
v h ( x ) = I = 1 N v ¯ I Ψ I ( x ) + λ 1 [ Q 11 ( x ) Q 12 ( x ) ] + λ 2 [ Q 21 ( x ) Q 22 ( x ) ]
with Laplace interpolation functions ϕ J ( x ) represented in Figure 3b. Here, ψ I ( x ) are constant-strain FE basis functions, which are defined over Delaunay triangles.
In addition, k 1 and k 2 are two unknown constants associated with a crack, and Q 1 α ( x ) and Q 2 α ( x ) represent the near-tip singular displacement fields given by
Q 11 ( x ) = 1 2 μ r 2 π c o s ( θ 2 ) [ κ 1 + 2 s i n 2 ( θ 2 ) ]
Q 12 ( x ) = 1 2 μ r 2 π s i n ( θ 2 ) [ κ + 1 2 c o s 2 ( θ 2 ) ]
Q 21 ( x ) = 1 2 μ r 2 π s i n ( θ 2 ) [ κ + 1 + 2 c o s 2 ( θ 2 ) ]
Q 22 ( x ) = 1 2 μ r 2 π c o s ( θ 2 ) [ κ 1 2 s i n 2 ( θ 2 ) ]
Referring to Figure 1, r is the radial distance from the crack tip, while θ is the angle from the x axis. Meanwhile, μ indicates the shear modulus, and κ is the Kolosov constant given by
κ = { 3 4 ν for   plane   strain ( 3 ν ) / ( 1 + ν ) for   plane   stress
with ν being the Poisson’s ratio. From Equation (23), the overall nodal coefficients u ¯ J are defined by
u ¯ J ( x ) = u h ( x J ) k 1 [ Q 11 ( x J ) Q 12 ( x J ) ] k 2 [ Q 21 ( x J ) Q 22 ( x J ) ]
and the overall stress and strain fields corresponding to u ¯ J are recovered by the global recovery technique. Meanwhile, the original essential boundary condition (2) is modified as
u ^ = u ^ k 1 [ Q 11 ( x ^ ) Q 12 ( x ^ ) ] k 2 [ Q 21 ( x ^ ) Q 22 ( x ^ ) ] , x ^ o n Γ D
Substituting Equations (23) and (24) into Equation (22), together with Equations (4) and (5), ends up with
[ K I J K I 1 K I 2 K I 1 T K 11 K 12 K I 2 T K 12 T K 22 ] { u ¯ J k 1 k 2 } = { F I f 1 f 2 }
In which the global stiffness matrices K I J and load vectors F I are defined by
K I J = Ω B I T D B J d Ω , F I = Γ N Φ I T t ^ d s ,    I , J = 1 , 2 , , N
where
B I T = [ ϕ I . x 0 ϕ I , y 0 ϕ I , y ϕ I , x ]
D = E ( 1 + ν ) ( 1 2 ν ) [ 1 ν ν 0 ν 1 ν 0 0 0 1 2 ν 2 ]   for   plane   strain
Meanwhile, the enriched stiffness matrices K α β , load vectors f α , and the interface matrices K I α are defined by, respectively
K α β = Ω B ^ α T D B ^ β d Ω ,   f α = Γ N ϕ ^ α T t ^ d s , α , β = 1 , 2
K I α = Ω B I T D B ^ α d Ω
with
ϕ ^ α T = [ Q α 1 ( x ) , Q α 2 ( x ) ]
B ^ α = [ Q α α , x ( x ) Q α β , y ( x ) Q α α , y ( x ) + Q α β , y ( x ) ] , α α = no   sum

3.2. Numerical Implementation of M ˜ ( 1 , 2 )

Figure 4 schematically represents the identification of domain of integration A and the definition of weighting function q ( x ) , which are prerequisite for the numerical implementation of M ˜ ( 1 , 2 ) in Eq. (6). For these two things, a circle of the radius r i n t , which is originated at the crack tip, is first imaginarily drawn on the NEM grid. Then, the nodal values of q ( x ) are assigned based on this circle, such that unity is specified to all the interior nodes while zero to the outer remaining nodes. Then, referring to Figure 2, the boundary of a square composed of eight darkened Delaunay triangles serves as an internal path Γ . On the other hand, the grayed region outside the internal path Γ automatically becomes a domain of integration A . According to the properties of Laplace interpolation functions [29], the defined weighting function q ( x ) becomes unity within a set of darkened Delaunay triangles while it has the values between zero and unity within the grayed domain of integration A .
Both the weighting function q ( x ) and the derivative q / x j become zero outside the domain of integration A . Furthermore, the derivative q / x j vanishes within the darkened squares. So, only the grayed Delaunay triangles are taken for the modified interaction integral M ˜ ( 1 , 2 ) such that
M ˜ ( 1 , 2 ) = K = 1 M G M ˜ K ( 1 , 2 )
Here, M G indicates the total number of grayed Delaunay triangles, while M ˜ K ( 1 , 2 ) stand for the triangle-wise modified interaction integrals, which are defined by
M ˜ K ( 1 , 2 ) = = 1 I N T { ( [ σ i j ( 1 ) u i ( 2 ) x 1 + σ i j ( 2 ) u i ( 1 ) x 1 W ( 1 , 2 ) δ 1 j ] x q x j | x w | J | x   + 1 2 [ σ i j ( 1 ) ε i j ( 2 ) x 1 σ i j ( 2 ) x 1 ε i j ( 1 ) + σ i j ( 2 ) ε i j ( 1 ) x 1 σ i j ( 1 ) x 1 ε i j ( 2 ) ] q ( x ) w | J | x
with I N T , x , and w being the total number, co-ordinates, and weights of sampling points for Gauss numerical integration for triangles. Meanwhile, T K in Figure 4 indicates the geometry transformation between the master element Ω ^ and the grayed triangle Ω K and the master element Ω ^ , which is defined by
T K :   x = i = 1 3 x i ψ i ( ξ , η ) , y = i = 1 3 y i ψ i ( ξ , η )
where ψ i stand for the constant-strain FE basis functions, which are employed for the test displacement field v ( x ) .

4. Numerical Experiments

The enriched NEM module was developed in Fortran and combined into the previously developed PG-NEM code [30] having the stress recovery module [36]. To demonstrate and validate the present method, the crack problem of an infinite plate with collinear cracks depicted in Figure 5a is considered. Differing from the cracks in plates of finite size which are of great practical interest, the present infinite plate with periodic cracks provides the closed form solution. Meanwhile, for the numerical study, one may take two kinds of periodic finite crack problems. One is obtained when the plate is cut along EF and GH, while the other is obtained if the plate is cut along AB and CD. Since the first one has been widely taken for the numerical studies (i.e., Ryicki and Kanninen [37] and Chow and Atluri [38]), the second periodic analysis model is taken for the current numerical experiments. Figure 5b represents the detailed second periodic model, as well as a NEM grid having higher grid density toward the crack tip for an upper left darkened quarter. The SIF for this infinite plate with collinear cracks is expressed in the following closed form (the correction factor C of 1.02) [39]:
K I = σ π a ( 2 b π a t a n π a 2 b ) C
with a being the half crack length, where a circular path Γ around the left crack tip has the relative radius r / a = 0.0139 , and the near-tip stress distributions are to be evaluated along this circular path. The Young’s modulus E is 200 GPa and Poisson’s ratio ν is 0.3.
The problem was solved with 13 Gauss integration points, and the stress field approximated by PG-NEM was smoothened by the stress recovery technique [36] with the Laplace interpolation functions. The stress distributions were evaluated along the circular path Γ and are compared in Figure 6. It can be observed that the case without enrichment shows similar stress distributions to the exact ones, but it provides a significantly low stress level. Thus, one can realize that only the fine gradient grid cannot sufficiently capture the near-tip stress distribution. Meanwhile, it is clear that the case with enrichment remarkably improves the stress level. Therefore, it has been justified that the enrichment of interpolation functions is effective in enhancing the near-tip stress approximation.
Figure 7a shows a rectangular FGM plate with an edge crack which is in the plane strain condition, where W and H / W are 1.0 and 8.0 units, while a / W is set variable for the parametric investigation as follows: 0.2, 0.3, 0.4, 0.5, and 0.6. The Young’s modulus E varies along the x axis in a form of
E ( x ) = E 1 e x p ( β x ) , 0 x W
but Poisson’s ratio is uniform as 0.3. When E 1 and E 2 are specified to the elastic moduli at x = 0 and x = L , respectively, the parameter β is calculated according to β = l n ( E 2 / E 1 ) , with E 1 being 1.0 unit. Two material variation cases of E 2 / E 1 are taken as follows: E 2 / E 1 = 0.1 and 10.0, and two loading cases of tension in Figure 7b and bending in Figure 7c are considered. This example was initially studied by Erdogan and Wu [40], who also provided an analytical solution. The external loading is set by as follows: for tension and σ y y ( x , ± 4 ) = ± 1.0 for tension and σ y y ( x , ± 4 ) = ± ( 2 x + 1 ) for bending.
An upper grayed half is taken for the crack analysis from the problem symmetry, and the following symmetric constraint is specified to the symmetric line. The non-cracked symmetric line is subjected to the vertical displacement constraint of u y = 0 , and the right end point is enforced by the lateral displacement constrain of u x = 0 . A 11 × 41 uniform NEM grid is used, and all the NEM analyses, the stress recovery, and the modified interaction integrals M ˜ ( 1 , 2 ) were carried out with 13 Gauss points. The radius r i n t for determining the domain of integration A was chosen using twice the side length of a square consisting of two Delaunay triangles [41].
Table 1 comparatively represents the normalized mode-I SIFs for uniform tensile loading and bending loading, respectively, where the values in parenthesis indicate the relative differences (%). It is confirmed that the proposed enriched PG-NEM provides the SIFs, which are consistent with the analytical solutions [40], for all the combinations of E 2 / E 1 and a / W of two types of loading, except for the relatively remarkable relative difference 6.807% in the tension case with E 2 / E 1 = 10 and a / W = 0.3 . For the further comparison, the numerical results by Kim and Paulino [15] and by Rao and Rahman [34] are also given in Table 1. It is observed that the SIFs of proposed method are correlated satisfactorily with the FEM results by the J k * -integral and the EFGM results by the modified M ˜ ( 1 , 2 ) -interaction integral. Table 2 represents the comparison with the results obtained without using the enrichment. In the uniform tensile loading, the case without enrichment provides the SIFs, which are lower than those with enrichment, such that the range of relative differences is from 27.59% and 51.04%. Meanwhile, in the bending, it is observed that the difference between with and without enrichment is not significant. Furthermore, the case without using enrichment provides the SIFs which are higher than those obtained by the enrichment. It is inferred that the near-tip stress field produced by bending is totally different from that by uniform tensile loading. The relative difference becomes smaller as the relative crack length a / W increases. It is observed from the detailed numerical values that the relative difference ranges from 2.41% to 24.68%. Thus, it has been found that the proposed enrichment method is influenced by the loading type.
The third example is a slant edge crack in a 2-D FGM plate in plane stress condition with the height H = 2 units and the width W = 1 unit. Referring to Figure 8a, the crack angle α is 45 ° and the relative crack length is a / W is 0.4 2 , where the enrichment and the crack evaluation are performed within the inclined Cartesian co-ordinate system { x , y } originated at the crack tip O . The Poison’s ratio is kept uniform as ν = 0.3 , but the elastic modulus varies in an exponential function along the x direction
E ( x ) = E ¯ e x p [ η ( x 0.5 ) ] , 0 x W
Here, E ¯ = 1 unit and η is variable for the parametric study such as η a = 0 , 0.1 , 0.25 , 0.5 , 0.75, and 1.0. As external load, an exponentially varying distributed load is applied to the upper edge with σ y y ( x , 1 ) = ε ¯ E ¯ e x p [ η ( x 0.5 ) ] , with ε ¯ being 1. The whole bottom edge is constrained in the vertical direction (i.e., u y = 0 ). At the same time, the right end is constrained in the horizontal direction. A uniform NEM grid, given in Figure 8b, is employed which is discretized by 11 × 21 , where a darkened area indicates the domain A for the interaction integral. The total number of grid points is 235 because four additional grid points are needed along the crack to split a crack line into upper and lower faces. All the NEM analyses, patch recovery, and the modified interaction integrals were performed with 13 Gauss points.
Table 3 comparatively represents the predicted normalized mixed-mode SIFs K I / ε ¯ E ¯ π a and K I I / ε ¯ E ¯ π a evaluated by three different methods, the present PG-NEM, J k * integral using FEM, and EFGM for seven different η a values. It is confirmed that the present method is in good agreement with FEM and EFGM, such that the biggest relative difference is 4.37% at K I I / ε ¯ E ¯ π a for η a = 0.25 . The values in [*] are the normalized mixed-mode SIFs obtained by PG-NEM without using enrichment. It is observed that the values are significantly smaller than those obtained using enrichment, such that the relative difference ranges from 23.32% to 65.63%. The relative difference increases in proportional to the exponential index η a . Hence, it was confirmed that the enrichment successfully and significantly improves the prediction accuracy of mixed-mode SIFs of the highly heterogeneous FGM, even for the coarse NEM grid.

5. Conclusions

In this paper, an enriched PG-NEM was introduced for the reliable crack analysis of inhomogeneous functionally graded materials (FGMs). The Laplace interpolation function was enhanced by the crack-tip singular displacement field and the essential boundary condition was modified accordingly. The unreached global displacement field was extracted, and its stress field was smoothened by the stress recovery. The validity and effectiveness of proposed method was illustrated and justified through three representative examples.
Through the numerical experiments, it was found that enriched PG-NEM is more effective to represent the near-tip stress singularity. The enrichment provides the stress level, which is more close to the exact solution when compared to the use of locally refined NEM grid. In addition, it was observed that enriched PG-NEM accurately predicts the mixed-mode SIFs of inhomogeneous functionally graded materials with exponentially varying elastic modulus. When compared to the case without using enrichment, the prediction accuracy was improved up to four times. Meanwhile, with respect to the analytical solution, the maximum relative difference was found to be less than 6.5%.

Funding

This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (Grant No. NRF-2017R1D1A1B03028879).

Conflicts of Interest

The author declares no conflicts of interest.

References

  1. Miyamoto, Y.; Kaysser, W.W.; Rabin, B.H.; Kawasaki, A.; Ford, R.G. Functionally Graded Materials: Design, Processing and Applications; Springer Science Business Media: New York, NY, USA, 1999. [Google Scholar]
  2. Giannakopolus, A.E.; Suresh, S.; Olsson, M. Elastoplastic analysis of thermal cycling: Layered materials with compositional gradients. Acta Metall. Mater. 1995, 43, 1335–1354. [Google Scholar] [CrossRef]
  3. Cho, J.R.; Oden, J.T. Functionally graded material: A Parametric study on thermal stress characteristics using the Crank-Nicolson-Galerkin scheme. Comput. Method. Appl. Mech. Eng. 2000, 188, 17–38. [Google Scholar] [CrossRef]
  4. Reiter, T.; Dvorak, G.J. Micromechanical models for graded composite materials: II. Thermomechanical loading. J. Phys. Solids 1998, 46, 1655–1673. [Google Scholar] [CrossRef]
  5. Cho, J.R.; Ha, D.Y. Volume fraction optimization for minimizing thermal stress in Ni-Al2O3 functionally graded materials. Mater. Sci. Eng. A 2002, 334, 147–155. [Google Scholar] [CrossRef]
  6. Apalak, M.K. Functionally graded adhesively bonded joints. Rev. Adhesion Adhesives 2014, 1, 56–84. [Google Scholar] [CrossRef]
  7. Zhang, C.; Sladek, J.; Sladek, V. Crack analysis in unidirectionally and bidirectionally functionally graded materials. Int. J. Fract. 2004, 129, 385–406. [Google Scholar] [CrossRef]
  8. Yahia, S.A.; Atmane, H.A.; Houari, M.S.A.; Tounsi, T. Wave propagation in functionally graded plates with porosities using various higher-order shear deformation plate theories. Struct. Eng. Mech. 2015, 53, 1143–1165. [Google Scholar] [CrossRef]
  9. Abdelaziz, H.H.; Meziane, M.A.A.; Bousahla, A.A.; Tounsi, A.; Mahmoud, S.R.; Alwabli, A.S. An efficient hyperbolic shear deformation theory for bending, buckling and free vibration of FGM sandwich plates with various boundary conditions. Steel Comp. Struct. 2017, 25, 693–704. [Google Scholar]
  10. Bourada, F.; Bousahla, A.A.; Bourada, M.; Azzaz, A.; Zinata, A.; Tounsi, A. Dynamic investigation of porous functionally graded beam using a sinusoidal shear deformation theory. Wind Struct. 2019, 28, 19–30. [Google Scholar]
  11. Cho, J.R.; Ha, D.Y. Averaging and finite-element discretization approaches in the numerical analysis of functionally graded materials. Mater. Sci. Eng. A 2001, 302, 187–196. [Google Scholar] [CrossRef]
  12. Birman, V.; Byrd, L.W. Modeling and analysis of functionally graded materials and structures. Appl. Mech. Rev. 2007, 60, 195–216. [Google Scholar] [CrossRef]
  13. Mahamood, R.M.; Akinlabi, E.T.; Shukla, M.; Pityana, S. Functionally graded material: An overview. In Proceedings of the World Congress on Engineering 2012 (WCE 2012), London, UK, 4–6 July 2012. [Google Scholar]
  14. Ivanov, I.V.; Sadowski, T.; Pietras, D. Crack propagation in functionally raded strip under thermal shock. Eur. Phys. J. Spec. Top. 2013, 222, 1587–1595. [Google Scholar] [CrossRef]
  15. Kim, J.H.; Paulino, G.J. Finite element evaluation of mixed mode stress intensity factors in functionally graded materials. Int. J. Numer. Methods Eng. 2002, 53, 1903–1935. [Google Scholar] [CrossRef]
  16. Tilbrook, M.T.; Moon, R.J.; Hoffman, M. Crack propagation in graded composites. Compos. Sci. Technol. 2005, 65, 201–220. [Google Scholar] [CrossRef]
  17. Delale, F.; Erdogan, F. The crack problem for a nonhomogeneous plane. J. Appl. Mech. 1983, 50, 609–614. [Google Scholar] [CrossRef]
  18. Eischen, J.W. Fracture of nonhomogeneous materials. Int. J. Fract. 1987, 34, 3–22. [Google Scholar]
  19. Anlas, G.; Lambros, J.; Santare, M.H. Dominance of asymptotic crack tip fields in elastic functionally graded materials. Int. J. Fract. 2002, 115(2), 193–204. [Google Scholar] [CrossRef]
  20. Ayhan, A.O. Stress intensity factors for three-dimensional cracks in functionally graded materials using enriched finite elements. Int. J. Solids Struct. 2007, 44, 8579–8599. [Google Scholar] [CrossRef] [Green Version]
  21. Atkinson, C.; List, R.D. Steady state crack propagation into media with spatially varying elastic properties. Int. J. Eng. Sci. 1978, 16, 717–730. [Google Scholar] [CrossRef]
  22. Gu, P.; Dao, M.; Asaro, R.J. A simplified method for calculating the crack tip field of functionally graded materials using the domain integral. J. Appl. Mech. 1999, 66, 101–108. [Google Scholar] [CrossRef]
  23. Anlas, G.; Santare, M.H.; Lambros, J. Numerical calculation of stress intensity factors in functionally graded materials. Int. J. Fract. 2000, 104, 131–143. [Google Scholar] [CrossRef]
  24. Dolbow, J.E.; Gosz, M. On the computation of mixed-mode stress intensity factors in functionally graded materials. Int. J Solids Struct. 2002, 39, 2557–2574. [Google Scholar] [CrossRef]
  25. Kaczmarczyk, J.; Grajcar, A. Numerical simulation and experimental investigation of cold-rolled steel cutting. Materials 2018, 11, 1263. [Google Scholar] [CrossRef] [PubMed]
  26. Rao, B.N.; Rahman, S. Mesh-free analysis of cracks in isotropic functionally graded materials. Eng. Fract. Mech. 2003, 70, 1–27. [Google Scholar] [CrossRef]
  27. Liu, K.Y.; Long, S.Y.; Li, G.Y. A meshless local Petrov-Galerkin method for the analysis of cracks in the isotropic functionally graded material. Comp. Model. Eng. Sci. 2008, 7, 43–57. [Google Scholar]
  28. Cho, J.R.; Lee, H.W. Calculation of stress intensity factors in 2-D linear fracture mechanics by Petrov–Galerkin natural element method. Int. J. Numer. Methods Eng. 2014, 98, 819–839. [Google Scholar] [CrossRef]
  29. Sukumar, N.; Moran, A.; Belytschko, T. The natural element method in solid mechanics. Int. J. Numer. Methods Eng. 1998, 43, 839–887. [Google Scholar] [CrossRef]
  30. Cho, J.R.; Lee, H.W. A Petrov–Galerkin natural element method securing the numerical integration accuracy. J. Mech. Sci. Technol. 2006, 20, 94–109. [Google Scholar] [CrossRef]
  31. Fleming, M.; Chu, Y.A.; Moran, B.; Belytschko, T. Enriched element-free Galerkin methods for crack tip fields. Int. J. Numer. Methods Eng. 1997, 40, 1483–1504. [Google Scholar] [CrossRef]
  32. Pant, M.; Singh, I.V.; Mishra, B.K. A novel enrichment criterion for modeling kinked cracks using element free Galerkin method. Int. J. Mech. Sci. 2013, 68, 140–149. [Google Scholar] [CrossRef]
  33. Moran, B.; Shih, F. Crack tip and associated domain integrals form momentum and energy balance. Eng. Fract. Mech. 1987, 27, 615–642. [Google Scholar] [CrossRef]
  34. Rao, B.N.; Rahman, S. An efficient meshless method for fracture analysis of cracks. Comput. Mech. 2000, 26, 398–408. [Google Scholar] [CrossRef]
  35. Anderson, T.L. Fracture Mechanics: Fundamentals and Applications, 1st ed.; CRC Press: Boca Raton, FL, USA, 1991. [Google Scholar]
  36. Cho, J.R. Stress recovery techniques for natural element method in 2-D solid mechanics. J. Mech. Sci. Technol. 2016, 30, 5083–5091. [Google Scholar] [CrossRef]
  37. Ryicki, E.F.; Kanninen, M.F. A finite element calculation of stress intensity factors by a modified crack closure integral. Eng. Fract. Mech. 1977, 9, 931–938. [Google Scholar] [CrossRef]
  38. Chow, W.T.; Atluri, S.N. Finite element calculation of stress intensity factors for interfacial crack using virtual crack closure integral. Comput. Mech. 1995, 16, 417–425. [Google Scholar] [CrossRef]
  39. ASTM. Fracture Toughness Testing and Its Applications; ASTM International: West Conshohocken, PA, USA, 1965; Volume 381, pp. 43–51. [Google Scholar]
  40. Erdogan, F.; Wu, B.T. The surface crack problem for a plate with functionally graded properties. ASME J. Appl. Mech. 1997, 64, 449–456. [Google Scholar] [CrossRef]
  41. Moës, N.; Dolbow, J.; Belytschko, T. A finite element method for crack growth without remeshing. Int. J. Numer. Meth. Eng. 1999, 46, 131–150. [Google Scholar] [CrossRef]
Figure 1. An inhomogeneous isotropic body with an edge crack.
Figure 1. An inhomogeneous isotropic body with an edge crack.
Applsci 09 03581 g001
Figure 2. An extension of contour Γ and a donut-type domain of integration A .
Figure 2. An extension of contour Γ and a donut-type domain of integration A .
Applsci 09 03581 g002
Figure 3. (a) Non-convex NEM grid N E M , (b) Laplace interpolation function ϕ J ( x ) at the crack tip.
Figure 3. (a) Non-convex NEM grid N E M , (b) Laplace interpolation function ϕ J ( x ) at the crack tip.
Applsci 09 03581 g003
Figure 4. A crack tip-originated circle and a graded domain of integration A .
Figure 4. A crack tip-originated circle and a graded domain of integration A .
Applsci 09 03581 g004
Figure 5. (a) Infinite plate with collinear cracks in plane strain state [2], (b) periodic model with double edge cracks and a gradient natural element method (NEM) grid (unit: m, N = 2516 ).
Figure 5. (a) Infinite plate with collinear cracks in plane strain state [2], (b) periodic model with double edge cracks and a gradient natural element method (NEM) grid (unit: m, N = 2516 ).
Applsci 09 03581 g005
Figure 6. Comparison of the stress distributions along the path Γ ( r / a = 0.0139 ) : (a) σ x x , (b) σ y y , and (c) τ x y .
Figure 6. Comparison of the stress distributions along the path Γ ( r / a = 0.0139 ) : (a) σ x x , (b) σ y y , and (c) τ x y .
Applsci 09 03581 g006
Figure 7. A cracked rectangular functionally graded material (FGM) plate with a varying Young’s modulus: (a) Geometry dimensions, (b) tension loading, (c) bending.
Figure 7. A cracked rectangular functionally graded material (FGM) plate with a varying Young’s modulus: (a) Geometry dimensions, (b) tension loading, (c) bending.
Applsci 09 03581 g007
Figure 8. (a) An FGM plate with an inclined edge crack under the exponential distribution load, (b) a uniform NEM grid and the domain of integration (235 nodes).
Figure 8. (a) An FGM plate with an inclined edge crack under the exponential distribution load, (b) a uniform NEM grid and the domain of integration (235 nodes).
Applsci 09 03581 g008
Table 1. Normalized stress intensity factors K I / σ y y π a for an edged cracked plate.
Table 1. Normalized stress intensity factors K I / σ y y π a for an edged cracked plate.
Relative   Crack   Length   a / W
Method E 2 / E 1 0.2 0.3 0.4 0.5
Uniform tensionProposed Method0.11.302 1.8432.5573.504
100.9601.1501.5812.176
Erdogan and Wu [40]0.11.297 (0.386)1.858 (−0.807)2.570 (−0.506)3.570 (−1.849)
101.002 (−4.192)1.229 (−6.428)1.588 (−0.441)2.176 (0.000)
Kim and Paulino
[ J k * ] [15]
0.11.284 (1.402)1.846 (−0.163)2.544 (0.118)3.496 (0.229)
101.003 (−4.287)1.228 (−6.352)1.588 (−0.441)2.175 (0.046)
Rao & Rahman
[ M ˜ ( 1 , 2 ) ] [34]
0.11.337 (−2.618)1.898 (−2.898)2.594 (−1.426)3.547 (−1.212)
100.996 (−3.614)1.234 (−6.807)1.598 (−1.064)2.189 (−0.594)
BendingProposed Method0.11.8851.8721.9082.141
100.5830.6360.8121.069
Erdogan and Wu [40]0.11.904 (−0.472)1.886 (−0.742)1.978 (−3.539)2.215 (−3.341)
100.565 (3.186)0.659 (−3.490)0.804(0.995)1.035 (3.285)
Kim and Paulino
[ J k * ] [15]
0.11.888 (−0.159)1.864 (0.429)1.943 (−0.035)2.145 (−0.140)
100.565 (3.186)0.659 (−3.490)0.804 (1.801)1.035 (3.285)
Rao & Rahman
[ M ˜ ( 1 , 2 ) ] [34]
0.11.903 (−0.420)1.875 (−0.160)1.954 (−2.354)2.155 (−0.650)
100.564 (3.369)0.664 (−4.217)0.812 (0.000)1.045 (2.297)
Table 2. Comparison of K I / σ y y π a between with and without enrichment.
Table 2. Comparison of K I / σ y y π a between with and without enrichment.
Relative   Crack   Length   a / W
Method E 2 / E 1 0.2 0.3 0.4 0.5
Uniform tensionWith enrichment0.11.3023 (46.902)1.8429 (29.041)2.5722 (37.785)3.5643 (27.587)
100.9960 (45.050)1.2225 (47.706)1.5210 (47.843)2.1473 (51.041)
Without enrichment0.10.69151.30771.60032.5810
100.54730.63930.79331.0513
BendingWith enrichment0.11.8847 (−24.683)1.8719 (−12.132)1.9083 (−4.910)2.1406 (−2.406)
100.5826(18.641)0.6361(7.436)0.8122 (9.554)1.0686 (15.319)
Without enrichment0.12.34992.09902.00202.1921
100.47400.58880.73460.9049
Table 3. Normalized stress intensity factors (SIFs) for an inclined crack in a functionally graded plate.
Table 3. Normalized stress intensity factors (SIFs) for an inclined crack in a functionally graded plate.
η a Present MethodKim and Paulino [ J k * ] [15]Rao and Rahman [ M ˜ ( 1 , 2 ) ] [34]
K I ε ¯ E ¯ π a K I I ε ¯ E ¯ π a K I ε ¯ E ¯ π a K I I ε ¯ E ¯ π a K I ε ¯ E ¯ π a K I I ε ¯ E ¯ π a
0.01.449 [0.951]0.606 [0.354]1.451 (−0.138)0.604 (0.331)1.448 (0.069)0.610 (−0.656)
0.11.377 [0.831]0.589 [0.335]1.396 (−1.361)0.579 (1.727)1.391 (−1.006)0.585 (0.684)
0.251.277 [0.673]0.525 [0.266]1.316 (−2.964)0.544 (−3.493)1.312 (−2.667)0.549 (−4.372)
0.51.163 [0.460]0.488 [0.241]1.196 (−2.759)0.491 (−0.611)1.190 (−2.269)0.495 (−1.414)
0.751.087 [0.357]0.434 [0.159]1.089 (−0.184)0.443 (−2.032)1.082 (0.462)0.446 (−2.691)
1.01.029 [0.240]0.411 [0.147]0.993 (3.625)0.402 (2.239)0.986 (4.361)0.404 (1.733)
[*] indicate the SIFs obtained without using the enrichment, while (*) indicate the relative differences (%).

Share and Cite

MDPI and ACS Style

Cho, J.-R. A Numerical Evaluation of SIFs of 2-D Functionally Graded Materials by Enriched Natural Element Method. Appl. Sci. 2019, 9, 3581. https://doi.org/10.3390/app9173581

AMA Style

Cho J-R. A Numerical Evaluation of SIFs of 2-D Functionally Graded Materials by Enriched Natural Element Method. Applied Sciences. 2019; 9(17):3581. https://doi.org/10.3390/app9173581

Chicago/Turabian Style

Cho, Jin-Rae. 2019. "A Numerical Evaluation of SIFs of 2-D Functionally Graded Materials by Enriched Natural Element Method" Applied Sciences 9, no. 17: 3581. https://doi.org/10.3390/app9173581

APA Style

Cho, J. -R. (2019). A Numerical Evaluation of SIFs of 2-D Functionally Graded Materials by Enriched Natural Element Method. Applied Sciences, 9(17), 3581. https://doi.org/10.3390/app9173581

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop