Full thermo-mechanical coupling using eXtended finite element method in quasi-transient crack propagation

This work aims to present a complete full coupling eXtended finite element formulation of the thermo-mechanical problem of cracked bodies. The basic concept of the extended finite element method is discussed in the context of mechanical and thermal discontinuities. Benchmarks are presented to validate at the same time the implementation of stress intensity factors and numerical mechanical and thermal responses. A quasi-transient crack propagation model, subjected to transient thermal load combined with a quasi-static crack growth was presented and implemented into a home-made object-oriented code. The developed eXtended finite element tool for modeling two-dimensional thermo-mechanical problem involving multiple cracks and defects are confirmed through selected examples by estimating the stress intensity factors with remarkable accuracy and robustness.


Background
The interest in fracture mechanics and its applications has gained considerable importance in recent years in various industries: aerospace engineering, automobile industry, civil engineering, etc. This attention is due to the high cost caused by the presence of cracks and defects, which require more energy, time, substantial efforts and dedicated strategies regarding intervention, maintenance, repair, etc. Practically, taking into account the real environmental conditions in service has become essential, when the material is subject to a significant gradient of temperature. For instance, temperature change in real structures, where the deformation are constrained, can engender a mechanical load and a high-stress concentration around crack tips. Subsequently, crack can propagate with a, a priori, not known orientation, direction, intensity etc. Since, cracks cannot be eliminated under any circumstances; this prompts engineers to guide our efforts towards winning strategies in prevention, design and especially analysis that can be provided by the tool of numerical modeling. The numerical modeling of cracked domains using finite element method (FEM) has clearly stood aside for the eXtended finite element method (XFEM) in the last two decades. XFEM has been able to provide essential answers for several situations, where the FEM method becomes numerically very expensive to have an optimal convergence, such as singularities, strong discontinuities, high gradient, moving surfaces, etc. This technique allows, by prior knowledge of the physical behavior of the problem, to enrich the space of the solutions by non-polynomial asymptotic functions when it is a singularity and a jump-function when it comes to a discontinuity or a combination of both of them. The resulting approximation space has to reproduce the Partition of Unity (PU), Babuška and Melenk [1]. The first work that introduced enriched FEM was Belytchko and Black's paper [2] which presented an implicit description of the crack with minimal remeshing. Moës et al. [3] improved this technique by incorporating a more suitable way to consider the discontinuities throughout the crack faces away from the crack tip by the generalized Heaviside function and branching functions for the near crack tip. Daux et al. [4] later extend the approach for multiple cracks and holes for the mechanical problem. Sukumar et al. [5] used the XFEM to model fracture in three-dimensional by using the PU concept, where the two-dimensional asymptotic crack tip displacement fields were added to the FE approximation to account for the crack. The XFEM for non-planar cracks in three dimensions illustrating the crack geometry using two signed distance functions was presented by Moës et al. [6]. Sukumar and Prévost [7] extended XFEM for two-dimensional crack modeling in isotropic and bimaterial media and later to demonstrate the numerical modeling of stress intensity factors in crack growth problems in Sukumar, and Prévost [8]. Lee et al. [9] exposed a combination of the XFEM and the mesh superposition method for modeling of stationary and growing cracks, where a step function implicitly described the discontinuity on the PU, and the crack tip was modeled by superimposed quarter point elements on an overlaid mesh. Budyn et al. [10] displayed a model for multiple crack growth considering the junction of cracks in brittle materials using XFEM, which does not require remeshing as the cracks grow. Other XFEM aspects have been addressed: In contact, Khoei et al. [11] used XFEM to model the frictional contact problem using the penalty method. Nistor et al. [12] developed an approach to couple the XFEM with the Lagrangian large sliding frictionless contact algorithm. An algorithm based on node-to-segment XFEM contact was presented by Khoei et al. [13] based on the XFEM to model the large deformation-large sliding contact problem using the penalty approach. In stabilization aspect, an XFEM pre-conditioner which stabilizes the enrichments by applying Cholesky decompositions to certain submatrices of the stiffness matrix was proposed by Béchet et al. [14]. Menk et al. [15] expose another pre-conditioning method suited for parallel computation. Also, another approach initially developed by Hansbo et al. [16] to simulate strong and weak discontinuities in solid mechanics. A similar method was used by Song et al. [17], named phantom nodes, for shear modeling dynamic crack and shear band propagation. Rabczuk et al. [18] developed a new crack tip element for the phantom node method suited for one-point quadrature scheme and can be used with other general quadrature schemes. XFEM numerical integration aspect is performed by Dolbow et al. [19] by using a sub-triangulation for computing the element area below and above the crack and to set criteria for node enrichment with discontinuity function. Laborde et al. [20] used a singular mapping for each sub-triangle and a bidirectional Gauss quadrature in each direction. In Ventura [21], the constructing sub-cells in the numerical integration of discontinuity functions is removed by defining an equivalent polynomial function. Schwarz-Christoffel conformal mapping was used to map an arbitrary polygon onto a unit disk by Natarajan et al. [22]. A fairly comprehensive review of the different aspects of XFEM was presented by Khoei [23]. All these advances in XFEM mentioned before are in the field of solid mechanics. In this paper, the approach taken is based on a semi-implicit thermo-mechanical-crackgrowth algorithm in which the combined full coupling thermal and mechanical responses have to be estimated beforehand. Then, the developed numerical fracture mechanics module takes those responses as inputs to evaluate the stress intensity factors, J-integral, the update of the crack in growth, etc. This actualization is done by an implicit description of the crack, using the Level-Set Method (LSM) presented firstly by Osher and Sethian [24]. The LSM provides a fundamental complementary to know when, where and how to enrich the crack by determining its relative position. Stolarska et al. [25] introduced an algorithm that combines the XFEM and LSM to model mechanical crack growth, where the LSM was used to model the crack surface and crack tip locations. Moreover, stress intensity factors (SIFs) computation, as the prime parameter of prediction, makes it possible to obtain an essential knowledge of the behavior of the crack. This evaluation enables to predict whether the structure becomes unsafe in service conditions, especially when it is in a thermo-mechanical context, where the spatial distribution of the mechanical stresses induced by the thermal field is unpredictable. Interest in thermo-mechanical applications appeared later with Michlik and Berndt [26] presented an approach of thermo-mechanical XFEM analysis to account for the existence of cracks in thermal barrier coating for predicting an effective thermal conductivity and Young's moduli of multi-layered. Duflot [27] used the XFEM for the analysis of steadystate thermally stressed, cracked solids in thermo-elastic problems, where he enriched both thermal and mechanical fields to represent the discontinuous temperature and displacement. Fagerstöm and Larsson [28] presented a thermo-mechanical fracture formulation based on discontinuous representation for temperature and displacements fields applicable to the fracture process zone into a cohesive zone. Zamani et al. [29] proposed a higher order XFEM to predict the SIFs for thermo-elastic with stationary cracks, The computation of SIF is extracted directly from the XFEM degrees of freedom. Zamani et al. [30], in a later work, implemented the XFEM to model the effect of the mechanical and thermal shocks on a body with a stationary crack. Lee et al. [31] presented an XFEM method for the analysis of heat conduction at submicron scales of geometrically complex nanostructured heterogeneous materials. Fan et al. [32] used XFEM to investigate the effect of thermally grown oxide on multiple surfaces cracking behavior in an air plasma sprayed thermal barrier coating system. Hosseini et al. [33] introduced a computational method based on the XFEM for fracture analysis of isotropic and orthotropic functionally graded materials (FGM) under mechanical and steady-state thermal loadings. Yu et al. [34] exploited XFEM for modeling the temperature field in heterogeneous materials, where the standard temperature field was enriched by using the level-set-based enrichment functions which model the interfaces. Macri et al. [35] presented a multiscale technique for modeling heterogeneous materials based on an enriched partition of unity that incorporates the thermal effects occurring on the microstructure into the global model for simulation. In Sapora et al. [36], an analogy between fracture and contact mechanics is proposed to investigate debonding phenomena at imperfect interfaces due to thermomechanical loading and thermal fields in bodies with cohesive cracks. From fracture mechanics point of view, Goli et al. [37] implemented the path-independent interaction integral in the context of the partition of unity for mixed mode adiabatic cracks under thermo-mechanical loadings particularly in orthotropic non-homogenous materials for a steady-state thermal problem. Bayesteh et al. [38] study a thermo-mechanical fracture of inhomogeneous cracked solids by the extended isogeometric analysis method, crack faces, and tip XFEM enrichment are incorporated into the non-uniform rational basis splines functions of isogeometric analysis (IGA) for static crack and steady-state thermal problem. Jia and Nie [39] adopted XFEM to analyze the interaction between a single or multiple macroscopic or microscopic inclusion and cracks for static crack and under the steady-state thermal problem. The work of Jaśkowiec [40] is concerned with modeling the heat flow through cracks in three-dimensional thermo-mechanical problems, the model for crack heat flow is combined with cohesive crack model. He et al. [41] established an XFEM thermo-elastic fracture problem for aluminium alloy metal inert gas welding, which includes a variable heat source with the initial and boundary conditions for a cracked plate structure. Li and Fish [42] developed a thermo-mechanical extended layerwise method for the laminated composite plates with delaminations and transverse cracks; transverse cracks are modeled using classical XFEM under pure mode-I. Recently, Zarmehri et al. [43] implemented XFEM to extract stress intensity factors for a stationary crack in an isotropic 2D finite domain under thermal shock, the coupled generalized thermo-elasticity theory employed is based on Green-Lindsay model. Although the plethora of works has treated numerical thermo-mechanical analysis using classical XFEM recently, few works have employed the enhanced version of XFEM named XFEM-f.a. in order to ensure an optimal convergence through a geometrical enrichment regardless of the mesh size. This work aims at developing the complete full thermomechanical coupling using XFEM in adiabatic cracked media adopting a geometrical enrichment. The implementation was firstly validated for a single crack from the existing examples in the literature. Then, validation of the case of the combination between a hole and cracks, and the influence of crack size and a single hole size on the stress intensity factors, i.e., on the behavior of the rupture in a given structure, is performed. This case was investigated with the work of Prasad et al. based on the dual boundary element method for thermo-elastic crack problems [44]. Notably, the case of transient thermal loading and its impact on the SIFs profiles was treated, then a situation of the growth in mode-I was analyzed. Moreover, this work study the case of the thermo-mechanical propagation of multiple cracks in the presence of multiple holes in mixed mode. This paper consists of six sections. The second one sketch the mathematical, physical and variational framework of the two-dimensional plane strain thermo-mechanical problem studied in a cracked medium. The third section intended for approximation spaces and the XFEM discretized forms of displacement and temperature fields as well as the full coupling XFEM matrices for each sub-problem part and the integration technique employed. Section four deals with the crack growth criterion assumed in this study, the form of the thermo-mechanical J-integral and the extraction of SIFs. Section five describes the specific numerical approach of a modified XFEM version involved. Several validation models from the literature are then considered for validation purpose; another benchmark with cracks and a manufacturing flaw idealized by a hole; an example of crack growth in mode-I, then in transient thermal loading; and lastly a mixed-mode crack growth model designed carefully for a specimen with multiple holes and cracks in the thermo-mechanical case. Finally, we conclude by a summary and some proposed extensions of this work.

