Local radial basis function collocation method for bending analyses of quasicrystal plates

Accepted Manuscript

Local radial basis function collocation method for bending analyses of quasicrystal plates Y.C. Chiang , D.L. Young , J. Sladek , V. Sladek PII: DOI: Reference:

S0307-904X(17)30391-8 10.1016/j.apm.2017.05.051 APM 11803

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

6 November 2016 16 May 2017 31 May 2017

Please cite this article as: Y.C. Chiang , D.L. Young , J. Sladek , V. Sladek , Local radial basis function collocation method for bending analyses of quasicrystal plates, Applied Mathematical Modelling (2017), doi: 10.1016/j.apm.2017.05.051

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights

AC

CE

PT

ED

M

AN US

  

Meshless formulations are developed for static and dynamic bending of QC plates. The phonon-phason coupling effects are studied numerically. The influence of elastic foundation and variable plate thickness is investigated. The response of QC plates to impact loading is simulated.

CR IP T



1

ACCEPTED MANUSCRIPT

Local radial basis function collocation method for bending analyses of quasicrystal plates

CR IP T

Y.C. Chiang1,2, D. L. Young*1,2, J. Sladek3 and V.Sladek3

AN US

Abstract The local radial basis function collocation method (LRBFCM) is proposed for plate bending analysis in orthorhombic quasicrystals (QCs) under static and transient dynamic loads. Three common types of the plate bending problems are considered: (1) QC plates resting on Winkler foundation (2) QC plates with variable thickness and (3) QC plates under a transient dynamic load. According to the Reissner–Mindlin plate bending theory, there is allowed to simulate the behavior of the two excitations in QC plates, phonon and phason, by 2D strong formulations for the system of governing equations. The governing equations, which describe the phason displacements, are based on Agiasofitou and Lazar elastodynamic model. Numerical results demonstrate the effect of the elastic foundation, as well as plate thickness on the phonon and phason characteristics in this paper. For the transient dynamic analysis, the influence of the phason friction coefficients on the responses of QC plate to transient dynamic loads is also studied.

Introduction

ED

1

M

Keywords: local radial basis function collocation method, orthorhombic quasicrystal, Reissner-Mindlin theory, phonon and phason displacements

AC

CE

PT

The quasicrystals (QCs) are special solid structures with long-range orientational order but no translational symmetry. That is, the quasicrystalline patterns of its atomic arrangement have the ability to fill all available space, but the arrangement is aperiodic or never repeats. The icosahedral QCs were first discovered and proposed by Shechtman et al. [1]. However, people did not believe that these controversial atomic arrangements could exist at that time. Until recent years, the extraordinary properties of QCs have caught more and more scientists’ attentions, and they eventually confirmed the existence of QCs. Shechtman was thereby awarded the Nobel Prize in Chemistry in 2011 for his remarkable works on QCs. Generally speaking, QCs exist in many polymers and metallic alloys and most often are discovered in aluminum alloys. These materials are widely applied on large numbers of engineering devices. Nowadays, QCs have become an important and potential field for researchers from several areas. Elasticity is one of the 1

Department of Civil Engineering, National Taiwan University, 10617 Taipei, Taiwan. CoreTech (Moldex3D), 30265 Chupei, Taiwan. 3 Institute of Construction and Architecture, Slovak Academy of Sciences, 84503 Bratislava, Slovakia. * Corresponding author, Fax/phone: +886-2-2362-6114, e-mail: [email protected]. 2

2

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

important properties that are concerned in this study. The elastic behavior of QCs is significantly different from the conventional crystals. Numerous physicists have made many efforts and progresses on the researches about the elasticity of QCs. The Landau density wave theory [2] is recognized as a physical basis of the elastic theory for QCs. Unlike conventional elastic theory for crystals which describes only phonon degrees of freedom, the quasi-periodicity in the QCs leads to additional degrees of freedom via phason displacements. Bak [3] pointed out that the phason structure disorders are described by the fluctuations in QCs. In view of the Bak’s argument, the generalized theories of elasticity of QCs are proposed [4-5]. Both phonon and phason play similar roles on the dynamics and the governing equations of the motion are both described as wave propagations. Lubensky et al. [6] considered that the phason displacements are insensitive to spatial translation, so that they claimed that the governing equations of the phason should be described by diffusive type with large diffusion time. Rochal and Lorman [7] proposed the minimal model of the phonon-phason dynamics in icosahedral QCs. Based on Rochal and Lorman works, the elasto-hydrodynamic model of QCs provided by Fan et al. [8] has been applied on many researches. Recently, Agiasofitou and Lazar [9] introduced the elastodynamic model of wave-telegraph type equations for the description of dynamics of quasicrystals. Phonons are represented by waves, and phasons by waves damped in time and propagating with finite velocity; that means the equations of motion for the phonons are partial differential equations of the wave type, and for the phasons partial differential equations of the telegraph type. This model considers the effect of the dynamic behavior of phasons through the phason friction. Some studies on the mathematical theory of elasticity of QCs and some relative applications could be found in the book by Fan [10]. Fan [11] has reviewed recent progress on the mathematical theory and mechanics foundation on QCs, such as elasticity, plasticity, dynamics and other subjects. The orthorhombic QCs are classified as one dimensional (1D) QCs. On the other words, the atom arrangements of 1D QCs are quasi-periodic only in one direction and the other two are periodic. There are a few papers theoretically analyzed this kind of materials in recent years. Li and Liu [12] presented the physical property tensors for 1D QCs. Chen et al. [13] analyzed the theory of 3D elastic problems with 1D hexagonal QC bodies. Wang and Pan [14] derived the analytical solutions for some defect problems in 1D hexagonal and 2D octagonal QCs. Gao et al. [15] developed the general solutions for the 1D QCs plane elasticity problems. Researchers also paid attention to the theory of QCs beam or plate structures. Analysis of plate bending in 1D hexagonal QCs was proposed by Gao et al. [15], and then Gao and Ricoeur [16] refined the theory of 1D QCs in thick plate structures. Moreover, Gao [17] pointed out the exact theory of 1D QC deep beams. In recent years, researchers have turned to solve the beam and plate bending problems numerically. Reliable numerical methods are required to analyze problems with more complicated geometry and boundary conditions. Mesh-dependent numerical methods, such as finite difference method (FDM) [18], finite element method (FEM), and boundary element method (BEM) have been widely developed for solving boundary value problems in several different engineering areas. However, mesh-dependent numerical methods still remain inherent challenges; it is time-consuming for the process of mesh generation or too tedious for the numerical integration when researchers conduct

3

ACCEPTED MANUSCRIPT

M

AN US

CR IP T

with complex geometry problems. Recently, meshless methods become potential alternatives to the conventional mesh-dependent methods [19]. Due to the high flexibility and low cost, meshless numerical methods spring up for researchers to deal with various complex engineering problems [20-28]. There are several papers in which the plate bending problems by meshless methods are analyzed. Sladek et al. [29] applied the moving least-squares (MLS) method with local boundary integral equations (LBIEs) to the simply-supported and clamped plate bending problems. The MLS was also employed to analyze the plates resting on elastic foundation [30]. Sator et al. [31-33] developed the unified formulation for the elastic plates which includes the Kirchhoff-Love theory as well as the 1st and 3rd order shear deformation plate theory for studying the coupling effects in elastic analysis of functionally graded plates and plates with variable thickness by the mesh-free methods. In terms of QC materials, the meshless local Petrov-Galerkin method (MLPG) based on the MLS scheme was applied to analyze orthorhombic QC plates and shells under static and transient dynamic loads [34, 35]. Another promising meshless method is the method of fundamental solutions (MFS). Recently Yaslan [36] used the method of time-dependent fundamental solutions to solve the 3D quasicrystal problems. The influence of the coupling parameter of QC plate and shell curvature of QC shell on phonon and phason displacements are investigated therein, respectively. Because the QCs received much more attention in recent years, there are worthy to develop and apply advanced numerical methods to analyze boundary value problems in QC. In order to broaden the field of meshless numerical methods for QC problems, this paper first attempts to utilize the LRBFCM for the plate bending analysis in orthorhombic QC problems.

AC

CE

