Thermal shock analysis of 2D cracked solids using the numerical manifold method and precise time integration

Engineering Analysis with Boundary Elements 75 (2017) 46–56

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

Thermal shock analysis of 2D cracked solids using the numerical manifold method and precise time integration ⁎

H.H. Zhanga, , G.W. Mab,c, L.F. Fanc, a b c

MARK

⁎

School of Civil Engineering and Architecture, Nanchang Hangkong University, Nanchang, Jiangxi 330063, China School of Civil, Environmental and Mining Engineering, The University of Western Australia, 35 Stirling Highway, Crawley, WA 6009, Australia College of Architecture and Civil Engineering, Beijing University of Technology, Beijing 100084, China

A R T I C L E I N F O

A BS T RAC T

Keywords: Crack Thermal shock Numerical manifold method Precise time integration method Transient temperature field Stress intensity factors

The numerical manifold method (NMM), combined with the precise time integration method (PTIM), is proposed for thermal shock fracture analysis. The temperature and displacement discontinuity across crack faces is naturally portrayed attributing to the cover systems in the NMM. The crack tip singularities are characterized through the use of asymptotic bases in the approximations. The discrete equations for transient thermal analysis are firstly solved with the PTIM and then the thermoelastic study is performed. With the interaction integral, the stress intensity factors are computed. Several examples are tested and the nice consistency between the present and existing results is found.

1. Introduction Engineering equipment such as the aero-engines, gas turbines and pressure vessels is frequently exposed to thermal shock. Under certain conditions, thermally induced stresses may cause the fracture and failure of cracked facilities. Consequently, the study of the behavior of cracked solids under thermal shock is of great importance. In light of the significance, a lot of work has been carried out during the past several decades. Many researchers focused on the analytical solutions. Oliveira and Wu [1] determined the thermal stress intensity factors (TSIFs) of axial cracks in hollow cylinders under thermal shock using the closed form weight function formula. Noda and Ashida [2] adopted the successive approximation, Fourier integral and Bessel series to solve transient thermal annular crack problem in an infinite transversely isotropic cylinder. Lee and Kim [3] calculated the TSIFs of elliptical surface cracks in thin-walled and thick-walled cylinders with the modified Vainshtok's weight function method. Noda and Wang [4] applied the approach of singular integration equation to tackle the thermal shock responses of the functionally graded materials (FGMs) with collinear cracks. Shahani and Nabavi [5] used the finite Hankel transform and the weight function method to compute the TSIFs of internal semi-elliptical crack in a thick-walled cylinder subjected to transient thermal stresses. Wang and Li [6] analyzed the transient thermoelastic fracture of piezoelectric material with the Laplace transformation and integral equation method. Although widely applied, analytical solutions are generally limited

⁎

to simple configurations (e.g., single or periodic cracks in infinite or semi-infinite mediums under regular initial or boundary conditions). For problems with finite dimensions, arbitrary cracks or under complex loadings, numerical tools such as the finite element method (FEM), the boundary element method (BEM) and the extended finite element method (XFEM) are much more popular. Emmel and Stamm [7] computed the transient TSIFs for cracked rectangular plate and hollow cylinder with the FEM. Magalhaes and Emery [8] inspected the effects of transient thermal loads on crack propagation in brittle film-substrate structure by the FEM. Through the three-dimensional elastic-plastic FEM, Kim et al. [9] evaluated the integrity of vessel with subclad crack under pressurized thermal shock. Prasad et al. [10] applied the dual BEM and path-independent J-integral to investigate the two-dimensional (2D) transient thermoelastic crack problems. Considering crack closure conditions, Giannopoulos and Anifantis [11] obtained both steady and transient TSIFs of 2D bimaterial interfacial cracks by the BEM. With the uncoupled thermoelastic theory, Zamani and Eslami [12] implemented the XFEM to model the effect of both mechanical and thermal shocks on 2D cracked solids. Rokhi and Shariati [13] studied the response of cracked FGMs under thermal shock with the XFEM in the framework of coupled thermoelasticity. In the past two decades, considerable efforts have been put on to the development of the numerical manifold method (NMM) proposed by Shi [14]. The soul of the NMM lies in the use of finite cover concept, which has been adopted by Bathe and his coauthors [15–17] to improve the FEM solutions most recently. Benefiting from the use of

Corresponding authors. E-mail addresses: [email protected] (H.H. Zhang), [email protected] (L.F. Fan).

http://dx.doi.org/10.1016/j.enganabound.2016.11.012 Received 26 September 2016; Received in revised form 29 November 2016; Accepted 30 November 2016 0955-7997/ © 2016 Published by Elsevier Ltd.

Engineering Analysis with Boundary Elements 75 (2017) 46–56

H.H. Zhang et al.

