Peridynamics Model with Surface Correction Near Insulated Cracks for Transient Heat Conduction in Functionally Graded Materials

Abstract

A peridynamic (PD) model of functionally graded materials (FGMs) is presented to simulate transient heat conduction in the FGM plate with insulated cracks. The surface correction is considered in the model to reduce the surface effect near the domain boundary and insulated cracks. In order to verify the proposed model, a numerical example for the FGM plate is carried out. The results show good agreement with the analytical solution. The convergence of the model with the surface correction for FGMs without cracks is then investigated. The results reveal that our model converges to the classical solutions in the limit of the horizon going to zero. The effects of two material points discretization schemes on the accuracy of numerical results are investigated. For transient heat conduction of FGMs with a static crack, the results obtained from the proposed PD model agree well with that from the finite element method. Finally, transient heat conduction of the FGM plate with a dynamic horizontal crack and intersecting cracks is simulated and discussed.

1. Introduction

Owing to their excellent mechanical and thermal properties, FGMs have attracted great interest in the field of materials and engineering [1]. The physical properties of FGMs are non-homogeneous and the material parameters present a continuously gradient change in the direction of thickness. Compared with traditional composite materials, FGMs have the advantage of no interfacial stress. In general, due to their superior physical properties, FGMs are used as structural components in aerospace and nuclear reactors in extreme high-temperature environments. Hence, transient heat conduction analysis is important for the design, optimization and engineering applications of FGMs [2].

Heat conduction of FGMs has been investigated by many researchers, and they mainly used analytical and numerical methods. Hosseini et al. [3] studied the transient thermal conduction of a functionally graded cylindrical shell assumed to be an axisymmetric boundary condition by using the Bessel functions. Zhao et al. [4] used perturbation method to calculate the transient temperature field of a functionally gradient ceramic plate under a convective boundary condition. Kayhani et al. [5] presented an exact analytical solution for the steady state heat conduction problem of cylindrical composite laminates. The solution can be directly applied to cylindrical composite pipes and reservoirs. Cinefra et al. [6,7] proposed models for the thermo-mechanical analysis of FGM shells and FEM thermal analysis of laminated shells by referring to Carrera’s Unified Formulation (CUF).

Since the non-homogeneous material properties, it is very difficult to obtain the analytical solution of transient heat conduction equation for FGMs. Therefore, numerical methods are widely used in the heat conduction problem of FGMs. Wang et al. [8] simulated temperature field and associated thermal stresses in FGMs with temperature-dependent material properties by a finite element/finite difference (FE/FD) method. Liu and Ming [9] developed a new three-dimensional control volume finite element method (CVFEM) for transient heat conduction in FGMs. Li et al. [2] used a multiple reciprocity boundary face method to solve the transient heat conduction problem of FGMs. Yu et al. [10] solved the transient heat conduction problem of FGMs by a new method consisting of differential transformation method and radial integral boundary element method. Zhang et al. [11] solved the two-dimensional transient heat conduction problems in FGMs by the numerical manifold method (NMM). Zhou et al. [12] employed the meshless weighted least-square (MWLS) method to solve the heat conduction problem for irregular FGMs with temperature-dependent material properties. Xi et al. [13] proposed a semi-analytical boundary collocation solver for inverse Cauchy problems in heat conduction under 3D FGMs with a heat source. Krahulec et al. [14] used the meshless local radial basis function method to simulate stationary and transient heat conduction problems in FGMs.

Under extreme thermal shock conditions, crack propagation is easily induced inside the FGM [15]. The inherent difficulty of this problem is that the spatial derivatives required for partial differential equations do not exist on a crack tip or surface. Therefore, any numerical method derived from these equations inherits this difficulty in modeling cracks [16]. The PD method is a nonlocal formulation of continuum mechanics, firstly proposed by Silling [17]. This theory discretizes the continuum into material points and assumes that the material points interact with each other in a finite range. Based on the level of interaction between the two points, damage is incorporated into the theory, so localization and fracture occur naturally [16]. After nearly 20 years of development, the peridynamic theory has been widely accepted to simulate and predict the damage and fracture of materials and structures [18,19,20,21,22].

A peridynamic approach to heat conduction has advantages over other continuum theory-based methods, as it not only accounts for nonlocality but it also allows for the determination of the temperature field including the discontinuities. In the context of peridynamic theory, Bobaru and Duangpanya [23] proposed 1-D transient heat conduction. They examined the effects of constant and triangular micro-conductivity functions on the heat transfer. Actually, they attributed micro-conductivity as a feature inherited in the peridynamic model of heat conduction. They extended their work to 2-D space later, and numerically proved that when the nonlocal parameter in the peridynamic model goes to zero, the results are converged to the classical solution. The numerical calculation considered the self-transfer and volume correction of the material points. The temperature load was applied to the material point closest to the boundary, and a two-dimensional transient heat conduction problem with adiabatic cracks was simulated [24]. Agwai [25] gave three different PD heat transfer kernel functions, and the corresponding constant micro-conductivity is obtained by equalizing the thermal potential energy function. The 1-D numerical calculation shown that the surface correction factor reduces the PD numerical calculation error. Based on the Lagrangian formula and PD theory, Oterkus et al. [26] developed the state-based PD heat conduction equation. They presented a thermal response function, a method for determining the micro-conductivity and the method of the virtual layer temperature boundary condition, but their examples still using the bond-based formulation. In another work, Chen and Bobaru [27] focused on the convergence issue regarding numerical solution of the peridynamic heat transfer equation in 1-D. They employed a one-point Gauss quadrature method to solve the equation numerically and showed that type of the peridynamic kernel affects the convergence. Shojaei et al. [28] developed a simple switching technique to couple PD with methods based on the classical local theory near the boundaries to remove the surface effect from the Peridynamic solution. Bazazzadeh et al. [29] used adaptive grid refinement techniques close to the boundaries to reduce the surface effect for thermal shock problems in ceramics.

