Three-Dimensional Fracture Analysis in Functionally Graded Materials Using the Finite Block Method in Strong Form

In this paper, the application of the strong-form finite block method (FBM) to three-dimensional fracture analysis with functionally graded materials is presented. The main idea of the strong-form FBM is that it transforms the arbitrary physical domain into a normalized domain and utilizes the direct collocation method to form a linear system. Using the mapping technique, partial differential matrices of any order can be constructed directly. Frameworks of the strong-form FBM for three-dimensional problems based on Lagrange polynomial interpolation and Chebyshev polynomial interpolation were developed. As the dominant parameters in linear elastic fracture mechanics, the stress intensity factors with functionally graded materials (FGMs) were determined according to the crack opening displacement criteria. Several numerical examples are presented using a few blocks to demonstrate the accuracy and efficiency of the strong-form FBM.


Introduction
In engineering structures such as reinforced concrete column steel beams, landing gear, and aircraft skin, the presence of cracks has a significant impact on their health.To improve the reliability of engineering structures and gain a deeper understanding of the mechanisms of cracks, many researchers have studied fracture mechanics in various fields, such as civil engineering [1][2][3], aerospace [4][5][6], and materials science engineering [7][8][9][10].In fracture mechanics, the asymptotic stress field surrounding a crack tip is a critical research issue.And, for linear elasticity, stress intensity factors (SIFs) are essential parameters, characterizing the strength of a stress field in the vicinity of the crack tip.Also, the SIFs are the dominant terms in the William's series expansion.Therefore, the primary object of linear fracture mechanics is to determine the corresponding SIFs.Up to now, many theoretical or numerical methods have been developed for determining SIFs [11][12][13][14][15].However, due to the difficulty of the mathematics involved, the number of models with which SIFs can be obtained analytically is limited, especially for FGMs.With the development of computers, numerical methods have become the most commonly used to solve general crack problems.The well-known numerical methods are the mesh-dependent methods and the meshless methods.
Regarding the mesh-dependent methods, the representative methods are the finite element method (FEM) [16,17] and the boundary element method (BEM) [18,19].The FEM discretizes the analyzed domain using small elements with different shapes, including triangular elements and quadrilateral elements for two-dimensional problems and tetrahedral elements and hexahedral elements for three-dimensional problems, granting the FEM very strong adaptability for complicate boundaries.However, for crack problems, the meshes near the crack tip have to be refined severely in order to obtain reliable results, which is usually very time consuming.Unlike the FEM, the BEM only divides the meshes on the boundary.And due to dimensional reduction and its semi-analytical character, it can obtain high accuracy while also providing timesaving properties.However, the fundamental solutions of the BEM are not easy to obtain in general problems.Furthermore, for FGMs, the boundary integral equation always needs an extra domain integral term, which robs the BEM of the advantage of dividing meshes only on the boundary.Although many methods have been proposed to transform the domain integral into a boundary integral, such as the double reciprocity method [20] and the radial integration method [21], additionally errors were introduced into the BEM.In contrast to the mesh-dependent methods, in recent decades, the meshless methods have gained much attention due to their advantage of not requiring mesh division.In meshless methods, a set of scattered nodes are distributed to discretize the physical domain and further form a linear system through corresponding kernel functions.The typically used meshless methods are the radial basis function collocation method (RBFCM) [22][23][24][25], the finite integration method (FIM) [26,27], the finite block method (FBM) [28], the method of fundamental solution (MFS) [29][30][31], the local radial basis function collocation method (LRBFCM) [32][33][34], and so on.Up to now, those meshless methods have been successfully applied in various problems, such as band structure computation in phononic crystals [35,36], elastic wave propagation [37], the simulation of water pollution [38], geometric modeling [39,40], etc.
Among the meshless methods, the FBM is a more suitable strong-form method for solving crack problems.The FBM was first proposed by Wen et al. [28] to solve elastic problems with FGMs.The main operation involved in the FBM is to transform the physical domain into the normalized domain and form the linear system by using mapping technique and a collocation method in the normalized domain, allowing the form of the node distribution in the physical domain to be neglected.And any-order partial differential matrices can be obtained by using a first-order partial differential matrix.Furthermore, due to the continuity of stresses and displacements between two adjacent blocks, the FBM can achieve higher accuracy with fewer blocks than other collocation methods.Therefore, it has attracted a great deal of research interest in recent years.Afterwards, the FBM was successfully applied to heat conduction problems [41], contact problems [42], non-linear elasticity [43], and fracture mechanics [44].Also, many novel numerical methods have been proposed based on the idea of the FBM, such as the finite-and infinite-block Petrov-Galerkin method [45] and the element differential method [46][47][48].The associated studies all demonstrated the superiority of the FBM.In the current paper, the strong-form FBM is further extended to analyze three-dimensional crack problems.The frameworks of the FBM based on Lagrange polynomial interpolation (FBML) and Chebyshev polynomial interpolation (FBMC) are established.The FEM solutions were selected as the benchmark and obtained using ABAQUS with the subroutine UMAT developed for FGMs.In ABAQUS, various utility subroutines are available in Abaqus/Standard, including UMAT.The UMAT user subroutine interface in ABAQUS is very powerful and efficient in numerical simulation in solid mechanics coded in Fortran and C++.The subroutine UMAT can be used to define the mechanical constitutive behavior of a material and can be used with any procedure that includes mechanical behavior.In addition, users can use solution-dependent state variables and update the stresses and solution-dependent state variables with respect to their values at the end of the increment for which they are called.All the corresponding details are described in the ABAQUS documentation [49].The crack opening displacement (COD) criteria were used to evaluate the SIFs with FGMs.
The structure of this paper is organized as follows.In Section 2, the mapping technique for three-dimensional problems is introduced briefly.In Section 3, the frameworks of the strong-form FBM are established.In Section 4, a description of three-dimensional elasticity problems is presented.In Section 5, several numerical examples are given to demonstrate the robustness of the strong-form FBM.Some conclusions are drawn in Section 6.