dual cover systems (i.e., the mathematical cover and the physical cover) [14,18], the NMM is very powerful in discontinuous analysis (e.g., in solving crack or inclusion problems). The major highlights of the NMM for crack modeling can be summarized in four aspects: (1) the mathematical cover can be independent of all domain boundaries including cracks; (2) the discontinuity of physical field across crack faces can be manifested in essence; (3) the local property at crack tip zone can be well captured through the use of associated local functions in the approximation, and (4) higher-order approximations can be achieved through the use of higher-order local functions on a fixed mathematical cover. To date, the NMM has been successfully improved to solve various stationary cracks and crack growth problems in homogeneous [19–30] and heterogeneous materials [31–35]. In the present paper, the NMM is further extended to tackle 2D stationary crack problems under thermal shock. The solution procedure is generally divided into two parts: Firstly, the transient heat diffusion problem is analyzed and the corresponding NMM discrete equations are solved with the precise time integration method (PTIM) [36], which has absolute stability, immunity to oscillations and timestep-independent precision. Subsequently, the calculated temperature fields at selected instants are imported into the thermoelastic part to extract the TSIFs. The remaining of the paper is addressed as follows. In Section 2, the governing equations and associated boundary and/or initial conditions are provided. The NMM formulations for both transient heat conduction and thermoelastic analysis are derived in Section 3. Section 4 presents the details of the PTIM for transient heat conduction analysis and the spatial integral scheme as well. To verify the proposed method, several numerical examples are tested in Section 5. Finally, the concluding remarks are addressed in Section 6.

q = − k ∇T with k the thermal conductivity. The stress tensor σ is expressed by the generalized Hooke's law as σ = C: (ε − εT ) with C the fourth-order Hooke tensor, ε = ∇s u the total strain tensor and εT = α(T − TR )I the thermal strain tensor. u is the displacement vector, α is the thermal conductivity and TR is the reference temperature.∇,∇s and I are, respectively, the gradient operator, the symmetric gradient operator and the second order identity tensor. The associated boundary conditions for Eqs. (1) and (2) are

2. Statement of problems

3.1. NMM approximations

As shown in Fig. 1, consider a cracked isotropic homogenous physical domain Ω enclosed by the boundary Γ in the unidirectional coupling transient linear thermoelasticity (i.e., the thermal loading affects the displacement, strain and stress fields, but not vice versa). Ignoring both the heat source and the body force, the governing equations for this problem are [10]

In the NMM, to solve a crack problem, the mathematical cover (MC), composed of a series of mathematical patches (MPs), is firstly built. Broadly speaking, the MP, formed by mathematical elements, can be of any shape and overlapped. The MC may be independent of all domain boundary (including cracks) but must be large enough to cover the whole cracked domain. On each MP, a partition of unity (PU) [37] weight function is defined. Next, the physical patches (PPs), the collection of which is termed as the physical cover (PC), are formed by the intersection of MPs and physical domain. On each PP, the local function is constructed to represent the local physical property. Then, the manifold elements (MEs) are generated through the shared region of PPs. To make the above concepts and procedures clear, an illustrated example in Fig. 2 is further adopted. Consider the physical domain Ω = Ω1 ∪ Ω2 in Fig. 2a. The MC consisting of two rectangular MPs, i.e., M1 and M2 in Fig. 2b, is chosen to cover the whole domain. The intersection of the physical domain and the MPs gives 4 PPs in Fig. 2c, that is, P1 and P2 from M1, and P3 and P4 from M2. Above 4 PPs finally generate 7 MEs, numbered e1, e2…, and e7 in Fig. 2d, where the bracketed contents represent the associated PPs. Following the aforementioned process, we can obtain the NMM approximation on each ME by pasting the local functions using the associated weight functions. Accordingly, for the present problem, the temperature and displacement field on any ME e is approximately expressed as

ρc

∂T (x , t ) + ∇q(x , t ) = 0 ∂t

∇⋅σ = 0

T (x , t ) = T (x , t (x ∈ ΓT ))

(3)

q(x , t )⋅n = q (x , t ) (x ∈ Γq )

(4)

u(x) = u(x) (x ∈ Γu )

(5)

σ(x)⋅n = t(x) (x ∈ Γt )

(6)

