Adaptive phase field modelling of crack propagation in orthotropic functionally graded materials

Hirshikesh .Emilio Mrtínez-P~ned .Sundrrjn Ntrjn ,*

a Department of Mechanical Engineering.Indian Institute of Technology Madras.Chennai.600036.India

b Department of Civil and Environmental Engineering.Imperial College London.London.SW7 2AZ.UK

Keywords: Functionally graded materials Phase field fracture Polygonal finite element method Orthotropic materials Recovery based error indicator

ABSTRACT In this work.we extend the recently proposed adaptive phase field method to model fracture in orthotropic functionally graded materials (FGMs).A recovery type error indicator combined with quadtree decomposition is employed for adaptive mesh refinement.The proposed approach is capable of capturing the fracture process with a localized mesh refinement that provides notable gains in computational efficiency.The implementation is validated against experimental data and other numerical experiments on orthotropic materials with different material orientations.The results reveal an increase in the stiffness and the maximum force with increasing material orientation angle.The study is then extended to the analysis of orthotropic FGMs.It is observed that.if the gradation in fracture properties is neglected.the material gradient plays a secondary role.with the fracture behaviour being dominated by the orthotropy of the material.However.when the toughness increases along the crack propagation path.a substantial gain in fracture resistance is observed.

1.Introduction

Functionally graded materials (FGMs) are a special class of composites with spatially varying microstructure - volume fractions of the constituent elements.These characteristics of FGMs allow the designer to develop ad hoc microstructures for specific,non-uniform service conditions.In addition.the continuous variation of material properties alleviates weak junctions within the system (for example in layered materials).i.e..avoiding the bimaterial interface.which could be a potential site for crack nucleation.The potential advantages in using the FGMs include: (a)enhanced thermal and fracture resistance [1,2].(b) reduced residual stresses [3].and (c) the smoothening of interfaces [4,5].Ceramic-based FGMs enjoy great popularity [6].However.these materials exhibit brittle fracture and complex fracture behaviour[7].particularly when a preferential direction of orthotropy develops.The preferential direction of orthotropy can arise due to the manufacturing process utilized for the synthesis.This is.for example.the case in FGMs manufactured with plasma spray techniques or electron beam physical vapor deposition.In the former,the outcome is a material with a lamellar structure with higher stiffness and weak cleavage planes parallel to the boundary.In FGMs manufactured via electron beam physical vapor deposition one observes a columnar structure.a higher stiffness in the thickness direction and weak fracture planes perpendicular to the boundary [8,9].

Several numerical techniques have been proposed in the literature to analyse the fracture processes in orthotropic FGMs[8,10-14].The vast majority of the works are based on discrete approaches; for example.the conventional finite element with displacement correlation technique(DCT)[15],the extended finite element method(XFEM)[8,10,11,14],and the scaled boundary finite element method (SBFEM) [12,13].However.predicting crack initiation and subsequent crack growth requires an ad hoc criterion,with crack trajectories being sensitive to this choice [16].Variational approaches based on energy minimization constitute a promising tool to overcome this limitation [17,18].Specifically.the phase field method (PFM) has proven to be efficient technique in modelling brittle fracture [19-21].ductile damage [22,23].dynamic fracture [24].fracture properties prediction ofnanocomposites [25].fiber cracking and composites delamination[26-28].plates and shells [29,30]and hydrogen embrittlement[31,32].among other phenomena.Recently.the success of phase field fracture methods has been extended to modelling cracking in isotropic FGMs by Hirshikesh et al.[33].Here.we extend the framework to deal with orthotropic FGMs and include an adaptive mesh refinement strategy to boost computational efficiency.

