A New Definition of Peridynamic Damage for Thermo-Mechanical Fracture in Brittle Materials

TL;DR

This paper introduces a new way to model thermal fractures in brittle materials by improving how damage is calculated in peridynamics.

Contribution

A novel peridynamic damage definition is proposed, incorporating both bond count and spatial distribution for more accurate thermal fracture simulations.

Findings

The new damage model captures directional damage, improving thermal fracture simulation realism.

Numerical examples confirm the model's accuracy against analytical solutions.

Thermal shock simulations demonstrate the model's effectiveness in analyzing fracture mechanisms.

Abstract

A thermo-mechanical fracture modeling is proposed to address thermal failure issues, where the temperature field is calculated by a heat conduction model based on classical continuum mechanics (CCM), while the deformation field with discontinuities is calculated using the peridynamic (PD) model. The model is calculated using a CCM/PD alternating solution based on finite element discretization, which ensures the calculation accuracy and facilitates engineering applications. The original PD model defines damage solely based on the number of broken bonds in the vicinity of the material point, neglecting the distribution of these bonds. To address this limitation, a new definition of the PD damage accounting for both the number of broken bonds and their specific distribution is proposed. As a result, damage in various directions can be captured, enabling more realistic thermal fracture simulations based on a unified mesh discretization. The effectiveness of the proposed model is validated by comparing numerical examples with analytical solutions. Moreover, simulation results, including a thermal shock case with a transient temperature field, demonstrate the model’s ability to aid in understanding the initiation and propagation mechanisms of complex thermal fractures.

Figures

Full text

1. Introduction

With the rapid development of industry, more and more high-temperature concrete and metal materials are used. However, unpredictable thermal deformation and stress, often resulting from uneven temperature distributions or inconsistent thermal expansion coefficients, can ultimately lead to structural failure. In order to ensure the safety and reliability of the structures, it is necessary to analyze the thermal deformation and thermal stress that may occur in the structures. However, experiments are complex, and it is costly to reproduce the high temperature and high pressure environment. Therefore, numerical simulation is an alternative approach to help understand the mechanisms of the thermo-mechanical coupling response of a structure. Therefore, an efficient and accurate numerical simulation method is necessary.

CCM is widely used for continuous thermo-mechanical problems, but its partial differential governing equation cannot handle discontinuous issues such as thermal fracture. To address this, several numerical methods have been developed, including the extended finite element method (XFEM) [1], phase-field fracture method (PFM) [2], and discrete element method (DEM) [3]. XFEM is effective for simulating discontinuous problems such as interfaces and crack propagation. Jaskowiec et al. used XFEM for three-dimensional numerical thermo-mechanical modeling of a laminated structure with a very thin inner layer [4]. Kumar et al. employed XFEM for a thermo-mechanical fracture analysis of porous functionally graded cracked plates [5]. PFM is used for simulating structural damage. Badnava et al. used PFM to simulate brittle fracture and thermal cracks in two-dimensional (2D) and three-dimensional (3D) continua [6]. Zhou et al. presented a novel coupled thermo-mechanical PFM for concrete at high temperatures [7]. DEM can reproduce macroscopic behavior comparable to laboratory tests and monitor microscopic variations in the failure process. Sun et al. presented a low-temperature thermo-mechanical coupling modeling framework to simulate frost crack evolution in rock masses using the finite-discrete element method (FDEM) [8]. Although the above methods can solve the thermo-mechanical coupling problem, it is still a challenge to deal with the initiation and propagation of complex multi-cracks.

