Engineering Analysis with Boundary Elements 60 (2015) 154–161
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements
journal homepage: www.elsevier.com/locate/enganabound
A novel meshless local Petrov–Galerkin method for dynamic coupled thermoelasticity analysis under thermal and mechanical shock loading
Bao-Jing Zheng a,b, Xiao-Wei Gao a,c,n, Kai Yang a, Chuan-Zeng Zhang b
a School of Aeronautics and Astronautics, Dalian University of Technology, Dalian 116024, PR China
b Department of Civil Engineering, University of Siegen, D-57068 Siegen, Germany
c State Key Laboratory of Structural Analysis for Industrial Equipment, Dalian University of Technology, Dalian 116024, PR China
article info
Article history:
Received 30 September 2014 Received in revised form
15 December 2014
Accepted 15 December 2014 Available online 6 January 2015
Keywords:
Dynamic coupled thermoelasticity Thermal and mechanical shock Meshless local Petrov–Galerkin method Moving Kriging interpolation method Backward difference method
1. Introduction
Coupled thermoelasticity is relevant to many problems in engi- neering and sciences. It encompasses the phenomena that describe the elastic and thermal behavior of thermoelastic solids and their interactions under mechanical and thermal loadings. In recently years, considerable attention has been paid to analyze the coupled thermoelasticity problems since the thermal stresses play an impor- tant role in the structures. Analytical solutions for some dynamic problems in uncoupled thermoelasticity were obtained by Danilovs- kaya [1,2]. For some quasi-static problems in coupled thermoelas- ticity, Boley and Tolins [3] obtained the analytical solutions. The class of problems which has analytical solutions in dynamic thermoelas- ticity is extremely limited. Therefore several computational methods have been proposed for the numerical analysis of the thermoelas- ticity problems with complex geometries and boundary conditions. Many of them have been directed to uncoupled problems in steady or transient states. Few investigations have been done successfully for coupled thermoelasticity.
Numerical solutions for the coupled thermoelasticity using the finite element method (FEM) have been reported by Keramidas and
n Corresponding author at: School of Aeronautics and Astronautics, Dalian University of Technology, Dalian, 116024, P.R. China. Tel.: +86 411 84706332;
fax: +86 411 84706332.
E-mail address: xwgao@dlut.edu.cn (X.-W. Gao).
http://dx.doi.org/10.1016/j.enganabound.2014.12.001
0955-7997/& 2015 Published by Elsevier Ltd.
abstract
A meshless local Petrov–Galerkin method (MLPG) based on the moving Kriging interpolation is further developed for two-dimensional linear dynamic coupled thermoelasticity problems under thermal and mechanical shock loading. In this method the moving Kriging interpolation is employed to construct the shape functions, and the Heaviside step function is applied in the local weak-form as the test functions. The equations of motion and transient heat conduction equations of the coupled thermoelasticity interact on each other. So these equations must be solved simultaneously. In this paper the backward difference method is applied to approximate the time derivatives. The main idea of treating the time- dependent problem is based on the reduction procedure of the equations of motion, in which the displacements, velocities and temperature are taken as the field variables. Then a unified explicit expression of the system of differential algebraic equations is obtained. Numerical examples are presented to illustrate the accuracy of this method.
& 2015 Published by Elsevier Ltd.
Ting [4], Ting and Chen [5], Prevost and Tao [6], Tamma and Railkar [7], Cannarozzi and Ubertini [8]. The boundary element method (BEM) has been developed and applied to analyze thermoelasticity problems by Rizzo and Shippy [9]. The radial-integration BEM for 2-D and 3-D uncoupled thermoelasticity analysis was presented by Gao [10]. The BEM also has been successfully applied to coupled thermoelasticity problems [11–17].
In recent years, meshless or element-free methods have been developed as alternative numerical approaches with efforts to eliminate some known shortcomings of the mesh-based methods [18–23]. The main objective of these methods is to reduce the difficulties of meshing and remeshing of complex structural problem domains. The solution is entirely built in terms of a set of distributed nodes, thus no element connectivity is required. The meshless local Petrov–Galerkin method (MLPG) is one of the most widely used meshless methods, as it does not need a mesh, either for the purpose of interpolation of the problem variables, or for the integration of the weak-form. A variety of meshless methods have been developed and applied to thermoelastic problems [24–31].
Most of the methods mentioned above for transient coupled thermoelasticity were based on the Laplace-transform technique to treat the time domain. However, the Laplace-transform of some boundary conditions may not exist. In this paper, the MLPG method based on the moving Kriging interpolation and with a Heaviside step function as the test functions [32–37] is further applied to solve two- dimensional dynamic coupled thermoelasticity problems. Since the
moving Kriging shape function possess the delta function property, the present method can easily enforce the essential boundary conditions as done in FEM. Then, the backward difference method is used to approximate the time derivatives instead of the Laplace-transform technique. For the dynamic coupled thermoelasticity, the equations of motion involve the second derivative with respect to the time, while the heat conduction equation has only the first time derivative. In order to use the backward difference method to discretize the time derivatives, velocities as extra field variables are used to obtain the reduction procedure for the equations of motion. A unified explicit expression of the system of differential algebraic equations is obtained. The accuracy and the efficiency of the present MPLG method are verified by two numerical examples.
2. Moving Kriging interpolation
Similar to the moving least-squares (MLS) approximation, the moving Kriging approach can be applied to any sub-domain Ωx D Ω. According to Ref. [34], the field variable uðxÞ in the problem domain Ωcan be approximated by uhðxÞ. For any sub-domain, the local approximation can be defined as follows:
[33]. In addition, the correlation matrix R1⁄2Rðxi;xjÞ is also given in an explicit form
B.-J. Zheng et al. / Engineering Analysis with Boundary Elements 60 (2015) 154–161 155
R1⁄2Rðxi;xjÞ1⁄46 ⋮ 4
⋮ ⋱⋮7 ð11Þ 5
2 1 Rðx1;x2Þ ⋯ Rðx1;xnÞ3 6Rðx2;x1Þ 1 ⋯ Rðx2;xnÞ7
Rðxn; x2Þ ⋯1
The partial derivatives of the shape functions can be easily
Rðxn; x1Þ obtained as follows:
φk;xl 1⁄4 ∂φkðxÞ 1⁄4 ∑m ∂pjðxÞAjk þ ∑n ∂riðxÞBik ð12Þ ∂xl j∂xl i∂xl
An important property of the moving Kriging shape function is
that the shape functions possess the delta function property [34],
namely
(
1 ðk1⁄4j;k;j1⁄41;2;⋯;nÞ
φkðxjÞ1⁄4 0 ðkaj;k;j1⁄41;2;⋯;nÞ ð13Þ
Another important property of the moving Kriging shape functions is the consistency property [34]
hTT
n
∑ φkðxÞ1⁄41
k1⁄41
ð14aÞ ð14bÞ
u ðxÞ 1⁄4 1⁄2p ðxÞAþr ðxÞBu ð1Þ
or n
n
uh ðxÞ 1⁄4 ∑ φk ðxÞuk ð2Þ
k
where the moving Kriging shape function φkðxÞ is defined by the following equation:
mn
φkðxÞ 1⁄4 ∑ pjðxÞAjk þ ∑ riBik ð3Þ
ji
The matrices A and B are given by the following equations: A1⁄4ðPTR1PÞ1PTR1 ð4Þ
B1⁄4R1ðIPAÞ ð5Þ where I is an identity matrix and the vector pðxÞ is the polynomial
with m basis functions
pðxÞ 1⁄4 ðp1ðxÞ; p2ðxÞ ;⋯ ; pmðxÞÞT ð6Þ
For the matrix P with the size nm, the values of the polynomial basis functions (6) at the given set of nodes are collected, i.e.,
∑ φkðxÞxk1⁄4x k1⁄41
3. MLPG discretization
3.1. Governing equations
26 p1ðx1Þ
p2ðx1Þ ⋯
pmðx1Þ 37
pmðx2Þ 7 ð7Þ
rðxÞ 1⁄4 ðRðx1 ; xÞ; Rðx2 ; xÞ; ⋯; Rðxn ; xÞÞT ð8Þ
Here, Rðxi;xjÞ is the correlation function between any pair of nodes located at xi and xj. Although many functions can be used as the correlation function Rðxi;xjÞ, however, a simple and frequently used correlation function is the Gaussian function
Rðxi;xjÞ1⁄4expðξr2ijÞ ð9Þ
in which
rij 1⁄4‖xixj‖=h ð10Þ
and h is the average distance of the nodes in the support domain, ξ represents the value of the correlation parameter used to fit the model, and it is a good choice to take ξ 1⁄4 0:03 0:2 according to
p2ðx2Þ ⋯
64 ⋮ ⋮ ⋱ ⋮ 75 p1ðxnÞ p2ðxnÞ ⋯ pmðxnÞ
P 1⁄4 6 p1ðx2Þ
and the vector rðxÞ in Eq. (1) is also given by the following
equation:
We consider the dynamic coupled thermoelasticity problem of a homogeneous, isotropic and linear elastic heat conductor based on the classical thermoelasticity theory. Let Ω be a finite domain with boundary Γ, and T1⁄41⁄20;t be the time interval. The basic equations of linear coupled thermoelasticity can be written as follows:
εij 1⁄42ðui;jþuj;iÞ ð15Þ Constitutive equation
σij 1⁄4 2μεkl þλεkkδij βθδij ð16Þ Equation of motion
σij;jþbiρu€i 1⁄40 ð17Þ Heat conduction equation
kθ;ii ρcθ_ βθ0u_i;i þQ 1⁄4 0 ð18Þ
In Eqs. (15)–(18), εij, ui, σij, β, θ, θ0, ρ, c, k, bi and Q are the strain tensor, displacement vector, Cauchy stress tensor, stress- temperature modulus, temperature difference, reference tempera- ture, mass density, specific heat, thermal conductivity, body force vector and heat sources, respectively, and λ and μ are Lame’s constants. The corresponding boundary and initial conditions are given as follows:
Boundary conditions
ui 1⁄4ui onΓu T ti 1⁄4σijnj 1⁄4ti onΓt T
Kinematic relation 1
156
B.-J. Zheng et al. / Engineering Analysis with Boundary Elements 60 (2015) 154–161
θ1⁄4θ onΓθ T Initial conditions
the time. They can be expressed as follows:
q1⁄4kθ;ini 1⁄4q onΓq T; uiðx;0Þ1⁄4ui0ðxÞ u_iðx;0Þ1⁄4u_i0ðxÞ θðx;0Þ1⁄4θ0ðxÞ inΩ:
θðx; tÞ 1⁄4 uðx; tÞ 1⁄4
n
∑ φkðxÞθkðtÞ ð25Þ
k1⁄41
n
∑ φkðxÞukðtÞ ð26Þ
k1⁄41
where Γu is the part of the global boundary with prescribed displacements, while on Γt, Γθ and Γq the traction vector, the temperature and the heat flux are prescribed, respectively.
3.2. Local weak formulation
The local weak-form is constructed over a sub-domain Ωs bounded by Γs. The local sub-domains overlap each other, and cover the whole problem domain Ω. They could be of any geometric shape and size [20]. For simplicity, the sub-domains are taken to be of a circular shape or rectangular shape. Using the local weighted residual method, the generalized local weak-form of Eq. (17) can be written as follows:
Z
wik ðσij;jþbiρu€iÞdΩ1⁄40 ð19Þ Ωs
where wik is the weight function.
Applying the divergence theorem to Eq. (19), the following
Substituting the displacement and temperature expressions in Eqs. (25) and (26) into the local weak-form (22), the discrete equations for all nodes are given as follows:
Mu u€ðtÞþK uuðtÞ þHuθðtÞ 1⁄4 FuðtÞ ð27Þ
where the submatrices of the mass matrix Mu, the stiffness matrix
Ku, the matrix Hu and the subvector of the load vector Fu can be
βnΦJdΩ ð30Þ
ð31Þ The matrices ΦJ, N, BJ and D, and the vector n in the above
( 2φ0321v03 δik; xAΩs J;x
Using the Heaviside step function as the test function and considering the natural boundary conditions, the local weak-form
(20) can be rewritten as follows: ZZZZZ
ρu€i dΩ tidΓ tidΓ1⁄4 tidΓþ bi dΩ ð22Þ Ωs Γsu Γsi Γst Ωs
where Γsi is a part of the local boundary of Ωs over which no boundary conditions are specified, Γst is the intersection of Γt and the boundary Γs, Γsu is the intersection of Γu and the boundary Γs.
Similarly, the local weak-form of Eq. (18) can be written as follows:
Z
In order to simplify Eq. (20), the test function w
ik
Φ1⁄4 J ;N1⁄4 x y ;n1⁄4ðnx nyÞT; J0φJ 0nynx
ðbi ρ u€iÞ dΩ 1⁄4 0 ð20Þ
equations are defined by the following equation: “φ0#”n0n#
wik1⁄4 0; x2=Ωs ð21Þ BJ1⁄464 0 φJ;y75;D1⁄4 E 64v 1 0 75;forplanestress;
written as follows:
Z
MuIJ1⁄4 ρΦJdΩ Ωs
ZZ
ð28Þ
NDBJdΓ ð29Þ
KuIJ 1⁄4 HuIJ 1⁄4
Γsi βnΦJdΩþ
Γsu
NDBJdΓ ZZ
Γsi ZZ
form can be obtained as follows:
ZZZ
wik σijnjdΓ wik;j σij dΩþ wik Γs Ωs Ωs
is chosen in such a way that it eliminates or simplifies the domain integral over Ωs. This can be accomplished by using the Heaviside step function
Γsu FuI 1⁄4 tdΓþ bdΩ
Γqt Ωq
φ φ
1v2 0 0 1v 2
J;y J;x
21v03 1v
w ðkθ;iiρcθ_βθ0u_i;iþQÞdΩ1⁄40 ð23Þ Ωs
Eð1vÞ6v 1 07
D1⁄4ð1þvÞð12vÞ41v 5; forplanestrain ð32Þ
0 0 12v 2ð1vÞ
Similarly, substituting the displacement and temperature expressions in Eqs. (25) and (26) into the local weak-form (24), the discrete equations for all nodes are given as follows:
Cθ θ_ðtÞþKθ θðtÞþGθu_ðtÞ1⁄4FθðtÞ ð33Þ
where the submatrices of the heat capacity matrix Cθ, the stiffness
matrix Kθ, the matrix Gθ and the subvector of the load vector Fθ
can be written as follows:
Z
CθIJ1⁄4 ρcφJdΩ ð34Þ Ωs
ZZ
KθIJ1⁄4 knTPJdΓ knTPJdΓ ð35Þ
where w is the weight function.
Applying the Gauss divergence theorem to the local weak-form
(23) and considering the Heaviside step function for the weight function w and the natural boundary conditions, the following
Γsθ
βθ0PTJ dΩ
Γsi
form can be obtained:
ZZZZZZ
ρcθ_dΩþ βθ0u_i;idΩ qdΓ qdΓ1⁄4 qdΓþ
Ωs Ωs Γsθ Γsi Γsq Ωs
Z ZZ
Q dΩ ð24Þ
GθIJ 1⁄4
FθI 1⁄4 Q dΩþ
ð36Þ ð37Þ
where Γsθ is the intersection of Γθ and the boundary Γs, Γsq is the intersection of Γq and the boundary Γs.
3.3. Numerical implementation
In the coupled thermoelastodynamics, the displacements and the temperature are functions of both the spatial coordinates and
qdΓ
Ωs
Ωs Γsq in which
PJ 1⁄4ðφJ;x φJ;y ÞT
Eqs. (27) and (33) are formulated for the displacements and the temperature, and they must be solved simultaneously. In order to obtain a unified expression for the system of differential algebraic
B.-J. Zheng et al. / Engineering Analysis with Boundary Elements 60 (2015) 154–161 equations, the velocity terms are used as the field variables.
157
Eq. (27) can be rewritten as follows:
!_! !_
Mu 0 u_ þ 0 Ku u þHuθ1⁄4 Fu ð38Þ 0IuI0u 0
Combining Eqs. (33) and (38) we can obtain a unified expres-
sion as follows:
2 30 1 2 30 1 0 1
Mu 0 0 _ 0 Ku Hu u_ Fu u
640 I 075B@uCAþ64I 0 075B@uCA1⁄4B@0CA ð39Þ 0 0 Cθ θ Gθ 0 Kθ θ Fθ
or
C U_ þ K U 1⁄4 F ð40Þ
Using the backward difference technique for the approximation of the first time derivative, Eq. (40) can be written as follows: 1 1
Cþ2KΔt Unþ1 1⁄4Fnþ1Δtþ C2KΔt Un ð41Þ Eq. (41) is an explicit time-stepping scheme for computing the
vector U at the (nþ1)-th time-step.
4. Numerical examples and discussion
In order to demonstrate the versatility and the accuracy of the present meshless method, first we consider a quasi-static thermo- elasticity problem as a special case, and then a dynamic coupled thermoelasticity problem under several loading cases is consid- ered. The computed results by the present method are compared with that of the finite element method. Square sub-domain and square interpolation domain are used for the numerical imple- mentation of the meshless method, the parameters of the influ- ence domain and the sub-domain are taken as 4.0 and 0.5 respectively, the correlation parameter ξ is taken as 0.1, and two Gaussian points for each direction are used for integration of the domain-integrals.
For convenience we introduce the following dimensionless variables:
x tC1 σij θθ0 λþμ x^1⁄4l;t^1⁄4l;σ^ij1⁄4βθ;θ^1⁄4θ;u^i1⁄4aβui ð42Þ
00
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
where l 1⁄4 k=ρcC1 is the dimensionless length and C1 1⁄4 ðλþ2μÞ=ρ is the velocity of the longitudinal wave. In the following examples the dimensionless variables are used, we drop the hat for convenience.
To study the coupling effect, a measure of the thermoelastic coupling is defined by the dimensionless thermoelastic parameter [8]
θ0 β2
C1⁄4ρcðλþ2μÞ ð43Þ
Fig. 1. A suddenly heated unit square panel.
Fig. 2. Nodal distribution. are given by Carslaw and Jaeger [38]
8> 4 1 ð1Þn ð2nþ1Þπ2κt ð2nþ1Þπy >θðy;tÞ1⁄41π ∑ exp 2 cos
< n1⁄402nþ1 4L 2L
u ðy;tÞ1⁄4ð1þvÞαRy θðy;tÞdy; >2 ð1vÞ0
>:σ11ðy;tÞ1⁄4 αE θðy;tÞ: ð1vÞ
;
where C 1⁄4 0 corresponds to the materials, C ranges from 0.01 to 0.1.
uncoupled
case. For
traditional
where L is the diffusion length and κ is the diffusivity, both are being assigned unit values. The material constants used are α 1⁄4 0:02, v 1⁄4 0:3 and plane strain condition is assumed.
As shown in Fig. 2, 121 (11 11) uniform nodes are used to discretize the problem domain. Analysis is conducted for several non-dimensional time steps in order to investigate convergence characteristics. Numerical results for the temperature θ, the dis- placement u2 and the lateral stress σ11 at the mid-point of the y- axis are presented in Figs. 3–5, respectively. It can be seen that the results obtained by the present method have an excellent agree- ment with the exact solution for the time variation while the time step is small enough.
4.2. Dynamic coupled thermoelasticity problem
An isotropic and homogeneous square plate under three loading cases on the left side is considered in this section. The boundary element solution of the problem is given by Hosseini-Tehrani [15],
4.1. Quasi-static thermoelasticity problem
Consider a unit square panel under a sudden heating on the top side as shown in Fig. 1. The inertial terms ρu€ i are neglected and the reference temperature is assumed to be zero in this problem. The dimensionless thermoelastic parameter C is taken as zero for the uncoupled case. The analytical solutions to uncoupled thermoelasticity
ð44Þ
158 B.-J. Zheng et al. / Engineering Analysis with Boundary Elements 60 (2015) 154–161
Fig. 3. Temporal variation of the temperature θ at the mid-point of the y-axis.
Fig. 4. Temporal variation of the displacement u2 at the mid-point of the y-axis.
Fig. 5. Temporal variation of the lateral stress σ11 at the mid-point of the y-axis.
and the finite element solution is obtained by the FEM software of COMSOL Multiphysics.
The dimensionless size of the square plate l1⁄410 and a unit thickness are taken, as shown in Fig. 6. 169 (13 13) nodes with a regular node distribution as shown in Fig. 7 are used in the numerical computation, and the FEM mesh distribution is shown in Fig. 8. The time steps are taken as Δt 1⁄4 0:02, The thermal and the mechanical boundary conditions are prescribed as follows:
ui 1⁄40;q1⁄40; atx1⁄410;
ti 1⁄40;q1⁄40; aty1⁄475:
where ui, ti and q are the displacements, the traction components and the heat flux density. The left side of the plate is kinematically free and thermally exposed to a temperature or pressure shock described by the known function
fðtÞ1⁄45texpð2tÞ ð45Þ
where t is the dimensionless time. Three loading cases are considered in the following.
Fig. 8. FEM mesh distribution over the square plate.
1) Thermal shock alone: The thermal shock described by Eq. (45) is applied on the left side where ti 1⁄4 0. Figs. 9–11 show the comparisons of the temperature θ, the displacement u1, and the stress σx along the x-axis for the uncoupled (C1⁄40) and coupled (C1⁄40.1) cases at dimensionless time t1⁄43 and t1⁄46, respec- tively. Due to the coupling effect, the temperature distributions
Fig. 6. A square plate subjected to thermal shock loading on the left side.
Fig. 7. Node distribution over the square plate.
Fig. 9. Comparison of the temperature θ along x-axis for thermal shock.
Fig. 12. Comparison of the temperature θ along x-axis for pressure shock loading.
Fig. 10. Comparison of the displacement u1 along x-axis for thermal shock loading.
Fig. 11. Comparison of the stress σx along x-axis for thermal shock loading.
of the uncoupled and coupled cases are significantly different, and the thermal shock generating stress wave is clearly displayed. It is seen that the temperature θ, the displacement u1, and the stress σx results of the present method have good agreement with the finite element solutions for the uncoupled and coupled case.
2) Pressure shock alone: The pressure shock given by Eq. (45) is applied to the left side where q 1⁄4 0. Fig. 12 shows the compar- ison of the temperature distribution along the x-axis for the uncoupled (C1⁄40) and coupled (C1⁄40.1) cases at dimensionless time t 1⁄4 3 and t 1⁄4 6. The temperature is always zero for the uncoupled case. However, even though a pressure shock is applied, the temperature rises due to the mechanical energy being changed into the heat energy. Figs. 13 and 14 show the comparison of the displacement u1 and the stress σx along the x-axis for uncoupled (C 1⁄4 0) to coupled (C 1⁄4 0.1) cases at dimensionless time t 1⁄4 3 and t 1⁄4 6, respectively. The results of the present method agree well with the finite element solu- tions for the uncoupled and coupled cases.
Fig. 14. Comparison of the stress σx along x-axis for pressure shock loading.
3) Thermal and pressure shock in combination: The combined thermal and pressure shock loading determined by Eq. (45) is applied on the left side of the plate. The distribution of the temperature θ, the displacement u1, and the stress σx along the x-axis for the uncoupled (C 1⁄4 0) and coupled (C 1⁄4 0.1) thermo- elastic cases at dimensionless time t 1⁄4 3 and t 1⁄4 6 are shown in Figs. 15–17, respectively. It can be found that the results could be also superimposed from those of cases (1) and (2) as discussed above. The results of the present method have good agreement with the finite element solutions for the uncoupled and coupled cases too.
5. Conclusions
A meshless local Petrov–Galerkin method in conjunction with a backward difference method is successfully implemented and applied to two-dimensional linear dynamic coupled thermoelasticity pro- blems under thermal and mechanical shock loading. In the present
B.-J. Zheng et al. / Engineering Analysis with Boundary Elements 60 (2015) 154–161 159
Fig. 13. Comparison of the displacement u1 along x-axis for pressure shock loading.
160 B.-J. Zheng et al. / Engineering Analysis with Boundary Elements 60 (2015) 154–161 Acknowledgments
Fig. 15. Comparison of the temperature θ along x-axis for combined thermal and pressure shock.
The authors acknowledge the support by the National Natural Science Foundation of China (11172055, 11202045) and the Ger- man Research Foundation (DFG, Project-no. ZH 15/10-2). Baojing Zheng also acknowledges the support by the China Scholarship Council (CSC) for his PhD visit study at the Chair of Structural Mechanics, Faculty of Science and Technology, University of Sie- gen, Germany.
References
[1] Danilovskaya VI. Thermal stresses in an elastic half-space duo to a sudden heating of its boundary (in Russian). Prikl Mat Mekh 1950;14:316–8.
[2] Danilovskaya VI. On a dynamic problem of thermoelasticity (in Russian). Prikl Mat Mekh 1952;16:341–4.
[3] Boley BA, Tolins IS. Transient coupled thermoelastic boundary value problems in the half space. J Appl Mech Trans ASME 1963;29:637–46.
[4] Keramidas GA, Ting EC. A finite element formulation for thermal stress analysis. I. Variational formulation, II. Finite element formulation. Nucl Eng Des 1976;39:267–87.
[5] Ting EC, Chen HC. A unified numerical approach for thermal stress waves. Comput Struct 1982;15:165–75.
[6] Prevost JH, Tao D. Finite element analysis of dynamic coupled thermoelasticity problem with relaxation times. J Appl Mech Trans ASME 1983;50:817–22.
[7] Tamam KK, Railkar SB. On heat displacement based hybrid transfinite element
formulations for uncoupled/coupled thermally induced stress wave propaga-
tion. Comput Struct 1988;30:1025–36.
[8] Cannarozzi AA, Ubertini F. A mixed variational method for linear coupled
thermoelastic analysis. Int J Solids Struct 2001;38:717–39.
[9] Rizzo FJ, Shippy DJ. An advanced boundary integral equation method for three-
dimensional thermoelasticity. Int J Numer Methods Eng 1977;11:1753–68.
[10] Gao XW. Boundary element analysis in thermoelasticity with and without
internal cells. Int J Numer Methods Eng 2003;57:957–90.
[11] Dargush GF, Banerjee PK. A new boundary element method for three- dimensional coupled problems of conduction and thermoelasticity. J Appl
Mech Trans ASME 1991;58:28–36.
[12] Suh IG, Tosaka N. Application of the boundary element method to 3-D linear
coupled thermoelasticity problems. Theor Appl Mech 1989;38:169–75.
[13] Tosaka N, Suh IG. Boundary element analysis of dynamic coupled thermo-
elasticity problems. Comput Mech 1991;8:331–42.
[14] Chen J, Dargush GF. Boundary element method for dynamic poroelastic and
thermoelastic analyses. Int J Solids Struct 1995;32:2257–78.
[15] Hosseini-Tehrani P, Eslami MR. BEM analysis of thermal and mechanical shock in a two-dimensional finite domain considering coupled thermoelasticity. Eng
Anal Bound Elem 2000;24:249–57.
[16] Sladek V, Sladek J. Boundary integral equation method in thermoelasticity Part I:
general analysis. Appl Math Model 1983;7:241–53.
[17] Sladek J, Sladek V. Boundary integral equation method in thermoelasticity Part II:
crack analysis. Appl Math Model 1984;8:27–36.
[18] Belytschko T, Lu YY, Gu L. Element-free Galerkin methods. Int J Numer
Methods Eng 1994;37:229–56.
[19] Liu WK, Jun S, Zhang Y. Reproducing kernel particle methods. Int J Numer
Methods Fluids 1995;20:1081–106.
[20] Atluri SN, Zhu TL. A new meshless local Petrov–Galerkin (MLPG) approach in
computational mechanics. Comput Mech 1998;22:117–27.
[21] Liu GR, Gu YT. A point interpolation method for two dimensional solid.
Int J Numer Methods Eng 2001;50:937–51.
[22] Chen L, Cheng YM. The complex variable reproducing kernel particle method
for two-dimensional elastodynamics. Chin Phys B 2010;19:090204.
[23] Zhang X, Song KZ, Lu MW. Meshless methods based on collocation with radial
basis function. Comput Mech 2000;26:333–43.
[24] Sladek V, Sladek J, Atluri SN. A pure contour formulation for meshless local
boundary integral equation method in thermoelasticity. Comput Model Eng
Sci 2001;2:423–34.
[25] Bobaru F, Mukherjee S. Meshless approach to shape optimization of linear
thermoelastic solids. Int J Numer Methods Eng 2003;53:765–96.
[26] Qian LF, Batra RC. Transient thermoelastic deformations of thick functionally
graded plate. J Therm Stresses 2004;27:705–40.
[27] Sladek J, Sladek V, Ch Zhang, Tan CL. Meshless local Petrov–Galerkin method
for linear coupled thermoelastic analysis. Comput Model Eng Sci
2006;16:57–68.
[28] Sladek J, Sladek V, Solek P, Tan CL, Zhang Ch. Two- and three-dimensional
transient thermoelastic analysis by the MLPG method. Comput Model Eng Sci
2009;47:61–95.
[29] Hosseini SM, Sladek J, Sladek V. Meshless Petrov–Galerkin method for coupled
thermoelasticity analysis of a functionally graded thick hollow cylinder. Eng
Anal Bound Elem 2011;35:827–35.
[30] Hosseini SM, Sladek J, Sladek V. Application of meshless local integral
equations to two dimensional analysis of coupled non-Fick diffussion-elasti- city. Eng Anal Bound Elem 2013;37:603–15.
Fig. 16. Comparison of the displacement u1 along x-axis for combined thermal and pressure shock.
Fig. 17. Comparison of the stress σx along x-axis for combined thermal and pressure shock.
method, the moving Kriging interpolation is employed to construct the shape functions, and the Heaviside step function is applied in the local weak-form as the test functions. In the dynamic coupled thermoelasticity problems the equations of motion and the transient heat conduction equation interact with each other, hence these equations must be solved simultaneously. In this paper the backward difference method is used to discretize the time domain, while the velocities as extra field variables are used to obtain the reduction procedure for the equations of motion. A unified explicit expression of the system of differential algebraic equations is obtained. A quasi- static thermoelasticity problem and a dynamic coupled thermoelas- ticity problem are studied to verify the accuracy and the efficiency of the present MLPG method.
[31] Zhu P, Zhang LW, Liew KM. Geometrically nonlinear thermomechanical analysis of moderately thick functionally graded plates using a local Petrov– Galerkin approach with moving Kriging interpolation. Compos Struct 2014;107:298–314.
[32] Zheng BJ, Dai BD. A meshless local moving Kriging method for two- dimensional solids. Appl Math Comput 2011;218:563–73.
[33] Dai BD, Cheng J, Zheng BJ. A moving Kriging interpolation-based meshless local Petrov–Galerkin method for elastodynamic analysis. Int J Appl Mech 2013;5:1350011.
[34] Gu L. Moving Kriging interpolation and element-free Galerkin method. Int J Numer Methods Eng 2003;56:1–11.
[35] Li XG, Dai BD, Wang LH. A moving Kriging interpolation-based boundary node method for two-dimensional potential problems. Chin Phys B 2010;19:120202–7.
[36] Zheng BJ, Dai BD. Improved meshless local Petrov–Galerkin method for two- dimensional potential problems. Acta Phys Sin 2010;59:5182–9 (in Chinese).
[37] Dai BD, Zheng BJ. Numerical solution of transient heat conduction problems using improved meshless local Petrov–Galerkin method. Appl Math Comput 2013;219:10044–52.
[38] Carslaw HS, Jaeger JC. Conduction of heat in solids. Oxford: Clarendon; 1959.
B.-J. Zheng et al. / Engineering Analysis with Boundary Elements 60 (2015) 154–161 161