Governing equations
Consider a linear-elastic, isotropic and homogeneous body occupying a geometrical cracked domain bounded by = ∂ in Fig. 1. The boundary is composed of parts u , T , t , q and c . The equations of thermo-mechanical problem, assuming small displacements and small strains on \ c are The objective is to find u(x, t) kinematically admissible, where, T f is the end time. The fields are displacement vector u, temperature T , stress tensor σ , strain tensor ε, 'thermal strain' tensor ε T defined with respect to a reference temperature T 0 , heat flux vector q; the properties of materials are scalar thermal conductivity k, thermal expansion coefficient α, density ρ, specific heat capacity c, and the isotropic fourth-order Hooke tensor C, Q(x) and b(x) are respectively the imposed heat source and the body force on \ c . I is the second-order identity tensor and ∇ s is the symmetric gradient operator on a vector field. Prescribed displacements u and temperatures T are imposed respictively on u and T , while tractions t 0 and heat flux q are imposed on t and q as u = u on u , σ · n = t on t , σ · n = 0 on c , T = T on T , q · n = q on q , q · n = 0 on c , The crack surface c is assumed to be traction-free. The problem is well-defined with T ∪ q ⊂ , T ∩ q = ∅ and u ∪ t ⊂ , u ∩ t = ∅. Thermal problem is merely time-dependent, while the mechanical one is quasi-static by neglecting inertial effects, Khoei et al. [45]. The time appeared in σ (., t), ε(., t) is a pseudo-time induced by the real-time in the time-space ]0, T f ]. In the staggered thermomechanical problem, which is not the case for this present work, the transient thermal problem is solved first to compute the temperature field at the real-time t, then the quasistatic mechanical problem defined by the equation Eq. (3) is solved by taking T (t) as input. Consequently, the resolution of the mechanical problem, in pseudo-time 'marching', becomes conditioned continuously in real-time by the resolution of the thermal problem. Hence, by splitting the mechanical stress σ (., t) to a global stress σ g (., t) and 'thermal stress' σ th (., t), the continuous space of mechanical pseudo-times can be defined for a traction-free crack by For a given real-time t, space PT t will be reduced to a singleton. This definition explains that for each real-time a unique pseudo-time is defined implicitly and naturally with the temporal evolution of the system. This allows to simply considert ≡ t. Also, from the numerical discretized point of view, it is possible to study a direct steady problem and to avoid the pseudo-time stepping with an associated small-time step in the non-steady state. Since, in this study, the problem is strongly coupled, Eqs. (1)-(6), the overall dynamic of the system is driven by the transient problem. Thus, the causal real-time retrace the 'dynamic' of the mechanical problem which becomes pseudo-time-dependent. Subsequently, the entire problem is treated as a monolithic object, and all sub-problem parts progressed simultaneously. Note that in case of several cracks, the mother crack c can be decomposed into many n adiabatic cracks, c = n 1 c i such every ith crack remains adiabatic q · n = 0 and traction-free σ · n = 0 on each c i , for any i ∈ 1, n . Henceforth, we will present the XFEM developments for the mother crack, which remains valid for all sub-cracks.

