The concept of Representative Crack Elements (RCE) for phase-field fracture: transient thermo-mechanics

The phase-field formulation for fracture based on the framework of representative crack elements is extended to transient thermo-mechanics. The finite element formulation is derived starting from the variational principle of total virtual power. The intention of this manuscript is to demonstrate the potential of the framework for multi-physical fracture models and complex processes inside the crack. The present model at hand allows to predict realistic deformation kinematics and heat fluxes at cracks. At the application of fully coupled, transient thermo-elasticity to a pre-cracked plate, the opened crack yields thermal isolation between both parts of the plate. Inhomogeneous thermal strains result in a curved crack surface, inhomogeneous recontact and finally heat flow through the crack regions in contact. The novel phase-field framework further allows to study processes inside the crack, which is demonstrated by heat radiation between opened crack surfaces. Finally, numerically calculated crack paths at a disc subjected to thermal shock load are compared to experimental results from literature and a curved crack in a three-dimensional application are presented.


Introduction
A crack contact criterion and the decomposition into degraded and transferred forces through the crack are approximated by splits of the deformation energy potential in phase-field models for fracture. Unphysical predictions of deformations at cracks are reported for common splits like the volumetric-deviatoric, the spectral and similar splits. Furthermore, the crack driving force of the phase-field variable is significantly controlled by the split approach and, thus, may result in unreliable predictions for the fracture process.
Artificial splits are replaced by Representative Crack Elements (RCE) in [1], which allows to overcome the kinematical deficiencies. After first application to anisotropic elasticity [1], the framework is successfully applied to phasefield fracture for visco-elastic materials [2], for inelastic materials, crack surface friction, finite deformations [3] and Michael.Kaliske@tu-dresden.de 1 Institute for Structural Analysis, Technische Universität Dresden, 01062 Dresden, Germany cohesive interface failure [4]. In the paper at hand, representative crack elements are developed towards phase-field fracture for transient coupled thermo-mechanics. The variational formulation is presented and derived into a finite element model. Through academic examples, the potential of representative crack elements is demonstrated for multiphysical applications and for processes, which take place inside the crack, e.g. heat radiation through the crack. It is also shown in the theoretical discussion, that the numerical solution scheme, which is introduced to solve the RCE model with heat radiation, can be used to model the reduced heat conductivity in closed cracks as a consequence of micro cavities, and to model heat convection at the crack surface.

Theory
The principle of total virtual power δP tot = δP int − δP ext = 0 balances the virtual power, which is associated with a thermodynamic system (internal and kinetic energy) and the virtual power transferred to the system from the exterior (mechanical work and heat) for all admissible generalised virtual velocities ∀δv ∈ Var v and their generalised local strain rate actions D(δv). Internal and external virtual pow-ers are formed by dual products of the tensor valued virtual velocity components δv i and strain rate components D i (δv) with their thermodynamically conjugated forces f i and stresses i , which yield where H is the displacement gradient and the dot represents time derivative.

Constitutive relations
It is shown in [5] that the variational formulation in the form of Eq. (3) represents the weak form of the momentum balance, the entropy balance and von Neumann boundary conditions. Momentum rate u , 1. Piola-Kirchhoff stress H , entropy rate θ , negative entropy flux ∇θ and contributions to the phase-field driving force p , ∇ p can be identified without additional assumptions on the constitutive relations. The variation of entropy balance in the form of [6] is obtained by accepting the relation of entropy flux to heat flux θ h = q, and entropy source to heat source θ d = r . The thermodynamically conjugated thermal stresses and mechanical stresses can be derived evaluating the first and second law of thermodynamics using Helmholtz free energy ψ. Fourier's empirical law of heat conduction states the relations between temperature and heat c = −θ ∂ 2 ψ ∂θ 2 , q = −K · ∇θ (6) and is experimentally validated at various thermodynamic systems. Heat capacity c and thermal conductivity tensor K are material parameters. 1. Piola-Kirchhoff stress H is known from standard thermo-mechanics and p together with ∇ p form the driving force for the phase-field. The influence of the presence of a crack on the mechanical and thermal behaviour is considered through the postulation of a constitutive structure in the form The degradation function g( p) acts as interpolation weight between the material behaviour of intact material, represented by ψ m,0 and q 0 , and fully degraded material, ψ m,c and q c , as introduced in [1]. The regularised fracture energy is given by ψ ph . ψ th is the stored heat in the material. The reader is referred e.g. to [6] for further discussion on ψ th at dissipative deformation processes.

