The MFS–MPS for two-dimensional steady-state thermoelasticity problems

Engineering Analysis with Boundary Elements 37 (2013) 1004–1020

Contents lists available at SciVerse ScienceDirect

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

The MFS–MPS for two-dimensional steady-state thermoelasticity problems Liviu Marin a,b,n, Andreas Karageorghis c a

Institute of Solid Mechanics, Romanian Academy, 15 Constantin Mille, P.O. Box 1–863, 010141 Bucharest, Romania Centre for Continuum Mechanics, Faculty of Mathematics and Computer Science, University of Bucharest, 14 Academiei, 010014 Bucharest, Romania c Department of Mathematics and Statistics, University of Cyprus, P.O. Box 20537, 1678 Nicosia, Cyprus b

art ic l e i nf o

a b s t r a c t

Article history: Received 31 July 2012 Accepted 6 April 2013

We consider the numerical approximation of the boundary and internal thermoelastic fields in the case of two-dimensional isotropic linear thermoelastic solids by combining the method of fundamental solutions (MFS) with the method of particular solutions (MPS). A particular solution of the nonhomogeneous equations of equilibrium associated with a planar isotropic linear thermoelastic material is derived from the MFS approximation of the boundary value problem for the heat conduction equation. Moreover, such a particular solution enables one to easily develop analytical solutions corresponding to any two-dimensional domain occupied by an isotropic linear thermoelastic solid. The accuracy and convergence of the proposed MFS–MPS procedure are validated by considering three numerical examples. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Linear thermoelasticity Navier–Lamé system Method of fundamental solutions (MFS) Method of particular solutions (MPS)

1. Introduction Numerous stress analysis problems in engineering deal with structures which are simultaneously subject to thermal and mechanical loadings, i.e. thermoelastic loadings. This type of problems is encountered whenever a solid is subject to heating conditions that give rise to a temperature distribution throughout its volume. This temperature distribution produces thermal expansions in the object under consideration. In an isotropic material, at a uniform reference temperature, a small uniform increase in the temperature field can produce a pure volumetric expansion, provided that the solid body is not constrained against such a movement. This phenomenon can be expressed in terms of the so-called thermal strain, which is related to the difference between the temperature of the solid and the reference temperature via the coefficient of thermal expansion. It is important to mention that such a thermal expansion may also occur with no stresses present in the solid body, see e.g. [1]. Mathematical problems of isotropic thermoelasticity have already been addressed by many authors in the literature using various numerical methods such as the boundary element method (BEM) [2,9,11,24,25,28,29], the dual reciprocity BEM (DRBEM) [23],

n Corresponding author at: Institute of Solid Mechanics, Romanian Academy, 15 Constantin Mille, P.O. Box 1–863, 010141 Bucharest, Romania. Tel./fax: +40 21 312 6736. E-mail addresses: [email protected], [email protected] (L. Marin), [email protected] (A. Karageorghis).

0955-7997/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.enganabound.2013.04.002

the finite element method (FEM) [3,4], the moving least-squares method combined with the local boundary integral method [30], etc. The method of fundamental solutions (MFS) is a meshless boundary collocation method which is applicable to boundary value problems for which a fundamental solution of the operator in the governing equation is known. In spite of this restriction, the MFS has become very popular primarily because of the ease with which it can be implemented, in particular for problems in complex geometries. It is important to mention that the MFS can be employed in conjunction with the so-called singularity subtraction technique (SST) for problems with boundary and/or solution singularities, in the sense that the standard MFS solution of the partial differential equation investigated is augmented by suitable singular functions/eigenfunctions associated with the corresponding partial differential operator as given by the asymptotic expansion of the solution near the singular point. For details on the MFS-SST technique applied to the numerical solution of direct and inverse problems associated with the Laplace, biharmonic, Helmholtz and modified Helmholtz, and Cauchy–Navier equations, we refer the reader to [12,17–19,15]. Moreover, the MFS does not involve integrations which could be potentially troublesome and complicated as is the case with the BEM. Another advantage of the MFS over domain discretisation methods such as the FEM or the finite difference method is the fact that the MFS is a boundary collocation method and hence only the boundary of the solution domain needs to be discretised. The major disadvantages of the MFS are related to the optimal location of the singularities, which is an important issue especially from the numerical viewpoint, and its inability to be directly applied to

L. Marin, A. Karageorghis / Engineering Analysis with Boundary Elements 37 (2013) 1004–10205

non-homogeneous or nonlinear problems, for which a fundamental solution is not available. Since its introduction as a numerical method by Mathon and Johnston [22], the MFS has been successfully applied to a large variety of physical problems, an account of which may be found in the survey papers [5,6,13]. The MFS, in conjunction with the method of particular solutions (MPS) and the dual reciprocity method, was applied to the numerical solution of direct problems in three-dimensional isotropic linear thermoelasticity in [16,33], respectively. Here, we aim to provide, apparently for the first time, a comprehensive study of the application of the MFS–MPS to direct two-dimensional isotropic linear thermoelasticity problems. Some preliminary results of this work may be found in [20]. Consequently, it is the purpose of this study to investigate the numerical approximation of the boundary and internal thermoelastic fields in the case of two-dimensional isotropic linear thermoelastic solids by combining the MFS with the MPS. More precisely, from the standard MFS approximation of the boundary value problem for the heat conduction equation, we first develop a particular solution of the non-homogeneous equations of equilibrium associated with a planar isotropic linear thermoelastic material. It is noteworthy that such a particular solution enables one to easily develop analytical solutions corresponding to any two-dimensional domain occupied by an isotropic linear thermoelastic solid. The accuracy and convergence of the proposed MFS– MPS algorithm are validated by analysing three numerical examples in both simply and doubly connected domains.

2. Mathematical formulation Consider a domain Ω⊂R2 which is bounded by a smooth or piecewise smooth curve ∂Ω and occupied by an isotropic solid characterised by the thermal conductivity, κ, the coefficient of linear thermal expansion, αT , Poisson's ratio, ν, and the shear modulus, G. In the framework of isotropic linear thermoelasticity, the strain tensor, ϵ ¼ ½ϵij 1 ≤ i;j ≤ 2 , is related to the stress tensor, r ¼ ½sij 1 ≤ i;j ≤ 2 , by means of the constitutive law of thermoelasticity, namely
 1 ν rðxÞ− trðrðxÞÞI þ α T TðxÞ I; x∈Ω ¼ Ω∪∂Ω: ð1Þ ϵðxÞ ¼ 2G 1þν Here I ¼ ½δij 1 ≤ i;j ≤ 2 denotes the identity matrix in R22 , ν ¼ ν for a plane strain state and ν ¼ ν=ð1 þ νÞ for a plane stress state, while α T ¼ αT and α T ¼ αT ð1 þ νÞ=ð1 þ 2νÞ for the plane strain and plane stress states, respectively. We note from Eq. (1) that the shear strains are not affected by the temperature as the free thermal expansion does not produce any angular distortion in an isotropic material. The constitutive law of thermoelasticity (1) can also be conveniently expressed in terms of the stresses as
 ν trðϵðxÞÞI −γ TðxÞI; x∈Ω; rðxÞ ¼ 2G ϵðxÞ þ ð2Þ 1−2ν where γ ¼ 2Gα T ð1 þ νÞ=ð1−2νÞ. The equations of equilibrium are similar to those of isotropic linear elasticity since they are based on purely mechanical considerations. The stress tensor can be expressed in terms of the displacement derivatives by combining the constitutive law (2)

0.8

0.8

0.6

0.6

d2

1005

d2 0.4

0.4

0.2

0.2

0.2

0.4

0.6

0.8

1.0

0.2

0.4

d1

0.6

0.8

1.0

0.6

0.8

1.0

d1

0.8

0.8

0.6

0.6

d2

d2 0.4

0.4

0.2

0.2

0.2

0.4

0.6 d1

0.8

1.0

0.2

0.4 d1

Fig. 1. The maximum normalised error Errðu1 Þ on ∂Ω as a function of the distances d1 and d2 , obtained using various numbers of collocation points and sources, namely (a) M ¼ N ¼ 48, (b) M ¼ N ¼ 60, (c) M ¼ N ¼ 72 and (d) M ¼ N ¼ 84, for Example 1.

1006

L. Marin, A. Karageorghis / Engineering Analysis with Boundary Elements 37 (2013) 1004–1020

with the kinematic relation ϵðxÞ ¼

1 2ð∇uðxÞ

T

þ ∇uðxÞ Þ;

ð3Þ

x∈Ω;

to yield


 2ν ð∇  uðxÞÞI −γ TðxÞI; rðxÞ ¼ G ð∇uðxÞ þ ∇uðxÞT Þ þ 1−2ν