Silling introduced a novel non-local continuum model known as the peridynamic (PD) model [9]. The PD model has advantages for complex multi-cracks with unknown a priori locations. Unlike the partial differential governing equations used in CCM, the PD model employs integro-differential governing equations, making it suited for simulation of structural fractures. PD naturally simulates crack initiation and propagation without requiring any predefined crack growth criteria. Subsequently, a PMB constitutive model of PD was introduced by Silling [10]. This particular PD model is termed the bond-based peridynamic (BB-PD) model, wherein the Poisson’s ratio remains constrained to a fixed value. Recognizing this limitation, Silling introduced a mathematical construct known as ‘state’ in 2007, thereby presenting two variations: the ordinary state-based peridynamic (OSB-PD) model and the non-ordinary state-based peridynamic (NOSB-PD) model [11]. Notably, the BB-PD model can be viewed as a specialized instance within the broader framework of the state-based peridynamic model [12]. Additionally, PD model demonstrates applicability in addressing thermal conduction and thermal fracture problems. Bobaru and Duangpanya introduced a PD formulation for transient heat conduction in solids with discontinuities [13]. Oterkus et al. derived OSB-PD heat conduction equations [14]. Based on the above works, PD can be used to solve thermo-mechanical problems. Oterkus et al. presented a fully coupled PD thermo-mechanical framework [15]. PD’s inherent ability to simulate crack initiation and propagation has made it practical for engineering applications. Wang et al. developed a thermo-mechanical BB-PD model to simulate thermal cracking processes in concrete exposed to fire scenarios [16]. Zang et al. introduced a fully coupled thermo-mechanical PD model for rock fracturing under blast loading, accounting for initial pore damage [17]. Cheng et al. presented a thermo-mechanical PD model to investigate damage in engineered cementitious composite-concrete bonding specimens at high temperatures [18]. While these studies effectively simulated thermal fracture, selecting the appropriate micro-conductivity remains a challenge. Sun et al. presented a novel computational framework for analyzing thermal fracture in brittle solids by coupling PD and CCM [19]. However, the original PD model defines damage solely based on the number of broken bonds but neglects their spatial distribution, thereby introducing inaccuracies in damage quantification. To address this limitation, we introduce a novel PD damage formulation that considers both the quantity and spatial distribution of broken bonds. As a result, damage in various directions can be captured, enabling the simulation of thermal fracture based on a unified mesh discretization framework.

The structure of the remainder of this paper is organized as follows: Section 2 reviews the fundamental formulations of the BB-PD model and presents a new definition of PD damage within this framework. Section 3 revisits the thermo-mechanical PD model and establishes a framework for integrating the new definition of the PD damage into the model. Section 4 details the finite element spatial discretization, as well as the time discretization, of the proposed framework. Section 5 demonstrates the effectiveness of the proposed model through four examples. Section 6 concludes with remarks summarizing the findings and contributions of this paper.

2. A New Definition of PD Damage

In this section, we first review the BB-PD model and the original definition of PD damage. Then, we elucidate the necessity and specific formula of the new definition of PD damage. It should be noted that the new definition of PD damage presented in this paper is applicable to the NOSB-PD model and the OSB-PD model, with the BB-PD model being introduced in detail as a special case.

2.1. A Review of the BB-PD Model and PD Damage

The PD model assumes that each material point $[eqn]$ has its own neighborhood $[eqn]$ and interacts with points $[eqn]$ located in $[eqn]$ . The equilibrium equation can be expressed as follows [20]:

$$[eqn]$$

where $[eqn]$ denotes the neighborhood of $[eqn]$ with a horizon of $[eqn]$ ; $[eqn]$ signifies the external force acting on point $[eqn]$ , and $[eqn]$ represents a pairwise force function. A potential constitutive model relating force to relative displacement for linear elasticity and small deformations can be expressed as follows:

$$[eqn]$$

where $[eqn]$ and $[eqn]$ represent the projections of the displacement at point $[eqn]$ and $[eqn]$ onto the bond, respectively; $[eqn]$ is the unit vector of bond $[eqn]$ ; $[eqn]$ ; and $[eqn]$ denote micro-modulus functions of bond $[eqn]$ . In this paper, we focus on homogeneous materials (i.e., $[eqn]$ ). The elastic energy density can be expressed as follows [21]:

$$[eqn]$$