with ΓT the temperature boundary, Γq the flux boundary, Γu the Γt displacement boundary and the traction boundary (ΓT ∪ Γq = Γu ∪ Γt = Γ , ΓT ∩ Γq = Γu ∩ Γt = ∅. The crack boundary Γc is a part of natural boundary and is assumed to be adiabatic in heat conduction analysis and traction-free in elastic analysis. T , q , u and t are, respectively, the prescribed temperature, flux, displacement and traction on corresponding boundary. n is the outward unit normal to the domain. The initial condition for Eq. (1) is

T (x , 0) =T0

(x ∈ Ω )

(7)

3. Thermal shock fracture analysis by the NMM

(1) (2)

where ρ is the density and c is the specific heat at constant pressure. ∂ denotes partial derivative. T (x , t ) is the temperature with x ∈ Ω and t the time. The heat flux q is determined by the Fourier's law as

T h (x , t ) =

nt

∑ wi(x)Ti (x, t ) i =1

uh(x) =

(8)

nt

∑ wi(x)ui(x) i =1

(9)

where nt is the amount of PPs shared by e. wi(x) is the PU weight function defined on the MP containing the ith PP. Ti (x , t ) and ui(x) are,

Fig. 1. A cracked homogeneous body in 2D transient linear thermoelasticity.

47

Engineering Analysis with Boundary Elements 75 (2017) 46–56

H.H. Zhang et al.

Fig. 2. A schematic of basic concepts in the NMM: (a) physical domain, (b) mathematical cover, (c) generation of physical patches and (d) generation of manifold elements.

ui(x) = Pu(x)ci + Φu(x)di

respectively, the thermal and mechanical local functions defined on the ith PP. For PPs containing no crack tip, i.e., the conventional PPs,

Ti (x , t ) = PT (x)ai(x , t )

(10)

ui(x) = Pu(x)ci

(11)

where bi and di are, respectively, the vector of additional thermal and mechanical DOFs defined on the ith PP. ΦT (x) is the thermal crack tip asymptotic basis being [38]

where ai and ci are, respectively, the vector of conventional thermal and mechanical degrees of freedom (DOFs) defined on the ith PP. PT (x) is the thermal polynomial basis being

PT (x) = [1 x y ⋯]

ΦT (x) =

θ 2

⎡ Φ 0 Φ2 0 Φ3 0 Φ4 0 ⎤ Φu(x) = ⎢ 1 ⎥ ⎣ 0 Φ1 0 Φ2 0 Φ3 0 Φ4 ⎦

(13)

For PPs containing a crack tip, i.e., the singular PPs [21],

Ti (x , t ) = PT (x)ai(x , t ) + ΦT (x)bi(x, t )

r sin

(16)

with (r, θ) the polar coordinates located at the crack tip as plotted in Fig. 1. Φu(x) is the mechanical crack tip asymptotic basis as [21]

(12)

and Pu(x) is the mechanical polynomial basis as

⎡ 1 0 x 0 y 0 ⋯⎤ Pu(x) = ⎢ ⎥ ⎣ 0 1 0 x 0 y ⋯⎦

(15)

with (14) 48

(17)

Engineering Analysis with Boundary Elements 75 (2017) 46–56

H.H. Zhang et al.

⎡ θ [Φ1, Φ2, Φ3, Φ4] = ⎢ r sin , ⎣ 2

r cos

θ , 2

r sin θ sin

θ , 2

r sin θ cos

θ⎤ ⎥ 2⎦ (18)

for the concerned problem. Through Eqs. (14)–(18), the local property around crack tip (e.g., the flux and stress singularity) can be well represented. In the case of the temperature and displacement discontinuity across the crack face, as proved in [21], it is naturally reflected owing to the use of independent local functions for PPs completely separated by the crack geometry. Therefore, the existence of cracks in linear thermoelasticity is totally manifested by the NMM.

∫Γ

e T

λT (T h − T )δT hdΓ −

h

e

e q

h

σ(u ): ε(δ u )dΩ +

NTi = [ wi PT wi ΦT ]

(31)

⎡ (wi PT ), x (wi ΦT ), x ⎤ ⎥ BTi = ⎢ ⎢⎣ (wi PT ), y (wi ΦT ), y ⎥⎦

(32)

K uU = Fu + F0

∫Γ

h

e u

h

λu(u − u)⋅δ u dΓ −

∫Γ

e t

Ku =

qδT hdΓ

(20)

δT h(x , t ) =

∫Γ

F0 =

∫Ω

i =1

∑ wi(x)δ ui(x)

(21)

(22)

i =1

∫Ω

KT =

∫Ω ∫Γ

e T

e

(23)

e

NTTρc NT dΩ

BTT k BT dΩ +

NTTλT T dΓ −

∫Γ

e T

∫Γ

e q

NTTλT NT dΓ

NTTqdΓ

e u

∫Γ

e t

NuTλu NudΓ

NuTtdΓ

BuT DεT dΩ

(38)

⎡ w 0 wx 0 wy 0 ⋯⎤ i i Niu = ⎢ i ⎥ 0 wy ⋯⎦ ⎣ 0 wi 0 wx i i

(39)

⎡ wi, x 0 (wx 0 (wy 0 ⋯⎤ i ), x i ), x ⎢ ⎥ 0 (wx ) 0 ( ) wy Biu = ⎢ 0 wi, y i ,y i , y ⋯⎥ ⎢ w w (wx ) (wx ) (wy ) (wy ) ⋯⎥ i ,y i ,x i ,y i ,x ⎣ i, y i, x ⎦

(40)

(42)

D=

(26)

⎡ 1 ν0 ⎤ 0 ⎢ ⎥ 1 0 ν 0 ⎢ ⎥ 1 − ν02 ⎢ ⎣ 0 0 (1 − ν0 )/2 ⎥⎦ E0

(43)

where

where the superscript T denotes the matrix transpose. The entries of NT and BT are

(28)

(41)

⎡ wi, x 0 (wx 0 (wy 0 ⋯ (wΦ 0 ⋯⎤ i ), x i ), x i 1), x ⎢ ⎥ 0 (wx 0 (wy 0 (wΦ Biu = ⎢ 0 wi, y i ), y i ), y ⋯ i 1), y ⋯⎥ ⎢ w w (wx ) (wx ) (wy ) (wy ) ⋯ (wΦ ) (wΦ ) ⋯⎥ i ,y i ,x i ,y i ,x i 1 ,y i 1 ,x ⎣ i, y i, x ⎦

(25)

BT = [ B1T BT2 … BTi … BTnt ]

(36)

Bu = [ B1u Bu2 … Biu … Bunt ]

(24)

(27)

(35)

(37)

for singular PP. D is the elasticity matrix being

NT = [ N1T NT2 … NTi … NTnt ]

(34)

Nu = [ N1u Nu2 … Niu … Nunt ]

⎡ w 0 wx 0 wy 0 ⋯ wΦ 0 ⋯⎤ i i i 1 Niu = ⎢ i ⎥ 0 wy ⋯ 0 wΦ ⎣ 0 wi 0 wx i i i 1 ⋯⎦

where T and Ṫ are, respectively, the vector of thermal unknowns and their time derivatives. CT , KT and FT are, respectively, the heat capacity matrix, the thermal conductivity matrix and the equivalent thermal load vector as

CT =

NuTλu udΓ +

∫Γ

for conventional PP and

On substituting Eqs. (8) and (21) into Eq. (19) and considering the arbitrariness of the variation of DOFs, the NMM discrete equations for thermal part of the concerned problem are derived as

CT Ṫ + KT T = FT

e

BuT DBu dΩ +

with

nt

δ uh(x) =

e u

e

where the entries of Nu and Bu are

nt

∑ wi(x)δTi (x, t )

∫Ω

Fu =

h

t⋅δ u dΓ = 0

(33)

where U is the vector of mechanical unknowns. K u , Fu and F0 are, respectively, the stiffness matrix, the equivalent mechanical load vector due to mechanical and thermal loading as

where λT and λu are the penalty numbers adopted to enforce the essential boundary conditions due to the inconsistency of MC with the physical boundary. Ω e and Γme (m denotes T, q, u and t) are, respectively, the domain and/or boundary occupied or shared by the ME e. Through Eqs. (8) and (9), the test functions δT h and δuh are expressed as

FT =

(30)

(19)

=0

∫Ω

∫Γ

⎡ (wi PT ), x ⎤ ⎥ BTi = ⎢ ⎢⎣ (wi PT ), y ⎥⎦

for singular PP. Similarly, by substituting Eqs. (9) and (22) into Eq. (20), the NMM formulations for elastic part of the problem are obtained as

The NMM discrete equations for the present problem can be derived using the weighted residual method in Galerkin form [39]. Let T ∈ H1(Ω ),u ∈ H1(Ω ) be the temperature and displacement trial function, and δT ∈ H1(Ω ), δ u ∈ H1(Ω ) be the corresponding test function with H1 the first Hilbert space and δ the first order variation. A weak form of the discrete problem on an ME e is to find T h and uh in the finite dimensional subspace V h ∈ H1(Ω ),∀ δT h ∈ V h, δ uh ∈ V h so that

∫

(29)

for conventional PP and

3.2. Discrete equations

⎞ ⎛ ∂T h ⎜ρc δ T + q(T h )k⋅q(δ T h)⎟dΩ + e⎝ ⎠ Ω ∂t

NTi = [wi PT ]

⎧ E plane stress E0 = ⎨ 2 E /(1 − ν ) plane strain ⎩

(44)

⎧ ν plane stress ν0 = ⎨ ⎩ ν /(1 − ν ) plane strain

(45)

with E the Young's modulus and ν the Poisson's ratio.

with 49

Engineering Analysis with Boundary Elements 75 (2017) 46–56

H.H. Zhang et al.

Defining

⎧T⎫ −1 G = − C−1 T K T , H = CT FT , X = ⎨ ⎬ , ⎩1⎭

(47)

we have

⎡ ̇⎤ ⎡ ⎤ Ẋ = ⎢ T ⎥ = ⎢ G H ⎥X ⎣0⎦ ⎣ 0 0 ⎦

(48)

⎡ ⎤ Further setting Q = ⎢ G H ⎥, we can obtain ⎣0 0⎦ Ẋ = QX

(49)

The solution to the above linear differential equations with constant coefficients is

X(t ) = exp(Qt )X 0

(50)

with X 0 = X(t )|t =0 . In the PTIM, X(t ) is calculated step by step. Assuming that X(tn ) at the instant tn is known, from Eq. (50), we can achieve

X(tn +1) = exp(QΔt )X(tn )

Fig. 3. An edge crack in a rectangular plate under thermal shock: (a) physical model and (b) discretized domain when a =0.5 and h =0.09 with 324 PPs and 282 MEs.

with Δt = tn +1 − tn the time step. To accurately compute X(t ) at concerned instant, the key lies in the calculation of the exponential matrix exp(QΔt ) in Eq. (51). To this end, the following precise computation scheme is adopted. Further dividing Δt equally into 2L intervals (L is a positive integral), we have

4. Computer implementations To solve Eq. (23), difference methods such as the forward difference scheme (i.e., Euler scheme), the central difference scheme (i.e., CrankNicolson scheme) and the backward difference scheme are mostly adopted [40]. However, the forward scheme is conditionally stable and may also be oscillating. As for the central and backward scheme, they both are unconditionally stable, whereas the former is also nonimmune to oscillation although its accuracy is higher at a given time step and mesh configuration. Summarily, for difference methods, the stability and accuracy of the solution are generally determined by the time step. The precise time integration method, originally proposed by Zhong and Williams [36], is absolutely stable, non-oscillating and commonly independent of time step in accuracy. Presently, the PTIM with dimensional expanding developed by Zhang and Deng [41] and adopted in [42] is applied to perform the time integration in Eq. (23). The corresponding procedure is demonstrated as follows. Through simple mathematical manipulations, Eq. (23) can be reexpressed as −1 Ṫ = − C−1 T K T T + CT FT

(51)

exp(QΔt ) = [exp(Qτ )]2

L

(52)

L

with τ = Δt /2 . Since τ is extremely small (take L = 20 as suggested in [36], then τ = Δt /220 ), through the second-order Taylor expansion technique, we obtain

exp(Qτ ) ≈ I′ + Qτ + (Qτ )2 /2 = I′ + A0

(53)

with I′ the identity matrix and

A0 = Qτ + (Qτ )2 /2

(54)

Substituting Eq. (53) into Eq. (52), we arrive

exp(QΔt ) ≈ (I′ + A0)2

L

(55)

In consideration that the elements of matrix A0 are also very small compared to those of I′, to reduce the round-off error, the following algorithm is taken.

(46)

Fig. 4. Temperatures at different instant for various time steps: (a) point A and (b) point B.

50

Engineering Analysis with Boundary Elements 75 (2017) 46–56

H.H. Zhang et al.

Table 1 The CPU time used by the PTIM and BDM.

Table 3 Temperatures of point B at different element size h and selected instants.

time step

PTIM

BDM

Time

h =0.25

h =0.15

h =0.09

h =0.06

h =0.045

0.02 0.04 0.08 0.10 0.20

14 s 10 s 8s 7s 6s

10 s – – – –

0.2 0.4 0.6 0.8 1.0

−0.42551 −0.51686 −0.53973 −0.54598 −0.54770

−0.41859 −0.51672 −0.54152 −0.54836 −0.55028

−0.41772 −0.51651 −0.54177 −0.54884 −0.55084

−0.41712 −0.51640 −0.54218 −0.54941 −0.55149

−0.41711 −0.51637 −0.54224 −0.54948 −0.55156

Note: the steady-state value when h =0.02 is −0.55283.

Simulation platform information: Intel XEON E5-1650 V2 3.5GHZ CPU, 32 GB RAM, Windows 7 OS Table 2 Temperatures of point A at different element size h and selected instants. Time

h =0.25

h =0.15

h =0.09

h =0.06

h =0.045

0.2 0.4 0.6 0.8 1.0

0.43183 0.56084 0.60596 0.61697 0.61914

0.43313 0.56836 0.60631 0.61893 0.61995

0.43440 0.56946 0.60824 0.61928 0.62242

0.43488 0.56997 0.60920 0.62047 0.62408

0.43495 0.57004 0.60929 0.62058 0.62417

Note: the steady-state value when h =0.02 is 0.62633. L

exp(QΔt ) ≈ (I′ + A0)2 = (I′ + A1)2 = (I′ + Ai)

2L − i

L −1

= (I′ + A2)2

=⋯

L −2

=⋯ (56)

where

Ai = 2 × Ai −1 + Ai −1 × Ai −1 (i = 1, 2, …, L )

(57) Fig. 5. Normalized mode II TSIFs at different instant and crack lengths.

After L times loop using Eq. (57), the elements of AL are no longer very small and we can obtain

exp(QΔt ) ≈ I′ + AL

(the edge length of a mathematical element) is also applied. The TSIFs are calculated with the domain form of the interaction integral [26,43].

(58)

without severe round-off error. Following the above procedure, the thermal unknowns at t = tn+1 are finally calculated by Eqs. (51) and (58). What needs to be emphasized is about the computation of the matrix G and H in Eq. (47). Due to the symmetry of matrix C (see Eq. (24)), with the frequently used matrix decomposition method (e.g., LDLT decomposition) and back substitution method, the calculation of inverse matrix involved in G and H can be avoided, which may bring down the simulation cost to some extent. Besides, to evaluate the corresponding matrix in Eqs. (23) and (33), the spatial integral scheme should also be provided. Presently, the numerical integral scheme is adopted. For the MEs containing no crack tip, they are firstly triangulated; then the widely used Gauss quadrature rules are used on each triangle; finally, the integrals on all the triangles add up to the final results on this ME. As for the MEs containing the crack tip, the Duffy transformation strategy described in [27] is taken.

5.1. An edge crack in a rectangular plate In the first example, as shown in Fig. 3a, a rectangular plate of size W × 2H with an edge crack of length a is considered to validate the PTIM and accuracy of TSIFs as well. The initial temperature of the plate is zero. From t = 0 , the temperatures ΔTH (t ) and −ΔTH (t ) are, respectively, applied on the top and bottom edges of the plate, with H (t ) the Heaviside step function, while the other boundaries are assumed to be thermally insulated. The displacement boundary condition is enforced such that the horizontal displacements of the left side and the vertical displacement at the center of the right side are set to be zero. Concerning the present configurations, the crack is theoretically of pure mode II. In the simulation, the plane strain condition is taken and the following parameters are adopted: W = H = 1.0, ΔT = 1.0 . Firstly, the effect of the time step on the distribution of temperature field is investigated. Herein, the crack length a is taken to be 0.5. MC with element size h=0.09 is used to cover the plate (the discretized domain is illustrated in Fig. 3b). In the computation, three kinds of time step with Δt = 0.08, 0.04 and 0.02 are inspected. In Fig. 4a and b, the calculated temperatures at point A (0.25, 0.50) and B (0.75, −0.50) by the PTIM are plotted together with those by the backward difference method (BDM). Clearly seen from these figures, the results by the PTIM are independent of the time step and tend to the steady-state value obtained by the methodology described in [26]. As for the BDM solutions, they are obviously time-step-dependent and closer to those by the PTIM with the decreasing of Δt . Then, the efficiency of the PTIM and BDM is compared under fixed total time i.e., 1.0 and the MC in Fig. 3b. The CPU time cost by these two schemes is listed in Table 1. As we can see, for the same time step with similar accuracy, e.g., Δt = 0.02 (see Fig. 4), the CPU time used by the former is somewhat more; however, since the accuracy of the PTIM

5. Numerical examples In this part, to verify the proposed method, three numerical examples are conducted. Firstly, a pure mode II edge crack in a rectangular plate is tested; then, two embedded mixed mode I-II cracks in a square plate is considered; finally, a symmetrical branched crack in a square plate is studied. Throughout the computations, constant polynomial basis (see Eqs. (12) and (13)) is considered. The same material with ρ = 1.0 ,c = 1.0 , k = 1.0 ,E = 1.0 ,ν = 0.3 and α = 1.0 is investigated. The temperature boundary conditions are prescribed with regard to the reference temperature, i.e., ΔT = T − TR . In the modeling, regular MPs consisting of square mathematical elements are used. The widely used finite element shape functions for rectangular element are adopted as the weight functions in Eqs. (8) and (9). To quantify the discretization, besides the number of PPs and MEs, the element size h 51

Engineering Analysis with Boundary Elements 75 (2017) 46–56

H.H. Zhang et al.

Fig. 6. Two embedded cracks in a square plate under thermal shock: (a) physical model and (b) discretized domain when β=45° and h =0.09 with 2033 PPs and 1892 MEs.

Fig. 7. Temperature distribution at different instant: (a) t=0.125, (b) t=0.5, (c) t=1.0 and (d) t=4.0.

52

Engineering Analysis with Boundary Elements 75 (2017) 46–56

H.H. Zhang et al.

Fig. 8. Normalized TSIFs at different instant when β=45°: (a) mode I and (b) mode II.

Fig. 9. Normalized TSIFs at different instant when β=0°: (a) mode I and (b) mode II.

Fig. 10. An embedded branched crack in a square plate under thermal shock: (a) physical model and (b) discretized domain when a =0.5 and h =0.045 with 2010 PPs and 1887 MEs.

are computed at different MC with h=0.25, 0.15, 0.09, 0.06 and 0.045. The results of point A and B at selected instants are listed in Tables 2 and 3, respectively. From these tables, we can see that with the refinement of the MC, the temperatures at each instant tend to be convergent; what's more, the convergent values at t = 1.0 are very close

is insensitive to Δt , for larger time step, e.g., 0.08, 0.10 and 0.20, its efficiency is higher than that of the BDM, which infers that the PTIM is of great accuracy and efficiency as well. Next, we focus on the convergence test of thermal analysis with the PTIM in the framework of the NMM. For this purpose, temperatures

53

Engineering Analysis with Boundary Elements 75 (2017) 46–56

H.H. Zhang et al.

Fig. 11. Temperature distribution at selected instant when a =0.5: (a) t=0.12, (b) t=0.24, (c) t=0.40 and (d) t=1.0.

In the modeling, we take L=2.0, a=1.0, ΔT=1.0 and q=ΔT/L. Two inclinations with β=45° and β=0° are, respectively, tested. Throughout the computations, the MC with h=0.09is used (see Fig. 6b for the discretization when β=45°). Firstly, we care about the accuracy of thermal analysis. The calculated temperature fields for β=45° at four instants, i.e., t=0.125, 0.5, 1.0 and 4.0, are plotted in contour form in Fig. 7a~d. For comparison, the results provided by Chang and Ma in [44] are also displayed. It is obvious that the temperature distributions obtained by the NMM and PTIM match well with those by the analytical alternating method. What's more, the temperature jump across crack faces can also be inferred from these figures. Next, the stress intensity factors for both inclinations are computed. The mixed mode TSIFs are normalized by αΔTE a and plotted in Fig. 8 for the case when β=45° and in Fig. 9 for β=0°. Considering that there are no reference solutions for this problem presently, we only compare the stable results with those obtained by the steady-state thermoelastic NMM described in [26]. Apparently, the transient TSIFs are closer to the steady-state results with the increase of time. In addition, it's observed that at some instants (e.g., t=1.0 for crack tip D in Fig. 8a), the mode I TSIFs are negative, suggesting the occurrence of crack

to the steady-state solutions at a finer MC. Finally, we study the accuracy of the proposed method in the computation of TSIFs. The TSIFs at different crack length with a =0.25, 0.50 and 0.75 are calculated using the fixed MC shown in Fig. 3b. The mode II TSIFs are normalized by αΔTE W and depicted in Fig. 5, in conjunction with the available reference solutions by the FEM [7] and the BEM [10]. From this figure, a nice agreement among the present and existing results is found. 5.2. Two embedded cracks in a square plate In this case, a square plate with two embedded cracks is examined in the plane stress state. As shown in Fig. 6a, the side length of the plate is 2L. The two cracks are of equal length 2a and their center is located at (−0.5a, −0.5a) and (0.5a, 0.5a), respectively. The included angle between the crack geometry and the x-axis is β. The initial temperature of the domain is set to be zero. From the instant t = 0 , temperatures are prescribed on the top and bottom side of the plate with the value being ΔTH (t ) and −ΔTH (t ), respectively, while the heat flux qH (t ) is applied on the other sides. As for the displacement boundary condition, both the left and bottom edges are clamped. 54

Engineering Analysis with Boundary Elements 75 (2017) 46–56

H.H. Zhang et al.

Fig. 12. Normalized TSIFs at different instant: (a) a =0.3, (b) a =0.4 and (c) a =0.5.

identical). As we can see from these figures, with the increase of time, the TSIFs move stably to the steady-state results.