Block Mapping in Three-Dimensional Problems
For three-dimensional problems, a block in a physical domain can be mapped into a normalized domain with 20 seeds, as shown in Figure 1.The shape functions are expressed as The coordinates in the physical domain can be written as in which (x k , y k , z k ) is the coordinate of kth seed.The first-order partial differential operators in Cartesian coordinates of real space are where the coefficients are From Equation (3), the second-order partial differential operators can be obtained by using the derivative rule, as follows: 13 33 (11) Figure 1.Three-dimensional normalized domain and its seeds.

Strong-Form Finite Block Method
The strong-form numerical method is also called the collocation method.In the strong-form finite block method, the direct collocation method is used to form a differential matrix in the normalized domain.Using the mapping technique, a differential matrix in the physical domain can then be constructed.In the following section, the approximation of the strong-form finite block method is introduced.

A Brief Description of the FBML
In the normalized domain shown in Figure 1, the smooth function ( ) u ξ,η,ζ can be approximated using Lagrange polynomial as in which ( ) is the shape function, N is the number of collocation nodes along axis x (y or z), and ( ) , q u denotes the nodal value.The collocation nodes i ξ , j η , and k ζ are selected as uniform nodes ( ) or roots of the Chebyshev polynomial of the first kind ( ) and the polynomial functions are defined as follows:

Strong-Form Finite Block Method
The strong-form numerical method is also called the collocation method.In the strongform finite block method, the direct collocation method is used to form a differential matrix in the normalized domain.Using the mapping technique, a differential matrix in the physical domain can then be constructed.In the following section, the approximation of the strong-form finite block method is introduced.