Variational form
The space of admissible displacement and temperature is U × T , X = (u, T ) ∈ U × T , where variational spaces U and T are defined on Sobolev space H 1 ( ) by U = {u ∈ H 1 ( ) 2 : u = u on u and u is discontinuous on c }, T = {T ∈ H 1 ( ) : T = T on T and T is discontinuous on c }, and the spaces of homogeneous essential conditions are given by The weak form of the thermo-mechanical problem, with the test functions v and w, can be expressed as follow: Find u ∈ U and T ∈ T such

Full coupled eXtended Finite Element form
Considering a finite element mesh M without taking account of the crack which is treated separately by an implicit description using the level-set method. The XFEM shifted discrete form of each component, u ∈ {u, v}, of the displacement field on M takes the following form: where A represents the whole set of nodes forming the mesh including all enriched nodes, A M; A cr describes all the nodes building the elements crossed by the crack without tips, A cr ⊂ A; and A tip denotes the nodes constructing the tip elements, A tip ⊂ A. N A , N A cr and N A tip are the countable sets of the nodes, respectively, of A, A cr and A tip . The singular asymptotic basis functions are given in polar coordinates (r, θ ) by Similarly, the discrete form of the temperature field with a single asymptotic function, and with the same definitions mentioned above can be written as namely that the asymptotic expression of the temperature field of an adiabatic crack can be expressed on the point (r, θ ) in the polar reference centered on the corresponding crack-tip, Duflot [27], as Henceforth, we adopt the notations of u := u h and T := T h for numerical development, without making a difference between the continuous and the discrete form of both displacement and temperature.
Time discretization can be obtained by assuming that two displacement-temperature {X i } at time t i and {X i+1 } at time t i+1 , t i+1 = t i + t, are related by the generalized trapezoid rule, including a parameter to set β, To improve stability of the previous scheme, we choose the particular case of Crank-Nicolson, with β = 1 2 which is unconditionally stable and with no numeric dissipation into the numerical approximation. Let X be the global variable of the full XFEM thermo-mechanical coupling, such where {U std } and {U enr } are respectively the standard part and the enriched part of the displacement; the same goes for the temperature. The full coupling form of the XFEM stiffness matrix, damping matrix and force vector are given by where [K UU ] is the purely mechanical part; [K TT ] is the purely thermal part; [K UT ] is the coupling part describes the influence of the thermal problem on the mechanical one; and the zero term matrix explains that there is no influence of the mechanical problem on the thermal one.
• Mechanical part: The unknowns for mechanical problem are {U i , b u i , c u i } T ; the mechanical stiffness matrix can be written as where the expression of each sub-matrix component, with respect to the Lamé's constants μ and λ denoting in the indicial notation of the stress σ ij = 2με ij +[λε kk −α(3λ+2μ)(T − T 0 )]δ ij , can be given explicitly by and, and, and, and, and, The mechanical part of the force vector, • Thermal part: The unknowns for thermal problem are the thermal stiffness matrix can be written as where the expression of each sub-matrix component can be given explicitly by and The thermal part of the global force vector, • Coupled part: The thermo-mechanical coupled part can be expressed as And the damping non-zero matrix [C TT ] is given by where each component can be expresses by,