To capture crack initiation and propagation, the PD model requires a bond failure criterion to trigger bond failure. Bond failure can be tracked using a history-dependent function $[eqn]$ , defined as follows:

$$[eqn]$$

where s is the bond stretch; t is the computational step; and $[eqn]$ denotes the critical bond stretch. The relationship between the bond stretch s and displacement can be formulated as follows:

$$[eqn]$$

It’s worth noting that the PD model employs a scalar $[eqn]$ to describe damage at a point $[eqn]$ :

$$[eqn]$$

2.2. A New Definition of PD Damage Based on the Spatial Distribution of the Broken Bonds

As demonstrated in Equation (6), the original PD damage is calculated based on the proportion of broken bonds to total bonds. While the original PD damage can indeed characterize the degradation of materials, it is difficult to accurately portray the specific bond failure distribution within the neighborhood. To illustrate, consider the scenario depicted in Figure 1, where two differently oriented cracks traverse the neighborhood in a two-dimensional plane. The spatial distribution of the broken bonds differs between the vertical crack and the horizontal crack. However, when Equation (6) is employed to compute the damage values for material points $[eqn]$ and $[eqn]$ , these two distinct scenarios may have identical damage values ( $[eqn]$ ). It underscores the limitation of the original PD damage in capturing the specific spatial distribution of the broken bonds. On the other hand, the original PD damage establishes a homogenized mapping mechanism between microscopic bond failure and macroscopic damage fields through statistical averaging approaches. However, when addressing multi-physics coupling simulations (e.g., thermo-mechanical or hygro-mechanical interactions), the macroscopic damage field necessitates explicit consideration of anisotropic damage characteristics. This requirement exposes an intrinsic limitation of original PD damage—its inherent isotropy assumption fundamentally restricts the characterization of direction-dependent damage evolution patterns.

To address the aforementioned limitation in the PD damage, it is necessary to implement a new definition of PD damage based on bond distributions. In the new definition, the PD damage can be represented as $[eqn]$ rather than a scalar to accurately describe damage in various directions:

$$[eqn]$$

where $[eqn]$ , $[eqn]$ , and $[eqn]$ are components defined in the material coordinate system. For isotropic materials, $[eqn]$ , $[eqn]$ , and $[eqn]$ represent the damage in the directions of the x, y, and z axes, respectively. And, for anisotropic materials, $[eqn]$ , $[eqn]$ , and $[eqn]$ represent damage in three main directions. However, due to the symmetry of the PD domain, despite dividing the bond among different directions, for the same PD neighborhood, $[eqn]$ , $[eqn]$ , and $[eqn]$ are almost equal. Therefore, it remains challenging to distinguish damage among different directions.

To distinguish $[eqn]$ , $[eqn]$ , and $[eqn]$ , some improvements have been made to present a new definition of the PD damage. As shown in Figure 2, unit basis vectors $[eqn]$ along material principal directions are defined in both 2D and 3D cases. For 2D cases, $[eqn]$ , and for 3D cases, $[eqn]$ . To distinguish the positions of each bond $[eqn]$ in the neighborhood $[eqn]$ , a scalar function $[eqn]$ is defined as follows:

$$[eqn]$$

Based on Equations (6) and (8), the following formula can be obtained:

$$[eqn]$$

where $[eqn]$ ( $[eqn]$ ) denotes damage in each direction i, ∗. For 2D cases, $[eqn]$ , and for 3D cases, $[eqn]$ .

Having computed the directional damage components $[eqn]$ for each hemisphere ( $[eqn]$ ), the next step is to assemble them into a damage tensor $[eqn]$ . A critical choice is the aggregation rule for the two hemispheres along each principal direction i.

Physical and Mathematical Justification for the Maximum Rule

The tensor components are constructed using the maximum of the two hemispherical damage values:

$$[eqn]$$

This choice is motivated by both physical reasoning and mathematical pragmatism.