A thermo-mechanical representative crack element
Whereas, constitutive laws for intact mechanical material ψ m,0 , the regularised fracture energy ψ ph and the heat storage ability ψ th are defined using standard constitutive models, the behaviour for cracked material ψ m,c and q c is determined by means of a Representative Crack Element, see Fig. 1. A local coordinate system (N 1 , N 2 , N 3 ) is introduced, where N 1 is the normal direction of the crack surface. The crack is modelled as discrete crack by the RCE, hence, the state variables are displacements and absolute temperature. Coupling operators between the cracked part of the phase-field model and the representative crack model are defined bȳ following the concept of variational homogenisation of [7]. Quantities denoted by bar are associated with the RCE model, Bx dV is the geometric centre,ū andθ are the unknown displacement and temperature fluctuations, and V is the RCE volume. Homogenisation of the state variables is given by Once, a solution for the unknown fluctuation fieldsū and θ is determined, the thermodynamically conjugated forces, stresses and heat flux can be calculated as wheref u are volume forces andf θ volume heat sources.
These relations are a result of the principle of multi-scale virtual power [7]. A solution for the unknown displacement and temperature fluctuations of the RCE can be obtained numerically, for instance by application of the finite element method. The embedding of finite element problems into the material model of another finite element problem is called FE 2 and is computationally very expensive. That is why an analytical and a semi-analytical solution are derived for linear and nonlinear bulk materials in the following. The computational costs of the resultant phase-field model are similar to standard phase-field models. Both RCE solutions are based on the assumptions of • homogeneous displacement and temperature gradients in B 1 andB 2 , • weak thermo-mechanical coupling through homogeneous thermal expansion, and • homogeneous volume forcesf u ,f θ and transient stresses u ,¯ θ .
Based on these simplifications, the RCE can be divided into a quasi-static mechanical problem for the determination of the displacement discontinuity vū w and a steady-state thermal problem for the determination of the temperature discontinuity vθ w.
The displacement field of the mechanical RCE for small deformations reads The linearised strain tensor is E = 1 2 (H + H ). The normalised relative displacement of the crack surfaces is called crack deformation i = vū w i /l 1 . The solution for the mechanical Representative Crack Element with thermal strains is presented in [1] for the first time. Thermal strains at the RCE are calculated based on the average temperature θ | x in order to ensure homogeneous strain fields.
The thermal RCE for the heat conduction is determined with the absolute temperature field where A solution for the unknown temperature gradient through the crack θ = vθ w/l 1 follows from the principle of virtual power of the RCE The linearisation and discretisation of the RCE with the nodal vector, shape function and the derivative of the shape function result in the linear system of equations for the nodal updates where θ | x and ∇θ | x are the quantities at the material point of the phase-field model and given by Dirichlet boundary conditions to the RCE.

Solution for thermally isolated cracks
The first analytical solution for the RCE model is derived under the assumptions of perfect heat conduction through a closed crack and that an opened crack behaves as thermal isolator, i.e. no heat flux through the crack exists (¯ ∇θ (x) = 0, ∀x ∈∂ B , for 1 > 0). As introduced in [1], linear thermo-elastic material at small deformations with the Lamé constants λ, μ and the coefficient of thermal expansion α yield where the thermal strain tensor is E θ = αθ| x I. For the kine- The solution for the temperature discontinuity reads which results in the homogenised thermal stress and tangent for an opened crack ( 1 > 0). The thermal conductivity tensor for an isotropic material is K = k I. The temperature gradient of both solid blocks of the RCE reads∇θ = ∇θ | The solution of the mechanical and thermal RCE depends on the crack contact state, which is identified by 1 . A tolerance range is adopted from [3] in order to avoid contact oscillations. However, still some additional Newton iterations due to contact oscillations are observed for the thermo-mechanical formulation with a tolerance range in the subsequent examples. The reason is the bilinear character in both, mechanical and thermal behaviour, due to crack contact. Thus, the contact criterion is fixed in each time step after the first Newton iteration, which results in quick and robust convergence.

Solution for cracks with heat radiation
Heat transfer between bodies through radiation is discribed in [8]. This process depends on the absolute temperature at the body surfaces, the surface areas, the emissivity ε, the Stefan-Boltzmann constant σ B and the geometrical view factor F. The emissivity relates the radiation absorption ability of the physical concept grey body (can reflect radiation) to those of a black body (no reflection). The view factor is the fraction of energy which leaves the first surface and reaches the second surface.
Crack surfaces are considered to be equally parallel planes and the radiation is approximately related to the average surface temperature in the following model. This simplification allows to avoid a numerical integration on the crack surfaces in the RCE for the calculation of radiation heat flux. Then, the heat flux reads The average temperature of the crack surfaces follows from Eq. (20) Following [9], the view factor of identical parallel square facets of length l 2 and distance l 1 1 is given by Then the iterative Newton-Raphson scheme for the unknown temperature discontinuity reads The residuum is and the algorithmic tangent is A numerical derivative based on finite differences is implemented for the view factor derivative for the purpose of simplicity. The homogenised thermal stress and tangent read for an opened cracked ( 1 > 0). Note, that heat radiation and, thus, the solution of the RCE, depends on the value of the crack opening vū w, whereas previously published solutions of the RCE depend only on the dimensionless crack deformation i . Thus, in contrast to the former RCE solutions, the model is not independent of the RCE size and l 1 becomes a model parameter.