Numerical integration
In XFEM, the standard Gauss approximation approach cannot be used for the elements crossed by the crack. It is then necessary to modify it appropriately to evaluate the contribution of the weak form W u t and W T for the two compartments generated by the crack at the sub-domain element e level. Indeed, the XFEM numerical integration requires a particular treatment due to the complexity encountered when integrating elements traversed by the discontinuity ( c,e = e ∩ c ) or crack-tip-element where the approximation functions are non-polynomials. The enrichment remains local in the vicinity of the crack region, termed enriched-zone; therefore, the area affected by this special treatment is located at the level of the enriched-zone. Beyond this zone, the elements are considered as standards with 4 or 8 integration points per non-enriched element. In the elements crossed completely by the crack c e , split or vertex element, the resulting configuration yields to a convex domain C and a complement non-convex domain C c c e . This situation needs a suitable geometrical approach, to deal with all possible cases, by a sub-polygons subdivision process to form a convex disjoint partition. A set of m e sub-convex-elements K of the same dimension taking account the relative position of the crack at the subdomain element, such e = m e 1 K , Dolbow et al. [19] and Laborde et al. [20], is taken. The same procedure can be done for tip-elements, or partially crossed by the crack, with much more attention due to the non-polynomial aspect of the functions of approximation. Commonly for both cases, the whole subdivision at the local element level can be seen as a spider-web delimited by the element borders, Fig. 2, centered on the crack-tip for the tip-element case and on the iso-barycenter of C or C c for split or vertex element cases. In split element each sub-element take 7 integration points per K , in total we obtain 7 * m e points per split/vertex element. Also, we assume more integration points on the tip-elements to capture well numerically the singularity by 19 integration points per K ; in total, we get 18 * m e per tip element. To note that to refine the XFEM approximation on the elements of transition between fully enriched elements and standard ones, we keep the same treatment as tip element with a spider-web centered on the iso-barycenter of the element.

Enriched zone update in crack growth
In the initial state of propagation, the position of the crack is predefined 0 c ; the mesh M is properly generated, once and for all, without any change on it during the process of growth. Then, the relative position of the crack is implicitly identified by level-set functions, leaving aside its known Cartesian global position. Consequently, crack is recognized independently of the mesh definition, relatively, with respect to its nodal environment thanks to the signed-distance function. At this early stage, to write the discrete XFEM form of the displacement, Eq. (14), and the temperature, Eq. (16), fields, we ought to select the two kinds of enriched nodes families, Heaviside and tip enriched nodes. All the mesh nodes including those wholly enriched are roughly approximated by the standard shape functions. The nodes enriched by Heaviside function are described by the nodes forming the elements thoroughly crossed by the crack. The tip candidate nodes are such as those nodes composing the tip element for a topological enrichment; in geometrical enrichment on a given disk, D of radius R and centered on the tip, the discrete set of tip nodes is formed by the intersection of the mesh M and the disk D, M ∩ D. In geometrical enrichment instance, D crosses undoubtedly, for a large radius, the already selected Heaviside nodes. As a result, these nodes should be enriched simultaneously by the combined form of branching function and Heaviside, F + H. This selection configuration is performed for the initial increment and will identify the enriched-zone EZ 0 related to the position of the crack 0 c which is supposed to be unique. At the next increment, the same procedure is followed until a (k − 1)th incremental crack is reached, k−1 c . At the kth increment, we suppose that the chosen crack growth criterion is satisfied, which allows extending the crack in the appropriate direction. This progress creates a new geometrical configuration of the crack and a new enriched-zone EZ k to be identified. We proceed in the same way, by recurrence, as the (k − 1)th increment. This time we find ourselves with two different configurations where some nodes in the previous increment that they were, for example, of Heaviside type become tip nodes, or standard nodes convert to Heaviside/tip nodes, or vice versa. The transport of nodal fields corresponding to temperature and displacements computed at the (k − 1)th increment to the kth one will be performed by an L 2 -projection on the discrete space generated by the kth recent configuration by means of least squares method. This strategy is possible since the different quantities are square-integrable, which allows ensuring the stability and efficiency of the used scheme. To note that, the advancement of the crack generates a new geometrical, topological and numerical reality of the thermo-mechanical problem resolution which requires a specific treatment at each increment. This results in an extensible and flexible set of degrees of freedom with the crack evolution; therefore, the linear system of discrete equations is also extensible and changes in dimension depending on the crack growth state, it can enlarge or diminish. On the other hand, the selection rule of the EZ k nodes is independent of the previous configuration and related only to the relative position of the crack in its nodal environment at the actual increment. We are thus left with two different configurations of two successive increments. This process is repeated iteratively until the estimated or evaluated end of the propagation process.