Although.the PFM has shown advantages over discrete approaches.the finite element discretization requires resolving the length scale parameter as ?o.In brittle materials.?ocan be very small and the discrete crack in linear elastic fracture mechanics is recovered for the limiting case of ?o→0.The need to resolve this region of high gradients creates a computational burden.Local refinement techniques can reduce the computational cost; however.this requires the crack path to be known a priori.which is often not the case.An alternative is to use adaptive refinement algorithms based on error indicators.Several strategies have been proposed [34-39].being most of them based on post-error estimation such as goal-oriented,recovery.and residual.For example,Areias et al.[40,41]presented an adaptive mesh refinement strategy that combines the staggered algorithm with the screened Poisson equation.Goswami et al.[42]proposed an adaptive fourthorder phase field model based isogeometric analysis (IGA).Recently,Samaniego et al.[43]solved the phase-field equations via machine learning approach.In this paper.we aim to extend the recently developed adaptive PFM by Hirshikesh et al.[37]to model fracture in orthotropic FGMs.The adaptive PFM is based on the combination of quadtree decomposition and recovery based error indicators.allowing for an automatic tracking of the crack trajectory and local domain discretization.The hanging nodes that arise due to the quadtree decomposition are treated within the framework of the polygonal finite element method (PFEM) with mean -value coordinate basis function.

The rest of the paper is organized as follows.Section 2 presents the governing equations for the PFM and the corresponding weak form.The adaptive refinement strategy based on the quadtree decomposition and recovery based error indicator is presented in Section 3.The applicability of the adaptive refinement strategy for the fracture in orthotropic FGM is shown in Section 4.Concluding remarks end the manuscript.

2.A phase field fracture formulation for orthotropic FGMs

Consider an orthotropic functionally graded solid with primary orientation directed along the axis e1.making an angle θ with respect to the global frame ex,and secondary orientation e2,which is orthogonal to e1as shown in Fig.1.The boundary (Γ) is considered to admit the decomposition with the outward normal n into three disjoint sets.i.e.,Γ = ΓD∪ΓN∪Γcand ΓD∩ΓN∩Γc= ?.where Γcis the crack surface,Dirichlet boundary and Neumann boundary conditions are specified on ΓDand ΓNrespectively.The closure of the domain is≡Ω∪Γ.

2.1.Governing equations

The spatial variation of the elastic and fracture properties inherent to functionally graded materials (FGMs) can be incorporated following the pioneering work by Hirshikesh et al.[33].Variational phase field fracture methods are particularly suited to capture the complex crack trajectories that are observed in FGMs due to the inherent crack tip mode mixity[33,44,45].As described below.we introduce a history field H to prevent damage irreversibility and we adopt the so-called hybrid model[46]to reduce the computational cost by keeping the linear form of the elasticity equation.In addition,we decompose the strain energy density into tensile and compressive parts ψ = ψ++ ψ-.so as to prevent damage under compressive stresses.Consider a linear elastic solid with spatially varying toughness Gc(x) undergoing small strains.For the displacement u and phase field φ.the strong form of the governing equations in the absence of inertia and body forces is given by Ref.[46,47]:

Fig.1.Schematic representation of a orthotropic FGM domain with a geometric discontinuity in PFM framework.?o is the characteristic length scale.

These balance equations are supplemented with the following boundary conditions:

where ?pand ?eare the scalar and the vector differential operators,given by,

Here A = I+ β[I-n?n].with n= ｛cosθ，sinθ｝Tintroduced to account for the crack path based on the material orientation.In this work,the penalty parameter,β=20 is considered which constraint the propagation of crack in the direction perpendicular to the cleavage plane.The Cauchy stress tensor.σ for the functionally graded orthotropic material is defined as:

where kpis a small positive number introduced for numericalstability and

with

and

The components of the tensor Q(x)are calculated as:

Here.E1(x) and E2(x) are the longitudinal and the transverse Young’s modulus respectively.G12(x) is the shear modulus.ν12is the major Poisson’s ratio,and ν21is the minor Poisson’s ratio.Thus,material properties vary at the element level.in what is usually referred to as a graded finite element approach [48].The small strain tensor(ε) is computed from the displacement field (u) as,

The history variable.H+is defined as,

The introduction of H+in Eqn.(1b)helps to decouple Equations(1a-1b)and a robust staggered scheme can be used for computing(u.φ) [46,49].However.one should note that monolithic quasi-Newton methods have recently shown great promise for phase field fracture problems [50,51].Further.to prevent the crack faces from inter-penetration.Eqn.(1b) is supplemented with the following constraint:

