Temperature and Displacement Discontinuity Boundary Element Method for Analysis of Cracks in Three-Dimensional Isotropic Thermoelastic Media

Much research has been conducted in the analysis of cracks in thermoelastic materials using analytical or numerical methods [1-4]. Because solving complicated practical problems analytically is difficult, numerical methods are needed. Among the various numerical methods, the boundary element method is very convenient and efficient in such analyses, and many studies have adopted this method to analyze crack problems in thermoelastic media [2-4]. Furthermore, the subsequent displacement discontinuity boundary integral equation and boundary element method are more efficient in studying crack problems, as they grasp the basic characteristic of crack problems; specifically, fields are discontinuous across crack faces. The displacement discontinuity method was first proposed by Crouch [5] to study two-dimensional crack problems, and was then extended to three-dimensional elastic media, piezoelectric media, and magnetoelectroelastic media [6-9], where the electric and magnetic potential discontinuities were introduced across crack faces. Based on previous work, this paper develops the temperature and displacement discontinuity boundary integral equation and boundary element method for which the temperature distribution across the crack faces is assumed discontinuous. The Green functions for unit-point temperature and displacement discontinuities are derived, and the temperature and displacement discontinuity boundary integral equations are obtained for an arbitrarily shaped planar crack. The singular fields ahead of the crack front are discussed, and the stress intensity factors are obtained. For numerical simulations, the Green functions for a triangular element are obtained, and an elliptical crack is analyzed to validate the correctness of the analytical solution and the proposed numerical method. The influence of the temperature and different thermal conditions are also discussed.

Consider an arbitrarily shaped planar crack S located in the plane $z=0$. The upper and lower crack surfaces are denoted $\mathrm{S}^{+}$and $\mathrm{S}^{-}$, respectively. Their respective outer normal vectors of $\mathrm{S}^{+}$and $\mathrm{S}^{-}$are respectively given by

A crack problem can be regarded as a combination of two problems: a no-crack problem under applied loadings and a perturbed problem with loadings only applied on crack faces. The first problem is analyzed to obtain the tractions and heat flux on the crack faces in the perturbed problem, which usually satisfy

where

and the overbar denotes the prescribed boundary conditions.

Because of the presence of the crack, the displacement and temperature across crack surfaces are discontinuous and expressed as

where $\left\|u_i\right\|\left(u_1=u, u_2=v, u_2=w\right)$ and $\|\theta\|$ denote respectively the displacement discontinuity and temperature discontinuity across the crack faces. Using the integral equation method in Zhao et al. [7], the Green functions for unit-point discontinuities are obtained. For simplicity, the final expressions of the Green functions are omitted. Based on the Green functions and the superposition principle, the stresses and heat flux of an arbitrary field point can be obtained in terms of the temperature and displacement discontinuities across crack faces. As the field point approaches the crack face, the integrals become hyper-singular. After using the finite-part integrals and some algebraic mathematical manipulations, one can obtain the boundary integral equations for an insulated arbitrarily shaped planar crack [8]:

where,

Eq. (5) shows that the displacement discontinuity in the normal direction of the planar crack and the temperature discontinuity are decoupled from the in-plane displacement discontinuities in the boundary integral equations, Eqs. (5a)-(5d). Moreover, each is governed by its own equation and related only to its corresponding applied loading. However, the two in-plane displacement discontinuities are coupled with the temperature discontinuity in the other two boundary integral equations. In these singular equations, there are two types of singularities: the hyper-singularity $O\left(1 / r^2\right)$ associated with the two in-plane displacement discontinuities and the strong singularity $O(1 / r)$ associated with the temperature discontinuity.

The method to analyze the singular stress and heat flux fields near the crack front is similar to that in [8]; the details are omitted for brevity. After analysis, it is clear that the temperature and displacement near the crack front have a classical singularity of $r^{1 / 2}$. Consequently, the stresses and heat flux near the crack front are also obtained, after substituting into the definition for the intensity factors, one derives expressions for the stress and heat flux intensity factors in terms of the temperature and displacement discontinuities:

Consider an arbitrarily shaped triangular crack element on the $o x y$ plane in an infinite thermoelastic medium. The three apices of the triangle are denoted $\mathrm{A}\left(x_1, y_1\right), \mathrm{B}\left(x_2, y_2\right)$ and $\mathrm{C}\left(x_3, y_3\right)$, respectively.

On the element, uniformly distributed temperature discontinuity $\|\theta\|^{\mathrm{e}}$ and displacement discontinuities $\|u\|^{\mathrm{e}},\|v\|^{\mathrm{e}},\|w\|^{\mathrm{e}}$ are applied. After integrating the fundamental solutions for unit-point displacement and temperature discontinuities on the triangular element, the associated Green functions are:

Detailed expressions for the influence functions $F_{i j}$ and $F_3$ can be found in Zhao et al. [8].

An elliptical crack is studied for simplicity. It is first discretized into $N$ triangular elements (Figure 1). On each element the displacement and temperature discontinuities are assumed constant for convenience. Using the Green functions for a constant triangular element obtained above and the superposition principle, one can obtain the following algebraic equations

where $i=1,2,3, \ldots, N$. There are a total of $4 N$ equations and corresponding $4 N$ unknown discontinuities for the whole crack. After solving these linear algebraic equations, one can obtain the discontinuities for each triangular element. The displacements, stresses, temperature, and heat flux at any point can all be determined.

As for the determination of the stress intensity factors, three points in front of the crack tip are used to fit the approximate expressions,

where, $\chi$ with different subscripts denotes the fitting coefficients, and $r$ is the distance from the point to the crack tip. Substituting Eq. (10) into Eq. (7), the intensity factors are thereby determined.

The integral equation in Eq. (5d) is for an insulated crack. However, if the effect of the medium inside the crack is to be more realistic, the thermal semi-permeable boundary condition must be adopted, and Eq. (5d) becomes

with

where $\beta^c$ denotes the thermal conductivity of the medium inside the crack. Similar to Fan et al. [9], an iterative approach is proposed to determine the heat flux at an arbitrary point in the crack.

First, the crack model is assumed insulated; specifically $h_z^c=0$. After using the boundary element method, we obtain the corresponding crack opening displacement discontinuity $\|w\|^{(1)}(r)$ and temperature discontinuity $\|\theta\|^{(1)}(r)$ at an arbitrary point. Second, by substituting the value of the discontinuities into the determining Eq. (12), we obtain a new heat flux $h_z^{c(1)}(r)$ in the crack cavity,

Inserting the new $h_z^{c(1)}(r)$ into Eq. (12) and beginning a new calculation, one obtains new discontinuities. As was mentioned earlier, the convergence criterion for the iteration to stop is that the value of the heat flux in the crack cavity remains the same (specifically, the difference in the two new heat fluxes is sufficiently small). In numerical simulations, if the convergence criterion satisfies a preset precision threshold $\delta$, the iteration is stopped, i.e.,

The value of the heat flux in the crack cavity is the final result, $h_z^c(r)=h_z^{c(k)}(r)$.

In this section, the proposed method is first checked. The thermoelastic medium aluminum is chosen for the simulation, and its material properties being

Uniform loadings are applied to the crack face

If the elliptical crack is assumed to be insulated and is reduced to penny-shaped, the numerical results can be compared with the analytical solutions given by Florence and Goodier [1]. Figures 2 and 3 plot the distribution of the displacement in the radial direction and temperature discontinuities, respectively. The numerical results are in perfect agreement with the analytical solution. As for the stress intensity factor $K_{\text {II }}$, the relative error is only $0.45 \%$ between the numerical results and the analytical solution.

In this section, the influence of the thermal boundary condition is studied. The final value of the heat flux in the crack cavity should be determined at first. Figure 4 plots the variation of the heat flux in the cavity versus the iteration number. The value of the heat flux clearly converges and remains constant after 30 steps in the iteration.

Figures 5 and 6 plot the elliptical ratio dependence of the normalized mode II stress intensity factor and heat flux intensity factor, respectively, for the major and minor axes under different thermal conditions. The value of the intensity factors at the major axis tip are smaller than those at the minor axis tip; this indicates that the elliptical crack is more inclined to extend and grow circularly in the minor axis direction with the facture occurring more easily at the minor axis tip. The value of the intensity factors under thermal semi-permeable condition is smaller than those under thermal insulated conditions, regardless whether it is the major or minor axis tip. In summary, the intensity factors under thermal semi-permeable conditions are smaller than those under thermal insulated conditions. As thermal semi-permeable conditions are more realistic, engineering designs based on fracture criteria for thermally insulated condition tend to be conservative.

The displacement discontinuity method was adopted and extended by introducing the temperature discontinuity to analyze planar cracks in three-dimensional isotropic thermoelastic media. The Green functions for unit-point discontinuities are obtained and employed to construct the temperature and displacement discontinuity boundary integral equations for an arbitrarily shaped planar crack. The singularities of the fields near the crack front are analyzed, and the stress and heat flux intensity factors are expressed in terms of temperature and displacement discontinuities across crack faces. In numerical simulations, the Green functions for a constant triangular element are derived, based on which the temperature and displacement discontinuity boundary element method was proposed. An iterative approach is also proposed to determine the value of the heat flux in the crack cavity for thermal semi-permeable conditions. From the results, it is concluded that the proposed method is in perfect agreement with the analytical solution, and that the influence of the thermal conditions cannot be ignored.