Propagation criterion and crack update
It has been shown that the use of level-set function plays an essential role in the implicit description of the crack and evaluation of enriched fields, mechanical Moës et al. [6] and thermal in this work. The crack is representing the zero-level set of a given function. The crack tip positions can be found by considering the intersection between zero-level contour and a second orthogonal level-set function Stolarska et al. [25] using the signeddistance function. The signed-distance in the level-set method is represented by a finite element approximation with the same mesh used for the mechanical and thermal problems. Adopting this representation makes the task easier when it is necessary to evaluate the level-set at element level by interpolation and when we need to compute its derivative which is well-defined by the derivative of shape functions. To monitoring crack growth, we use the maximum hoop (circumferential) tensile stress theory introduced firstly by Erdogan and Sih [46]. In mixed-mode, the information is extracted in the vicinity of the crack tip by evaluating the stress state, written in polar coordinates. We assume that the crack extension starts at its tip in a radial direction, it is produced in the plane perpendicular to the direction of uttermost tension, i.e., at a critical angle θ c , and it begins when σ θθ reaches a critical given value. When K II = 0 then θ c = 0 also, in this case, we have a pure mode-I. By considering K II < 0 the critical crack growth θ c > 0, and if K II > 0 the angle θ c < 0. A handy expression of θ c was given by Sukumar [7], The extension of the crack path is determined by a constant increment of growth as an attractive approach. The selection of a is almost always made a priori as an input parameter of the numerical crack propagation model. Several settings affect the quality of crack propagation path; those factors are widely studied by [8] using many examples illustrating the impact of these choices on the path. Therefore, it is more judicious numerically to choose a a that takes into account those number of parameters Belytschko [2] to ensure convergence toward the appropriate path. Principally, three parameters influence the quality of the crack path: Firstly the crack growth magnitude (length of the crack incremental segment) which have to be considered within a range of l e ≤ a ≤ 3 2 l e , such, the element size l e = √ A e and A e is the average area of the elements. Secondly, the mesh size is important to have the best approximation of the field near crack with a finer one. Finally, the choice of J-integral domain is decisive to evaluate adequately the J-integral which allow extracting stress intensity factors in mode-I, II and in mixed mode and determining after that the value of crack growth orientation θ c .
On the other hand, different crack extension criteria exist in the literature and adequately ensures the crack progression, governed by fatigue law varieties. They are adapted to the crack progress when it is subjected to cyclic loading. The crack rate increment with respect to the loading cycle, i.e., speed growth, appeared in these laws and assumed to be, in general, a function which depends on the stress intensity factor range between two cycles and the stress ratio, Beden et al. [47]. The popular one is the classical law of Paris which is a version of the general law of fatigue, where the speed growth depends on the stress intensity factor range, and two constants, called constants of Paris law, that have to be identified for each specific material, Cherepanov et al. [48]. Its limitation lies in the fact that it requires a minimum stress intensity factor to ensure the propagation and does not take into account the stress ratio. Another version appeared later by Xiaoping et al. [49] that overcomes these limitations of the classical Paris law but requires three additional parameters more than the classical Paris law. All these models can 'better' capture the crack progress and monitor the history of the adapted crack increment for each promotion. They are more suited to fatigue propagation fashion and also require additional parameters related to the material that can be determined by fatigue tests. This last point may be a drawback for the attractiveness of these methods for the present work. However, the convergence of fixed crack increment method may be 'lower' in some cases, but with a suitable choice of the crack increment, which depends on the mesh and other parameters as cited previously, one can reach good results. Besides, fixed crack increment technique is more attractive; it needs less material parameters compared with the earlier mentioned laws. Its ability to obtain crack paths that coincide very well with reference solutions is investigated by Baydoun et al. [50].

Stress intensity factors evaluation
The J-integral, with free body force b, was introduced by Rice [51] as a way to compute the energy release rate G. Rice defined a line path independent integral, which keeps the same value for any path surrounding crack tip as where δ .. is the Kronecker operator, n i is the component in i-direction of the normal outward vector to the contour ε , t j = σ ij n i and u j are components of the interior traction and displacements, W is the strain energy density per unit volume defined in the thermo-mechanical state by where T = T − T 0 , ε m ij represents the mechanical part of strain and ε t ij denotes the total strain. The form of the integral (62) is not adapted for a finite element computation, in particular in XFEM, while preserving the same shape functions. An enclosed contour * is considered as a sum of piecewise lines as defined in Fig. 3, * = γ + ∪ γ 0 ∪ γ − ∪ γ 1 . Hence, the J-integral can be converted into a domain integral by introducing a weight function q in the expression of (62), that is unity on γ 0 , zero on γ 1 and varying monotonically in-between. In this work, we used a plateau truncated cone. By applying the divergence theorem, the equivalent domain integral (EDI) form of the J-integral is obtained as The J-integral of the superimposed of two equilibrium states: State 1 with u, σ and ε corresponds to the real state and state 2 with u aux , σ aux and ε aux corresponds to an auxiliary situation, is given by which can be explicitly written as By developing J s , the interaction integral I is obtained by By assuming crack faces to be traction free, using equilibrium (i.e., σ ij,j =0), straindisplacement equations, and after some handling, we obtain For general mixed mode problems and isotropic materials, the direct relationship between J-integral and the stress intensity factors, having dimensions of [stress. lenght ], in mode I and II is given by The extraction of individual mode-I and mode-II stress intensity factors can be done by the choice of K aux I = 1 and K aux II = 0 to find K I and K aux I = 0 and K aux II = 1 to find K II as