x∈Ω:

ð4Þ

By assuming the absence of body forces, one obtains the equilibrium equations of two-dimensional isotropic linear thermoelasticity in terms of the displacement vector and the temperature (also known as the Navier–Lamé system of two-dimensional isotropic linear thermoelasticity), namely −∇  rðxÞ≡LuðxÞ þ γ ∇TðxÞ ¼ 0;

ð5Þ

x∈Ω;

where L ¼ ðL1 ; L2 ÞT is the partial differential operator associated with the Navier–Lamé system of isotropic linear elasticity, i.e.
 2ν ∇ð∇  uðxÞÞ ; x∈Ω: ð6Þ LuðxÞ≡−G ∇  ð∇uðxÞ þ ∇uðxÞT Þ þ 1−2ν Further, we let nðxÞ ¼ ðn1 ðxÞ; n2 ðxÞÞT be the outward unit normal vector at ∂Ω and tðxÞ be the traction vector at x∈∂Ω given by tðxÞ≡rðxÞnðxÞ;

x∈Γ u ;

ð8aÞ

and ~ tðxÞ ¼ tðxÞ;

qðxÞ≡−ðκ ∇TðxÞÞ  nðxÞ;

ð9Þ

x∈∂Ω:

If the temperature and normal heat flux are prescribed on the boundary segments Γ T ⊂∂Ω and Γ q ⊂∂Ω, respectively, where Γ T ∩Γ q ¼ ∅ and Γ T ∪Γ q ¼ ∂Ω, then one obtains the mixed boundary value problem associated with the two-dimensional steady-state heat conduction. In the absence of heat sources, this problem can be recast as −∇  ðκ ∇TðxÞÞ ¼ 0;

ð10aÞ

x∈Ω;

TðxÞ ¼ T~ ðxÞ;

x∈Γ T ;

ð10bÞ

~ qðxÞ ¼ qðxÞ;

x∈Γ q :

ð10cÞ

3. Solution algorithm

ð7Þ

x∈∂Ω:

If the boundary segments Γ u ⊂∂Ω and Γ t ⊂∂Ω are such that Γ u ∩Γ t ¼ ∅ and Γ u ∪Γ t ¼ ∂Ω, and the displacement and traction vectors are prescribed on these boundaries, i.e. ~ uðxÞ ¼ uðxÞ;

then equations (5), (8a) and (8b) represent the mixed boundary value problem associated with the equations of equilibrium for two-dimensional steady-state isotropic linear thermoelasticity. Next, we let qðxÞ be the normal heat flux at a point x∈∂Ω defined by

x∈Γ t ;

ð8bÞ

We propose the numerical solution of boundary value problems associated with two-dimensional isotropic linear thermoelastic solids, i.e. Eqs. (5), (8a), (8b), and (10a)–(10c), using the MFS in conjunction with the MPS. Assuming that κ is constant, the heat conduction boundary value problem (10a)–(10c) is first solved numerically by applying the standard MFS. Next, based on this numerical approximation for the temperature field in the solid, we derive (apparently for the first time) a particular solution of the

0.8

0.8

0.6

0.6

d2

d2 0.4

0.4

0.2

0.2

0.2

0.4

0.6

0.8

1.0

0.2

0.4

d1

0.6

0.8

1.0

0.6

0.8

1.0

d1

0.8

0.8

0.6

0.6

d2

d2 0.4

0.4

0.2

0.2

0.2

0.4

0.6 d1

0.8

1.0

0.2

0.4 d1

Fig. 2. The maximum normalised error ErrðqÞ on ∂Ω as a function of the distances d1 and d2 , obtained using various numbers of collocation points and sources, namely (a) M ¼ N ¼ 48, (b) M ¼ N ¼ 60, (c) M ¼ N ¼ 72 and (d) M ¼ N ¼ 84, for Example 1.

L. Marin, A. Karageorghis / Engineering Analysis with Boundary Elements 37 (2013) 1004–10205

and traction vector

equilibrium equations of two-dimensional isotropic linear thermoelasticity (5) or, equivalently, the non-homogeneous equilibrium equations associated with two-dimensional isotropic linear elasticity. Finally, we apply the MFS to the resulting boundary value problem corresponding to the homogeneous equilibrium equations for a two-dimensional isotropic linear elastic material. The numerical procedure described above may be summarised as follows:

tðPÞ ðxÞ ¼ rðPÞ ðxÞ nðxÞ;

LuðHÞ ðxÞ ¼ 0;

stress tensor
rðPÞ ðxÞ ¼ 2G ϵðPÞ ðxÞ þ

 ν trðϵðPÞ ðxÞÞI ; 1−2ν

0.3

0.2

0.2

0.1

0.1

u1 0.0

u2 0.0 Analytical M = N = 48 M = N = 60 M = N = 72 M = N = 84

-0.2 -0.3 -1.0

-0.8

-0.6

-0.4

x∈Γ t ;

ð12cÞ

Analytical M = N = 48 M = N = 60 M = N = 72 M = N = 84

-0.1 -0.2

-0.2

ð12bÞ

using the MFS to determine the unknown boundary displacement uðHÞ jΓt and traction tðHÞ jΓu , as well as uðHÞ jΩ , ϵðHÞ jΩ and rðHÞ jΩ . Step 4. On applying the superposition principle, determine the unknown boundary displacement ujΓt ¼ uðHÞ jΓt þ uðPÞ jΓ t and boundary traction tjΓ u ¼ tðHÞ jΓu þ ðtðPÞ −γ T nÞjΓu , as well as the mechanical fields inside the domain, namely ujΩ ¼ uðHÞ jΩ þ uðPÞ jΩ , ϵjΩ ¼ ϵðHÞ jΩ þ ϵðPÞ jΩ and rjΩ ¼ rðHÞ jΩ þ ðrðPÞ −γ T IÞjΩ .

0.3

-0.1

x∈Γ u ;

ðPÞ ~ ðxÞ−γ TðxÞnðxÞ; tðHÞ ðxÞ ¼ tðxÞ−½t

ð11bÞ

x∈Ω;

ð12aÞ

x∈Ω;

ðPÞ ~ ðxÞ; uðHÞ ðxÞ ¼ uðxÞ−u

ð11aÞ

x∈Ω;

ð11cÞ

x∈∂Ω:

Step 3. Solve the resulting direct problem for the homogeneous equilibrium equations in two-dimensional isotropic linear elasticity, i.e.

Step 1. Solve the thermal problem (10a)–(10c) using the MFS to determine the unknown boundary temperature TjΓq and flux qjΓT , as well as the temperature distribution in the domain TjΩ . Step 2. Determine a particular solution uðPÞ of the non-homogeneous equilibrium equations for two-dimensional isotropic linear elasticity (5) in Ω, as well as the corresponding particular strain tensor: ϵðPÞ ðxÞ ¼ 12ð∇uðPÞ ðxÞ þ ∇uðPÞ ðxÞT Þ;

1007

0.0

-0.3 -1.0

-0.8

/(2)

-0.6

-0.4

-0.2

0.0

/(2)

-5.76

-5.77

-5.78 q -5.79

-5.8

-5.81 -1.0

-0.8

-0.6 -0.4 /(2) Analytical M = N = 60 M = N = 84

-0.2

0.0

M = N = 48 M = N = 72

Fig. 3. The analytical and numerical displacements (a) u1 jΓ int and (b) u2 jΓ int , and (c) normal heat flux qjΓ int , obtained using d1 ¼ 1:0, d2 ¼ 0:30 and M ¼ N∈f48; 60; 72; 84g, for Example 1.

1008

L. Marin, A. Karageorghis / Engineering Analysis with Boundary Elements 37 (2013) 1004–1020

It should be stressed that the solution uðHÞ of the boundary value problem for the homogeneous equilibrium equations in twodimensional isotropic linear elasticity (12a)–(12c) is the homogeneous term of the solution u of the non-homogeneous boundary value problem (5), (8a) and (8b) associated with the mechanical field for two-dimensional isotropic linear thermoelasticity.

Proof. The following relations hold for any x∈R2 \X K , 1 ≤ k ≤ K, K∈Zþ and i; j ¼ 1; 2: ∂ ðx −xðkÞ Þ ¼ δij ; ∂xj i i

ð14aÞ

3.1. Particular solutions

xj −xðkÞ ∂ j ð∥x−xðkÞ ∥Þ ¼ ; ∂xj ∥x−xðkÞ ∥

ð14bÞ

The proposed MFS–MPS algorithm relies on the existence of a particular solution of the non-homogeneous equilibrium equations (5) and this is justified by the following result:

xj −xðkÞ ∂ j ðlog∥x−xðkÞ ∥Þ ¼ ; ∂xj ∥x−xðkÞ ∥2

ð14cÞ

Proposition 1. Let Ω⊂R2 be a domain occupied by an isotropic solid characterised by the constant thermal conductivity, κ, the coefficient of linear thermal expansion, αT , Poisson's ratio, ν, and the shear modulus, G, respectively. K Then for any set X K ¼ fxðkÞ gk ¼ 1 ⊂R2 \Ω, K∈Zþ , the temperature field K

T ðPÞ ðxÞ ¼ ∑ T k log∥x−xðkÞ ∥;

ð13aÞ

x∈Ω;

where δij denotes the Kronecker delta. Clearly, the function T ðPÞ given by Eq. (13a) satisfies the Laplace equation (10a) in R2 \X K and hence in Ω. Moreover, by differentiating expression (13a) with respect to xi and using identity (14c), one obtains

k¼1

where T k ∈R, 1 ≤k ≤K, and displacement vector   K γ 1−2ν ∑ T ðx−xðkÞ Þ log∥x−xðkÞ ∥; uðPÞ ðxÞ ¼ 4G 1−ν k ¼ 1 k

γ

K xi −xðkÞ ∂ ðPÞ i T ðxÞ ¼ γ ∑ T k ; ðkÞ ∥2 ∂xi ∥x−x k¼1

where γ ¼ 2Gα T ð1 þ νÞ=ð1−2νÞ. If for 1 ≤k ≤K, K∈Zþ , we use the notation ðkÞ ðkÞ uðkÞ i ðxÞ ¼ ðx−x Þ log∥x−x ∥;

0.4

0.4

0.2

0.2

u1 0.0

-0.4

u2 Analytical M = N = 48 M = N = 60 M = N = 72 M = N = 84 0.0

0.2

0.4

0.6 /(2)

x∈R2 \X K ; i ¼ 1; 2;

Analytical M = N = 48 M = N = 60 M = N = 72 M = N = 84

-0.4 0.0

1.0

Analytical M = N = 60 M = N = 84

ð16aÞ

0.0

-0.2

0.8

-2.8918

ð15Þ

ð13bÞ

x∈Ω;

where γ ¼ 2G α T ð1 þ νÞ=ð1−2νÞ, represent a particular solution ðT ðPÞ ; uðPÞ Þ∈ðC ∞ ðΩÞ3 Þ of the governing thermo-mechanical equations of two-dimensional isotropic linear thermoelasticity (5) and (10a).

-0.2

x∈R2 \X K ; i ¼ 1; 2;

0.2

0.4

0.6 /(2)

0.8

1.0

M = N = 48 M = N = 72

-2.892

q -2.8922

-2.8924

-2.8926 0.0

0.2

0.4

0.6 /(2)

0.8

1.0

Fig. 4. The analytical and numerical displacements (a) u1 jΓ out and (b) u2 jΓ out , and (c) normal heat flux qjΓ out , obtained using d1 ¼ 1:0, d2 ¼ 0:30 and M ¼ N∈f48; 60; 72; 84g, for Example 1.

L. Marin, A. Karageorghis / Engineering Analysis with Boundary Elements 37 (2013) 1004–10205

then by differentiating expression (16a) with respect to xj and employing identities (14a) and (14c), we obtain xj −xðkÞ xi −xðkÞ ∂ ðkÞ j i ui ðxÞ ¼ log∥x−xðkÞ ∥δij þ ; ∂xj ∥x−xðkÞ ∥ ∥x−xðkÞ ∥

x∈R2 \X K ; i; j ¼ 1; 2: ð16bÞ

Further, taking the derivative of expression (16b) with respect to xj yields " # 2 xj −xðkÞ ðxj −xðkÞ xi −xðkÞ ∂2 ðkÞ j j Þ i u ðxÞ ¼ 2 δ þ 1−2 ; ij ∥x−xðkÞ ∥2 ∥x−xðkÞ ∥2 ∥x−xðkÞ ∥2 ∂x2j i x∈R2 \X K ; i; j ¼ 1; 2:

xi −xðkÞ i ; ∥x−xðkÞ ∥2

x∈R2 \X K ; i ¼ 1; 2:

ð18bÞ

Substituting relations (18a) and (18b) into the expression of the i-th component, Li ; i ¼ 1; 2, of the partial differential operator associated with the Navier–Lamé system of two-dimensional isotropic linear elasticity, L, see Eq. (6), yields " # 2 ∂2 ν ∂ 2 ∂ ðkÞ ðkÞ ðxÞ ¼ −2G ∑ u ðxÞ þ ∑ u ðxÞ Li uðkÞ i 2 i 1−2ν ∂xi j ¼ 1 ∂xj j j ¼ 1 ∂xj  ¼ −4G

ð17Þ

1−ν 1−2ν



xi −xðkÞ 2γ i ¼− αT ∥x−xðkÞ ∥2



1−ν 1þν



xi −xðkÞ i ; ∥x−xðkÞ ∥2

x∈R2 \X K ; i ¼ 1; 2;

Next, from Eqs. (16b) and (17), one obtains xi −xðkÞ ∂ 2 ∂ ðkÞ ∂ i ∑ u ðxÞ ¼ ð2 log∥x−xðkÞ ∥ þ 1Þ ¼ 2 ; ∂xi j ¼ 1 ∂xj j ∂xi ∥x−xðkÞ ∥2 x∈R2 \X K ; i ¼ 1; 2;

ð18aÞ

ð19Þ

and, consequently, one obtains "  #  K K xi −xðkÞ αT 1 þ ν ðPÞ ðkÞ i ∑ T k ui ðxÞ ¼ −γ ∑ T k ; Li ui ðxÞ ¼ Li 1−ν 2 ∥x−xðkÞ ∥2 k¼1 k¼1 x∈R2 \X K ; i ¼ 1; 2:

and " # 2 2 ðxj −xðkÞ xi −xðkÞ ∂2 ðkÞ j Þ i u ðxÞ ¼ ∑ 2δ þ 1−2 ij 2 i ∥x−xðkÞ ∥2 j ¼ 1 ∥x−xðkÞ ∥2 j ¼ 1 ∂xj

ð20Þ

Finally, from Eqs. (15) and (20), it follows that the displacement vector uðPÞ given by (13b) is a particular solution of the Navier– Lamé system of two-dimensional isotropic linear thermoelasticity

2

∑

1.0E-02

1.0E-02

1.0E-04

1.0E-04 err(u2(x))

err(u1(x))

¼2

1009

1.0E-06 1.0E-08

1.0E-06 1.0E-08

1.0E-10 1.0E-10 1.0E-12 -1.0

-0.8

-0.6 -0.4 /(2)

M = N = 48 M = N = 72

-0.2

0.0

1.0E-12 -1.0

-0.8

M = N = 60 M = N = 84

-0.6 -0.4 /(2)

M = N = 48 M = N = 72

-0.2

0.0

M = N = 60 M = N = 84

1.0E-02

err(q(x))

1.0E-03 1.0E-04 1.0E-05 1.0E-06 1.0E-07 -1.0

-0.8

-0.6

-0.4

-0.2

0.0

/(2) M = N = 48 M = N = 72

M = N = 60 M = N = 84

Fig. 5. The pointwise normalised errors (a) errðu1 ðxÞÞ, (b) errðu2 ðxÞÞ, and (c) errðqðxÞÞ, x∈Γ int , obtained using d1 ¼ 1:0, d2 ¼ 0:30 and M ¼ N∈f48; 60; 72; 84g, for Example 1.

1010

L. Marin, A. Karageorghis / Engineering Analysis with Boundary Elements 37 (2013) 1004–1020

1.0E-01

1.0E-01

1.0E-02

1.0E-02

1.0E-03

1.0E-03

1.0E-04

1.0E-04

1.0E-05

1.0E-05 M = N = 48 M = N = 72

M = N = 48 M = N = 72

M = N = 60 M = N = 84

1.0E-06

M = N = 60 M = N = 84

1.0E-06 0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

M = N = 48 M = N = 72

M = N = 60 M = N = 84

0.4

0.6

0.8

1.0

1.0E-03 5.0E-04

1.0E-04 5.0E-05

1.0E-05 5.0E-06

1.0E-06

0.0

0.2

0.4

0.6

0.8

1.0

Fig. 6. The pointwise normalised errors (a) errðu1 ðxÞÞ, (b) errðu2 ðxÞÞ, and (c) errðqðxÞÞ, x∈Γ out , obtained using d1 ¼ 1:0, d2 ¼ 0:30 and M ¼ N∈f48; 60; 72; 84g, for Example 1.