Remarks on further heat flux processes at cracks
Perfect heat conduction behaviour through a closed crack is assumed for both solutions above, i.e. θ = 0 for 1 = 0. However, the heat conduction in compressed cracks is typically less than in the bulk material due to spot contact, conduction in micro cavities and further effects. [10] has identified the following constitutive structure for the heat flux through a closed crack where q N is the heat flux normal to the crack surface and p N is the normal pressure of the crack surfaces in contact. Thus, the iterative solution scheme for the temperature discontinuity at the RCE model with heat radiation ( 1 > 0), given by Eqs. (33-35), (37) and (38), also applies to the proposed RCE model when Eq. (39) is considered for closed cracks ( 1 = 0). It is further possible to consider heat convection at the crack surface, which depends on the temperature of the crack surfacesθ 1 andθ 2 , and the temperature of the medium in the crackθ c . Then, the heat flux at the crack surface can be given in the form However, the medium temperature in the crack depends on the heat conduction and fluid mechanics of the medium in the whole crack in general, which is not part of the proposed phase-field model. Thus, heat convection at crack surfaces is restricted to slow processes, whereθ c = const. can be assumed, or very quick processes, whereθ c can be calculated from the heat capacity of the medium in the crack and the heat transferred through the crack surfaces.

Finite element formulation
The finite element method is applied to the phase-field model in Eq. (3). Thus, a solution is obtained iteratively starting at the state (u i , p i , θ i ) with and the linearisation by means of the first order Taylor expansion with respect to the state variables. For the discretisation into finite elements, the same shape functions N (x) and derivatives B(x) are used for the state variables and their corresponding virtual velocities The total virtual power is required to vanish for all admissible virtual velocities, which yields the linear system of equations for the nodal updates of the state variables A n e=1 is the assembly operator for the n finite elements towards the global system of equations.

Pre-cracked plate
The proposed thermo-mechanical phase-field model for fracture is applied to a pre-cracked plate. For a plate of 100 mm length with an initial middle crack, a displacement of u y = 4 mm is applied to open the crack, cf. Fig. 2a. The temperature at one edge parallel to the crack is increased by 6 K. The load history of the boundary conditions is given in Fig. 2b. The time is discretised into 100 time steps. The thermoelastic material parameters are λ = 19.6 MPa, μ = 2 MPa, α = 0.01 mm/K, ρ = 0.1 kg/mm 3 , c = 38.8 J/kg/K and k = 1000 W/m/K, while thermal isolation is assumed for opened crack states.
The temperature distribution on the deformed plate is shown in Fig. 3c, and the temperature and displacement profile at the vertical symmetry axis (x = 50 mm) are given in Fig. 3a, b. The heat propagates from the plate edge towards the crack and causes thermal expansion. The heat cannot pass the opened crack and yields a curved crack surface, which propagates towards the opposite, straight crack surface. The first crack contact takes place at the middle of the plate and allows the heat to propagate through the region in contact into the second half of the plate. The ongoing thermal expansion in both halves of the plate increases the contact area of the crack until contact along the entire crack surface.
The example demonstrates realistic deformations at the crack for thermo-mechanical material behaviour and the influence of the crack state on the thermal conduction, while the opened crack is considered as thermal isolator. Convergence is obtained in three or less Newton iterations for all time steps.
The simulation of the pre-cracked plate is repeated and heat radiation between the crack surfaces is considered, cf. Fig. 4a. The boundary condition loads are changed to u y = 7 mm and 6 K in order to prevent crack contact when only one plate half is thermally expanding. The emissivity is considered as ε = 1, which represents heat radiation for black bodies. The heat is propagating from the plate edge to the crack. However, the temperature difference between the crack surfaces increases slowly and, thus, only a small amount of heat is transferred by radiation. With increasing temperature at the crack surface, also the distance between the crack surfaces reduces, which increases heat flux from radiation. The slightly curved crack yields an inhomogeneous distribution of the heat flux at the crack. The second half of the plate is heated up faster at the middle than at the plate sides. The temperature distribution on the deformed plate is given in Fig. 4c for a simulation with and without heat radiation at the opened cracked. The thermal expansion of the lower half plate yields crack contact after a while, which cannot be achieved without taking heat radiation into account.
The solution scheme of the non-linear problem consists of two Newton-Raphson loops, one applied to the global system of equations and one at the RCE model. The convergence behaviour in terms of the number of iterations is presented in Fig. 5. Up to nine iterations per time step at the global level are required due to the strongly non-linear radiation relation. Corresponding relative residuum norms of the iterations are given for three time steps at all phases of the simulation. The local Newton loop for the radiation flux is plotted for an element in the centre of the plate in Fig. 5b and, thus, at the centre of the crack. Convergence is achieved within two to three iterations. Quadratic convergence is demonstrated by the relative residuum norms of three time steps. The local loop is not activated when the crack is in contact, which is the case from step 26 onwards for the investigated element.