Numerical examples of thermo-mechanical analysis
A set of thermo-mechanical examples are herein discussed by considering a strong material discontinuity; for a static adiabatic crack and in propagation state of an isotropic material. Validation of the results is fulfilled by a comparison with the computation of the stress intensity factors which allows validating both mechanical and thermal responses as well as the quantification of the linear elastic fracture mechanics (LEFM) parameters. The computation domains chosen for the benchmarks are extracted from the literature and meshes are generated using Gmsh [52]. A hybrid object-oriented code has been developed in a monolithic multi-physical philosophy treating each step starting from the mesh generated from Gmsh, the definition of the enrichment-zone, the XFEM matrix computation blocks associated to each physical segment and to each coupled part, the computation of fracture mechanics quantities and post-processing context. The rate of convergence of conventional XFEM, using a 'topological' enrichment, is not improved when the characteristic mesh length h goes to zero because of the presence of a singularity. Laborde et al. [20] proposed a modified version of XFEM by enriching a whole fixed area (f.a) around the crack-tip, named XFEM-f.a. In the standard XFEM, only the nodes of the crack tip element are enriched by branching functions, the support of the additional basis functions vanishes when h is going to zero. In two dimensions, the fixed enriched area of a radius E j R according to the jth crack-tip is giving by the disk The major drawback of 'topological' enrichment is that the size of the enriched zone depends linearly on the size of the mesh. However, 'geometrical' enrichment has an asset by enriching all the elements containing in the disk D . for a given radius E . R regardless of the mesh size. Therefore, XFEM-f.a. achieves the expected optimal rate of convergence of O(h). For a given configuration where several singularities (crack-tips) are apparently present, the fixed enriched area defined by gathering the multiple disks assigned to each singularity, ensuring that they remain disjointed by a judicious choice of the radius of each disk. Then, the global f.a. is given by where N tip is the discrete set of crack-tips. The discrete approximations of displacement, Eq. (14), and temperature, Eq. (16), by XFEM keep the same expression with a significant change in the topological enrichment of the crack-tip. Thus, the set N A tip of the nodes enriched by branching functions is transformed to N D which represents all the nodes forming the geometrical enrichment zone established by D. It is noteworthy that the effect of the blending elements decreases systematically with the increase of the enrichment area on the whole D. Also, there is no significant effect observed on the numerical solutions Fries [53]. Numerical results are performed for a full thermo-mechanical coupling problem in a cracked domain using a plane strain analysis, where the mechanical loading is induced by a pure thermal one under the prescribed temperatures and flux on boundaries. This case represents the most relevant situation, which can be easily combined with a pure mechanical load acted by external forces. One can simplify the analysis by considering = T − T 0 , with T 0 = 0 initially for the whole domain, and with no heat source Q and no body force b, which is the case for our analyzes. The crack surface is thermally insulated, so the flux lines have to circumvent the crack. Stress intensity factors computation is commonly normalized with respect to another choice of the triplet (E, k, α) material and with a fixed value of Poisson ration to 0.3. Normalized SIFs are presented for all the examples, including the negative values of K I , for a static crack, which represents an important indicator of the compressive effect at the crack lips. The contact between crack surfaces is not taken into account in the numerical model, leaving XFEM-f.a. to produce information that can predict an inter-penetration of crack faces for certain thermo-mechanical configurations/domains. This plight remains entirely true, valid and adequate from a conceptual point of view. The J-integral radius is considered as a function of enrichment radius and have to be greater than or equal to E j R to obtain good results of SIFs computations and   Table 1.
In this section, we present various examples to validate the thermo-mechanical model by XFEM-f.a. implementation, comparing with several benchmarks taken from the literature. The primary objective of all these cases is to investigate the accuracy and robustness of the numerical results. Then, we present an example of the cracked domain under transientthermal load and in crack growth governed by mode-I. Finally, we design a model of crack propagation in mixed mode for a pure thermal loading, with round holes and multiple cracks.

Rectangular plate with a centered slope crack
A rectangular plate specimen with a centered inclined crack subjected to a pure thermal load is analyzed, with the dimensions 2L, 2W , the crack is defined with the half-length a and the slope is characterized by the β angle Fig. 4a. The displacements along the e yaxis is fixed at the bottom extreme right corner, and the bottom left corner is clamped. Both right and left boards are completely insulated, an imposed temperatures of ± 0 are defined at the top and bottom sides, such 0 = 10 • C. We consider a uniform enrichment disks radius in the case where several crack-tips exist; hence, E R ≡ E j R . The radius of disks   Fig. 4c. The objective of this example is, firstly, to study the accuracy of J-integral computation independently of the choice of Rice integral contour. Secondly, to show the robustness of computation of stress intensity factors in the dominated mode-II for various horizontal crack lengths. The stress intensity factors are normalized by dividing K II by α 0 E √ W which gives K norm II . In Table 2, the numerical normalized SIFs results are presented for four selected paths centered on right crack tip, referred by numbers '1' to '4'. There is no difference related to the choice of the right or left tip. The physical domain is discretized with a structured quadrilateral mesh with a characteristic length of 0.016 m. The variation of the SIFs values remains in the range [0, 0.9482%] with a maximum variation of 0.94% with respect to the minimum value, which corresponds to Path 2. The results obtained with XFEM-f.a. show a good outcome for path independence. Next, we consider different crack lengths starting from 0.1 to 0.6 with a jump of 0.1 with the same specimen configuration. The normalized SIFs results obtained with XFEM-f.a. agree closely with those presented by Murakami [54], Prasad et al. [44] and Duflot [27] for each crack length as shown shortly in Table 3 with respect to the results presented by the previous cited references. Complete results of this example are presented in Appendix: Table 7, with 6 digits, including the negative values of SIFs illustrating, as an indicator, of an important compressive aspect in the vicinity of the two tips. Temperature distribution is illustrated in Fig. 5a, as well as the e x , in Fig. 5b, and e y , in Fig. 5c, displacements. An important concentration of the stresses at the crack-tips are observed by Fig. 6a-c. Thermal flux is perpendicular to the crack Fig. 7a, b, since the crack is adiabatic; one notes a gradual deviation of the flux lines to circumvent the geometry of the crack Fig. 7c. These results are not presented by the references cited above. Second, as a benchmark problem, we treat the general case of any choice of β ∈ [0, π 2 ] and various choices of crack lengths. Dimensions is chosen such L/W = 0.5, inclined crack is defined by the total crack length 2a Fig. 4b. The main objective of this example is to show  Table 4 summarizes the results of normalized SIFs for a fixed angle β = 30 • and various crack lengths varying from 0.2 to 0.6. In Table 5, we give the normalized SIFs results for a fixed crack length,    Tables 8 and 9 with 6 digits, including again the negative values of SIFs. Results are observed to be in good agreement with Murakami [54], Prasad et al. [44] and Duflot [27]. Temperature distribution influenced by the prescence of crack, e x and e y displacements are presented respectively in Fig. 8a-c. Horizontal, vertical and line flux are plotted respectively in Fig. 9a-c. Stresses, Fig. 10a-c, show the same behavior around the crack-tips like in the case of the square plate. Again, these results are not presented by the references cited above.