Physical Motivation: In brittle fracture, the macroscopic material degradation (e.g., loss of stiffness or thermal conductivity) along a given direction is dominantly governed by the most severe local damage within that directional neighborhood. A crack, which is a localized plane of broken bonds, will severely degrade properties in its normal direction regardless of the state of the opposite hemisphere. The maximum rule adopts a conservative engineering perspective by ensuring that the damage tensor reflects the envelope of directional damage severity, which is crucial for accurately driving coupled processes like anisotropic thermal conductivity reduction (see Equation (21)).

Mathematical Motivation: The goal is to map a discrete set of bond failures to a continuous, symmetric second-order tensor suitable for constitutive modeling. Alternative rules, such as taking the average ( $[eqn]$ ), can smooth out localized damage. For instance, if severe damage exists in one hemisphere ( $[eqn]$ ) while the opposite is nearly intact ( $[eqn]$ ), the average (≈0.5) would significantly underestimate the true degradation along that direction. The maximum rule preserves the monotonicity between microscopic bond failure and the macroscopic damage measure and naturally yields a symmetric, positive semi-definite damage tensor.

Behavior under Symmetric and Asymmetric Damage:

• Symmetric Damage: If damage is diffuse and approximately equal in both hemispheres ( $[eqn]$ ), then $[eqn]$ . The rule effectively reduces to a representative average for that direction.
• Asymmetric Damage: This is the typical case for a localized crack. The rule selects the hemisphere with the more severe damage as the representative value for direction i. Information about the less damaged hemisphere is not entirely lost, as it may influence the damage components in other directions. The full tensor $[eqn]$ , through the differences among its diagonal components, still captures the overall anisotropic damage pattern.

Therefore, the maximum rule provides a robust, conservative, and physically interpretable method for constructing a damage tensor from directional bond failure statistics.

For 2D cases, the new definition of the PD damage can be written as follows:

$$[eqn]$$

This means that the maximum vector should be selected among all possible damage vectors. Similarly, for 3D cases, the new definition of the PD damage can be written as follows:

$$[eqn]$$

As shown in Figure 1, for point A, the new definition of the PD damage $[eqn]$ , whereas for point B, $[eqn]$ . This means that the new definition of the PD damage varies depending on different bond failure distributions within the PD domain. In order to further clarify the effectiveness of the new definition of the PD damage in characterizing cracks, as depicted in Figure 3a, consider a crack surface that intersects the horizontal plane. Assume that all bonds crossing the crack surface are broken. The angle between this crack surface and the principal axes of the material is denoted as $[eqn]$ . Broken bonds are indicated by green dashed lines, while intact bonds are represented by blue solid lines. Furthermore, as illustrated in Figure 3b, the red line represents $[eqn]$ , the blue line represents $[eqn]$ , and the black line represents the original PD damage. As the angle $[eqn]$ changes, the new definition of the PD damage also varies accordingly. When the angle is $[eqn]$ , $[eqn]$ is the maximum and $[eqn]$ is the minimum. When the angle is $[eqn]$ , $[eqn]$ is the maximum and $[eqn]$ is the minimum.

2.3. Generalization to State-Based Peridynamic Models

The novel damage tensor $[eqn]$ proposed in Equations (11) and (12) is formulated within the BB-PD framework for clarity. However, its core concept—quantifying damage by counting broken bonds in different spatial directions—is general and can be directly extended to state-based peridynamics, including both OSB-PD and NOSB-PD models.

In state-based peridynamics, the pairwise force function is replaced by a more general force state, but the geometric concept of a “bond” as a connection between two material points persists. Therefore, the bond-failure indicator $[eqn]$ can still be defined according to a suitable failure criterion appropriate for the state-based model. Common choices include:

• A critical stretch criterion analogous to Equation (4) for certain material models.
• An energy-based criterion where a bond breaks when its contribution to the strain energy density reaches a critical fracture energy $[eqn]$ .
• A stress- or strain-invariant-based criterion, especially in the NOSB-PD model, where bonds associated with a point are considered broken when a local stress or strain measure exceeds the material strength.