0.03 0.025 0.02

0.0

t2/1010

t1/1010

0.02

-0.02

Analytical M = N = 32 M = N = 40 M = N = 48 M = N = 56

-0.04 0.0

0.1

0.2

0.015 Analytical M = N = 32 M = N = 40 M = N = 48 M = N = 56

0.01 0.005 0.3

0.4

0.0

0.5

0.0

0.1

0.2

0.3

0.4

0.5

60 Analytical M = N = 32 M = N = 40 M = N = 48 M = N = 56

50

T

40 30 20 10 0.0

0.1

0.2

0.3

0.4

0.5

Fig. 7. The analytical and numerical tractions (a) t 1 jΓ u , (b) t 2 jΓ u , and (c) Tjrq temperature qjΓT , where Γ u ¼ Γ q ¼ fx∈R2 j ∥x∥ ¼ R; θ∈½0; πÞg, obtained using d ¼ 0:50 and M ¼ N∈f32; 40; 48; 56g, for Example 2.

L. Marin, A. Karageorghis / Engineering Analysis with Boundary Elements 37 (2013) 1004–10205

(5) in R2 \X K and hence in Ω, namely ∂ ðPÞ T ðxÞ ¼ 0; ∂xi

Li uðPÞ i ðxÞ þ γ

x∈Ω; i ¼ 1; 2:

Moreover, Proposition 1 enables one to easily construct analytical solutions in two-dimensional isotropic linear thermoelasticity regardless of the geometry of the domain. This is achieved by simply setting K∈Zþ and choosing a set of points X K ¼ K fxðkÞ gk ¼ 1 ⊂R2 \Ω, as well as constants T k ∈R, 1 ≤k ≤K. The corresponding analytical solution of the governing thermo-mechanical equations of two-dimensional isotropic linear thermoelasticity (5) and (10a) is given by Eqs. (13b) and (13a), respectively.

ð21Þ

The regularity of the particular solution, i.e. ðT ðPÞ ; uðPÞ Þ∈ðC ∞ ðΩÞÞ3 , follows directly from the expressions of its components as given by Eqs. (13a) and (13b). □ As a direct consequence of Proposition 1, the expressions for the corresponding particular strain tensor ϵðPÞ , stress tensor rðPÞ , and traction vector tðPÞ are obtained by substituting the particular displacement vector given by (13b) into Eqs. (11a), (11b), and (11c), respectively, i.e. ϵðPÞ ðxÞ ¼

γ 4G



1−2ν 1−ν




 K x−xðkÞ x−xðkÞ ⊗ ; ∑ T k log ∥x−xðkÞ ∥ I þ ðkÞ ∥ ðkÞ ∥ ∥x−x ∥x−x k¼1

4. Method of fundamental solutions Step 1. For the first step of the algorithm described in Section 3, consider the fundamental solution, F, of the two-dimensional steady-state heat conduction equation (10a) in an isotropic homogeneous medium [5], namely 1 Fðx; ξÞ ¼ − log∥x−ξ∥; x∈Ω; ð23Þ 2πκ

x∈Ω;

ð22aÞ r

ðPÞ





where x ¼ ðx1 ; x2 Þ is a collocation point and ξ ¼ ðξ1 ; ξ2 Þ∈R2 \Ω is a singularity or source point. The main idea in the MFS is to approximate the temperature in the solution domain by a linear




 1  ∑ Tk log ∥x−xðkÞ ∥ þ ν I 1−2ν k¼1  x−xðkÞ x−xðkÞ þ ⊗ ; x∈Ω; ðkÞ ðkÞ ∥x−x ∥ ∥x−x ∥

γ ðxÞ ¼ 2

1−2ν 1−ν

K

combination of fundamental solutions with respect to NLs singula-

ð22bÞ

NL

rities, fξðnÞ gn s¼ 1 , in the form NLs

and tðPÞ ðxÞ ¼

1011

γ 2



1−2ν 1−ν



K

TðxÞ≈T NL ðcð1Þ ; ξ; xÞ ¼ ∑ cnð1Þ Fðx; ξðnÞ Þ;




s

1 ðlog∥x−xðkÞ ∥ þ νÞnðxÞ 1−2ν  ðx−xðkÞ Þ  nðxÞ þ ðx−xðkÞ Þ ; x∈∂Ω: ðkÞ 2 ∥x−x ∥ ∑ Tk

L

NL

coordinates of the singularities fξðnÞ gn s¼ 1 . From Eqs. (9) and (23) it follows that the normal heat flux, through a curve defined by the

ð22cÞ

*10-4 Analytical M = N = 32 M = N = 40 M = N = 48 M = N = 56

Analytical M = N = 32 M = N = 40 M = N = 48 M = N = 56

-5

u2

u1

L

s

*10-3 0.0

-1.0

ð24Þ

ð1Þ T Ns and ξ∈R2Ns is a vector containing the where cð1Þ ¼ ½cð1Þ 1 ; …; cNL  ∈R

k¼1

-0.5

x∈Ω;

n¼1

-1.5

-10

-2.0 -15 -2.5 -3.0 0.5

0.6

0.7

0.8

0.9

1.0

0.5

0.6

0.7

0.8

0.9

1.0

Analytical M = N = 32 M = N = 40 M = N = 48 M = N = 56

q

50

0

-50 0.5

0.6

0.7

0.8

0.9

1.0

Fig. 8. The analytical and numerical displacements (a) u1 jΓ t , (b) u2 jΓ t , and (c) qjrT normal heat flux qjΓq , where Γ t ¼ Γ T ¼ fx∈R2 j∥x∥ ¼ R; θ∈½π; 2πÞg, obtained using d ¼ 0:50 and M ¼ N∈f32; 40; 48; 56g, for Example 2.

1012

L. Marin, A. Karageorghis / Engineering Analysis with Boundary Elements 37 (2013) 1004–1020

outward unit normal vector nðxÞ, can be approximated on the boundary ∂Ω by NLs

ðnÞ qðxÞ≈qNL ðcð1Þ ; ξ; xÞ ¼ − ∑ cð1Þ n ½κ ∇x Fðx; ξ Þ  nðxÞ; s

 ∑

n¼1

ð25Þ

x∈∂Ω: NT

T

Nqc

collocation points, fxðNc þnÞ gn ¼ 1 ⊂Γ q , where NTc þ N qc ¼ N Lc , and collocate the boundary conditions (10b) and (10c) to obtain the following system of linear equations with respect to the unknown L

coefficients cð1Þ ∈RNs : ð1Þ

:

ð26Þ

Here Að11Þ ∈RNc Ns is the matrix whose elements are calculated L

L

ð1Þ

NLc

from Eqs. (24) and (25), while f ∈R is the right-hand side vector containing the corresponding discretised Dirichlet and Neumann data as given by Eqs. (10b) and (10c), respectively. Step 2. The MFS approximation for the particular solution of the non-homogeneous equilibrium equations (5) in R2 is given by   L γ 1−2ν Ns ð1Þ ð1Þ ∑ c ðy−ξðnÞ Þlog∥y−ξðnÞ ∥; ; ξ; yÞ ¼ − uðPÞ ðyÞ≈uðPÞ L ðc Ns 8πκG 1−ν n ¼ 1 n NLs

y∈R2 \ ⋃ fξðnÞ g:

ð27Þ

n¼1

" cð1Þ n

# 1  ðy−ξðnÞ Þ  nðyÞ ðnÞ ðnÞ log ∥y−ξ ∥ þ ν nðyÞ þ ðy−ξ Þ ; 1−2ν ∥y−ξðnÞ ∥2 ð28Þ

y∈∂Ω;

n¼1

Next, we select N Tc collocation points, fxðnÞ gn c¼ 1 ⊂Γ T , and N qc

Að11Þ cð1Þ ¼ f

NLs

ðPÞ

−γ TnÞ is approximated on ∂Ω by   γ 1−2ν ð1Þ ; ξ; yÞ−γ T tðPÞ ðyÞ−γ TðyÞnðyÞ≈tðPÞ L ðyÞnðyÞ ¼ L ðc Ns Ns 4πκ 1−ν " #  NLs ν ðy−ξðnÞ Þ  nðyÞ ðnÞ ðnÞ ð1Þ nðyÞ þ  ∑ cn log∥y−ξ ∥− ðy−ξ Þ ; 1−2ν ∥y−ξðnÞ ∥2 n¼1

while the term ðt

ð29Þ

y∈∂Ω: ð1Þ

NLs