PT

ED

The LRBFCM is a local type meshless method, which employs the interpolation of radial basis functions (RBFs) associated to the local supported nodes to approximate the unknown physical variables. In contrast to the conventional mesh-dependent numerical methods, LRBFCM focuses only on nodes within arbitrary computational domain; it is no requirement for the implementation of mesh and integration of governing equations. The RBFs have been widely applied on numerical approximation for the partial differential equations. One of the most popular RBFs, multiquadric (MQ) RBF, was proposed by Hardy [37]. Frankle [38] made efforts on investigating kinds of RBFs by scattered data interpolation and verified the sensitivity and accuracy of the MQRBF. Kansa [39, 40] first introduced the radial basis function collocation method (RBFCM), which belongs to global type meshless method and based on the MQRBF. Due to the easy implementation, it is convenient for the RBFCM to solve the partial differential equations (PDEs) numerically [41, 42]. However, the global type meshless methods, which take account of all computational nodes within whole analyzed domain, tend to consume time for solving the dense interpolation matrix. Additionally, it is more probable to cause the interpolation matrix ill-conditioning and large memory loading for the computation. In order to overcome these drawbacks, the localized scheme is applied to improve the efficiency of RBFCM, which is the so-called LRBFCM. The selection of the supported local nodes and the optimal shape parameter c play important roles on LRBFCM. Shen et al.[43- 45] discussed the detail about the two major selection issues by the local radial basis function differential quadrature method (LRBFDQ), which is also a

4

ACCEPTED MANUSCRIPT

CR IP T

local type RBF method. Several papers have made efforts on the selection of the shape parameter. Rippa [46] proposed a well-known principle, which is called leave-one-out cross validation (LOOCV). According to the criteria of LOOCV, Uddin [47] proposed an extended algorithm for time-dependent PDEs. Gherlone et al. [48] provided a novel algorithm which is on the basis of convergence analysis and demonstrated the effectiveness on static analysis of laminated composite and sandwich plate problems. Iurlaro et al. [49] proposed a selection method which is based on the principle of the minimum of the total potential energy. Recent developments to avoid the difficulty to choose optimal shape parameter is to consider the condition number results of the interpolation matrix, the concept has been successfully implemented by Mavrič et al.[50] on thermo-elasticity problems.

PT

ED

M

AN US

With the promotion of the accuracy and the stability, the LRBFCM has been broadly applied on solving some kinds of engineering problems in recent years, such as diffusion problems [51, 52], shallow water equations [53], thermo-driven fluid flow problems [54], thermo-mechanical problems [55], hyperbolic partial differential equations [56], timefractional heat diffusion in functionally graded materials [57] and functionally graded piezoelectric semi-conductor analysis problems [58]. In this paper, the LRBFCM is first applied on plate bending analysis in orthorhombic QC problems. According to the Reissner–Mindlin plate bending theory [59, 60], we could utilize 2D governing equations of the motion for the phonon and phason in the mid-plane to analyze the phonon and phason displacements which are independent on the transversal coordinate. The normalized technique for the shape parameter of the MQRBF will be introduced in order to improve the stability of the numerical solutions. Two static study cases for the QC plates under static uniform load; (1) the plates resting on elastic foundation and (2) plates with variable thickness are investigated, respectively. Moreover, the case of the QC plates under transient dynamic loads is also analyzed. The time scheme of the second order resulting from the equations of motion is solved by the Houbolt higher order finite difference method [61]. Numerical results are presented for different values of the elastic foundation parameters, variable plate thickness, and phason friction coefficients for each case in order to analyze the behavior of the phonon and phason displacements.

AC

CE

The organization of this paper is as follows: In the section 2, we describe the governing equations based on the Reissner-Mindlin theory and two types of boundary conditions for the orthorhombic QC plates. Numerical implementations of the LRBFCM are introduced thoroughly in the section 3. Then, we present numerical examples for the proposed numerical method in the section 4 and conclude the present study in the section 5. 2

Governing equations and boundary conditions for the orthorhombic QC plates

2.1 Governing equations Let us consider a QC plate with constant thickness h and homogeneous material properties, as shown in Fig.1. The mean surface of the QC plate within the domain  in

5

ACCEPTED MANUSCRIPT

q(x ,t)

x3 x 2 x1

h

Q1

CR IP T

the Cartesian plane x  ( x1, x2 ) and the third coordinate x3  z is perpendicular to the mid-plane where placed x3  0 of the domain. Thus, x3  h / 2 and x3  h / 2 denotes the top and bottom surfaces of the QC plate, respectively.

AN US

M11

M22

M21 Q2

M12

M

Figure 1: Sign convention of bending moments and forces for a QC plate

ED

The Reissner–Mindlin plate bending theory [59, 60] is used to take the QC plate deformation into account. By the concept of the Reissner–Mindlin plate bending theory, we can transfer the original 3D plate into the 2D one in the mid plane of the plate with constant shear strain throughout the plate thickness. The phonon displacements components (u1 , u2 , u3 ) of the QC plate can be expressed as

PT

u1 (x, x3 , t )  u10 (x, t )  x3 1 (x, t )

u2 (x, x3 , t )  u20 (x, t )  x3 2 (x, t )

(1)

u3 (x, t )  u30 (x, t )

CE

where t is the time variable,  1 (x, t ) ,  2 (x, t ) represent the rotations in x1 direction and

AC

x2 direction around the in-plane axes. u10 (x, t ) , u20 (x, t ) denote in-plane displacements in x1 direction and x2 direction, respectively. u30 (x, t ) represents the out-of-plane displacements by the assumption that u3 is independent on x3 . The linear phonon strains  ij (x, x3 , t ) are symmetric and relative with the phonon displacements as can be seen in equation (2) 1  ui u j  2  x j xi

 ij (x, x3 , t )   

6

  . 

(2)

ACCEPTED MANUSCRIPT

According to the Reissner–Mindlin plate bending theory in equation (1), we can obtain the linear phonon strains by the following derivations. u (x, x3 , t ) u10 (x, t )  1 (x, t ) 11 (x, x3 , t )  1   x3 (3) x1 x1 x1 u2 (x, x3 , t ) u20 (x, t )  2 (x, t )   x3 x2 x2 x2

 33 (x, x3 , t ) 

u3 (x, t ) 0 x3

(4)

CR IP T

 22 (x, x3 , t ) 

(5)

1  u1 (x, x3 , t ) u2 ( x, x3 , t )    2 x2 x1 

12 (x, x3 , t )  

AN US

1  u (x, t ) u20 (x, t )  1   1 (x, t )  2 (x, t )    10     x3   2  x2 x1  2  x2 x1 

(6)

u30 (x, t )  1  u1 (x, x3 , t ) u3 (x, t )  1    =  1 (x, t )   2 x3 x1  2  x1 

(7)

u30 (x, t )  1  u2 (x, x3 , t ) u3 (x, t )  1      2 (x, t )   2 x3 x2  2  x2 

(8)

13 (x, t )  

M

 23 (x, t )  

The generalized Hooke’s law for plane elasticity of orthorhombic QC is given as [10] (9)

 22  c1211  c22 22  c23 33  R2 w33

(10)

 31  2c55 31  R6 w31

(11)

 32  2c44 32  R5 w32

(12)

12   21  2c6612

(13)

H31  2R6 31  K1w31

(14)

H32  2R5 32  K2 w32

(15)

H33  R111  R2 22  R3 33  K3 w33

(16)

CE

PT

ED

11  c1111  c12 22  c13 33  R1w33

AC

where  ij , and H ij represent the stress and the generalized stress, respectively. The material coefficients cij denote the classical phonon elastic coefficients, Ri represent the phonon-phason coupling parameters, and K i are the phason elastic coefficients. wij (x, x3 , t ) denotes the phason strains, which is not-symmetric and defined as

7

ACCEPTED MANUSCRIPT

w31 (x, x3 , t ) 

w3 (x, x3 , t ) w (x, x3 , t ) w (x, x3 , t ) , w32 (x, x3 , t )  3 , w33 (x, x3 , t )  3 x1 x2 x3

(17)

According to the plate bending assumptions, the phonon displacements u3 is independent on x3 , consistently we can also assume the phason displacements w3 is to be independent