A Brief Description of the FBML
In the normalized domain shown in Figure 1, the smooth function u(ξ, η, ζ) can be approximated using Lagrange polynomial as in which Φ(ξ, η, ζ) is the shape function, N is the number of collocation nodes along axis x (y or z), and q = (k − 1)N 2 + (j − 1)N + i, u q denotes the nodal value.The collocation nodes ξ i , η j , and ζ k are selected as uniform nodes or roots of the Chebyshev polynomial of the first kind and the polynomial functions are defined as follows: Then, the first-order partial derivatives of shape function Φ with respect to ξ, η, and ζ can be determined straightforwardly as follows ∂G η, η j ∂η Substituting all the collocation nodes into partial derivatives of the shape function, the matrix form of the first-order partial derivatives of u can be formulated as in which u = [u 1 , u 2 , . . . ,u N 3 ] is the vector of nodal values with sizes of N 3 × 1, while D ξ , D η , and D ζ are all the first-order differential matrices with dimensions of N 3 × N 3 .Thus, we have where D 0 is the first-order differential matrix for a one-dimensional case with dimensions of N × N, which can be obtained from Equation (19).By using the transformation matrix T η and T ζ , one can obtain in which For any-order partial derivatives with respect to ξ, η, and ζ, we have Therefore, from Equation ( 3), one can obtain where D x , D y , and D z constitute the first-order differential matrix in the physical domain, and B γµ denotes the diagonal matrix, defined as Similarly, any-order partial derivatives in the physical domain with respect to x, y, and z can be expressed as

The Brief of the FBMC
Besides Lagrange polynomial interpolation, the framework of the FBM can also be constructed using Chebyshev polynomial interpolation.The unknown function u(ξ, η, ζ) can be written in the following form where α ijk is the unknown coefficient and T i (ξ), T j (η), and T k (ζ) are i-th, j-th, and k-th order Chebyshev polynomials of the first kind, respectively.According to the properties of the Chebyshev polynomial, the first-order and second-order partial derivatives of u with respect to ξ, η, and ζ can be obtained analytically as follows in which U p (w), w = ξ, η, ζ, p = i, j, k, denotes the p-th order Chebyshev polynomial of the second kind.Similar to the Lagrange polynomial interpolation, by substituting all uniform nodes or roots of the Chebyshev polynomial of the first kind into Equations ( 34)-( 42), the matrix form of different orders of partial derivatives is expressed as where α = [α 1 , α 2 , . . . ,α N 3 ] T is the vector of unknown coefficients with dimensions of N 3 × 1, while D κ and D κχ , for which κ, χ = ξ, η, ζ, are differential matrices in the normalized domain with dimensions of N 3 × N 3 .From Equations ( 3)-( 11), the matrix form of the first-order partial derivatives in the physical domain can be formulated as (48) in which Θ x , Θ y , and Θ z are the first-order differential matrices.For second-order partial derivatives, we have in which the coefficient matrices C (bd) t , t = 1, 2, . . ., 9, are presented in Appendix A.

Three-Dimensional Elasticity 4.1. Governing Equations
For 3D static elasticity, the governing equations of the FGMs are given as where σ xx , σ yy , σ zz , σ xy , σ xz , and σ yz denote stress components, while f x , f y , and f z are the body force components, respectively.The stress components are defined in the following way in which u x , u y , and u z indicate the displacement components; λ = Eν/(1 + ν)/(1 − 2ν), G = E/2/(1 + ν), and E = E(x, y, z) denote the Young's modulus with material gradation; and ν indicates the Poisson's ratio.The governing equation can be further formulated as where

Boundary Conditions
The displacement condition and traction condition are defined on boundaries Γ 1 and Γ 2 as where u x , u y , and u z are given displacements; n = n x , n y , n z is the outward unit normal vector; and Γ 1 and Γ 2 indicate the Dirichlet boundary and the Neumann boundary, respectively.

Numerical Discretization
In this section, the governing equations of FGMs are discretized using the FBML and the FBMC.For the FBML, by substituting N 3 collocation nodes into partial derivatives, Equations ( 57)-( 59) can be formulated in matrix form as where u d = [u d1 , u d2 , . . . ,u dN 3 ] T , d = x, y, z, denotes the displacement component vector with dimensions of N 3 × 1; f d = [ f d1 , f d2 , . . . ,f dN 3 ] T is the body force component vector with dimensions of N 3 × 1; and k γ , λ d , and G d denote the diagonal coefficient matrix, expressed as The displacement components can be obtained directly by solving the linear system of Equations ( 67)-(69).For the FBMC, the matrix form of the governing equations can be expressed as y, z, denotes the unknown coefficient vector.Once the coefficient vectors are determined, the displacement components can be acquired from Equation (33).