Note that once the coefficients, c ∈R , corresponding to the thermal problem (10a)–(10c) are retrieved by solving Eq. (26), the particular solutions for the boundary displacement vector on Γ t and boundary traction vector on Γ u are expressed via Eqs. (27) and (28), respectively. Step 3. In the case of the Cauchy–Navier system associated with the two-dimensional isotropic linear elasticity, the fundamental solution matrix U ¼ ½Uij 1 ≤ i;j ≤ 2 , for the displacement vector is given by [1]
 1 y −η yj −ηj −ð3−4νÞ log; ∥y−η∥δij þ i i ; Uij ðy; ηÞ ¼ 8πGð1−νÞ ∥x−η∥ ∥x−η∥ y∈Ω; i; j ¼ 1; 2;

From approximation (27), the corresponding particular traction vector on the boundary ∂Ω is approximated as   γ 1−2ν tðPÞ ðyÞ≈tðPÞ ðcð1Þ ; ξ; yÞ ¼ − NLs 4πκ 1−ν

ð30Þ

where y ¼ ðy1 ; y2 Þ∈Ω is a collocation point and η ¼ ðη1 ; η2 Þ∈R2 \Ω is a singularity or source point. On differentiating Eq. (26) with respect to yk , k ¼ 1; 2, one obtains the derivatives of the fundamental solution for the displacement vector, denoted by ∂yk Uij ðy; ηÞ, where

1.0E-03

1.0E-03

1.0E-04

1.0E-04

1.0E-05 1.0E-05 1.0E-06 1.0E-06 1.0E-07 1.0E-08 1.0E-09 0.0

1.0E-07 M = N = 32 M = N = 48 0.1

0.2

M = N = 40 M = N = 56 0.3

M = N = 32 M = N = 48

0.4

1.0E-08 0.0

0.5

0.1

0.2

M = N = 40 M = N = 56 0.3

0.4

0.5

1.0E-04 1.0E-05 1.0E-06 1.0E-07 1.0E-08 1.0E-09

M = N = 32 M = N = 48

M = N = 40 M = N = 56

1.0E-10 0.0

0.1

0.2

0.3

0.4

0.5

Fig. 9. The pointwise normalised errors (a) errðt 1 ðxÞÞ, (b) errðt 2 ðxÞÞ, and (c) errðTðxÞÞ, x∈Γ q x∈fx∈R2 j∥x∥ ¼ R; θ∈½0; πÞg, obtained using d ¼ 0:50 and M ¼ N∈f32; 40; 48; 56g, for Example 2.

L. Marin, A. Karageorghis / Engineering Analysis with Boundary Elements 37 (2013) 1004–10205

∂yk ≡∂=∂yk . By combining Eqs. (8b) and (30), the fundamental h solution matrix T ¼ Tij 1 ≤ i;j ≤ 2 , for the traction vector in the case of two-dimensional isotropic linear elasticity is obtained as [1] T1k ðy; ηÞ ¼

2G ½ð1−νÞ ∂y1 U1k ðy; ηÞ þ ν ∂y2 U2k ðy; ηÞn1 ðyÞ 1−2ν

þG½∂y2 U1k ðy; ηÞ þ ∂y1 U2k ðy; ηÞn2 ðyÞ;

1013

solutions (31), namely NEs

tðHÞ ðyÞ≈tðHÞ ðcð2Þ ; η; yÞ ¼ ∑ Tðy; ηðnÞ Þcð2Þ n ; NE

ð33Þ

y∈∂Ω:

n¼1

s

By collocating the boundary conditions (12b) and (12c) at the Nu

y∈∂Ω; k ¼ 1; 2;

ð31aÞ

collocation points fyðnÞ gn c¼ 1 ⊂Γ u and fyðN

u

t þnÞ Nc gn ¼ 1 ⊂Γ t ,

respectively,

where N u þ N t ¼ N Ec , one obtains the following system of linear E

equations with respect to the unknown coefficients cð2Þ ∈R2Ns :

and T2k ðy; ηÞ ¼ G½∂y2 U1k ðy; ηÞ þ ∂y1 U2k ðy; ηÞn1 ðyÞ

A

2G ½ν∂y1 U1k ðy; ηÞ þ ð1−νÞ∂y2 U2k ðy; ηÞn2 ðyÞ; þ 1−2ν

ð32Þ

y∈Ω;

n¼1

s

ð2Þ ð2Þ T E 2 ð2Þ where cð2Þ ¼ ½ðc1ð2Þ ÞT ; …; ðcð2Þ ÞT T ∈R2Ns n ¼ ½cn;1 ; cn;2  ∈R , 1 ≤n ≤N s , c NE

E

s

E

¼f

and η∈R2Ns is a vector containing the coordinates of the singulaNE

rities fηðnÞ gn s¼ 1 . In a similar manner, the traction vector, tðHÞ , associated with the homogeneous equilibrium equation (12a) is approximated by a linear combination of the traction fundamental

c :

ð34Þ

E

ð2Þ

the matrix Að21Þ ∈R2Nc Ns are determined from those elements of the matrices that approximate uðPÞ ðyðnÞ Þ and ðtðPÞ −γ TnÞðyðnÞ Þ, E

L

1 ≤n ≤N Ec , according to Eqs. (12a), (12b), (27) and (29). E Step 4. Having determined the coefficients cð2Þ ∈R2Ns , the approximations of the boundary displacement, ujΓ t , and traction vectors, tjΓu , are obtained via the superposition principle and Eqs. (27), (29), (32) and (33). It is important to mention that, in order to uniquely determine L E the solutions cð1Þ ∈RNs and cð2Þ ∈R2Ns , the corresponding numbers of boundary collocation points and singularities must satisfy the inequalities N Ls ≤N Lc and N Es ≤N Ec , respectively. For exact boundary data, these systems can be solved by a direct method, such as the least-squares method or, equivalently, the inversion of the normal equation. However, in the case of perturbed boundary data, such

1.0E-04

1.0E-06

1.0E-06

1.0E-08

1.0E-08 M = N = 40 M = N = 56

M = N = 32 M = N = 48

1.0E-10

E

from Eqs. (32) and (33), f ∈R2Nc is the right-hand side vector containing the corresponding discretised Dirichlet and Neumann data as given by Eqs. (12b) and (12c), respectively. The elements of

1.0E-04

M = N = 32 M = N = 48

−A

ð21Þ ð1Þ

E

Analogously to the MFS approach for E the thermal problem, we N now consider N Es singularities, fηðnÞ gn s¼ 1 , and approximate the displacement vector, uðHÞ , associated with the homogeneous equilibrium equation (12a) in the solution domain by a linear combination of the displacement fundamental solutions (30) with respect to these singularities, i.e. NEs

c

ð2Þ

Here Að22Þ ∈R2Nc 2Ns is the matrix whose elements are calculated

y∈∂Ω; k ¼ 1; 2: ð31bÞ

uðHÞ ðyÞ≈uðHÞ ðcð2Þ ; η; yÞ ¼ ∑ Uðy; ηðnÞ Þcð2Þ n ; NE

ð22Þ ð2Þ

M = N = 40 M = N = 56

1.0E-10 0.5

0.6

0.7

0.8

0.9

1.0

0.5

0.6

M = N = 32 M = N = 48

M = N = 40 M = N = 56

0.7

0.8

0.9

1.0

1.0E-04 1.0E-05 1.0E-06 1.0E-07 1.0E-08 1.0E-09 1.0E-10 0.5

0.6

0.7

0.8

0.9

1.0

Fig. 10. The pointwise normalised errors (a) errðu1 ðxÞÞ, (b) errðu2 ðxÞÞ and (c) errðqðxÞÞ, x∈Γ t x∈fx∈R2 j∥x∥ ¼ R; θ∈½π; 2πÞg, obtained using d ¼ 0:50 and M ¼ N∈f32; 40; 48; 56g, for Example 2.

1014

L. Marin, A. Karageorghis / Engineering Analysis with Boundary Elements 37 (2013) 1004–1020

an approach would produce a highly unstable solution. Hence, in order to alleviate this problem, systems (26) and (34) should be solved by using e.g. the Tikhonov regularization method [31] or the singular value decomposition (SVD) [8].

5. Numerical results

qðanÞ ðxÞ ¼ −κ

1.0

1.0

0.5

0.5

0.0

-0.5

-1.0

-1.0 -0.5

0.0 x1

0.5

1.0 1.0

0.5

0.5

x2

1.0

0.0

logð∥x∥=Rint Þ logðRout =∥x∥Þ þ T int ; log ðRout =Rint Þ log ðRout =Rint Þ

T out −T int x  nðxÞ ; logðRout =Rint Þ ∥x∥2

x∈∂Ω;

x∈Ω;

ð35aÞ

ð35bÞ

0.0

-0.5

-1.0

x2

