A Fully Coupled Thermomechanical Phase Field Method for Modeling Cracks with Frictional Contact

: In this paper, a thermomechanical coupled phase ﬁeld method is developed to model cracks with frictional contact. Compared to discrete methods, the phase ﬁeld method can represent arbitrary crack geometry without an explicit representation of the crack surface. The two distinguishable features of the proposed phase ﬁeld method are: (1) for the mechanical phase, no speciﬁc algorithm is needed for imposing contact constraints on the fracture surfaces; (2) for the thermal phase, formulations are proposed for incorporating the phase ﬁeld damage parameter so that different thermal conductance conditions are accommodated. While the stress is updated explicitly in the regularized interface regions under different contact conditions, the thermal conductivity is determined under different conductance conditions. In particular, we consider a pressure-dependent thermal conductance model (PDM) that is fully coupled with the mechanical phase, along with the other three thermal conductance models, i.e., the fully conductive model (FCM), the adiabatic model (ACM), and the uncoupled model (UCM). The potential of this formulation is showcased by several benchmark problems. We gain insights into the role of the temperature ﬁeld affecting the mechanical ﬁeld. Several 2D boundary value problems are addressed, demonstrating the model’s ability to capture cracking phenomena with the effect of the thermal ﬁeld. We compare our results with the discrete methods as well as other phase ﬁeld methods, and a very good agreement is achieved.


