Bending analysis of FG plates using a general third-order plate theory with modified couple stress effect and MLPG method

Engineering Analysis with Boundary Elements 94 (2018) 159–171

Contents lists available at ScienceDirect

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

Bending analysis of FG plates using a general third-order plate theory with modified couple stress effect and MLPG method V.S. Khorasani a,∗, M. Bayat b a b

Department of Mechanical Engineering, Faculty of Engineering, Kharazmi University, Mofatteh Avenue, P.O. Box 15719-14911, Tehran, Iran Department of Mechanical Engineering, Islamic Azad University, Yadegar-e-Imam Khomeini (RAH) Shahre Rey Branch, P.O. Box 18155-144, Tehran, Iran

a r t i c l e

i n f o

Keywords: Meshless local Petrov–Galerkin method General third-order plate theory Modified couple stress effect Functionally graded plate Moving least-square approximation

a b s t r a c t Meshless local Petrov–Galerkin analysis of functionally graded plates based on a general third-order shear deformation plate theory with a modified couple stress effect is presented. Governing equations of problem are a fourth-order partial differential equations system which derived in terms of eleven generalized displacement variable, by applying the principle of virtual displacements. The moving least-squares approach is used for approximation of unknown variables and the Gauss weight function is employed as test function for obtaining local weak form. The Gauss–Legendre quadrature method is utilized for numerical integration of weak equations. Static bending results of a simply-supported plate is obtained for various power law index and length scale parameter, and is compared to analytical solutions that shows high accuracy in results.

1. Introduction In the mechanical problems, either time-dependent or timeindependent, it is recommended to reduce three-dimensional governing equations of problems to a two-dimensional formulations by a proper method. One of the best approaches for reducing three-dimensional elasticity is representing a dimension of problem through power series and multiplying this series to the functions that describe roles of other two dimensions. For analyzing plate structures, it is preferred to define the direction of thickness by proper order of power series to account for the kinematic of deformation and derive constitutive relations. This is due to fact that thickness of plates is quite small compared to their inplane dimensions. Types of two-dimensional plate theories consist of two categories, displacement-based and stress-based. These two categories are similar in expansion of the fields by increasing powers of the thickness coordinate, and different in their strain/stress compatibility conditions. Displacement-based theories are strain/stress compatible, therefore these theories have preferred in literature and in this research too. To find the optimal order of power series for expanding thickness coordinate in formulation, it should be considered that what degree of variation of strains and stresses is recommended to have better results. By expanding displacement field up to three order for in-plane displacements and two order for out-plane displacements in most of cases, the transverse shear strains have quadratic variation through the plate thickness. This approach is the third-order shear deformation theory (TSDT) and if all terms in power series remain in expansion, the theory is called

∗

the general third-order plate theory (GTPT). In addition to having proper variation for the transverse shear strains, TSDT or GTPT do not need to any shear correction factor unlike the first-order shear deformation theory (FSDT). In conventional continuum mechanics, the effects of micro- and nano-scale interactions are not considered. For analyzing size-dependent behavior of micro-scale plate structures, the modified couple stress theory can be implemented in deriving governing equations. This theory has only one length scale parameter for representing microstructural effect. This is an advantage of the modified couple stress approach in comparison to the classical couple stress theories, due to complexity of determining two length scale parameter that each of them related together and to size-dependent effect. A meshless method is a method used to establish system algebraic equations for the whole problem domain without the use of a predefined mesh for the domain discretization [1]. Meshless methods use a set of nodes scattered within the problem domain as well as sets of nodes scattered on the boundaries of the domain to represent (not discretize) the problem domain and its boundaries. These scattered nodes do not form a mesh, so it does not need to any a priori information on the relationship between scattered nodes for the interpolation or approximation of the unknown functions of field variables [2]. There are two basic forms of meshless methods: the global form and the local form. The meshless methods based on the global form can be applied easily for solving partial differential equations (PDE) and integral equations. But these methods are not advantageous for solv-

Corresponding author. E-mail address: [email protected] (V.S. Khorasani).

https://doi.org/10.1016/j.enganabound.2018.06.015 Received 23 February 2018; Received in revised form 23 May 2018; Accepted 14 June 2018 0955-7997/© 2018 Elsevier Ltd. All rights reserved.

V.S. Khorasani, M. Bayat

Engineering Analysis with Boundary Elements 94 (2018) 159–171

ing all categories of PDEs. To overcome these disadvantages, the meshless methods based on local approach have been introduced. These approaches consist of two classes [3]:

moving least-squares (MLS) approach for approximation of unknown variables.

• Local meshless methods based on the variational weak form • Local meshless methods based on the strong form.

2. Displacements and strains of GTPT If in-plane displacement and out-of-plane displacements extend up to third power and second power of thickness direction, respectively, the displacement field of GTPT can be formulated as [25,26]:

In the first class, it is necessary to use a numerical integration procedure for solving the local weak equations. This procedure can be meshless or use background cells for discretizing integration domain. In the numerous meshless methods that have applied for solving partial differential equations until today, the meshless local Petrov– Galerkin method (MLPG) is one of the few methods that is truly meshless, because this approach do not need to meshing either for the interpolation of unknown variables or for the numerical integration of weak equations. In other words, no domain or boundary element is required for discretizing the domain of problem. The MLPG method is based on local weak form that means for every node in the domain, weak form equations are satisfied in local sub-domains around nodes. In view of the fact that no global integration is involved, the MLPG method is addressed as a truly meshless method [4]. This method has been successfully applied for solving a wide range of problems in engineering [5].

𝑢1 (𝑥, 𝑦, 𝑧, 𝑡) = 𝑢(𝑥, 𝑦, 𝑡) + 𝑧𝜃𝑥 (𝑥, 𝑦, 𝑡) + 𝑧2 𝜑𝑥 (𝑥, 𝑦, 𝑡) + 𝑧3 𝜓𝑥 (𝑥, 𝑦, 𝑡) 𝑢2 (𝑥, 𝑦, 𝑧, 𝑡) = 𝑣(𝑥, 𝑦, 𝑡) + 𝑧𝜃𝑦 (𝑥, 𝑦, 𝑡) + 𝑧2 𝜑𝑦 (𝑥, 𝑦, 𝑡) + 𝑧3 𝜓𝑦 (𝑥, 𝑦, 𝑡) 𝑢3 (𝑥, 𝑦, 𝑧, 𝑡) = 𝑤(𝑥, 𝑦, 𝑡) + 𝑧𝜃𝑧 (𝑥, 𝑦, 𝑡) + 𝑧2 𝜑𝑧 (𝑥, 𝑦, 𝑡)

(1)

where u, v, w, 𝜃 x , 𝜃 y , 𝜃 z , 𝜙x , 𝜙y , 𝜙z , 𝜓 x , and 𝜓 y are unknown generalized displacements. In this expansion, normality and straight condition of the transverse normal lines under the Kirchhoff assumptions are released because of cubic variation of in-plane displacements. Also, inextensibility of transverse normal lines is released due to quadratic variation of outof-plane displacements. Here linearized strains are considered and their relations to the generalized displacements take the form: 𝜕 𝜃𝑥 𝜕 𝜙𝑥 ⎧ 𝜀 ⎫ ⎧ 𝜕𝑢 ⎫ ⎧ ⎫ ⎧ ⎫ ⎧ 𝜕 𝜓𝑥 ⎫ 𝜕𝑥 𝜕𝑥 𝜕𝑥 ⎪ 𝑥𝑥 ⎪ ⎪ 𝜕𝑥 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ 𝜕 𝜃 𝜕 𝜙 𝑦 𝑦 𝜕𝑣 𝜕 𝜓𝑦 ⎪ 𝜀𝑦𝑦 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ 𝜕𝑦 𝜕𝑦 𝜕𝑦 ⎪ ⎪ ⎪ 𝜕𝑦 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ 0 0 ⎪ 𝜀𝑧𝑧 ⎪ ⎪ 𝜃𝑧 ⎪ ⎪ 2𝜙𝑧 ⎪ ⎪ ⎪ ⎪ ⎪ 2 3 ⎨ ⎬ = ⎨ 𝜕𝑢 𝜕𝑣 ⎬ + 𝑧⎨ 𝜕 𝜃𝑥 𝜕 𝜃𝑦 ⎬ + 𝑧 ⎨ 𝜕 𝜙𝑥 𝜕 𝜙𝑦 ⎬ + 𝑧 ⎨ 𝜕 𝜓 𝜕 𝜓𝑦 ⎬ + 𝑥 2 𝜀 + + ⎪ 𝑥𝑦 ⎪ ⎪ 𝜕𝑦 𝜕𝑥 ⎪ ⎪ 𝜕𝑦 ⎪ 𝜕𝑦 ⎪ 𝜕𝑦 + 𝜕𝑥 ⎪ 𝜕𝑥 ⎪ 𝜕𝑥 ⎪ 𝜕 𝜃𝑧 ⎪ 𝜕 𝜙𝑧 ⎪ ⎪2𝜀 ⎪ ⎪𝜃 + 𝜕𝑤 ⎪ ⎪ ⎪ ⎪ ⎪ 0 ⎪ 𝑥𝑧 ⎪ ⎪ 𝑥 𝜕𝑥 ⎪ ⎪2𝜙𝑥 + 𝜕𝑥 ⎪ ⎪3𝜓𝑥 + 𝜕𝑥 ⎪ ⎪ ⎪ ⎪2𝜀𝑦𝑧 ⎪ ⎪𝜃𝑦 + 𝜕𝑤 ⎪ ⎪2𝜙 + 𝜕 𝜃𝑧 ⎪ ⎪3𝜓 + 𝜕 𝜙𝑧 ⎪ ⎪ ⎪ 0 𝜕𝑦 ⎩ ⎭ ⎩ ⎭ ⎩ 𝑦 ⎩ 𝑦 ⎩ ⎭ 𝜕𝑦 ⎭ 𝜕𝑦 ⎭

1.1. Present study After the MLPG method introduced by Atluri and Zhu [6] in 1998, and discussed in detail by Atluri and his co-authors [7–9], several researches is published until today which employed the MLPG approach in higher-order plate theories. Qian et al. [10] analyzed deformation of a homogenous and isotropic thick plate with a higher-order shear and normal deformable plate theory (HOSNDPT) by using MLPG method that displacements according to thickness derived by assuming the Legendre polynomial as basis function. They also computed static deformations, and free and forced vibrations of a thick functionally graded plate by applying HOSNDPT and MLPG method [11]. Qian and Batra [12,13] studied transient thermoelastic deformations of a thick FG plate by HOSNDPT and MLPG approach. Comparison of the MLPG method and the finite element method (FEM) presented in this research. They also designed a FG plate for optimal natural frequencies by employing MLPG method and HOSNDPT [14]. Xiao et al. [15] analyzed thick FG plates by using HOSNDPT and MLPG approach. They expanded displacement field by the Legendre polynomial basis and utilized the Radial Basis Functions (RBFs) for interpolation of variables. Gilhooley et al. [16] and Xiao et al. [17] attempted this concept for analyzing thick FG plates and thick composite laminates, respectively. In addition to MLPG method, several researches investigated solutions of higher-order plate theories by employing various methods and compared results with MLPG solutions. Ferreira et al. [18,19] solved TSDT governing equations of FG plates by using the collocation RBF method and compared their results with Qian et al. [11] solutions. Natural frequencies of a FG plate are achieved by employing the TSDT and the collocation RBF approach by Ferreira et al. [20] and compared with Qian et al. [11] results. Sheikholeslami and Saidi [21] obtained solutions of vibration of a FG plate with HOSNDPT and analytical approach, and compared their results with MLPG method [11]. Reddy and Kim [22] proposed the GTPT with the modified couple stress effect and derived governing equations by employing the principle of virtual displacements and the fundamental lemma of variational calculus. For first time, analytical solutions for bending, buckling and vibration of FG plates obtained by Kim and Reddy [23] by using Navier solution technique. They also computed solutions for bending of FG plates with von Kármán nonlinearity by implementing finite element method (FEM) with conforming element which has four degree of freedom per node [24]. In this paper, the GTPT governing equations with modified couple stress effect of FG plates are solved by using the MLPG method with the

(2) 3. Governing equations of GTPT For deriving governing equations of GTPT based on the modified couple stress theory of FG plates, it is needed to represent modified couple stress model and constitutive relations of considered FG plate. 3.1. Modified couple stress model Yang et al. [27] proposed a modification of the classical couple stress theory (see [28]), which established that the couple stress tensor is symmetric and the symmetric curvature tensor is the only proper conjugate strain criterion that enters in equation of the total strain energy of the body. In the modified couple stress theory, strain energy density function depends only on the strain and the symmetric part of the curvature tensor, therefore only one length scale parameter is involved. This fact and inclusion of a symmetric couple stress tensor are two main advantages of the modified couple stress model over the classical couple stress theory. Virtual strain energy 𝛿 using modified couple stress model is defined as [29]: 𝛿 =

∫𝑉

(𝛿 𝜺 ∶ 𝝈 + 𝛿 𝝌 ∶ 𝐦) 𝑑𝑣

(3)