T ðanÞ ðxÞ ¼ T out

x2

x2

In this section, we investigate the performance of the MFS–MPS algorithm described in Section 3. This is done by solving the boundary value problem governed by the partial differential equations (5) and (10a), and subject to various choices of boundary conditions (8a), (8b), (10b) and (10c), for an isotropic linear thermoelastic material (copper alloy) characterised by the material constants G ¼ 4:80  1010 N=m2 , ν ¼ 0:34, κ ¼ 4:01 W m−1 K−1 , αT ¼ 16:5  10−6 1C−1 and γ ¼ 2Gα T ð1 þ νÞ=ð1−2νÞ.

Example 1. We consider the annular domain Ω ¼ fx∈R2 j Rint o ∥x∥ o Rout g, where Rint ¼ 1:0 and Rout ¼ 2:0, which is bounded by the inner and outer boundaries Γ int ¼ fx∈R2 j ∥x∥ ¼ Rint g and Γ out ¼ fx∈R2 j ∥x∥ ¼ Rout g, respectively. We also assume that the thermoelastic fields associated with the material occupying the domain Ω correspond to constant inner and outer temperatures, T int ¼ 1 1C and T out ¼ 2 1C, as well as constant inner and outer radial pressures, sint ¼ 1:0  1010 N=m2 and sout ¼ 2:0 1010 N=m2 , respectively, which describe a plane strain state. The boundary value problem (5), (8a), (8b) and (10a)–(10c) is solved for Γ T ¼ Γ t ¼ Γ int ∪Γ out and Γ q ¼ Γ u ¼ ∅, and admits the following analytical solution:

-1.0

-0.5

0.0 x1

0.5

-1.0

-0.5

0.0 x1

0.5

1.0

0.0

-0.5

-0.5

-1.0

-1.0 -1.0

-0.5

0.0 x1

0.5

1.0

1.0

1.0

x2

0.5

0.0

-0.5

-1.0 -1.0

-0.5

0.0 x1

0.5

1.0

Fig. 11. The analytical and numerical stress ðs11  10−10 ÞjΩ , obtained using d ¼ 0:50 and various numbers of boundary collocation points and singularities, namely (a) M ¼N ¼32, (b) M ¼N ¼ 40, (c) M ¼ N¼ 48, (d) M ¼ N¼ 56, (e) Analytical. M ¼ N∈f32; 40; 48; 56g, for Example 2.

L. Marin, A. Karageorghis / Engineering Analysis with Boundary Elements 37 (2013) 1004–10205


  γ 1−2ν T out −T int log∥x∥ 2 1−ν logðRout =Rint Þ    1−ν 1 x ; x∈Ω; −W þV 2 1þν ∥x∥ 2G

and

uðanÞ ðxÞ ¼

