# Effect of T-stress on Fracture of a Long Cracked Plate in Unsteady Heat Transfer

Authors: Shaoqing Zhou; Limin Li
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

In order to accurately predict the structure fracture caused by thermal load, a modified maximum tensile stress (MTS) criterion combined withT-stress is proposed. The modified MTS uses a two-parameter model (stress intensity factor K and T-stress) to describe the fracture behavior under thermal load. The T-stress and stress intensity factor at the crack tip are solved by using J-integral in the theoretical calculation of a cracked strip with temperature difference. The results show that T-stress can affect the fracture toughness and the stress at the crack tip of the cracked strip with temperature difference. This provides a basis for the simulation of structural fracture under thermal load.


---

## Body

## 1. Introduction

In practical engineering, thermal stress analysis is an essential aspect because some structures are often affected by temperature changes. There are many engineering structure fractures due to thermal load. In recent years, the international rock mechanics field has considered temperature as an important factor affecting rock mechanical properties [1, 2]. Therefore, some research has been done on the influence of temperature on rock strength and fracture toughness [3, 4]. Dölling et al. [5] predicted thermally induced edge-crack initiation using finite fracture mechanics. Yang and Wang [6] studied the fracture mechanism of crack under the interaction of constraints. Because the variation of J-resistance curve is complex, it is necessary to consider the interaction of constraints. Qian et al. [7] analyzed the effect of temperature on the model parameters in local approaches to cleavage fracture. Kim et al. [8] analyzed the effect of rapid thermal cooling on mechanical rock properties. Liu and Xu [9] carried out the dynamic mechanical experiments on marble under different temperatures and different strain rates by using the high-temperature split Hopkinson pressure bar (SHPB) experimental system. In addition, Le and Gardin [10] discussed the stress intensity factor of the cracked body for different crack lengths obtained by a closed-form integration of the stress field. The traditional theory usually uses a single parameter to describe the stress and strain fields at the crack tip. That is, the stress function only considers the singular part of Williams function (stress intensity factor K) but ignores the nonsingular part (T-stress). Tang et al. [11] studied the effect of T-stress on crack growth path in rock and fracture strength. The traditional theory usually uses a single parameter to describe the stress and strain fields at the crack tip. Therefore, there are few methods to study thermal stress fracture with T-stress.To study the fracture of similar rock or pavement under the action of temperature difference, combined withT-stress, the fracture of an infinite plate with cracks under unsteady heat transfer is studied theoretically.


## 2. Theoretical Model of Thermal Stress Fracture Calculation

As shown in Figure1, the height of the strip is L, the initial temperature inside the plate is T0, the insulation bottom layer is set to y = 0, the height of the upper surface is y = L, and the crack length at the upper surface is a. When the time is t = 0, the surface is suddenly heated or cooled to Ta.Figure 1
Initial temperatureT0 and external temperature Ta of crack strip.For convection heat transfer problems, it is convenient to use the Biot number to calculate the results. The Biot number is expressed as follows:(1)β=hLκ,where h is the heat transfer coefficient, L is the thickness of the plate, and κ is the thermal conductivity of the plate.Using the linear elasticity of materials and the superposition principle, the problem is divided into two parts. As shown in Figure2, the stress produced by thermal shock in the long elastic plate without crack is calculated in Figure 2(a). Then, the effect of cracks in the long plate is considered in Figures 2(b)–4, and the thermal elastic stress in the opposite direction is taken according to the stress results of the long elastic plate without crack.Figure 2
Calculation of crack by superposition method: (a) stress distribution of long plate without crack; (b) phase reaction force of cracked plate at crack location.As shown in Figure2, the stress field of the strip with crack is obtained.


## 3. Finite Element Calculation Theory of Heat Conduction of Long Elastic Plate without Crack