u10 u  1  2  x3 )  c22 ( 20  x3 ) x1 x1 x2 x2

 31   13  c55 ( 1 

u30 w )  R6 3 x1 x1

 32   23  c44 ( 2 

u30 w )  R5 3 x2 x2

 x2



  1  2   u20     x3   x1  x1    x2

u30 w )  K1 3 x1 x1

H 32  R5 ( 2 

u30 w )  K2 3 x2 x2

(20) (21) (22) (23)

u10 u  1  2  x3 )  R2 ( 20  x3 ) x1 x1 x2 x2

(25)

PT

H 33  R1 (

(19)

(24)

ED

H 31  R6 ( 1 

M

 u10

 12   21  c66 

AN US

 22  c12 (

CR IP T

on x3 . As a result, we can know that the phason strain w33 (x, x3 , t )  0 . Then, the generalized Hooke’s law for plate elasticity of orthorhombic QC plate can be written as u u   2  11  c11 ( 10  x3 1 )  c12 ( 20  x3 ) (18) x1 x1 x2 x2

Bending moments M  , normal forces N , shear forces Q , and generalized shear force P are defined as [34]

M   

h /2

(26)

  dx3

(27)

  3 dx3

(28)

P   h /2 H3 dx3

(29)

CE

  x3 dx3

 h /2

N  

h /2

AC

 h /2

Q   

h /2

 h /2

h /2

where  ,  take values 1,2 and   5 / 6 according to the Reissner-Mindlin plate theory.

8

ACCEPTED MANUSCRIPT

We can substitute the equations (18)-(25) into the moment and force (26)-(29), the bending moments M  , normal forces N , shear forces Q and generalized shear forces will be combined with the rotations, displacements, in-plane phonon displacements and phason displacements as follows:

 u10 u  1  2   x3 )  c12 ( 20  x3 )  x3 dx3 c11 ( x1 x1 x2 x2  h/2  h3  1 h3  2  c11  c12 12 x1 12 x2



 u10 u  1  2   x3 )  c22 ( 20  x3 )  x3 dx3 c12 ( x1 x1 x2 x2   h /2  h3  1 h3  2  c12  c22 12 x1 12 x2 h /2

AN US



M 22 (x, t ) 

CR IP T

h/2

M 11 (x, t ) 

   1  2    h3   1  2    u10 u20   c   x  x dx  c            66 3 3 3 66  x1  x1    12  x2 x1   x2  h /2    x2  h /2

(31)

(32)

N11 (x, t ) 

 u10 u u u  1  2   x3 )  c12 ( 20  x3 )  dx3  c11h 10  c12 h 20 c11 ( x1 x1 x2 x2  x1 x2 h/ 2 

(33)

N 22 (x, t ) 

 u10 u u u  1  2   x3 )  c22 ( 20  x3 )  dx3  c12 h 10  c22 h 20 c12 ( x1 x1 x2 x2  x1 x2 h/ 2 

(34)

N12 (x, t ) 

  u10 u20     1  2    u10 u20     ) c66    x3     dx3  c66 h( x1  x1    x2 x1  x2   x2   h /2 

(35)

Q1 (x, t )  

PT

M12 (x, t ) 

(30)

(36)

   u30 w  u  w  )  R5 3 )  dx3   h c44  2  30   R5 3  c44 ( 2  x2 x2  x2  x2  h/ 2   

(37)

h /2

M



h /2

h /2



ED



   u30 w  u  w  )  R6 3 )  dx3   h c55  1  30   R6 3  c55 ( 1  x1 x1  x1  x1   h /2    h /2

CE



Q2 (x, t )  

h /2





   u30 w  u  w  )  K1 3  dx3  h  R6  1  30   K1 3   R6 ( 1  x1 x1  x1  x1   h /2   

(38)

P2 (x, t ) 

   u30 w  u  w  )  K 2 3  dx3  h  R5  2  30   K 2 3   R5 ( 2  x2 x2  x2  x2   h /2   

(39)

AC

P1 (x, t ) 

h /2

h /2



9

ACCEPTED MANUSCRIPT

It is well-known that there has not yet been an agreement, concerning the description of the dynamic behavior of QCs. Agiasofitou and Lazar [9] published the elastodynamic model in which phonons are represented by waves, and phasons by waves damped in time and propagating with finite velocity. In the rest of the paper we will use this model described by the following system of governing equations x j

(x, x3 , t )  X i (x, x3 , t )  

 2ui (x, x3 , t ) t 2

(x, x3 , t )  Gi (x, x3 , t )  

 2 wi wi (x, x3 , t )  D f (x, x3 , t ) 2 t t

H ij x j

(40)

CR IP T

 ij

(41)

where  is the mass density, X i , Gi represent the conventional (phonon) body force

density and a generalized (phason) body force density, respectively [9]. D f  1 is the so-

AN US

w

called phason friction coefficient, where  w is the kinematic coefficient of the phason  2 ui  2 wi , denote the acceleration t 2 t 2 w of the phonon and phason displacements, respectively, and i represents the velocity of t the phason displacements. Moreover, according to the Reissner-Mindlin plate bending theory, the following governing equations will correspond with the elastodynamic model provided by Agiasofitou and Lazar [9]. M  (x, t )  2   Q (x, t )  I M (x, t ) (42) x t 2

ED

M

field of the materials defined by Lubensky et al. [6].

Q (x, t )  2u30  q( x, t )  I Q ( x, t ) x t 2  q (x, t )  I Q

 2u 0 (x, t ) t 2

(44)

CE

x

PT

N (x, t )

(43)

where I M 

 h3 12

, IQ   h

AC

Then, integrating equation (41) in the x3 direction, we can obtain the last governing equation as following: h /2



 h /2



H 3 j x j

h /2

(x, x3 , t )dx3 



 h /2

  2 w3  w3 ( x, x3 , t )  dx3   2 (x, x3 , t )  D f  t  t   h /2  h /2

G3 (x, x3 , t )dx3 



P (x, t ) 2 w w  g (x, t )  I Q 2 3 (x, t )  I D 3 (x, t ) x t t

(45)

10

ACCEPTED MANUSCRIPT

where g (x, t )  h  R1  1 (x, t )  R2  2 (x, t )    G3 (x, x3 , t )dx3 , I D  D f h  h /2 h /2



x1

x2



CR IP T

For the LRBFCM, it is convenient to derive governing equations for primary fields, phonon displacements and rotations, and phason displacements. Substituting equations (30)-(39) into governing equations (42)-(45) we get a system of partial differential equations

u30 w3  2 1  2 1  2 2  2 1 + d  d   d  d  d  I 6 7 1 3 7 8 M x12 x2 2 x1x2 x1 x1 t 2

(46)

d3

u30 w3  2 1  2 2  2 2  2 2 +d6 + d  d   d  d  I 2 4 2 4 5 M x1x2 x12 x2 2 x2 x2 t 2

(47)

d7

 2u30  2u30  2 w3  2 w3  2u30  1  2 +d 4  d7 + d  d  d  q  I 4 8 5 Q x1 x2 x12 x2 2 x12 x2 2 t 2

(48)

c11

 2u10  2u10  2u20 1 1  2u10  c  c  c  q  IQ 2  12 66  66 1 x12 x2 2 x1x2 h h t

(49)

 2u10  2u20  2u20 1 1  2u20 +c66  c  q  IQ 22 2 x1x2 x12 x2 2 h h t 2

(50)

 c12  c66 

 2u30  2u30  1    R2  R5  2  R6  R 5 x1 x2 x12 x2 2

(51)

h3 h3 h3 h3 , d 2  c22 , d3   c12  c66  , d4   hc44 , d5   hR5 , d 6  c66 , 12 12 12 12

PT

where d1  c11

 2 w3  2 w3 1 1  2 w 1 w  K2  G(x, t )  I Q 2 3  I D 3 2 2 x1 x2 h h t h t

ED

 K1

M

 R1  R6 

AN US

d1

d7   hc55 , d8   hR6 , G(x, )  

h /2

 h /2

G3 (x, x3 , )dx3 .

AC

CE

Similarly, one can analyse the QC plate resting on an elastic foundation. If the Winklertype elastic foundation is considered the reaction of the foundation on the plate is proportional to the plate deflection. Then, it is sufficient to modify only the governing equation (43) as Q (x, t )  2u30  q(x, t )  ku30 (x, t )  I Q ( x, t ) x t 2

(52)

where k is the parameter characterising the Winkler foundation. It creates a change only in the partial differential equation (48)

11

ACCEPTED MANUSCRIPT

d7

 2u30  2u30  2 w3  2 w3  2u30  1  2 +d 4  d7 + d  d  d  ku  q  I 4 8 5 30 Q x1 x2 x12 x2 2 x12 x2 2 t 2

(53)

u w  2 1  2 1  2 2  1  2  1 +d 6  d 7 1  d3  d 7 30  d8 3  b1  b12  b6 2 2 x1 x2 x1x2 x1 x1 x1 x2 x2  2  2 1  IM x1 t 2

b6

u w  2 1  2 2  2 2  2   +d 6 +d 2  d 4 2  d 4 30  d5 3  b2 +b21 1 +b61 1 2 x1x2 x1 x2 2 x2 x2 x2 x1 x2

+b61

d7

AN US

d3

 2  2 2  IM x1 t 2

 2u30  2u30  2 w3  2 w3 u u  1  2 +d 4  d7 + d  d  d  b7 1  b4 2  b4 30  b7 30 4 8 5 x1 x2 x12 x2 2 x12 x2 2 x2 x1

w w  2u30 b8 3  b5 3  ku30  q  I Q x1 x2 t 2

 2u10  2u10  2u20 h u10 h u u h  hc  h c  c  c11  c12 20  c66 10   66 12 66 2 2 x1 x2 x1x2 x1 x1 x1 x2 x2 x2

ED

hc11

PT

u  2u10 h  c66 20  q1  I Q x2 x1 t 2

 2u10  2u20  2u20 h u u u h h +hc66  hc  c66 20  c66 10  c12 10 22 x1x2 x12 x2 2 x1 x1 x1 x2 x2 x1

CE

h  c12  c66 

u  2u20 h  c22 20  q2  I Q x2 x2 t 2

AC

h  R1  R6   hK 2 

(54)

(55)

(56)

M

d1

CR IP T

and others are unchanged. It is also interesting to analyse plates with continuously varying plate thickness. The plate thickness occurs in expressions for bending moments (30)-(31), normal forces (33)-(35), shear forces and generalized shear forces (36)-(39). If bending moments, normal forces, shear forces and generalized shear forces with variable plate thickness are substituted into governing equations (42)-(44) we get a system of partial differential equations as

(57)

(58)

 2u30  2u30  2 w3  1   h  R2  R5  2  hR6  hR  hK 5 1 x1 x2 x12 x2 2 x12

 2 w3 h u30 h u h h  R6  R6 1  R5 30  R5 2 2 x2 x1 x1 x1 x2 x2 x2

w w 2w w h h K1 3   K 2 3  G(x, t )  I Q 2 3  I D 3 x1 x1 x2 x2 t t

12

(59)

ACCEPTED MANUSCRIPT

2 2 2 h 2 h h 2 h where b1  c11 h h , b2  c22 h h , b6  c66 h h , b21  c12 , b21  c12 , 4 x2 4 x2 4 x1 4 x2 4 x2

b61  c66

h 2 h h h h h , b  c44 , b7   c55 , b5   R5 , b8   R6 4 x1 4 x2 x1 x2 x1

AN US

CR IP T

The system of partial differential equations (54)-(59) describe the QC plate resting on the Winkler-type elastic foundation and the plates with continuously varying plate thickness. The coupling system consists of 6 equations for 6 unknowns  1 , 2 , u30 , u10 , u20 , w3 . The system of PDEs will be solved numerically by the LRBFCM for given boundary conditions which are introduced as follows. Two types of boundary conditions are discussed in this paper. One is simply-supported and the other is clamped. Two boundary conditions with uniform or transient load on the QC plate demonstrate the possible effect of the coupling parameters on the induced phason and phonon displacements. 2.2 Boundary conditions

2.2.1 Simply-supported boundary conditions

M

In the simply supported boundary conditions, normal traction and bending moment are vanishing. In addition, the phonon and phason displacements u30 , w3 are also vanishing. The simply supported boundary conditions can be expressed as

n1 N11  n2 N12  0 u30  0

ED

n1M11  n2 M12  0

(60)

n1M12  n2 M22  0

(61)

(62)

n1 N12  n2 N22  0

(63)

(64)

w3  0

(65)

AC

CE

PT

n1 , n2 represent the x1 -component and x2 -component of the unit normal vector of the boundary, respectively. The simply supported boundary conditions consist of both Neumann type and Dirichlet type boundary conditions. For the Neumann type boundary conditions, we can substitute the LRBFCM spatial partial differential operators into the derivatives for approximation of the differential term. For the Dirichlet type boundary conditions, we can directly substitute the boundary conditions into the LRBFCM global system. The LRBFCM spatial partial differential operators will be derived in the section 3.

2.2.2 Clamped boundary conditions For the clamped boundary conditions, all primary field variables are vanishing on the boundary edge.

13

ACCEPTED MANUSCRIPT

(66)

2  0

(67)

u10  0

(68)

u20  0

(69)

u30  0

(70)

w3  0

(71)

CR IP T

1  0

All the clamped boundary conditions are Dirichlet type so that we can easily substitute all clamped boundary conditions into the LRBFCM global system. 3

Approach of local radial basis function collocation method

In the LRBFCM, the approximation of a physical field variable  ( x, t ) in the near vicinity of the reference node xgi is supported by N L supported local nodes from the local

NL



AN US

domain  gi as



  x , t    gi ,lj f x  xgi ,lj , xgi ,lj  gi , lj 1

where gi , gj are the global indices, lj are the local indices ,  gi ,lj for (lj  1,2,

(72) , N L ) are

M

the expansion coefficients which will be determined. f is the radial basis function (RBF). (see Fig.2) In this study, we utilize the multiquadrics (MQ) RBF. The trial function associated with the supported local node xk ,m is given as xk  xk ,m

2

(73)

 c2 ,

ED

f ( xk  xk ,m ) 

CE

PT

where c denotes the shape parameter. The shape parameter c determines the shape of RBFs. It is adopted to avoid the singularity when the reference node overlaps its own supported local node. Thus, the accuracy and stability of the LRBFCM solutions are highly relative with the selection of the shape parameter. The value of the optimal shape parameter can be chosen by the criteria elaborating in the followings and also in the proposed studies [43-50].

AC

In this way, we normalize the shape parameter directly as cgi  cconst .  Lgi

(74)

where cgi is the normalized shape parameter of the gi th supported local nodes, cconst . is the given constant shape parameter for each case, and Lgi is the maximum distance between the reference node and its supported local nodes, which corresponds to the scale

14

ACCEPTED MANUSCRIPT

NL

of the problem domain. Moreover, we add an additional constraint



gi ,lj

 0 [45],

lj 1

cgi

into MQRBF.

AN US

proper value of the normalized shape parameter

CR IP T

which states the summation of the expansion coefficients of the supported local nodes equals zero. The given constant shape parameter cconst . is also needed to satisfy this constraint. From equation (74), each reference node has its normalized shape c parameter gi which takes account the local distance Lgi . Thus, we could utilize the

M

Figure 2: The illustration of one reference node with its local supported nodes within the local domain  gi

(75)

PT

Φ gi =f gi α gi

ED

In order to evaluate the expansion coefficients, one should enforce equation (72) at each node from  gi , with getting the system of algebraic equations

CE

where Φgi    xgi ,1 , t  ,  xgi ,2 , t  ,,  xgi , NL , t  column vectors, while

fi

 



 f x

AC

and α gi   gi ,1 , gi ,2 , gi , NL 

is N L  N L matrix

 

 f xgi ,1  xgi ,1   f x x gi ,2 gi ,1 fgi      f xgi , NL  xgi ,1 

T

x f x f

gi ,1

 xgi ,2

gi ,2

 xgi ,2

gi , NL

 xgi ,2

 

x f x



f

f

x

15

gi ,1

 xgi , NL

gi ,2

 xgi , NL

gi , NL

 xgi , NL

    

    NL NL

T

are

ACCEPTED MANUSCRIPT

Since fgi is invertible, the expansion weighting coefficients can be determined by the inverting matrix fgi 1 : α gi =f gi 1Φ gi

(76)

NL

 {  x , t }   gi ,lj { f lj 1

 x  x },

CR IP T

In view of equation (72), we may write

xgi ,lj  gi

gi ,lj

(77)

 {  xgi , t }=  {f gi ( x )}

x  xgi

α gi

AN US

where    = 2 / xi x j   represents any spatial differential operator for the physical field variable  ( x, t ) .Hence the spatial differential value of field variable at each reference node can be expressed in matrix form as

where  {f gi ( x )} is the vector defined as











(78)



 {f gi ( x )}   { f x  x gi ,1 },  { f x  x gi ,2 },,  { f x  x gi , NL } 

M

Substitute α i from equation (76) into equation (78), we can obtain

ED

 {  xgi , t }   {f gi ( x gi )}f gi -1Φgi  s gi Φgi where s gi = {f gi ( xgi )}f gi -1   sgi ,1 sgi ,2

(79) sgi , NL



PT

Then, we can transfer the local system into the global system by the following concept

CE

0, Sgi =  S gi ,1 , S gi ,2, , S gi, N  , S gi , gj  

 sgi ,lj ,

xgj  gi xgj  xgi ,lj  gi

.

(80)

As a result, the global system of the global computational domain  becomes

AC

ζ{Φ}=SΦ

where S  S1 ,S2 ,

(81) , S N  , Φ    x1 , t  ,  x2 , t ,  xN , t  , N is the total number of T

T

nodes. The linear system matrix S is named as the LRBFCM spatial differential operator. All spatial differential operators are derived by the above procedure and utilized to approximate the spatial derivatives in the discretized system of PDEs (54)-(59).

16

ACCEPTED MANUSCRIPT

Moreover, the second partial derivative terms with respect to time t can be discretized by the Houbolt finite difference method [61] Φt t =

2Φt t  5Φt  4Φt t  Φt  2t t 2

(82)

Φt t =

CR IP T

and the first partial derivative terms with respect to time t can be approximated by the backward finite difference method Φt t  Φt t

(83)

AN US

where t is the time-step and the dots over the column-vector Φ indicate differentiations with respect to time t. Then, the discretized system of the PDEs for the Agiasofitou and Lazar model could be rearranged into the LRBFCM system completed with the discretized boundary conditions. The LRBFCM system becomes

t t

 ψ1   B1  ψ  B   2  2 u 30  B 3      u  10  B4  u  B   20   5  w 3  6 N 1 B6  6 N 1

(84)

ED

M

2I M   A12 A13 A14 A15 A16  A11  t 2 I    2I   A 21 A 22  M2 I A 23 A 24 A 25 A 26   t   2IQ   A 31 A 32 A 33  2 I A 34 A 35 A 36   t   2 I Q   A 41 A 42 A 43 A 44  2 I A 45 A 46   t   2IQ   A 51 A 52 A 53 A 54 A 55  2 I A 56   t   2 I I Q D  A 61 A 62 A 63 A 64 A 65 A 66  2 I  I   t t  6 N 6 N

Numerical examples

CE

4

PT

and the details are given in Appendix A. For the localized scheme of the LRBFCM, the global matrix systems are transformed into sparse form which can reduce the risk of illconditioned problem and save computational time.

We consider a square QC plate whose side-length a is taken as 0.254 m. The center of the QC

AC

plate locates at the origin of the Cartesian coordinate system.

21 21 computational nodes with

an orthogonal uniform distribution (Fig.3) are utilized to the approximation of the phonon field variables ( u30 ,  1 , 2 ) and the phason displacements w3 in the three types of plate bending problems. Moreover, 25  25 non-uniform nodes distribution (Fig.4) is used to densify the clamped side of the QC plate resting on Winkler foundation. Both the x1 and

x2 directions of the non-uniform nodes distribution are generated by Chebyshev nodes.

17

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 3: The illustration of the uniform nodal distribution of the mid-plane of the QC

AC

CE

PT

ED

M

plate

Figure 4: The illustration of the non-uniform nodal distribution of the mid-plane of the QC plate

The material properties for the QC plate, which correspond to Al-Ni-Co QC, were measured and collected by Fan [10]. . The relevant elastic constants are c11  c22  23.43 1010 Nm2 ,

c12  5.74 1010 Nm2 ,

18

c44  c55  7.19 1010 Nm2 ,

ACCEPTED MANUSCRIPT

K1  12.2 1010 Nm2 ,

K2  2.4 1010 Nm2 ,

  4180 kg / m3 ,

w 

c66  (c11  c12 ) / 2 ,

1  4.8  1019 m3 s / kg . Df

AN US

4.1 Results for QC plates resting on Winkler foundation

CR IP T

We assume the coupling parameter R1  R2  0 , and R5  R6  R .The coupling parameter R is normalized by the stiffness parameter M  c66 . The phonon and phason displacements are dependent on the ratio of the coupling parameter R to the stiffness parameter M . If the coupling parameter equals zero, the property of the QC plate could be easily transferred to the conventional crystal plate with vanishing phason displacements. In this paper, the ratio of the parameter to the stiffness parameter for the QC plates is chosen as 0.4, R / M  0.4 . The supported local nodes are selected by the nearest 13 nodes from the reference node for all cases. ( N L  13 ).

AC

CE

PT

ED

M

At first, a simply-supported QC plate with constant thickness h  0.012 m resting on Winkler foundation is considered here. On the top of the plate a uniformly distributed transversal load is applied with intensity q  1107 N / m2 . For the LRBFCM model, it is convenient to analyze this static state by eliminating the velocity and acceleration terms. The variation of the phonon and phason displacements in the QC plate along normalized coordinate x1 / a with x2  0 are presented in Fig.5 and Fig.6, respectively. Different values of the foundation parameter k are given in order to observe the reaction of the Winkler foundation. For the case of parameter k equals zero, the QC plate could be reduced to the one without resting on any elastic foundation, and was analyzed early by Sladek et al. [34]. The LRBFCM results for the case k  0 show good agreement with the MLPG results, so that we could further validate our LRBFCM computational model. As can be seen in Fig.5 and Fig.6, both the phonon and phason displacements will decrease with the increasing value of the foundation parameter k.

19

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

Figure 5: Variation of the phonon displacements along normalized coordinate x1 / a for a simply-supported QC plate with constant thickness h  0.012 m resting on a Winkler foundation.

20

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

CE

PT

Figure 6: Variation of the phason displacements along normalized coordinate x1 / a for a simply-supported QC plate with constant thickness h  0.012 m resting on a Winkler foundation.

AC

The variation of the bending moment M11 along normalized coordinate x1 / a for the simply-supported plate is shown in Fig.7. Likewise, we also investigate the same simplysupported QC plate with constant thickness equals 0.012 m to analyse the behavior of the bending moment. The values of the bending moment are also decreasing with the increasing value of the foundation parameter k. For the case of the foundation parameter k equals 1010 Pa/m, the foundation stiffness is quite high so that the QC plate will have less bending moments.

21

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

PT

Figure 7: Variation of the bending moment along normalized coordinate x1 / a for a simply- supported QC plate with constant thickness h  0.012 m resting on a Winkler foundation.

AC

CE

Then, a clamped QC plate resting on Winkler foundation subjected to a uniform distribution load is also analyzed here. The geometric, material parameters and the intensity of the load are all the same as for the former simply-supported QC plate. Numerical results for the phonon and phason displacements given for different values of the foundation parameter k are presented in Fig.8 and Fig.9, respectively. The variation of the bending moment M 11 for the clamped QC plate with various foundation parameter k is presented in Fig.10. The non-uniform nodal distribution is also analyzed in this case, so that we could compare the results between uniform and non-uniform nodal distribution. Numerical results for the phonon and phason displacements and bending moments (Fig.8 to Fig.10) all reveal that good agreement between uniform and non-uniform nodal distribution, which validate our LRBFCM is robust on computational domain utilized by both uniform and non-uniform nodal distributions. In Fig.8 and Fig.9, we could observe

22

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

that the phonon and phason displacements are also decreasing with the increasing value of the foundation parameter k. The slope of the phonon and phason displacement close to the clamped side can be better estimated by using densified non-uniform distribution. The absolute values of the phonon displacements for this clamped plate are lower than the former simply-supported plate due to the stronger restriction of the clamped boundary condition. Moreover, the influence of the foundation parameter k to the phason displacements is more obvious for the simply-supported plate than for the clamped plate in comparison with Fig.6 and Fig.9, respectively. Fig.10 reveals that the stronger foundation stiffness leads to the smaller gradient variation of the bending moment of the plate.

Figure 8: Variation of the phonon displacements along normalized coordinate x1 / a for a clamped QC plate with constant thickness h  0.012 m resting on a Winkler foundation.

23

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

Figure 9: Variation of the phason displacements along normalized coordinate x1 / a for a clamped QC plate with constant thickness h  0.012 m resting on a Winkler foundation.

24

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

PT

Figure 10: Variation of the bending moment along normalized coordinate x1 / a for a clamped QC plate with constant thickness h  0.012 m resting on a Winkler foundation.

AC

CE

In order to verify our selection of the optimal shape parameter in the LRBFCM, we demonstrate one sensitive analysis for the uniform nodal distribution case of the QC plate without resting on any elastic foundation (k equals zero), which case was proposed in MLPG meshless method by Sladek et al. [34]. We estimate the value of the optimal given constant shape parameter cconst . by considering the L2 norm difference value which denotes as the following form, L2 norm difference 

2

N



| ugiLRBF

 ugiMLPG

|

(85)

gi 1

where u giLRBF represents our LRBFCM solutions of the gi th node, u giMLPG represents the proposed MLPG solutions [34]. N is the number of nodes in our computational domain. In this case N equals 441. In Fig.11, we take the results of the phonon displacements for

25

ACCEPTED MANUSCRIPT

PT

ED

M

AN US

CR IP T

two chosen numbers (9 and 13) of supported local nodes N L to analyze here. We could observe that the L2 norm difference value is stable within a wider range of the given constant shape parameter cconst . , when the number of chosen supported local nodes equals 9. However, if we choose proper given constant shape parameter cconst . , we could obtain better agreement between LRBFCM solutions and MLPG solutions when the number of chosen supported local nodes is 13. In this case, the reliable value of given shape parameter cconst . locates within the range of 10 to 20.

CE

Figure 11: L2 norm difference of phonon displacements between the LRBFCM solutions and MLPG solutions with different values of the given shape parameter cconst .

AC

4.2 Results for QC plates with variable thickness For the second numerical examples, a clamped QC plate with variable thickness is 3 x analyzed here. The variable thickness is chosen as a linear function, h( x1 )  h0 (  1 ) , 2 a where h0 equals 0.012 m. That is, this thickness is linearly various for x1 direction and uniform for x2 direction. Moreover, we consider two QC plates with constant thickness.

26

ACCEPTED MANUSCRIPT

CE

PT

ED

M

AN US

CR IP T

One is the same thickness as the first example, 0.012 m; the other is the double thickness, 0.024 m. The material properties for these three QC plates are all the same as in the first example, and all of these QC plates are subjected to the same uniform load q  1107 N / m2 without resting on any elastic foundation. Fig.12 shows the comparison of the phonon displacements along normalized coordinate x1 / a for these three clamped square plate. For the case of the two plates with constant thickness, we could display the robustness of our LRBFCM model by the nice agreements with the MLPG results [34]. The comparison of phason displacements is presented in Fig.13. Both the results of phonon and phason displacements for QC plates with constant thickness are symmetric. Moreover, the absolute values of both the phonon and phason displacements for the QC plate with constant thickness which equals 0.024 m are smaller than displacements of the plate with thickness 0.012 m by the knowledge of the stronger flexural rigidity. For the case of the QC plate with linearly variable thickness, the results of the phonon and phason displacements both are skewed distribution.

AC

Figure 12: Comparison of the phonon displacements along normalized coordinate x1 / a for clamped QC plates with various thickness

27

AN US

CR IP T

ACCEPTED MANUSCRIPT

M

Figure 13: Comparison of the phason displacements along normalized coordinate x1 / a for clamped QC plates with various thickness

AC

CE

PT

ED

Likewise, we also calculated the bending moment M11 for the QC plate with these three thicknesses by equation (30). As can be seen in Fig.14, the results for the QC plate with the both constant thickness are symmetric and almost the same. For the case of QC plate with variable thickness, we could observe the effect to the bending moment due to the gradation of the plate stiffness.

28

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

ED

Figure 14: Comparison of the bending moments along normalized coordinate x1 / a for clamped QC plates with various thickness

PT

4.3 Response of QC plates to transient dynamic loads

AC

CE

The QC plate under transient dynamic loads is analysed too. In this part, we considered a clamped QC square plate with constant thickness, 0.012 m, under a transient dynamic load. The load is expressed as a Heaviside time variation qH (t  0) . We applied the elastodynamic model introduced by Agiasofitou and Lazar to investigate this problem. The time-step t was taken as 1.25 105 sec for our LRBFCM numerical analysis. Numerical results for the phonon and phason displacements at the center of the QC plate stat are normalized by the corresponding static LRBFCM value u30 (0,0)  1.6776 103 m , which have been validated with published MLPG results [34]. The LRBFCM numerical results are presented for various values of the phason friction. When the phason friction coefficient equals zero, the employed elastodynamic model becomes identical with the model provided by Bak [3]. The temporal variations of the normalized phonon and phason displacements at the plate center can be seen in Fig.15 and Fig.16, respectively. We could observe strong influence of the phason friction coefficient on the phason displacements. The bending moments are normalized by the corresponding static

29

ACCEPTED MANUSCRIPT

PT

ED

M

AN US

CR IP T

LRBFCM value, M11stat (0,0)  1.44 105 Nm , which has been validated with published MLPG results [34]. The normalized bending moments at the center of the clamped QC plate with various phason frictions, M11 (0,0) / M11stat (0,0) are shown in Fig.17. With respect to the Agiasofitou and Lazar model, Fig.15 and Fig.17 reveal that the phason friction coefficients have little influence on the phonon displacements and bending moments.

AC

CE

Figure 15: Comparison of the temporal variations of the normalized phonon displacements for the clamped QC plate under a transient dynamic load with various phason friction coefficients by Agiasofitou and Lazar model

30

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Figure 16: Comparison of the temporal variations of the normalized phason displacements for the clamped QC plate under a transient dynamic load with various phason friction coefficients by Agiasofitou and Lazar model

31

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

PT

5. Conclusions

ED

Figure 17: Comparison of the temporal variations of the normalized bending moment for the clamped QC plate under a transient dynamic load with various phason friction coefficients by Agiasofitou and Lazar model

AC

CE

The researches of the QCs are becoming more popular in recent years. The fields in QC problems are coupled and it is needed to have a reliable, advanced numerical method to analyse general boundary value problems. However, only a few papers are devoted to the QC problems by numerical methods. Moreover, the commercial computer codes which are based on the FEM or BEM cannot be utilized to analyse QC problems due to lacking of adequate modifications to these multi-field problems. As the first attempt, the LRBFCM is developed for bending problems of QC plates in this paper. The ReissnerMindlin plate bending theory reduces the actual 3D thick plate problem to 2D governing equations in mid-plane of plate. The LRBFCM scheme is the strong meshless formulation which is applied for the plate bending problem. This paper succeeds to use the LRBFCM to solve some bending problems in QC plates and partially verifies the LRBFCM results by nice agreements with MLPG results. The numerical results reveal the influence of the elastic foundation parameters and variable thickness of the plate on the phonon and phason displacements and bending

32

ACCEPTED MANUSCRIPT

CR IP T

moments in QC plates under a static uniform load. Furthermore, the transient dynamic analysis of the QC plates under Heaviside transient dynamic load is also presented in this study. The phason friction coefficient shows strong influence on the phason displacements but little influence on phonon displacements and bending moments. The peak value of the phason displacements is reduced with the increasing phason friction. Due to the ease implementation and flexibility of adjusting the governing equations, this paper shows that LRBFCM is suitable to analyse static and transient dynamic plate bending problems in QC plates. In the future, we may extend more complex bending problems by the LRBFCM.

AN US

Acknowledgement The authors gratefully acknowledge the supports by the Ministry of Science and Technology (MOST) of Taiwan under the grant number: NSC 103-2923-E002-004-MY2, as well as Slovak Research and Development Agency under the contract No. APVV-14-0440. It is greatly appreciated. References:

AC

CE

PT

ED

M

[1] D. Shechtman, I. Blech, D. Gratias, J. W. Cahn, Metallic phase with long-range orientational order and no translational symmetry. Phys. Rev. Lett. 53 (1984)1951-1953. [2] L. D. Landau, E. M. Lifshitz, Theoretical Physics V: Statistical Physics, third ed., Pergamon Press, New York, 1980. [3] P. Bak, Phenomenological theory of icosahedral incommensurate (quasiperiodic) order in Mn-Al alloys. Phys. Rev. Lett. 54(14) (1985) 1517-1519. [4] D. H. Ding, W. G. Yang, C. Z. Hu, R. H. Wang, Generalized elasticity theory of quasicrystals. Phys. Rev. B 48(10) (1993) 7003-7010. [5] D. H. Ding, W. G. Yang, C. Z. Hu, R. H. Wang, Linear elasticity of quasicrystals and defects in quasicrystals. Mater. Sci. Forum. 150 (1994) 345-354. [6] T. C. Lubensky, S. Ramaswamy, J. Joner, Hydrodynamics of icosahedral quasicrystals. Phys. Rev. B 32 (1985) 7444-7452. [7] S. B. Rochal, V. L. Lorman, Minimal model of the phonon-phason dynamics in icosahedral quasicrystals and its application to the problem of internal friction in the iAlPdMn alloy. Phys. Rev. B 66 (2002) 144204. [8] T. Y. Fan, X. F. Wang, W. Li, A. Y. Zhu, Elasto-hydrodynamics of quasicrystals. Philos. Mag. 89(6) (2009) 501-512. [9] E. Agiasofitou, M. Lazar, The elastodynamic model of wave-telegraph type for quasicrystals, I. J. Solids and Struct. 51(5) (2014) 923–929. [10] T. Y. Fan, Mathematical Theory of Elasticity of Quasicrystals and its Applications. Springer, Beijing, 2011. [11] T. Y. Fan, Mathematical theory and methods of mechanics of quasicrystalline materials. Eng. 5 (2013) 407-448. [12] C. L. Li, L.Y. Liu, The physical property tensors of one-dimensional quasicrystals. Chin. Phys. 13 (2004) 924-931.

33

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

[13] W. Q. Chen, Y. L. Ma, H. J. Ding, On three-dimensional elastic problems of onedimensional hexagonal quasicrystal bodies. Mech. Res. Commun. 31(6) (2004) 633-641. [14] X. F. Wang, E. Pan, Analytical solutions for some defect problems in 1D hexagonal and 2D octagonal quasicrystals. Pranama J. Phys. 70 (2008) 911-933. [15] Y. Gao, B. X. Xu, B. S. Zhao, T. C. Chang, New general solutions of plane elasticity of one-dimensional quasicrystals. Phys. Stat. Sol. b 245(1) (2008) 20–27. [16] Y. Gao, A. Ricoeur, The refined theory of one-dimensional quasi-crystals in thick plate structures. J. Appl. Mech. ASME. 78 (2001) 031021. [17] Y. Gao, The exact theory of one-dimensional quasicrystal deep beams. Acta Mech. 212 (2010) 283-292. [18] X. F. Wang, F. T. Fan, Study on the dynamics of the double cantilever-beam specimen of decagonal Al-Ni-Co quasicrystals. Appl. Math. Comput. 211 (2009) 336-346. [19] G. R. Liu, Meshfree Methods: Moving Beyond the Finite Element Method. CRC Press. (2002) [20] J. Sladek, P. Stanak, Z. D. Han, V. Sladek, S. N. Atluri, Applications of the MLPG method in engineering & sciences: A review. Comput. Model. Eng. Sci. 92 (2013) 423475. [21] W. H. Chen, M. H. Lee, A novel meshless analysis procedure for three-dimensional structural problems with complicated geometry. Comput. Model. Eng. Sci. 93 (2013) 149-166. [22] Q. G. Liu, B. Šarler, Non-singular method of fundamental solutions for twodimensional isotropic elasticity problems. Comput. Model. Eng. Sci. 91 (2013) 235-268. [23] L. Dong, S. N. Atluri, SGBEM Voronoi cells (SVCs), with embedded arbitraryshaped inclusions, voids, and/or cracks, for micromechanical modeling of heterogeneous materials. Comput. Mater. Contin. 33 (2013) 111-154. [24] Z. D. Han, S. N. Atluri, On the (Meshless Local Petrov-Galerkin) MLPG-Eshelby method in computational finite deformation solid mechanics - Part II. Comput. Model. Eng. Sci. 97 (2014) 199-237. [25] L. Dong, A. Alotaibi, S. A. Mohiuddine, S.N. Atluri, Computational methods in engineering: A variety of primal & mixed methods, with global & local Interpolations, for well-posed or ill-posed BCs. Comput. Model. Eng. Sci. 99 (2014) 1-85. [26] A. Mazzia, G. Pini, F. Sartoretto, MLPG refinement techniques for 2D and 3D diffusion problems. Comput. Model. Eng. Sci. 102 (2014) 475-497. [27] M. Uddin, H. Ali, A. Ali, Kernel-based local meshless method for solving multidimensional wave equations in irregular domain. Comput. Model. Eng. Sci. 107 (2015) 463-479. [28] L. H. Kuo, M. H. Gu, D. L. Young, C. Y. Lin, Domain type kernel-based meshless methods for solving wave equations. Comput. Mater. Contin. 33 (2013) 213-28. [29] J. Sladek, V. Sladek, H. A. Mang, Meshless formulations for simply supported and clamped plate problems. I. J. Numer. Methods Eng. 55(3) (2002) 359-375.

34

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

[30] J. Sladek, V. Sladek, H. A. Mang, Meshless local boundary integral equation method for simply supported and clamped plates resting on elastic foundation. Comput. Methods Appl. Mech. Eng. 191 (2002) 5943-5959. [31] L. Sator, V. Sladek, J. Sladek, Coupling effects in elastic analysis of FGM composite plates by mesh-free methods. Compos. Struct. 115 (2014) 100-110. [32] V. Sladek, J. Sladek, L. Sator, Physical decomposition of thin plate bending problems and their solution by mesh-free methods. Eng. Anal. Bound. Elem. 37 (2013) 348-365. [33] L. Sator, V. Sladek, J. Sladek, D.L. Young, Elastodynamics of FGM plates by meshfree method, Compos. Struct. 140 (2016) 309-322. [34] J. Sladek, V. Sladek, E. Pan, Bending analyses of 1D orthorhombic QC plates. I. J. Solids Struct. 50 (2013) 3975-3983. [35] J. Sladek, V. Sladek, C. Zhang, M. Wünsche, Modelling of orthorhombic quasicrystal shallow shells. Euro. J. Mech.-A/Solids, 49 (2015) 518-530. [36] H. C. Yaslan, Equations of anisotropic elastodynamics in 3D quasicrystals as a symmetric hyperbolic system: deriving the time-dependent fundamental solutions. Appl. Math. Modell. 17 (2013) 8409-8418. [37] R. L. Hardy, Multiquadric equations of topography and other irregular surfaces. J. Geophys. Res. 76(8) (1971) 1905-1915. [38] R. Franke, Scattered data interpolation: Tests of some methods. Math. Comput. 38 (1982) 181-200. [39] E. J. Kansa, Multiquadrics-a scattered data approximation scheme with applications to computational fluid-dynamics-I surface approximations and partial derivative estimates. Comput. Math. Applic. 19(8-9) (1990) 127-145. [40] E. J. Kansa, Multiquadrics-a scattered data approximation scheme with applications to computational fluid-dynamics-II solutions to parabolic, hyperbolic and elliptic partial differential equations. Comput. Math. Applic. 19(8-9) (1990) 147-161. [41] M. A. Golberg, C. S. Chen, S. R. Karur, Improved multiquadric approximation for partial differential equations. Eng. Anal. Bound. Elem. 18(1) (1996) 9-17. [42] B. Šarler, A radial basis function collocation approach in computational fluid dynamics. Comput. Model. Eng. Sci.7 (2005) 185-193. [43] L. H. Shen, K. H. Tseng, D. L. Young, Evaluation of multi-order derivatives by local radial basis function differential quadrature method. J. Mech. 29(1) (2013) 67-78. [44] Y. L. Chan, C. T. Wu, L. H. Shen, D. L. Young, The interpolation technique for scattered data by local radial basis function differential quadrature method. Int. J. Comput. Math. 10(2) (2013) 1341011. [45] K. H. Tseng, Evaluation of Multi-Order Derivatives by Local Radial Basis Function Differential Quadrature Method. Master Thesis, National Taiwan Univ., Taipei, Taiwan. [46] S. Rippa, An algorithm for selecting a good value for the parameter c in radial basis function interpolation. Adv. Comput. Math. 11(2) (1999), 193-210.

35

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

[47] M. Uddin, On the selection of a good value of shape parameter in solving timedependent partial differential equations using RBF approximation method. Appl. Math. Model. 38(1) (2014), 135-144. [48] M. Gherlone, L. Iurlaro, M. Di Sciuva, A novel algorithm for shape parameter selection in radial basis functions collocation method. Compos. Struct. 94(2) (2012), 453461. [49] L. Iurlaro, M. Gherlone, M. Di Sciuva, Energy based approach for shape parameter selection in radial basis functions collocation method. Compos. Struct. 107 (2014), 70-78. [50] B. Mavrič, B. Šarler, Local radial basis function collocation method for linear thermoelasticity in two dimensions. Int. J. Heat Fluid Flow. 25(6) (2015) 1488-1510 [51] B. Šarler, R. Vertnik, Meshfree explicit local radial basis function collocation method for diffusion problems. Comput. Math. Applic. 51(8) (2006) 1269-1282. [52] G. Kosec, B. Sarler, Local RBF collocation method for Darcy flow. Comput. Model. Eng. Sci. 25(3) (2008) 197-208. [53] C.K. Chou, C. P. Sun, D. L. Young, J. Sladek, V. Sladek, Extrapolated local radial basis function collocation method for shallow water problems. Eng. Anal. Bound. Elem. 50 (2015) 275-290. [54] Y.C. Hon, B. Šarler, D. F. Yun, Local radial basis function collocation method for solving thermo-driven fluid-flow problems with free surface. Eng. Anal. Bound. Elem. 57 (2015) 2-8. [55] U. Hanoglu, B. Šarler. Simulation of hot shape rolling of steel in continuous rolling mill by local radial basis function collocation method. Comput. Model. Eng. Sci. 109(5) (2015) 447-479. [56] R. Vertnik, B. Šarler, Local radial basis function collocation method along with explicit time stepping for hyperbolic partial differential equations. Appl. Numer. Math. 67 (2013) 136-151. [57] S. Krahulec, J. Sladek, V. Sladek, Y. C. Hon, Meshless analyses for time-fractional heat diffusion in functionally graded materials. Eng. Anal. Bound. Elem. 62 (2016) 57-64. [58] H. H.-H. Lu, D.L. Young, J. Sladek, V. Sladek, Three-dimensional analysis for functionally graded piezoelectric semiconductors, J. Intell. Mater. Syst. Struct. (2016) DOI:10.1177/1045389x16672566. [59] R. D. Mindlin, Influence of rotary inertia and shear on flexural motions of isotropic, elastic plates. J. Appl. Mech. ASME 18 (1951) 31-38. [60] J. N. Reddy, Mechanics of Laminated Composite Plates: Theory and Analysis. CRC Press, Boca Raton, 1997. [61] J. C. Houbolt, A recurrence matric solution for the dynamic response of aircraft. J. Aeronaut. Sci. 17 (1950) 371-376.

36

ACCEPTED MANUSCRIPT

Appendix A For the LRBFCM, we use the LRBFCM spatial differential operators S i and S ij to

CR IP T

approximate the partial spatial derivatives where i,j represent the Cartesian coordinate indices. The application of the LRBFCM to the governing equations yields the following system of algebraic equations.

where I denotes the identity matrix.

A21   d3S12  b21S1  b61S2 N  N A41  [0]N  N A51  [0]N  N

ED

A31   d7S1  b7 I N  N

M

A11   d1S11  d6S22  d7 I  b1S1  b6S2 N  N

AN US

2I M   A12 A13 A14 A15 A16  A11  t 2 I    2I   A 21 A 22  M2 I A 23 A 24 A 25 A 26   t   2IQ   A 31 A 32 A 33  2 I A 34 A 35 A 36   t   2IQ   A 41 A 42 A 43 A 44  2 I A 45 A 46   t   2 I Q   A 51 A 52 A 53 A 54 A 55  2 I A 56   t   2 I I Q  A 61 A 62 A 63 A 64 A 65 A 66  2 I  D I   t t  6 N 6 N

CE

PT

  h A61   h( R1  R6 )S1  R6 I  x1   N N

A12   d3S12  b12S2  b6S1 N  N

A22   d6S11  d2S22  d4I  b2S2  b61S1 N  N

A32   d4S2  b4 I N  N A42  [0]N  N A52  [0]N  N   h A62   h( R2  R5 )S 2  R5I   x 2   N N

A13   d7S1 N  N

A14  0N  N

A23   d4S2 N  N

A24  0N  N

A33   d7S11  d4S22  b4S2  b7S1  kI N  N

A34  0N  N

A43  [0]N  N

A44  [0]N  N

A53  [0]N  N

A54  [0]N  N

  h h A63   hR6S11  hR5S 22  R6S1  R5S 2   x  x 1 2   N N

A64  [0]N  N

AC

t t

 ψ1   B1  ψ  B   2  2 u 30  B 3      u10  B4  u  B  20    5  w 3  6 N 1 B 6  6 N 1

37

ACCEPTED MANUSCRIPT

A16   d8S1 N  N

A25  0N  N

A26   d5S2 N  N

A35  0N  N

A36   d8S11  d5S22  b8S2  b7S1 N  N

A45  [0]N  N

A46  [0]N  N

A55  [0]N  N

A56  [0]N  N

A65  0N  N

  h h A66   hK1S11  hK 2S 22  K1S1  K 2S 2  x1 x2   N N

where 0 is noted as null matrix. ψ1   1,1 , 1,2 ,..., 1, N  N 1

CR IP T

A15  0N  N

AN US

u10  u10,1 , u10,2 ,..., u10, N  N 1

ψ 2   2,1 , 2,2 ,..., 2, N  N 1

u20  u20,1 , u20,2 ,..., u20, N  N 1

u30  u30,1 , u30,2 ,..., u30, N  N 1

w 3   w3,1 , w3,2 ,..., w3, N  N 1

4 1  5  B1    2 I M ψ1t  2 I M ψ1t-Δt  2 I M ψ1t-2Δt  t t  t  N 1

M

4 1  5  B2    2 I M ψ 2t  2 I M ψ 2t-Δt  2 I M ψ 2t-2Δt   t  t  t   N 1

ED

4 1  5  B3    2 I Q u 30 t  2 I Q u 30 t-Δt  2 I Q u 30 t-2Δt  q   t  t  t   N 1

4 1  5  B4    2 I Q u10 t  2 I Q u10 t-Δt  2 I Q u10 t-2Δt  t t  t  N 1

PT

4 1  5  B5    2 I Q u 20t  2 I Q u 20 t-Δt  2 I Q u 20 t-2Δt   t  t  t   N 1

CE

1 4 1  5  B6    2 I Q w 3t  I D w 3t  2 I Q w 3t-Δt  2 I Q w 3t-2Δt  t t t  t  N 1

AC

The discretized equations should be solved numerically.

38