Numerical Examples and Discussion
According to the COD criteria, once the displacement components on the crack surface are determined using the linear algebraic Equations ( 67)-(69) or Equations ( 73)-( 75), the approximated SIFs can be evaluated directly.The relationships between the SIFs and the displacement components near the crack-tip are defined as where K I , K II , and K III are the mode I, mode II, and mode III stress intensity factors, and r is the small distance from the crack front to the node on the crack surface, respectively.∆u x , ∆u y , and ∆u z are the CODs in the local coordinate system, which is located at the front of the crack.

Central-Crack Plate under Tension Load
In this example, consider a central-crack plate for which T/W = 1 and H/W = 1 subjected to unit uniform loads σ 0 in the z-direction on the top boundary, as shown in Figure 2a.

Numerical Examples and Discussion
According to the COD criteria, once the displacement components on the crack surface are determined using the linear algebraic Equations ( 67)-(69) or Equations ( 73)-( 75), the approximated SIFs can be evaluated directly.The relationships between the SIFs and the displacement components near the crack-tip are defined as where K I , K II , and K III are the mode I, mode II, and mode III stress intensity factors, and r is the small distance from the crack front to the node on the crack surface, respectively.
x u ' Δ , y u ' Δ , and z u ' Δ are the CODs in the local coordinate system, which is located at the front of the crack.

Central-Crack Plate under Tension Load
In this example, consider a central-crack plate for which T W = 1 and H W = 1 subjected to unit uniform loads σ 0 in the z-direction on the top boundary, as shown in Figure 2a.The rigid constraint is imposed on the bottom boundary.The crack length is 2a, and the Poisson's ratio was set as 0.3.For the normalized Young's modulus E, consider the following three cases: Case 1-Homogeneous material: The rigid constraint is imposed on the bottom boundary.The crack length is 2a, and the Poisson's ratio was set as 0.3.For the normalized Young's modulus E, consider the following three cases: Case 1-Homogeneous material: E = 1.
Case 2-Material gradation in the y-direction: E 2 and E 1 are the normalized Young's modulus values at y = T (or z = H) and y = 0 (or z = 0), respectively.The ratio of the Young's modulus is E 1 /E 2 = 2. Due to the symmetry, only half of the plate is analyzed in this example.In the numerical procedure, for the FBM, the computational domain is divided into four blocks, as depicted in Figure 2b.The 15 × 15 × 15 roots of the Chebyshev polynomial of the first kind are distributed in each block.For the FEM, 63,600 20-node quadratic hexahedral elements and 15-node wedge elements with a total 269,377 nodes are used, as shown in Figure 3, while a UMAT subroutine was developed for FGMs.By utilizing the COD criteria, the results of the comparison of the normalized mode I SIF with respect to a/W = 0.25, 0.5, 0.75 are given in Figure 4; in this section, we select the displacements of the third point in front of the crack tip for the SIF calculation in the FBM and FEM.Obviously, the numerical results of the FBML, the FBMC, and the FEM show very satisfactory agreement both in homogeneous material and FGM.In Figure 4, it can be observed that the SIFs of the homogeneous material are nearly uniform throughout most of the crack front.There are some differences in the SIFs at the two ends of the crack front and in the interior of the crack front due to the near-plain stress conditions.In addition, a difference in the SIF with a different number of collocation points was observed.This shows that convergent results pertaining to the SIF can be obtained when the number of collocation points along each side N > 11.Furthermore, the normalized mode I SIF becomes larger with the increase in crack length.A comparison of the results also implies that the FBML and the FBMC are almost equivalent.4; in this section, we select the displacements of the third po of the crack tip for the SIF calculation in the FBM and FEM.Obviously, th results of the FBML, the FBMC, and the FEM show very satisfactory agreem homogeneous material and FGM.In Figure 4, it can be observed that the SIF mogeneous material are nearly uniform throughout most of the crack fron some differences in the SIFs at the two ends of the crack front and in the int crack front due to the near-plain stress conditions.In addition, a differenc with a different number of collocation points was observed.This shows that results pertaining to the SIF can be obtained when the number of collocation p each side > N 11 .Furthermore, the normalized mode I SIF becomes larger w crease in crack length.A comparison of the results also implies that the FBM FBMC are almost equivalent.