The PD heat conduction method does not include spatial derivatives and uses instead an integral over a region around a material point. The new model is suitable for modeling, for example, heat flow in bodies with evolving discontinuities such as growing insulated cracks. However, the effect of the surface correction factor on the accuracy and convergence of numerical calculation results is still unclear, especially for 2-D transient heat conduction problem. For the problem of heat conduction with a crack, no one has considered the surface correction near the crack.

To address this concern, the aim of this paper is focused to analyze the effect of surface correction on transient heat conduction of the FGM plates. This work is organized into four sections. In Section 2, we review the peridynamic formulation for transient heat conduction and describe the PD heat transfer kernel function formulas with constant and conical micro-conductivity. In Section 3, based on the kernel function of constant micro-conductivity, a discrete model for transient heat conduction in FGMs with cracks is given. Near a domain boundary and insulated cracks, we propose a surface correction factor to reduce surface effect. In Section 4, the effects of surface correction on the accuracy of PD transient thermal conduction model of FGMs are studied. Moreover, the convergence of the FGM model after surface correction is analysed and the effects of two material points discretization schemes on the accuracy of numerical results are investigated. In Section 5, the transient heat conduction problem of FGMs with static and dynamic cracks are simulated. The influence of the surface correction about the material points close to an insulated crack on the numerical results is analysed. In the last section, some conclusions are presented.

2. Bond-Based Peridynamic Thermal Transfer

As illustrated in Figure 1a, considering a body that occupies a region $\Omega$ consisting of material points, thermal transfer exists between each material point $\mathit{x}$ and the associated neighboring points $\widehat{\mathit{x}}$ within a finite domain, which is named horizon H. Each material point has the associated mass and volume. The thermal conductor between two interacting material points $\mathit{x}$ and $\widehat{\mathit{x}}$ is called ‘thermal bond’ or in short: t-bond (see Figure 1b). The peridynamic heat flux per unit volume transferring along the t-bond is assumed to depend on the relative location and the temperature difference between material points $\mathit{x}$ and $\widehat{\mathit{x}}$. The transient heat transfer equation in the form of a bond can be written as [26]:

$$\rho c_v\dot{\Theta }\left(\mathit{x},t\right)={{\displaystyle \int }}_H^{}{f}_h\left(\widehat{\mathit{x}},\mathit{x},t\right)d{V}_{\widehat{\mathit{x}}}+h_s\left(\mathit{x},t\right)$$

(1)

where $\rho$ and $c_v$ are the density and specific heat capacity of the material point $\mathit{x}$, respectively, and $h_s$ represents the heat generated per unit volume of the internal heat source. The kernel function (the function under the integral sign) $f_h\left(\widehat{\mathit{x}},\mathit{x},t\right)$ represents the heat exchange per unit volume of the material point $\widehat{\mathit{x}}$ to the material point $\mathit{x}$, usually expressed in the following form [27]:

