Application of a hybrid mesh-free method for shock-induced thermoelastic wave propagation analysis in a layered functionally graded thick hollow cylinder with nonlinear grading patterns

Engineering Analysis with Boundary Elements 43 (2014) 56–66

Contents lists available at ScienceDirect

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

Application of a hybrid mesh-free method for shock-induced thermoelastic wave propagation analysis in a layered functionally graded thick hollow cylinder with nonlinear grading patterns Seyed Mahmoud Hosseini Industrial Engineering Department, Faculty of Engineering, Ferdowsi University of Mashhad, PO Box 91175-1111, Mashhad, Iran

art ic l e i nf o

a b s t r a c t

Article history: Received 15 January 2014 Accepted 17 March 2014

This article exploits a hybrid mesh-free method for coupled thermoelasticity analysis (without energy dissipation) and thermoelastic wave propagation analysis in layered FGMs subjected to shock loading. The presented hybrid mesh-free method is based on generalized finite difference (GFD) and Newmark finite difference (NFD) methods. The Green–Naghdi (GN) theory of coupled thermoelasticity is assumed to derive the governing equations for FG thick hollow cylinder. The layered FG cylinder is assumed to be under thermal shock loading. The mechanical properties of layered FG cylinder are considered to vary along the radial direction as nonlinear functions in terms of volume fraction. Thermoelastic wave propagations are studied in details at various time instants for various grading patterns of mechanical properties. The effects of nonlinear grading patterns on thermoelastic wave propagations are obtained and discussed using the presented effective mesh-free method. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Mesh-free method Generalized coupled thermoelasticity Generalized finite difference (GFD) method Thick hollow cylinder Thermal shock Functionally graded materials (FGMs)

1. Introduction A very well-known kind of composites, which are widely used in thermal applications, is functionally graded materials (FGMs). In these kinds of materials, the mechanical properties are varied along a certain direction as linear or nonlinear functions. In the recent years, some research works were carried out to study on dynamic analysis, transient heat conduction and wave propagation of FGMs. Some accurate results were presented by employing Carrera unified formulation (CUF) [1,2] for dynamic analysis of composite and FG structures. The Carrera unified formulation (CUF) and the Reissner's mixed variational theorem (RMVT) were extended for bending analysis of functionally graded plates by Brischetto and Carrera [3]. The ‘thermal shock’ analysis for a thin, cylindrical plate made of functionally graded ceramics was presented by Nakonieczny and Sadowski [4] using a meshfree, semidiscrete finite element method. In their work, the unsteady heat conduction equation was discretized in space with the partition of unity FEM and the temporal discretization was realized by using an explicit finite difference method. Youssef [5] presented an analytical method to study a problem of thermoelastic interactions in an elastic infinite medium with

E-mail addresses: [email protected], [email protected] http://dx.doi.org/10.1016/j.enganabound.2014.03.007 0955-7997/& 2014 Elsevier Ltd. All rights reserved.

cylindrical cavity subjected to thermal shock and moving heat source with constant velocity using a direct approach in the Laplace transforms domain. The formulation based on the statespace approach was derived and then applied to solve a boundary value problem of thermoelastic interactions with energy dissipation on the basis of the theory of generalized thermoelasticity type-III (Green–Naghdi theory) of an isotropic elastic half space by Mukhopadhyay et al. [6]. Although the analytical methods are very valuable even for coupled problems in engineering (for example coupled thermoelasticity problem [7]), some numerical methods such as finite element (FE) [8–10] have been successfully developed in this area of engineering because of mathematical limitations of analytical methods. Also, some other numerical methods were developed especially for coupled thermoelasticity problem. The locally transversal linearization (LTL) technique with the numerical inverse Laplace transform method were successfully used to solve GN coupled thermoelasticity in an annulus by Taheri et al. [11]. The lasergeneration of ultrasound induced-coupled thermoelasticity was studied by Veres et al. [12] using numerical solution based on the finite difference method. The two sets of equations were decoupled and therefore a combination of implicit and explicit finite difference techniques was applied. In spite of the great success of the finite element and boundary element methods as the most efficient numerical tools for the

S.M. Hosseini / Engineering Analysis with Boundary Elements 43 (2014) 56–66

solution of coupled problems in complex domains, there has been a growing interest in the so-called meshless and mesh-free methods over the past decade. One of the most efficient meshless techniques is the meshless local Petrov–Galerkin (MLPG) method [13–16]. Recently, the author employed this method successfully for coupled thermoelasticity analysis [17] and coupled non-Fick diffusion–elasticity [18]. There is other mesh-free method based on generalized finite difference (GFD) method, which was successfully used for some engineering problems. Benito et al. [19,20] developed the GFD method and in another work, they analyzed the possibility of employing the GFD method over adaptive clouds of points progressively increasing the number of nodes. A procedure in the GFD method was given that easily assures the quality of numerical results by obtaining the residual at each point by Gavete et al. [21]. Also, they compared the GFD method with another meshless method the, so-called, element free Galerkin method (EFG). Benito et al. [22] applied the GFD method to solve the parabolic and hyperbolic differential equations. In their research, some examples from two kinds of differential equations were solved. The application of generalized finite difference (GFD) method to the problem of seismic wave propagation was presented by Ureña et al. [23] in which independent stability conditions and star dispersion of the phase velocity were obtained for the P and S waves. Dynamic analysis of beams and plates was successfully carried out using the generalized finite difference (GFD) method with irregular clouds of nodes by Gavete et al. [24]. Recently, the author used a hybrid mesh-free method based on generalized finite difference (GFD) method for elastic wave propagation analysis (elastodynamic analysis) in a functionally graded (FG) thick hollow cylinder [25]. Also, the application of GFD based mesh-free method was developed for coupled thermoelasticity analysis considering GN theory for homogenous isotropic cylinder by the author in [26]. To develop the application of GFD based mesh-free method for nonhomogenous materials and thermoelasticity analysis, the presented article employs the above hybrid mesh-free method to study thermoelastic wave propagation in a layered FG thick hollow cylinder subjected to thermal shock loading. The presented hybrid mesh-free method is based on the generalized finite difference (GFD) method and Newmark finite difference (NFD) method for temporal section. The thermoelastic wave propagation through both displacement and temperature fields are studied in details. The influence of various grading patterns on dynamic behaviors of displacement and temperature fields and also on propagation of thermoelastic wave fronts are obtained.

