Functionally Graded Plate Fracture by Field Boundary Element

The paper describes the Field Boundary Element Method applied to the fracture analysis of a 2D rectangular plate made of Functionally Graded Material to calculate Mode I Stress Intensity Factor. The object of the Field Boundary Element Method is the transversely isotropic plane plate. Its material presents an exponential variation of the elasticity tensor depending on a scalar function of position, i.e., the elastic tensor results from multiplying a scalar function by a constant taken as a reference. Several examples using a parametric representation of the structural response show the suitability of the method that constitutes a sight of Stress Intensity Factor evaluation of Functionally Graded Materials plane plates even in the case of more complex geometries.


Introduction
In the recent past, a new class of composite materials, namely Functionally Graded Materials (FGM), have been attracted researchers and designed to satisfy Multiphysics requests, by combining different materials into one, resulting in continuous variable material. The main aspect of FGM with respect to traditional composites is that component properties are combined to get intermediate behavior between the constituents. When the FGM is used at the interface between two materials, it yields to continuous transition from one material to the other deleting the jump discontinuity of the response at the interface level. In other words, FGMs have been presented as an alternative to laminated composite materials that exhibit a discrepancy in properties at material interfaces. This material discontinuity in laminated composite materials leads to large interlaminar stresses and the possibility of crack initiation and propagation. Continuous transition greatly influences crack propagation at the interface level and avoids or reduces debonding and slip of layers [1,2,3,4].
Moreover, they are widely used in traditional laminate composites to overcome interface problems [5,6].
As a common example, we mention metal/ceramic which incorporates the excellent corrosion resistance and insulating properties of ceramics with the high strength and bonding ability of metals without requiring a high thermal stress rate. However, ceramics, in contrast to metals, have a fragile nature so crack-like defects are usually induced under service load conditions. From the above statement, FGMs are widely used in the research and manufacturing industry to optimize and design [7,8,9].
For the evaluation of cracks induced by several loading conditions, many theories and numerical applications have grown over the years and as well as many parameters to describe the fracture. The so-known SIF as the acronyms of Stress Intensity Factor in the cracked tip is an important parameter to determine the safety of a cracked structure [10,11].
As previously mentioned, the increasing interest for the cracks in the environment of Functionally Graded Materials (FGM) gave rise to a lot of methods as the Weight Function approach [12,13], Energy Release Method [14], and others.
The present work has aimed to apply the field boundary element method to heterogeneous material, in the case of simple scalar variability of properties, to model the graded bi-material interface.
The formulation enriches the Boundary Element Method with field terms that account for the variation of material elasticity [15,16,17,18]. Using the formulation, the application of Elastic Fracture Mechanics (EFM) is shown to calculate the Stress Intensity Factor under Mode I crack propagation.
The calculations follow the traditional formulation for EFM where it is supposed that around the crack tip the elasticity of the material is constant within the small region of stress concentration. The interface, the joints, or the surfaces, have great importance in graded material because the stress is usually limited to these locations. In these cases, the specific zone is the so-called "graded zone" and, if the volume fractions of the constituent materials change, by a certain smooth function, the material is defined as a Functionally Graded Material, FGM [19,20,21]. Furthermore, several benchmarks and reports calving models using Linear Elastic Fracture Mechanics (LEFM) are presented in [22]. Moreover, it is possible to evaluate other aspects of fractures in a microscale approach through anisotropic 2D models.
For cracked FGMs with general geometry and loading conditions, advanced numerical methods must be applied, because of the high mathematical complex nature of the governing partial differential equations with variable coefficients, and because the most available analytical methods can only be successfully applied to cracked FGMs with very simple geometry and loading conditions.
In this framework, we only cite the singular integral equation method [23,24,25,26,27,28], the classical finite element method (FEM) [29,30,31,32,33,34], the extended finite element method (XFEM) [32], the boundary integral equation method (BIEM) or boundary element method (BEM) [36,37,38,39,40,41,42,18]. This paper is organized as follows: In Section 2, the mathematical framework of the applied model, geometry, variation of mechanical properties, applied load, and fundamentals of integral equations for FGMs. Section 3, shows the result of FEM analysis and the present FBEM respectively, highlighting the matching between them. In section 4, all the results in terms of displacement, stress, and stress intensity factors (SIF) of mode I were discussed. Finally, some observations and further research developments are provided in Section 5.

