Study on thermal fracture modeling by the Scaled boundary finite element polygons

Thermal load is an important cause of structural fracture. A model based on the Scaled boundary finite element method (SBFEM) is presented for the fracture modeling due to thermal load. The SBFEM is a semi-analytical method, which has the advantages of both the Finite element method (FEM) and the Boundary element method (BEM) while possessing its own special features. The structure is discretized only at the boundary, thus reducing the dimension of solution. In addition, the stress field and displacement field are obtained analytically in the radial direction, so the stress intensity factors can be obtained directly from its definition and this works for homogeneous material, bi-material and orthotropic material. In the presented model, the temperature load is transformed from the temperature field of the structure, which is included in the structural analysis to get the displacement field, stress field and therefore the stress intensity factors. By introducing the polygon technique, the remeshing process due to crack propagation can be modeled conveniently. Two numerical examples including a rectangular plate and a cross plate are employed for verification of the model.


Introduction
For high concrete dams, the fracture problem due to temperature load is a great concern. The heat of hydration, air temperature, casting temperature, solar radiation all contributes to temperature change in the dam body, inducing thermal stress which may exceed the tensile strength of concrete. In this situation, crack may occur at the surface of or inside the dam, which does harm to the integrity and durability of the dam [1]. Researchers have performed quite some research on this issue and made progress in the calculation of temperature field and fracture modeling. Wilson developed the code DOT-DICE based on the FEM for two-dimensional analysis of mass concrete structures [2]. Subsequently the engineers from the United States Army Corps of Engineers made some revisions to the code and applied it to the analysis of the temperature field of the Willow Creak dam [3]. The temperature data for different time period was obtained and accorded with the data measured very well. Cheng and Chang et al. applied the Meshless Galerkin method to compute the stress intensity factors and model the crack propagation of hydraulic structures [4]. No remeshing was needed during the crack propagation and the distribution of nodes could be adjusted freely, the crack path could be traced dynamically and was not influenced by the element boundary. Liu and Zhou et al presented a model based on the Meshless method for modeling of crack propagation of concrete block due to temperature stress. The results were found to be consistent with crack distribution of mass concrete block due to uniform temperature drop and strong constraint at the base [5] The study on the thermal fracture includes two aspects, i.e., the modeling of the temperature field and the crack propagation due to the thermal load. The Finite element method falls short of the latter due to poor precision at the crack tip and great computational effort involved in remeshing due to crack propagation. Yet the SBFEM has great advantages for this. It is a semi-analytical method, which has the advantages of both the Finite element method (FEM) and the Boundary element method (BEM). The structure is discretized only at the boundary as in the BEM, but as compared to BEM no fundamental solution is needed [6]. When applied to fracture problems, the stress and displacement in the radial direction is obtained analytically, and the stress intensity factor is obtained directly from its definition with high precision. Song developed the definition and computation method for thermal load, by which the temperature field, displacement field and stress field can be obtained with high precision [7]. Concerning the modeling of crack propagation, the polygon technique is becoming more and more popular these years, which can reduce the effort involved in remeshing to the minimum [8][9][10][11].
In the presented research, a procedure based on the Scaled boundary finite element polygons is developed for the computation of thermal fracture process modeling of structures. By discretizing the structure with polygons and solved by the SBFEM for both thermal analysis and structural analysis, the fracture process of the structure due to thermal load can be modeled and two examples are investigated for verification.

Solution of temperature field
The SBFEM (the Scaled Boundary Finite Element Method) is employed for thermal analysis and structural analysis. The formulation of the fundamental equations for the SBFEM is referred to the literature [7] and only the main equations are presented here. The fundamental equation of temperature T() with  stands for the local radial coordinate for two-dimensional case can be derived by the Principle of virtual work or the Matrix function method [13] as follows.
Here [E 0t ], [E 1t ] and [E 2t ] are the coefficient matrices of the polygons at the boundary as follows. The superscript t denotes temperature in order to distinguish from the matrices in structural analysis. The coefficient matrices are decided only by the geometry and material properties of the polygons. The matrices of each polygon are then assembled as in standard Finite element method.
Here k denotes the coefficient of thermal conductivity. The inner heat source in the radial direction is The following first-order ordinary differential equations are obtained by transforming Eqs. (1) and (5): in which [Z] is a coefficient matrix which is a Hamiltonian matrix as follows.
By performing the block-diagonalization Schur decomposition of the matrix [Z], the following relation holds.
Then by substituting Eq. (8) into (6), the temperature solution can be obtained as follows.
Here c t is the integration constant.
By interpolating T() by the shape function N()，the temperature field T(,) of each polygon can be arrived at.