57

The coupled thermoelsticity (without energy dissipation) based on the Green–Naghdi theory (type II) is used for the problem formulations.

2. Layered functionally graded materials In layered functionally graded cylinder, the cylinder is divided in to many homogenous isotropic sub-cylinders along the radial direction of thickness. By suitable arrangement of layers, the grading patterns in mechanical properties of FGMs can be created. The mechanical properties of layered FG cylinder are varied from one bounding surface to other bounding surface through the radial direction. The layered FG cylinders are usually made by arrangement of some isotropic layers, which create the gradation of two distinct material phases where one of them is ceramic and the other one usually is metal alloy phases, and the volume fractions of the constituents vary as a nonlinear function. The layered FG cylinder is considered to be made of a combined metal–ceramic material for which the mixing ratio is varied in the radial direction from pure ceramic to pure metal. In the present problem, the bounding surfaces of the layered FG cylinder are made of materials with mechanical properties P(rin) for inner surface and P(rout) for the outer surface. Also, it is assumed that the terms “r in ” and “r out ” are inner and outer radii layered FG cylinder, respectively. To show the variation of mechanical properties through the radial direction, the nonlinear volume fraction is employed to simulate the gradation of all mechanical properties. The volume fraction nonlinear functions can be explained as follows: V in ¼ ½1  λnr r 

ð1Þ

V out ¼ ½λnr r 

ð2Þ

where
 r  r in nr λnr r ¼ r out  r in

ð3Þ

The material properties at an arbitrary point in the layered FG cylinder can be simulated using the rule of mixtures as follows: PðrÞ ¼ Pðr in Þ V in þ Pðr out ÞV out

ð4Þ

PðrÞ ¼ Pðr in Þ ½1  λnr r  þ Pðr out Þ ½λnr r 

ð5Þ

0.7 Volume fraction exponent nr = 0.1

0.68

Volume fraction exponent nr = 0.25 Volume fraction exponent nr = 0.5

0.66

Volume fraction exponent nr = 1

0.64

Volume fraction exponent nr = 1.5 Volume fraction exponent nr = 3

Cp

0.62

Volume fraction exponent nr = 6

0.6 0.58 0.56 0.54 0.52 0.5

1

1.05

1.1

1.15

1.2

1.25

1.3

1.35

1.4

1.45

Non-dimensional radius Fig. 1. Various grading patterns for a sample mechanical property in a layered FG cylinder.

1.5

58

S.M. Hosseini / Engineering Analysis with Boundary Elements 43 (2014) 56–66

The mechanical properties of each homogenous isotropic layer for example Jth layer can be obtained as follows:  n ðJ  1ÞH 0 r þ Pðr in Þ ð6Þ P J ¼ ½Pðr out Þ  Pðr in Þ ðr out r in Þ where the term nr stands for volume fraction exponent. H 0 ¼ ðr out r in Þ=M

ð7Þ

The term M stands for number of homogenous isotropic subcylinders. The obtained nonlinear grading patterns using Eq. (6) for a sample mechanical property (cp ) and various values of volume fraction exponents are presented in Fig. 1.

3. Coupled thermoelasticity governing equations To find the dynamic response of displacement and temperature fields, the coupled thermo-elasticity governing equations are employed in the problem. In coupled thermo-elasticity, the time rate of the first invariant change of strain tensor is employed in the first law of thermodynamics. Consequently, the dependency between the temperature and displacement fields is caused and thus the coupling between elasticity and energy equations. To find the dynamic behaviors of displacement field, the temperature field and elastic field should be solved together as a coupled system of partial differential equations. One of the most important theories in thermo-elasticity is Green and Naghdi (GN) theory of generalized coupled thermo-elasticity [27]. The governing coupled thermo-elasticity equations based on GN theory without energy dissipation are given for every homogenous isotropic layer as € ∇:r þρ F ¼ ρ u

ð8Þ

n cT€ þγT 0 ∇:u€ ¼ ρg_ þ ∇:ðk ∇TÞ

ð9Þ

where “u” is the displacement vector, “T” is the temperature change above the uniform reference temperature “T 0 ”, “F” is the _ is the external rate of supply of heat. Both external force and “g” _ are assumed to be absent in this work. Furtherthe “F” and “g” more, “ρ” is the mass density, “c” is the specific heat, “λ” and “μ” are the Lame constants and γ ¼ ð3λ þ 2μÞβn

ð10Þ n

n

where “β ” is the coefficient of linear thermal expansion and “k ” is a material parameter of the GN theory. The dot over the symbol denotes the time differentiation. The stress field can be reached using following equation:

sij ¼ δij ðλ uk;k  γTÞ þ μ ðui;j þuj;i Þ

ð11Þ

The axi-symmetry and plane strain conditions are assumed for the problem. Consequently, the following relations are taken into account to calculate the parameters: uθ ¼ 0; uz ¼ 0; srr ¼ 2 μ ur;r þ ðλ e γTÞ; sθθ ¼ 2 μ ur =r þ ðλ e γTÞ; szz ¼ λ e γT; srθ ¼ 0 ¼ srz ¼ szθ ; e ¼ ur;r þur =r: The governing Eqs. (8) and (9) are reduced to the equations μ∇2 u þ ðλ þ μÞ∇divu  γ∇T þ ρF ¼ ρu€

ð12Þ

cT€ þγ T 0 divu€ ¼ ρ g_ þ k ∇2 T

ð13Þ

n

To analyze the problem, the following non-dimensional parameters are employed for the problem. r v 1 ðλ þ 2μÞ T sr sθ r ¼ ; t ¼ t; u ¼ u; T ¼ ; sr ¼ ; sθ ¼ l l l γ T0 T0 γ T0 γ T0

ð14Þ

where “l” is a standard length and “v” is a standard speed. The governing equation and heat transfer can be rewritten for each

sub-cylinder in layered FG cylinder by using nondimensional parameters. C 2s ∇2 u þ ðC 2p  C 2s Þ∇divu  C 2p ∇T ¼ u€

ð15Þ

C 2T ∇2 T ¼ T€ þ εn divu€

ð16Þ

where C 2p ¼

n