$$f_h\left(\widehat{\mathit{x}},\mathit{x},t\right)=\{\begin{array}{c}\kappa \left(\mathit{x},\widehat{\mathit{x}}\right)\frac{\Theta \left(\widehat{\mathit{x}},t\right)-\Theta \left(\mathit{x},t\right)}{{\Vert \widehat{\mathit{x}}-\mathit{x}\Vert }^n}\\ 0, \Vert \widehat{\mathit{x}}-\mathit{x}\Vert >\delta \end{array}, \Vert \widehat{\mathit{x}}-\mathit{x}\Vert \le \delta$$

(2)

where $\kappa \left(\mathit{x},\widehat{\mathit{x}}\right)$ represents the micro-conductivity of the t-bond between material points $\widehat{\mathit{x}}$ and $\mathit{x}$. $\Theta \left(\widehat{\mathit{x}},t\right)-\Theta \left(\mathit{x},t\right)$ is the temperature difference between the two material points $\widehat{\mathit{x}}$ and $\mathit{x}$. Here n is an integer, normally selected to be 0, 1, or 2 [27]. $\delta$ is the horizon size.

As described by Bobaru and Duangpanya [24], the PD heat flux vector between the material points $\mathit{x}$ and $\widehat{\mathit{x}}$ can be written as:

$$\mathit{q}\left(\mathit{x},t\right)=K\left(\mathit{x},\widehat{\mathit{x}}\right)\frac{\Theta \left(\widehat{\mathit{x}},t\right)-\Theta \left(\mathit{x},t\right)}{\Vert \widehat{\mathit{x}}-\mathit{x}\Vert }\mathit{e}$$

(3)

where $K\left(\mathit{x},\widehat{\mathit{x}}\right)$ is the heat conductivity of the t-bond between material points $\mathit{x}$ and $\widehat{\mathit{x}}$, and where $\mathit{e}=\mathit{i}cos\varphi +\mathit{j}sin\varphi$ is the unit vector along the corresponding t-bond (see Figure 1b).

In the 2-D case, the classical heat flux vector between $\mathit{x}$ and $\widehat{\mathit{x}}$ along the t-bond because of their temperature difference can be written as [24]:

$${\mathit{q}}^{classical}\left(\mathit{x},t\right)=K_x\left(\mathit{x},\widehat{\mathit{x}}\right)\frac{\Theta \left(\widehat{\mathit{x}},t\right)-\Theta \left(\mathit{x},t\right)}{\Vert \widehat{\mathit{x}}-\mathit{x}\Vert }cos\varphi \mathit{i}+K_y\left(\mathit{x},\widehat{\mathit{x}}\right)\frac{\Theta \left(\widehat{\mathit{x}},t\right)-\Theta \left(\mathit{x},t\right)}{\Vert \widehat{\mathit{x}}-\mathit{x}\Vert }sin\varphi \mathit{j}$$

(4)

where $K_x\left(\mathit{x},\widehat{\mathit{x}}\right)$ and $K_y\left(\mathit{x},\widehat{\mathit{x}}\right)$ are the thermal conductivities of the t-bond ($\mathit{x},\widehat{\mathit{x}}$) respectively in x and y directions.

In the bond-based peridynamics transient heat transfer theory, assuming that the thickness of the plate is h, the PD heat flux between $\mathit{x}$ and $\widehat{\mathit{x}}$ along the t-bond can be obtained from [24]:

$${\mathit{q}}^{Peri}\left(\mathit{x},t\right)={{\displaystyle \int }}_{H^{+}}^{}\Vert \widehat{\mathit{x}}-\mathit{x}\Vert f_h\left(\widehat{\mathit{x}},\mathit{x},t\right)d{V}_{\widehat{\mathit{x}}}={{\displaystyle \int }}_{H^{+}}^{}\kappa \left(\mathit{x},\widehat{\mathit{x}}\right)\frac{\Theta \left(\widehat{\mathit{x}},t\right)-\Theta \left(\mathit{x},t\right)}{{\Vert \widehat{\mathit{x}}-\mathit{x}\Vert }^{n-1}}d{V}_{\widehat{\mathit{x}}}$$

(5)

The natural value for the parameter n in Equation (5) was 2 in the literature [24]. Here we let n be 0, 1 and 2 to get the different types of kernels mentioned before in the 2-D case.

Two kinds of micro-conductivity functions are considered. The first one is a constant:

$$\kappa \left(\mathit{x},\widehat{\mathit{x}}\right)={\kappa }_0\left(\mathit{x},\widehat{\mathit{x}}\right)$$

(6)

Another one varies linearly over the horizon size called conical micro-conductivity function:

$$\kappa \left(\mathit{x},\widehat{\mathit{x}}\right)={\kappa }_1\left(\mathit{x},\widehat{\mathit{x}}\right)(1-\frac{\Vert \widehat{\mathit{x}}-\mathit{x}\Vert }{\delta })$$

(7)

The micro-conductivity is obtained by equating the peridynamic heat flux with the classical heat flux. The temperature field is assumed to be linearly distribution of the following form $\Theta \left(x,y,t\right)=ax+b$, where $a$ and $b$ are constants, as shown in Figure 2. The classical heat flux and peridynamic heat flux perpendicular to the isothermal line at the point x are respectively expressed as:

$${\mathit{q}}^{classical}\left(\mathit{x},t\right)=K\left(\mathit{x},\widehat{\mathit{x}}\right)\frac{\partial \Theta \left(x,y\right)}{\partial x}=K\left(\mathit{x},\widehat{\mathit{x}}\right)a$$

(8)

.

With the constant shape factor, substituting Equation (6) into Equation (5), the PD heat flux can be written as:

$${\mathit{q}}^{Peri}\left(\mathit{x},t\right)=\frac{{\delta }^{4-n}}{2\left(4-n\right)}{\kappa }_0\left(\mathit{x},\widehat{\mathit{x}}\right)ah\pi$$

(9)

With the ‘‘triangular’’ shape factor, substituting Equation (7) into Equation (5), the PD heat flux can be written as:

$${\mathit{q}}^{Peri}\left(\mathit{x},t\right)=\frac{{\delta }^{4-n}}{2\left(4-n\right)\left(5-n\right)}{\kappa }_1\left(\mathit{x},\widehat{\mathit{x}}\right)ah\pi$$

(10)

Equating (9), (10) with the classical flux in (8), gives:

$${\kappa }_0\left(\mathit{x},\widehat{\mathit{x}}\right)=\frac{2\left(4-n\right)}{\pi h{\delta }^{4-n}}K\left(\mathit{x},\widehat{\mathit{x}}\right)$$

(11)

$${\kappa }_1\left(\mathit{x},\widehat{\mathit{x}}\right)=\frac{2\left(4-n\right)\left(5-n\right)}{\pi h{\delta }^{4-n}}K\left(\mathit{x},\widehat{\mathit{x}}\right)$$

(12)

For the conical micro-conductivity, $\kappa \left(\mathit{x},\widehat{\mathit{x}}\right)$ decreases linearly with the bond length, making the PD model closer to the classical localized model, reducing the integral error caused by the PD surface effect. In this paper, in order to make the integral error more obvious, the constant micro-conductivity and the integer $n=1$ is adopted. In the next section, a peridynamic model for transient heat conduction in the FGM with cracks is given.

3. A Peridynamic Transient Heat Conduction Model for FGMs

3.1. Numerical Discretization of FGMs in 2D

Functionally graded materials are basically heterogeneous materials, which are considered as composites composed of two-phase materials, the volume fraction of which gradually changes from one end to the other. For the convenience of analysis, we assume that the material parameters such as density, specific heat capacity, and heat transfer coefficient of FGMs are functions of spatial coordinates, and the local material parameters can be regarded as fixed constants. They can be expressed as:

$$\rho \left(x,y\right)=P\left({\rho }_0,x,y\right)$$

(13)

$$c\left(x,y\right)=P\left(c_0,x,y\right)$$

(14)

$$K\left(x,y\right)=P\left(K_0,x,y\right)$$

(15)

where P represents the gradient function of the material parameter. Cheng et al. proposed a peridynamic model of FGMs, which calculated the peridynamic micro-modulus of the bond by the mean values of the material properties at the two interacting material points [30]. The equivalent micro-conductivity ${\kappa }_{eff}$ of the t-bond between the two material points $\mathit{x}$ and $\widehat{\mathit{x}}$ can be approximated as:

$${\kappa }_{eff}\left(\mathit{x},\widehat{\mathit{x}}\right)=\frac{6}{\pi h{\delta }^3}{K}_{eff}\left(\mathit{x},\widehat{\mathit{x}}\right)$$

(16)

where $K_{eff}\left(\mathit{x},\widehat{\mathit{x}}\right)=0.5\left(K\left(x,y\right)+K\left(\widehat{x},\widehat{y}\right)\right)$. The bond-based PD transient heat conduction equation of FGMs can be discretely expressed as:

$$\rho \left({\mathit{x}}_{\mathit{i}}\right)c_v\left({\mathit{x}}_{\mathit{i}}\right)\dot{\Theta }\left({\mathit{x}}_{\mathit{i}},t\right)={\displaystyle \sum }_{j=1}^n{\kappa }_{eff}\left({\mathit{x}}_{\mathit{i}},{\mathit{x}}_{\mathit{j}}\right)\frac{\Theta \left({\mathit{x}}_{\mathit{j}},t\right)-\Theta \left({\mathit{x}}_{\mathit{i}},t\right)}{\Vert {\mathit{x}}_{\mathit{j}}-{\mathit{x}}_{\mathit{i}}\Vert }{V}_{{\mathit{x}}_{\mathit{j}}}+h_s\left({\mathit{x}}_{\mathit{i}},t\right)$$

(17)

where $i$ represents the point of interest, $j$ represents the family members of point $i$ and $n$ is the total number of all points $j$.

3.2. Dirichlet Boundary Condition

In general, the PD method causes the boundary effect by applying a fixed temperature directly to the material points at the boundary [23,24]. This boundary condition application method will result in incomplete neighborhood of the material points near the boundary, thereby weakening the boundary condition. This paper refers to the method of the fictitious nodes [26] when dealing with the Dirichlet boundary conditions. Based on his research, the extent of the virtual boundary layer to be equal to the horizon, δ in order to ensure that the prescribed temperatures sufficiently reflected in the actual material region. As shown in Figure 3, Dirichlet boundary condition, $\Theta \left({\mathit{x}}^{*},t\right)$ can be applied in the virtual boundary layer, ${ℛ}_t$ along the boundary of the body surface, as:

$$\Theta \left({\mathit{y}}_{\mathit{i}},t+\mathrm{∆}t\right)=2\Theta \left({\mathit{x}}^{*},t+\mathrm{∆}t\right)-\Theta \left({\mathit{z}}_{\mathit{i}},t\right),{\mathit{x}}^{*}\in S_t, \mathit{y}\in {ℛ}_t,\mathit{z}\in ℛ, i=1,2,3$$

(18)

where ${\mathit{z}}_{\mathit{i}}$ represents the location of a material point in $ℛ$, and ${\mathit{x}}^{*}$ is the position of a point on the body surface, $S_t$. $\Theta \left({\mathit{z}}_{\mathit{i}},t\right)$ represents the temperature of the material point in the body, which is obtained by solving the PD transient heat conduction equation. $\Theta \left({\mathit{y}}_{\mathit{i}},t\right)$ represents the temperature of the material point in the virtual layer. When time equals $t+\mathrm{∆}t$, $\Theta \left({\mathit{x}}^{*},t+\mathrm{∆}t\right)$ is known from the analytical boundary condition, so we can get the $\Theta \left({\mathit{y}}_{\mathit{i}},t+\mathrm{∆}t\right)$ by the Equation (1). During numerical simulation, $\Theta \left({\mathit{y}}_{\mathit{i}},t\right)$ can be calculated at each time step in the virtual boundary layer.

3.3. Modeling of Insulated Crack

To simulate transient heat conduction for FGMs with cracks, cracks need to be modelled by PD. As shown in Figure 4, the partial t-bonds in the horizon of material point $\mathit{x}$ were broken. Referring to Bobaru and Duangpanya [24], the kernel function $f_n\left(\widehat{\mathit{x}},\mathit{x},t\right)$ in Equation (2) is set as a “history-dependent” function:

$${\tilde{f}}_n\left(\widehat{\mathit{x}},\mathit{x},t\right)=f_n\left(\widehat{\mathit{x}},\mathit{x},t\right)\mu (t,\Vert \widehat{\mathit{x}}-\mathit{x}\Vert )$$

(19)

where $\mu$ is a history-dependent scalar valued function, which is defined as:

$$\mu (t,\Vert \widehat{\mathit{x}}-\mathit{x}\Vert )=\{\begin{array}{c}1 if bond is intact \\ \alpha \left(t,\xi ,\eta \right) otherwise \end{array}$$

(20)

.

In this paper, it is assumed that the crack is an insulated crack, and the heat conduction of the t-bonds crossing the crack were completely suppressed. The scalar valued function $\alpha$ is zero if the t-bond crosses the insulated crack.

3.4. Surface Correction

The micro-conductivity parameters are determined by equating the peridynamic heat flux values with the classical heat flux values. The micro-conductivity parameters are calculated by assuming integration over a full horizon. Near a domain boundary, the peridynamic horizon is not complete. Therefore, the value of $\kappa$ needs to be corrected if the material points near the free surfaces or material interface.

A number of methods/algorithms have been proposed recently for surface correction. Le and Bobaru [30] discuss in detail the accuracy and computational efficiency of multiple surface correction methods for elasticity and fracture. These methods include the volume method [31], the force density method [32], the energy method [32], the force normalization [33], the position-aware method [34], and fictitious nodes method [26,32]. For transient heat conduction problem, we reference the surface correction factor $g\left(\mathit{x}\right)$ proposed in [26] to reduce the PD surface effect based on the energy method. It is obtained by comparing the PD thermal potential with the corresponding classical thermal potential energy for simple temperature distribution.

The classical thermal potential formulation at point $\mathit{x}$ can be written as [26]:

$$Z_{\infty }\left(\mathit{x}\right)=\frac{1}{2}\mathit{G}·K\left(\mathit{x}\right)\mathit{G}$$

(21)

where $K\left(\mathit{x}\right)$ is the thermal conductivity and $\mathit{G}=\nabla \Theta$.

The PD thermal potential is the integral over all microthermal potentials associated with this point, and is expressed as [26]:

$$Z\left(\mathit{x}\right)=\frac{1}{2}{{\displaystyle \int }}_H^{}{\kappa }_{eff}\left(\mathit{x},\widehat{\mathit{x}}\right)\frac{{\left(\Theta \left(\widehat{\mathit{x}},t\right)-\Theta \left(\mathit{x},t\right)\right)}^2}{2\Vert \widehat{\mathit{x}}-\mathit{x}\Vert }d{V}_{\widehat{\mathit{x}}}$$

(22)

where ${\kappa }_{eff}$ is the equivalent micro-conductivity of the t-bond between the two material points $\mathit{x}$ and $\widehat{\mathit{x}}$. The correction factor of the material points in the body can be determined by [26]:

$$g\left(\mathit{x}\right)=\frac{Z_{\infty }\left(\mathit{x}\right)}{Z\left(\mathit{x}\right)}$$

(23)

Therefore, the discretized thermal diffusion Equation (17) containing the surface correction factor for point $\mathit{x}$ can be written:

$$\rho \left({\mathit{x}}_{\mathit{i}}\right)c_v\left({\mathit{x}}_{\mathit{i}}\right)\dot{\Theta }\left({\mathit{x}}_{\mathit{i}},t\right)={\displaystyle \sum }_{j=1}^{n}g\left({\mathit{x}}_{\mathit{i}},{\mathit{x}}_{\mathit{j}}\right){\kappa }_{eff}\left({\mathit{x}}_{\mathit{i}},{\mathit{x}}_{\mathit{j}}\right)\frac{\Theta \left({\mathit{x}}_{\mathit{j}},t\right)-\Theta \left({\mathit{x}}_{\mathit{i}},t\right)}{\Vert {\mathit{x}}_{\mathit{j}}-{\mathit{x}}_{\mathit{i}}\Vert }{V}_{{\mathit{x}}_{\mathit{j}}}+h_s\left({\mathit{x}}_{\mathit{i}},t\right)$$

(24)

where $g\left(\mathit{x},\widehat{\mathit{x}}\right)=\left(g\left(\mathit{x}\right)+g\left(\widehat{\mathit{x}}\right)\right)/2$.

For the material points close to insulated cracks, the peridynamic horizon is incomplete due to the partial t-bond disconnection (see Figure 4). In principle, in order to match the conductivity parameter of a material point with complete horizon while integrating over a smaller area (due to the presence of an insulating crack), one has to increase the micro-conductivity value of the t-bond to compensate for the reduction of integration area.

Considering the partial t-bond of the material point $\mathit{x}$ in the horizon is broken, the thermal potential energy expressed by the Equation (22) can be written as:

$$Z\left(\mathit{x}\right)=\frac{1}{2}{{\displaystyle \int }}_H^{}\mu (t,\Vert \widehat{\mathit{x}}-\mathit{x}\Vert ){\mathsf{\kappa }}_{eff}\left(\mathit{x},\widehat{\mathit{x}}\right)\frac{{\left(\Theta \left(\widehat{\mathit{x}},t\right)-\Theta \left(\mathit{x},t\right)\right)}^2}{2\Vert \widehat{\mathit{x}}-\mathit{x}\Vert }d{V}_{\widehat{\mathit{x}}}$$

(25)

The classical thermal potential is shown in Equation (21), so the surface correction factor can also be expressed by Equation (23). When dealing with transient heat conduction problems with dynamic cracks, new t-bonds will break near the crack as the crack grows. The PD thermal potential energy expressed by Equation (25) will change. The correction factor needs to be calculated at each step according to Equation (23). However, for transient thermal conduction problems with static cracks or without cracks, this parameter only needs to be calculated in the first time step and remains constant after the first step.

4. Numerical Convergence Studies

We analyse the convergence of the proposed peridynamic model for FGMs by the following example of a 2-D FGM plate without any crack. The convergence tests in the limit of nonlocal parameter going to zero are performed by comparing the peridynamic solutions with the classical analytic solution.

4.1. The Peridynamic Model for FGMs

In the literature [19], the bond-based PD model for dynamic fractures in FGMs was verified by elastic wave propagation and comparisons with analytical results for the classical model. In order to verify the reliability of the PD heat conduction model of FGMs proposed in Section 3, the temperature field and heat flux of the following example are calculated and compared with the analytical solution.

As shown in Figure 5a, the FGM plate has a length and width of 1.0 m, a thickness of 0.01 m, and an initial temperature of 0 °C. The top edge is fixed at a constant temperature of 100 °C. The bottom edge is 0 °C. The left and right edges are well insulated (zero normal flux). The geometry details and the boundary conditions are shown in Figure 5a. Figure 5b is the PD discretization model of the FGM plate.

The thermal conductivity and the specific heat of FGMs are defined by

$$K\left(x,y\right)=5.0e^{\lambda y} W/\left(m·K\right)$$

(26)

$$c\left(x,y\right)=1.0e^{\lambda y} J/\left(kg·K\right)$$

(27)

where exponential parameters $\lambda =3$ m⁻¹, and the density $\rho =1.0 \mathrm{kg}/{\mathrm{m}}^3$.

The analytical solution for temperature filed is given in [2] by

$$\begin{array}{c}\Theta \left(x,y,t\right)=T\frac{1-e^{-2\beta y}}{1-e^{-2\beta L}}+{\displaystyle {\displaystyle \sum }_{n=1}^{\infty }}{B}_{n}sin\frac{n\pi y}{L}{e}^{-\beta y}{e}^{-(n^2{\pi }^2/L^2+{\beta }^2)\alpha t}\\ B_n=\frac{2T{e}^{\beta L}n\pi cos\left(n\pi \right)}{{\beta }^2{L}^2+n^2{\pi }^2}\end{array}$$

(28)

where $L$ is the length of the plate in the y-direction, $\beta =\lambda /2=1.5$ m⁻¹, $\alpha =5.0$. The distribution of heat flux can be obtained by the derivative of temperature with respect to time

$$\begin{array}{c}\mathit{q}\left(x,y,t\right)=2K\left(x,y\right)\beta T\frac{e^{-2\beta y}}{1-e^{-2\beta L}}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} +{\displaystyle {\displaystyle \sum }_{n=1}^{\infty }}K\left(x,y\right)B_n(\frac{n\pi }{L}cos\frac{n\pi y}{L}-\beta sin\frac{n\pi y}{L})e^{-\beta y}{e}^{-(n^2{\pi }^2/L^2+{\beta }^2)\alpha t}\hfill \end{array}$$

(29)

The x and y directions of the model are discretized into 100 material points, the point spacing $\mathrm{∆}=0.01 \mathrm{m}$, the horizon size $\delta =3\mathrm{∆}$, and the time step $\mathrm{∆}t={10}^{-5}$ s. In Figure 6, the temperature distribution along the line x $=0.495$ m is compared with the solution obtained by analytical method. The solution is plotted at different times: 0.002, 0.005, 0.01, 0.02, 0.05, and 0.1 s, to observe its transient heat transfer characteristics. The vertical heat flux comparison at t = 0.05 s along the line x $=0.495$ m is plotted in Figure 7. The numerical results agree very well with the analytical solution.

4.2. Convergence Study

In order to investigate the effect of surface correction on the numerical accuracy of peridynamic heat conduction model for FGMs, the relative difference of our nonlocal solutions with the classical solution introduced in [24] can be defined as:

$$\frac{\Vert {\theta }_{classical}-{\theta }_{peri}\Vert }{\Vert {\theta }_{classical}\Vert }$$

(30)

We compared the relative differences in temperature between the corrected peridynamic model and the uncorrected peridynamic model. As we can see from Figure 8, the relative difference after correction is obviously smaller than uncorrected PD. Therefore, it can be said that the correction of the surface effect improves calculation accuracy.

In order to observe the convergence behavior of the PD transient heat conduction model after surface correction to the classical analytical solution of the FGM plate, the m convergence and δ convergence (see [35]) of this model are studied respectively. The temperature of the material point with the coordinates of (0.495 m, 0.495 m) in the middle of the FGM plate at t = 0.01 s is used to evaluate convergence of the PD solutions to the classical solution.

In terms of m-convergence, as shown in Figure 9, the fixed horizon $\delta$ is equal to 0.04, 0.02, and 0.01 m, the results show that when $\delta$ is constant, the relative difference tends to decrease with increasing the nonlocal ratio, $m=\delta /\mathrm{∆}$. In terms of δ-convergence, the nonlocal ratio m is equal to 2, 4, and 8, it can be seen that the relative difference tends to decrease with decreasing horizon size δ in Figure 10. When δ = 0.01 m and m = 4, 8 and δ = 0.02 m and m = 8, the numerical calculation has higher precision, and the relative difference is less than 1%. Another type of m-convergence is also used for analysis, the point spacing $\mathrm{∆}$ is fixed and the nonlocal ratio m is changed. As shown in Figure 11, we analyse the m-convergence for three different point spacing: $\mathrm{∆}$ = 0.02, 0.01, and 0.005 m. The relative difference of the numerical calculation changes very little with increasing m, while the point spacing is fixed. The relative difference when m = 3 is slightly less than when m is equal to some other values. The relative difference decreases as the material point spacing decreases.

In PD analysis, the computational time increases with the point spacing decreasing and the nonlocal ratio increasing. Thus, it is important to use the optimal point spacing for the desired computational efficiency and accuracy. Based on the convergence study, the family of each point is established by the horizon size of 3 associated with the grid spacing of 0.01 m.

4.3. Discretization Schemes for PD Material Points

The domain is divided into square grids, with each material point at the center of the grid. There are two ways to assign the material points positions. One of these is shown in Figure 5b, the center of the material points closest to the boundary is not on the boundary of the FGM plate, but inside the body. As shown in Figure 12a, it is called discretization scheme Ⅰ. When Dirichlet boundary condition is applied to the boundary of the FGMs plate, the temperature is actually applied to the material point with a certain volume and it has an error with the true Dirichlet boundary condition. Based on the above reason, we use discretization scheme Ⅱ to keep the center of the material point on the sample boundary, as shown in Figure 12b. The other parameters are the same as in Section 4.1. and the calculated result is compared with the scheme Ⅰ.

As we can see from Figure 13, the PD scheme Ⅱ has a lower relative difference than the scheme Ⅰ, and the relative error of each position point is smaller.

Table 1. Comparison of the temperature of the discretization scheme Ⅱ with the analytical solution along the line x = 0.5 m at t = 0.01 s, t = 0.02 s.

y (m)	t = 0.01 s			t = 0.02 s
	Analytical (°C)	PD (°C)	Relative Difference	Analytical (°C)	PD (°C)	Relative Difference
0.0	0.0000	0.0000	0.0000%	0.0000	0.0000	0.0000%
0.1	1.3779	1.3832	0.3857%	9.9293	9.9394	0.1015%
0.2	3.4151	3.4245	0.2764%	18.8305	18.8416	0.0589%
0.3	7.0249	7.0388	0.1974%	28.0408	28.0538	0.0464%
0.4	13.0901	13.1088	0.1427%	38.1238	38.1393	0.0406%
0.5	22.3446	22.3686	0.1077%	49.0987	49.1169	0.0370%
0.6	35.0706	35.1008	0.0861%	60.6112	60.6320	0.0344%
0.7	50.7914	50.8282	0.0724%	72.0826	72.1056	0.0319%
0.8	68.1695	68.2119	0.0622%	82.8518	82.8759	0.0291%
0.9	85.2493	85.2942	0.0527%	92.3084	92.3321	0.0258%
1.0	100.0000	100.0000	0.0000%	100.0000	100.0000	0.0000%

is the comparison of the temperature between the discretization scheme Ⅱ and the analytical solution at t = 0.01 s, t = 0.02 s along x = 0.5 m line, and the relative differences of the scheme Ⅱ are all below 0.5%. Therefore, the numerical solution obtained by the scheme Ⅱ which align the material points on the sample’s boundary is closer to the analytical solution.

5. Transient Heat Conduction of the FGM Plate with Cracks

5.1. Transient Heat Conduction of the FGM Plate with Static Crack

As shown in Figure 14, the FGM plate has a horizontal static insulated crack with a length of 0.5 m. Material properties, the geometry details and the boundary conditions are the same as the previous example.

The example is discretized by the scheme Ⅱ and the model (PD2) considering the micro-conductivity correction of the material points close to the crack is adopted. Results from the PD1 (the PD model without considering the micro-conductivity correction of material points close to the crack) and the classical finite element method (FEM) are used for comparison.

Figure 15 shows the temperature field of FEM result at t = 0.03 s. Figure 16 shows the temperature distribution of PD1 and PD2 at t = 0.03 s. At this time, the temperature in the FGM plate has not reached a steady state. It can be observed from Figure 15 and Figure 16 that the temperature contours of PD2 and FEM are almost identical. But, the temperature of PD1 at the bottom of the crack is lower than the temperature obtained by FEM. This phenomenon can also be observed by comparing the temperature results of the FGM plate along the line x = 0.5 m at t = 0.03 s (see Figure 17). The reason for the above results can be analyzed as follows. The micro-conductivity of the material points close to crack from the PD1 method is lower than the true value, which leads to a decrease in the heat transfer efficiency. Hence, our proposed PD model which considers the micro-conductivity correction of the material points close to the crack can simulate transient heat conduction of FGMs plate with insulated crack exactly.

5.2. Transient Heat Conduction of the FGM Plate with Dynamic Horizontal Crack

In the practical coupled thermo-mechanical engineering problems, the initiation and propagation of cracks is often abrupt and uncontrollable. The classical numerical method is difficult to deal with this sudden spatial discontinuity problem because spatial derivatives don’t exist over discontinuities. Peridynamics fundamentally eliminates these difficulties since spatial derivatives do not appear in the formulation.

To demonstrate the advantage of PD, we consider the problem of transient heat conduction of the FGM plate with a horizontal dynamic crack. The growth of the crack is imposed by us here. As shown in Figure 18, consider the pre-existing adiabatic horizontal crack in the middle of the FGM plate. The initial length of the crack is 0.1 m. The crack can be dynamically extended along the x direction. The rate of 80 m/s is extended to both sides, and the final crack length is 0.5 m at t = 0.05 s. The material parameters of the FGM plate are the same as example 1. The geometry details and the boundary conditions are the same as the first example. It is assumed that the extended crack can cut the heat conduction of the t-bond passing through the crack and we consider the micro-conductivity correction of the material points close to crack introduced in Section 3.4.

The solutions are shown in Figure 19, the temperature distributions of the FGM plate with growing crack are shown at 0.01 s, 0.02 s, 0.03 s, 0.04 s and 0.05 s, respectively. With the increase of time, the dynamic crack become longer and hinders the heat conduction of the t-bond passing through the crack. As observed in these Figures, it is obvious that the presence of the crack causes discontinuous temperature fields above and below the crack surface. The temperature values below the crack drop suddenly because the insulated cracks behave as a thermal barrier in the FGM plate. It is also observed that the temperature gradient near the end of the crack is higher than its neighborhood. According to the thermodynamics knowledge, the position with larger temperature gradient always generates larger thermal stress, then leads to crack growth. It is worth noting that the shielding effect induced by the dynamic horizontal crack is correctly captured by the peridynamic simulation.

5.3. Transient Heat Conduction of the FGM Plate with Dynamic Intersecting Cracks

Under thermal shock, FGMs will generate complex cracks, including some interacting cracks. As shown in Figure 20, we consider the problem of transient heat conduction for the FGM plate in which two insulating cracks initiation, propagation and intersect. The geometry of the plate is the same as in the first example. The growth of two cracks is imposed by us. The two cracks begin to grow at time t = 0 s with a constant extension velocity 0.2828 m/s. A symmetric x-like crack is obtained at the end of the simulation with the spacing $a=0.25 \mathrm{m}$. The initial temperature of the FGM plate is 0 °C. The left and right boundaries are insulated and temperatures of ±100 °C are imposed on the top and bottom boundaries, respectively (see Figure 20). To investigate the influence of the material gradation on the temperature variation, three different exponential parameters $\lambda =$ 0, 1 and 3 ${\mathrm{m}}^{-1}$ are chosen for comparison in numerical calculation of the FGMs.

The temperature distributions at different times for three values of the exponential parameter $\lambda$ are shown in Figure 21, Figure 22, Figure 23 and Figure 24. The steady-state results are given at t = 5 s. As the cracks grow and the crack gap closes, the temperature rises above the V-shaped crack and drops below it. Once the cracks intersect with each other and continue to grow, the temperature values jump around the intersection point of x-like crack in FGM plate. Before the steady state is reached, the temperature difference between the top and bottom of the intersection point increases with the calculation time. When $\lambda =$ 0, the material becomes homogeneous and the temperature is symmetrically up and down when t = 5 s.

Comparing the results of three different exponential parameters, it is found that the temperature increases along with an increase in $\lambda$-values (or equivalently in thermal conductivity) at the same moment. It can also be found that as the $\lambda$-values increase, the temperature gradient at the crack tip of the FGM plate decreases before the cracks intersect, but it increases after the cracks intersect. From these pictures, difference between the homogeneous material and the FGMs can be observed. The shielding effect induced by the dynamic intersecting cracks is correctly captured by the peridynamic solution.

6. Conclusions and Future Work

In this paper, a peridynamics model for transient heat conduction in functionally graded materials with insulated cracks was introduced. To eliminate the peridynamic surface effect near a domain boundary and insulated cracks, the improved surface correction factor containing the crack was presented by referring to the original surface correction factor.

Subsequently, the present model was verified by comparison with analytical solutions for FGMs without cracks. We compared the relative differences of temperature between the corrected PD model and the uncorrected PD model, it was observed that the numerical calculation accuracy was improved by surface correction. The convergence of the model was studied. The influences of the horizon and the nonlocal ratio on the corrected PD model were investigated by measuring the relative difference. We showed that as the horizon decreases, the PD solutions converge to the analytical solution. However, the relative difference of the numerical calculation changed very little with an increasing nonlocal ratio m, while the point spacing is fixed. By comparing two different discretization schemes for PD material points, it was found that when the center of the boundary material points is located on the boundary, the PD solution is more accurate.

For the transient heat conduction of FGMs with cracks, we introduced a method of considering the surface correction of the material points close to an insulated crack to improve the accuracy of numerical calculation near the crack. We then solved the transient heat conduction of FGMs with horizontal static cracks. The results showed better agreement with results from FEM than the results from the uncorrected PD model. The transient heat conduction of the FGM plate with dynamic horizontal crack and dynamic intersecting cracks was simulated. The shielding effect of the thermal conduction of the temperature field caused by the dynamic crack propagation was correctly captured. The developed corrected PD heat conduction model can accurately simulate the transient heat conduction problem of functionally graded materials with cracks.

We plan to study coupling response between the thermal and mechanical of FGMs under thermal shock loading, and simulate the crack initiation, propagation and failure process of FGMs.