Phase Field Models for Thermal Fracturing and Their Variational Structures

It is often observed that thermal stress enhances crack propagation in materials, and, conversely, crack propagation can contribute to temperature shifts in materials. In this study, we first consider the thermoelasticity model proposed by M. A. Biot and study its energy dissipation property. The Biot thermoelasticity model takes into account the following effects. Thermal expansion and contraction are caused by temperature changes, and, conversely, temperatures decrease in expanding areas but increase in contracting areas. In addition, we examine its thermomechanical properties through several numerical examples and observe that the stress near a singular point is enhanced by the thermoelastic effect. In the second part, we propose two crack propagation models under thermal stress by coupling a phase field model for crack propagation and the Biot thermoelasticity model and show their variational structures. In our numerical experiments, we investigate how thermal coupling affects the crack speed and shape. In particular, we observe that the lowest temperature appears near the crack tip, and the crack propagation is accelerated by the enhanced thermal stress.


Introduction
Cracking is a phenomenon that occurs everywhere in our lives, but, if it is allowed to continue, it can cause fatal damage. A crack in a material occurs when the material experiences a continuous overload. However, several other factors, such as thermal expansion and contraction due to temperature changes [1][2][3], fluid pressure (e.g., in hydraulic fracturing) [4], the diffusion of hydrogen (or hydrogen embrittlement) [5,6], chemical reactions [7], and humidity [2], cause cracks in materials. In particular, among these phenomena, cracks due to thermal expansion are interesting to study from the viewpoint of the energy balance between elastic, thermal, and surface energies. M. A. Biot proposed a theoretical framework for coupled thermoelasticity based on the principle of minimum entropy production [8]. Biot's model is now widely known as the traditional coupled thermoelasticity model, and it has been extended to dynamical theory [9] and to various other situations [10][11][12][13][14][15]. As shown in Section 2.2, it satisfies an energy balance equality between the elastic and thermal energies.
In fracture mechanics, especially in the modeling and simulation of crack propagation, a phase field approach has been recently recognized as a powerful tool. The phase field model (PFM) for fractures was first proposed by Bourdin et al. [16] and Karma et al. [17]. Then, based on the framework of variational fracture theory [18,19], the techniques and applications of PFM have been extensively developed, for example [20][21][22][23][24][25]. We refer to [26] for further information on the development of PFM for fracture mechanics. PFM for fracture mechanics is derived as a gradient flow of the total energy, which consists of The organization of this paper is as follows: in Section 2, we introduce the linear thermoelasticity model by M.A. Biot and derive its variational principle and energy dissipation property. In addition, we numerically investigate the effect of the thermal coupling term on the elastic and thermoelastic energies in an expanding region. Section 3 is devoted to PFMs for crack propagation under thermal stress. In Section 3.1, we give a brief review of the irreversible fracturing phase field model (F-PFM) and its energy equality, which guarantees the energy dissipation property (Theorem 2) and follows the works [22,26]. In Sections 3.2 and 3.3, we propose two types of thermal fracturing phase field models (TF-PFMs). The first model, TF-PFM1, is a straightforward coupling of F-PFM and the Biot thermoelasticity model. Based on the variational principle of the Biot model (Proposition 2), we show a partial energy equality for a fixed temperature (Theorem 3). However, it does not satisfy the energy equality for the total energy, which consists of the elastic, thermal, and surface energies.
The second model, TF-PFM2, presented in Section 3.3 is another natural coupling of F-PFM and the Biot thermoelasticity model based on the energy equality of the Biot model (Theorem 1). We prove an energy equality for TF-PFM2 in Theorem 4. Since we consider several models (Biot's model, F-PFM, and TF-PFMs) and their energy qualities, for the readers' convenience, we list the energies and energy equalities for each model in Tables 1 and 2.
In Section 4, we show some numerical comparisons between two TF-PFMs using non-dimensionalized equations. We investigate the effects of the thermal coupling in TF-PFM1 and TF-PFM2 on the crack speed and the crack path by changing a dimensionless coupling parameter δ. As noted, although the temperature influences material properties micro-structurally [33], it is not considered in the present study for simplicity. Generally, the effect of micro-structure of material gives the typical crack path, such as: curving, branching, kinking, etc. A clear study of this is addressed by [34]. The last section shows some conclusions and comments on further topics.
To easily understand the relevant notation and symbols in this paper, we introduce them in this section. Let Ω be a bounded domain in R d (d = 2 or 3). The position in R d is denoted by x = (x 1 . . . . , x d ) T ∈ R d , where T denotes the transposition of a vector or matrix. Let ∇, div, and ∆ be the gradient, divergence, and Laplacian operators with respect to x, respectively. For simplicity, we writeu,Θ, andż as the partial derivatives of u, Θ and z with respect to t, respectively. For simplicity, we often denote u(t) := u(·, t), etc. The space of the real-valued (symmetric) d × d matrix is denoted by R d×d (R d×d sym ). The inner product of square matrices A, B ∈ R d×d is denoted by A : B := ∑ d i,j=1 A ij B ij . Using L 2 (Ω), we refer to the Lebesgue space on Ω, while H 1 (Ω, R d ) and H 1 2 (Γ u D , R d ) represent the Sobolev space on Ω and its trace space on the boundary Γ u D , respectively. For more details on Sobolev spaces, we refer to the review in [35]. In addition, we summarize the physical properties used in this paper in Table 3. Table 1. List of energies.

Formulation of the Problem
M.A. Biot [8] proposed the following mathematical model for coupled thermoelasticity: where Ω is a bounded domain in R d (d = 2 or 3). We suppose that Ω is an isotropic elastic body and consider the thermoelastic coupling between the mechanical deformation and the thermal expansion in Ω. The constant β is defined by β := a L (dλ + 2µ) with a L > 0 as the coefficient of linear thermal expansion and µ(> 0); λ(> − 2µ d ) are Lamé's constants. The unknown functions in (1) and (2) are the displacement u(x, t) = (u 1 (x, t), . . . , u d (x, t)) T ∈ R d and the temperature Θ(x, t) ∈ R. In addition, the constant Θ 0 > 0 is a fixed reference temperature. Similarly, strain e[u] and stress tensors σ [u] are defined as is an isotropic elastic tensor and I is the identity matrix of size d. From (3b), (1) is also written in the form The term β∇Θ in (1) and the term Θ 0 β ∂ ∂t (divu) in (2) represent the body force due to thermal expansion and the heat source due to the volume change rate, respectively. We remark that, when a L = 0, (1) and (2) are decoupled.
It is convenient to introduce the following strain and stress tensors, including the thermal effect: Using the thermal stress tensor σ * [u, Θ], (1) can be written in the following form: This means that the force σ * [u, Θ] is in equilibrium in Ω. In the preceding, Equations (1) and (2) represent the force balance and the thermal diffusion in Ω, respectively.
The system in (1) and (2) is complemented by the following boundary and initial conditions: in Ω, where n is the outward unit normal vector along the boundary, The boundaries Γ u D and Γ u N (Γ Θ D and Γ Θ N ) are the Dirichlet and Neumann boundaries for u (for Θ), respectively. We suppose that the (d − 1)-dimensional volume of Γ u D is positive for the solvability of u.

Remark 1.
Instead of boundary conditions (5a) and (5b), we can also consider the following mixed-type condition. When d = 2, on a part of the boundary (which we denote by Γ u DN ), u = (u 1 , u 2 ) T and where u Di := Γ u DN → R is a given horizontal or vertical displacement and e 1 = (1, 0) T , e 2 = (0, 1) T . These types of mixed boundary conditions are considered in Sections 2.3.3 and 4.4.1. Even for these mixed-type boundary conditions, we can easily extend the following arguments on weak solutions, variational principles, and energy equalities.

Variational Principle and Energy Equality
This section aims to show a variational principle and provide an energy equality that implies the energy dissipation property for the system (1) and (2). In linear elasticity theory, a weak form of the boundary value problem for which is given by A weak solution uniquely exists and is given by is an elastic energy. This is known as a variational principle [37,38]. For a fixed Θ(x), a weak form for u of (1) and its variational principle are derived as follows: is equivalent to the following weak form: The equivalency immediately follows from this equation: there exists a unique weak solution u ∈ H 1 (Ω; R d ) that satisfies (9). Furthermore, the solution u is a unique minimizer of the variational problem: We remark that E * el (v, Θ) represents thermoelastic energy.
Proof. The unique existence of a weak solution for u is shown by the Lax-Milgram theorem [38] since (9) is written as The coercivity of the above weak form is known as Korn's second inequality [38]: For a weak solution u and any v ∈ V u (0), using the equalities we have This shows that u is a minimizer of E * el (u, Θ) among V u (u D ).
On the other hand, if u is a minimizer, the first variation of E * el vanishes at u; i.e., for all v ∈ V u (0), we have Hence, u is a weak solution. Summarizing the above, there exists a unique weak solution to (8), and u is a weak solution if and only if it is a minimizer of E * el among V u (u D ).
The next theorem represents a dissipation of the sum of the elastic and thermal energies during the thermomechanical process. We define thermal energy as Theorem 1 (Energy equality for Biot's model). Let (u(x, t), Θ(x, t)) be a sufficiently smooth solution to (1), (2) and (5). In addition, we suppose that u D does not depend on t and Θ D = Θ 0 . Then, Proof. Since d dt we obtain Substituting (2) into (14) and using the boundary conditions (5c) and (5d) for Θ, we obtain This gives the energy equality for (5).
As shown in Proposition 2 and Theorem 1, Biot's thermoelasticity model is related to both energies E el (u) and E * el (u, Θ). We denote their energy densities as follows: where W(u) and W * (u, Θ) are the elastic and thermoelastic energy densities, respectively.

Non-Dimensional Setting
In the following numerical examples, we introduce a non-dimensional form of Biot's model. We consider the following scaling for x, t, u, C (or λ, µ), and Θ: where c x , c t , c u , c e , and c Θ > 0 are the scaling parameters. Let c x [m], c e [Pa], and c Θ [K] be characteristic scales for the length of the domain, the size of the elastic tensor, and the temperature, respectively. The parameters c t and c u are defined as where (1) and (2) are written in the following non-dimensional form: The system (19) has only three parameters,λ,μ, and δ. The parameter δ is a nondimensional thermoelastic coupling parameter defined by and δ > 0. If we choose δ = 0, (19b) is decoupled from (19a), and the temperature fieldΘ in (19a) is essentially a given function. In the following example, the case δ = 0 is referred to as the uncoupled case. Under the above scaling, we denote the (thermo)elastic strain, stress tensors, and (thermo)elastic energy densities as follows: In the following section, we apply these non-dimensional forms and omit ∼for simplicity.

Numerical Setup and Time Discretization
In the following examples, we set Young's modulus E Y = 1, Poisson's ratio ν P = 0.32, the coefficient of linear thermal expansion a L = 0.475 and the thermoelasticity coupling parameter δ = 0.0, 0.1, 0.5 in the non-dimensional form of (19). We consider two numerical examples for (19), an L-shaped cantilever domain and a square domain with a crack (more precisely, a very sharp notch), as illustrated in Figure 2.  Figure 2. An L-shaped cantilever (left) and a cracked domain (right) with the subdomain A as an observation area.
We apply the following implicit time discretization for (19): where u k and Θ k are approximations to u and Θ at t = k∆t (k = 0, 1, 2, . . .). At each time step k = 1, 2, . . ., we solve (21) with given boundary and initial conditions (5) using the finite element method. The details of the weak forms for (21) and their unique solvability are described in Appendix A.
In observation area A illustrated in Figure 2, we define the average of (thermo)elastic energy densities in A as follows: and the differences between W (A) and W * (A) for each δ > 0 and for δ = 0 are defined by In the following examples, we use the software FreeFEM [39] with P2 elements and unstructured meshes. For the time interval and time step, we use 0 ≤ t ≤ 0.1 and ∆t = 1 × 10 −4 , respectively.

L-Shape Cantilever
Here, we consider the L-shaped cantilever whose left side is fixed, and the vertical displacement u 2 is given on the right side, as illustrated in Figure 2 left. We denote the left and right boundaries by Γ u D and Γ u DN , respectively, and define Γ u . The boundary conditions for u are For Θ, we suppose ∂Θ ∂n = 0 on Γ and the initial temperature Θ * = 0. Although we adopt the above slightly modified boundary conditions in this example, the previous arguments are valid with small modifications, and we omit their details.
We apply the finite element method to (21). The total number of triangular meshes = 18,215 and the number of nodes (the vertices of the triangles) = 9301.
As shown in the lower part of Figure 3, we observe that the highest temperature is in the contracting area and the lowest is in the expanding area. Furthermore, there exists a contribution δ for each δ > 0 during the loading process. Although the disparity is small, the thermoelastic coupling parameter δ contributes to the variations in W (u) and W * (u), as shown in Figure 4a    In the L-shape cantilever case for each δ > 0, we conclude that the thermal coupling parameter enhances the singularity of (thermo)elastic energy in the expanding area. The (thermo)elastic energy plays a role in the driving force in the phase field model [23], which means that the parameter δ can accelerate crack growth in the expanding area.

Cracked Domain
Here, we consider a cracked domain with vertical displacements on the top and bottom sides, and the other sides are free traction, as shown in Figure 2 right. The boundary conditions for u are where Γ u +D and Γ u −D denote the top and bottom boundaries of Ω, respectively, and . For Θ, we suppose ∂Θ ∂n = 0 on Γ Θ N = Γ and the initial temperature Θ * = 0.
We use the finite element method to solve (21). Therefore, the total number of triangular meshes and the number of nodes (the vertices of the triangles) are 11,176 and 5722, respectively.
From Figure 5 left, we conclude that the area that expands the most (i.e., divu is largest) appears near the crack tip. This can be compared with the analytical solution for the linear elasticity in a cracked domain in Appendix B. We also observe that the region with the lowest temperature appears to the right of the crack tip in Figure 5 right. From the temporal change in the temperature along the x 1 axis plotted in Figure 6 right, we also observe that the lowest temperature region appears in 0.5 < x 1 < 0.6 and that the temperature decreases over time. This is shown in Figure 6 left, where the value of divu is plotted along the x 1 axis and divu is increasing over time; i.e., the heat source term divu in (2) is positive. Experimentally, the lowest temperature around the crack tips does not match with the studies of of Zehnder et al. [40], Rusinek et al. [41], and Wang et al. [42]. They record that the highest temperature occurs around the crack tips, which result from a plastic zone around the crack tips. In the present study, we do not consider the plastic zone. However, it would be more possible to use a thermo-viscous-elasticity condition. We will consider and study it using the thermo-viscous-elasticity equation in the future work [43].  Similar to Section 2.3.3, for each δ > 0, we obtain variations of W (A) and W * (A) in subdomain A (Figure 7), where the subdomain A corresponds to the area that expands the most. From Figure 7, it is observed that W * (A) is larger than W (A). This suggests that the thermoelastic energy density W * (u, Θ) has a higher value than the elastic energy density W(u). These observations are confirmed by the comparison of our thermal fracturing phase field models.

Crack Propagation under Thermal Stress
This section is devoted to the phase field models for thermal fracturing, which are the main purpose of this paper.

Fracturing Phase Field Model (F-PFM)
According to the works [22,26], we introduce fracturing PFM (we call it F-PFM) in this section. Let Ω be a bounded (uncracked) domain in R d and Γ := ∂Ω = Γ u D ∪ Γ u N , similar to Section 2. In F-PFM, a crack in Ω at time t is described by a damage variable z(x, t) ∈ [0, 1] for x ∈ Ω with space regularization. The cracked and uncracked regions are represented by z ≈ 1 and z ≈ 0, respectively, and z ∈ (0, 1) indicates slight damage. A typical example of a straight crack in a square domain is illustrated in Figure 8.
The F-PFM is described as: with the following boundary and initial conditions: In (22b), the term W works as a driving force for z. The symbol ( ) + on the right-hand side in (22b) denoted the positive part (s) + := max(s, 0), and it represents the irreversible property of crack growth.
F-PFM is derived as a unidirectional gradient flow of the total energy E el (u, z) More precisely, u(t) obeys the following variational principle: and (22b) becomes a gradient flow of the energy min u E el (u, z) + E s (z).
We remark that E el (u, z) is a modified elastic energy, which corresponds to the elastic energy with a damaged Young's modulusẼ Y = (1 − z) 2 E Y . The energy E s (z) is regularized surface energy, which approximates the crack area (d = 3) or length (d = 2) as → 0. Please see [26] for more details. The following energy equality for F-PFM is shown in [26] ( [22] for the antiplane setting).
Theorem 2 (Energy equality for F-PFM). Let (u(x, t), z(x, t)) be a sufficiently smooth solution to (22) and (23). If u D is independent of t, then we have Proof. Differentiating the total energy in t and applying integration by parts, we obtain where we define H := div(γ * ∇z) − γ * z + (1 − z)W(u). Since (22b) is written as αż = (H) + , using the equality H(H) + = (H) 2 + , we conclude that

Thermal Fracturing Phase Field Model 1 (TF-PFM1)
To combine the Biot model in (1) and (2) and F-PFM in (22), their variational principles for u, Proposition 2, and (26) suggest that we consider the following modified thermoelastic energy: and a variational principle: From the definition of the modified thermoelastic energy (29), it is natural to replace the driving force term W(u) = σ[u] : e[u] in (22b) by the thermoelastic energy density For heat Equation (2), since β = a L (dλ + 2µ) and Lamè's constants (λ, µ) are replaced by damaged constants ((1 − z) 2 λ, (1 − z) 2 µ), β should also be replaced by damaged constant (1 − z) 2 β. The thermal conductivity κ 0 is also considered to be modified by z because the heat is usually insulated across the crack. We suppose κ = κ(z) > 0 in this section, and we set it as κ(z) = (1 − z) 2 κ 0 in Section 4.
Summarizing the above statements, we obtain the following thermal fracturing model, PFM 1 (TF-PFM1): Similar to (1), (2) and (22), the boundary and the initial conditions to solve (31) are presented as follows: in Ω.
In the following, for simplicity, we define As a natural extension of Proposition 2 and Theorem 1, we obtain the following "partial" energy equality for TF-PFM1.

Thermal Fracturing Phase Field Model 2 (TF-PFM2)
In the previous section, we proposed TF-PFM1 based on the thermoelastic energy E * el (u, Θ). We proved a variational principle but proved only partial energy equality. As shown in Section 2.2, the Biot model is related to both energies E * el (u, Θ) and E el (u). The variational principle holds for E * el (u, Θ) (Proposition 2), and the energy equality holds for E el (u) (Theorem 1). This motivates us to consider another type of thermal fracturing PFM based on elastic energy E el (u). We call the following thermal fracturing model TF-PFM2: The associated boundary and initial conditions are given by (32). For this model, we can show the following energy equality.
Theorem 4 (Energy equality for TF-PFM2). We suppose that (u(x, t), Θ(x, t), z(x, t)) is a sufficiently smooth solution for (35) and (32). If u D is independent of t and Θ D = Θ 0 , then the following energy equality holds: Proof. Since the relation in (13) is written as d dt Hence, we have On the other hand, Taking a sum of these equalities (37) and (38), then we obtain the energy equality (36).

Numerical Experiments
In this section, we conduct numerical experiments to test F-PFM, TF-PFM1, and TF-PFM2, which were derived in Section 3, and report the numerical results. Through the numerical experiments, we observe the effect of thermal coupling on the crack speed and the crack path during its growth process.

Time Discretization
To solve problem (39), we adopt the following semi-implicit time discretization scheme [22,26]: where u k , z k , and Θ k are the approximations of u, z, Θ, respectively, at time t k := k∆t (k = 1, 2, 3, . . .). Since the adaptive mesh technique in the FEM is often effective and accurate in numerical experiments with phase field models, problems (41) and (42) are calculated using adaptive finite elements with P2 elements with a minimum mesh size of h min = 2 × 10 −3 and a maximum mesh size of h max = 0.1. The adaptive mesh control at each time step is performed by the adaptmesh() command in FreeFEM based on the variable z. An example of the adaptive mesh is illustrated in Figure 9 right. In addition, the code for the following numerical experiments in the current study is written on FreeFEM [39] and executed on a desktop with an Intel(R) Core i7−7820X CPU@3.60 GHz, 16 core processor, and 64 GB RAM.
x 1 x 3 Figure 9. Domain for Section 4.3 with z * (x) as the initial crack (left) and the adaptive mesh for the initial crack (right).

Thermoelastic Effect on the Crack Speed
We set a square domain Ω := (−1, 1) 2 ⊂ R 2 with the initial crack z * (x) := e (−(x 2 /η) 2 ) /(1 + e (x 1 /η) ) and η = 1.5 × 10 −2 . The initial mesh is adapted to z * (x), as illustrated in Figure 9 right. The material constants for the following examples in the non-dimensional form are listed in Table 4. The boundary conditions for u and Θ are illustrated in Figure 9 left. For z, we set ∂z ∂n = 0 on Γ.
In Figure 10, the numerical results obtained by F-PFM, TF-PFM1, and TF-PFM2 are shown in the upper, middle, and bottom parts, respectively, where we set δ = 0.5 for TF-PFM1 and TF-PFM2. In addition, the profile of z on line x 2 = 0 is shown in Figure 11. From Figures 10 and 11, we observe that the crack propagation rate obtained by F-PFM is slower than that obtained by the others, and that the crack propagation rate obtained by TF-PFM1 is slightly faster than that obtained by TF-PFM2.
The temperature distributions obtained by TF-PFM1 and TF-PFM2 are shown in Figure 12. In the equation for Θ, the heat resource is given by −(1 − z) 2 δ d dt (divu). During crack propagation (0.4 ≤ t ≤ 0.8), the areas near the crack tip, the upper-right corner, and lower-right corner are continuously expanding when divu > 0 and ∂ ∂t (divu) > 0. Therefore, due to the negative source − ∂ ∂t (divu), lower temperatures are observed in those areas. On the other hand, at t = 1, due to the sudden compression caused by the total fracture, positive heat is generated, and a higher temperature is observed, especially near the upper-right and lower-right corners. In this condition, it does not allow for temperature discontinuities along the crack even if we set κ(z) = (1 − z) 2 κ 0 .  To see how the thermoelastic coupling parameters contribute to enhanced crack propagation, we consider δ = 0, 0.1, 0.2, 0.5 for TF-PFM1 and TF-PFM2, and their elastic and surface energies are plotted in Figure 13. From Figure 13, we observe that faster crack propagation occurs with a larger coupling parameter. The figure also shows that crack propagation using TF-PFM1 is faster than that using TF-PFM2.

Thermoelastic Effect on the Crack Path
In this section, we investigate the effect of the thermoelastic coupling parameter on crack path selection using our proposed models. Under a given temperature gradient, we consider crack propagation of an opening mode (Mode I) and a mixed mode (Mode I + II). In the following numerical examples, we also use the parameters in Table 4.

Mode I
We use an edge-cracked square domain, which is shown in Figure 14 left. We set the domain as follows: , and we define The boundary conditions for u and Θ are given as follows: The initial condition for Θ is given as Θ * = 0. For z, similar to the previous example (Section 4.3) , we set ∂z ∂n = 0 on Γ and choose the initial value as z * (x) := e (−(x 2 /η) 2 ) /(1 + e ((x 1 +0.2)/η) ) with η = 1.5 × 10 −2 . In this numerical experiment, we apply the thermoelastic coupling parameter δ = 0.5. Figure 15 shows the different crack paths obtained by the three models when Θ D = 10. Straight cracks occur in the F-PFM path since the thermal effect is ignored there. On the other hand, crack curves occur in the TF-PFM1 and TF-PFM2 paths. Here, the crack path is more curved in the TF-PFM2 path than in the TF-PFM1 path. These results show good qualitative agreement with the results reported in [44]. Figure 16 shows the crack paths for different temperature gradients Θ D = 0, 3, 5, 7, 10 obtained by TF-PFM1 (left) and TF-PFM2 (right). A larger temperature gradient generates a more curved crack path, and TF-PFM2 obtains a more curved crack path than TF-PFM1. Both have significant differences in the magnitude of angle deviation but have the same crack path directions. Therefore, it is clear that thermal expansion changes the crack path.
The temperature distributions during crack growth are shown in Figure 17. There exists a temperature discontinuity along the crack path, which is caused by κ(z) = (1 − z) 2 κ 0 . It approximately represents a thermal insulation condition across the crack. Different from the previous condition in Section 4.3, although we involve ∂ ∂t (divu), its contribution is small.

Mode I + II
According to the numerical experiment in [26], we consider the following setting for mixed mode crack propagation under a thermal gradient. Let Ω := (−1, 1) 2 ∈ R 2 , as shown in Figure 14 right, and Γ := ∂Ω. We set The boundary conditions for u are given as follows: The boundary conditions for Θ and z are the same as those in Section 4.4.1. The initial crack profile is given as z * (x) := e (−(x 2 /η) 2 ) /(1 + e ((x 1 −0.5)/η) ) − e (−(x 2 /η) 2 ) /(1 + e ((x 1 +0.5)/η) ) with η = 1.5 × 10 −2 . We fix the thermoelastic coupling parameter δ = 0.15 and change the temperature gradient to Θ D = 0, 2, 3, 5, 6.   Figure 18 shows the crack paths obtained by TF-PFM1 and TF-PFM2. The cracks are kinked, and the kink angle becomes larger when the thermal gradient Θ D increases. The two models provide similar results, but the kink angle in the TF-PFM2 crack is larger than that in the TF-PFM1 crack, as shown in Figure 19. Therefore, we conclude that thermal expansion changes the crack path.
Here, we do not show the temperature distribution during thermal expansion. We observe that the temperature distribution is quite similar to that of Mode I in Section 4.4.1, and a temperature discontinuity exists along the crack path during temperature injection. As mentioned, since it is relatively difficult to find the available experimental result for the thermal fracturing under mode I + II, we do not compare our result with the experimental result.
At the end of this section, we give a remark on the extendability of our TF-PFM to anisotropic material. When the material has a strong anisotropy, we have to take into account the anisotropies on the elasticity tensor C, a coefficient of linear thermal expansion a L , and the critical energy release rate γ * , especially among many material properties. For C and a L , we can easily include an anisotropic effect by using an anisotropic tensor in (3b) and replace the matrix a L I in (4a) by anisotropic one. On the other hand, for γ * , it has not well succeeded to include the anisotropy, which means that the dependency of the crack surface direction, even in the standard PFM, as far as the authors' knowledge.

Conclusions and Future Works
We proposed two thermal fracturing phase field models, TF-PFM1 and TF-PFM2, by coupling the Biot thermoelasticity model [8] and the fracturing phase field model (F-PFM) by Takaishi-Kimura [22,26].
For the Biot model, we studied a variational principle (Proposition 2) and energy equality (Theorem 1), which were related to different energies E * el (u, Θ) and E el (u) + E th (Θ), respectively (see Tables 1 and 2).
On the other hand, F-PFM has a gradient flow structure with respect to the total energy E el (u, z) + E s (z) and admits energy equality (Theorem 2).
As the first model, TF-PFM1 was derived based on the variational principle of the Biot model and the gradient flow structure of F-PFM, while TF-PFM2 is based on the energy equalities of the Biot model and F-PFM. The difference between them is the driving force term for the crack: W * (u, Θ) in TF-PFM1 (31b) and W(u) in TF-PFM2 (35b).
Consequently, we established partial energy equality for TF-PFM1 (Theorem 3) and energy equality for TF-PFM2 (Theorem 4). From the viewpoint of energy consistency, both models are satisfactory, but TF-PFM2 is more energetically consistent than TF-PFM1.
Based on the obtained numerical experiments, the following conclusions can be drawn.

2.
TF-PFM1 accelerates the crack speed more than TF-PFM2 ( Figure 11). On the other hand, the effect of the temperature gradient on the crack path in TF-PFM2 is larger than that in TF-PFM1 (Figures 16-19).
The analytical and numerical comparisons between the two models are briefly summarized in Table 5.   In this study, we did not consider the unilateral contact condition along the crack for the sake of simplicity. To further improve TF-PFM, the ideal unilateral condition for fracturing PFM [21,26] should be introduced in our PFM. In addition, there are many other effects that should be included in the model. For example, although we assumed that the critical energy release rate γ * (x) is a priory given, it may depend on the temperature in the real material. A possible extension of our model is to suppose that γ * depends on Θ linearly as γ * (x, Θ) =γ(1 − α 0 (Θ − Θ 0 )) for someγ > 0 and α > 0 [45].
Such relatively easy extendability is one of the advantages of PFM. However, we should remark that the energy equalities which we derived in this paper may not be valid for all extended models.  We define a weak form for (A4) as where V Θ := {ψ ∈ H 1 (Ω); ψ| Γ Θ D = 0}.