Thermal shock of a circular specimen
The second example examines the capability of the present phase-field model applied to a practical engineering problem. A recent experiment of [11] attempts to investigate crack evolution patterns of ceramic materials exposed to a thermal shock (quench treatment), which obtains a quasi periodical and hierarchical crack profile within a thin layer of circular specimen. Thereafter, a numerical method [12] using a thermal coupled phase-field fracture model studies the cracking pattern, showing good agreement with the experimental investigation in [11]. The work at hand reproduces this model problem using the proposed thermo-phase-field approach within an RCE framework based on the setup shown in Fig. 6.
In particular, a two-dimensional boundary value problem, which is characterized as homogeneous and isotropically linear elastic, is taken into account. For the initial phase, the specimen is uniformly heated up by 200 K, where the whole heating procedure yields a complete stress-free state of the specimen. In the sequel, the most outside surface is supposed to be cooling down in a significantly fast manner. Naturally, this process leads to material shrinkage at the outside layer and results in tensile stress along the tangential direction.
As long as the elastic strain potential sufficiently exceeds the critical energy release rate of the material, e.g. G c = 8 × 10 −3 GPa/m 2 , multiple cracks are nucleated and start to propagate inwards along the radial direction. Due to the symmetry characteristics of the geometric setup and the temperature loading condition, the model problem is simplified as a quarter of the original one with proper boundary conditions, depicted in Fig. 6a. The simulation is performed by a standard parallel computation technique, where the finite element discretization is partitioned by 80 individual domains, see Fig. 6b. The material parameters are given as κ = 196 GPa, μ = 20 GPa, ρ = 1000 kg/m 3 , k = 50 W/m/K, c = 3.88 × 10 −3 kJ/kg/K,

Thermally induced cracking of a notched beam
This example studies the crack evolution induced by thermal shrinking within a beam with an initial notch. Similar to the aforementioned thermal shock cracking mechanism, sudden cooling is applied to the surface of the structure. Thus, the thermal shrinking yields a stress concentration at the notch tip, where the crack is supposed to be initiated.
In order to demonstrate the propagation of curved cracks in a three-dimensional specimen, an inclined notch is taken into account. The geometrical setup of the boundary value problems is depicted in Fig. 8. The material parameters are given as κ = 196 GPa, μ = 20 GPa, ρ = 3000 kg·m −3 , k = 520 W/m/K, c = 1.66e −5 kJ/kg/K, α = 1.2e −4 mm/K, G c = 5e −2 GPa/m 2 , l = 2 mm. The maximum principal stress direction is used to approximate the crack normal orientation [1,3]. The temperature loading is applied to the bottom surface, where the initial notch is located at, as a Dirichlet boundary condition. In the sequel, due to continuous heat conduction, the lower temperature propagates upwards reaching the notch tip. As a result, a concentrated elastic deformation occurs at the notch tip to initiate the phase-field crack. The inclined notch yields a curved crack surface, which reduces the crack surface area while dividing the beam into two parts, see Fig. 9. For the sake of better visualization, a post-processing technique generates the isosurface with p = 0.95 to represent the evolved crack surfaces.

Conclusions and outlook
The abstract variational principle of total virtual power is applied to phase-field fracture for coupled thermo-mechanics. The state variables are chosen to be the displacement vector, the phase-field variable and the absolute temperature. The obtained variational problem represents the weak form of momentum balance, entropy balance and phase-field balance. A representative model for a portion of a crack is introduced in order to determine crack contact, force degradation, degradation of the heat conductivity and further processes inside the crack. This Representative Crack Element (RCE) is coupled to the phase-field formulation by means of variational homogenisation.
The presented framework is used to derive a phase-field fracture model for transient thermo-elasticity, where the crack acts either as thermal isolator or where heat radiation through the crack is considered. An iterative solution scheme to solve the non-linear RCE model is introduced. The possible application of the solution scheme to considered further thermal effects, e.g. spot contact, micro-cavity conduction or heat convection, is discussed. The derived models are applied to a plate with an initial crack, demonstrating the interaction of thermal strain and heat conduction with the crack with and without considering heat radiation. A thermal shock simulation is presented and the results are compared to experimental findings. The three-dimensional application is performed at a clamped beam with an inclined notch, which yields a curved crack surface. The presented procedure to derive multi-physical models of phase-field fracture based on representative crack elements can be adopted to further formulations in future works, e.g. nonlocal material models, higher order continuum models, multi-scale problems and other multi-physical couplings.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.