where summation on repeated indices is implied: here 𝜎 ij are Cartesian components of the symmetric part of the stress tensor, 𝜀ij denotes the strain components, mij are the components of the deviatoric part of the symmetric couple stress tensor, and 𝜒 ij donates the symmetric curvature tensor components and can be written as [30]: 𝝌=

] 1[ ∇𝝎 + (∇𝝎)𝑇 , 2

or 𝜒𝑖𝑗 = 160

1 2

(

) 𝜕 𝜔𝑖 𝜕 𝜔𝑗 + , 𝜕 𝑥𝑗 𝜕 𝑥𝑖

𝝎=

1 ∇×𝐮 2

𝑖, 𝑗 = 1, 2, 3

(4)

(5)

V.S. Khorasani, M. Bayat

Engineering Analysis with Boundary Elements 94 (2018) 159–171

Components of the rotation vector, 𝝎 (𝜔x , 𝜔y , and 𝜔z ) in terms of generalized displacements are: 𝜕𝜃

𝜕𝜙

⎧𝜔 ⎫ ⎧𝜔 ⎫ ⎧ 𝜕𝑤 − 𝜃𝑦 ⎫ ⎧ 𝜕𝑦𝑧 − 2𝜙𝑦 ⎫ ⎧ 𝜕𝑦𝑧 − 3𝜓𝑦 ⎫ ⎪ 𝑥⎪ ⎪ 1 ⎪ ⎪ 𝜕𝑦 𝜕𝑤 ⎪ ⎪ ⎪ 𝜕 𝜃𝑧 ⎪ 𝜕𝜙 ⎪ 2 2⎨𝜔𝑦 ⎬ = 2⎨𝜔2 ⎬ = ⎨𝜃𝑥 − 𝜕𝑥 ⎬ + 𝑧⎨2𝜙𝑥 − 𝜕𝑥 ⎬ + 𝑧 ⎨3𝜓𝑥 − 𝜕𝑥𝑧 ⎬ ⎪ 𝜔𝑧 ⎪ ⎪𝜔3 ⎪ ⎪ 𝜕𝑣 − 𝜕𝑢 ⎪ ⎪ 𝜕 𝜃𝑦 𝜕 𝜃𝑥 ⎪ ⎪ 𝜕 𝜙𝑦 𝜕 𝜙𝑥 ⎪ ⎩ ⎭ ⎩ ⎭ ⎩ 𝜕𝑥 𝜕𝑦 ⎭ ⎩ 𝜕𝑥 − 𝜕𝑦 ⎭ ⎩ 𝜕𝑥 − 𝜕𝑦 ⎭ ⎧ 0 ⎪ 0 + 𝑧3 ⎨ ⎪ 𝜕 𝜓𝑦 − ⎩ 𝜕𝑥

⎫ ⎪ ⎬ 𝜕 𝜓𝑥 ⎪ 𝜕𝑦 ⎭

(6)

thus curvature tensor components take the form: ( ) ⎧ ⎫ 𝜕 𝜕𝑤 − 𝜃𝑦 ⎧𝜒 ⎫ ⎧𝜒 ⎫ ⎪ ⎪ 𝜕𝑥 ( 𝜕𝑦 ) 𝑥𝑥 11 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ 𝜕 𝜕𝑤 𝜃 − ⎪ 𝜒𝑦𝑦 ⎪ ⎪ 𝜒22 ⎪ ⎪ ⎪ 𝑥 𝜕𝑦 𝜕𝑥 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ 𝜕 𝜃𝑦 𝜕 𝜃𝑥 − ⎪ 𝜒𝑧𝑧 ⎪ ⎪ 𝜒33 ⎪ ⎪ ( ⎪ 𝜕𝑥) 𝜕𝑦( ) 2⎨ ⎬ = 2⎨ ⎬ = ⎨ 𝜕 𝜕𝑤 𝜕 𝜕𝑤 ⎬ 𝜒 2 𝜒 − 𝜃 − + 𝜃 ⎪ 𝑥𝑦 ⎪ ⎪ 12 ⎪ ⎪ 𝜕𝑦 𝜕𝑦 𝑦 𝑥 𝜕𝑥 𝜕𝑥 ⎪ ⎪𝜒 ⎪ ⎪2𝜒 ⎪ ⎪ 𝜕 ( 𝜕𝑣 𝜕𝑢 ) ⎪ 𝜕 𝜃𝑧 ⎪ 𝑥𝑧 ⎪ ⎪ 13 ⎪ ⎪ 𝜕𝑥 𝜕𝑥 − 𝜕𝑦 + 𝜕𝑦 − 2𝜙𝑦 ⎪ ⎪ 𝜒𝑦𝑧 ⎪ ⎪2𝜒23 ⎪ ⎪ 𝜕 ( 𝜕𝑣 𝜕𝑢 ) ⎪ 𝜕𝜃 ⎩ ⎭ ⎩ ⎭ ⎪ − 𝜕𝑦 + 2𝜙𝑥 − 𝜕𝑥𝑧 ⎪ 𝜕𝑦 𝜕𝑥 ⎩ ⎭ ( ) ⎫ ⎧ 𝜕 𝜃𝑧 𝜕 − 2𝜙𝑦 ⎪ ⎪ 𝜕𝑥 ( 𝜕𝑦 ) ⎪ ⎪ 𝜕 𝜃𝑧 𝜕 𝜙 − 2 𝑥 ⎪ ⎪ 𝜕𝑦( 𝜕𝑥 ) ⎪ ⎪ 𝜕𝜙 𝜕𝜙 2 𝜕𝑥𝑦 − 𝜕𝑦𝑥 ⎪ ⎪ ) ( )⎬ + 𝑧⎨ 𝜕 ( 𝜕 𝜃 𝜕 𝜃𝑧 𝜕 𝑧 ⎪ 𝜕𝑦 𝜕𝑦 − 2𝜙𝑦 + 𝜕𝑥 2𝜙𝑥 − 𝜕𝑥 ⎪ ( )⎪ ⎪ 𝜕 ( 𝜕 𝜃𝑦 𝜕 𝜃𝑥 ) 𝜕 𝜙𝑧 ⎪ 𝜕𝑥 𝜕𝑥 − 𝜕𝑦 + 2 𝜕𝑦 − 3𝜓𝑦 ⎪ ( )⎪ ⎪ 𝜕 ( 𝜕 𝜃𝑦 𝜕 𝜃𝑥 ) 𝜕𝜙 ⎪ 𝜕𝑦 𝜕𝑥 − 𝜕𝑦 + 2 3𝜓𝑥 − 𝜕𝑥𝑧 ⎪ ⎩ ⎭ ( ) ⎧ ⎫ 𝜕 𝜙𝑧 𝜕 − 3𝜓𝑦 ⎧ ⎪ ⎪ 0 𝜕𝑥 ( 𝜕𝑦 ) ⎪ ⎪ ⎪ 𝜕 𝜙𝑧 𝜕 3𝜓𝑥 − 𝜕𝑥 0 ⎪ ⎪ ⎪ 𝜕𝑦( ) ⎪ ⎪ ⎪ 𝜕 𝜓𝑦 𝜕 𝜓𝑥 0 3 − ⎪ ⎪ ⎪ 𝜕𝑥 𝜕𝑦( ) ) + 𝑧3 ⎨ + 𝑧2 ⎨ 𝜕 ( 𝜕 𝜙 𝜕 𝜙𝑧 ⎬ 0 𝜕 𝑧 ⎪ 𝜕𝑦 𝜕𝑦 − 3𝜓𝑦 + 𝜕𝑥 3𝜓𝑥 − 𝜕𝑥 ⎪ ⎪ ( 𝜕𝜓 ( 𝜕𝜙 ) 𝜕 𝑦 ⎪ ⎪ ⎪ − 𝜕 𝜙 𝜕 𝑦 − 𝜕𝑦𝑥 ⎪ ⎪ ⎪ 𝜕𝑥 ( 𝜕𝜕𝑥 𝜕𝑥 ( 𝜕𝑥 𝜓𝑦 ) 𝜕 ⎪ ⎪ ⎪ 𝜕𝜙 𝜕 𝜕 𝜙𝑦 ⎩ 𝜕𝑦 𝜕𝑥 − − 𝜕𝑦𝑥 ⎪ ⎪ 𝜕𝑦 𝜕𝑥 ⎩ ⎭

Fig. 1. The non-dimensional resultants of modulus.

⎫ ⎪ ⎪ ⎪ ⎪ ⎬ )⎪ 𝜕 𝜓𝑥 ⎪ 𝜕𝑦 )⎪ 𝜕 𝜓𝑥 ⎪ 𝜕𝑦 ⎭ (7)

metal, respectively and n denotes the volume fraction exponent, called power-law index [22]. The function f(z) has a key role for defining material properties of FGMs and consists of one or more parameters for representing proper volume fraction of two constituents [35–37]. In this research, the FG plate is assumed isotropic and isothermal, with variation of two constituents (ceramic and metal) through thickness. Thus, the linear constitutive equations can be written as [23,38]: ⎧𝜎𝑥𝑥 ⎫ ⎡1 − 𝜈 ⎪ ⎪ ⎢ 𝜈 𝜎 ⎪ 𝑦𝑦 ⎪ ⎢ ⎪ 𝜎𝑧𝑧 ⎪ ⎢ 𝜈 = Λ ⎨𝜎 ⎬ ⎢ 0 ⎪ 𝑥𝑦 ⎪ ⎢ ⎪𝜎𝑥𝑧 ⎪ ⎢ 0 ⎪𝜎 ⎪ ⎢ ⎣ 0 ⎩ 𝑦𝑧 ⎭

𝜈 1−𝜈 𝜈 0 0 0

⎧ 𝜀𝑥𝑥 ⎫ ⎪ ⎪ ⎪ 𝜀𝑦𝑦 ⎪ ⎪ 𝜀𝑧𝑧 ⎪ ×⎨ ⎬ ⎪2𝜀𝑥𝑦 ⎪ ⎪2𝜀𝑥𝑧 ⎪ ⎪ 2𝜀 ⎪ ⎩ 𝑦𝑧 ⎭

3.2. FG plate constitutive equations Functionally graded materials (FGMs) are a class of materials that have a predetermined variation of material properties from point surface to another (see [31,32]). FGMs are nonhomogeneous but often isotropic or isothermal. A typical FGM is one in which two materials are mixed to create a composition that provides desired functionality. For example, thermal-barrier structures are made of a mixture of ceramics and metals. The gradation in properties of the material reduces thermal stresses, residual stresses, and stress concentrations [22]. The material properties of FGMs vary continuously and smoothly in the thickness direction z and are functions of volume fractions of constituent materials [33]. Each mechanical quantity related to the first constituent, ceramic, is denoted by the subscript c, whereas the subscript m is used to specify the properties of the metal [34]. A most common material property of the FGM through the thickness direction is assumed to be represented by a power-law: ( ) [ ] 1 𝑧 𝑛 𝑃 (𝑧) = 𝑃𝑐 − 𝑃𝑚 𝑓 (𝑧) + 𝑃𝑚 , 𝑓 (𝑧) = + (8) 2 ℎ where Pc and Pm are the values of a typical material property, such as the modulus, density, and conductivity, of the ceramic material and

𝜈 𝜈 1−𝜈 0 0 0

0 0 0 1 (1 − 2𝜈) 2 0 0

0 0 0 0 1 1 − 2𝜈) ( 2 0

0 ⎤ ⎥ 0 ⎥ 0 ⎥ ⎥ 0 ⎥ 0 ⎥ ⎥ 1 1 − 2 𝜈) ( ⎦ 2

(9)

where Λ = 𝐸∕[(1 + 𝜈)(1 − 2𝜈)]. Poisson’s ratio, 𝜈, is assumed as constant but Young’s modulus, E, varies through thickness direction by power law as: ( ) 1 𝑧 𝑛 𝐸(𝑧) = (𝐸𝑡 − 𝐸𝑏 )𝑓 (𝑧) + 𝐸𝑏 , 𝑓 (𝑧) = + (10) 2 ℎ where the subscripts t and b indicate top (𝑧 = ℎ2 ) and bottom (𝑧 = − ℎ2 ) surfaces and h is thickness of plate. The resultants of Young’s modulus are: 𝐸𝑘 =

ℎ 2 ∫− ℎ 2

𝑧(𝑘) 𝐸 (𝑧)𝑑 𝑧

(11)

According to Eqs. (10) and (11), Ek decreases for larger power law index, n, due to the assumption that material on top surface is stiffer than one on bottom surface. The variation of non-dimensional Ek in terms of power law index is showed in Fig. 1.

161

V.S. Khorasani, M. Bayat

Engineering Analysis with Boundary Elements 94 (2018) 159–171

The generalized forces and couples can be written as: (𝑖) ⎧𝑀𝑥𝑥 ⎫ ℎ ⎪ (𝑖) ⎪ 2 ⎨𝑀𝑦𝑦 ⎬ = ∫ ℎ − ⎪ (𝑖) ⎪ 2 ⎩𝑀𝑧𝑧 ⎭

⎧𝜎𝑥𝑥 ⎫ ⎡ (𝑘) 3+𝑖 𝐴11 ∑ ⎪ ⎪ 𝑖 ⎢ (𝑘) ⎨ 𝜎𝑦𝑦 ⎬(𝑧) 𝑑𝑧 = ⎢𝐴12 ⎪𝜎 ⎪ 𝑘=𝑖 ⎢ (𝑘) ⎣𝐴 ⎩ 𝑧𝑧 ⎭

⎧𝑀 ( 𝑖 ) ⎫ ℎ ⎪ 𝑥𝑦 2 (𝑖) ⎪ ⎨𝑀𝑥𝑧 ⎬ = ∫ ℎ − ⎪𝑀 ( 𝑖 ) ⎪ 2 ⎩ 𝑦𝑧 ⎭

⎧𝜎 ⎫ ⎡ (𝑘) 3+𝑖 𝐵11 ∑ ⎪ 𝑥𝑦 ⎪ 𝑖 ⎢ ⎨𝜎𝑥𝑧 ⎬(𝑧) 𝑑𝑧 = ⎢ 0 ⎪𝜎𝑦𝑧 ⎪ 𝑘=𝑖 ⎢ 0 ⎣ ⎩ ⎭

⎧ ( 𝑖 ) ⎫ ℎ ⎪ 𝑥𝑥 2 (𝑖) ⎪ ⎨ 𝑦𝑦 ⎬ = ∫ ℎ − ⎪ (𝑖) ⎪ 2 ⎩ 𝑧𝑧 ⎭

⎧𝑚 ⎫ ⎡ (𝑘) 3+𝑖 11 ∑ ⎪ 𝑥𝑥 ⎪ 𝑖 ⎢ ⎨ 𝑚𝑦𝑦 ⎬(𝑧) 𝑑𝑧 = ⎢ 0 ⎪ 𝑚𝑧𝑧 ⎪ 𝑘=𝑖 ⎢ 0 ⎣ ⎩ ⎭

0 (11𝑘) 0

(𝑘−𝑖) ⎫ 0 ⎤⎧𝜒𝑥𝑥 ⎥⎪ (𝑘−𝑖) ⎪ 0 ⎥⎨𝜒𝑦𝑦 ⎬ (𝑘−𝑖) ⎪ (11𝑘) ⎥⎦⎪ ⎩𝜒𝑧𝑧 ⎭

(14)

⎧ ( 𝑖 ) ⎫ ℎ ⎪ (𝑥𝑦 2 𝑖) ⎪ =  ⎨ 𝑥𝑧 ⎬ ∫ ℎ − ⎪ (𝑖) ⎪ 2 ⎩ 𝑦𝑧 ⎭

⎧𝑚 ⎫ ⎡ (𝑘) 3+𝑖 11 ∑ ⎪ 𝑥𝑦 ⎪ 𝑖 ⎢ ⎨𝑚𝑥𝑧 ⎬(𝑧) 𝑑𝑧 = ⎢ 0 ⎪𝑚𝑦𝑧 ⎪ 𝑘=𝑖 ⎢ 0 ⎣ ⎩ ⎭

0 (11𝑘) 0

(𝑘−𝑖) ⎫ 0 ⎤⎧𝜒𝑥𝑦 ⎥⎪ (𝑘−𝑖) ⎪ 0 ⎥⎨𝜒𝑥𝑧 ⎬ (𝑘−𝑖) ⎪ (11𝑘) ⎥⎦⎪ ⎩𝜒𝑦𝑧 ⎭

(15)

12

𝐴(12𝑘) 𝐴(11𝑘) 𝐴(12𝑘) 0 (𝑘) 𝐵11 0

𝑘−𝑖) ⎫ 𝐴(12𝑘) ⎤⎧𝜀(𝑥𝑥 (𝑘) ⎥⎪ (𝑘−𝑖) ⎪ 𝐴12 ⎥⎨𝜀𝑦𝑦 ⎬ ⎥⎪ ⎪ 𝐴(11𝑘) ⎦⎩𝜀(𝑧𝑧𝑘−𝑖) ⎭

(12)

(𝑘−𝑖) ⎫ 0 ⎤⎧𝛾𝑥𝑦 ⎥⎪ (𝑘−𝑖) ⎪ 0 ⎥⎨𝛾𝑥𝑧 ⎬ (𝑘) ⎥⎪ (𝑘−𝑖) ⎪ 𝐵11 ⎦⎩𝛾𝑦𝑧 ⎭

(13)

where couple stresses (mij ) and plate stiffnesses (𝐴11 , 𝐴12 , 𝐵11 , and 11 ) are following: 𝑚𝑖𝑗 = 2𝜇𝑙2 𝜒𝑖𝑗 ℎ

𝐴(11𝑘) =

2 (𝑘) 1−𝜈 𝑧 𝐸 (𝑧)𝑑 𝑧, (1 + 𝜈)(1 − 2𝜈) ∫− ℎ 2 ℎ 2 ∫− ℎ 2

𝐴(12𝑘) =

𝜈 (1 + 𝜈)(1 − 2𝜈)

(𝑘) 𝐵11 =

2 (𝑘) 1 𝑧 𝐸 (𝑧)𝑑 𝑧 , 2(1 + 𝜈) ∫− ℎ

𝑧(𝑘) 𝐸 (𝑧)𝑑 𝑧

ℎ

ℎ

(11𝑘) =

2

2 (𝑘) 𝑙2 𝑧 𝐸 (𝑧)𝑑 𝑧 (1 + 𝜈) ∫− ℎ

(16)

2

3.3. The partial differential equations system of motion By implementing the principle of virtual displacements, equation of motion of the GTPT can be derived that based on linear kinematic relation, a modified couple stress model, and variation of two constituents through thickness. The dynamic version of the principle of virtual displacements is generalized to the Hamilton’s principle (17). Equation of motion of the GTPT in terms of the generalized forces derived by Reddy and Kim [22]. Hamilton’s principle is: 𝑇

∫0

(𝛿 − 𝛿 − 𝛿 )𝑑𝑡 = 0

(17)

where 𝛿 is the virtual kinematic energy, 𝛿 is the virtual strain energy, and 𝛿 is the virtual work done by external forces, and defined as [22]: (

𝛿 =

ℎ 2 ∫Ω ∫− ℎ 2

𝜌

𝛿 =

ℎ 2 ∫Ω ∫− ℎ 2

(𝛿 𝜺 ∶ 𝝈 + 𝛿 𝝌 ∶ 𝐦)𝑑 𝑧𝑑 𝑥𝑑 𝑦

𝜕 𝑢1 𝜕𝛿𝑢1 𝜕 𝑢2 𝜕𝛿𝑢2 𝜕 𝑢3 𝜕𝛿𝑢3 + + 𝜕𝑡 𝜕𝑡 𝜕𝑡 𝜕𝑡 𝜕𝑡 𝜕𝑡

) (18)

(19)

ℎ ⎡ ( ) ) 2 ( ⎢ 𝛿 = − 𝑞𝑥𝑡 𝛿𝑢1 + 𝑞𝑦𝑡 𝛿𝑢2 + 𝑞𝑧𝑡 𝛿𝑢3 𝑑 𝑥𝑑 𝑦 𝑓𝑥 𝛿𝑢1 + 𝑓𝑦 𝛿𝑢2 + 𝑓𝑧 𝛿𝑢3 𝑑 𝑧𝑑 𝑥𝑑 𝑦 + ℎ ∫Ω+ ⎢ ∫Ω ∫− ⎣ 2 ] ( ) ( ) + 𝑡𝑥 𝛿𝑢1 + 𝑡𝑥 𝛿𝑢2 + 𝑡𝑥 𝛿𝑢3 𝑑𝑆 𝑞𝑥𝑏 𝛿𝑢1 + 𝑞𝑦𝑏 𝛿𝑢2 + 𝑞𝑧𝑏 𝛿𝑢3 𝑑 𝑥𝑑 𝑦 + ∫Ω− ∫𝑆