Materials and Methods
Once the crack is formed, the crack behavior is determined by the direction of the forces or moments as depicted in Figure 1. Mode I is a tensile mode and Mode II is a sliding mode where the crack surfaces slip over each other. Finally, Mode III is a shear mode where this is primarily due to torsion moments.

Figure 1 Crack Modes
This section aims to recall the fundamentals of integral equations applied to FGM through a mixed field and boundary approach. After a brief description of the mechanical behavior of FGM within the framework of continuum mechanics, the integral formulation of equilibrium is described. The proof originates from the Somigliana identity specifically to materials for which there is no analytic fundamental solution. Betti's reciprocity theorem is transformed using a sample material whose fundamental equation is known. The resultant equation includes, in addition to the ordinary boundary variables, internal variables that cannot be simplified. This yields a combined field and boundary integral equation method.
Elastic constitutive law of the material is governed by fourth-order elastic tensor that, in a Cartesian coordinate system, has the form: where ( ) is the Cauchy stress and = , , is the infinitesimal strain.
Constitutive tensor has the usual symmetries: The structure, covering volume V, is subjected to body forces ( ) traction ( ) on the loaded boundary , and prescribed displacements on = \ Equilibrium equation of the structure, in terms of displacement, ( ) is: Let us consider * ( , ) the displacement at a point x of the elastic infinite space with constant elastic tensor , due to a one-point force acting at the point y in Fig.2.
where ( ) is a numerical coefficient depending on the boundary regularity at y. Moreover, eq. (4) can be simplified by considering the actual form of ( ) tensor. For scalar graded materials it depends on a function ( ) and on a reference constant array . Here we assume that Kelvin's solution is described by the same tensor , then eq. (4) becomes [21,41]: In particular, ( ) and ( ) are the displacements and tractions on the boundary S. The terms * ( , ) and * ( , ) represent the displacements and tractions of Kelvin's solution. Numerical solution of eq. (5) gives the displacement field of the structure, which is used to calculate the displacement at the crack tip. The proposed formulation has been applied to symmetrical plates where the crack has been modeled using a single line, so no dual equation is needed, nor do hypersingular integrals occur. Another kind of strategy for crack modeling would consist of multi-region discretization that avoids hypersingular integration as well [21]. In fact, in the present paper, to overcome the hypersingular integration, a symmetric quarter plate geometry, as depicted in Figure 3, has been considered varying material properties and increasing crack length also in the case of where ( ) is the Mode I displacement of a nodal point and is the distance from the crack tip. Coefficients and are coefficients to be determined by minimizing the difference between (6) and numerical solution given by (5). The numerical value of SIF is calculated assuming, that within a small neighborhood of the crack tip, the elastic modulus of the material is and is the Poisson ratio both assumed to be constant. By these considerations the SIF has the usual displacement dependency form: The described method has been applied to a quarter of a rectangular plate due to symmetry. The plate has length and height ℎ = 2 , with a crack at the center and aligned to the short plate side. The crack length is .
shows the mesh and boundary conditions of the analyzed structure. The plate has null displacement on the nodes with y=0, outside the crack. The x=0 side has symmetry constraints. The structure is loaded by a traction = 1 along y-direction on the top side. Three different crack lengths have been investigated, namely = 0.11; = 0.33, and = 0.55. The material property variation has been ascribed to the Elastic Modulus only, i.e. the Poisson ratio has been fixed to = 0.3. In particular, the elastic tensor varies through a multiplicative scalar field ( , ), and fulfilling the following equation: In the proposed example, only the case of variability along with x and y-directions have been considered. The reference of the coordinate system has the origin coincident with the center of the crack, consequently, the analyzed specimen has the origin at the leftmost side. Two variabilities have been accounted: 1.
: Young Modulus decreasing in the y-direction with ( , ) = ( ) spanning from 8 to 1: 4. : Young Modulus increasing in the y-direction with ( , ) = ( ) spanning from 1 to 8: The numerical calculation has been performed assuming plane strain condition with reference Young Modulus equal to = 10 , Poisson ratio = 0.3 , and crack length variation equal to = 0.11, 0.33, 0.55. Figure 4 and Figure 5 show the variation of Young Modulus along x-direction and ydirection. Furthermore, the dark blue and light blue curves indicate the decreasing Young Modulus ( ) and increasing Young Modulus ( ) respectively. In the same manner, dark green and light green curves indicate the increasing Young Modulus ( ) and decreasing Young Modulus ( ) respectively.