The nondimensional equations governing the structural response are the plain-strain equations(2)∂σx∂x+∂σxy∂y=0,(3)∂σxy∂x+∂σy∂y=0.The heat conduction differential equation that the field variableϕx,y,t of the transient temperature field should satisfy in rectangular coordinates is as follows:(4)ρc∂θ∂t−D∂2θ∂x2+∂2θ∂y2−ρQ=0withinΩ,where θ is the temperature field, D is the thermal diffusivity, t is the time, Q is the heat source density inside the object, c is the specific heat, and ρ is the density.The boundary conditions that the field variables should satisfy at the same time are as follows:(5)θ=θ0S,ton the boundary ofS1,(6)D∂θ∂xnx+∂θ∂yny=qS,ton the boundary ofS2,(7)D∂θ∂xnx+∂θ∂yny=hθa−θon the boundary ofS3,where nx¸ ny refer to the direction cosine of the normal outside the boundary; θ0S,t is the given temperature on the S1 boundary; qS,t is the given heat flow on the S2 boundary; h is the heating coefficient; and θaS,t is the ambient temperature.The boundary should meet the following requirements:(8)S1+S2+S3=S,where S is the entire boundary of Ω domain.When timet = 0, the initial conditions are as follows:(9)θx,y,t=0=θ0x,y.The corresponding variational formulation is as follows: in the allowable temperature field satisfying the boundary conditionS1, S2, S3 and initial condition equation (7), the real temperature field should make the following function take the minimum value, i.e.,(10)Iθmin=12∫Ωα∂θ∂x2+∂θ∂y2+2ρQ−ρc∂θ∂tdΩ.For the transient heat transfer problem, the temperature field of the unit will change with time:(11)θex,y,t=Nx,y⋅qθet.Here the node temperatureqθet changes with time, i.e.,(12)qθet=θ1t,θ2t⋯θntT.Substituting equation (12) into equation (11) and finding the variational extremum for qθet, we get(13)Cθeq˙θe+Kθeqθe=Pθe.Equation (13) is a set of linear ordinary differential equations with time t as the independent variable, where Cθe is the heat capacity matrix, Kθe is the heat conduction matrix, Pθe is the temperature load matrix, qθe is the node temperature matrix, and q˙θe is the derivative array of node temperature to time, q˙θe=dqθe/dt=dθ1t/dt,dθ2t/dt,…,dθnt/dtT.The basic concept of using the numerical integration method to solve ordinary differential equations is to discretize the time domain. The backward difference scheme is adopted:(14)Δθn=θn+1−θn=Δtn∂θ∂tn+1.Substituting equation (14) into equation (13),(15)Kθe+1ΔtnCθe⋅θn+1−1ΔtnCθe⋅θn=Pθn+1e,n=0,1,2,…,where Kθe is the heat conduction matrix, Δtn is the time step, Cθe is the heat capacity matrix, θn is the temperature distribution, and Pθn+1e is the temperature load matrix.In equation (15), θn, Pθn+1e are known and θn+1 is unknown. The temperature distribution θn+1 at time t=tn+1 can be obtained by the linear equations. Finally, the instantaneous node temperature can be obtained by recursion.Because of thermal expansion, only linear strain is produced, but the shear strain is zero. The strain caused by thermal deformation can be regarded as the initial strain of the object. In the calculation of thermal stress, the initial strainε0 caused by thermal deformation is calculated, and the equivalent nodal load Pε0 caused by the corresponding initial strain is obtained. Then, the node displacement S caused by thermal deformation can be solved as usual, and then the thermal stress σ can be obtained from S.Thermal expansion will be caused by the thermal load; then, the physical equation of the object becomes(16)σij=Dijklεkl−εkl0.According to the general virtual work principle of elasticity, its functional expression is as follows:(17)Hpu=∫ΩσijδεijdΩ−∫ΩfδuidΩ+∫SppδuidA.Substituting equation (15) into equation (16), we get that(18)Hpu=∫ΩDijklεkl−εkl0δεijdΩ−∫ΩfδuidΩ+∫SppδuidA.Equation (18) is the virtual work principle of thermal stress problem.The solution domainΩ is discretized by finite element method, and the finite element solution equation can be obtained from δHp=0:(19)Keqe=P,where P=Pf+Pp+Pθ, in which Pf and Pp are the load terms caused by volume load and surface load and Pθ is the load term caused by temperature strain, i.e., Pθ=∑e∫ΩeBTDε0dΩ..It can be seen from the above formulas that compared with the stress analysis problem without thermal load, the structure thermal stress problem is the same except for adding a temperature load.


## 4. Finite Element Calculation Theory of J-Integral and T-Stress for Cracked Elastic Long Plates