( tðanÞ ðxÞ ¼

sðHÞ ≡sint −γ T int þ int

−sout nðxÞ;

x∈Γ out ≡fx∈∂Ω j∥x∥ ¼ Rout g x∈Γ int ≡fx∈∂Ω j∥x∥ ¼ Rint g;

ð35dÞ

W≡

ðHÞ 2 2 ðsðHÞ out −sint Þ Rout Rint

R2out −R2int

γ T out −T int 2 logðRout =Rint Þ

x2

sðHÞ out ≡sout −γ T out þ

;

ð36aÞ

 1 log Rout þ 1 ; 1−ν

ð36bÞ

Example 3. We consider the same geometry as that given in Example 1, whilst the analytical solution (plane strain state) is

1.0

1.0

0.5

0.5

0.0

0.0

-0.5

-0.5

-1.0

-1.0 -0.5

0.0 x1

0.5

1.0

-1.0

1.0

1.0

0.5

0.5

x2

x2

;



-1.0

ð36cÞ

Example 2. We consider the unit disk Ω ¼ fx∈R2 j ∥x∥ oRg, R ¼ 1:0, and the analytical solution (plane strain state) given by Eqs. (13a), (13b) and (10a)–(10c), where K ¼ 1, xð1Þ ¼ ð2:0; 1:0Þ and T 1 ¼ 100 1C. Here, we solve the direct problem (5), (8a), (8b) and (10a)–(10c) with Γ q ¼ Γ u ¼ fx∈R2 j∥x∥ ¼ R; θ∈½0; πÞg and Γ T ¼ Γ t ¼ fx∈R2 j∥x∥ ¼ R; θ∈½π; 2πÞg, where θ ¼ θðxÞ is the radial angular polar coordinate associated with the point x∈R2 .

x2

R2out −R2int

 1 log Rint þ 1 : 1−ν

Some preliminary results for Example 1 have already been presented in [20].

−sint nðxÞ;

ðHÞ 2 2 sðHÞ out Rout −sint Rint



ð35cÞ

where V≡−

γ T out −T int 2 logðRout =Rint Þ

1015

0.0

-0.5

-0.5

0.0 x1

0.5

1.0

0.0

-0.5

-1.0

-1.0 -1.0

-0.5

0.0 x1

0.5

1.0

-1.0

-0.5

-1.0

-0.5

0.0 x1

0.5

1.0

1.0

x2

0.5

0.0

-0.5

-1.0 0.0 x1

0.5

1.0

Fig. 12. The analytical and numerical displacement u2 jΩ , obtained using d ¼ 0:50 and various numbers of boundary collocation points and singularities, namely (a) M ¼N ¼32, (b) M ¼ N¼ 40, (c) M ¼ N¼ 48, (d) M ¼ N¼ 56, (e) Analytical. M ¼ N∈f32; 40; 48; 56g, for Example 2.

1016

L. Marin, A. Karageorghis / Engineering Analysis with Boundary Elements 37 (2013) 1004–1020

described by Eqs. (13a), (13b) and (10a)–(10c), where K ¼ 2, xð1Þ ¼ ð2:5; 2:5Þ, xð2Þ ¼ ð−0:2; 0:2Þ, T 1 ¼ 100 1C and T 2 ¼ 50 1C. Here, we solve the direct problem (5), (8a), (8a) and (10a)– (10c) with Γ T ¼ Γ t ¼ Γ int and Γ q ¼ Γ t ¼ Γ out . For all examples analysed herein, we have considered N Lc ¼ N Ec ¼ N uniformly distributed collocation points on ∂Ω, as well as N Ls ¼ N Es ¼ M uniformly distributed singularities, which are preassigned and kept fixed throughout the solution process (i.e. the socalled static approach has been employed) on a pseudo-boundary ~ of a similar shape to that of ∂Ω such that distð∂Ω; ~ ∂ΩÞ is a fixed ∂Ω constant, see e.g. [7]. According to the notations used in Section 4, the corresponding MFS parameters have been set as follows:

ðnumÞ

respectively, where f ðxÞ denotes an approximate numerical value for f ðxÞ, x∈∂Ω. Analogously, one can define the pointwise normalised error of f at x∈Ω and the maximum normalised error of f inside the domain Ω.

0.5

0.0

-0.5

-0.5

-1.0

-1.0 0.0 x1

0.5

1.0

1.0

1.0

0.5

0.5

0.0

-0.5

-1.0

-1.0 -0.5

0.0 x1

0.5

-1.0

-0.5

0.0 x1

0.5

1.0

-1.0

-0.5

0.0 x1

0.5

1.0

0.0

-0.5

-1.0

ð37bÞ

x∈∂Ω

x2

x2

Errðf Þ ¼ max errðf ðxÞÞ;

0.5

-0.5

~ ¼ ∂Ω

and

1.0

-1.0

on

y∈∂Ω

1.0

0.0

M¼N

In order to assess the accuracy and convergence of the proposed MFS–MPS approach, for any real-valued function f : ∂Ω⟶R, we define the corresponding pointwise normalised error of f at x∈∂Ω and the maximum normalised error of f on ∂Ω by ðnumÞ jf ðxÞ−f ðxÞj errðf ðxÞÞ ¼ ; x∈∂Ω; ð37aÞ maxjf ðyÞj

x2

x2

(i) Examples 1 and 3: N∈f48; 56; 72; 84g with ð2N=3Þ on Γ out and ðN=3Þ on Γ int ; M ¼ N with ð2M=3Þ on Γ~ out ¼ fx∈R2 j ∥x∥ ¼ Rout þ ~ ¼ d1 g and ðM=3Þ on Γ~ int ¼ fx∈R2 j ∥x∥ ¼ Rint −d2 g, where ∂Ω Γ~ out ∪Γ~ int , d1 ∈½0:05; 1:0 and d2 ∈½0:02; 0:90.

(ii) Example2: N∈f32; 40; 48; 56g on ∂Ω; fx∈R2 j ∥x∥ ¼ R þ dg, where d∈½0:05; 5:0.

1.0

1.0

x2

0.5

0.0

-0.5

-1.0 -1.0

-0.5

0.0 x1

0.5

1.0

Fig. 13. The analytical and numerical temperature TjΩ , obtained using d ¼ 0:50 and various numbers of boundary collocation points and singularities, namely (a) M¼ N¼ 32, (b) M¼ N¼ 40, (c) M¼ N¼ 48, (d) M¼ N¼ 56, (e) Analytical. M ¼ N∈f32; 40; 48; 56g, for Example 2.

L. Marin, A. Karageorghis / Engineering Analysis with Boundary Elements 37 (2013) 1004–10205

Fig. 1(a)–(d) presents the numerical results for the x1  component of the displacement on the boundary ∂Ω, in terms of the maximum normalised error Errðu1 Þ defined by Eq. (37b), obtained using the proposed MFS–MPS approach, various numbers of boundary collocation points and singularities, i.e. M ¼ N∈f48; 60; 72; 84g, d1 ¼ dist ðΓ~ out ; Γ out Þ∈½0:05; 1:0 and d2 ¼ distðΓ~ int ; Γ int Þ∈½0:02; 0:90, for Example 1. From these figures we observe that the numerical solution for the displacement u1 j∂Ω is accurate and converges to its corresponding exact solution as we increase the number of collocation points, as well as the distances of the boundaries Γ out and Γ int to the pseudo-boundaries Γ~ out and Γ~ int , respectively. Although not presented, similar results were obtained for the x2 component of the displacement and the traction vector on ∂Ω, as well as the displacement vector and stress tensor in the domain Ω. Convergent and accurate numerical results have also been retrieved for the thermal fields in Ω with M ¼ N∈f48; 60; 72; 84g, d1 ∈½0:05; 1:0 and d2 ∈½0:02; 0:9, as can be seen from Fig. 2(a)– (d) which displays the maximum normalised error of the normal heat flux on ∂Ω, ErrðqÞ. Furthermore, it can be observed that, as the numbers of boundary collocation points M and singularities N increase, d1 increases and d2 ⟶0:90, the distribution of the maximum normalised error ErrðqÞ becomes irregular. This is a direct consequence of the fact that the corresponding MFS matrix becomes highly ill-conditioned since not only the dimension of this matrix increases (i.e. M ¼ N increase), but also the singularities on the pseudo-boundary Γ~ int associated with the inner boundary Γ int cluster around the origin. Although not illustrated, it is worth mentioning that similar results were obtained in Examples 2 and 3.

Figs. 3 and 4 present the analytical and numerical results for the unknown boundary displacements and normal heat fluxes on the inner (u1 jΓint , u2 jΓint and qjΓint ) and outer (u1 jΓout , u2 jΓout and qjΓ out ) boundaries, respectively, obtained using d1 ¼ 1:0, d2 ¼ 0:30 and various numbers of boundary collocation points and singularities, namely M ¼ N∈f48; 60; 72; 84g, for Example 1. From these figures we observe that the numerical results retrieved for the unknown data on both boundaries Γ int and Γ out are very accurate approximations to their corresponding exact values. Moreover, these numerical approximations are, as expected, convergent with respect to increasing the numbers of boundary collocation points and singularities, as can be observed from Figs. 5 and 6 which illustrate the convergence of the MFS–MPS algorithm in terms of the associated pointwise normalised error (37a) on Γ int and Γ out , respectively, for Example 1. Similar conclusions regarding the convergence of the proposed MFS–MPS procedure can also be drawn from Figs. 7 and 8, which present the numerically retrieved boundary data on Γ u ¼ Γ q ¼ fx∈R2 j ∥x∥ ¼ R; θ∈½0; πÞg and Γ t ¼ Γ T ¼ fx∈R2 j ∥x∥ ¼ R; θ∈½π; 2πÞg, respectively, using d ¼ 0:50 and M ¼ N∈f32; 40; 48; 56g, in comparison with their corresponding analytical values, for Example 2. To get a better insight into the convergence process in the case of Example 2, the pointwise normalised error (37a) associated with the boundary traction vector and boundary temperature on Γ u ¼ Γ q , and boundary displacement vector and normal heat flux on Γ t ¼ Γ T , obtained using the same values of the MFS–MPS approach, are displayed in Figs. 9(a)–(c) and 10(a)–(c), respectively. Furthermore, a comparison between the exact and numerically retrieved for the stress component s11 , displacement u2 and temperature T inside the solution domain Ω,

0.1

0.1 Analytical M = N = 48 M = N = 60 M = N = 72 M = N = 84

0.0

0.05

t2/1010

0.05

t1/1010

1017

0.0

-0.05

-0.05

-0.1

-0.1

-0.15

Analytical M = N = 48 M = N = 60 M = N = 72 M = N = 84

-0.15 -1.0

-0.8

-0.6

-0.4

-0.2

0.0

-1.0

-0.8

-0.6

-0.4

-0.2

0.0

0 Analytical M = N = 48 M = N = 60 M = N = 72 M = N = 84

q

-100

-200

-300

-400 -1.0

-0.8

-0.6

-0.4

-0.2

0.0

Fig. 14. The analytical and numerical tractions (a) t 1 jΓ int , (b) t 2 jΓ int , and (c) normal heat flux qjΓ int , obtained using d1 ¼ 1:0, d2 ¼ 0:30 and M ¼ N∈f48; 60; 72; 84g, for Example 3.

1018

L. Marin, A. Karageorghis / Engineering Analysis with Boundary Elements 37 (2013) 1004–1020

0.0

Analytical M = N = 48 M = N = 60 M = N = 72 M = N = 84

-0.005

u2

u1

0.0

-0.01

-0.005

Analytical M = N = 48 M = N = 60 M = N = 72 M = N = 84

-0.01

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

200

150 T Analytical M = N = 48 M = N = 60 M = N = 72 M = N = 84

100

0.0

0.2

0.4

0.6

0.8

1.0

Fig. 15. The analytical and numerical displacements (a) Example 1: u1 jrout t 1 jΓ out , (b) Example 1: u2 jrout t 2 jΓ out , and (c) Example 1: temperature TjΓ out , obtained using d1 ¼ 1:0, d2 ¼ 0:30 and M ¼ N∈f48; 60; 72; 84g, for Example 3.

1.0E-02

1.0E-02

1.0E-03

1.0E-03

1.0E-04 1.0E-04 1.0E-05 1.0E-05 1.0E-06 1.0E-06

1.0E-07

M = N = 48 M = N = 72

M = N = 60 M = N = 84

M = N = 48 M = N = 72

1.0E-08

M = N = 60 M = N = 84

1.0E-07 -1.0

-0.8

-0.6

-0.4

-0.2

0.0

-1.0

-0.8

M = N = 48 M = N = 72

M = N = 60 M = N = 84

-0.6

-0.4

-0.2

0.0

1.0E-02 1.0E-03 1.0E-04 1.0E-05 1.0E-06 1.0E-07 -1.0

-0.8

-0.6

-0.4

-0.2

0.0

Fig. 16. The pointwise normalised errors Example 1: (a) errðt 1 ðxÞÞ, (b) errðt 2 ðxÞÞ, and (c) errðqðxÞÞ, x∈Γ int , obtained using d1 ¼ 1:0, d2 ¼ 0:30 and M ¼ N∈f48; 60; 72; 84g, for Example 3.

L. Marin, A. Karageorghis / Engineering Analysis with Boundary Elements 37 (2013) 1004–10205

1019

1.0E-04

1.0E-04

1.0E-05

1.0E-05

1.0E-06 1.0E-06 1.0E-07 1.0E-07 1.0E-08 1.0E-08 M = N = 48 M = N = 72

1.0E-09 0.0

0.2

0.4

1.0E-09

M = N = 60 M = N = 84 0.6

0.8

1.0E-10

1.0

M = N = 48 M = N = 72 0.0

0.2

0.4

M = N = 60 M = N = 84 0.6

0.8

1.0

1.0E-03 M = N = 48 M = N = 72

M = N = 60 M = N = 84

1.0E-04

1.0E-05

1.0E-06

1.0E-07

0.0

0.2

0.4

0.6

0.8

1.0

Fig. 17. The pointwise normalised errors Example 1: (a) errðu1 ðxÞÞ, (b) errðu2 ðxÞÞ, and (c) errðTðxÞÞ, x∈Γ out , obtained using d1 ¼ 1:0, d2 ¼ 0:30 and M ¼ N∈f48; 60; 72; 84g, for Example 3.

for Example 2, are presented in Figs. 11, 12 and 13, respectively. In addition, it can be seen from these figures that very accurate approximations for the internal stress, displacement vector and temperature field are obtained even for relatively small numbers of boundary collocation points and singularities, say M ¼ N ¼ 32. Finally, a comparison between the analytical and numerical results for the unknown boundary data, obtained using d1 ¼ 1:0, d2 ¼ 0:30 and M ¼ N∈f48; 60; 72; 84g, in the case of Example 3, is presented in Figs. 14(a)–(c) and 15(a)–(c). The corresponding pointwise normalised errors on inner (errðu1 ðxÞÞ, errðu2 ðxÞÞ and errðqðxÞÞ, x∈Γ int ) and outer (errðt 1 ðxÞÞ, errðt 2 ðxÞÞ and errðTðxÞÞ, x∈Γ out ) boundaries are displayed in Figs. 16(a)–(c) and 17(a)–(c), respectively. Overall, from the examples investigated, we can conclude that the proposed MFS–MPS algorithm provides accurate and convergent numerical approximations with respect to the numbers of collocation points and singularities, for both the unknown boundary data and the internal thermal and mechanical fields, for various direct problems associated with two-dimensional planar thermoelastic materials.

homogeneous equations of equilibrium which only depends on the MFS approximation of the boundary value problem for the Laplace equation. Furthermore, the development of such a particular solution enables one to easily develop analytical solutions in the framework of two-dimensional isotropic linear thermoelasticity, regardless of the geometry of the problem in question. The accuracy and convergence properties of the proposed MFS–MPS procedure were investigated by considering three numerical examples in both simply and doubly connected domains. Future work is related to the development of fast MFS– MPS algorithms for two-dimensional isotropic thermoelasticity [14], as well as the application of the proposed MFS–MPS procedure for the stable numerical solution of inverse boundary value problems in planar isotropic thermoelasticity, see [21] for preliminary results. In addition, the suitability of the proposed approach in relation to the recent work of Trevelyan and his co-workers on the enriched BEM [27,32] as well as on the isogeometric BEM [26], which is based on the corresponding work using the FEM presented in [10], could be investigated.

6. Conclusions Acknowledgements The solution of boundary value problems in two-dimensional linear isotropic thermoelasticity was investigated using the MFS in conjunction with the MPS. This approach is based on the construction of a novel particular solution of the non-

The financial support received from the Romanian National Authority for Scientific Research (CNCS-UEFISCDI), Project no. PNII-ID-PCE-2011-3-0521, is gratefully acknowledged.

1020

L. Marin, A. Karageorghis / Engineering Analysis with Boundary Elements 37 (2013) 1004–1020

References [1] Aliabadi MH. The boundary element method. Applications in solids and structures, vol. 2. London: John Wiley & Sons; 2002. [2] Cheng AH-D, Chen CS, Golberg MA, Rashed YF. An adaptive boundary element scheme for steady thermoelastic analysis. Eng Anal Boundary Elem 2001;25: 377–87. [3] Dennis BH, Dulikravich GS. A finite element formulation for the detection of boundary conditions in elasticity and heat conduction. In: Tanaka M, Dulikravich GS, editors. Inverse problems in engineering mechanics. UK: Elsevier Science; 1998. p. 61–70. [4] Dennis BH, Dulikravich GS. Simultaneous determination of temperatures, heat fluxes, deformations, and tractions on inaccessible boundaries. J Heat Transf Trans ASME 1999;121:537–45. [5] Fairweather G, Karageorghis A. The method of fundamental solutions for elliptic boundary value problems. Adv Comput Math 1998;9:69–95. [6] Fairweather G, Karageorghis A, Martin PA. The method of fundamental solutions for scattering and radiation problems. Eng Anal Boundary Elem 2003;27:759–69. [7] Gorzelańczyk P, Kołodziej JA. Some remarks concerning the shape of the shape contour with application of the method of fundamental solutions to elastic torsion of prismatic rods. Eng Anal Boundary Elem 2008;32:64–75. [8] Hansen PC. Rank-deficient and discrete ill-posed problems: numerical aspects of linear inversion. Philadelphia: SIAM; 1998. [9] Henry Jr DP, Banerjee PK. A new boundary element formulation for two- and three-dimensional thermoelasticity using particular integrals. Int J Numer Methods Eng 1988;26:2061–77. [10] Hughes TJR, Cottrell JA, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Comput Methods Appl Mech Eng 2005;194:4135–95. [11] Kamiya N, Aikawa Y, Kawaguchi K. An adaptive boundary element scheme for steady thermoelastic analysis. Comput Methods Appl Mech Eng 1994;119: 311–24. [12] Karageorghis A. Modified methods of fundamental solutions for harmonic and biharmonic problems with boundary singularities. Numer Methods Partial Differential Equations 1992;8:1–19. [13] Karageorghis A, Lesnic D, Marin L. A survey of applications of the MFS to inverse problems. Inverse Probl Sci Eng 2011;19:309–36. [14] Karageorghis A, Marin L, Efficient MFS algorithms for problems in thermoelasticity. J Sci Comput, in press; http://dx.doi.org/10.1007/s10915-012-9664-x. [15] Karageorghis A, Poullikkas A, Berger JR. Stress intensity factor computation using the method of fundamental solutions. Comput Mech 2006;37:445–54. [16] Karageorghis A, Smyrlis Y-S. Matrix decomposition MFS algorithms for elasticity and thermo-elasticity problems in axisymmetric domains. J Comput Appl Math 2007;206:774–95. [17] Marin L. Stable MFS solution to singular direct and inverse problems associated with the Laplace equation subjected to noisy data. CMES: Comput Model Eng Sci 2008;37:203–42.

[18] Marin L. Treatment of singularities in the method of fundamental solutions for two-dimensional Helmholtz-type equations. Appl Math Model 2010;34: 1615–33. [19] Marin L. A meshless method for the stable solution of singular inverse problems for two-dimensional Helmholtz-type equations. Eng Anal Boundary Elem 2010;34:274–88. [20] Marin L, Karageorghis A. MFS-based solution to two-dimensional linear thermoelasticity problems. In: Brebbia CA, Poljak D, editors. Boundary elements and other mesh reduction methods XXXIV (BEM/MRM 2012). Southampton, UK: WIT Press; 2012. p. 39–49. [21] Marin L, Karageorghis A. MFS solution of inverse boundary value problems in two-dimensional linear thermoelasticity. In: Prochazka P, Aliabadi MH, editors. Advances in boundary element techniques XIII: international conference on boundary element and meshless XII. Techniques. UK: EC Ltd.; 2012. p. 141–6. [22] Mathon R, Johnston RL. The approximate solution of elliptic boundary value problems by fundamental solutions. SIAM J Numer Anal 1977;14:638–50. [23] Partridge PW, Brebbia CA, Wrobel LC. The dual reciprocity boundary element method. Southampton: Computational Mechanics Publications; 1992. [24] Rizzo FJ, Shippy DJ. An advanced boundary integral equation method for three-dimensional thermoelasticity. Int J Numer Methods Eng 1977;11: 1753–68. [25] Rizzo FJ, Shippy DJ. The boundary element method in thermoelasticity. In: Banerjee PK, Butterfield R, editors. Developments in boundary element methods—I. London: Applied Science Publishers; 1979. [26] Simpson RN, Bordas SPA, Trevelyan J, Rabczuk T. A two-dimensional isogeometric boundary element method for elastostatic analysis. Comput Methods Appl Mech Eng 2012;209–212:87–100. [27] Simpson R, Trevelyan J. Evaluation of J1 and J2 integrals for curved cracks using an enriched boundary element method. Eng Fract Mech 2011;78:623–37. [28] Sladek V, Sladek J. Boundary integral equation in thermoelasticity. Part I: general analysis. Appl Math Model 1983;7:413–8. [29] Sladek V, Sladek J. Boundary integral equation in thermoelasticity. Part III: uncoupled thermoelasticity. Appl Math Model 1984;8:413–8. [30] Sladek J, Sladek V, Atluri SN. A pure contour formulation for the meshless local boundary integral equation method in thermoelasticity. CMES—Comput Model Eng Sci 2001;2:423–33. [31] Tikhonov AN, Arsenin VY. Methods for solving ill-posed problems. Moscow: Nauka; 1986. [32] Trevelyan J, Simpson R. Enriched dual BEM for curved cracks. In: Lesnic D, editor. Proceedings of the 8th UK conference on boundary integral methods. Leeds, UK: Leeds University Press; 2011. p. 1–8. [33] Tsai CC. The method of fundamental solutions with dual reciprocity for threedimensional thermoelasticity under arbitrary forces. Int J Comput Aided Eng Softw 2009;26:229–44.