Results
In this section, the solution of the FBEM algorithm has been discussed. The Field Boundary Element calculation has been compared with a commercial FEM code result (Ansys ©). In the following, Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 21 July 2021 doi:10.20944/preprints202107.0487.v1 Figure 8, Figure 8 the displacement components along the x-direction and the y-direction have been shown for the three different crack lengths and together with the FBEM results. It can be seen that FBEM and FEM produce the same displacements. In Figure 10, Figure 11, Figure 11 the stress components along the x-direction and along the y-direction were compared too.    To investigate the goodness of the implemented FBEM results in comparison to the FEM modeling results, numerical calculations with different variations of Young's Modulus were also performed and reported in Table 1 regarding KI fracture parameter only: Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 21 July 2021 doi:10.20944/preprints202107.0487.v1 For completeness, summary plots of the Stress Intensity Factor KI for all Young Modulus variations and the comparison between FEM and FBEM for all normalized crack lengths were shown in Figure 12

Discussion
As depicted in the previous sections, it has been studied an isotropic 2D case where the comparison of the differences between the SIF in FEM versus other literature equations is negligible. Conversely, in the case of increasing anisotropy ratio and constant Poisson's ratio, differences increase as in the case of metals where ≅ 0.3. In this paper, it has been analyzed the fracture on a large-scale simulation, but as previously introduced, the majority of materials are anisotropic, then is recommended to investigate the same model with anisotropic behavior where the crack propagation is not homogenized and the thickness influences the results getting closer to the center where the crack force is maximum. Moreover, the actual method resulted to be versatile because it is possible to use the same approach, with few changes, to heterogeneous and anisotropic ones, even for more complex geometry.
From the above figures of displacements and stresses, it can be seen that the present FBEM formulation fits well with the standard FEM. In particular, it can be seen that the present FBEM method doesn't need mesh refinement to appreciate the convergence that Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 21 July 2021 doi:10.20944/preprints202107.0487.v1 can be seen even using a sparse mesh allowing to obtain good results without heavy computational time and cost. In Table 1, is listed as a function of . The results are also summarized in Figure   12, Figure 13, Figure 14, where the Stress Intensity Factor, , is plotted for all the computed numerical cases. It can be seen that for increasing , the SIF increases even in case of increasing . When approaches 1, the plate tends to be homogeneous. Conversely when decreases, the stiffness of the crack tip decreases and so also does the calculated SIF.
In Figure 12 for = 0.11 it is possible to see a pseudo-convenience in the sense of stress to make an FGM plate with an increasing material gradation along the x-direction with a SIF of 240.3 % lower than the homogeneous case.
In Figure 13 and Figure 14, where = 0.33 and = 0.55 the pseudo-convenience in manufacturing an FGM plate with an increasing material gradation along the x-direction with SIF values around 152.2% and 86.3% lower than the homogeneous cases as well.
The analysis highlight that the stress concentration is particularly effective where the material has increasing stiffness. It has to be expected that cracks start from the stiffer side toward the softer one. It should be considered that in common FGM, where the stiffer side is metallic and the softer one is ceramic, the ductility of the system plays a favorable role concerning the growth of the crack.

Conclusions
The numerical calculation shows the feasibility of the method that allows the application of Integral Equations to fracture mechanics even in the cases of heterogeneous materials; moreover, it is shown that the variability of the mechanical properties of the structure affects the value of the SIF. The procedure proposed here can be used to design composite structures in which the variation of constitutive parameters is prescribed to minimize the stress intensity around the crack tip. The method could also be used to evaluate crack path propagation since it is greatly influenced by the variability of the material near the crack tip. It should be considered that although many authors assume a constant behavior of the material around the crack tip, the influence of the variability can be studied in detail using the integral formulation which can be realized with a slight modification of the present FBEM.
The present method showed with good approximation, matching the standard Finite Element Method, even in the case of Young's Modulus variation in both x and y directions, predicting crack propagation and manufacturing convenience also due to the lower rank of computational costs compared to FEM.
Finally, the formulation can be extended to anisotropic variable material, the application to anisotropic or Mode II material gradient dependence of FGM through present FBEM formulation will be the subject of further research activity.