closure. From the theory of fracture mechanics, such kind of results is not valid and can only be served in a superposition which may yield a positive solution [45]. 5.3. An embedded branched crack in a square plate

6. Concluding remarks In the last example, to show the power of the proposed approach for more complex crack problems, a square plate with an embedded branched crack as shown in Fig. 10a is considered. The plate dimension, crack geometries, temperature and displacement boundary conditions are all denoted in Fig. 10a. Similarly to the above two examples, the initial temperature is also taken to be zero. To the authors’ knowledge, there are no existing results in literature for both thermal and elastic analysis. In the computation, the plane stress state is supposed. L and ΔT are set to be unity. Three kinds of crack length with a =0.3, 0.4 and 0.5 are, respectively, investigated on a fixed MC with h =0.045 (the discretized domain when a =0.5 is presented in Fig. 10b). The calculated temperature fields at selected instants with t=0.12, 0.24, 0.4 and 1.0 when a =0.5 are graphically provided in Fig. 11a~d, in which the temperature discontinuity across crack faces is clearly exhibited. The computed TSIFs at crack tips A and B are normalized by αΔTE W and plotted in Fig. 12a~c for different crack lengths (theoretically, due to the symmetry, the mode I TSIFs at crack tip B are equal in size and opposite in sign with those at crack tip C, while their mode II TSIFs are

In this work, the numerical manifold method, combined with the precise time integration scheme, has been developed to tackle 2D crack problems subjected to thermal shock. The temperature and displacement discontinuity across crack surfaces were naturally represented due to the introduction of dual cover systems in the NMM. The singularity of field gradients at crack tip was described through the use of asymptotic basis in the local functions on singular PPs. The time integration in transient heat conduction analysis was realized with the PTIM and then the thermal results were transported into the elastic part to calculate the TSIFs. The proposed method was validated through three numerical examples, including single, double and multi-branched cracks. The advantage of the PTIM (e.g., time-stepindependence, high accuracy and efficiency) in the framework of the NMM was numerically proved. The calculated TSIFs by the present method are in good accordance with previously published results. The present approach can be extended to the dynamic fracture analysis of homogeneous and nonhomogeneous materials under thermal shock and the results are awaited. 55

Engineering Analysis with Boundary Elements 75 (2017) 46–56

H.H. Zhang et al.

crack modeling. Theor Appl Fract Mec 2005;44(3):234–48. [21] Ma GW, An XM, Zhang HH, Li LX. Modeling complex crack problems using the numerical manifold method. Int J Fracture 2009;156(1):21–35. [22] Zhang HH, Li LX, An XM, Ma GW. Numerical analysis of 2-D crack propagation problems using the numerical manifold method. Eng Anal Bound Elem 2010;34(1):41–50. [23] Gao HF, Cheng YM. A complex variable meshless manifold method for fracture problems. Int J Comp Meth-Sing 2010;7(1):55–81. [24] Zhang HH, Zhang SQ. Extract of stress intensity factors on honeycomb elements by the numerical manifold method. Finite Elem Anal Des 2012;59:55–65. [25] Wu ZJ, Wong LNY. Elastic-plastic cracking analysis for brittle-ductile rocks using manifold method. Int J Fracture 2013;180(1):71–91. [26] Zhang HH, Ma GW, Ren F. Implementation of the numerical manifold method for thermo-mechanical fracture of planar solids. Eng Anal Bound Elem 2014;44:45–54. [27] Zheng H, Xu DD. New strategies for some issues of numerical manifold method in simulation of crack propagation. Int J Numer Meth Eng 2014;97(13):986–1010. [28] Zheng H, Liu F, Du XL. Complementarity problem arising from static growth of multiple cracks and MLS-based numerical manifold method. Comput Method Appl M 2015;295:150–71. [29] He J, Liu QS, Ma GW, Zeng B. An improved numerical manifold method incorporating hybrid crack element for crack propagation simulation. Int J Fracture 2016;199(1):21–38. [30] Yang YT, Zheng H. A three-node triangular element fitted to numerical manifold method with continuous nodal stress for crack analysis. Eng Fract Mech 2016;162:51–75. [31] Terada K, Ishii T, Kyoya T, Kishino Y. Finite cover method for progressive failure with cohesive zone fracture in heterogeneous solids and structures. Comput Mech 2007;39(2):191–210. [32] Kurumatani M, Terada K. Finite cover method with multi-cover layers for the analysis of evolving discontinuities in heterogeneous media. Int J Numer Meth Eng 2009;79(1):1–24. [33] An XM, Zhao ZY, Zhang HH, He L. Modeling bimaterial interface cracks using the numerical manifold method. Eng Anal Bound Elem 2013;37(2):464–74. [34] Zhang HH, Ma GW. Fracture modeling of isotropic functionally graded materials by the numerical manifold method. Eng Anal Bound Elem 2014;38:61–71. [35] Wu ZJ, Wong LNY. Modeling cracking behavior of rock mass containing inclusions using the enriched numerical manifold method. Eng Geol 2013;162:1–13. [36] Zhong WX, Williams FW. A precise time-step integration method. P I Mech Eng CJ Mec 1994;208(6):427–30. [37] Melenk JM, Babuska I. The partition of unity finite element method: Basic theory and applications. Comput Method Appl M 1996;139(1–4):289–314. [38] Duflot M. The extended finite element method in thermoelastic fracture mechanics. Int J Numer Meth Eng 2008;74(5):827–47. [39] Lin JS. A mesh-based partition of unity method for discontinuity modeling. Comput Method Appl M 2003;192(11–12):1515–32. [40] Cebeci T. Convective heat transfer, 2nd Revised edition. Horizons Publishing and Springer; 2002. [41] Zhang SY, Deng ZC. Increment-dimensional precise integration method for nonlinear dynamic equation. Chin J Comput Mech 2003;20(4):423–6. [42] Yu TT, Gong ZW. Numerical simulation of temperature field in heterogeneous material with the XFEM. Arch Civ Mech Eng 2013;13(2):199–208. [43] Zamani A, Gracie R, Eslami MR. Higher order tip enrichment of eXtended Finite Element Method in thermoelasticity. Comput Mech 2010;46(6):851–66. [44] Chang CY, Ma CC. Transient thermal conduction analysis of a rectangular plate with multiple insulated cracks by the alternating method. Int J Heat Mass Transf 2001;44(13):2423–37. [45] Konda N, Erdogan F. The mixed-mode crack problem in a nonhomogeneous elastic medium. Eng Fract Mech 1994;47(4):533–45.

Acknowledgements The present work was supported by the National Natural Science Foundation of China (Grant Nos. 11462014, 11572282), the Provincial Natural Science Foundation of Jiangxi, China (Grant No. 20151BAB202003), the Science and Technology Program of Educational Committee of Jiangxi Province of China (Grant No. GJJ14526). References [1] Oliveira R, Wu XR. Stress intensity factors for axial cracks in hollow cylinders subjected to thermal-shock. Eng Fract Mech 1987;27(2):185–97. [2] Noda N, Ashida F. Thermal-shock in a transversely isotropic cylinder containing an annular crack. Int J Solids Struct 1993;30(3):427–40. [3] Lee KY, Kim JS. Determination of thermal shock stress intensity factor for elliptical crack by modified Vainshtok's weight function method. Eng Fract Mech 1997;56(3):423–35. [4] Noda N, Wang BL. Transient thermoelastic responses of functionally graded materials containing collinear cracks. Eng Fract Mech 2002;69(14–16):1791–809. [5] Shahani AR, Nabavi SM. Transient thermal stress intensity factors for an internal longitudinal semi-elliptical crack in a thick-walled cylinder. Eng Fract Mech 2007;74(16):2585–602. [6] Wang BL, Li JE. Hyperbolic heat conduction and associated transient thermal fracture for a piezoelectric material layer. Int J Solids Struct 2013;50(9):1415–24. [7] Emmel E, Stamm H. Calculation of stress intensity factors of thermally loaded cracks using the finite-element method. Int J Pres Ves Pip 1985;19(1):1–17. [8] Magalhaes JF, Emery AF. Transient thermoelastic fracture of brittle substrates bonded to brittle films. J Therm Stresses 1997;20(1):35–45. [9] Kim JS, Koo BG, Choi JB, Kim YJ. A finite element study on the integrity evaluation method of subclad cracks under pressurized thermal shock transients. J Press VessT Asme 2003;125(1):46–51. [10] Prasad NNV, Aliabadi MH, Rooke DP. The dual boundary element method for transient thermoelastic crack problems. Int J Solids Struct 1996;33(19):2695–718. [11] Giannopoulos GI, Anifantis NK. Interfacial steady-state and transient thermal fracture of dissimilar media using the boundary element contact analysis. Int J Numer Meth Eng 2005;62(10):1399–420. [12] Zamani A, Eslami MR. Implementation of the extended finite element method for dynamic thermoelastic fracture initiation. Int J Solids Struct 2010;47(10):1392–404. [13] Rokhi MM, Shariati M. Implementation of the extended finite element method for coupled dynamic thermoelastic fracture of a functionally graded cracked layer. J Braz Soc Mech Sci 2013;35(2):69–81. [14] ShiG.H.. Manifold method of material analysis. Transcations of the 9th Army Confernece on Applied Mathematics and Computing 1991;57-76. [15] Kim J, Bathe KJ. The finite element method enriched by interpolation covers. Comput Struct 2013;116:35–49. [16] Kim J, Bathe KJ. Towards a procedure to automatically improve finite element solutions by interpolation covers. Comput Struct 2014;131:81–97. [17] Jeon HM, Lee PS, Bathe KJ. The MITC3 shell finite element enriched by interpolation covers. Comput Struct 2014;134:128–42. [18] Ma GW, An XM, He L. The numerical manifold method: a review. Int J Comp Meth-Sing 2010;7(1):1–32. [19] Tsay RJ, Chiou YJ, Chuang WL. Crack growth prediction by manifold method. J Eng Mech-Asce 1999;125(8):884–90. [20] Li SC, Li SC, Cheng YM. Enriched meshless manifold method for two-dimensional

56