(20) where 𝜌 is mass density, 𝑞𝜉𝑡 and 𝑞𝜉𝑏 are distributed forces on top and bottom surfaces, f𝜉 is the body forces, and 𝑡𝜉 are surface forces on lateral surfaces of plate; subscript 𝜉 can taking the values of x, y, or z. For deriving equation of motion of GTPT, the length scale parameter, l, has important mathematical role, because if l = 0, modified couple stress model suspends and equation of motion will be a PDE system of order two. But if l ≠ 0, modified couple stress installs and PDE system of motion will be of order four. In this research, solutions of both of these two cases for various values of power law index are presented. Substituting kinematic relation (2) and constitutive equations (9) into generalized forces (12)–(15), and then substituting obtained relations into equation of motion in terms of generalized forces [22], we acquire equation of motion in terms of generalized displacements. Linear static PDE system of motion of the GTPT in terms of generalized displacements are following: 162

V.S. Khorasani, M. Bayat

Engineering Analysis with Boundary Elements 94 (2018) 159–171

( ) 2 2 2 𝜕 2 𝜃𝑥 𝜕 4 𝜃𝑥 𝜕2 𝑢 1 𝜕4 𝑢 1 𝜕4 𝑢 1 𝜕4 𝑣 1 𝜕4 𝑣 1 (0) 𝜕 𝑢 (0) 𝜕 𝑣 (1) 𝜕 𝜃𝑥 + 𝐵11 − (0) − (0) + 𝐴(0) + 𝐵11 + (0) + (0) + 𝐴(1) + 𝐵11 − (1) 11 11 12 11 11 11 11 2 2 2 2 4 3 3 2 2 4 4 𝜕 𝑥𝜕 𝑦 4 4 4 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕 𝑥𝜕 𝑦 𝜕𝑥 𝜕𝑦 𝜕 𝑥2 𝜕 𝑦2 2 4 4 ( ) ( ) 4 2 2 4 𝜕 𝜃𝑦 𝜕 𝜃𝑦 𝜕 𝜃𝑥 𝜕 𝜃𝑧 𝜕 𝜙𝑥 𝜕 𝜙𝑥 1 (2) 𝜕 𝜙𝑥 1 1 1 1 (1) 𝜕 𝜃𝑦 (2) − (1) + 𝐴(1) + 𝐵11 + (1) + (1) + 𝐴(0) + 𝐴(2) + 𝐵11 + (0) − 11 12 12 𝜕𝑥 11 𝜕 𝑥2 4 11 𝜕 𝑦4 𝜕 𝑥𝜕 𝑦 4 11 𝜕 𝑥3 𝜕𝑦 4 11 𝜕 𝑥𝜕 𝑦3 2 11 𝜕 𝑦2 4 𝜕 𝑥2 𝜕 𝑦2 2 4 4 ) 𝜕 𝜙𝑦 𝜕 𝜙𝑦 𝜕 𝜙𝑦 𝜕 4 𝜙𝑥 ( (2) 𝜕 𝜙𝑧 𝜕 2 𝜓𝑥 ( (3) 3 (1) ) 𝜕 2 𝜓𝑥 1 1 1 1 (2) − (2) + 𝐴12 + 𝐵11 − (0) + (2) + (2) + 2𝐴(1) + 𝐴(3) + 𝐵11 + 11 12 𝜕𝑥 11 𝜕 𝑥2 4 11 𝜕 𝑦4 2 11 𝜕 𝑥𝜕 𝑦 4 11 𝜕 𝑥3 𝜕𝑦 4 11 𝜕 𝑥𝜕 𝑦3 2 𝜕 𝑦2 (0) 2𝜓 4𝜓 4𝜓 ( ) 4 4 𝜕 𝜕 𝜕 𝜕 𝜓𝑥 𝜕 𝜓𝑥 𝑦 𝑦 𝑦 1 1 3 1 1 1 𝜕𝑐𝑧 (3) − (3) − (3) + 𝐴(3) + 𝐵11 − (1) + (3) + (3) + 𝐹𝑥(0) + =0 12 4 11 𝜕 𝑥2 𝜕 𝑦2 4 11 𝜕 𝑦4 2 11 𝜕 𝑥𝜕 𝑦 4 11 𝜕 𝑥3 𝜕𝑦 4 11 𝜕 𝑥𝜕 𝑦3 2 𝜕𝑦

𝛿𝑢 ∶ 𝐴(0) 11

(

( ) 2 2 𝜕 4 𝜃𝑥 1 𝜕4 𝑢 1 𝜕4 𝑢 𝜕 2 𝑣 1 (0) 𝜕 4 𝑣 1 (0) 𝜕 4 𝑣 1 𝜕2 𝑢 (0) 𝜕 𝑣 (1) 𝜕 𝜃𝑥 + (0) + (0) + 𝐵11 + 𝐴(0) − 11 −  + 𝐴(1) + 𝐵11 + (1) 11 𝜕 𝑦2 12 𝜕 𝑥𝜕 𝑦 4 11 𝜕 𝑥3 𝜕 𝑦 4 11 𝜕 𝑥𝜕 𝑦3 4 𝜕 𝑥𝜕 𝑦 4 11 𝜕 𝑥3 𝜕 𝑦 𝜕 𝑥2 𝜕 𝑥4 4 11 𝜕 𝑥2 𝜕 𝑦2 2 ) 𝜕2 𝜙 𝜕 2 𝜃𝑦 1 (1) 𝜕 4 𝜃𝑦 1 (1) 𝜕 4 𝜃𝑦 𝜕 4 𝜃𝑥 𝜕 𝜃𝑧 ( (2) 𝜕 4 𝜙𝑥 𝜕 4 𝜙𝑥 1 1 1 1 (1) 𝜕 𝜃𝑦 (2) 𝑥 + (1) + 𝐵11 + 𝐴(1) − 11 − 11 + 𝐴(0) + 𝐴12 + 𝐵11 − (0) + (2) + (2) 11 11 12 11 11 11 3 2 2 4 2 2 3 4 4 4 𝜕𝑦 2 𝜕 𝑥𝜕 𝑦 4 𝜕 𝑥𝜕 𝑦 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕 𝑥 𝜕𝑦 4 𝜕 𝑥𝜕 𝑦3 2 4 4 ( ) 𝜕 2 𝜙𝑦 ( ) 2 4 𝜕 𝜙𝑦 1 (2) 𝜕 𝜙𝑦 1 (2) 𝜕 𝜙𝑦 𝜕 𝜙𝑧 𝜕 4 𝜓𝑥 𝜕 𝜓𝑥 1 (3) 𝜕 𝜓𝑥 1 3 1 (2) (3) + 𝐵11 + (0) + 𝐴(2) − 11 − 11 + 2𝐴(1) + 𝐴(3) + 𝐵11 − (1) + 11 + (3) 11 11 12 12 11 11 2 2 4 2 2 3 2 4 4 𝜕𝑦 2 𝜕 𝑥𝜕 𝑦 4 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕 𝑥 𝜕𝑦 4 𝜕 𝑥𝜕 𝑦3 2 2 4 4 (0) ( ) 𝜕 𝜓𝑦 𝜕 𝜓𝑦 1 (3) 𝜕 𝜓𝑦 1 (3) 𝜕 𝜓𝑦 3 1 𝜕𝑐𝑧 (3) + 𝐵11 + (1) + 𝐴(3) − 11 − 11 + 𝐹𝑦(0) + =0 11 𝜕 𝑦2 2 11 𝜕 𝑥2 4 4 2 𝜕𝑥 𝜕 𝑥4 𝜕 𝑥2 𝜕 𝑦2

𝛿𝑣 ∶

(0) + 𝐵11 𝐴(0) 12

(0) 𝛿𝑤 ∶ 𝐵11

(21)

)

(22)

3 3 2 2 2 𝜕 4 𝜃𝑧 𝜕2 𝑤 1 𝜕 4 𝑤 1 (0) 𝜕 4 𝑤 1 (0) 𝜕 𝜃𝑥 1 (0) 𝜕 𝜃𝑦 1 (0) 𝜕 𝑤 (0) 𝜕 𝜃𝑥 (0) 𝜕 𝜃𝑦 (1) 𝜕 𝜃𝑧 (1) 𝜕 𝜃𝑧 + 𝐵11 − (0) −  + 𝐵 +  + 𝐵 +  + 𝐵11 + 𝐵11 − (1) 11 𝜕 𝑥 11 𝜕 𝑦 4 11 𝜕 𝑥4 4 11 𝜕 𝑦4 4 11 𝜕 𝑥3 4 11 𝜕 𝑦3 4 11 𝜕 𝑥4 𝜕 𝑥2 𝜕 𝑦2 𝜕 𝑥2 𝜕 𝑦2

2 2 𝜕 3 𝜙𝑦 𝜕 4 𝜃𝑧 𝜕 3 𝜙𝑥 𝜕 4 𝜙𝑧 1 (2) 𝜕 4 𝜙𝑧 1 1 1 1 (1) 𝜕 𝜙𝑥 (1) 𝜕 𝜙𝑦 (2) 𝜕 𝜙𝑧 (2) 𝜕 𝜙𝑧 (2) 𝜕 𝜓𝑥 − (1) + 2𝐵11 + (1) + 2𝐵11 + (1) + 𝐵11 + 𝐵11 − (2) − 11 + 3𝐵11 11 11 11 11 4 3 3 2 2 4 𝜕𝑥 2 𝜕𝑦 2 4 4 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕 𝑥4 𝜕 𝑦4

+

(0) (0) 3 𝜕 3 𝜓𝑦 3 (2) 𝜕 𝜓𝑥 3 1 𝜕𝑐𝑦 1 𝜕𝑐𝑥 (2) 𝜕 𝜓𝑦 (0) 11 + 3𝐵11 + (2) + 𝐹 + + =0 𝑧 4 𝜕𝑦 4 11 𝜕 𝑦3 2 𝜕𝑥 2 𝜕𝑦 𝜕 𝑥3

(23)

( ) 2 2 𝜕2 𝑢 1 𝜕4 𝑢 1 𝜕4 𝑢 1 𝜕4 𝑣 1 𝜕4 𝑣 1 𝜕3 𝑤 (1) 𝜕 𝑢 (1) 𝜕 𝑣 (0) 𝜕𝑤 (0) + 𝐵11 − (1) − (1) + 𝐴(1) + 𝐵11 + (1) + (1) − 𝐵11 − (0) − 𝐵11 𝜃𝑥 11 11 12 11 11 11 2 2 2 2 4 3 3 4 4 𝜕 𝑥𝜕 𝑦 4 4 𝜕𝑥 4 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕 𝑥 𝜕𝑦 𝜕 𝑥𝜕 𝑦 𝜕 𝑥3 ( ) 𝜕2 𝜃 ( ) 𝜕2 𝜃 ) 𝜕 2 𝜃𝑦 𝜕 4 𝜃𝑦 𝜕 4 𝜃𝑦 𝜕 4 𝜃𝑥 𝜕 4 𝜃𝑥 ( (2) 1 1 1 1 1 1 1 (2) (2) 𝑥 𝑥 + 𝐴(2) + (0) + 𝐵11 + (0) − (2) − (2) + 𝐴12 + 𝐵11 − (0) + (2) + (2) 11 11 11 11 11 11 11 11 2 2 2 2 4 3 4 2 4 4 2 𝜕 𝑥𝜕 𝑦 4 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕 𝑥 𝜕𝑦 4 𝜕 𝑥𝜕 𝑦3 ( ) ( ) 2 ( ) 2 ( ) 2 3 4 4 1 (1) 𝜕 𝜃𝑧 3 (1) 𝜕 𝜙𝑥 1 (3) 𝜕 𝜙𝑥 1 (3) 𝜕 𝜙𝑥 3 (1) 𝜕 𝜙𝑦 (1) 𝜕 𝜃𝑧 (1) (3) 1 (1) 𝜕 𝜙𝑥 (3) (3) (3) + 𝐴(1) − 𝐵 −  − 2 𝐵 𝜙 + +  + +  −  −  + + 𝐵 −  𝐴 𝐵 𝐴 𝑥 12 11 11 11 2 11 11 12 11 𝜕𝑥 4 11 𝜕 𝑥3 2 11 𝜕 𝑦2 4 11 𝜕 𝑥2 𝜕 𝑦2 4 11 𝜕 𝑦4 2 11 𝜕 𝑥𝜕 𝑦 𝜕 𝑥2 4𝜙 4𝜙 ( ) ( ) ( ) 3 2 2 4 𝜕 𝜕 𝜕 𝜙𝑧 𝜕 𝜓𝑥 𝜕 𝜓𝑥 1 (4) 𝜕 𝜓𝑥 𝑦 𝑦 1 1 1 3 (2) 𝜕 𝜙𝑧 (2) (4) + (3) + (3) + 2𝐴(2) − 𝐵11 − (2) − 3𝐵11 𝜓𝑥 + 𝐴(4) + (2) + 𝐵11 + 3(2) − 11 12 11 11 4 11 𝜕 𝑥3 𝜕𝑦 4 11 𝜕 𝑥𝜕 𝑦3 𝜕𝑥 4 11 𝜕 𝑥3 4 11 𝜕 𝑥2 4 𝜕 𝑦2 𝜕 𝑥2 𝜕 𝑦2 2 4 4 (1) ( ) 4 𝜕 𝜓𝑦 𝜕 𝜓𝑦 𝜕 𝜓𝑦 𝜕 𝜓𝑥 1 1 1 1 1 𝜕𝑐𝑧 (4) − (4) + 𝐴(4) + 𝐵11 − 3(2) + (4) + (4) + 𝐹𝑥(1) + 𝑐𝑦(0) + =0 (24) 12 11 𝜕 𝑥𝜕 𝑦 4 11 𝜕 𝑦4 4 11 𝜕 𝑥3 𝜕𝑦 4 11 𝜕 𝑥𝜕 𝑦3 2 2 𝜕𝑦

𝛿𝜃𝑥 ∶ 𝐴(1) 11

(

)

2 1 𝜕4 𝑢 1 𝜕4 𝑢 𝜕 2 𝑣 1 (1) 𝜕 4 𝑣 1 (1) 𝜕 4 𝑣 1 𝜕3 𝑤 𝜕2 𝑢 (1) 𝜕 𝑣 (0) 𝜕𝑤 + (1) + (1) + 𝐵11 + 𝐴(1) − 11 −  − 𝐵11 − (0) 11 𝜕 𝑦2 𝜕 𝑥𝜕 𝑦 4 11 𝜕 𝑥3 𝜕𝑦 4 11 𝜕 𝑥𝜕 𝑦3 4 𝜕𝑦 4 11 𝜕 𝑦3 𝜕 𝑥2 𝜕 𝑥4 4 11 𝜕 𝑥2 𝜕 𝑦2 ( ) 𝜕2 𝜃 ( ) 𝜕 2 𝜃𝑦 ( ) 𝜕 2 𝜃𝑦 𝜕 4 𝜃𝑦 𝜕 4 𝜃𝑥 𝜕 4 𝜃𝑥 1 1 1 1 1 1 (2) (0) (2) 𝑥 + 𝐴(2) + 𝐵11 − (0) + (2) + (2) − 𝐵11 𝜃𝑦 + 𝐵11 + (0) + 𝐴(2) + (0) − (2) 12 11 11 11 11 11 11 11 3 3 2 2 2 𝜕 𝑥𝜕 𝑦 4 2 4 4 𝜕 𝑥 𝜕𝑦 4 𝜕 𝑥𝜕 𝑦 𝜕𝑥 𝜕𝑦 𝜕 𝑥4 4 ( ) ( ) ( ) 𝜕 2 𝜙𝑦 3 2 4 4 𝜕 𝜃𝑦 𝜕 𝜃𝑧 𝜕 𝜙𝑥 𝜕 𝜙𝑥 1 (3) 𝜕 𝜙𝑥 1 1 3 1 3 (1) 𝜕 𝜃𝑧 (3) (1) (3) − (2) + 𝐴(1) − 𝐵11 − (1) + 𝐴(3) + 𝐵11 − (1) + 11 + (3) − 2𝐵11 𝜙𝑦 + 𝐵11 + (1) 11 12 11 12 11 11 11 4 𝜕𝑦 4 2 𝜕 𝑥𝜕 𝑦 4 2 𝜕 𝑥2 𝜕 𝑦2 𝜕 𝑦3 𝜕 𝑥3 𝜕𝑦 4 𝜕 𝑥𝜕 𝑦3 𝜕 𝑥2 4𝜙 4𝜙 ( ) 𝜕 2 𝜙𝑦 ( ) ( ) 3𝜙 2𝜓 4𝜓 𝜕 𝜕 𝜙 𝜕 𝜕 𝜕 𝜕 𝑦 𝑦 1 1 1 1 1 (2) (4) 𝑧 𝑧 𝑥 𝑥 + 𝐴(3) + (1) − (3) − (3) + 2𝐴(2) − 𝐵11 − (2) + 𝐴(4) + 𝐵11 − 3(2) + (4) 11 12 12 11 𝜕 𝑥𝜕 𝑦 2 11 𝜕 𝑦2 4 11 𝜕 𝑥4 4 11 𝜕 𝑥2 𝜕 𝑦2 𝜕𝑦 4 11 𝜕 𝑦3 4 11 𝜕 𝑥3 𝜕𝑦 (1) ( ) 𝜕 2 𝜓𝑦 ( ) 𝜕 2 𝜓𝑦 𝜕 4 𝜓𝑦 1 (4) 𝜕 4 𝜓𝑦 𝜕 4 𝜓𝑥 1 3 1 1 1 𝜕𝑐𝑧 (2) (4) + (4) − 3𝐵11 𝜓𝑦 + 𝐵11 + 3(2) + 𝐴(4) + (2) − (4) − 11 + 𝐹𝑦(1) + 𝑐𝑥(0) + =0 (25) 11 11 11 11 11 3 2 2 4 2 2 4 4 4 4 2 2 𝜕𝑥 𝜕 𝑥𝜕 𝑦 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑥 𝜕𝑦

𝛿𝜃𝑦 ∶

(1) 𝐴(1) + 𝐵11 12

( ) ) 𝜕 𝜃𝑦 2 2 𝜕 3 𝜃𝑦 𝜕 𝜃𝑥 𝜕 𝜃𝑥 𝜕 3 𝜃𝑥 ( (1) 1 𝜕 4 𝑤 1 (1) 𝜕 4 𝑤 1 1 (1) 𝜕 𝑤 (1) 𝜕 𝑤 (1) 𝜕 𝜃𝑥 (1) − 𝐴(0) + 𝐵11 + 𝐵11 − (1) − 11 4 + −𝐴(1) + 𝐵11 + (1) + −𝐴12 + 𝐵11 + (1) − 𝐴(0) 𝜃 12 𝜕𝑥 12 11 𝑧 𝜕𝑥 4 11 𝜕 𝑥4 4 𝜕𝑥 4 11 𝜕 𝑥3 𝜕𝑦 4 11 𝜕 𝑦3 𝜕 𝑥2 𝜕 𝑦2 𝜕𝑦 ( ) 𝜕2 𝜃 ( ) 𝜕2 𝜃 ) 𝜕𝜙 ) 𝜕 𝜙𝑦 𝜕 4 𝜃𝑧 1 (2) 𝜕 4 𝜃𝑧 ( (2) 𝜕 3 𝜙𝑥 ( (2) 1 1 1 1 1 1 𝑧 𝑧 𝑥 (2) (2) (2) (2) + 𝐵11 + (0) + 𝐵11 + (0) − (2) − 11 + −𝐴12 + 2𝐵11 − (0) + (2) + −𝐴12 + 2𝐵11 − (0) 4 11 𝜕 𝑥2 4 11 𝜕 𝑦2 4 11 𝜕 𝑥4 4 2 11 𝜕𝑥 2 11 𝜕 𝑥3 2 11 𝜕𝑦 𝜕 𝑦4 ( ) 𝜕2 𝜙 ( ) 𝜕2 𝜙 ) 𝜕𝜓 𝜕 3 𝜙𝑦 𝜕 4 𝜙𝑧 1 (3) 𝜕 4 𝜙𝑧 ( (3) 𝜕 3 𝜓𝑥 1 1 1 1 3 3 𝑧 𝑧 𝑥 (3) (3) (3) + (2) − 2𝐴(1) 𝜙𝑧 + 𝐵11 + (1) + 𝐵11 + (1) − (3) − 11 + −𝐴12 + 3𝐵11 − (1) + (3) 11 11 11 11 11 11 11 3 2 2 4 4 2 2 2 4 4 2 𝜕𝑥 4 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕 𝑥3 (1) 3 (1) ( ) 𝜕 𝜓𝑦 𝜕 𝜓𝑦 3 3 1 𝜕𝑐𝑦 1 𝜕𝑐𝑥 (3) + −𝐴(3) + 3𝐵11 − (1) + (3) + 𝐹𝑧(1) + + =0 (26) 12 2 11 𝜕𝑦 4 11 𝜕 𝑦3 2 𝜕𝑥 2 𝜕𝑦

𝛿𝜃𝑧 ∶ −𝐴(0) 12

163

V.S. Khorasani, M. Bayat

Engineering Analysis with Boundary Elements 94 (2018) 159–171

( ) 2 ( ) 2 𝜕2 𝑢 1 1 𝜕4 𝑢 1 1 𝜕4 𝑣 1 𝜕4 𝑣 1 𝜕3 𝑤 𝜕 𝑢 1 (2) 𝜕 4 𝑢 𝜕 𝑣 (2) (2) (1) 𝜕𝑤 + 𝐵11 + (0) − 11 − (2) + 𝐴(2) + 𝐵11 − (0) + (2) + (2) − 2𝐵11 − (1) 11 11 12 11 11 11 11 2 2 2 2 4 3 3 2 4 4 2 𝜕 𝑥𝜕 𝑦 4 𝜕𝑥 4 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕 𝑥 𝜕𝑦 4 𝜕 𝑥𝜕 𝑦 𝜕 𝑥3 2 4 4 ( ) 𝜕2 𝜃 ( ) 𝜕2 𝜃 ( ) 4 4 𝜕 𝜃𝑦 𝜕 𝜃𝑦 𝜕 𝜃𝑦 𝜕 𝜃𝑥 𝜕 𝜃𝑥 1 3 1 1 3 1 1 (1) (3) (3) 𝑥 𝑥 −2𝐵11 𝜃𝑥 + 𝐴(3) + (1) + 𝐵11 + (1) − (3) − (3) + 𝐴(3) + 𝐵11 − (1) + (3) + (3) 11 12 2 11 𝜕 𝑥2 2 11 𝜕 𝑦2 4 11 𝜕 𝑥2 𝜕 𝑦2 4 11 𝜕 𝑦4 2 11 𝜕 𝑥𝜕 𝑦 4 11 𝜕 𝑥3 𝜕𝑦 4 11 𝜕 𝑥𝜕 𝑦3 ( ) 𝜕𝜃 ) ( ) 𝜕2 𝜙 ( ) 𝜕2 𝜙 𝜕 3 𝜃𝑧 ( 𝜕 4 𝜙𝑥 𝜕 4 𝜙𝑥 1 1 1 1 (2) (2) (4) 𝑧 𝑥 𝑥 𝜙𝑥 + 𝐴(4) + 𝐴(2) − 2𝐵11 + (0) − (2) + −4𝐵11 − (0) + (2) + 𝐵11 + 3(2) − (4) − (4) 12 11 11 11 11 2 11 𝜕𝑥 2 11 𝜕 𝑥3 4 11 𝜕 𝑥2 𝜕 𝑦2 4 11 𝜕 𝑦4 𝜕 𝑥2 𝜕 𝑦2 ( ) 𝜕 2 𝜙𝑦 ( ) 𝜕𝜙 ) ( ) 𝜕2 𝜓 𝜕 4 𝜙𝑦 𝜕 4 𝜙𝑦 𝜕 3 𝜙𝑧 ( 1 1 1 3 (4) (3) (3) 𝑧 𝑥 + 𝐴(4) + 𝐵11 − 3(2) + (4) + (4) + 2𝐴(3) − 2𝐵11 + (1) − (3) + −6𝐵11 − 3(1) + (3) 𝜓𝑥 + 𝐴(5) 12 11 𝜕 𝑥𝜕 𝑦 11 11 12 11 11 11 11 11 3 3 3 4 4 𝜕𝑥 2 2 𝜕 𝑥 𝜕𝑦 𝜕 𝑥𝜕 𝑦 𝜕𝑥 𝜕 𝑥2 (2) ( ) 𝜕2 𝜓 ( ) 𝜕 2 𝜓𝑦 𝜕 4 𝜓𝑦 𝜕 4 𝜓𝑦 𝜕 4 𝜓𝑥 𝜕 4 𝜓𝑥 1 1 1 1 1 1 𝜕𝑐𝑧 (5) (5) 𝑥 + 𝐵11 + 5(3) − (5) − (5) + 𝐴(5) + 𝐵11 − 5(3) + (5) + (5) + 𝐹𝑥(2) + 𝑐𝑦(1) + =0 (27) 11 11 11 12 11 11 11 4 𝜕 𝑥𝜕 𝑦 4 2 2 𝜕𝑦 𝜕 𝑦2 𝜕 𝑥2 𝜕 𝑦2 4 𝜕 𝑦4 𝜕 𝑥3 𝜕𝑦 4 𝜕 𝑥𝜕 𝑦3

𝛿𝜙𝑥 ∶ 𝐴(2) 11

) ( ) 2 𝜕 𝑣 1 (0) 𝜕 2 𝑢 1 𝜕4 𝑢 1 𝜕4 𝑢 1 𝜕 2 𝑣 1 (2) 𝜕 4 𝑣 1 (2) 𝜕 4 𝑣 1 𝜕3 𝑤 (2) (1) 𝜕𝑤 11 + (2) + (2) + 𝐵11 + (0) + 𝐴(2) − 11 −  − 2𝐵11 − (1) 11 11 11 11 𝜕 𝑦2 2 𝜕 𝑥𝜕 𝑦 4 2 4 𝜕𝑦 2 11 𝜕 𝑦3 𝜕 𝑥3 𝜕𝑦 4 𝜕 𝑥𝜕 𝑦3 𝜕 𝑥2 𝜕 𝑥4 4 11 𝜕 𝑥2 𝜕 𝑦2 ( ) 𝜕2 𝜃 ( ) 𝜕 2 𝜃𝑦 ( ) 𝜕 2 𝜃𝑦 𝜕 4 𝜃𝑦 1 (3) 𝜕 4 𝜃𝑦 𝜕 4 𝜃𝑥 𝜕 4 𝜃𝑥 3 1 1 3 1 1 (3) (1) (3) 𝑥 + 𝐴(3) + 𝐵11 − (1) + (3) + (3) − 2𝐵11 𝜃𝑦 + 𝐵11 + (1) + 𝐴(3) + (1) − (3) − 11 12 11 2 11 𝜕 𝑥𝜕 𝑦 4 11 𝜕 𝑥3 𝜕𝑦 4 11 𝜕 𝑥𝜕 𝑦3 2 11 𝜕 𝑥2 2 11 𝜕 𝑦2 4 11 𝜕 𝑥4 4 𝜕 𝑥2 𝜕 𝑦2 ( ) 𝜕𝜃 ( ) ( ) ( ) 𝜕 2 𝜙𝑦 3 2 4 4 𝜕 𝜃𝑧 𝜕 𝜙𝑥 𝜕 𝜙𝑥 1 (4) 𝜕 𝜙𝑥 1 1 1 (2) (4) (2) (4) 𝑧 𝜙𝑦 + 𝐵11 + 𝐴(2) − 2𝐵11 + (0) − (2) + 𝐴(4) + 𝐵11 − 3(2) + 11 + (4) + −4𝐵11 − (0) + 3(2) 12 11 11 12 11 11 11 11 3 3 3 2 𝜕𝑦 4 𝜕 𝑥𝜕 𝑦 4 𝜕𝑦 𝜕 𝑥 𝜕𝑦 4 𝜕 𝑥𝜕 𝑦 𝜕 𝑥2 2 4 4 ( ) 𝜕 𝜙𝑦 ( ) 𝜕𝜙 ) 𝜕2 𝜓 𝜕 𝜙𝑦 1 (4) 𝜕 𝜙𝑦 𝜕 3 𝜙𝑧 ( (5) 𝜕 4 𝜓𝑥 1 1 1 (3) (5) 𝑧 𝑥 + 𝐴(4) + (2) − (4) − 11 + 2𝐴(3) − 2𝐵11 + (1) − (3) + 𝐴12 + 𝐵11 − 5(3) + (5) 11 11 12 11 11 𝜕 𝑥𝜕 𝑦 4 11 𝜕 𝑥4 4 𝜕𝑦 4 11 𝜕 𝑦3 4 11 𝜕 𝑥3 𝜕𝑦 𝜕 𝑦2 𝜕 𝑥2 𝜕 𝑦2 (2) ( ) ( ) 𝜕 2 𝜓𝑦 ( ) 𝜕 2 𝜓𝑦 𝜕 4 𝜓𝑦 1 (5) 𝜕 4 𝜓𝑦 𝜕 4 𝜓𝑥 1 3 1 1 1 𝜕𝑐𝑧 (3) (5) 𝜓𝑦 + 𝐵11 + (5) + −6𝐵11 − 3(1) + 5(3) + 𝐴(5) + (3) − (5) − 11 + 𝐹𝑦(2) + 𝑐𝑥(1) + =0 (28) 11 11 11 11 11 2 2 4 2 2 4 11 𝜕 𝑥𝜕 𝑦3 2 4 4 2 2 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑥 𝜕𝑦

𝛿𝜙𝑦 ∶

(

(2) + 𝐵11 − 𝐴(2) 12

( ) ) 2 2 𝜕 𝜃𝑥 𝜕 𝜃𝑥 𝜕 3 𝜃𝑥 ( 1 𝜕 4 𝑤 1 (2) 𝜕 4 𝑤 1 (2) 𝜕 𝑤 (2) 𝜕 𝑤 (2) 𝜕 𝜃𝑥 (2) 𝜕 𝜃𝑦 − 2𝐴(1) + 𝐵11 + 𝐵11 − (2) − 11 + −2𝐴(2) + 𝐵11 + (2) + −2𝐴(2) + 𝐵11 12 11 12 11 12 2 2 4 4 3 𝜕𝑥 𝜕𝑥 4 4 𝜕𝑥 4 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕𝑥 3𝜃 ( ) ( ) ( ) 2𝜃 2𝜃 4𝜃 4𝜃 𝜕 𝜕 𝜕 𝜙 𝜕 3 𝜙𝑥 𝜕 𝜕 𝜕 𝑦 1 1 1 1 1 1 (3) (3) (3) 𝑧 𝑧 𝑧 𝑧 𝑥 + (2) − 2𝐴(1) 𝜃 + 𝐵11 + (1) + 𝐵11 + (1) − (3) − (3) + −2𝐴(3) + 2𝐵11 − (1) + (3) 11 𝑧 12 11 4 11 𝜕 𝑦3 2 11 𝜕 𝑥2 2 11 𝜕 𝑦2 4 11 𝜕 𝑥4 4 11 𝜕 𝑦4 𝜕𝑥 2 11 𝜕 𝑥3 ( ) 𝜕 𝜙𝑦 ( ) 𝜕2 𝜙 ( ) 𝜕2 𝜙 𝜕 3 𝜙𝑦 𝜕 4 𝜙𝑧 1 (4) 𝜕 4 𝜙𝑧 1 1 (3) (4) (4) 𝑧 𝑧 + −2𝐴(3) + 2𝐵11 − (1) + (3) − 4𝐴(2) 𝜙𝑧 + 𝐵11 + (2) + 𝐵11 + (2) − (4) − 11 12 11 11 11 11 11 11 3 2 2 𝜕𝑦 2 4 4 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕 𝑥4 𝜕 𝑦4 (2) 3 (2) ( ) 𝜕𝜓 ) 𝜕 𝜓𝑦 𝜕 𝜓𝑦 𝜕 3 𝜓𝑥 ( 3 3 1 𝜕𝑐𝑦 1 𝜕𝑐𝑥 (4) (4) 𝑥 + −2𝐴(4) + 3𝐵11 − 3(2) + (4) + −2𝐴(4) + 3𝐵11 − 3(2) + (4) + 𝐹𝑧(2) + + =0 12 11 12 11 𝜕𝑥 4 11 𝜕 𝑥3 𝜕𝑦 4 11 𝜕 𝑦3 2 𝜕𝑥 2 𝜕𝑦

𝛿𝜙𝑧 ∶ −2𝐴(1) 12

(29)

( ) 2 ( ) 2 𝜕2 𝑢 3 1 𝜕4 𝑢 3 1 𝜕4 𝑣 1 𝜕4 𝑣 3 𝜕3 𝑤 𝜕 𝑢 1 (3) 𝜕 4 𝑢 𝜕 𝑣 (3) (3) (2) 𝜕𝑤 + 𝐵11 + (1) −  − (3) + 𝐴(3) + 𝐵11 − (1) + (3) + (3) − 3𝐵11 − (2) 12 2 11 𝜕 𝑦2 4 11 𝜕 𝑥2 𝜕 𝑦2 4 11 𝜕 𝑦4 2 11 𝜕 𝑥𝜕 𝑦 4 11 𝜕 𝑥3 𝜕𝑦 4 11 𝜕 𝑥𝜕 𝑦3 𝜕𝑥 4 11 𝜕 𝑥3 𝜕 𝑥2 ( ) 𝜕2 𝜃 ( ) 𝜕2 𝜃 ) 𝜕 2 𝜃𝑦 𝜕 4 𝜃𝑦 𝜕 4 𝜃𝑦 𝜕 4 𝜃𝑥 𝜕 4 𝜃𝑥 ( (4) 3 1 1 1 1 (2) (4) (4) 𝑥 𝑥 − 3𝐵11 𝜃𝑥 + 𝐴(4) + (2) + 𝐵11 + 3(2) − (4) − (4) + 𝐴12 + 𝐵11 − 3(2) + (4) + (4) 11 11 11 11 11 11 11 11 2 2 2 2 4 3 4 4 4 𝜕 𝑥𝜕 𝑦 4 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕 𝑥 𝜕𝑦 4 𝜕 𝑥𝜕 𝑦3 ( ) 𝜕𝜃 ( ) ( ) ( ) 3 2 2 4 4 𝜕 𝜃𝑧 𝜕 𝜙𝑥 𝜕 𝜙𝑥 𝜕 𝜙𝑥 1 (5) 𝜕 𝜙𝑥 3 3 3 1 (3) (3) (5) 𝑧 + 𝐴(3) − 3𝐵11 + (1) − (2) + −6𝐵11 − 3(1) + (3) + 𝐵11 + 5(3) − 11 − (5) 𝜙𝑥 + 𝐴(5) 12 11 11 11 2 11 𝜕𝑥 4 11 𝜕 𝑥3 2 11 𝜕 𝑥2 4 𝜕 𝑦2 𝜕 𝑥2 𝜕 𝑦2 4 11 𝜕 𝑦4 4 4 ( ) 2 ( ) 𝜕𝜙 ) ( ) 𝜕2 𝜓 𝜕 3 𝜙𝑧 ( 1 (5) 𝜕 𝜙𝑦 1 (5) 𝜕 𝜙𝑦 3 9 (5) (3) 𝜕 𝜙𝑦 (4) (4) 𝑧 𝑥 + 𝐴(5) + 𝐵 − 5  +  +  + 2𝐴(4) − 3𝐵11 + 3(2) − (3) + −9𝐵11 − 9(2) + (4) 𝜓𝑥 + 𝐴(6) 12 11 11 𝜕 𝑥𝜕 𝑦 12 11 11 11 4 11 𝜕 𝑥3 𝜕𝑦 4 11 𝜕 𝑥𝜕 𝑦3 𝜕𝑥 4 11 𝜕 𝑥3 4 11 𝜕 𝑥2 (3) ( ) 2 ) 2 4 𝜕 4 𝜓𝑦 𝜕 4 𝜓𝑦 𝜕 4 𝜓𝑥 ( (6) 15 (4) 𝜕 𝜓𝑥 1 (6) 𝜕 𝜓𝑥 1 15 (4) 𝜕 𝜓𝑦 1 1 3 1 𝜕𝑐𝑧 (6) (6) + 𝐵11 + 11 − 11 − (6) + 𝐴12 + 𝐵11 − 11 + (6) + (6) + 𝐹𝑥(3) + 𝑐𝑦(2) + =0 11 11 11 2 2 2 4 3 3 2 4 4 2 𝜕 𝑥𝜕 𝑦 4 2 2 𝜕𝑦 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕 𝑥 𝜕𝑦 4 𝜕 𝑥𝜕 𝑦 (30)

𝛿𝜓𝑥 ∶ 𝐴(3) 11

) ( ) 2 3 (1) 𝜕 2 𝑢 1 𝜕4 𝑢 1 𝜕4 𝑢 3 𝜕 2 𝑣 1 (3) 𝜕 4 𝑣 1 (3) 𝜕 4 𝑣 3 𝜕3 𝑤 𝜕 𝑣 (3) (2) 𝜕𝑤 11 + (3) + (3) + 𝐵11 + (1) + 𝐴(3) − 11 −  − 3𝐵11 − (2) 11 11 11 11 𝜕 𝑦2 2 𝜕 𝑥𝜕 𝑦 4 2 4 𝜕𝑦 4 11 𝜕 𝑦3 𝜕 𝑥3 𝜕𝑦 4 𝜕 𝑥𝜕 𝑦3 𝜕 𝑥2 𝜕 𝑥4 4 11 𝜕 𝑥2 𝜕 𝑦2 4 4 ( ) 2 ( ) 2 ( ) 2 4 4 1 (4) 𝜕 𝜃𝑥 1 (4) 𝜕 𝜃𝑥 3 (2) 𝜕 𝜃𝑦 1 (4) 𝜕 𝜃𝑦 1 (4) 𝜕 𝜃𝑦 (4) (2) 𝜕 𝜃𝑥 (2) (4) (2) 𝜕 𝜃𝑦 (4) + 𝐴(4) + 𝐵 − 3  +  +  − 3 𝐵 𝜃 + + 3  + +  −  −  𝐵 𝐴 𝑦 12 11 11 𝜕 𝑥𝜕 𝑦 11 11 11 𝜕 𝑥2 11 4 11 𝜕 𝑥3 𝜕𝑦 4 11 𝜕 𝑥𝜕 𝑦3 4 11 𝜕 𝑦2 4 11 𝜕 𝑥4 4 11 𝜕 𝑥2 𝜕 𝑦2 ( ) 𝜕𝜃 ( ) ( ) ( ) 𝜕 2 𝜙𝑦 3 2 4 4 𝜕 𝜃𝑧 𝜕 𝜙𝑥 𝜕 𝜙𝑥 1 (5) 𝜕 𝜙𝑥 3 3 1 (3) (5) (3) (5) 𝑧 + 𝐴(3) − 3𝐵11 + (1) − (3) + 𝐴(5) + 𝐵11 − 5(3) + 11 + (5) + −6𝐵11 − 3(1) + 5(3) 𝜙𝑦 + 𝐵11 12 11 11 12 11 11 11 11 3 3 3 2 𝜕𝑦 4 𝜕 𝑥𝜕 𝑦 4 𝜕𝑦 𝜕 𝑥 𝜕𝑦 4 𝜕 𝑥𝜕 𝑦 𝜕 𝑥2 4 4 ( ) 𝜕 2 𝜙𝑦 ( ) ( ) 3 2 4 𝜕 𝜙𝑦 1 (5) 𝜕 𝜙𝑦 𝜕 𝜙𝑧 3 (4) 𝜕 𝜙𝑧 3 1 15 (4) 𝜕 𝜓𝑥 1 (6) 𝜕 𝜓𝑥 (4) (6) + 𝐴(5) + (3) − (5) − 11 + 2𝐴(4) − 3𝐵11 + 3(2) − 11 + 𝐴(6) + 𝐵11 −  +  11 12 11 12 2 11 𝜕 𝑦2 4 11 𝜕 𝑥4 4 𝜕𝑦 4 2 11 𝜕 𝑥𝜕 𝑦 4 11 𝜕 𝑥3 𝜕𝑦 𝜕 𝑥2 𝜕 𝑦2 𝜕 𝑦3 (3) ( ) ( ) 2 ( ) 𝜕 2 𝜓𝑦 𝜕 4 𝜓𝑦 1 (6) 𝜕 4 𝜓𝑦 𝜕 4 𝜓𝑥 1 15 (4) 𝜕 𝜓𝑦 9 1 1 1 𝜕𝑐𝑧 (4) (6) + (6) + −9𝐵11 − 9(2) +  + 𝐴(6) + (4) − (6) − 11 + 𝐹𝑦(3) + 𝑐𝑥(2) + = 0 (31) 𝜓𝑦 + 𝐵11 11 11 2 2 4 11 𝜕 𝑥𝜕 𝑦3 2 11 𝜕 𝑥2 4 11 𝜕 𝑦2 4 11 𝜕 𝑥4 4 2 2 𝜕𝑥 𝜕𝑥 𝜕𝑦

𝛿𝜓𝑦 ∶

(

(3) + 𝐵11 − 𝐴(3) 12

164

V.S. Khorasani, M. Bayat

Engineering Analysis with Boundary Elements 94 (2018) 159–171

where: (𝑖)

𝐹𝜉(𝑖) = 𝑓 𝜉 +

𝑐𝜉(𝑖)

=

𝑐 (𝜉𝑖)

( )𝑖 [ ] ℎ 𝑞𝜉𝑡 + (−1)𝑖 𝑞𝜉𝑏 , 2

( )𝑖 [ ] + ℎ2 𝑝𝑡𝜉 + (−1)𝑖 𝑝𝑏𝜉 ,

𝑓

𝑐

(𝑖)

(𝑖)

=

=

ℎ 2 ∫− ℎ 2

ℎ 2 ∫− ℎ 2

(𝑧)𝑖 𝑓𝜉 𝑑𝑧

𝑖

(𝑧) 𝑐𝜉 𝑑𝑧

[ 𝑥 − 𝑧 𝑦 − 𝑧 (𝑥 − 𝑧)2 (𝑥 − 𝑧)(𝑦 − 𝑧) (𝑦 − 𝑧)2 𝐩𝑇 (𝐱) = 1, , , , , ℎ ℎ ℎ2 ℎ2 ℎ2

(32)

×

(33)

] (𝑥 − 𝑧)3 (𝑥 − 𝑧)2 (𝑦 − 𝑧) (𝑥 − 𝑧)(𝑦 − 𝑧)2 (𝑦 − 𝑧)3 , , , , for 𝑙 ≠ 0 ℎ3 ℎ3 ℎ3 ℎ3 (39)

The coefficient vector a(k) (x) can be determined by minimizing a weight discrete L2 norm [45,46]:

4. The MLPG solving procedure For employing the MLPG method as a local approximation, we need to choose trial and test functions and proper subdomain geometry for each of them, to write local weak form of governing PDEs. Then, by applying suitable numerical integration scheme with adequate numbers of integration points, stiffness matrix can be obtained. Finally, by solving linear equations system that consisting of the stiffness matrix and the applied forces vector, we reach to the solution of the bending problem.

𝐽 (𝑘) (𝐱) =

𝑁 ∑ 𝑖=1

[ ]2 [ ]𝑇 𝐰𝑖 (𝐱) 𝐩𝑇 (𝐱𝑖 )𝐚(𝑘) (𝐱) − 𝑈𝑖(𝑘) = 𝑃 ⋅ 𝐚(𝑘) (𝐱) − 𝐔(𝑘)

[ ] ⋅𝑊 ⋅ 𝑃 ⋅ 𝐚(𝑘) (𝐱) − 𝐔(𝑘)

(40)

where wi (x) is the weight function of node i, N is number of nodes in Ωx for which the weight function wi (x) > 0, and 𝑈𝑖(𝑘) are the fictitious nodal values for unknown trial function. The matrices P and W are following:

4.1. The MLS approximation ⎡𝐩𝑇 (𝐱1 )⎤ ⎢ 𝑇 2⎥ 𝐩 (𝐱 )⎥ 𝑃 =⎢ ⎢ ⋮ ⎥ ⎢𝐩𝑇 (𝐱𝑛 )⎥ ⎣ ⎦

In present work, the MLS approximation is used to represent the trial functions for estimation of unknown variables at nodes in the problem domain. Let consider that Ωx is a subdomain of Ω in the neighborhood of a point x for defining the MLS approximation of the trial function around x. The MLS approach is based on the approximation of polynomial (or monomial) of order m (power basis function) and on the choice of a set of weight functions. These functions depend on the distance between the considered node x and neighborhood node xj as following: ‖ ‖ ‖𝑥 − 𝑥 𝑗 ‖ ‖, 𝑟𝑗 (𝑥) = ‖ 𝐷

[ ] 𝐿 𝑟 ∈ 0, 𝐷

∀𝐱 ∈ Ω𝑥

(34)

(41)

𝑁×𝑚

(42)

𝐴(𝐱) = 𝑃 𝑇 𝑊 𝑃 = 𝐵(𝐱)𝑃 =

𝑁 ∑ 𝑖=1

𝐰𝑖 (𝐱)𝐩(𝐱𝑖 )𝐩𝑇 (𝐱𝑖 ),

[ ] 𝐵(𝐱) = 𝑃 𝑇 𝑊 = 𝐰1 (𝐱)𝐩(𝐱1 ), 𝐰2 (𝐱)𝐩(𝐱2 ), … , 𝐰𝑛 (𝐱)𝐩(𝐱𝑛 ) .

(43) (44)

The matrix A that called the moment matrix, is of size m × m. When A is non-singular, the MLS approximation is well-defined. In other word, rank of P must be equal to m, that means at least m weight functions should be non-zero for a well-defined MLS approximation. By substituting a(k) (x) from Eq. (42) into Eq. (35), we have:

(35)

𝑈ℎ(𝑘) (𝐱) = Φ𝑇 (𝐱) ⋅ 𝐔(𝑘) =

𝑁 ∑ 𝑖=1

𝜑𝑖 (𝐱)𝑈𝑖(𝑘) ,

𝐱 ∈ Ω𝑥

(45)

where Φ𝑇 (𝐱) = 𝐩𝑇 (𝐱)𝐴−1 (𝐱)𝐵(𝐱),

= [1, 𝑥, 𝑦, 𝑥𝑦, for 𝑙 = 0 𝐩𝑇 (𝐱) = [1, 𝑥, 𝑦, 𝑥2 , 𝑥𝑦, 𝑦2 , 𝑥3 , 𝑥2 𝑦, 𝑥𝑦2 , 𝑦3 ], for 𝑙 ≠ 0

𝑦2 ],

(46)

or 𝜑𝑖 (𝐱) =

The MLS approximation can be implemented in a more stable fashion, if a shifted and scaled polynomial basis functions is used as basis for approximation [42]. In this paper, we use the following basis: ⎧ ⎫ ⎪ (𝐱 − 𝑧)𝛼 ⎪ ⎨ |𝛼| ⎬ ⎪ ℎ𝑋,Ω ⎪ ⎩ ⎭0≤|𝛼|≤𝑚

0 ⎤ ⋮ ⎥ ⎥ 𝐰𝑛 (𝐱)⎦

where matrices A(x) and B(x) are:

where 1 (x), p2 (x),…, pm (x)] is a complete monomial basis or order m, and a(k) (x) is a vector with components 𝑎(𝑗𝑘) (𝐱), 𝑗 = 1, 2, … , 𝑚, which are functions of the space coordinates x = (x, y). U can be each of the unknown generalized displacements. Here, we use two different m-order monomial depending on value of length scale parameter, l, as: 𝑥2 ,

… ⋱ ⋯

𝐴(𝐱)𝐚(𝑘) (𝐱) = 𝐵(𝐱)𝐔(𝑘)

pT (x) = [ p

𝐩𝑇 (𝐱)

⎡𝐰1 (𝐱) 𝑊 =⎢ ⋮ ⎢ ⎣ 0

The stationary point of J in Eq. (40), with respect to a(k) (x), form the linear relation between a(k) (x) and U(k) as:

where D is the influence radius of domain and L is the domain length [39]. To approximate the distribution of each eleven trial functions for unknown generalized displacements in Ωx over a number of nodes {xi }, i = 1, 2, …, N, the MLS approximation 𝑈ℎ(𝑘) (𝐱) of U(k) is defined as [40,41]: 𝑈ℎ(𝑘) (𝐱) = 𝐩𝑇 (𝐱)𝐚(𝑘) (𝐱),

,

𝑚 ∑ 𝑗=1

[ ] 𝑝𝑗 (𝐱) 𝐴−1 (𝐱)𝐵(𝐱) 𝑗𝑖

(47)

For each nodal points xi , 𝜑i (x) is the shape function of the MLS approximation. In this paper, the Gaussian weight function is applied for computing shape functions around nodes: { exp[−(𝑑𝑖 ∕𝑐𝑖 )2 ]−exp[−(𝑟𝑖 ∕𝑐𝑖 )2 ] 0 ≤ 𝑑𝑖 ≤ 𝑟𝑖 1−exp[−(𝑟𝑖 ∕𝑐𝑖 )2 ] 𝐰𝑖 (𝐱) = (48) 0 𝑑𝑖 > 𝑟𝑖

(36)

where z is fixed and depends on the evaluation point to be considered and hX,Ω is the fill distance that for a set of points X = {x1 ,x2 ,…, xn } in a bounded domain Ω is defined as: ‖ ‖ ℎ𝑋,Ω = sup min ‖𝐱 − 𝑥𝑗 ‖ (37) ‖2 𝑥∈Ω 1≤𝑗≤𝑁 ‖

where ri is the size of the support domain, ci is the shape parameter of the weight function wi , and di = ‖x − xi ‖. To ensure that matrix A is non-singular, ri should be large enough to collect adequate number of nodes in the domain of definition of every sample point (N > m). Here, the radius of support domain, ri , is selected as following:

By using scaled and shifted basis functions to overcome the inherent instability of MLS approximation, the relations for pT (x) are updated as following [43,44]: [ ] 𝑥 − 𝑧 𝑦 − 𝑧 (𝑥 − 𝑧)2 (𝑥 − 𝑧)(𝑦 − 𝑧) (𝑥 − 𝑧)2 𝐩𝑇 (𝐱) = 1, , , , , , for 𝑙 = 0 (38) ℎ ℎ ℎ2 ℎ2 ℎ2

𝑟𝑖 = 2𝒎ℎ

(49)

where m is polynomial degree of basis function and h is fill distance (m = 2 for l = 0 and m = 3 for l≠0). The reason for this selection is that 165

V.S. Khorasani, M. Bayat

Engineering Analysis with Boundary Elements 94 (2018) 159–171

equations of GTPT earlier, to avoid a lengthy paper, we introduce only first local weak equation of 121 local weak equations, other equations can be obtained in same manner according to Eqs. (21)–(31). Thus first local weak equation (first submatrix of K) is: 𝐾𝑖𝑗0101 = −

(1) (1) 𝜕𝑤(1) ⎛ 𝜕𝜑(1) 𝜕𝑤𝑗 𝑗 (0) 𝜕𝜑𝑖 ⎜𝐴(0) 𝑖 + 𝐵 11 𝜕𝑦 ∫Ω ⎜ 11 𝜕𝑥 𝜕𝑥 𝜕𝑦 ⎝

⎛ 2 (1) 𝜕 2 𝑤(1) ⎞ ⎞ 𝜕 2 𝑤(1) 𝜕 2 𝜑(1) 𝑗 𝑗 ⎟⎟ 1 (0) 𝜕 𝜑𝑖 𝑖 + ⎜(0) +  𝑑Ω 11 𝜕 𝑦2 4 ⎜ 11 𝜕 𝑥𝜕 𝑦 𝜕 𝑥𝜕 𝑦 𝜕 𝑦2 ⎟ ⎟ ⎝ ⎠⎠

(51)

The boundary conditions are imposed directly by the MLS approximation for boundary nodes. The boundary conditions for simplysupported rectangular plate in this problem are following: 𝑢(𝑥, 0) = 𝑢(𝑥, 𝑏) = 0,

𝜃𝑥 (𝑥, 0) = 𝜃𝑥 (𝑥, 𝑏) = 0,

𝜙𝑥 (𝑥, 0) = 𝜙𝑥 (𝑥, 𝑏) = 0, 𝑣(0, 𝑦) = 𝑣(𝑎, 𝑦) = 0,

𝜓𝑥 (𝑥, 0) = 𝜓𝑥 (𝑥, 𝑏) = 0, 𝜃𝑦 (0, 𝑦) = 𝜃𝑦 (𝑎, 𝑦) = 0,

𝜙𝑦 (0, 𝑦) = 𝜙𝑦 (𝑎, 𝑦) = 0,

𝜓𝑦 (0, 𝑦) = 𝜓𝑦 (𝑎, 𝑦) = 0,

𝑤(𝑥, 0) = 𝑤(𝑥, 𝑏) = 0,

𝑤(0, 𝑦) = 𝑤(𝑎, 𝑦) = 0,

𝜃𝑧 (𝑥, 0) = 𝜃𝑧 (𝑥, 𝑏) = 0,

𝜃𝑧 (0, 𝑦) = 𝜃𝑧 (𝑎, 𝑦) = 0,

𝜙𝑧 (𝑥, 0) = 𝜙𝑧 (𝑥, 𝑏) = 0,

𝜙𝑧 (0, 𝑦) = 𝜙𝑧 (𝑎, 𝑦) = 0,

(𝑖) 𝑀𝑥𝑥 (0, 𝑦)

=

(𝑖) 𝑀𝑥𝑥 (𝑎, 𝑦)

= 0,

𝑗) 𝑗) (𝑥𝑦 (0, 𝑦) = (𝑥𝑦 (𝑎, 𝑦) = 0,

Fig. 2. Comparison between results of various types of weight functions as basis of the MLS approximation.

𝑗) 𝑗) (𝑥𝑦 (𝑥, 0) = (𝑥𝑦 (𝑥, 𝑏) = 0,

(52)

where i = 0, 1, 2, 3 and j = 0, 1, 2. For applying the boundary conditions, we first separated boundary test nodes from inner test nodes. Then we calculated MLPG solutions for all trial nodes and inner test nodes (trial nodes coordinates are considered as like as test nodes). For example, for 11 × 11 distribution of test nodes in a rectangular domain, we separated 40 boundary test nodes from 81 inner test nodes, and then achieved MLPG solutions of 121 trial nodes and 81 test nodes. In this step we got a 81 × 121 matrix. The governing equations are a system of PDEs with 11 dependent displacement variables that mean we have boundary conditions about 11 variables. According to relation (52) we have two kind of boundary conditions for bending of considered plate, zero for all four edges (40 nodes) for 3 variables (w, 𝜃 z ,𝜙z ), and zero for only two edges (22 nodes) for remaining 8 variables. In first case, we defined a 40 × 121 zero matrix and then attached it to 81 × 121 solution matrix to get 121 × 121 matrix for all nodes. In second case we used the MLS approximation procedure for non-zero boundary test nodes (18 nodes) and 121 trial nodes to get a 18 × 121 matrix, and then attached this matrix to a 22 × 121 zero matrix to get a 40 × 121 matrix for joining to related MLPG solution matrix. We repeated this algorithm for 121 local weak equations to gain the complete 1331 × 1331 stiffness matrix.

by increasing the polynomial degree, it is needed to more trial nodes for approximation, so the radius of support domain must be increased. A comparison between results of the Gauss function, the cubic spline, and quartic spline as basis for weight function for present problem is showed in Fig. 2. Differences between these three functions are little but as shown, the Gauss function yields to smoother results and thus is selected for present study. 4.2. The test function (MLPG type) Atluri and Shen [7] categorized the MLPG method with respect to type of the applied test function in discretizing equations. Among six types of the MLPG approach, the MLPG1 and the MLPG5 has received more attention in the literature due to their symmetry of local weak form and regularity of their integration. As mentioned earlier, the governing PDE of motion is of fourth order when l ≠ 0, thus we cannot use the MLPG5 for discretizing equations because in the MLPG5, test function is the Heaviside step function, which cannot be differentiated two times to reach proper weak equations. For this reason, we employed the MLPG1 for discretizing governing equations either l = 0 or l ≠ 0 for consistency in results. Thus in present study the test function is Gauss function as the MLS weight function.

4.4. The Gauss–Legendre quadrature rule The Gauss–Legendre quadrature rule for arbitrary interval [a, b] is: 𝑏

4.3. The local weak form

∫𝑎

According to the governing equations of GTPT (21–31), we have a PDE system with eleven equations for eleven unknown generalized displacements, thus the stiffness matrix, K, has 121 square submatrices which their row index in K is number of their equation in system, and their column index in K is number of the generalized displacement component in respect to other. So matrix form of the MLPG model for equation of motion is: [𝐊](11∗𝑁 )×(11∗𝑁 ) {𝐔}11∗𝑁 = {𝐅}11∗𝑁

(𝑖) (𝑖) 𝑀𝑦𝑦 (𝑥, 0) = 𝑀𝑦𝑦 (𝑥, 𝑏) = 0,

𝑓 (𝑦)𝑑𝑦 =

𝑛 𝑏−𝑎 ∑ 𝜔 𝑓 (𝑦𝑖 ) + 𝑅𝑛 2 𝑖=1 𝑖

(53)

where the nodes are 𝑦𝑖 = ( 𝑏−2 𝑎 )𝑥𝑖 + ( 𝑏+2 𝑎 ), and nodes xi are the ith zeros of Legendre polynomials Pn (x). The weights 𝜔i and the remainder Rn are following: 𝜔𝑖 = ( 𝑅𝑛 =

(50)

where U is the displacements vector and F is the force vector; both of them consist of 11 subvectors. Because we presented the governing

2 )[ ]2 , 1 − 𝑥2𝑖 𝑃 ′ 𝑛 (𝑥𝑖 )

(𝑏 − 𝑎)2𝑛+1 22𝑛+1 (𝑛!)4 (2𝑛 + 1)[(2𝑛)!]3

(54)

The Gauss–Legendre rules are mostly used in computer subroutine, because they have the highest possible degree of precision [47]. So 166

V.S. Khorasani, M. Bayat

Engineering Analysis with Boundary Elements 94 (2018) 159–171

Fig. 4. Non-dimensional deflection 𝑤(𝑥, 2𝑏 , 0)versus distance for various values of n and l = 0.

Fig. 3. The geometry of plate and distribution of nodes.

we choose this quadrature rule for numerical integration of local weak equations in problem domain. 5. Numerical results In this section, we present the results of static bending of a simplysupported FG plate by using the MLPG method. The plate thickness, h, is considered to be 17.6 × 10−6 m, and length of square plate, a, is assumed to be 20 h. The 11 by 11 node distribution with equal distances is used for discretizing the domain of micro-plate. The coordinates of nodes and domain of assumed plate is showed in Fig. 3. For both materials of FG plate, Poisson’s ratio, 𝜈, is assumed to be 0.38 and the Young modulus of top (metal) and bottom (ceramic) surfaces are considered as Et = 14.4 × 109 N/m2 and Eb = 1.44 × 109 N/m2 , respectively. In addition to the boundary conditions (52), the out-of-plane rigid body motion is constrained. For purpose of comparison the present results with the analytical solution [23], we assume uniformly distributed load value to be 𝑞𝑧𝑡 = 1N∕m2 on top surface, and no point load is considered thus 𝑐𝜉(𝑖) = 0. Therefore for right hand side of PDE system for the bending problem, we have: [ ]𝑇 2 𝑓 = 0 0 (55) 𝑞𝑧𝑡 0 0 𝑞𝑧𝑡 ℎ 0 0 𝑞𝑧𝑡 ℎ 0 0 2

4

The Gauss–Legendre quadrature approach is used as numerical integration method in present study with 15 integration points for computing both sides of the equations system. Fig. 4 show variation of non-dimensional deflections 𝑤(𝑥, 2𝑏 , 0) versus distances for various values of power law index, n, when l = 0. In Fig. 5, 𝑤(𝑥, 2𝑏 , 0) are shown when micro-structural effects are included, i.e. l = 17.6 × 10 − 6 m. When n increases, the volume fraction of top surface material decreases, so because we assumed Et > Eb , the stiffness of FG plate decreases and in both cases (l = 0orl ≠ 0), larger deflections values are obtained for larger values of n. Fig. 6 shows comparison between values of 𝑤(𝑥, 2𝑏 , 0) for same power law indexes, but different l values.

Fig. 5. Non-dimensional deflection 𝑤(𝑥, 2𝑏 , 0)versus distance for various values of n and l = 17.6 × 10 − 6 .

theory to the governing equations of GTPT. To have a better understanding of the variation of deflections on all domain of plate, Figs. 7 and 8 show 𝑤 in three-dimensional for all nodes which considered in problem. Figs. 9 and 10 show the bending stress, 𝜎 xx , and the transverse 𝑎 𝑏 shear stress, 𝜎 xz , which are computed for left corner node at ( 11 , 11 , 𝑧), respectively. Since we do not force the transverse shear strain to be zero

By considering microstructural size effects, additional plate stiffness (11𝑖) is included, thus stiffness matrix, K, becomes stiffer. Therefore for same values of n, smaller deflections are obtained by applying couple stress

167

V.S. Khorasani, M. Bayat

Engineering Analysis with Boundary Elements 94 (2018) 159–171

Fig. 8. Non-dimensional deflection 𝑤on all domain of plate for n = 1 and l = 17.6 × 10 − 6 .

Fig. 6. Non-dimensional deflection 𝑤(𝑥, 2𝑏 , 0)versus distance for various values of n and l.

Fig. 9. Bending stress through thickness for various values of n and l = 0. Fig. 7. Non-dimensional deflection 𝑤 on all domain of plate for n = 1 and l = 0.

that use the Heaviside step function and the Kronecker delta function as test function, respectively. Fig. 13 shows a comparison between results of MLPG1, MLPG2, and MLPG5 for a rectangular plate when l = 0 and n = 1. The local weak form in MLPG1 and MLPG5 is symmetric but in MLPG2 is asymmetric. To show the efficiency of the solving method that used in this research, we considered Poisson equation as a test problem in a square domain with the Dirichlet boundary conditions for all edges. The maximum error of obtained results with exact solutions that calculated with TM the Frank’s function is acquired about 0.0025. By using an Intel® Core i5-2410 M CPU with 2.30 GHz frequency, the calculation time for solving Poisson equation as test problem is obtained about 1.8 s. The CPU times for main problem in this research are presented in Table 1 for various values of n and l. In next section, we will compare

on top and bottom surfaces, thus the transverse shear stresses on both surfaces are not zero. In this research, it is focused on rectangular domain because the results can be compared to analytical solutions in the literature. But to ensure that the solving procedure (algorithms and codes) is implemented in the general manner, it is necessary to consider an irregular domain. For this purpose we define an L-shape domain by dimension close to considered rectangular domain. Fig. 11 shows nodes coordinates in the L-shape plate, and the results for the bending of simply-supported conditions for all edges are shown in Fig. 12. As it mentioned earlier, MLPG1 method is selected because when the length scale parameter, l, is not zero, the governing PDEs becomes of fourth order. But when l is zero, we can apply MLPG5 and MLPG2 168

V.S. Khorasani, M. Bayat

Engineering Analysis with Boundary Elements 94 (2018) 159–171

Fig. 12. Non-dimensional deflection 𝑤 on all domain of L-shape plate for n = 1 and l = 0. Fig. 10. Transverse shear stress through thickness for various values of n and l = 17.6 × 10 − 6 .

Fig. 13. Comparison between results of MLPG1, MLPG2, and MLPG5 for n = 1 and l = 0. Fig. 11. The geometry of L-shape plate and distribution of nodes.

the results of MLPG method with the Navier (analytical) solutions and investigate the sources of errors.

Table 1 The CPU time of solving method, MLPG1, for various values of n and l. n

0 0.5 1 2 5 10

CPU time (s) l=0

l≠0

27.68 27.06 28.61 28.96 25.89 31.71

139.09 112.87 156.45 114.30 124.59 135.16

6. Error analysis The results of the MLPG method in this research are compared to results of analytical approach [23] in Table 2. Although it is traditional in the literature that diagrams are drawn for non-dimensionalized variables, it is needed to compare the actual results for error analysis. According to Table 2, it is clear that the exact errors between present solutions and analytical results [23] are very little, the NRMS errors are 169

V.S. Khorasani, M. Bayat

Engineering Analysis with Boundary Elements 94 (2018) 159–171

Table 2 Comparison on center deflection 𝑤( 𝑎2 , 2𝑏 , 0) for various values of l and n. l/h

n

𝑤( 𝑎2 , 2𝑏 , 0) Present

Analytical [23]

0

0 0.5 1 5 10 0 5 10

6.2501 × 10−13 8.9287 × 10−13 1.1363 × 10−12 2.5000 × 10−12 3.4375 × 10−12 1.3677 × 10−13 5.4711 × 10−13 7.5227 × 10−13

8.2732 × 10−12 1.3663 × 10−11 1.9304 × 10−11 3.7731 × 10−11 4.1742 × 10−11 1.2964 × 10−12 5.3704 × 10−12 7.2222 × 10−12

1

Fig. 15. Values of stiffness matrix, K, arrays for a 11-by-11 nodes distribution.

Fig. 14. Normalized rank of stiffness matrix, K, versus number of nodes.

acceptable but the relative errors are not negligible. In this section, two sources of this error is presented. 6.1. The rank of stiffness matrix, K By employing 121 nodes for domain discretization in present problem, stiffness matrix, K, becomes 1331 × 1331. The rank of K is obtained 522, which means that we cannot use regular matrix inverting techniques, therefore for inverting K, the Moore–Penrose method is applied for solving linear equations system. It is founded that by increasing number of nodes in discretizing domain, the rank of K will increase as shown in Fig. 14. Until we reach to a full rank K, accuracy will not increase noticeably because of using pseudoinverse techniques for solving the system of equations. For obtaining a full rank K, it is needed to consider more than 200 × 200 nodes. This nodes distribution require high RAM memory capacity, which could not be granted in this research.

Fig. 16. Comparison of two sample arrays of first submatrix of K, 𝐾𝑖𝑗0101 .

values, leads to a matrix with high condition number and therefore reduction of accuracy for inverting this matrix for solving linear system of equations after numerical integration procedure. 7. Conclusion In this paper, bending behavior of a FG plate based on a general third-order plate theory and the modified couple stress model is analyzed by using the meshless local Petrov–Galerkin method. The moving least-squares approximation is employed to represent the unknown generalized displacements in governing equations which are a fourth-order PDE system. The Gauss–Legendre quadrature method is used for numerical integration of the local weak form. The variation of Young’s modulus of two constituents through the thickness of FG plate is considered by applying a power-law model and microstructure-dependent size effects are captured using a length scale parameter. The MLPG solutions show that for larger values of power-law index, n, larger deflections is obtained due to reduction of plate stiffness in all cases. Also it is clearly

6.2. The condition number of stiffness matrix, K The governing equations of GTPT is a linear PDE system with the constant coefficients. These coefficients which are plate stiffnesses vary from 106 to 10−36 . This wide range of coefficients increases the condition number of K so highly that makes K so sensitive to variation of its arrays. Fig. 15 shows values of K arrays for n = 1 and l = 0. In all cases, arrays of K are near to zero except for some arrays of submatrices 𝐾𝑖𝑗0101 and 𝐾𝑖𝑗0202 . Fig. 16 shows high difference between arrays even in submatrix 𝐾𝑖𝑗0101 of K. This inconsistency and irregular difference between arrays of stiffness matrix, which is consequence of the plate stiffnesses

170

V.S. Khorasani, M. Bayat

Engineering Analysis with Boundary Elements 94 (2018) 159–171

shown that by considering microstructural size effects through modified couple stress theory, stiffness of plate is increased and smaller deflections are obtained, therefore length scale parameter makes plate stiffer. By employing the GTPT, no shear correction factor are needed, and the transverse shear stress are not obtained zero on both surfaces of plate. The MLPG results are compared to the analytical solutions that shows very small exact errors, and sources of relative errors is fully investigated.

[19] Ferreira AJ, Roque C, Jorge R, Fasshauer G, Batra R. Analysis of functionally graded plates by a robust meshless method. Mech Adv Mater Struct 2007;14:577–87. [20] Ferreira A, Batra R, Roque C, Qian L, Jorge R. Natural frequencies of functionally graded plates by a meshless method. Compos Struct 2006;75:593–600. [21] Sheikholeslami S, Saidi A. Vibration analysis of functionally graded rectangular plates resting on elastic foundation using higher-order shear and normal deformable plate theory. Compos Struct 2013;106:350–61. [22] Reddy J, Kim J. A nonlinear modified couple stress-based third-order theory of functionally graded plates. Compos Struct 2012;94:1128–43. [23] Kim J, Reddy J. Analytical solutions for bending, vibration, and buckling of FGM plates using a couple stress-based third-order theory. Compos Struct 2013;103:86–98. [24] Kim J, Reddy J. A general third-order theory of functionally graded plates with modified couple stress effect and the von Kármán nonlinearity: theory and finite element analysis. Acta Mech 2015;226:2973. [25] Reddy J. A general non-linear third-order theory of plates with moderate thickness. Int J Non-Linear Mech 1990;25:677–86. [26] Reddy J. A general nonlinear third-order theory of functionally graded plates. Int J Aerosp Lightweight Struct (IJALS) 2011;1:1–21. [27] Yang F, Chong A, Lam DCC, Tong P. Couple stress based strain gradient theory for elasticity. Int J Solids Struct 2002;39:2731–43. [28] Kolter W. Couple stresses in the theory of elasticity. In: Proceedings of the 1964 koninklijke nederlandse akademie van wetenschappen, 67; 1964. [29] Ma H, Gao X-L, Reddy J. A microstructure-dependent Timoshenko beam model based on a modified couple stress theory. J Mech Phys Solids 2008;56:3379–91. [30] Reddy J. Microstructure-dependent couple stress theories of functionally graded beams. J Mech Phys Solids 2011;59:2382–99. [31] Koizumi M. The Concept of FGM, Ceramic Transactions. Functionally Gradient Materials 1993;34:3–10. [32] Miyamoto Y, Kaysser W, Rabin B, Kawasaki A, Ford RG. Functionally graded materials: Design, processing and applications. Springer Science & Business Media; 2013. [33] Tornabene F. Free vibration analysis of functionally graded conical, cylindrical shell and annular plate structures with a four-parameter power-law distribution. Comput Methods Appl Mech Eng 2009;198:2911–35. [34] Tornabene F, Fantuzzi N, Bacciocchi M, Viola E, Reddy JN. A numerical investigation on the natural frequencies of FGM sandwich shells with variable thickness by the local generalized differential quadrature method. Appl Sci 2017;7:131. [35] Tornabene F, Viola E. Free vibration analysis of functionally graded panels and shells of revolution. Meccanica 2009;44:255–81. [36] Viola E, Tornabene F. Free vibrations of three parameter functionally graded parabolic panels of revolution. Mech Res Commun 2009;36:587–94. [37] Tornabene F, Viola E. Free vibrations of four-parameter functionally graded parabolic panels and shells of revolution. Eur J Mech-A/Solids 2009;28:991–1013. [38] Reddy J. Analysis of functionally graded plates. Int J Numer Methods Eng 2000;47:663–84. [39] Tornabene F, Fantuzzi N, Bacciocchi M, Neves AM, Ferreira AJ. MLSDQ based on RBFs for the free vibrations of laminated composite doubly-curved shells. Compos Part B Eng 2016;99:30–47. [40] Lancaster P, Salkauskas K. Surfaces generated by moving least squares methods. Math Comput 1981;37:141–58. [41] Levin D. The approximation power of moving least-squares. Math Comput Am Math Soc 1998;67:1517–31. [42] Mirzaei D. Analysis of moving least squares approximation revisited. J Comput Appl Math 2015;282:237–50. [43] Li X, Wang Q. Analysis of the inherent instability of the interpolating moving least squares method when using improper polynomial bases. Eng Anal Bound Elem 2016;73:21–34. [44] Abbaszadeh M, Dehghan M. The two-grid interpolating element free Galerkin (TG-IEFG) method for solving Rosenau-regularized long wave (RRLW) equation with error analysis. Appl Anal 2017;97(7):1129–53. [45] Dehghan M, Mirzaei D. The meshless local Petrov–Galerkin (MLPG) method for the generalized two-dimensional non-linear Schrödinger equation. Eng Anal Bound Elem 2008;32:747–56. [46] Mirzaei D, Dehghan M. Meshless local Petrov–Galerkin (MLPG) approximation to the two dimensional sine-Gordon equation. J Comput Appl Math 2010;233:2737–54. [47] Kythe PK, Schäferkotter MR. Handbook of computational methods for integration. CRC Press; 2004.

Acknowledgment Valuable comments and suggestions from Dr. Jinseok Kim (Western Michigan University) and Dr. Davoud Mirzaei (University of Isfahan) are gratefully acknowledged. References [1] Liu G-R, Gu Y-T. An introduction to meshfree methods and their programming. Springer Science & Business Media; 2005. [2] Dehghan M, Abbaszadeh M, Mohebbi A. Numerical solution of system of n-coupled nonlinear schrodinger equations via two variants of the meshless local Petrov–Galerkin (MLPG) method. Comput Model Eng Sci 2014;100:399–444. [3] Dehghan M, Abbaszadeh M. Numerical investigation based on direct meshless local Petrov Galerkin (direct MLPG) method for solving generalized Zakharov system in one and two dimensions and generalized Gross–Pitaevskii equation. Eng Comput 2017;33:983–96. [4] Dehghan M, Salehi R. A meshless local Petrov–Galerkin method for the time-dependent Maxwell equations. J Comput Appl Math 2014;268:93–110. [5] Sladek J, Stanak P, Han Z, Sladek V, Atluri S. Applications of the MLPG method in engineering & sciences: a review. Comput Model Eng Sci 2013;92:423–75. [6] Atluri SN, Zhu T. A new meshless local Petrov–Galerkin (MLPG) approach in computational mechanics. Comput Mech 1998;22:117–27. [7] Atluri SN, Shen SP. The meshless local Petrov–Galerkin (MLPG) method: a simple & less-costly alternative to the finite element and boundary element methods. Comput Model Eng Sci 2002;3:11–51. [8] Atluri SN, Shen S. The meshless local Petrov–Galerkin (MLPG) method. Crest; 2002. [9] Atluri SN. The meshless method (MLPG) for domain & BIE discretizations. Tech Science Press Forsyth; 2004. [10] Qian L, Batra R, Chen L. Elastostatic deformations of a thick plate by using a higher-order shear and normal deformable plate theory and two meshless local Petrov— Galerkin (MLPG) methods. Comput Model Eng Sci 2003;4:161–76. [11] Qian L, Batra R, Chen L. Static and dynamic deformations of thick functionally graded elastic plates by using higher-order shear and normal deformable plate theory and meshless local Petrov–Galerkin method. Compos Part B Eng 2004;35:685–97. [12] Qian L, Batra R. Transient thermoelastic deformations of a thick functionally graded plate. J Therm Stress 2004;27:705–40. [13] Qian L, Batra R. Three-dimensional transient heat conduction in a functionally graded thick plate with a higher-order plate theory and a meshless local Petrov— Galerkin method. Comput Mech 2005;35:214–26. [14] Qian L, Batra R. Design of bidirectional functionally graded plate for optimal natural frequencies. J Sound Vib 2005;280:415–24. [15] Xiao J, Batra R, Gilhooley D, Gillespie J, McCarthy M. Analysis of thick plates by using a higher-order shear and normal deformable plate theory and MLPG method with radial basis functions. Comput Methods Appl Mech Eng 2007;196:979–87. [16] Gilhooley D, Batra R, Xiao J, McCarthy M, Gillespie J. Analysis of thick functionally graded plates by using higher-order shear and normal deformable plate theory and MLPG method with radial basis functions. Compos Struct 2007;80:539–52. [17] Xiao J, Gilhooley D, Batra R, Gillespie J, McCarthy M. Analysis of thick composite laminates using a higher-order shear and normal deformable plate theory (HOSNDPT) and a meshless method. Compos Part B Eng 2008;39:414–27. [18] Ferreira A, Batra R, Roque C, Qian L, Martins P. Static analysis of functionally graded plates using third-order shear deformation theory and a meshless method. Compos Struct 2005;69:449–57.

171