λ þ 2μ 2 k μ γ2 T 0 ; ; CT ¼ 2; Cs2 ¼ ; εn ¼ 2 2 c ðλ þ 2μÞ ρv ρv cv

ð17Þ

The axisymmetry and plane strain conditions are assumed for the problem. The governing equations for axisymmetry and plane strain conditions can be obtained as follows: ∂2 u

1 ∂u u ∂T  ðC 2p  C s 2 Þ 2  C 2p ¼ u€ þ Cp2 ∂r r ∂r r ∂r 2 " # " # ∂2 T 1 ∂T ∂u€ u€ þ C 2T þ ¼ T€ þ εn 2 ∂r r r ∂r ∂r

C 2p

ð18Þ

ð19Þ

The following continuity conditions should be satisfied in interfaces between each two layers such as Jth and (J þ1)th layers in layered FG cylinder. uJ ¼ uJ þ 1 J

T ¼T

J þ1

;

sJr ¼ sJr þ 1

ð20Þ

;

J Jþ1 q_ ¼ q_

ð21Þ

where q_ is the heat flux rate, which is defined as follows: ∂T q_ ¼ C T ∂r

ð22Þ

To solve the above coupled thermo-elasticity governing Eqs. (18) and (19), a hybrid mesh-free method based on generalized finite difference (GFD) method is developed in the present article. 4. Generalized finite difference (GFD) method In GFD method, the partial derivatives are linearly approximated by Taylor series expansion on some nodes (center nodes) in the analyzed domain that each center node is surrounded by some other nodes. The nodes can be regular or randomly distributed on analyzed domain. Consequently, partial derivatives are obtained at the rest of each center nodes and the group of nodes with a center node and surrounding other nodes is called a star in this method. Consider the non-dimensional radial displacement at a center node to be “u0 ” and non-dimensional temperature to be “T 0 ” and the terms “ui ” and “T i ” are the values of non-dimensional radial displacement and temperature at the rest of surrounding nodes. The function values “ui ” and “T i ” can be approximated using Taylor expansion as
 ∂ u0 1 2 ∂ 2 u0 þ hi ui ¼ u0 þ hi þ⋯ ð23Þ 2 ∂r 2 ∂r and T i ¼ T 0 þ hi

! ∂ T 0 1 2 ∂2 T 0 hi þ þ⋯ ∂r 2 ∂ r2

ð24Þ

The term “i” is number of surrounding nodes. The analyzed domain in the problem is linear and through the radial direction of thickness of cylinder. Consequently, the term “hi ” can be calculated as hi ¼ r i  r o

ð25Þ

The terms over second order are ignored in Eqs. (23) and (24) and the linear approximation of second order can be obtained for

S.M. Hosseini / Engineering Analysis with Boundary Elements 43 (2014) 56–66

radial displacement and temperature. To minimize the error in this method, the function of norm should be minimized. The functions of norm for radial displacement and temperature are
 
2 ∂ u0 1 2 ∂ 2 u0 þ hi wðh NormðuÞ ¼ ∑ u0  ui þ hi Þ i ∂r 2 ∂ r2 i¼1 N

ð26Þ

(

∂2 T 0 ∂r

¼ BT1 2

)

N

2

∑ ð  T 0 þ T i Þhi w ðhi Þ

i¼1

N

NormðTÞ ¼ ∑

T 0  T i þ hi

i¼1

∂ T 0 1 2 ∂2 T 0 h þ ∂ r 2 i ∂ r2

!!

where the coefficients “BT2 ” are obtained as follows: N

#2 ð27Þ

Au1 ¼

wðhi Þ ¼

ðdistÞ3

¼

i¼1

ð28Þ

3

hi

¼ ξu2

ð30Þ



∂ u0 ∂2 u0 ; Q u2 ¼ ∂ r ∂ r2 2

T

∂ T0 ∂ T0 ; ∂ r ∂ r2 2

N

ð31Þ

∑

2

i¼1 N

2 hi w2 ðhi Þ

N

3

hi 2 w ðhi Þ 7 7 i¼12 7 7 N h4 5 ∑ i w2 ðhi Þ i¼14

∑

i¼1

∑

N

2 hi w2 ðhi Þ

ð33Þ

∑

i¼1

2 hi w2 ðhi Þ

AT2 ¼

∑

ð34Þ

i¼1

N

BT1 ¼ ð35Þ

∂T 0 ¼ AT1 ∂r

)

(

∑ ð  u0 þ ui Þhi w2 ðhi Þ  Bu2

N

2

∑ ð u0 þ ui Þ

i¼1

(

N

)

∑ ð T 0 þ T i Þhi w2 ðhi Þ  AT2

i¼1

N

∑

i¼1

∑

i¼1

(

N

∑ ð  T 0 þ T iÞ

i¼1

N

∑

)

hi 2 w ðhi Þ 2

) 2 hi 2 w ðhi Þ 2 ð38Þ

! 

2 hi w2 ðhi Þ

i¼1

! 2 hi w2 ðhi Þ

!

N

!2

ð44Þ

∑

!2

ð45Þ

!2

ð46Þ

!2

ð47Þ

3

hi 2 2 w ðhi Þ



N

∑

i¼1

3

hi 2 2 w ðhi Þ

3

hi 2 2 w ðhi Þ

!

hi 2 4 w ðhi Þ



N

∑

i¼1

3

hi 2 2 w ðhi Þ

!

∑

i¼1

ð43Þ

!

i¼1 N

∑

i¼1

hi 2 4 w ðhi Þ

N

∑



4

∑

!2 3

hi 2 2 w ðhi Þ

h3i 2 2 w ðhi Þ

N

i¼1

N

i¼1

!

∑

N

ð42Þ

hi 2 4 w ðhi Þ

i¼1

!

∑

!2

2 hi w2 ðhi Þ

4

i¼1

 BT2 ¼

ð37Þ and

2 hi w2 ðhi Þ

N

h3i 2 2 w ðhi Þ

!

i¼1

!



ð36Þ ¼ Bu1

N

N

i¼1

4

∑

∑

!

N