Square plate with round hole and two cracks
In this example, we consider the case of a square plate with a hole and two cracks. The objective is: firstly, to show the influence of the presence of a hole, due to a manufacturing defect or willingly introduced into the material, on the stress intensity factors. Secondly, to examine the influence of radius of the fixed enriched zone on the SIFs. Thirdly, to show the influence of the characteristic length (h) on the convergence of SIFs when h goes to  Fig. 11a. The two cracks are defined at the two ends, right and left, of the hole are centered (right and left cracks), with a length l. Half-length of the apparent crack is defined by a = l + R. The bottom left corner is clamped and displacements along e y -axis is fixed. The heat flux is zero at the right and left edges, an imposed temperature of ± 0 are defined at the top and bottom sides, such 0 = 10 • C. We investigate the influence of various fractions of hole size R/L and cracks sizes l/L on the SIFs computations. We choose a set of R/L resp., l/L between 0.0 and 0.3, resp., 0.1 and 0.6 with a jump of 0.1. The structured mesh is used Fig. 11b, such the characteristic length is 0.011 m. It is worth noting that when R/L and l/L become too small, we refine sufficiently close to the two cracks to ensure a good approximation of the SIFs for the J-integral domain. The stress intensity factors are dominated by mode-II; we normalize it by α 0 E √ W which gives K norm II . The normalized SIFs for several choices of the two ratios are illustrated in Fig. 12. Results obtained by the XFEM-f.a. are close to those given by Prasad et al. [44]. Complete results of the two cracks are given in Appendix: Table 10. Influence of the radius of the fixed enriched area on the SIFs computation is inspected for several values of E R . A typical case are chosen for a hole with R/L = 0.1 and l/L = {0.5, 0.6} for both left and right cracks, Table 6. It can be seen that there is no significant difference in the computation of SIFs, for both left and right cracks, with respect to the choice of the radius value of the enrichment disk. A relatively large radius related (and independently) to the characteristic length of the mesh is desirable. As mentioned before, when the two a b

Fig. 11
Square plate with round hole and two cracks: a boundary conditions representation, b structured mesh used  ratios have become small, we tend to refine the mesh in the vicinity of the crack. Therefore, the selection of radius must be adapted to the size of the crack and the characteristic length. Convergence of SIFs computation has been demonstrated with respect to the size of the mesh, by comparing between a range of a coarse mesh and a sufficiently finer one. As illustrated in Fig. 13, XFEM-f.a. guarantees a significant convergence of the thermomechanical model and the SIFs computation. Distribution of the temperature generated by the effect of the crack and hole is presented in Fig. 14a. Displacements in e x and e y Directions are given respectively in Fig. 14b, c. We note an important concentration of the stress in the two crack-tips and around the perimeter of the hole in the vertical direction, as showed respectively in Fig. 15a-c for σ xx , σ yy and σ xy . The heat flux in e x and e y directions are figured respectively in 16a-c represents the spatial distribution of the flux lines that bypasses both the cracks and the hole. is linear in space and in time horizontally and remains uniform along the e y axis, so a typical cut over the line-section {x ∈ \ c ; y = 0} is presented for different increments until reaching the thermal equilibrium, as shown in Fig. 18a. On the left border, the plate tends to expand. Additionally, the displacement is fixed along e y -axis throughout the top and bottom sides; this generates a significant displacement, with respect to the pseudotime, of the upper crack surface towards e y and symmetrical displacement of the lower crack surface towards −e y . This behavior leads to a gradual opening of the crack with the continuous transient load until the achievement of the equilibrium state. Figure 18b depict the Euclidean norm of displacement combining the horizontal and vertical one over the line-section {y ∈ \ c ; x = 1 5 }. The stress, with respect to the pseudo-time, also becomes important near the crack tip after the incremental thermal loading.

Crack growth
This section shows an application of XFEM-f.a. in the thermo-mechanical growth of an edge crack governed by mode-I. The configuration is taken with the same considerations and initial crack as mentioned above. Crack growth is monitoring using hoop stress, by introducing an additional virtual crack extension (VCE) based on VCE-method Hellen [55], Millwater et al. [56] after each equilibrium step. The direction of discrete crack propagation is determined by the orientation in which the maximum energy is released from the system. Crack propagation was simulated for a total of 14 steps, with each step size of length a = 0.01 m. Convergence of SIFs is proved by Fig. 13; simply choose a sufficiently finer mesh to guarantee an optimal convergence, here characteristic length is 0.01 m. Determination of crack path is tested for several crack magnitudes a, 2 a, 3 a and 5 a and converges to the same path. Stress intensity factors K I , K II in Fig. 21a show a crack growth driven by mode-I, K II keeps a zero value for all steps. Stability of crack growth is illustrated by Fig. 21b, where the variation of energy release rate G remains negative for the total of increments. Crack progresses in a parallel direction to e x and at y = 0,