Edge-Crack Plate under Tension Load
An edge-crack plate with a fixed bottom subjected to unit uniform tensile loads σ 0 on the top boundary is discussed in the current example, as shown in Figure 5a.The scale parameters of the plate are T W = 2 and H W = 2 .The crack length is a = 0.5, and the cracks form an angle θ with the x axis.The Poisson's ratio is 0.3.Herein, the following three different situations are still considered.

Edge-Crack Plate under Tension Load
An edge-crack plate with a fixed bottom subjected to unit uniform tensile loads σ 0 on the top boundary is discussed in the current example, as shown in Figure 5a.The scale parameters of the plate are T/W = 2 and H/W = 2.The crack length is a = 0.5, and the cracks form an angle θ with the x axis.The Poisson's ratio is 0.3.Herein, the following three different situations are still considered.In the computational process, for the FBM, the computational domain is divided into four blocks, as depicted in Figure 5b.The × × 15 15 15 roots of the Chebyshev polynomial of the first kind are distributed in each block.In ABAQUS with the UMAT subroutine, 67,360 20-node quadratic hexahedral elements and 15-node wedge elements with a total of 284,699 nodes were employed, as shown in Figure 6.Hence, we adopted the displacements of the third collocation node in front of the crack front to carry out an evaluation of the COD criterion; the comparison of the results of the normalized mode I and mode II SIFs with respect to    θ = , , 0 30 45 are shown in Figures 7-9.It is clear that the normalized mixed-mode SIFs of the FBML, the FBMC, and the FEM yield excellent consistencies in both homogeneous material and FGM.
In Figure 7, it can be observed that the SIFs of the homogeneous material have a nearly uniform distribution along most of the crack front with respect to Actually, the mode II SIFs are nearly zero due to the mode I edge-crack when  θ = 0 .And when ≠  θ 0 , this edge-crack model transforms into a mixed-mode problem.The influence of material gradation along the z-direction is very small due to the isotropy on the crack surface.Furthermore, the tilt angle θ of the crack surface directly affects the distribution of mixed-mode SIFs.Similar to the central-crack, the normalized SIFs become larger with the increase in the tilt angle.The comparison of the results suggests that the FBML and the FBMC can achieve very high accuracy in three-dimensional edge-crack analysis.
E 2 and E 1 represent the normalized Young's modulus at y = T (or z = H) and y = 0 (or z = 0), and assume that E 1 /E 2 = 2.The local coordinate system x y z located on the crack surface is used to evaluate the CODs, namely, ∆u x , ∆u y , and ∆u z , as depicted in Figure 5b.
In the computational process, for the FBM, the computational domain is divided into four blocks, as depicted in Figure 5b.The 15 × 15 × 15 roots of the Chebyshev polynomial of the first kind are distributed in each block.In ABAQUS with the UMAT subroutine, 67,360 20-node quadratic hexahedral elements and 15-node wedge elements with a total of 284,699 nodes were employed, as shown in Figure 6.Hence, we adopted the displacements of the third collocation node in front of the crack front to carry out an evaluation of the COD criterion; the comparison of the results of the normalized mode I and mode II SIFs with respect to θ = 0 • , 30 • , 45 • are shown in Figures 7-9.It is clear that the normalized mixed-mode SIFs of the FBML, the FBMC, and the FEM yield excellent consistencies in both homogeneous material and FGM.
In Figure 7, it can be observed that the SIFs of the homogeneous material have a nearly uniform distribution along most of the crack front with respect to θ = 0 • , 30 • , 45 • .Actually, the mode II SIFs are nearly zero due to the mode I edge-crack when θ = 0 • .And when θ = 0 • , this edge-crack model transforms into a mixed-mode problem.The influence of material gradation along the z-direction is very small due to the isotropy on the crack surface.Furthermore, the tilt angle θ of the crack surface directly affects the distribution of mixed-mode SIFs.Similar to the central-crack, the normalized SIFs become larger with the increase in the tilt angle.The comparison of the results suggests that the FBML and the FBMC can achieve very high accuracy in three-dimensional edge-crack analysis.