(



hi 2 4 w ðhi Þ

i¼1

ð41Þ

4

∑

The first and second derivatives are calculated as ( ) ( ) 2 N N ∂u0 hi 2 u u 2 ¼ A1 ∑ ð  u0 þ ui Þhi w ðhi Þ  A2 ∑ ð  u0 þ ui Þ w ðhi Þ ∂r 2 i¼1 i¼1

!

4

N

!2 3

hi 2 2 w ðhi Þ

i¼1

hi 2 4 w ðhi Þ

i¼1

!



N

hi 2 2 w ðhi Þ

and also

N

ð40Þ

!

∑

∑

∑

!2 3

hi 2 2 w ðhi Þ

!

i¼1

i¼1

and

3

N

N

N

i¼1

h4i 2 4 w ðhi Þ

∑

!



!

∑

i¼1

“BT1 ”,

h3i 2 2 w ðhi Þ

N

N

2 hi w2 ðhi Þ

“AT2 ”,

! !

i¼1

!

“AT1 ”,

hi 2 4 w ðhi Þ

4

i¼1

 Bu2 ¼

AT1 ¼

2

“Bu2 ”,

hi 2 4 w ðhi Þ

N

∑

“Bu1 ”,

N

!

∑

3 N 2 6 ∑ ð T 0 þ T i Þhi w ðhi Þ 7 6 i¼1 7 7 ξT2 ¼ 6 2 6 N 7 h 4 ∑ ð  T þ T Þ i w2 ðh Þ 5 0 i i 2 i¼1

∂r

∑

!

“Au2 ”,

hi 2 4 w ðhi Þ

i¼1

3

2

2

N

∑

ð32Þ

) 2 hi 2 ∑ ð  T 0 þ T i Þ w ðhi Þ 2 i¼1

4

 Bu1 ¼

i¼1

3 N 2 6 ∑ ð u0 þ ui Þhi w ðhi Þ 7 6 i¼1 7 7 ξu2 ¼ 6 2 6 N 7 4 ∑ ð  u þ u Þhi w2 ðh Þ 5 0 i i 2 i¼1

∂2 u0

N

i¼1

)T

2 6 ∑ hi w ðhi Þ 6 i ¼ 1 ψ u2 ¼ ψ T2 ¼ 6 6 N h3 4 ∑ i w2 ðhi Þ i¼12

and

2 hi w2 ðhi Þ

∑

Au2 ¼

where the terms “ψ u2 ” and “ψ T2 ” stand for 2  2 matrices in displacement and temperature fields, respectively. The vectors “Q u2 ” and “Q T2 ” are given, respectively, by

Q T2 ¼

N

ð29Þ

ψ T2 Q T2 ¼ ξT2