Rectangular plate with two circular holes and multiple cracks
A numerical example of a H × L rectangular plate, containing two circular holes and two cracks was designed. Physical and numerical framework ensuring a thermo-mechanical crack propagation is well defined and sketched in Fig. 23a. The width H = 0.5 m, the length L = 1.0 m and the radius of each circle is R = 0.07 m. Initial cracks are set to start from the limit of the left hole; the first crack termed 'crack 1' is a straight crack with a length a 1 = 0.05 m, the second crack termed 'crack 2' is an inclined crack with an angle β = 60 • and a length a 2 = 0.1 m. Displacements of the specimen are fixed vertically throughout the upper and lower part, excluding the two right corners which are embedded. The right edge is cold at − 0 and hot on the left one at 0 temperature, where 0 = 20 • C. The plate is completely insulated on the top and bottom borders. The computational domain is outlined in Fig. 23b, characteristic length of the mesh used is 0.012 m; we refine a bit more throughout the borders of the two circles. Material properties are defined in the Table 1, with the consideration of Young's modulus E and a specific choice of thermal expansion coefficient α = 10α. The radius E R of disks related to each crack tip is taken uniform and equal to 0.1m. No reference found in the literature dealing with this crack growth example for comparison. The configuration of specimen defined in this case, presenting two material defects (holes) and multiple cracks, allows simulating a thermo-mechanical crack growth by a b Fig. 24 Crack growth in plate with two circular holes and multiple cracks: stress intensity factors K I and K II ; a crack 1, b crack 2 XFEM-f.a. driven by mixed mode. Similarly, the VCE method is assumed for the progression of cracks. Furthermore, considering the two cracks on the border of the left hole is not an arbitrary choice. The idea is to identify the zone which has an important stress concentration enabling a progression 'opening'. Cracks move forward simultaneously with a combined propagation criterion for both of them. The growth was simulated for a total of 13 extension steps, with each step size of length a = 0.01 m. Several crack magnitudes and sufficient characteristic lengths of the mesh are chosen and converge to the appropriate crack path. Stress intensity factors in mode-I and II, K I and K II , of the two cracks '1' and '2' are presented respectively in Fig. 24a, b. Crack 1 is driven by a dominated mode-I from the beginning with a range of ∼ 10 6 that is sufficient to preserve a progressive growth. Whereas, the intensity of the driven magnitude at the vicinity of the initial crack-tip 2 is 10 times less than the crack-tip 1 and 2 times less than the intensity required to evolve the crack 2. This case brings to a crack 'initiation' controlled by both mode-I and II for the first extension and described by a drop, of ∼ 6 times less, of K II and rise, of ∼ 2 times more, of K I . Energy release rate G 1 of crack 1 and G 2 of crack 2, Fig. 25a, keep the same behavior as the stress intensity factors. The G 1 remains stable throughout the incremental progression; while crack 2 seeks to reach local stability for the first two extensions by getting the required energy to progress crack and subsequently maintain the stability. Crack paths are depicted in Fig. 25b including a representation of the initial preexisting cracks and the left round hole.
The crack arrest is taken for a state related to the configuration defined by 13 a. This represents the state where the crack 1 develops locally an important compression at the crack surfaces; it produced by the new material configuration caused by the displacement of the lower surface of crack 2 along the −e y axis. The specimen is thermally loaded throughout the process of crack propagation. Consequently, the temperature profile is affected by the new configuration of the cracks which causes an incremental change in the spatial redistribution of the temperature; the initial state and the final one are shown respectively in Fig. 26a, b.  Displacements are more pronounced over e y direction; Fig. 27a, b stand for respectively e y -displacement at initial state and at 13 a state. Stress combining the information of different directional stresses is given by the von Mises stress, σ VM , at the initial state in Fig. 28a and at the final one in Fig. 28b. Spatial redistribution of σ VM , induced by the temperature profile, varies with the crack growth. Stress becomes maximal near the crack tips and over the outer borders of the two holes at the initial state; it increases for the holes and decreases slightly at the tip points thanks to the release of energy caused by the crack progression.

Conclusions
A new thermo-mechanical crack propagation model in a cracked body was presented which can be applied, for instance, to ensure the safety of structures subjected to thermal loading. The developed geometrical eXtended finite element method was successfully applied to model crack growth and achieving the expected optimal rate of convergence by confirming the benefit of the fixed enrichment area approach on the computation of stress intensity factor profile. Numerical development and various matrices in full coupling were presented for each sub-problems, mechanical and thermal, and for the full coupled XFEM part. The criteria for crack growth, as well as for the direction of the virtual crack extension are described, and their performance in the context of the XFEM is discussed. From three examples, various benchmarks result in a cracked domain are examined and validated from the existing results in the literature. The robustness and the accuracy of the model implementation to extract the thermo-mechanical responses and to compute the associated stress intensity factors for stationary crack, with and without holes, as well as the effect of crack length and hole position on the SIFs are proved. Furthermore, a quasi-transient load example governed by mode-I is presented and the contribution of this loading on the profile of the SIFs until reaching thermal equilibrium is analyzed. Finally, an example of multiple mixed-mode cracks growth and multiple holes that may be present as small flaws in the material manufacturing stage is examined; only the limiting cases of stable crack are discussed. When the heat flow is distributed by the presence of the cracks, we observe a high local intensification of thermal gradients followed by an intensification of thermo-mechanical stress around them, which may lead to the crack growth or inevitable collapse of the structure.
As outlook of future works, possible improvement of this study can be made by taking into consideration the mechanical contact aspect between the crack surfaces; this will be important to extend to study of the last example to simulate the complete process of growth. Another point can be viewed by holding the crack propagation in the overall dynamic of the whole problem and admitting a crack-pseudo-time-dependant; which make it possible to control the evolution of the crack with the transient loading. This case requires a sophisticated treatment of the stiffness matrix; K needs to be evaluated at two different times for two different configurations of the crack, K i and K i+1 .