Curved Edge-Crack under Shearing Load
In this example, consider a curved edge-crack plate with a fixed bottom subjected to unit uniform shearing load τ 0 on the top boundary, which is shown in Figure 10a.The crack front is a curved line whose geometric parameters are shown in Figure 10b.The front end of the crack surface is an arc of radius r with its center at point O, where and O is the midpoint of line segment AB, a is the crack length.And we have W T = 2 and W H =1 .The Poisson's ratio adopted is 0.3.Herein, the material gradation in the y-direction is considered, and the normalized Young's modulus is In the numerical computation, for the FBM, the computational domain is divided into four blocks, as depicted in Figure 10c.The × × 15 15 15 roots of the Chebyshev polynomial of the first kind are used in each block.In ABAQUS with the UMAT subroutine, as shown in Figure 11, 20-node quadratic hexahedral elements and 15-node wedge elements are employed, with a total number of 30,080 elements and 129,523 nodes.Then, the displacements of the third point in front of the crack tip are selected to calculate the

Curved Edge-Crack under Shearing Load
In this example, consider a curved edge-crack plate with a fixed bottom subjected to unit uniform shearing load τ 0 on the top boundary, which is shown in Figure 10a.The crack front is a curved line whose geometric parameters are shown in Figure 10b.The front end of the crack surface is an arc of radius r with its center at point O, where r = (0.5T) 2 + a 2 and O is the midpoint of line segment AB, a is the crack length.And we have T/W = 2 and H/W = 1.The Poisson's ratio adopted is 0.3.Herein, the material gradation in the y-direction is considered, and the normalized Young's modulus is E = E 1 e y/Tln(E 2 /E 1 ) , in which E 1 /E 2 = 2. E 2 and E 1 are the normalized Young's modulus values at y = T and y = 0, respectively.

Curved Edge-Crack under Shearing Load
In this example, consider a curved edge-crack plate with a fixed bottom subjected to unit uniform shearing load τ 0 on the top boundary, which is shown in Figure 10a.The crack front is a curved line whose geometric parameters are shown in Figure 10b.The front end of the crack surface is an arc of radius r with its center at point O, where and O is the midpoint of line segment AB, a is the crack length.And we have W T = 2 and W H =1 .The Poisson's ratio adopted is 0.3.Herein, the material gradation in the y-direction is considered, and the normalized Young's modulus is In the numerical computation, for the FBM, the computational domain is divided into four blocks, as depicted in Figure 10c.The × × 15 15 15 roots of the Chebyshev polynomial of the first kind are used in each block.In ABAQUS with the UMAT subroutine, as shown in Figure 11, 20-node quadratic hexahedral elements and 15-node wedge elements are employed, with a total number of 30,080 elements and 129,523 nodes.Then,  In the numerical computation, for the FBM, the computational domain is divided into four blocks, as depicted in Figure 10c.The 15 × 15 × 15 roots of the Chebyshev polynomial of the first kind are used in each block.In ABAQUS with the UMAT subroutine, as shown in Figure 11, 20-node quadratic hexahedral elements and 15-node wedge elements are employed, with a total number of 30,080 elements and 129,523 nodes.Then, the displacements of the third point in front of the crack tip are selected to calculate the SIFs; the comparison of the results of the normalized mode I, mode II, and mode III SIFs with respect to a/W = 0.25, 0.5, 0.75 are depicted in Figure 12.There is little difference in the numerical results regarding the SIFs of the FBML, the FBMC and the FEM.This further validates that the FBML and the FBMC are equivalent numerically.
SIFs; the comparison of the results of the normalized mode I, mode II, and with respect to a / W = 0.25,0.5,0.75 are depicted in Figure 12.There is l in the numerical results regarding the SIFs of the FBML, the FBMC and further validates that the FBML and the FBMC are equivalent numerically.
In Figure 12, it can be observed that the related error among the so FBML, the FBMC, and the FEM is approximately less than 10%.Comp straight crack front, the fluctuation of the SIFs of the curved crack is la mode I and mode II SIFs possess a zero-point due to the distribution of dis addition, the mode III SIFs increase with the increase in crack length comparisons of the SIFs presented above, the FBML and the FBMC can ac curacy in the analysis of three-dimensional mixed-mode fracture mechanic SIFs; the comparison of the results of the normalized mode I, mode II, and mode III SIFs with respect to a / W = 0.25,0.5,0.75 are depicted in Figure 12.There is little difference in the numerical results regarding the SIFs of the FBML, the FBMC and the FEM.This further validates that the FBML and the FBMC are equivalent numerically.
In Figure 12, it can be observed that the related error among the solutions of the FBML, the FBMC, and the FEM is approximately less than 10%.Compared with the straight crack front, the fluctuation of the SIFs of the curved crack is larger.And the mode I and mode II SIFs possess a zero-point due to the distribution of displacements.In addition, the mode III SIFs increase with the increase in crack length a W . Through comparisons of the SIFs presented above, the FBML and the FBMC can achieve high accuracy in the analysis of three-dimensional mixed-mode fracture mechanics.= in which the partial derivatives of β γµ , γ, µ = 1, 2, 3, can be obtained from Equation (4) by using derivative rule.