where,

2.2.Weak form

Let W (Ω) include the linear displacement field and the phase field variable.and let (U ，P) and (V ，Q) be the trial and the test function spaces:

Let the domain be partitioned into elements Ωhand on using shape functions N that span at least the linear space,we substitute the trial and the test functions:and ｛v，q｝ =into Eqn.(13).The system of equations can be readily obtained upon applying the standard Bubnov-Galerkin procedure.Find uh∈U and φh∈P such that， for all v∈V 0 and q∈Q0,

which leads to the following system of linear equations:

where

Here.B= ?eN is the strain-displacement matrix and Bφ= ?pN is the scalar gradient of the shape function matrix N.The above system of equations are solved by the staggered approach [46,49].The present framework is implemented in Matlab.The reader is referred to Ref.[52]for a FEniCS-based implementation.Ref.[53]for a COMSOL-based implementation,Ref.[54]for an Abaqus-based implementation.and to Ref.[43]for machine learning solution scheme.

Fig.2.Diffraction method for calculating weight function considering discontinuities in moving least-squares approximations (red line shows the discontinuity).

3.Recovery based error indicator and quadtree decomposition

In this section,we present a brief overview of the recovery based error indicator proposed by Bordas and Duflot[55,56]for the XFEM.This is done to assess the error and identify the elements/regions which have to be refined.Later.the process of quadtree decomposition is discussed.

3.1.Recovery based error indicator

In this method,the enhanced strain field is computed using the standard nodal solution through the eXtended Moving Least Square(XMLS) derivative recovery process.This is then further used as error indicator.Let x be a point in the domain,and nxXMLS points contain x in their domain of influence.Then.using the displacement values at these nxpoints.the enhanced displacement field and the strain field at x can be written as,

where Ψk(x)is the MLS shape function value associated with node k at x,D is the derivative operator and D is the MLS shape function derivative matrix.The matrix form of the MLS shape function is given by:

where p(x) denotes the m reproducing polynomial used for the MLS shape function.For two dimensions.p(x) = [1 x y].and,

Here.A is a m×m matrix and B is a m×nxmatrix.For the matrix A to be invertible.we need nx＞m.i.e..we need more number of points whose domains of influence contains x that the basis functions in p.However.note that this is not a sufficient condition.The weight function wIassociated with a node xIis calculated by the diffraction method with a circular domain of influence.The domain of influence also changes if it intersects with the discontinuity.In this work,a fourth order spline is taken as the weighting function [56]:

The enhanced derivatives of the shape functions are computed by finding the derivatives of the MLS shape functions,see Eqn.(18).The enhanced derivatives of the displacements and the enhanced small strain εscan then be found.The error is computed by comparing the enhanced strain field to the standard compatible

Fig.3.Quadtree decomposition: (a) representative quadtree mesh and (b) tree structure employed to store the mesh details.

strain field.Thus.the total error of domain Ω is,

while the error within the element i with area Ωiis,

The tolerance is chosen based on the maximum error criteria.Thus,the elements with high individual error are discretized in the next level;all elements whose individual error is higher than given tolerance value will be sent for discretization in the next level.

3.2.Quadtree decomposition

The quadtree decomposition is used for local refinement once the error is quantified.The quadtree decomposition entails several features; namely.(a) is easy to implement (b) it requires less degrees of freedom (Dofs) and (c) retains hierarchical meshstructures.The hierarchical mesh structure facilitates efficient computations.particularly efficient storage and data retrieval.In this decomposition.the so-called stopping criterion is used to decide which element requires to be further refined.This criterion can be a geometry based factor or any error indicator.The criteria for an element to be refined based on the error indicator could be based on either equal distribution criterion or Min-number criteria[57].In this work.equal distribution criterion is used to minimizethe global error and balance the local error throughout the domain.If the given element does not satisfy the stopping criterion within the user specified tolerance limit.it will be divided into four child elements as shown in Fig.3.This process can be repeated several times until the criteria is met.The tolerance in all the examples is chosen to be 1× 10-5.

Fig.4.Schematic representation of an element with hanging node and the construction of mean value shape function.