(

i¼1

!

N

1

If the norms (26) and (27) are minimized with respect to the partial derivatives, a set of linear equations system is obtained as follows: ψ u2 Q u2

N

∑

where “wðhi Þ” is the weight function. In this article, the weight function is defined by 1

(

4

∑

wðhi Þ

 BT2

ð39Þ “Au1 ”,

and "

59

2 hi w2 ðhi Þ

!

4

hi 2 4 w ðhi Þ



N

∑

i¼1

3

hi 2 2 w ðhi Þ

The derivatives of radial displacement and temperature can be also rewritten in star forms as follows: N ∂ u0 ¼  α0 u0 þ ∑ αi ui ∂r i¼1

ð48Þ

where N

2

N

i¼1

2

hi 2 w ðhi Þ i¼1 2

α0 ¼ Au1 ∑ hi w2 ðhi Þ  Au2 ∑

60

S.M. Hosseini / Engineering Analysis with Boundary Elements 43 (2014) 56–66

2 hi

αi ¼ Au1

w

2

(

2 h ðhi Þ  Au2 i w2 ðhi Þ

2

N

α0 ¼ ∑ αi

ð49Þ

i¼1

For second derivative of radial displacement, it is 2

∂ u0

N

¼  β 0 u0 þ ∑ β i ui

∂r 2

N

þ ∑ f  ðC Jp Þ2 γ i gT i ¼ u€ 0 and also " # " # ∂2 T 0 1 ∂T 0 € þ ðεn ÞJ ∂u€ 0 þ u€ 0 ðC JT Þ2 þ T ¼ 0 ∂r r 0 ∂r 2 r 0 ∂r

where N

2

hi 2 w ðhi Þ i¼1 2

2

β0 ¼ Bu1 ∑ hi w2 ðhi Þ  Bu2 ∑

i¼1 2 h u 2 2 βi ¼ B1 hi w ðhi Þ Bu2 i w2 ðhi Þ N

β0 ¼ ∑ βi

(

ð51Þ

n J

þ ðε Þ

i¼1

The temperature derivatives can be obtained using the similar method. N ∂T 0 ¼  γ0 T 0 þ ∑ γi T i ∂r i¼1



ð52Þ

where N

2

γ 0 ¼ AT1 ∑ hi w2 ðhi Þ  AT2 ∑ i¼1

hi 2 w ðhi Þ 2

¼ T€ 0

) N ðεn ÞJ € € €  α0 u0 þ ∑ αi ui þ u0 r0 i¼1

ð61Þ

   N 1 1 T 0 þ ∑ ðC JT Þ2 ψ i þ ðC JT Þ2 γ i T i ¼ T€ 0  ψ 0 ðC JT Þ2  γ 0 ðC JT Þ2 r0 r0 i¼1 ( ) n J N ðε Þ þ  ðεn ÞJ α0 þ u€ 0 þ ∑ fðεn ÞJ αi gu€ i r0 i¼1

ð53Þ

i¼1

fϕgT ¼

For second derivative of temperature, it is N

¼  ψ 0T0 þ ∑ ψ iT i

∂r 2

!

ð62Þ

ð63Þ

where

N

γ0 ¼ ∑ γi

∂2 T 0

N 1  γ0 T 0 þ ∑ γi T i r0 i¼1

€ ðN þ 1Þn1 þ ½KðN þ 1ÞnðN þ 1Þ fϕgðN þ 1Þn1 ¼ ½f ðN þ 1Þn1 ½MðN þ 1ÞnðN þ 1Þ fϕg

2

2

γ i ¼ AT1 hi w2 ðhi Þ AT2

i¼1

ð60Þ

The following system of linear equations is obtained for the distributed nodes on the analyzed domain.

hi 2 w ðhi Þ i¼1 2

2

!

N

ðC JT Þ2  ψ 0 T 0 þ ∑ ψ i T i þ ðC JT Þ2

2

N

ð59Þ

i¼1

ð50Þ

i¼1

N

) i 1 h J 2 J 2 1  ðC p Þ  ðC s Þ 2 u0 r0 r0   n o N 1 þ ∑ ðC Jp Þ2 βi þ ðC Jp Þ2 αi ui þ ðC Jp Þ2 γ 0 T 0 r 0 i¼1

 β0 ðC Jp Þ2  α0 ðC Jp Þ2

n

u0

T0

u1

T1

:

:

:

uN

TN

u€ 0

T€ 0

u€ 1

T€ 1

:

:

:

u€ N

T€ N

oT

ð64Þ

and ð54Þ

i¼1

€ T¼ fϕg



T ð65Þ

where N

N

2

hi 2 w ðhi Þ i¼1 2

2

ψ 0 ¼ BT1 ∑ hi w2 ðhi Þ  BT2 ∑ i¼1

5. Time domain analysis

2 h 2 ψ i ¼ BT1 hi w2 ðhi Þ BT2 i w2 ðhi Þ

2

N

ψ0 ¼ ∑ ψi

ð55Þ

i¼1

Also, the second derivative of radial displacement with respect to time can be approximated for first derivative with respect to radius as follows: N ∂u€ 0 ¼ α0 u0 þ ∑ αi ui ∂r i¼1

ð56Þ

where the terms “α0 ” and “αi ” were introduced in Eq. (49). By substituting the obtained relations in star forms for first and second derivatives in governing Eqs. (18) and (19) at a center node, the coupled thermo-elasticity governing equations can be obtained in new form based on the GFD method. In other words, the governing equations should be valid at every center node located on each sub-cylinder for example Jth layer on layered FG cylinder. iu ∂ 2 u0 1 ∂u0 h J 2 ∂T 0 € 0 ðC Jp Þ2 2 þ ðC Jp Þ2  ðC p Þ  ðC Js Þ2 2  ðC Jp Þ2 ¼ u0 ð57Þ r 0 ∂r ∂r ∂r r0 ðC Jp Þ2 h

!

N

 β 0 u0 þ ∑ β i ui i¼1

 ðC Jp Þ2  ðC Js Þ2

iu

0

r0 2

! N 1  α0 u0 þ ∑ αi ui r0 i¼1 ! N  γ u0 þ ∑ γ ui ¼ u€ 0 0

i¼1

i

tp

€ g þ ½Kfϕtp g ¼ ½f p  ½Mfϕ t

ð66Þ 0

0

€ g ¼ ff g  ½Kfϕ0  ½Mfϕ 0

The matrices ½K m and ½K m  ¼ ½K þ n

ð67Þ t ff mp g

are defined as follows:

1 ½M λ1 Δt 2

o n o t t f mp ¼ f p þ

ð68Þ

 1 ½M ð ϕt p  1 λ1 Δt 2

t

tp  1

gÞ

ð69Þ

tp

_ , and ½ϕ€  can be computed using the The matrices of ½ϕtp , ½ϕ following equations: tp

ð58Þ

0

Using the initial conditions ff g and fϕ g, the following equation can be obtained:

þ Δtfϕ_ p  1 g þ ð0:5  λ1 ÞΔt 2 fϕ€

þ ðC Jp Þ2

 ðC Jp Þ2

There are a number of numerical methods to solve the system differential Eq. (63) resulting from the application of the GFD method. In this article, the Newmark time approximation scheme with suitable time step is used, and the time-dependent temperature and displacement fields are obtained for the cylinder. Consider the system to be expressed in terms of nondimensional time t ¼ t p in which the governing equation of system takes the form

fϕtp g ¼ ½K m   1 ff mp g t

ð70Þ

S.M. Hosseini / Engineering Analysis with Boundary Elements 43 (2014) 56–66

n

€ tp ϕ

o

¼

n t o n t o 1  t p  t p  1 € p1 ϕ  ϕ  Δt ϕ_ p  1  Δt 2 ð0:5  λ1 Þ ϕ 2 λ1 Δt ð71Þ

_ g ¼ fϕ _ p  1 g þ Δt ½ð1  λ2 Þfϕ€ fϕ tp

t

tp  1

tp

g þ λ2 fϕ€ g

ð72Þ

tp Using aforementioned equations, the matrices of ½ϕt p , ½ϕ_ , tp €  can be obtained for an arbitrary time. The best converand ½ϕ gence rate can be achieved in this method by choosing λ1 ¼ 14 and λ2 ¼ 12.

61

uðr out ; tÞ ¼ 0

ð78Þ

It is assumed that N nodes are distributed on thickness of cylinder through the radial direction. The first node (i ¼ 1) is located on inner surface and the last one (i ¼ N) is located on outer surface of cylinder. The following constitutive properties are assumed for FG thick hollow cylinder. C P ðr in Þ ¼ 0:7; C T ðr in Þ ¼ 1:2; C S ðr in Þ ¼ 0:28; εn ðr in Þ ¼ 0:08 C P ðr out Þ ¼ 0:5; C T ðr out Þ ¼ 1; C S ðr out Þ ¼ 0:26; εn ðr out Þ ¼ 0:07 ð79Þ

6. Numerical example and discussions To study the thermoelastic wave propagation in thick hollow cylinder, an axisymmetric hollow cylinder in plane strain conditions is considered with inner and outer radii r in andr out , respectively. Required nodes are distributed regularly through the radial directions, which is considered to be N r ¼ 50. The number of star nodes is assumed to cover all distributed nodes. The term “H” is defined as “H ¼ r out r in ”. The following well-known thermal shock loading is considered as the first example. The inner surface of cylinder is assumed to be under sudden heat flux, the outer surface is isolated and radial displacement of nodes on both inner and outer surfaces are considered to be prescribed as follows: q_ ðr in ; tÞ ¼ HðtÞ; ∂Tðr out ; tÞ ¼ 0; ∂r

uðr in ; tÞ ¼ 0

ð73Þ

uðr out ; tÞ ¼ 0

ð74Þ

where q_ is the heat flux rate that is defined in the GN theory of coupled thermo-elasticity as q_ ¼  C 2T ∂T and the term “HðtÞ” is ∂r Heaviside step function. The thermal boundary conditions can be also written in the form of generalized finite difference (GFD) as follows: ! N ∂Tðr in ; tÞ _ in ; tÞ ¼  C T 2 ðr in Þ ¼  C T 2 ðr in Þ  γ 1 T 1 þ ∑ γ i T i ¼  HðtÞ qðr ∂r i¼2 ð75Þ uðr in ; tÞ ¼ 0

ð76Þ

N 1 ∂T ðr out ; tÞ ¼  γ N T N þ ∑ γ i T i ¼ 0 ∂r i¼1

ð77Þ

To verify the present computational method, the volume fraction exponent is assumed to be zero as nr ¼ 0. Consequently, the FG cylinder becomes a homogenous isotropic thick hollow cylinder with following properties, which are the same with Ref. [11] and Ref. [17]. C p ¼ 0:5;

C T ¼ 1:0;

C s ¼ 0:26;

εn ¼ 0:073

Using the above properties (80), the propagation of thermal wave are obtained and compared with those obtained by using the finite element (FE) method [9] and meshless local Petrov–Galerkin (MLPG) method [17]. The thermal wave front is observed using three methods in Fig. 2. Also, the non-dimensional radial displacement wave propagation is shown in Fig. 3 at various times using FE, MLPG and presented hybrid mesh-free methods. A very good agreement can be observed in both Figs. 2 and 3 that means the presented hybrid mesh-free method has a high capability and accuracy to solve the coupled thermoelasticity problems. The influences of volume fraction exponent on radial displacement and temperature wave fronts can be observed in Figs. 4 and 5. It is evidence that the wave fronts in FGM with big value of volume fraction exponent propagate faster than those with small value of exponent. It means that the mean value of thermoelastic wave propagation speed for grading patterns with big values of n r is bigger than mean value of speed for small values of nr . To better assess of the presented method, some results in form of tables are presented for temperature and displacement fields of first example in Tables 1 and 2, respectively. To show more ability of the presented hybrid mesh-free method in coupled thermo-elasticity case, the following thermal shock loading is considered as second example. The inner surface of cylinder is assumed to be under sudden temperature rise and simply supported conditions for displacements as follows: Tðr in ; tÞ ¼ HðtÞ;

uðr in ; tÞ ¼ 0

ð81Þ

0.4 Present method

0.35

Non-dimensional temperature

MLPG method [17]

0.3

LTL method [11]

0.25 Non-dimensional time t = 0.4

0.2

Non-dimensional time t = 0.3

0.15 Non-dimensional time t = 0.2

0.1 0.05 0 −0.05

1

1.05

1.1

1.15

1.2

ð80Þ

1.25

1.3

1.35

1.4

1.45

1.5

Non-dimensional radius Fig. 2. Thermal wave propagation along the thickness of cylinder at various times using three solution techniques for the first example.

62

S.M. Hosseini / Engineering Analysis with Boundary Elements 43 (2014) 56–66

10

x 10−3

Non-dimensional radial displacement

Present method

8

MLPG method [17] LTL method [11]

6 Non-dimensional time t = 0.4

4

Non-dimensional time t = 0.3 Non-dimensional time t = 0.2

2

0

−2

1

1.05

1.1

1.15

1.2

1.25

1.3

1.35

1.4

1.45

1.5

Non-dimensional radius Fig. 3. Non-dimensional radial displacement wave propagation along the thickness of cylinder at various times using three solution techniques for the first example.

0.6 Volume fraction exponent nr = 0.5

Non-dimensional temperature

0.5

Volume fraction exponent nr = 1 Volume fraction exponent nr = 2

0.4

Non-dimensional time t = 0.4

0.3

Non-dimensional time t = 0.3

0.2 Non-dimensional time t = 0.2

0.1 0 −0.1

1

1.05

1.1

1.15

1.2

1.25

1.3

1.35

1.4

1.45

1.5

Non-dimensional radius Fig. 4. Thermal wave propagation along the thickness of cylinder for various grading patterns of FGM in the first example.

20

x 10−3 Volume fraction exponent nr = 0.5

Non-dimensional radial displacement

Volume fraction exponent nr = 1

15

Volume fraction exponent nr = 2

10

Non-dimensional time t = 0.4

Non-dimensional time t = 0.3

5

Non- dimensional time t = 0.2

0

−5

1

1.05

1.1

1.15

1.2

1.25

1.3

1.35

1.4

1.45

1.5

Non-dimensional radius Fig. 5. Radial displacement wave propagation along the thickness of cylinder for various grading patterns of FGM in the first example.

S.M. Hosseini / Engineering Analysis with Boundary Elements 43 (2014) 56–66

63

Table 1 Some results of temperature field in the first example. t ¼ 0:2

nr nr nr nr

¼ 0:1 ¼ 0:5 ¼1 ¼3

t ¼ 0:3

t ¼ 0:4

r ¼ 1:1

r ¼ 1:25

r ¼ 1:4

r ¼ 1:1

r ¼ 1:25

r ¼ 1:4

r ¼ 1:1

r ¼ 1:25

r ¼ 1:4

0.1154 0.1532 0.1721 0.1848

0.0016 0.0088 0.0212 0.0419

0 0 0 0

0.2091 0.2564 0.2822 0.3075

0.0697 0.1117 0.1369 0.1642

0 0.0022 0.0056 0.0217

0.2955 0.3493 0.3806 0.4178

0.1543 0.206 0.2353 0.2721

0.0182 0.0675 0.0966 0.1324

Table 2 Some results of displacement field in the first example. t ¼ 0:2

nr nr nr nr

¼ 0:1 ¼ 0:5 ¼1 ¼3

t ¼ 0:3

t ¼ 0:4

r ¼ 1:1

r ¼ 1:25

r ¼ 1:4

r ¼ 1:1

r ¼ 1:25

r ¼ 1:4

r ¼ 1:1

r ¼ 1:25

r ¼ 1:4

0.0023 0.0039 0.0046 0.005

0 0 0 0

0 0 0 0

0.0055 0.0075 0.0086 0.0095

0 0.0024 0.004 0.0064

0 0 0 0

0.0086 0.011 0.0125 0.0139

0.0047 0.0095 0.0131 0.0174

0 0 0.0017 0.0041

1.4 Volume fraction exponent nr = 0.5

Non-dimensional temperature

1.2

Volume fraction exponent nr = 1 Volume fraction exponent nr = 2

1 0.8 Non-dimensional time t = 0.4

0.6

Non-dimensional time t = 0.3

0.4 Non-dimensional time t = 0.2

0.2 0

1

1.05

1.1

1.15

1.2

1.25

1.3

1.35

1.4

1.45

1.5

Non-dimensional radius Fig. 6. Thermal wave propagation along the thickness of cylinder for various grading patterns of FGM in the second example.

Tðr out ; tÞ ¼ 0;

uðr out ; tÞ ¼ 0

ð82Þ

where HðtÞ is the Heaviside unit step function. The propagations of radial displacement and thermal waves along the radial direction of FG cylinder are illustrated in Figs. 6 and 7, respectively. The variation in values of volume fraction exponent effects on propagation of wave fronts is clear. Similar to the previous example, it is concluded from both Figs. 6 and 7 that the mean values of thermoelastic wave propagation speed for FGM with big values of n r are bigger than those with small values of n r . Also, various grading patterns of mechanical properties in FG cylinder effect directly on time histories of radial displacement and temperature fields. In Figs. 8 and 9, the time histories of radial displacement and temperature for middle point on thickness of FG cylinder are drawn, respectively. The frequency of variation in time histories diagrams is decreased when the exponent value of volume fraction is increased. This phenomenon depends to thermoelastic wave propagation speed when the mean value of thermoelastic wave propagation is increased, the frequency of variation in time histories diagrams are increased. The time histories of radial displacement and temperature fields for various points on

thickness of FG cylinder are presented in Figs. 10 and 11 for nr ¼ 0:5. The points, which are located close to inner surface, start to oscillate faster than points located far from inner surface. This is because the thermoelastic wave starts to propagate from inner surface of FG cylinder, which are excited by a thermal shock loading. The effects of number of nodes for regular distribution of nodes through the radial direction on temperature field can be observed for the first example in Fig. 12. It can be concluded that the diagrams for values of N ¼ 20 are smoother than diagrams for small ones. Also, there is no any significant difference between diagrams for N ¼ 20 and N ¼ 15. It means that the optimum values of N should be more than 20. In other word, the selected value as N ¼ 50 is effective value with high convergence rate. A similar study can be carried out for displacement field and it is concluded that the selected value as N ¼ 50 is proper for analysis. By using the random irregular distribution of nodes, the number of nodes should be increased to achieve the same convergence rate and stability by regular distribution. The effects of number of nodes on stability of results in temperature field are presented in Fig. 13. As it is shown in Fig. 13, the number of nodes

64

S.M. Hosseini / Engineering Analysis with Boundary Elements 43 (2014) 56–66

0.09 Volume fraction exponent nr = 0.5

Non-dimensional radial displacement

0.08

Volume fraction exponent nr = 1 Non-dimensional time t = 0.3

Volume fraction exponent nr = 2

0.07 0.06

Non-dimensional time t = 0.4

0.05 0.04 0.03 0.02 0.01 0

Non-dimensional time t = 0.2

1

1.05

1.1

1.15

1.2

1.25

1.3

1.35

1.4

1.45

1.5

Non-dimensional radial distance Fig. 7. Radial displacement wave propagation along the thickness of cylinder for various grading patterns of FGM in the second example.

0.25 Volume fraction exponent nr=0.5 r = rin + H/2

Non-dimensional radial displacement

0.2

Volume fraction exponent nr=1 Volume fraction exponent nr=2

0.15 0.1 0.05 0 -0.05 -0.1

0

1

2

3

4

5

6

Non-dimensional time Fig. 8. Time history of radial displacement for various grading patterns of FGM in the second example.

r = rin + H/2

Non-dimensional temperature

1 0.8 0.6 0.4 0.2 0 Volume fraction exponent nr = 2

-0.2

Volume fraction exponent nr = 1 Volume fraction exponent nr = 0.5

-0.4 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

Non-dimensional time Fig. 9. Time history of temperature for various grading patterns of FGM in the second example.

7. Conclusion

the Green–Naghdi theory (without energy dissipation). The thick hollow cylinder is considered as the analyzed domain for the problem, which is under thermal shock loading. The axi-symmetry and plane strain conditions are assumed in this article. The main outputs of this paper can be outlined as follows:

The application of a hybrid mesh-free method based on the generalized finite difference (GFD) method has been developed to solve the coupled thermo-elasticity governing equations based on

 The governing coupled thermo-elasticity equations of thick hollow cylinder are derived in new mesh-free forms. The propagations of thermal and non-dimensional radial displacement

should be considered more than 150 to achieve stable results in temperature field.

S.M. Hosseini / Engineering Analysis with Boundary Elements 43 (2014) 56–66

65

0.25 r = r + H/4

Non-dimensional radial displacement

Volume fraction exponent nr = 0.5

in

r = r + H/2

0.2

in

r = r in + 3H/4

0.15 0.1 0.05 0 -0.05 -0.1 -0.15

0

1

2

3

4

5

6

Non-dimensional time Fig. 10. Time history of radial displacement for several points located on the thickness of FG cylinder in the second example.

Non-dimensional temperature

r = r + H/4 in

Volume fraction exponent nr = 0.5

1.2

r = r in + H/2 r = r + 3H/4 in

1 0.8 0.6 0.4 0.2 0 -0.2 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

Non-dimensional time Fig. 11. Time history of temperature for several points located on the thickness of FG cylinder in the second example.

0.3

Non-dimensional temperature

0.25

Number of nodes = 5 Number of nodes = 10 Number of nodes = 15 Number of nodes = 20 and more

Non-dimensional time t = 0.2

0.2 0.15 0.1 0.05 0 -0.05

1

1.05

1.1

1.15

1.2

1.25

1.3

1.35

1.4

1.45

1.5

Non-dimensional radius Fig. 12. The effects of number of nodes on stability of temperature field for regular distribution of nodes.

wave through the radial direction of cylinder are simulated and discussed in details for two kinds of thermal shock loadings.  Employing the presented hybrid mesh-free method, the propagation thermoelastic wave fronts are observed and it is

discussed in detail for various kinds of grading patterns in FGMs.  The influence of various grading patterns in FGMs is studied in time histories of radial displacement and temperature fields.

66

S.M. Hosseini / Engineering Analysis with Boundary Elements 43 (2014) 56–66

0.3

Non-dimensional temperature

0.25

Number of nodes = 25 Number of nodes = 50 Number of nodes = 75 Number of nodes = 100 Number of nodes = 125 Number of nodes = more than 150

Non-dimensional time t = 0.2

0.2 0.15 0.1 0.05 0 -0.05

1

1.05

1.1

1.15

1.2

1.25

1.3

1.35

1.4

1.45

1.5

Non-dimensional radius Fig. 13. The effects of number of nodes on stability of temperature field for random irregular distribution of nodes.

The presented hybrid mesh-free method is successfully used to achieve this purpose.  The proposed method is a truly mesh-free method, which requires neither domain elements nor background cells in either the interpolation or the integration.  The mathematical formulations and also the obtained results show that the presented hybrid mesh-free method can be used to solve the coupled system of differential equations. In general, the paper develops the application of a hybrid mesh-free method based on the GFD in coupled thermoelasticity and thermoelastic wave propagation analysis in FG thick hollow cylinder. The obtained results confirm that the presented method in this article is an effective method for analysis of heat conduction, thermal sciences and also coupled systems.

Acknowledgments This work is supported by the Ferdowsi University of Mashhad as a research project with No. 2/22445 and date (10-07-2012). References [1] Carrera E. A class of two-dimensional theories for multilayered plates analysis. Accademia delle Scienze Torino, Memorie Scienze Fisiche 1995:1–39. [2] Carrera E. Theories and finite elements for multilayered plates and shells: a unified compact formulation with numerical assessments and benchmarking. Arch Comput Methods Eng 2003;10(3):215–96. [3] Brischetto S, Carrera E. Coupled thermo-mechanical analysis of one-layered and multilayered plates. Compos Struct 2010;92:1793–812. [4] Nakonieczny K, Sadowski T. Modelling of ‘thermal shocks’ in composite materials using a meshfree FEM. Comput Mater Sci 2009;44:1307–11. [5] Youssef HM. Two-temperature generalized thermoelastic infinite medium with cylindrical cavity subjected to moving heat source. Arch Appl Mech 2010;80:1213–24. [6] Mukhopadhyay S, Kumar R. State-space approach to thermoelastic interactions in generalized thermoelasticity type III. Arch Appl Mech 2010;80: 869–81. [7] Hosseini SM, Abolbashari MH. Analytical solution for thermoelastic waves propagation analysis in thick hollow cylinder based on Green–Naghdi model of coupled thermoelasticity. J Therm Stresses 2012;35:363–76. [8] Hosseini Zad SK, Komeili A, Eslami MR, Fariborz S. Classical and generalized coupled thermoelasticity analysis in one-dimensional layered media. Arch Appl Mech 2012;82:267–82.

[9] Hosseini SM. Coupled thermoelasticity and second sound in finite length functionally graded thick hollow cylinders (without energy dissipation). Mater Des 2009;30:2011–23. [10] Hosseini SM, Akhlaghi M, Shakeri M. Heat conduction and heat wave propagation in functionally graded thick hollow cylinder base on coupled thermoelasticity without energy dissipation. Heat Mass Transfer 2008;44: 1477–84. [11] Taheri H, Fariborz S, Eslami MR. Thermoelastic analysis of an annulus using the Green–Naghdi model. J Therm Stresses 2005;28(9):911–27. [12] Veres IA, Berer T, Burgholzer P. Numerical modeling of thermoelastic generation of ultrasound by laser irradiation in the coupled thermoelasticity. Ultrasonics 2013;53:141–9. [13] Sladek J, Sladek V, Zhang CH. Application of meshless local Petrov–Galerkin (MLPG) method to elastodynamic problems in continuously nonhomogeneous solids. CMES: Comput Model Eng Sci 2003;4:637–48. [14] Sladek J, Sladek V, Hon YC. Inverse heat conduction problems by meshless local Petrov–Galerkin method. Eng Anal Bound Elem 2005;30:650–61. [15] Sladek J, Sladek V, Tan CL, Atluri SN. Analysis of transient heat conduction in 3D anisotropic functionally graded solids, by the MLPG method. CMES 2008;32(3):161–74. [16] Sladek J, Sladek V, Ch Zhang, Tan CL. Linear coupled thermoelastic analysis for 2-d orthotropic solids by MLPG. ICCES 2007;3(2):87–92. [17] Hosseini SM, Sladek J, Sladek V. Meshless local Petrov–Galerkin method for coupled thermoelasticity analysis of a functionally graded thick hollow cylinder. Eng Anal Bound Elem 2011;35(6):827–35. [18] Hosseini SM, Sladek J, Sladek V. Application of meshless local integral equations to two dimensional analysis of coupled non-Fick diffusion-elasticity. Eng Anal Bound Elem 2013;37:603–15. [19] Benito JJ, Ureña F, Gavete L. Influence of several factors in the generalized finite difference method. Appl Math Model 2001;25(12):1039–53. [20] Benito JJ, Ureña F, Gavete L. An h-adaptive method in the generalized finite differences. Comput Methods Appl Mech Eng 2003;192:735–9. [21] Gavete L, Benito JJ, Gavete ML. Improvements of generalized finite difference method and comparison with other meshless method. Appl Math Model 2003;27(10):831–47. [22] Benito JJ, Ureña F, Gavete L. Solving parabolic and hyperbolic equations by the generalized finite difference method. J Comput Appl Math 2007;209:208–33. [23] Ureña F, Benito JJ, Salete E, Gavete L. A note on the application of the generalized finite difference method to seismic wave propagation in 2D. J Comput Appl Math 2012;236(12):3016–25. [24] Gavete L, Ureña F, Benito JJ, Salete E. A note on the dynamic analysis using the generalized finite difference method. J Comput Appl Math 2012. http://dx.doi. org/10.1016/j.cam.2012.06.035. [25] Hosseini SM. Analysis of elastic wave propagation in a functionally graded thick hollow cylinder using a hybrid mesh-free method. Eng Anal Bound Elem 2012;36:1536–45. [26] Hosseini SM. Shock-induced thermoelastic wave propagation analysis in a thick hollow cylinder without energy dissipation using mesh-free generalized finite difference (GFD) method. Acta Mech 2013;224(3):465–78. [27] Green AE, Naghdi PM. Thermoelasticity without energy dissipation. J Elast 1993;31:189–208.