Figure 1 .
Figure 1.Three-dimensional normalized domain and its seeds.

Figure 2 .
Figure 2. Central-cracked plate: (a) geometry of the plate; (b) half of the plate divided by four blocks.

Figure 2 .
Figure 2. Central-cracked plate: (a) geometry of the plate; (b) half of the plate divided by four blocks.

erials 2023 ,
16, x FOR PEER REVIEW wedge elements with a total 269,377 nodes are used, as shown in Figure UMAT subroutine was developed for FGMs.By utilizing the COD criteria, th the comparison of the normalized mode I SIF with respect to a W = , 0.25 0 given in Figure

Figure 3 . 5 FEM
Figure 3.The FEM model: (a) mesh distribution of the half plate; (b) results regardin half of the plate.

Figure 3 .Figure 3 .Figure 4 .
Figure 3.The FEM model: (a) mesh distribution of the half plate; (b) results regarding u z of the half of the plate.

FEM 5 FEMFigure 4 .
Figure 4. Comparison of the results obtained for the FEM, FBML, and FBMC: (a) normalized mode I SIFs of homogeneous material; (b) normalized mode I SIFs of material gradation in y-direction; (c) normalized mode I SIFs of material gradation in z-direction.

Figure 5 .
Figure 5. Edge-crack plate: (a) geometry of the plate; (b) the whole plate divided by four blocks.

E 2 and 1 E 1 .
Figure 5b.In the computational process, for the FBM, the computational domain is divided into four blocks, as depicted in Figure5b.The × × 1515 15 roots of the Chebyshev polynomial of the first kind are distributed in each block.In ABAQUS with the UMAT subroutine, 67,360 20-node quadratic hexahedral elements and 15-node wedge elements with a total of 284,699 nodes were employed, as shown in Figure6.Hence, we adopted the displacements of the third collocation node in front of the crack front to carry out an evaluation of the COD criterion; the comparison of the results of the normalized mode I and mode II SIFs with respect to

Figure 6 .
Figure 6.The FEM model: (a) mesh distribution of the whole plate with a slant edg sults regarding z u of the whole plate.

Figure 6 .Figure 6 .Figure 7 .Figure 7 .
Figure 6.The FEM model: (a) mesh distribution of the whole plate with a slant edge crack; (b) results regarding u z of the whole plate.

1 . E 2 and 1 EFigure 10 .
Figure 10.Curved edge-crack plate: (a) geometry of the plate; (b) geometry of crack front; (c) the whole plate divided by four blocks.

Figure 9 .
Figure 9. Result comparisons for FGM with material gradation in z-direction: (a) normalized mode I SIFs; (b) normalized mode II SIFs.

Figure 9 .
Figure 9. Result comparisons for FGM with material gradation in z-direction: (a) normalized mode I SIFs; (b) normalized mode II SIFs.

1 . E 2 and 1 EFigure 10 .
Figure 10.Curved edge-crack plate: (a) geometry of the plate; (b) geometry of crack front; (c) the whole plate divided by four blocks.

Figure 10 .
Figure 10.Curved edge-crack plate: (a) geometry of the plate; (b) geometry of crack front; (c) the whole plate divided by four blocks.

Figure 11 .FEMFigure 11 .
Figure 11.The FEM model: (a) mesh distribution of the whole plate with a curved results regarding z u of the whole plate.

Figure 11 .
Figure 11.The FEM model: (a) mesh distribution of the whole plate with a curved edge crack; (b) results regarding z u of the whole plate.