Fig.5.Edge crack orthotropic specimen: (a)geometry,material properties and boundary conditions (b)crack propagation direction for different material orientation.[where L=1 mm and a= 0.5 mm].

Fig.6.Domain discretization of edge crack specimen with θ=30° at (a) 0.035 (b) 0.038 and (c) 0.039 mm.

Table 1 Computational time (in seconds) comparison for adaptive PFM and PFM with uniform refinement.

Fig.7.Final domain discretization for (a) θ=0° (b) θ=30° (c) θ=45° and (d) θ=60° and the corresponding crack trajectory in (e,f,g,h).respectively.

Fig.8.Quadtree evolution as a function of load step.

Table 2 Crack propagation angle compared with the experiments and the XFEM.

The aforementioned decomposition leads to elements with hanging nodes; see.Fig.3.The conventional finite element approach cannot handle such elements without additional work.This is because of lack of compatibility between the elements.Torestrict the number of hanging nodes per edge.a general practice 2:1 rule is applied,in which the mesh is constructed in such a way that two neighboring elements do not differ by more than one level.A number of techniques to handle these hanging nodes have been proposed,such as triangulation[58],transformation of the hanging degrees of freedom to corner degrees of freedom using constraint equations [59].use of special conforming shape functions [60],considering the element with hanging nodes as a polygon [61].or the use of other advanced methods like Scaled Boundary Finite Element Method (SBFEM) or Smoothed Finite Element Method(SFEM) [62].

In this work,the elements with hanging nodes are considered as n-sided polygons (see Fig.4a).The mean - value shape functions proposed by Floater [63]are used to approximate the unknown fields.The reasons behind this choice are the facts that some angles of elements with hanging nodes are 180°and the mean value shape-functions work efficiently for non-convex polygons.

The mean - value coordinates for a point P(x) in an arbitrary polygon are given by:

where n is the number of nodes in an element.xiare the coordinates of point Piand αi’s are the internal angles.Fig.4b shows the mean value shape function for the polygon with hanging node.The numerical integration for the polygonal elements is performed by subdividing the polygon into triangles and employing standard quadrature rule.

4.Results

In this section.the performance and the robustness of the adaptive PFM for fracture of orthotropic FGMs is investigated.We first validate the adaptive PFM results against experimental and numerical results for failure of orthotropic materials.Then,cracking of orthotropic FGMs is investigated for different material grading possibilities.The numerical stability parameter kpis assumed to be 1×10-6in all the numerical examples,unless specified otherwise.The proposed adaptive PFM is implemented in MATLAB R2014b and the simulations were performed on Intel quad Core i5-4590CPU@3.30 GHz with 8 GiB RAM.

4.1.Validation: fracture of orthotropic materials

The framework developed is validated first with the experimental and numerical work(XFEM)of Cahill et al.,[64].In order to study the fracture processes in an orthotropic material.an edge crack specimen subjected to tensile loading is considered.see Fig.5a.The material properties are chosen as:E1=114.8 GPa,E2=11.7 GPa.G12=9.66 GPa.ν12=0.21 and the critical toughness G c =2.7 MPa mm.

The simulation starts with a coarse mesh and an assumed characteristic length scale.?o.For each load step.the domain is discretized as explained in Section 3.which allows to track the crack trajectory continuously.Fig.6 shows the domain discretization for the evolving crack trajectory;the combination of quadtree decomposition and post-errori error estimator strategies leads to a fine discretization in the vicinity of the propagating crack tip.This feature substantially reduces the size of the global stiffness matrix,resulting in reduced CPU memory requirement and efficient computations.

We proceed to quantify the computational gains by comparing to the case of uniform refinement.Computation times are compared to those obtained from a solution with uniform mesh,where the characteristic element size equals the element size in the crack region in the adaptive re-meshing case.Results are shown in Table 1 for the initial load step Δu.after convergence has been achieved.The number of degrees of freedom(DOFs)is also shown.Very substantial computationally gains are reported relative to the uniform mesh scenario,due to the significantly smaller number of DOFs required when using an adaptive mesh refinement strategy.The error indicator is the most time consuming operation in the adaptive PFM.