Once $[eqn]$ is defined, the subsequent formulas for calculating the directional damage components $[eqn]$ (Equation (9)) and assembling the damage tensor $[eqn]$ (Equations (11) and (12)) remain identical to those presented in the bond-based formulation. The key step is the directional counting of broken bonds via the half-space indicator $[eqn]$ (Equation (8)), which is purely geometric and independent of the constitutive law.

Thus, the proposed damage tensor provides a unified measure for anisotropic damage evolution that can be coupled with multi-physics processes (e.g., anisotropic thermal conductivity reduction as in Equation (21) within both BB-PD and state-based peridynamic frameworks.

The subsequent sections introduce the applications and benefits of the new definition of PD damage in modeling thermal fracture.

3. A New Definition of PD Damage for Modeling Thermal Fracture

3.1. An Improved Thermo-Mechanical PD Model

The PD equilibrium equation with temperature can be expressed as follows:

$$[eqn]$$

The bond force can be divided into two parts:

$$[eqn]$$

where $[eqn]$ and $[eqn]$ are the bond forces of point $[eqn]$ over point $[eqn]$ and $[eqn]$ over point $[eqn]$ ; $[eqn]$ is the temperature variation of the bond. A possible constitutive equation can be written as follows [22]:

$$[eqn]$$

where $[eqn]$ denotes micro-expansivity at point $[eqn]$ ; for homogeneous materials, $[eqn]$ . Substitute Equation (15) into Equation (14):

$$[eqn]$$

where $[eqn]$ is the thermal modulus, written as follows:

$$[eqn]$$

In this paper, we focus on homogeneous materials. (i.e., $[eqn]$ ).

Consider a uniform and isotropic solid containing a heat source and undergoing heat exchange with its surrounding medium. We study the distribution and variations of temperature within the solid. This analysis is grounded in the principles of the energy conservation equation:

$$[eqn]$$

where c is the specific heat capacity; $[eqn]$ is the density; $[eqn]$ is the temperature; $[eqn]$ is the heat flux; t is time, $[eqn]$ is the heat source. The above equation can be simplified as follows:

$$[eqn]$$

According to the Fourier law,

$$[eqn]$$

where $[eqn]$ is the thermal conductivity.

3.2. Anisotropic Thermal Conductivity for Modeling Thermal Fracture

According to Section 3.1, in the thermo-mechanical model, the temperature field and deformation field can interact with each other. Changes in temperature, whether increasing or decreasing, affect the material’s thermal deformation. The bond failure can be calculated through deformation fields. Progressive accumulation of microscopic bond failures induces macroscopic damage, which can locally characterize the degradation of thermal conductivity. As shown in Figure 4, a fixed temperature $[eqn]$ is applied to the left boundary of the solid. When considering heat flow through both horizontal and vertical cracks, it is observed that the vertical crack significantly hinders heat flow, while the horizontal crack has a weaker impact. This demonstrates that cracks oriented in different directions can hinder heat flow to varying degrees. Therefore, the reduction in thermal conductivity cannot solely be attributed to the original PD damage, as described by Equation (6), and must also consider the direction and extent of cracking.

A suitable way to define thermal conductivity $[eqn]$ based on PD damage is as follows:

$$[eqn]$$

where $[eqn]$ is the reference thermal conductivity; $[eqn]$ is the unit matrix; $[eqn]$ is defined in Equations (11) and (12). This formula reflects the anisotropic effect of damage on thermal conductivity, with different degrees of degradation in different directions.

4. Numerical Algorithm

4.1. A Unified Finite Element Discretization for Thermo-Mechanical Crack Propagation

In a previous study [19], temperature fields were computed using a local model based on the finite element discretization, while deformation fields with discontinuities were computed using a PD model based on an element-free discretization. In this paper, we unify the discretization of both model through a shared mesh system. Figure 5 illustrates the proposed computational framework, where finite element discretization of the computational domain is implemented during the initialization phase. The computational procedure initiates with the computation of the temperature field. Sequentially, the temperature field is incorporated into the PD model to obtain the deformation field. The deformation field is used to determine bond failures. Progressive accumulation of bond failures induces micro-crack nucleation. Micro-cracks coalesce through damaged bands to form macro-cracks. The PD damage defined in this paper enables characterization of degradation in thermal conductivity, affecting the distribution of the temperature field.

4.2. Time and Spatial Discretization

The time discretization of the temperature field is as follows:

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ is the time increment; $[eqn]$ is the integration parameter; different values of $[eqn]$ correspond to different differential formulas. $[eqn]$ and $[eqn]$ are temperatures at time step $[eqn]$ and n.

The governing equation for temperature field computation can be expressed in the following canonical finite element form:

$$[eqn]$$

where $[eqn]$ is the heat capacity matrix; $[eqn]$ is the heat conduction matrix; $[eqn]$ is the temperature vector; $[eqn]$ is the temperature load vector; $[eqn]$ is the derivative vector of node temperature with respect to time.

The heat conduction matrix $[eqn]$ , the heat capacity matrix $[eqn]$ , and the temperature load vector $[eqn]$ can be written as follows:

$$[eqn]$$

where n is the number of total finite elements; h is the heat convection coefficient corresponding to the boundary $[eqn]$ ; $[eqn]$ is the temperature on boundary $[eqn]$ ; $[eqn]$ is the heat flux density corresponding to the boundary $[eqn]$ ; $[eqn]$ denotes the matrix of shape function; $[eqn]$ denotes the matrix of differential operators; $[eqn]$ is the heat source.

The finite element spatial discretization of PD model can be written as the following formula [23]:

$$[eqn]$$

where $[eqn]$ is the total stiffness matrix; $[eqn]$ is the displacement vector; $[eqn]$ is the external load force vector.

The total stiffness matrix $[eqn]$ and the external load force vector $[eqn]$ can be written as follows:

$$[eqn]$$

where $[eqn]$ is the amount of relative elements of point $[eqn]$ ; $[eqn]$ is the body force.

4.3. Flowchart of the Proposed Numerical Algorithm

As depicted in the flowchart in Figure 6, to address thermal fracture problems based on a shared mesh system between temperature and deformation computation, the temperature field and deformation field are computed independently. At the beginning of the algorithm, the geometric model only needs to be discretized once. Subsequently, for each time step, the new definition of the PD damage is computed to characterize the degradation of the thermal conductivity. We solve the linear Equation (24) to obtain the temperature field, which is subsequently utilized to compute the equilibrium Equation (26), incorporating temperature terms and obtain displacement field. The displacement field allows us to identify bond failures. Progressive accumulation of bond failures induces micro-crack nucleation. For quasi-static problems, if no new bond failure is detected, we proceed to the calculation for the subsequent time step, continuing this process until the end of the computation.

5. Numerical Examples

5.1. Thermal Deformation Without Damage

To validate the efficacy of the proposed model, we consider a scenario involving thermal deformation without damage. The condition of plane strain is adopted for this example. Timoshenko et al. [24] analyzed a square plate with three edges thermally insulated and mechanically restrained against normal displacement. As illustrated in Figure 7, the top edge of the plate is subject to a Dirichlet boundary condition with a temperature $[eqn]$ °C, while the initial temperature $[eqn]$ of the whole plate is 0 °C. The material properties are detailed in Table 1.

Timoshenko et al. [24] and Carslaw and Jaeger [25] derived the analytical solution for this problem as follows:

$$[eqn]$$

and

$$[eqn]$$

In this numerical example, the PD micromodulus coefficient is assumed to be an exponential function: $[eqn]$ , where $[eqn]$ is a constant coefficient that is calculated according to the given Poisson’s ratio and Young’s modulus, and l is a characteristic length. In this example, $[eqn]$ . The model is discretized into 10,000 uniform quadrilateral finite elements with dimensions of $[eqn]$ mm. Figure 8 shows the vertical displacement of points a and b with different horizons, i.e., $[eqn]$ , where $[eqn]$ is the element size. By comparing with the analytical solution, the thermo-mechanical PD model is validated. In order to improve computational efficiency and obtain accurate calculation results, $[eqn]$ is a suitable choice of horizon.

5.2. Heat Flow in a Plate with Two Thermal Insulation Cracks

To validate the proposed new definition of the PD damage, heat flow in a plate containing thermal insulation cracks is considered. The geometry configuration and boundary conditions are shown in Figure 9. There are two pre-existing thermal insulation cracks, one oriented horizontally and the other vertically. The top edge of the plate is subjected to a Dirichlet boundary condition with a temperature of $[eqn]$ °C. The initial temperature $[eqn]$ of the whole plate is 0 °C. The material properties are shown in Table 1. The model is discretized into 10,000 uniform quadrilateral finite elements with dimensions of $[eqn]$ mm. The horizon is $[eqn]$ .

For comparison purposes, a classical approach to model thermal insulation cracks adopts an isotropic degradation of thermal conductivity based on a scalar damage variable. In the previous model, the thermal conductivity is expressed as:

$$[eqn]$$

where $[eqn]$ is thresholded scalar damage, defined as:

$$[eqn]$$

Here, d is the original scalar damage defined by Equation (6), and $[eqn]$ and $[eqn]$ are two threshold values. This model assumes that the influence of cracks on thermal conductivity is the same in all directions (isotropic).

To account for directional effects, we propose an anisotropic thermal conductivity model based on the new damage tensor $[eqn]$ . A suitable way to define thermal conductivity $[eqn]$ has been given in Equation (21).

The calculation results based on Equations (30) and (31) and the proposed model based on Equation (21) for thermal insulation cracks are compared. Figure 10 shows the heat flux calculation results between the classical and the proposed model for thermal insulation cracks. Figure 11 shows the comparison of thermal conductivity $[eqn]$ and $[eqn]$ between the two models. Notably, the classical model calculates thermal conductivity based on a scalar damage value, resulting in equal thermal conductivity in both directions. In contrast, the new definition of the PD damage captures directional differences in thermal conductivity. Specifically, the pre-existing crack in the y-direction significantly affects the thermal conductivity in the x-direction but has minimal impact on the thermal conductivity in the y-direction. Conversely, the effect of cracks in the x-direction on thermal conductivity is opposite. Figure 12 compares the temperature field contours between the two models, further demonstrating the necessity and advantages of the proposed model.

5.3. Thermal Fracture of a Cruciform Plate with a Corner Crack

This example is a quasi-static crack propagation problem under thermo-mechanical conditions. Figure 13 depicts a cruciform plate with a corner crack, where the crack length is 10 mm and forms a $[eqn]$ angle with the vertical axis. The numerical computations are conducted under the assumption of plane stress. We consider crack propagation paths under three distinct mechanical and thermal boundary conditions. In all three conditions, the initial temperature is set to 0 °C, and the three edges of the plate are mechanically restrained against normal displacement. The material properties are shown in Table 2, while the boundary conditions for Figure 13 are detailed in Table 3. To accelerate the calculations, we introduce the PD domain only in the vicinity of the corner crack and employ the CCM model in the remaining domain. The PD domain is discretized into quadrilateral finite element meshes with dimensions of $[eqn]$ mm, while the CCM domain is discretized into finite element meshes with dimensions of $[eqn]$ mm. The PD horizon is set to $[eqn]$ .

Figure 14 displays the final crack paths under conditions 1 and 2 based on the present model (represented by colored contours), compared with those predicted using PFM proposed by Mandal et al. [26] (represented by purple dashed lines), showing good agreement between the two models. For condition 3, as illustrated in Figure 15, we compare the crack paths obtained from various methods: PFM by Mandal et al. [26], XFEM by Duflot et al. [27], the adaptive mesh refinement (AMR) method by Pham et al. [28], the boundary element method (BEM) by Prasad et al. [29], the gradient-enhanced damage method by Sarkar et al. [30], and the present method (represented by colored contours). As shown, the crack paths predicted using these methods are in good agreement, validating the capability of the present model for simulating thermo-mechanically coupled crack propagation. The new definition of PD damage, which can distinguish damage in two directions, enables a more accurate description of the effect of thermal insulation cracks on thermal conductivity. Figure 16a,b show the degradation of thermal conductivity $[eqn]$ and $[eqn]$ due to the thermal insulation crack. Figure 16c shows the temperature field.

5.4. Thermal Shock Fractures in Ceramics

Ceramic materials exhibit excellent mechanical properties at high temperatures, but they may break when subjected to sudden temperature changes. Consequently, thermal shock resistance is a crucial metric for assessing the suitability of ceramic materials for high-temperature engineering applications. The quenching test of ceramics is a widely adopted method to investigate the failure mechanisms associated with thermal shock-induced crack patterns. In the previous experiments, thin specimens measuring $[eqn]$ mm were heated to temperature $[eqn]$ and then rapidly quenched in a water bath maintained at $[eqn]$ °C [31]. Due to the central symmetry of the boundary conditions and geometric model, only a quarter of the numerical model, with dimensions of $[eqn]$ mm, needs to be simulated. Figure 17 presents a simplified two-dimensional numerical simulation sketch, illustrating the geometry and boundary conditions of the computational domain. The upper and right edges are constrained, while convection boundaries are applied to the left and bottom edges. The numerical computations are based on the plane stress assumption. Table 4 lists the material properties reported in [19,32]. A uniform finite element discretization with a size of $[eqn]$ mm is employed to calculate both the temperature and displacement fields. The PD horizon is set to $[eqn]$ . Figure 18 compares the final crack patterns obtained through numerical simulations and experiments [32]. Both methods demonstrate crack propagation under thermal shock conditions. The initial temperature $[eqn]$ and convective heat transform coefficient $[eqn]$ of the ceramic plates influence the final crack patterns. The left and middle columns of the figure show the new definition of the PD damage $[eqn]$ proposed in this paper for different $[eqn]$ values. The right column displays the final crack patterns obtained from experiments [32]. Figure 19 presents the final temperature field contours obtained through numerical simulations for various initial temperatures. To compare numerical simulations with experiments, the dimensionless crack length is selected as a reference index. The dimensionless crack length is defined as the ratio of the crack length to the total height (5 mm). Cracks with a dimensionless crack length greater than 0.6 are classified as long cracks. Figure 20 compares the numerical results from the proposed model with experimental results. The crack length in numerical simulations is slightly shorter than that observed in experiments. This discrepancy may be attributed to unpredictable manufacturing defects or material heterogeneities in the ceramic specimens used in the experiments.

6. Conclusions

A new definition of PD damage has been developed to model thermal fractures. A comprehensive solution framework has been proposed to compute both the temperature and displacement fields. Specifically, the temperature field is calculated via CCM, whereas the displacement field is computed via PD. To ensure accuracy, unified finite element discretization is employed for both methods. To elucidate the impact of thermal cracks or defects on the temperature field, the insulation crack is modeled by reducing thermal conductivity through the new definition of the PD damage, which indicates the influence of various bond distributions on thermal conductivity.

The proposed thermo-mechanical model utilizes finite element discretization. The shared mesh system for both temperature and displacement computation eliminates the need for remeshing, making it highly efficient and convenient for engineering applications. Furthermore, the new definition of PD damage differs from the original definition. The original definition can only represent the number of bond failures within the domain and cannot depict specific bond failure distributions. In contrast, the new definition of PD damage can characterize damage in different directions. Notably, the model presented in this paper is applicable not only to isotropic materials but also to anisotropic materials. For isotropic materials, the new PD damage is defined along the three coordinate axis directions, whereas for anisotropic materials, it is defined along the material’s three principal directions. Subsequent studies will provide a discussion.