First,J-integral is calculated, and then T-stress is calculated by J-integral:(20)J=∫Γ0Wn1−Fiui,1dΓ+∬Ωσiiαθ,1dΩ,where θ is the temperature value and α is the coefficient of thermal expansion.According to Duhamel–Neumann relation, for an isotropic linear elastic body, the following stress expression can be obtained:(21)σkk=3λ+2μuk,k−3αθ,where λ and μ are Lamé constants.Using equation (21), the strain energy density is obtained:(22)W=∫0εijσijdεij=μεijεij+12λεkk2−2μ+3λαθεkk+322μ+3λα2θ2.By using Gauss divergence theorem for the region integral of the second formula on the right in equation (20), it is obtained that(23)∬Ωσiiαθ,kdΩ=α3λ+2μ∫Γui,iθnkdΓ−∬Ωui,ikθdΩ−3α23λ+2μ∬Ωθθ,kdΩ.By substituting equations (22) and (23) into equation (20), the final J-integral expression is obtained:(24)J=∫Γ0μεijεij+12λεkk2−Fiui,1−3λ+2μ2λ+μα2θ2n1+μ3λ+2μλ+μαθμ1,i−θ,iu1nidΓ.To calculate the T-stress, the paper uses the auxiliary solution, namely, a unit load is set at the crack tip and the load is parallel to the crack surface, as shown in Figure3.Figure 3
Integration path and definition coordinates.Consider two independent equilibrium statesσijA,εijA,uijA,θA and σijB,εijB,uijB. The second case is the auxiliary solution, which only considers the case of unit force f at the crack tip.We useJA+B to express the J-integral after superposition of two cases (i.e., A and B). According to the theory of small displacement:(25)JA+B=∫Γ012σijA+σijBεijA+εijBn1−σijA+σijBnjui,1A+ui,1BdΓ+∬ΩσiiA+σiiBαθ1AdΩ.The following integrals can be obtained, respectively:(26)JA=∫Γ012σijAεijAn1−σijAnjui,1AdΓ+∬ΩσiiAαθ1AdΩ,(27)JB=∫Γ012σijBεijBn1−σijBnjui,1BdΓ.IntroduceM, and its expression is as follows:(28)M=JA+B−JA−JB.M integral can be obtained:(29)M=∫Γ012σijAεijB+σijBεijAn1−σijAnjui,1B−σijBnjui,1AdΓ+∬ΩσiiBαθ1AdΩ.The integral path is optional, so the radius ofΓε can be reduced to zero. In this special case, the singular term will not affect the integral. Thus, only the nonsingular term (T-stress) contributes to this integral, so equation (29) becomes(30)M=limε⟶0∫ΓεσijTεijBn1−σijTui,1Bnj−σijBnjui,1TdΓ,where σijT=Tδi1δj1.Under plane strain condition:(31)εijT=1+υEδi1δj1−υδijT+1+υαθδij,εi3T=ε33T=0.Therefore,(32)ui,1T=u1,1Tδi1=ε11Tδi1.Combining (30) with (32), the M integral is obtained:(33)M=−limε⟶0∫Γεσ1jBnjε11TdΓ=1−ν2ET⋅f+αθf⋅1+ν.Then, T-stress is calculated by equation (33).(34)T=E1−ν2fM−αθf⋅1+ν.


## 5. Modified MTS Criterion for Mixed-Mode Fracture

The elastic T-stress represents the constant tensile stress acting parallel to the crack flanks, which is related to the second term (the first nonsingular term) in the Williams expansion of the stress field:(35)σθθ=122πrcosθ2KΙ1+cosθ−3KΙΙsinθ+Tsin2θ+or1/2,where KI, KII and T denote the mixed-mode SIF and T-stress, respectively.In the present paper, a two-parameter modification of the maximum tensile stress (MTS) criterion is applied to determine the stress intensity factor and theθ value. The MTS criterion is one of the most widely used theories for mixed mode I and mixed mode II crack growth. It states that a crack propagates in the direction for which the tangential stress is maximum. The two-parameter modification of the MTS criterion takes into account the existence of the second term in the Williams expansion of the stress field and has the following form:(36)KΙcθ=cosθc2KΙcos2θc2−32KΙΙsinθc+2πrcTsin2θc.


## 6. Finite Element Simulation Analysis