Fig.9.Load displacement response for the orthotropic material with different material orientation.θ.

Fig.10.Orthotropic FGM specimen: domain and boundary conditions.(α，β，γ) are indices that control the material variation.

Fig.5b shows the crack propagation trajectory for four selected values of material orientation i.e.,θ = 0°，30°，45°，60°.The corresponding domain discretization is shown in Fig.7.In agreement with expectations.the crack propagation path strictly follows the material orientations.see Fig.5b.

Fig.8 shows how the number of quadtree elements increases as the crack propagates.It is shown that the material orientation angles that translate into a larger crack deflection require a larger number of quadtree elements.The comparison with the results obtained in experiments and XFEM-based calculations [64]are shown in Table 2.A very good agreement is attained for all values of material orientation.

Finally.Fig.9 shows the load-displacement response of the orthotropic specimen with different material orientations.The stiffness of the response increases with the material orientation angle θ.

4.2.Fracture of orthotropic functionally graded materials

Next,we examine the fracture processes in an orthotropic FGM specimen as shown in Fig.10.In terms of material gradation.we consider the following representative case studies:

· plate with a crack parallel to the material gradation i.e..xdirection grading,

· plate with a crack perpendicular to the material gradation i.e.,y- direction grading.

For simplicity.we assume that the material property variation follows an exponential gradation,as characterized by the indices α，β，γ - see Fig.10.As shown in Fig.11.two different scenarios have been considered: (i) proportional and (ii) non-proportional gradation strategies.in terms of the material indices.In the former.the indices,α，β，γ are set to 0.2.whilst for non-proportional variation,we choose.(α，β，γ) = (0.5，0.4，0.3).The material constants arechosen as: E01= 114.8 GPa.E02=11.7 GPa.G012= 9.66 GPa.and the critical toughness G c = 2.7 MPa mm.With respect to the critical energy release rate Gc.two cases are considered; one.where no material gradation is assumed and a graded one that follows the material gradation depicted in Fig.10.For all results,the orthotropic material orientation is take to be θ = 0°.

Fig.11.FGM orthotropic material gradation in (a) x- direction with proportional material gradation (b) x- direction with non-proportional material gradation (c) x- direction with proportional material gradation.and (d) y- direction with non-proportional material gradation.

Fig.12.Load-displacement response for FGM orthotropic specimen with (a) with constant toughness G c.and (b) with varying toughness G c(x).

In all cases.given that the same value of θ is considered.the predicted crack trajectories follow an almost identical path.However.differences can be seen in the load-displacement curves.as shown in Fig.12.Consider first the effect of a proportional or nonproportional material gradation.the same qualitative trends are seen in both Fig.12a and b.For the case of material gradation in xdirection.the non-proportional material gradation shows stiffer response than the proportional material gradation.This trend is reversed for the case of material gradation in y- direction,where differences are minimal.Differences are due to the crack tip nonhomogeneity.which affects the mode mixity.Consider now the influence of spatially varying the material critical energy release rate G c toughness.i.e.Fig.12a versus Fig.12b.It can be seen that differences are substantial in the case of material grading along the x-direction.In agreement with expectations,the propagating crack encounters an increasing resistance to fracture as the magnitude of G c at the crack tip raises with crack advance.

5.Conclusions

We have presented a novel framework for modelling fracture problems in orthotropic functionally graded materials(FGMs).The framework builds upon the phase field fracture method for FGMs and an adaptive mesh technique based on a recovery based error indicator and quadtree decomposition.Results show the capability of the model in capturing complex crack trajectories,not known a priori.while minimising the computational cost.The numerical framework is validated by comparing to experimental and numerical results on non-graded orthotropic materials.A good agreement is observed.Then.calculations are shown for orthotropic FGMs and the role of the material gradation indeces explored.Topics of interest for future work involve extending the present framework to dynamic crack growth.three dimensions problems and enabling mesh coarsening behind the crack.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

E.Martínez-Pa~neda acknowledges financial support from the Royal Commission for the 1851 Exhibition through their Research Fellowship programme (RF496/2018).