Solution of displacement field
The displacement can be obtained analytically in the radial direction by the SBFEM. The displacement in the radial direction is denoted as {u()}with  stands for the local radial coordinate. Based on that, for a point with local coordinate (,) where  is the circumferential coordinate, its displacement is obtained by interpolation through the shape functions in the circumferential direction as in standard FEM.
And the stress field in which [D] is the elasticity matrix of the material, B 1 () and B 2 () are the matrices relating the strain with the displacement [12]. The fundamental equation of the SBFEM with the displacement function as the basic unknown is [15][16]: Here [E 0 ], [E 1 ] and [E 2 ] are the coefficient matrices of the polygons at the boundary, which are decided only by the geometry and material properties of the polygons. The matrices of each polygon are then assembled as in conventional Finite element method. {F(ξ)} is the external load vector. If the external load involves only the temperature load, {F(ξ)} can be written as follows [14]: where d is the opposite of the eigenvalues when computing the temperature field and () is the determinant of the Jacobi matrix. The internal nodal force is  (19) in which {q0()}is the internal force induced by the initial temperature stress.
The initial stress {σ0} induced by temperature can be obtained as follows.
Here T is the temperature solved from Eq. (11) and {β} is a coefficient vector related to the material properties. For plane strain problem, And for plain stress problem, in which ν is the Poisson's ratio and α is the coefficient of thermal expansion. By solving Eq. (14), the displacement can be obtained.

Computation of stress intensity factors and crack propagation angle
A polygon containing a crack is shown in Figure 1. The displacement and the stress intensity factors in the local coordinate system has the following relation [17].  Figure 1. A polygon containing a crack tip Figure 2. Plate with a crack For calculation of the stress intensity factors using the SBFEM, only the terms related to singular stress are involved, which leads to high precision of the stress intensity factors and with less computational efforts. The crack propagation angle C  is decided by the maximum circumferential stress criteria. By substituting the angle into the following equation   2 3 cos sin cos 2 2 2 The equivalent stress intensity factor   C K  can be obtained. And if it is equivalent to the fracture toughness KIC of the material, the crack propagates.

Numerical examples
A rectangular plate and a cross plate are investigated for verification of the model presented.

A central-cracked isotropic isothermal plate
Considering a central-cracked isotropic isothermal plate as shown in Figure 2. Its dimension is 2W×2W, the length of the central crack is 2a. The modulus of elasticity E=2.184x10 5 Pa, Poisson's ratio ν=0.3, and the coefficient of thermal expansion α=1.67×10 -5 /℃. The boundary condition of the plate is also given in the figure, where T denotes the temperature and h q means adiabatic condition. In order to compare with results in literature, plain strain condition is assumed. Considering the symmetry of the structure and the load, one half of the model is analyzed and discretized with 45 polygons (Figure 3). The contour of temperature distribution is shown in Figure 4. As it can be seen that the temperature increases from the crack to the boundary of the plate, which agrees well with the result obtained by the XFEM. Shown in Figure 5 is the contour of displacement. The horizontal displacement increases from the left to the two corners at the right side. The vertical displacement reaches the maximum and minimum respectively at the upper left corner and lower left corner of the plate.
The normalized stress intensity factors K * =K/Eα△TW 1/2 are listed in Table 1. As it can be seen that the current results are consistent with those obtained by XFEG [18] and FFEM [19]. With the increase of the crack length, the stress intensity factor increases gradually.

A cross plate
Consider a cross plate with a crack as shown in Figure 6. The contour of temperature distribution is shown in Figure 9. As it can be seen that the temperature increases from the bottom of the crack to the top and the temperature contour agrees quite well with that obtained by XEFG [18]. Shown in Figure 10   The crack path is shown in Figure 8 and it can be seen that the crack path agrees well with that obtained by the Double boundary element method (DBEM). The crack propagates towards the upper right direction until the plate fails, which is consistent with the actual failure mode of a cross plate.
The normalized stress intensity factors are also compared with those derived by the DBEM [20] as shown in Figure 11. It can be seen that the two series of results are quite similar. With the growth of the crack, K I * increases and then decreases continuously, K II * decreases sharply first and then fluctuates until decreasing gradually.

Conclusion
A procedure is presented for the thermal fracture modeling based on the Scaled boundary finite element polygons. The fundamental equations for both the thermal analysis and structural analysis based on the SBFEM are brought forward. After solution of the temperature field, the thermal load is approximated by the function {F(ξ)}, which is then substituted into the fundamental equation of structural analysis to obtain the displacement and stress field. The SBFEM enables the computation of the stress intensity factors directly obtained from its definition with high precision. And the polygon technique greatly facilitates the remeshing process and reduces the computational efforts involved in remeshing to the minimum. The procedure is verified by two examples of plates.