Let the elastic modulusE=2×105MPa, Poisson’s ratio ν=0.3, and the coefficient of thermal expansion α=1.67×10−5deg−1 of the long plate. The nondimensional stress intensity factor is K∗=K1−ν/αEΛθa, where Δθ is the temperature difference, a is the crack length, t∗ is nondimensional time parameter, t∗=Dt/L2, and D is the thermal diffusivity. The dimensionless T-stress is set to T∗=T/αEΔθ.Figure4 shows that the stress intensity factor also changes with the crack length changing with the fixed Biot value β=10. For a given crack length ratio a/L, the dimensionless stress intensity factor K∗ increases first and then decreases after reaching its max value. It can be seen from Figure 4 that for very small crack length ratio a/L, the stress intensity factor is relatively much larger than that for long cracks. With the increase of crack length a/L, the stress intensity factor decreases.Figure 4
Variation of stress intensity factor with different crack lengths (K∗=K/K0,t∗=Dt/L2,β=10).Figure5 compares the dimensionless stress intensity factor changes of different Biot values under the condition of constant crack length a/L=0.1. It can be seen from Figure 5 that in thermal shock, not only the Biot value can affect the maximum stress intensity factor but also the maximum stress intensity factor is related to the dimensionless time t∗. Therefore, when the Biot number is reduced, its effect on the stress intensity factor is similar to that of increasing the crack length.Figure 5
Variation of stress intensity factor with different Biot numbers (K∗=K/K0,t∗=Dt/L2,a/L=0.1).Figure6 shows the relationship between dimensionless T-stress and crack length. For a given crack length ratio a/L, as the dimensionless time function t∗ increases, the dimensionless T-stress increases first and then decreases gradually after reaching its maximum value. T-stress decreases with increasing crack length. The time required for T-stress to reach its maximum is different for different crack lengths.Figure 6
T-stress variations for different crack lengths (T∗=T/αET0−Ta,t∗=Dt/L2,β=10).Figure7 compares the dimensionless T-stress changes of different Biot values with the same crack length a/L=0.1. It can be seen from Figure 7 that not only the Biot value can affect the maximum T-stress but also the maximum T-stress is related to the dimensionless time t∗. Therefore, combined with Figures 6 and 7, T-stress at the crack tip is related to the Biot number, crack length, and time.Figure 7
T-stress variations with different Biot numbers (T∗=T/αET0−Ta,t∗=Dt/L2,a/L=0.1).


## 7. Comparison between the Experimental Results and the FEM results

Kokini and Long [12] have studied the fracture of PMMA under thermal loading. They heated the specimen first and then tried to cool one end of the specimen with cold air. The material parameters are shown in Table 1.Table 1
Thermal parameters of PMMA materials.
Elastic modulus (GPa)Specific heat capacity (J/kg⋅°C)Density (kg/m3)Thermal conductivity (W/m⋅°C)Coefficient of thermal expansion (°C)Poisson ratio2.391104.611690.3651.46 × 10−40.38As can be seen from Figure8, because T-stress is considered, the predicted value of fracture toughness is consistent with the experimental results. For mixed-mode loading, the effect of singular stress is represented by an effective stress intensity factor Keff=KI2+KII2. Zhou and Guo [13] showed that the mixed-mode fracture toughness, i.e., Keff, depends on the magnitude and the sign of T-stress. In general, a negative T-stress increases the mixed mode I/II fracture toughness of cracked specimens, and conversely the mixed-mode fracture toughness decreases for specimens having a positive T-stress. Zhou and Guo [13] indicated that the T-stress term in Williams series expansion has a great influence on the calculation results of stress intensity factor, and the calculation results of stress intensity factor considering the influence of T-stress are more consistent with the experimental results.Figure 8
Experimental results and finite element simulation predictions.


## 8. Conclusion

(1)
In a cracked strip with a temperature difference, if the Biot value is fixed, the stress intensity factor changes with the length of the crack. For a given crack length ratio, as the dimensionless time function increases, the dimensionless stress intensity factor increases first and then decreases after reaching its maximum value. For very small crack length ratios, the stress intensity factor is relatively much larger than that for long cracks. The stress intensity factor decreases with increasing crack length.(2)
In a cracked strip with temperature field difference, not only the Biot value can affect the maximum stress intensity factor but also the maximum stress intensity factor is related to the dimensionless time. Therefore, when reducing the Biot number, its effect on the stress intensity factor is similar to that of increasing the crack length.(3)
For a given crack length ratio, the T-stress increases first and then decreases gradually as the time function increases. Different crack lengths take different times for the T-stress to reach its maximum. The shorter the crack, the shorter the time required. T-stress at the crack tip is related to Biot number, crack length, and time.



---

*Source: 2700064-2021-11-13.xml*