Introduction
Thermomechanical contact cracks are ubiquitous in engineering applications, ranging from metal forming and powder compaction processes [1] to heat conduction across the contact interfaces in cooling systems of microelectronics [2]. It is essential yet challenging to accurately model the frictional crack phenomena, including fracture onset, propagation, and branching, particularly in a multi-physical environment [1,3,4]. To tackle the aforementioned complex physical phenomena with regard to cracks and fractures that cannot be resolved using analytical methods, developments within the context of computational mechanics and numerical simulations have been a matter of intensive research in recent years.
In general, thermomechanical contact constraints are treated using methods such as the Lagrange multiplier method and augmented Lagrange multiplier method [5][6][7][8], or the penalty method [9], which are inherited from the discontinuous Galerkin (DG) method. In these methods, the contact constraints are either enforced through the Lagrange multipliers or by penalty regularization. When using the Lagrange multiplier method, an inf-sup stability requirement needs to be satisfied with the additional unknown degree of freedom [6,7]. On the other hand, for the penalty method, while it remains primal and is easy to implement, a carefully chosen penalty parameter is crucial to avoiding the loss of coercivity [10] and ill-conditioned systems with too many constraints [11]. The Nitschetype method, originally proposed to enforce element boundaries in a weak sense [12,13], is another method used to model contact and frictional conditions [2,10,14,15]. Different from the Penalty-type method, the Nitsche-type method is generally consistent without penetration. One advantage of the Nitsche-type method over the Lagrange multiplier method is that it contains appropriate and consistent terms that only involve the primary fields, i.e., there is no need to fulfill the inf-sup condition [10]. Despite all the advantages of the Nitsche-type method, the interface weighting parameters and stability parameters are user-defined with special treatments. To resolve this issue in Nitsche's method with regard to user-defined parameters, the variational multiscale discontinuous Galerkin (VMDG) method [16,17] has been proposed to provide a primal formulation with interfacial stabilization for both weak and strong discontinuities using the variational multiscale (VMS) technique. In the VMDG method, the formulation is naturally derived, and no ad hoc parameters are appended. Nevertheless, in all these interfacial element formulations, continuity restrictions are imposed at the junction of overlapping finite element grids, and all methods exhibit element-wise conservative approximations [14,18] for cracks. Moreover, the discontinuities are aligned with element boundaries. Due to the complexity of the crack patterns, meshes should be fine enough around the crack region to ensure arbitrary crack branching [18]. On the other hand, allow the discontinuity to be embedded inside elements, local remeshing techniques, such as the advanced methods based on X-FEM [19][20][21] or E-FEM [22][23][24], are used near crack regions for sharp crack propagation. However, these local enrichment strategies suffer in three-dimensional (3D) applications [3,25,26] and present limitations when predicting crack initiation, branching, and coalescence for multiple crack fronts [27]. In addition, in the context of fracture mechanics, modeling the structural behavior of multiaxial fatigue problems is challenging. The experimental results show that it is essential to model the fracture surface based on qualitative fractography [28] and to understand the fracture surface response based on fracture surface topography parameters [29]. In the work recently published by Lanwer et al. [30], a combined experimental and computational method is proposed to investigate the fatigue and degradation behavior of the fiber-matrix composite. In their work, a mechanism-oriented bonding model on the inserted fiber-matrix interfaces is developed for FE analysis. However, to better understand the fiber-matrix bond mechanism and the possibility of cracks propagating in the matrix, it is desired to have a numerical method where the crack path is not defined a priori.
The phase field method is a feasible algorithm proposed in recent years that provides a new perspective toward simulating complex crack topologies [3]. Different from the discrete methods such as the above-mentioned VMDG, Nitsche's, penalty, or Lagrange multiplier methods, the phase field method is rooted in classical energy-based Griffith's theory [26,31]. The sharp crack surface is regularized within a thin band of a diffusive crack zone governed by a scale variable, called the damage parameter, that distinguishes the damaged and undamaged domains [3,25]. This nonlocal damage variable is governed by a Poisson-type partial differential equation (PDE) and evolves with crack propagation. Compared to discrete methods, the phase field method has decisive advantages; the explicit crack interface tracking is unveiled by this scale variable [32], and the numerical implementation is straightforward [3]. This method has gained extensive attention in recent years, making remarkable contributions to modeling quasi-static and dynamic brittle fractures [4,26,33,34], ductile fractures [35,36], fatigue [37][38][39][40][41][42][43], and multi-physics applications [32,[44][45][46]. For modeling frictional cracks, because preexisting phase field models for fractures only consider non-contact or stick contact conditions without incorporating a frictional slip condition, Fei and Choo [47,48] proposed a stress decomposition scheme within the context of the phase field method to model frictional slip condition. In their approach, unlike discrete methods, there is no need to implement specific tracking strategies for the interface/cracking surface and for the evolution of the gap function along the contact surface. Their approach, although very promising, only works for pure mechanical prob-lems. Other phase field methods for multi-physics applications [32,[44][45][46] only degrade or maintain the effect of thermal conductivity based on the damage variable, which is not sufficient to accommodate different thermal conductance conditions at the contact surfaces. In discrete modeling, the contact heat flux is characterized by a pressure-dependent thermal contact conductance as well as by the temperature jump across the contact interface [49,50]. A similar treatment should be developed within the context of the phase field method. To the best of the authors' knowledge, there is no such phase field model that accommodates different scenarios for thermal conductance.
In the present work, following the work done by Fei and Choo [47,48], we extend the pure mechanical phase field formulation to a thermomechanically coupled formulation for frictional contact problems. The proposed framework has the following features: (1) the contact constraints at the interfaces are represented without any explicit functions or algorithms, and (2) different thermal conductance conditions, including the pressure-dependent thermal conductance, are accommodated via the existence of the contact pressure and the phase field damage parameter. In other words, the key idea of the proposed method is to incorporate the thermal conductance constraints through the suitable calculation of the thermal conductivity in the fracture surface region. In addition to the proper calculations of the stress tensor for different contact behaviors and constraints for the mechanical phase [47], a fully coupled thermomechanical formulation for frictional contact problems is derived in the current work. The proposed phase field method is able to handle the complex crack geometry and the contact condition in the mechanical phase while coping with the different thermal conductance conditions in the thermal phase.
In the following sections, we first briefly introduce of the phase field method for coupled thermomechanical problems in Section 2. In Section 3, we develop the governing equations and corresponding weak forms for coupled thermomechanical problems that involve frictional contact. In Section 4, following a summary of the approach used to explicitly generate the stress tensor at the contact interfaces based on the different contact conditions, we introduce a new strategy to calculate the thermal conductivity within the interface region. A set of numerical simulations verifying the proposed algorithm is shown in Section 5, followed by concluding remarks in Section 6.

Phase Field Method for Coupled Thermomechanical Problems
The purpose of this section is to develop a new phase field formulation that accommodates both thermal and mechanical constraints within the fracture surface region. The formulation is derived along the same lines as [48] to model general frictional interfaces in solids, which we extend to incorporate thermal conductance constraints. This formulation is useful for modeling frictional contact problems in a thermomechanical setting.

Regularized Variational Framework
We first establish the regularized variational framework for thermomechanical problems in developing governing equations for the phase field problem. As shown in Figure 1, let Ω be the open domain and Γ I the curve for the possible cracks. In a regularized framework, the scalar parameter d in the fracture surface energy equation is used to approximate the crack geometry, while l c is interpolated to control the thickness of the crack in its regularized form. The scalar parameter d ∈ [0, 1] is incorporated to denote the damage crack region (d = 1) and the intact bulk domain (d = 0).
The variational approach for fracture mechanics proposed by Francfort and Marigo [18] states the energy function for the cracked body [25]. Based on the energy functions from [18], the regularized form of the energy describing the cracked structure induced by the thermomechanical effects can be expressed as follows: where E e (u, θ, d) is the elastic energy stored in the cracked body, E d (d, ∇d) is the fracture surface energy, u is the displacement field, and θ = T − T ref is the relative temperature field of the current temperature T with respect to the reference temperature T ref . In the standard linear theory of elasticity for isotropic solids, the global energy storage functional is the sum of the elastic energy and the fracture energy: where Ψ(ε, θ, d) is the energy density of the bulk and G c is the Griffth's critical energy release. Assuming the quadratic free energy functional, the undamaged strain energy functional Ψ(ε, θ) can be expressed as where ε is the strain, C is the fourth-order tensor representing the material moduli for the bulk domain, c v = ρ 0 c p is the volumetric heat capacity (with ρ 0 being the density and c p the specific heat capacity), and m is the second-order thermomechanical coupling tensor, which is discussed in detail in Section 3. In Equation (2), G c is the Griffth's critical energy release, with Γ I G c dΓ I = Ω G c γ(d, ∇d) dΩ measuring the critical fracture energy; γ(d, ∇d) is the crack surface density functional, and we adopt the most common one, which is written as the convex function

Phase Field Approximation
It is noted that the crack surface behavior is governed by the unilateral conditions of stress carrying capacity, depending on whether the solid is in tension or compression [1,51]. In the present work, we follow the idea from Miehe [3]. Considering the unilateral contact, the damage is assumed to only modify the positive part of the strain energy, with the positive part defined in terms of the principal strains at which spectral decomposition of the strain tensor is required. Based on this, the free energy density function in Equation (3) can be rewritten as Ψ(ε, θ, d) in the bulk to accommodate the damage parameter d, then decomposed to the positive part and negative part, as shown in Equation (5).
The existence of the small positive variable k is used to ensure that the partly broken systems are well-posed when numerical discretization is applied [3]. In the following applications, we set k = 10 −8 . A simple degradation function that describes the degradation of the tensile part of the stored energy evolving damage is interpolated as It is assumed to have the following properties: (1) g(d = 0) = 1 represents the undamaged phase; (2) g(d = 1) = 0 represents the fully damaged phase; (3) g (d = 1) = 0 guarantees that the energetic fracture force converges to a finite value when the fully damaged phase is reached [3]. Here, Ψ + (ε) and Ψ − (ε) are the strain energies decomposed from the positive and negative components of the strain tensor.
In order to define the positive and negative parts of the stored energy, the strain field is decomposed into the tensile and compressive modes, as follows: with ε i defined as the principal strain, and n i as the principal strain direction, and where x ± = (x ± |x|)/2. Using the definitions above, the positive (tensile) part Ψ + (ε) and the negative (compressive) part Ψ − (ε) of the stored energy due to the mechanical loadings are decomposed as follows: By substituting Equation (4) and Equation (5) into Equation (2) and taking the variational derivative with respect to the damage parameter d, the governing equation for the phase field problem to be solved by evaluating the phase field d(x, t) at time t is summarized as follows: ∇d(x, t) · n = 0 on ∂Ω (10) where H is the strain history function, G c is the fracture energy, and n is the outward unit normal to ∂Ω. As stated in [48], H can be arbitrarily defined as a constant for models with stationary cracks, meaning that the phase field equation only needs to be solved once in order to find the distribution of the phase field parameter. When dealing with crack propagation, H needs to be calculated based on the updated stress tensor and the phase field parameter. For the calculation of the strain history functional H, we adopt the idea from Miehe et al. [3], as follows:

Governing Equations and Corresponding Weak Forms
The balance of the linear momentum, as well as the balance of the energy of the problem, are provided below: where σ is the Cauchy stress tensor, ∇ · (·) = tr[grad(·)] is the divergence operator, ρ > 0 is the density, b is the mass-specific body force, c p is the specific heat capacity, θ = T − T ref is the temperature difference of the current temperature T with respect to the reference temperature T ref ,θ is the time derivative of the temperature field, q is the heat flux vector, v is the velocity of the mechanical field, r is the heat source, and m = 3κβI is the second-order thermomechanical coupling tensor, where κ is the bulk modulus and β is the coefficient of thermal expansion. Using Fourier's relation, the heat flux vector is defined as q = −k∇θ; by assuming isotropic, k = kI is defined as the thermal conductivity.

Remark 1.
It should be noted that Equations (12) and (13) are fully coupled through the existence of the second-order thermomechanical coupling tensor m.
The boundary conditions of this coupled problem are provided by whereū,t,θ, andq are the prescribed displacement, traction, temperature, and heat flux, respectively, and n is the unit normal vector pointing outwards at the boundaries shown in Figure 1.
The corresponding weak forms are summarized as follows by multiplying the test functions {η, δθ, δd}, integrating by part, and applying the divergence theorem. The formal statement of the weak forms is as follows: (20) are satisfied for all {η, δθ, δd}.
Weak form of the phase field: Weak form of the balance of momentum: Weak form of the balance of energy: In the context of phase field modeling with the expression of g(d) mentioned in the previous section, in Equation (19) we can express the stress tensor in the domain as a combination of the stress for the bulk and the interface [48], as follows: where σ bulk and σ interface are the stresses in the interior and on the interface, respectively. Equation (21) is considered a generalized expression of the stress in the phase field model for the frictional contact problem [47]. The bulk Cauchy stress tensor σ bulk is represented as σ = C : − mθ, where C is the fourth-order tensor representing the elastic modulus of the bulk material and is the strain as defined in Section 2.1. The calculation of the interface stress σ interface is discussed in Section 4 under different contact conditions.
In a similar fashion, we propose a generalized expression of the thermal conductivity k in Equation (20) for the coupled thermomechanical problem as where k bulk = k bulk I is the thermal conductivity of the interior domain, k interface = k interface I is the thermal conductivity of the interface domain, and A is a triggering value that distinguishes the bulk from the interface domain. We discuss how to decide this triggering parameter A and provide the expressions for the bulk and interface thermal conductivity k bulk and k interface in Section 4.

Stress Tensor σ Updates in the Phase Field Modeling
For the mechanical phase, following [47], updates to the Cauchy stress are made to incorporate the contact-dependent mechanical response of the interface based on different contact conditions.
We first summarize the stress tensor updates of σ interface in Equation (21) based on the different contact conditions. The normal strain is calculated to determine the following debonding/contact conditions: For the non-contact condition with ε N > 0, the stress free state results in the interface stress tensor σ interface = 0, and there is an open crack between two surfaces. Therefore, the stress tensor in Equation (21) is expressed as σ = g(d)σ bulk .
For the in contact condition with ε N ≤ 0, the yield function f (t T , t N ) is introduced to distinguish the stick and slip conditions. The yield function f (t T , t N ) with a frictional coefficient µ f for the frictional Coulomb model is defined as where t T = (I − n ⊗ n)σ bulk and t N = σ bulk : (n ⊗ n).
When the interface is considered in the stick condition, there is no relative motion between the bulk and interface domains [47], and σ bulk = σ interface yields the stress tensor in Equation (21) as σ = σ bulk . On the other hand, when the interface is considered in the slip condition with f (t T , t N ) ≤ 0, the interface stress tensor is calculated by decomposing into a frictional part σ friction and a no-penetration part σ no-penetration , as proposed by Fei and Choo [47], with the expressions as follows: When substituting Equations (25) and (26) into Equation (21), the overall stress tensor under the slip condition is obtained as The procedures for updating the interface stress tensor σ interface and the overall stress tensor σ are summarized in Box 1, where x int represents the geometry description of an integration point. Box 1. Stress update algorithm for phase field modeling of the frictional contact problem.
Step 1: Inputs at x int : strain ε, phase field variable d, the normal unit vector n, and the tangential unit vector m.
Step 2: If the phase field variable d = 0, THEN Bulk domain,

Thermal Conductivity Tensor k Updates in Phase Field Modeling
To the best of our knowledge, there is no specific formulation of thermal conductivity in the existing literature on the phase field method for coupled thermomechanical problems when considering frictional contact. In most existing phase field models, the thermal conductivity k is degraded as g(d)k when considering the damage effect. To accommodate frictional contact in the thermal field, we consider the heat transfer associated with the fracture damage as well as different types of thermal contact conditions. To this end, we present thermal conductivity k updates based on four thermal conductance contact models.
As shown in the Figure 2, we enlarge the domain near the crack region in Figure 1, and the crack is shown in Figure 2a. Using the phase field method, the crack is modeled by the thin diffusive damage band shown in Figure 2b. As provided in Equation (22), the trigger parameter A is used to distinguish the bulk and interface domains. By setting the threshold to the damage value, the distribution of the trigger parameter A is shown in Figure 2c. The procedure for updating the thermal conductivity for different regions is explained in Box 2; the value of A is determined by comparing the phase field variable d with the threshold, then jumping into the corresponding thermal conductivity updates.

Remark 2.
Recalling the phase field governing equation Equation (8), the analytical solution can be written as d = e −|x|/l c , where x is the distance from the current location to the nearest d = 1 point. This argument is strong evidence for defining the damage threshold. From this equation, the width of the diffusive band is determined by the length parameter l c . We set the threshold as d = e −1 ≈ 0.378, where x = l c , to divide the bulk and interface regions. When the damage value at the integration point surpasses this value, we assume that the interface area is entered. With Box 2, different thermal situations can be monitored in the phase field method. In this section, we introduce four different thermal models and demonstrate how to define k bulk and k interface for the respective models.

Box 2.
Thermal conductivity update algorithm for phase field modeling.
Step 1: Inputs at x int : phase field variable d, thermal conductivity k bulk and k interface .

Four Thermal Conductance Contact Models
Four different thermal conductance models are introduced in this section and implemented in the following applications. The definitions of k bulk and k interface in different thermal scenarios are displayed in Table 1. For the fully conductive model (FCM), the thermal conductivity is kept as a constant, assuming the whole model is in the undamaged region with k = k 0 . For the adiabatic model (ADM), a small enough value for the k interface avoids heat interaction through the crack region. Normal contact stress is utilized for both the pressure-dependent (PDM) and uncoupled (UCM) models, meaning that the capability of heat transfer is dominated by both the mechanical loading and the damage effect. It is noted that the thermal expansion coefficient β is set to be very small in the UCM in order to eliminate the coupling effect from the mechanical and thermal field.
Here, we use the pressure-dependent thermal contact model as an example to show how the process of thermal conductivity updating works. Instead of having a thermal contact conductance formulation at the contact surface, in the phase field formulation, for the first time, we develop a formulation for the interfacial thermal conductivity that is a function of the phase-field variable and the contact pressure, which accounts for the normal stress-related heat exchange through the spots. Inspired by the thermal conductance expression from discrete methods [52], we propose the following formulation for the pressure-dependent thermal conductivity: where t N,bulk is the contact traction in the normal direction, H e is the Vickers hardness coefficient, and ε c is the thermal constant. The updates for the thermal conductivity based on the pressure-dependent contact conditions are summarized in Box 3. The bulk and interface thermal domains are separated by the phase field variable d with the threshold setting. When d < threshold, we consider the bulk domain with A = 0, and employ the expression of k bulk for the thermal conductivity. Otherwise, we consider the interface domain with A = 1 and k interface . For the pressure-dependent model, we use the norm of the normal strain to differentiate the contact and debonding conditions. When ε N ≤ 0, we apply the pressure-dependent model to calculate k interface .
Box 3. Thermal conductivity update algorithm for phase field modeling of the pressuredependent contact model.
Step 1: Inputs at x int : the strain ε, the thermal conductivity k 0 , the phase field variable d, the normal unit vector n, and the tangential unit vector m, the threshold value: e −1 ≈ 0.378. Step 2: Set k = (1 − A)k bulk + Ak interface Step 3: If the phase field variable d < threshold, Bulk domain A = 0, THEN k = k bulk Set k bulk = g(d)k 0 EXIT Step 4: ELSE, the phase field variable d ≥ threshold, Interface domain A = 1, THEN k = k interface Calculate the normal strain ε N = ε : (n ⊗ n) If ε N > 0, THEN Non−contact condition, k interface = g(d)k 0 , EXIT ELSE, Contact condition,

Numerical Examples
In this numerical simulation section, we investigate the performance of the proposed method by means of representative numerical examples. Here, we aim to: (1) demonstrate our novel thermal formulation for the thermomechanical contact problem within the context of the phase field framework; (2)  In the following sections, we demonstrate the feasibility of the thermal formulation through two pure heat conduction problems, then compare our results with the numerical solutions from the X-FEM method in Section 5.1. Two different boundary conditions are explored, namely, the isothermal and the adiabatic situations. Then, the fully coupled thermomechanical frictional contact model is investigated, with an inner crack in Section 5.2 and a prescribed interface in Section 5.3 used to illustrate the accuracy of the proposed method. Comparisons with the literature for the mechanical phase (Section 5.2) and the coupled thermomechanical phase (Section 5.3) are presented under four different thermal contact conductance conditions. Lastly, we study a frictional crack propagation problem in Section 5.4 to show the ability of the proposed method to capture cracking phenomena using the thermomechanical coupled setting.

Square Plate with (a) Central Horizontal and (b) Inclined Cracks
In the first numerical example, we demonstrate the feasibility of the proposed thermomechanical coupled phase field formulation in a pure thermal setting. Two different inner cracks are utilized, namely, a horizontal crack and a slanting crack. As shown in Figure 3a,b, we first consider a square domain with a horizontal center crack that is subjected to two different thermal boundary conditions [53] for the adiabatic (Figure 3a) and isothermal (Figure 3b) cases. For the adiabatic case, temperatures of T = 1 • C and T = −1 • C are imposed on the top and bottom surfaces, respectively, with no heat flux transferring across the crack [49]. For the isothermal case, the temperature is maintained at T = 1 • C for all edge surfaces, while the temperature field of the inner crack is different from the boundaries, with T = −1 • C. Figure 4 depicts the results for the model with a slant center crack using the same thermal boundary settings.
For the horizontal center crack model shown in Figure 3, the block size is 2 m × 2 m, with the red line depicting the inner crack. The inner crack is located in the middle, with crack length a of 1.2 m, meaning that a/L = 0.6, where L is the side length of the domain. The mesh size of this simulation is 50 × 50, as shown in Figure 5a, with the length parameter set to l c = 0.16 m. For the adiabatic case, in order to guarantee that the heat transferring through the crack is small enough to be ignored, the interface thermal conductance k interface is set to 10 −7 W/mm • C.  For the slant center crack, as illustrated in Figure 4, the domain size is 1 m × 1 m, with a stationary crack starting from (0.3 m, 0.33 m) to (0.7 m, 0.68 m). The problem domain is discretized by unstructured meshes, as shown in Figure 5b, with an approximated mesh size of h = 0.05 m and with l c = 0.2 m. Detailed material properties for these two crack models are summarized in Table 2. The resulting contour plots of the temperature distribution for the models with horizontal and slant center cracks under two different thermal conditions (the adiabatic and isothermal cases) are exhibited in Figure 6. For the adiabatic case, shown in Figure 6a,c, the discontinuities in the thermal field are clearly shown near the fracture region, as the two sides of the crack are supposed to be insulated from each other. The width of the observed crack is controlled by the length parameter l c , which is the width of the diffusive damage band. On the other hand, for the isothermal case, shown in Figure 6b,d, the thermal contours are uniformly distributed. Compared to the analysis in [49,53], where the X-FEM method is used, we observe that the results match remarkably well with the reference. This agreement demonstrates that the proposed thermomechanical coupled phase field formulation is able to correctly identify and reproduce the temperature discontinuity and continuity profile near the crack regions for both the isothermal and adiabatic cases.

Remark 3.
For the adiabatic case, as shown in Table 2, k interface should be set with a sufficiently small value rather than 0 in order to avoid any singularity in the matrix calculation.

Squared Domain with an Internal Slant Frictional Crack
After validating the pure thermal problem, we extend the model with the slant center crack in Section 5.1 to model a thermomechanical coupled problem with a frictional contact crack. The purpose of the second example is to investigate the ability of the proposed phase field method for thermomechanical problems with frictional contact. The problem was initially proposed by Dolbow et al. [20], and was later elaborated in [48] to verify the phase field method for frictional contact cracks in purely mechanical problems. As shown in Figure 7a, the geometry of the domain is a 1 m × 1 m square with an existing crack located from (0.3 m, 0.33 m) to (0.7 m, 0.68 m) [48]. A uniform displacement loading d = −0.1 m in total is applied on the top surface in 10 steps. We adopt the same material properties used in [48], with elastic modulus E = 10, 000 MPa, Poisson's ratio ν = 0.3, and the frictional coefficient of the crack surface µ f = 0.1. The two fractional mesh distributions are shown in Figure 7, in (b) with coarse mesh size h = 0.008 m, l c = 0.032 m and in (c) with the finest mesh size h = 0.002 m, l c = 0.008 m are investigated in this section. In this simulation, the crack driving force parameter H is calculated following the same strategy as in Appendix A of Borden et al. [4], and the critical fracture energy is set to G c = 50 N/m 2 . As the main focus of this paper, we investigate the effect of the different thermal conductance conditions on the mechanical phase. In addition to mechanical loading, we apply the thermal boundary conditions shown in Figure 7a, with 10 • C prescribed to the top surface and 0 • C to the bottom surface, for further exploration.
In Table 1, we display different thermal conductance models with different thermal conductivity k as provided in Equation (22). The thermal expansion coefficients are presented in Table 1. The four different thermal conductivity models discussed in Section 4.2, namely, the pressure-dependent model (PDM), uncoupled model (UCM), fully conductive model (FCM), and adiabatic contact model (ACM), are employed. First, we compare the normal contact stress along the crack with [47] using the uncoupled model ( Figure 8). The comparison shows that with mesh refinement, the proposed method approaches the reference solution with good agreement.
The results of the displacement and temperature contours at the last step using the fine mesh are shown in Figure 9. It can be noted that for the uncoupled case (UCM), we set the thermal expansion coefficient to be small enough that the coupling effect can be ignored. Comparing our results with the purely mechanical simulation presented in Fei and Choo's work [48], excellent agreement on the displacement contours is achieved. For the other three cases with thermal effects, displacements in the x direction are more significant compared to the UCM case due to coupled thermomechanical effects. Because the displacement in the y direction is fixed, there is not much difference shown on contour plots in Figure 9.
As for the thermal field, the temperature distribution is uniform in the fully conductive model with the same thermal conductivity for the bulk (k bulk ) and the interface (k interface ) domains. Thus, in the fully conductive case shown in Figure 9 (3), the location of the inclined crack is not visible. For the other three models, the fracture crack in the temperature contours can be clearly observed. As shown in Table 1, sharing the same k interface , which is determined by the contact pressure and damage parameter, the difference in the temperature distribution for the PDM and UCM case is dominated by the coupling term between the thermal and mechanical field. As shown in Figure 9 (3), the PDM case with a stronger coupling effect exhibits a more continuous temperature distribution around the crack region. In the adiabatic model (ACM), with k interface set to be a small enough constant, there is almost no heat transference in that area, and strong discontinuity in the temperature field is visible.

The Frictional Sliding Problem under Different Thermomechanical Coupled Conditions
In this section, we use the proposed frictional contact model to solve a frictional sliding problem. The problem was simulated in our previous paper, where a discrete VMDG method was used [50]. The detailed geometry is shown in Figure 10, where the 1 m × 1 m domain is subjected to two prescribed displacement loadings on the top surface. In addition, a heat fluxq n = 1200 W/m 2 is applied on the top surface pointing inward, with 0 • C fixed on the bottom surface. The total displacement is applied in 10 steps, while the time increment ∆t is set to be ∆t = 2.8 × 10 16 s to ensure that the thermal phase is in the quasi-static state. As stated in Section 4.2, four different thermal conductivity conditions are considered, as provided in Table 1. For the adiabatic contact model (ACM), as mentioned in [52], it should be noted that the heat flux boundary condition must be replaced by the temperature boundary condition T = 10 • C on the top surface to avoid the singularity in the equations of the thermal system. The detailed material properties are provided in Table 3.  In Figure 11, the contours of temperature distribution are shown in their deformed configurations. A comparison is performed in this figure among the different thermal contact conductivity models, i.e., PDM, UCM, FCM, and ADM. In addition, we compare the results with our previous work [50], where the VMDG method was used. Very good agreement is achieved between the proposed phase field method and the previous VMDG method. It is noted that a slight difference in the temperature contour between the VMDG and the proposed phase field method can be observed in the adiabatic case, as shown in Figure 11. Because the discontinuity is controlled by the length parameter l c in the phase field method, when s smaller length parameter is applied the sharper discontinuity is captured in the temperature as well as the displacement. The VMDG method is a discrete method with duplicated nodes at the fracture surface; thus, the sharp discontinuity in the displacement is captured exactly. In contrast, in the proposed phase field method, discontinuity is captured in a diffuse fashion influenced by the mesh size and the length parameter l c . To further assess the performance and accuracy of the proposed phase field model for the frictional sliding problem, we compare the normal contact stress results with Khoei and Bahmani's [52] analysis, where the X-FEM method was used. As in [52], we use a uniform structured FE mesh of 75 × 75 quadrilateral elements, and the penalty parameters are used to constrain the tangential and normal stick conditions. As shown in Figure 12, a strong agreement is achieved when comparing the proposed phase field method to the X-FEM method. As stated in the previous sections, the proposed phase field method does not require a specific treatment or algorithm at the contact surface and is easier to implement compared to discrete methods.

Propagation of an Inclined Crack in a Rectangular Plate
To enrich the applications of the proposed method, in this example we simulate the crack propagation to demonstrate the capacity of the framework following the verification of the stationary interface problems. As shown in Figure 13a Following the algorithm proposed by Amor [54], the strain history energy H should be updated and employed in the calculation at each load step, rather than arbitrarily defined as for stationary cracks, owing to the dynamic crack propagation. This means that the evolution of the phase field damage parameter d is dependent on the displacement from previous steps. The thermal expansion coefficient β is chosen to be a relatively small number in this simulation.
The evolution of the phase field parameter d and the temperature contour plots are shown in Figure 14 at three different load steps and with three different thermal contact conductance models. The propagation of the crack and the evolution of the temperature are clearly captured using the proposed phase field method. For the fully conductive case (FCM), the temperature is uniformly distributed, as the heat is free to transfer between the upper and lower parts of the domain. For the adiabatic case (ACM), the temperature field exhibits a sharp discontinuity across the damaged crack. As the crack propagates, the temperature distribution shows a stronger discontinuity, while the damage region expands. At the last load step, when the crack is fully developed, and the whole domain is separated into two parts, and the temperature fields of the upper and lower parts are completely discontinuous. For the pressure-dependent model (PDM), when compared to the ACM case, heat transfer is more obvious at the beginning of the crack propagation process, as it is dependent on the magnitude of the contact stress at the crack region. During the crack initiation and propagation stage, the temperature discontinuity becomes less distinguishable as the crack propagates, because k interface is dominated by increasing displacement loading rather than the phase field damage variable d. After the crack is fully developed, the damage variable is dominated by the conductivity, and there is no thermal interaction between the two subdomains.

Pressure dependent
Fully conductive Adiabatic contact Damage contour

Concluding Remarks
In this paper, a thermomechanical coupled phase field method has been developed for modeling cracks with frictional contact. In the mechanical phase, the proposed method follows the work of Fei and Choo [48], regularizing the stress tensor at contact regions under different contact conditions, with a specific focus on frictional contact to accommodate the no-penetration constraints. For the thermal phase, we propose a novel formulation for thermal conductivity in a regularized contact region by accommodating different thermal conductance conditions and the phase field damage parameter. We investigate several benchmark examples to assess the proposed computational method. First, we verify the proposed thermal algorithm using different pure thermal models. Then, we simulate a problem with an internal slant frictional crack and extend the model into thermomechanical coupled fields. In the third example, we study the thermomechanical coupled frictional sliding problem and compare the results with the X-FEM and the VMDG methods, finding good agreement, which demonstrates the capacity of the proposed phase field framework in simulating coupled thermomechanical problems with frictional sliding. Lastly, we demonstrate that the proposed method can simulate crack propagation with thermal effects. Based on the above conclusions, the proposed method has potential applications in various areas, including modeling of thermomechanical coupled crack propagation of composite structures, deposition processes in friction stir welding, and metal forming processes.
There are several key features in this paper, including (i) the formulation accommodates no-penetration constraints of the contact behavior in the phase field method and (ii) the formulation accommodates different thermal conductance conditions without the need for a sophisticated algorithm to impose thermal contact constraints on crack surfaces. The proposed phase field method can be implemented easily in comparison to the existing discrete methods for frictional cracks. The proposed model can be applied in the future for complex geometry and crack branching, and merging in coupled multi-physics settings. Although the proposed thermomechanical coupled phase field method is promising, it has a number of deficiencies. The main limitations are threefold. First, the inertial effect in the mechanical field is not considered. Second, in terms of accuracy, the phase field method does not show superior performance for stationary cracks compared to discrete methods. Third, without duplicating nodes at the cracking interfaces it is difficult to capture the sharp discontinuity along the cracking surfaces. With all these considerations, our future work will focus on advancing the model in order to (1) accommodate dynamic effects in the mechanical field, (2) adopt a discrete method for predefined crack regions to capture the sharp discontinuities, and (3) track both stationary interface cracks and propagating cracks simultaneously. Work is underway to adopt the VMDG method [50] to model stationary interfaces and combine it with the present phase field model for crack propagation.
Funding: This research was supported in part by a Seed Grant award from the Institute for Computational and Data Sciences at the Pennsylvania State University.
Institutional Review Board Statement: Ethical review and approval were waived for this study, due to the reason that the studies are not involving humans or animals.

Conflicts of Interest:
The authors declare no conflict of interest.