RBF-FD meshless optimization using direct search (GLODS) in the analysis of composite plates

ARTICLE IN PRESS

JID: EABE

[m5GeSdc;November 24, 2017;9:40]

Engineering Analysis with Boundary Elements 000 (2017) 1–10

Contents lists available at ScienceDirect

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

RBF-FD meshless optimization using direct search (GLODS) in the analysis of composite plates C.M.C. Roque a,∗, J.F.A. Madeira b,c a

INEGI, Faculdade de Engenharia da Universidade do Porto, Rua Dr. Roberto Frias, Porto 4200-465, Portugal IDMEC, Instituto Superior Tecnico, Universidade de Lisboa, Lisbon, Portugal c Department of Mathematics, ISEL, IPL, Lisbon, Portugal b

a b s t r a c t An optimization technique, Global and Local Optimization using Direct Search (GLODS) is used to optimize stencil size and shape parameter of the multiquadric function used in the Radial Basis Function-Finite Difference (RBF-FD) meshless numerical method. The mechanical problem to be solved by the meshless method is a composite plate in bending, under sinusoidal load. Results show GLODS is capable of finding values for stencil size and shape parameter that produce excellent numerical solutions when compared with analytical solutions. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction Composite structures such as composite plates are increasingly used in many areas of engineering, including civil and aeronautic engineering. Regarding the numerical simulation of composite structures in general, although the finite element method is still very popular, many authors are now using meshless numerical methods to study the behavior of composite structures. Meshless numerical methods used to study mechanical problems includes the SPH method [1], generalized finite difference method [2], element-free Galerkin (EFG) method [3], h–p cloud method [4,5], moving Kriging interpolation-based meshless, [6], reproducing kernel method [7], amongst many others. In particular, authors used the multiquadric global collocation method to study the static and dynamic behavior of composite beams, plates and shells [8–10]. The method was originally developed by Kansa to solve partial differential equations and uses a global collocation procedure. In its global version, the method can produce severe ill conditioning in matrices. More recently, a local version of the method having the finite difference method as a basis has emerged. The use of local numerical schemes, such as finite differences produces much better conditioned matrices. By combining finite differences and radial basis functions it is possible to obtain a versatile meshless method, the radial basis function-finite difference technique (RBF-FD) [11]. The method shows great potential to solve large engineering problems when compared to traditional global RBF collocation method, since the conditioning of the problem is greatly improved. Due to its excellent results, many authors are increasingly using local versions of the RBF method in diverse areas of physics and engineering [12–21].

∗

Still, the most favorable grid distribution, stencil size and shape parameter in RBF-FD method remains an open problem. A recent comparison between meshless weak and strong formulations for boundary value problems shows how sensitive the method is to the shape parameter, specially when Neumann boundary conditions are used [22]. Both global and local RBF collocation methods using the multiquadric radial basis function depend on a user defined shape parameter, here denoted as c, that has a big influence in the accuracy of solutions. For global collocation, a practical rule consists in choosing a shape parameter proportional to the distance between consecutive nodes, or proportional to the number of nodes in the grid [23,24]. This rule of thumb is used by many researches, producing not optimal, but quite accurate results for many applications, using very simple formulas [25,26]. In cases where these simple formulas do not produce enough accurate results, other methodologies can be used to find a good shape parameter. Most reported work in the literature is focused on global RBF collocation. Some techniques evolved around the idea of cross validation [27–30]. Authors considered multiobjective optimization strategies to choose an optimal shape parameter and grid distribution to solve boundary value problems in plates in bending [31,32]. More recently, a theoretical approach was developed by Luh for the generalized multiquadric and inverse multiquadric interpolation problems [33]. Regarding local collocation, in [34] RBF-FD with RBF polyharmonic spline and polynomials is applied to different grid distributions to solve the Navier–Stokes equations. Tested grids include cartesian, hexagonal, and quasi-uniform. Hexagonal, and quasi-uniform grids produced marginally better results than cartesian. In [35,36] constant and variable shape parameters are found based on analytical formulas derived for the local approximation error of

Corresponding author. E-mail address: [email protected] (C.M.C. Roque).

https://doi.org/10.1016/j.enganabound.2017.10.019 Received 20 July 2017; Received in revised form 12 October 2017; Accepted 14 October 2017 Available online xxx 0955-7997/© 2017 Elsevier Ltd. All rights reserved.

Please cite this article as: C.M.C. Roque, J.F.A. Madeira, Engineering Analysis with Boundary Elements (2017), https://doi.org/10.1016/j. enganabound.2017.10.019

ARTICLE IN PRESS

JID: EABE C.M.C. Roque, J.F.A. Madeira

Engineering Analysis with Boundary Elements 000 (2017) 1–10

RBF-FD formulas. Results show that a variable shape parameter (different for each stencil) produces in general more accurate results. This is in line with previous recommendations for global RBF method [37,38]. Also regarding shape parameter, Sarra applies to RBF-FD an algorithm to restrain the conditioning of matrices to a chosen interval [39]. This is done by choosing an appropriate shape parameter. The method produces excellent results but a good initial guess is needed. In this work, an optimization technique GLODS - Global and Local Optimization using Direct Search is used to optimize grid distribution and find a suitable shape parameter for each stencil in RBF-FD collocation method. Regarding stencil size, a compromise between accuracy and dimension has to be reached. Too small stencils can produce poor solutions but large stencils can sometimes introduce numerical instabilities. Shape parameter defining the radial basis function is a parameter chosen by the user and is of great importance to achieve stable and accurate solutions. In order to find the best stencil size and shape parameter for each grid point, an optimization algorithm (GLODS) is used. Three optimization problems are proposed. In the first, the shape parameter is optimized for each grid point, (which coincides with a center of the radial basis function). In the second problem, the radius defining the number of points in the RBF-FD stencil is optimized, for each stencil. In the third optimization problem both shape parameter and stencil size are chosen by GLODS, making the RBF-FD completely free of user defined parameters. GLODS is suited for global derivative-free constrained optimization [40]. Using direct search of directional type, the algorithm alternates between a search step, where potentially good regions are located, and a poll step where the previously located promising regions are explored. The method is here applied to the bending of plates and can be applied to any other engineering problems, regardless its complexity. As an example of application to a complex geometry, a simple Poisson problem is solved considering an irregular domain.

Fig. 1. Global set 𝜒 and local subset 𝜒 i .

combination of smooth radial basis functions, for N grid points 𝑢(𝑥) ≈ 𝑠(𝑥) =

≅

𝑛 ∑ 𝑖=1

𝐮 = 𝐀𝛂

where 𝐀 = 𝜙(𝑥 − 𝑥𝑖 ). By using the RBF function we guarantee the nonsingularity of matrix A for distinct data centers [41] (For certain radial basis functions this property is not guaranteed and a polynomial has to be added to interpolant (3). Using a global collocation method, we apply a linear operator  to Eq. (3) to obtain: 𝑠 =

𝑛 ∑

𝛼𝜙(𝑥 − 𝑥𝑖 )

(5)

The computational implementation of the global collocation method is straightforward and produces excellent results, if the shape parameter c is well chosen. Unfortunately, this approach produces dense, illconditioned, matrices. In order to improve the conditioning number of matrix A, a local collocation approach was formulated by Tolstykh et al. [42], Cecil et al. [43] and Wright and Fornberg [44], in the so-called RBF-FD technique. In this local RBF-FD collocation approach we compute weights pi such that:

(1)

where wk are weights. In classical finite difference formulas, nodes from i to n are equidistant and weights are computed using polynomial interpolation. This limitation of the spatial distribution of nodes is not desirable for more generic problems. However, using RBFs as an interpolant, a scattered distribution of nodes can be used. Radial basis functions basically depend on the distance between points. This property of RBFs adds flexibility to the method, since it is easy to compute distances between points, even for complex 2D/3D geometries. Consider a set of nodes 𝑥1 , 𝑥2 , … , 𝑥𝑁 ∈ Ω ⊂ ℝ𝑛 . The radial basis functions centered at xj are defined as ( ) 𝜙(𝐱) ≡ 𝜙 ‖𝐱 − 𝐱𝑖 ‖ ∈ ℝ𝑛 , 𝑖 = 1, … , 𝑁 (2)

𝑢(𝑥) =

𝑛𝑖 ∑ 𝑖=1

𝑝𝑖 𝑢(𝑥𝑖 )

(6)

For each node, weights pi are computed on a subset 𝜒𝑖 = [1, … , 𝑛𝑖 ] of the original set of points 𝜒 = [1, … , 𝑛] (see Fig. 1). Subsets 𝜒 i can have an irregular distribution of points and different subsets can have different number of points, increasing the flexibility of the present formulation. In the present paper we usually fix a support distance and consider the nodes within the circle, for each node (also denoted as center). In order to derive the RBF-FD formulation, an interpolant (3) is considered in a Lagrangean form as

where ‖𝐱 − 𝐱𝑖 ‖ is the Euclidian norm. Radial basis functions rely on the euclidean distance between nodes and in some cases on a user-defined shape parameter c. Although many RBF functions could be used, some of the most used RBFs are [37,38]:

𝑠(𝑥) =

𝑛 ∑ 𝑖=1

𝜓𝑖 (𝑥)𝑢(𝑥𝑖 )

(7)

where 𝜓𝑖 (𝑥𝑘 ) = 𝛿𝑖𝑘 , k=1, ..., n. A closed form solution exists for the case (𝜓𝑖 (𝑥𝑘 ) = 𝛿𝑖𝑘 ), and is given by [45]:

1

• Multiquadrics: 𝜙(𝐱) = (‖𝐱 − 𝐱𝑖 ‖ + 𝑐 2 ) 2

(3)

(4)

𝑖=1

𝑤𝑘𝑗,𝑖 𝑢(𝑥𝑖 )

𝛼𝑖 𝜙(𝑥 − 𝑥𝑖 ), 𝑖 = 1, … , 𝑁

or in matrix-vector notation:

The basic formulation of the RBF-FD method is presented. The finite difference method approximates derivatives of a function u(x) at point 𝑥 = 𝑥𝑗 by: 𝑑𝑥𝑘

𝑛 ∑ 𝑖=1

2. The RBF-FD method

𝑑 𝑘 𝑢(𝑥𝑗 )

[m5GeSdc;November 24, 2017;9:40]

1

• Inverse multiquadrics: 𝜙(𝐱) = (‖𝐱 − 𝐱𝑖 ‖ + 𝑐 2 )− 2 2 2 • Gaussians: 𝜙(𝐱) = 𝑒−𝑐 ‖𝐱−𝐱𝑖 ‖ • Thin Plate Splines: 𝜙(𝐱) = ‖𝐱 − 𝐱𝑖 ‖2 log ‖𝐱 − 𝐱𝑖 ‖

𝜓𝑖 (𝑥) =

𝑑𝑒𝑡(𝐴𝑖 (𝑥)) 𝑑𝑒𝑡(𝐀)

(8)

where Ai (x) is obtained by replacing rows i in matrix A (defined in (4)) by vector [ ( ) ( ) ( )] 𝐁 = 𝜙 𝑥 − 𝑥1 𝜙 𝑥 − 𝑥2 … 𝜙 𝑥 − 𝑥𝑛 (9)

In this paper we restrict the analysis to multiquadrics radial basis functions, as interpolants for the solution of partial differential equations. The approximate solution of a PDE can be represented by a linear 2

ARTICLE IN PRESS

JID: EABE C.M.C. Roque, J.F.A. Madeira

[m5GeSdc;November 24, 2017;9:40] Engineering Analysis with Boundary Elements 000 (2017) 1–10

one of the most popular is the third-order theory of Reddy [3,4] that proved to be quite adequate for the analysis of laminated plates and is adopted in this work. The third-order theory of Reddy [50,51] is based on the same assumptions than the classical and first-order plate theories, except that the assumption of straightness and normality of a transverse normal to the middle-surface of the plate after deformation is relaxed by expanding the displacements (u, v, w) as cubic functions of the thickness coordinate. 3.1. Displacements The displacement field for every material point of the plate is defined as ( ) 𝜕𝑤0 (𝑥, 𝑦) 4 3 𝑢(𝑥, 𝑦, 𝑧) = 𝑢0 (𝑥, 𝑦) + 𝑧𝜃𝑥 (𝑥, 𝑦) − 𝑧 𝜃𝑥 (𝑥, 𝑦) + 𝜕𝑥 3ℎ2 ( ) 𝜕𝑤0 (𝑥, 𝑦) 4 3 𝑣(𝑥, 𝑦, 𝑧) = 𝑣0 (𝑥, 𝑦) + 𝑧𝜃𝑦 − 𝑧 𝜃𝑦 (𝑥, 𝑦) + 𝜕𝑦 3ℎ2 𝑤(𝑥, 𝑦, 𝑧) = 𝑤0 (𝑥, 𝑦) (15)

y

0

where u and v are the inplane displacements at any point (x, y, z), u0 and v0 denote the inplane displacement of the point (x, y, 0) on the midplane, w is the transverse displacement (here considered to be constant across the thickness direction), 𝜃 x and 𝜃 y are the rotations of the normals to the midplane about the y and x axes, respectively.

1

x

3.2. Deformations of the plate

Fig. 2. Optimization problem 1. Optimized shape parameter for each grid point. 𝑐𝑚𝑖𝑛 = 1.1; 𝑐𝑚𝑎𝑥 = 9.1.

The deformations–displacement relationships are defined by: { } { } 𝜕𝑢 𝜕𝑣 𝜕𝑢 𝜕𝑣 𝜕𝑢 𝜕𝑤 𝜕𝑣 𝜕𝑤 , , + , + , + 𝜖𝑥𝑥 , 𝜖𝑦𝑦 , 𝛾𝑥𝑦 , 𝛾𝑥𝑧 , 𝛾𝑦𝑧 = 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕𝑥 𝜕𝑧 𝜕𝑥 𝜕𝑧 𝜕𝑦

The approximation of derivative 𝑢(𝑥𝑗 ) is now performed by applying the linear differential operator on the interpolant in (7) as 𝑢(𝑥𝑗 ) ≈ 𝑠(𝑥𝑗 ) =

𝑛 ∑ 𝑖=1

𝜓𝑖 (𝑥𝑗 )𝑢(𝑥𝑖 )

Deformations can be depicted as { } { } { } { (0) (0) (0) } (1) (1) (1) (3) (3) (3) , 𝜖𝑦𝑦 , 𝛾𝑥𝑦 + 𝑧 𝜖𝑥𝑥 , 𝜖𝑦𝑦 , 𝛾𝑥𝑦 + 𝑧3 𝜖𝑥𝑥 , 𝜖𝑦𝑦 , 𝛾𝑥𝑦 𝜖𝑥𝑥 , 𝜖𝑦𝑦 , 𝛾𝑥𝑦 = 𝜖𝑥𝑥

(10)

(17)

{ } { } { (0) (0) } (2) (2) , 𝛾𝑦𝑧 + 𝑧2 𝛾𝑥𝑧 , 𝛾𝑦𝑧 𝛾𝑥𝑧 , 𝛾𝑦𝑧 = 𝛾𝑥𝑧

From Eqs. (10) and (6) it is possible to conclude that 𝑝𝑖 = 𝜓𝑖 (𝑥𝑗 ). Although is possible to find a closed-form solution for 𝜓𝑖 (𝑥𝑗 ) (as in the interpolation problem), in practice weights are computed by solving the linear system: [ ( )] 𝚽𝐩 =  𝐁(𝑥𝑖 ) (11)

(18)

where the in-plane deformations are given by } { } { 𝜕𝑢 𝜕𝑣 𝜕𝑢 𝜕𝑣 0 (0) (0) (0) 𝜖𝑥𝑥 , 𝜖𝑦𝑦 , 𝛾𝑥𝑦 = , 0, 0 + 0 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕𝑥

where B(x) was computed in Eq. (9), 𝚽 is a matrix containing the RBF functions, and p is a matrix of weights. For the static analysis of plates in bending, a boundary value problem with domain Ω and boundary 𝛿Ω is considered. The PDE problem is defined as: { 𝐮(𝑥) = 𝐟 (𝑥), 𝑥 ∈ Ω (12) 𝐮(𝑥) = 𝐠(𝑥), 𝑥 ∈ 𝛿Ω

(19)

the flexural deformations are expressed as } { } { 𝜕𝜃 𝜕𝜃𝑦 𝜕𝜃 𝜕𝜃𝑦 𝑥 (1) (1) (1) 𝜖𝑥𝑥 , 𝜖𝑦𝑦 , 𝛾𝑥𝑦 = , , 𝑥 + 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕𝑥

(20)

the higher-order deformations are obtained as { } { } 𝜕𝜃𝑥 𝜕 2 𝑤0 𝜕𝜃𝑦 𝜕 2 𝑤0 𝜕𝜃𝑥 𝜕𝜃𝑦 𝜕 2 𝑤0 (3) (3) (3) 𝜖𝑥𝑥 , 𝜖𝑦𝑦 , 𝛾𝑥𝑦 = −𝑐1 + , + , + +2 𝜕𝑥 𝜕𝑥 𝜕 𝑥𝜕 𝑦 𝜕𝑥2 𝜕𝑦 𝜕𝑦2 𝜕𝑦 (21) The transverse shear deformations are expressed as } { } { 𝜕𝑤 𝜕𝑤0 0 (0) (0) 𝛾𝑥𝑧 , 𝛾𝑦𝑧 = + 𝜃𝑥 , + 𝜃𝑦 (22) 𝜕𝑥 𝜕𝑦

Here,  is a linear differential operator, and  is a linear differential operator related to boundary conditions. Weights p are computed for each node i by solving the linear system of equations [ ( )] 𝚽𝐩 =  𝐵(𝑥𝑖 ) (13)

while the higher-order transverse shear deformations are expressed as { } { } 𝜕𝑤0 𝜕𝑤0 (2) (2) 𝛾𝑥𝑧 , 𝛾𝑦𝑧 = −𝑐2 + 𝜃𝑥 , + 𝜃𝑦 (23) 𝜕𝑥 𝜕𝑦

where the operator  is defined by  =  in case i ∈ Ω or by  =  in case i ∈ 𝜕Ω. Solution u can now be found by assembling the p weights in global matrix M and then solving the linear algebraic collocation system: 𝐌𝐮 = 𝐅

(16)

and 𝑐1 = 3ℎ4 2 , 𝑐2 = 3𝑐1 . Due to the cubic evolution in the z-coordinate of the in-plane displacements and the quadratic evolution of the transverse deformations, this theory does not consider the shear correction factor, as in the case of first-order shear deformation theories.

(14)

[ 𝑝Ω ] [ 𝑢Ω ] where 𝐌 = 𝛿Ω is the matrix of weights; 𝐮 = 𝛿Ω is the vector of 𝑝 𝑢 [ 𝐟 (𝑥) ] global unknowns and 𝐅 = is the vector of external forces. In the 𝐠(𝑥) case of plates in bending, f(x) represent applied forces on the domain and g(x) the (essential and natural) boundary conditions.

3.3. Stresses Neglecting 𝜎 z for each layer (k), the stress-deformation relations in the layer local coordinate system can be expressed as ⎧ 𝜎1 ⎫ ⎪𝜎 ⎪ ⎪ 2⎪ ⎨𝜏12 ⎬ ⎪𝜏23 ⎪ ⎪ ⎪ ⎩𝜏31 ⎭

3. Third-order shear deformation theory The study of deformation theories to model composite laminates is today a very active field of research [46–49]. Amongst many theories, 3

(𝑘)

⎡𝑄11 ⎢𝑄 ⎢ 12 =⎢ 0 ⎢ 0 ⎢ ⎣ 0

𝑄12 𝑄22 0 0 0

0 0 𝑄33 0 0

0 0 0 𝑄44 0

(𝑘)

0 ⎤ 0 ⎥ ⎥ 0 ⎥ 0 ⎥ ⎥ 𝑄55 ⎦

⎧ 𝜀1 ⎫ ⎪𝜀 ⎪ ⎪ 2⎪ ⎨𝛾12 ⎬ ⎪𝛾23 ⎪ ⎪ ⎪ ⎩𝛾31 ⎭

(𝑘)

(24)

ARTICLE IN PRESS

JID: EABE

[m5GeSdc;November 24, 2017;9:40]

C.M.C. Roque, J.F.A. Madeira

Engineering Analysis with Boundary Elements 000 (2017) 1–10

where subscripts 1 and 2 are, respectively, the fiber and the normal to fiber inplane directions, 3 is the direction normal to the plate, and the reduced stiffness components, Qij are given by 𝑄(11𝑘) =

(𝑘) 𝐸1 (𝑘) (𝑘) 1−𝜈12 𝜈21

𝑄(22𝑘) =

(𝑘) 𝑄(33𝑘) = 𝐺12 (𝑘) (𝑘) 𝜈21 = 𝜈12

(𝑘) 𝐸2 (𝑘) (𝑘) 1−𝜈12 𝜈21

(𝑘) 𝑄(44𝑘) = 𝐺23

and 𝛿𝑉 = −

(𝑘) (𝑘) 𝑄(12𝑘) = 𝜈12 𝑄22

(𝑘)

(𝑘) 𝑄(55𝑘) = 𝐺13

(𝑘) 𝐸2 (𝑘) 𝐸1

⎡𝑄11 ⎢𝑄 ⎢ 12 = ⎢𝑄16 ⎢ 0 ⎢ ⎣ 0

𝑄12 𝑄22 𝑄26 0 0

𝑄16 𝑄26 𝑄66 0 0

0 0 0 𝑄44 𝑄45

(𝑘)

0 ⎤ 0 ⎥ ⎥ 0 ⎥ 𝑄45 ⎥ ⎥ 𝑄55 ⎦

⎧𝜀𝑥𝑥 ⎫ ⎪𝜀 ⎪ ⎪ 𝑦𝑦 ⎪ ⎨ 𝛾𝑥𝑦 ⎬ ⎪ 𝛾𝑦𝑧 ⎪ ⎪ ⎪ ⎩ 𝛾𝑧𝑥 ⎭

𝑞 𝛿𝑤0 𝑑 𝑥𝑑 𝑦

(27)

where Ω0 denotes the midplane of the laminate, q is the external distributed load. The stress resultants are defined as

(𝑘) (𝑘) (𝑘) (𝑘) in which 𝐸1(𝑘) , 𝐸2(𝑘) , 𝜈12 , 𝐺12 , 𝐺23 and 𝐺13 are independent material properties for each ayer. By performing adequate coordinate transformation, the stress–strain relations in the global x–y–z coordinate system can be obtained for each layer as

⎧𝜎𝑥𝑥 ⎫ ⎪𝜎 ⎪ ⎪ 𝑦𝑦 ⎪ ⎨ 𝜏𝑥𝑦 ⎬ ⎪ 𝜏𝑦𝑧 ⎪ ⎪ ⎪ ⎩ 𝜏𝑧𝑥 ⎭

∫Ω0

⎧𝑁 ⎫ ⎧1⎫ ⎧1⎫ 𝑁𝑙 ℎ∕2 𝑧𝑘+1 ∑ ⎪ 𝛼𝛽 ⎪ ⎪ ⎪ (𝑘) ⎪ ⎪ = 𝜎 𝑑𝑧 = 𝜎 𝑀 𝑧 𝑧 𝑑𝑧 ⎨ 𝛼𝛽 ⎬ ∫ 𝛼𝛽 ⎨ ⎬ 𝛼𝛽 ⎨ ⎬ ∫ −ℎ∕2 ⎪ 𝑃𝛼𝛽 ⎪ ⎪𝑧 3 ⎪ ⎪𝑧3 ⎪ 𝑘=1 𝑧𝑘 ⎩ ⎭ ⎩ ⎭ ⎩ ⎭

(28)

{ } { } { } 𝑁𝑙 ℎ∕2 𝑧𝑘+1 ∑ 𝑄𝛼 1 1 (𝑘) = 𝜎𝛼𝑧 2 𝑑𝑧 = 𝜎𝛼𝑧 2 𝑑𝑧 ∫−ℎ∕2 ∫𝑧𝑘 𝑅𝛼 𝑧 𝑧 𝑘=1

(29)

where 𝛼, 𝛽 take the symbols x, y, Nl represents the number of layers, and 𝑧𝑘 , 𝑧𝑘+1 the bottom and top z-coordinates of each layer k. Substituting for 𝛿U, 𝛿V, into the virtual work statement, noting that the virtual strains can be expressed in terms of the generalized displacements, integrating by parts to relieve from any derivatives of the generalized displacements and using the fundamental lemma of the calculus of variations, we obtain the equations of motion in terms of the stress resultants [50]:

(𝑘)

(25)

The third-order theory of Reddy satisfies zero transverse shear stresses on the bounding planes.

𝜕𝑁𝑥𝑥 𝜕𝑁𝑥𝑦 + =0 𝜕𝑥 𝜕𝑦

(30)

3.4. Equations of motion 𝜕𝑁𝑥𝑦

The equations of motion of the third-order theory are derived from the principle of virtual displacements [50,51]. The virtual strain energy (𝛿U) and the virtual work done by applied forces (𝛿V) are given by { 𝛿𝑈 =

∫Ω0

ℎ∕2

∫−ℎ∕2

[

𝜕𝑥

∫Ω0

𝜕𝑁𝑦𝑦 𝜕𝑦

=0

𝜕 𝑄𝑥 𝜕 𝑄𝑦 4 + + 𝜕𝑥 𝜕𝑦 3ℎ2

( (0) ) (1) (3) 𝜎𝑥𝑥 𝛿𝜖𝑥𝑥 + 𝑧𝛿𝜖𝑥𝑥 − 𝑐1 𝑧3 𝛿𝜖𝑥𝑥

( ) (0) (1) (3) +𝜎𝑦𝑦 𝛿𝜖𝑦𝑦 + 𝑧𝛿𝜖𝑦𝑦 − 𝑐1 𝑧3 𝛿𝜖𝑦𝑦 ( ) ( (0) ) (0) (1) (3) (2) +𝜏𝑥𝑦 𝛿𝛾𝑥𝑦 + 𝑧𝛿𝛾𝑥𝑦 − 𝑐1 𝑧3 𝛿𝛾𝑥𝑦 + 𝑧2 𝛿𝛾𝑥𝑧 + 𝜏𝑥𝑧 𝛿𝛾𝑥𝑧 ( )] } (0) (2) + 𝑧2 𝛿𝛾𝑦𝑧 𝑑𝑧 𝑑𝑥 𝑑𝑦 +𝜏𝑦𝑧 𝛿𝛾𝑦𝑧

=

+

(31)

(

𝜕 2 𝑃𝑥𝑥 𝜕𝑥2

+2

𝜕 2 𝑃𝑥𝑦 𝜕 𝑥𝜕 𝑦

+

𝜕 2 𝑃𝑦𝑦

)

𝜕𝑦2

+𝑞 =0

𝜕 𝑀 𝑥𝑥 𝜕 𝑀 𝑥𝑦 + − 𝑄𝑥 = 0 𝜕𝑥 𝜕𝑦

𝜕 𝑀 𝑥𝑦

( (0) (1) (3) (0) (1) 𝑁𝑥𝑥 𝛿𝜖𝑥𝑥 + 𝑀𝑥𝑥 𝛿𝜖𝑥𝑥 − 𝑐1 𝑃𝑥𝑥 𝛿𝜖𝑥𝑥 + 𝑁𝑦𝑦 𝛿𝜖𝑦𝑦 + 𝑀𝑦𝑦 𝛿𝜖𝑦𝑦

𝜕𝑥

+

𝜕 𝑀 𝑦𝑦 𝜕𝑦

(32)

(33)

− 𝑄𝑦 = 0

(34)

with

(3) (0) −𝑐1 𝑃𝑦𝑦 𝛿𝜖𝑦𝑦 + 𝑁𝑥𝑦 𝛿𝛾𝑥𝑦 (1) (3) (0) (2) (0) +𝑀𝑥𝑦 𝛿𝛾𝑥𝑦 − 𝑐1 𝑃𝑥𝑦 𝛿𝛾𝑥𝑦 + 𝑄𝑥 𝛿𝛾𝑥𝑧 + 𝑅𝑥 𝛿𝛾𝑥𝑧 + 𝑄𝑦 𝛿𝛾𝑦𝑧 ) (2) +𝑅𝑦 𝛿𝛾𝑦𝑧 𝑑𝑥 𝑑𝑦

𝑀 𝛼𝛽 = 𝑀𝛼𝛽 −

4 𝑃𝛼𝛽 ; 3ℎ2

𝑄𝛼 = 𝑄 𝛼 −

4 𝑅𝛼 ; ℎ2

(35)

(26) The stress resultants are related with deformations in the material reference axes by:

⎧⎧ ⎫⎫ ⎪ ⎪𝑁 1 ⎪ ⎪ ⎡ 𝐴12 𝐴16 ⎤ ⎪ ⎨𝑁2 ⎬ ⎪ ⎢⎡ 𝐴11 ⎪ ⎪ ⎪ ⎪ ⎢⎢ 𝐴22 𝐴26 ⎥ ⎥ ⎪ ⎩𝑁6 ⎭ ⎪ ⎢⎢⎣𝑠𝑦𝑚𝑚. 𝐴 66 ⎦ ⎪ ⎪ ⎢ ⎪ ⎪ ⎪⎧𝑀 ⎫⎪ ⎢⎢ 1 ⎪⎪ ⎪⎪ ⎢ ⎨⎨𝑀 2 ⎬⎬ = ⎢ ⎪⎪𝑀 6 ⎪⎪ ⎢ ⎪⎩ ⎭⎪ ⎢ ⎪ ⎪ ⎪ ⎧ ⎫ ⎪ ⎢⎢ ⎪ ⎪𝑃 1 ⎪ ⎪ ⎢ 𝑠𝑦𝑚𝑚𝑒𝑡𝑟𝑖𝑐 ⎪ ⎨𝑃 ⎬ ⎪ ⎢ ⎪ ⎪ 2⎪ ⎪ ⎣ ⎪ ⎩𝑃 6 ⎭ ⎪ ⎩ ⎭

4

⎡ 𝐵11 ⎢ ⎢ ⎣𝑠𝑦𝑚𝑚.

𝐵12 𝐵22

𝐵16 ⎤ 𝐵26 ⎥ ⎥ 𝐵66 ⎦

⎡ 𝐸11 ⎢ ⎢ ⎣𝑠𝑦𝑚𝑚.

𝐸12 𝐸22

⎡ 𝐷11 ⎢ ⎢ ⎣𝑠𝑦𝑚𝑚.

𝐷12 𝐷22

𝐷16 ⎤ 𝐷26 ⎥ ⎥ 𝐷66 ⎦

⎡ 𝐹11 ⎢ ⎢ ⎣𝑠𝑦𝑚𝑚.

𝐹12 𝐹22

⎡ 𝐻11 ⎢ ⎢ ⎣𝑠𝑦𝑚𝑚.

𝐻12 𝐻22

⎧⎧ 0 ⎫⎫ ⎪ 𝜖 ⎪ 𝐸16 ⎤ ⎤⎪⎪ 10 ⎪⎪ ⎨𝜖 ⎬ ⎥ 𝐸26 ⎥ ⎥⎪⎪ 20 ⎪⎪ ⎥ ⎪⎩𝜖6 ⎭⎪ 𝐸66 ⎦ ⎥⎪ ⎪ ⎥⎪ ⎥⎪⎧ 0 ⎫⎪ 𝐹16 ⎤ ⎥ 𝛾1 ⎪ ⎪⎪ ⎪⎪ 𝐹26 ⎥ ⎥⎨⎨𝛾20 ⎬⎬ ⎥ ⎥⎪⎪ 0 ⎪⎪ 𝐹66 ⎦ ⎥ ⎩𝛾6 ⎭ ⎪ ⎥⎪ ⎪ ⎥⎪ 𝐻16 ⎤⎥⎪⎧ 2 ⎫⎪ ⎪ 𝛾 ⎪ 𝐻26 ⎥⎥⎪⎪ 12 ⎪⎪ ⎥⎥ ⎨𝛾2 ⎬ 𝐻66 ⎦⎦⎪⎪ 2 ⎪⎪ ⎪⎩𝛾6 ⎭⎪ ⎩ ⎭

(36)

ARTICLE IN PRESS

JID: EABE C.M.C. Roque, J.F.A. Madeira

[m5GeSdc;November 24, 2017;9:40] Engineering Analysis with Boundary Elements 000 (2017) 1–10

1

0.9

0.8

0.7

d

0.6

0.5

i

0.4

0.3

0.2

0.1

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Fig. 3. Example of radial distance d defining stencil size. Points inside circle with radius d define stencil with center at point i.

{ } [ ] 𝑄2 𝐴45 ⎧ ⎫ ⎡ 𝐴44 ⎪ ⎪ ⎢ 𝐴 𝑄1 𝐴55 ⎪ ⎪ ⎢ 45 ⎨ { } ⎬=⎢ ⎪ ⎪ ⎢ 𝑅2 ⎪ ⎪ ⎢ 𝑠𝑦𝑚𝑚. 𝑅 ⎩ ⎭ ⎣ 1

[

] { } 𝐷45 ⎤⎧ 𝜖40 ⎫ 𝐷55 ⎥⎪ 𝜖50 ⎪ ⎥⎪ ⎪ ] ⎥⎨{ 2 }⎬ ⎥ ⎪ 𝐹45 𝛾 ⎪ ⎥⎪ 4 ⎪ 𝐹55 ⎦⎩ 𝛾52 ⎭

𝐷44 𝐷45

[ 𝐹45 𝐹45

(37)

where Aij , Bij , etc., are the stiffness components of the plate, obtained as (

y

𝑁

) ∑𝑙 𝐴𝑖𝑗 , 𝐵𝑖𝑗 , 𝐷𝑖𝑗 , 𝐸𝑖𝑗 , 𝐹𝑖𝑗 , 𝐻𝑖𝑗 =

𝑘=1

𝑧𝑘+1

∫𝑧𝑘

𝑁

( ) ∑𝑙 (𝑖, 𝑗 = 1, 2, 6), 𝐴𝑖𝑗 , 𝐷𝑖𝑗 , 𝐹𝑖𝑗 =

𝑘=1

( ) 𝑄̄ (𝑖𝑗𝑘) 1, 𝑧, 𝑧2 , 𝑧3 , 𝑧4 , 𝑧6 𝑑𝑧

𝑧𝑘+1

∫𝑧𝑘

( ) 𝑄̄ (𝑖𝑗𝑘) 1, 𝑧2 , 𝑧4 𝑑𝑧

(𝑖, 𝑗 = 4, 5), (38)

being 𝑄̄ (𝑖𝑗𝑘) the transformed elastic constants of each layer [50,51]. For boundary conditions, only simply supported plates are considered. As an example a simply-supported condition in 𝑥 = 𝑎 edge imposes five boundary conditions as follows: 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

𝑣0 = 0,

1

x

𝑤0 = 0,

𝜃𝑦 = 0 ,

𝑀 𝑥 = 0,

𝑁𝑥 = 0

(39)

Considering simply supported boundary conditions, a simple analytical solution is possible to be found, using the Navier method [52].

Fig. 4. Optimization problem 2. Optimized distance for each stencil. 𝑑𝑚𝑖𝑛 = 0.1; 𝑑𝑚𝑎𝑥 = 0.8.

5

ARTICLE IN PRESS

JID: EABE C.M.C. Roque, J.F.A. Madeira

[m5GeSdc;November 24, 2017;9:40] Engineering Analysis with Boundary Elements 000 (2017) 1–10

y

y

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0

1

0.1

0.2

0.3

0.5

0.6

0.7

0.8

0.9

1

Fig. 7. Optimization problem 3. Optimized distance for each stencil. 𝑑𝑚𝑖𝑛 = 0.1; 𝑑𝑚𝑎𝑥 = 0.8.

Fig. 5. Optimization problem 2. Correspondent shape parameter for each grid point, √ computed by 𝑐 = 2∕ 𝑛𝑖 .

where 𝑓 ∶ Ω ⊆ 𝐼 𝑅𝑛 → 𝐼𝑅 ∪ {+∞} represents a real-extended value function and Ω⊆IRn a compact set, defining the problem feasible region. Locating and identifying points as global minimizers is, in general, a hard and time-consuming task. Difficulties increase in the impossibility of using the derivatives of the functions defining the problem. GLODS [40] is an algorithm for single optimization, suited for bound constrained, derivative-free and global optimization. Using direct search of directional type, this optimization method alternates between a search step, where potentially good regions are located, and a poll step where the previously located regions are explored. This exploration is made through the launching of several pattern search methods, one in each of the regions of interest. Differently from a multi-start strategy, the several pattern search methods will merge when sufficiently close to each other. The goal of GLODS algorithm is to end with as many active pattern searches as the number of local minimizers, which would allow to easily locating the possible global extreme value. Similarly to other derivative-free optimization algorithms, it makes use of the extreme barrier function f defined as:

y

𝑓Ω (𝐱) =

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Fig. 6. Optimization problem 3. Optimized shape parameter for each grid point. 𝑐𝑚𝑖𝑛 = 1; 𝑐𝑚𝑎𝑥 = 10.

4. GLODS (Global and local optimization using direct search) method A constrained nonlinear optimization problem can be mathematically formulated as: Find n design variables: ( )𝑇 𝐱 = 𝑥1 , 𝑥2 , … , 𝑥𝑛 ∈ Ω ⊆ 𝐼 𝑅𝑛 (40)

𝑖𝑓 𝐱 ∈ Ω 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒

(42)

5. Numerical examples Two different problems are presented. The first problem is a composite plate in bending. The second is a Poisson problem with an irregular amoeba-like geometry. Most results are presented graphically, where circles with variable diameters represent the values of shape parameter or stencil distance for each point in the grid.

which minimize: min 𝑓 (𝐱)

{ 𝑓 (𝐱) +∞

Unfeasible points will not be evaluated, being the corresponding objective function value set equal to +∞. This approach allows dealing with black-box type constraints, where only a yes/no type of answer is returned. Several details are omitted in the present text and the reader is referred to Custodio et al. [40] for a complete description: the algorithmic structure considered, the convergence analysis and report numerical results.

1

x

𝐱∈ Ω

0.4

x

x

(41) 6

ARTICLE IN PRESS

JID: EABE C.M.C. Roque, J.F.A. Madeira

Engineering Analysis with Boundary Elements 000 (2017) 1–10

5.1. Composite plate

5.1.2. Optimization problem 2. Optimization of stencil size In the second optimization problem, the shape parameter is defined through a typical formula related to the number of points in each stencil √ ni , 𝑐 = 2∕ 𝑛𝑖 . The distance defining the number of points in the stencil is now optimized (see Fig. 3). Results corresponding to the best solution are presented in Figs. 4 and 5. Fig. 4 represents the optimized distance for each stencil. Distances range from 0.1 to 0.8, corresponding to stencils with 3 to 17 points, respectively. The error for the best solution in optimization problem 2 is error = 0.0160. Since the optimization produced a smooth variation in distances and therefore in the number of points for each stencil, the shape parameter also presents a smooth distribution in boundaries and domain, as shown in Fig. 5.

A square composite plate of side 𝑎 = 1 and thickness h is considered. The plate is composed of four equally thick layers with lay-up (0/90/90/0). The plate is subjected to sinusoidal load of the form ( ) ( 𝜋𝑦 ) 𝜋𝑥 𝑝𝑧 = 𝑃 sin sin 𝑎 𝑎 with the origin of the coordinate system being located at the lower left corner on the midplane. The material properties are as follows 𝐸1 = 25.0𝐸2 ;

𝐺12 = 𝐺13 = 0.5𝐸2 ;

𝐺23 = 0.2𝐸2 ;

𝜈12 = 0.25;

Since the plate has has symmetric lay-up, the formulation is simplified in order to consider only variables w, 𝜙x and 𝜙y . Three different optimization problems are defined. In the first problem, only the shape parameter is optimized. In the second problem, only the stencil size is optimized. Finally, in the third optimization problem both shape parameter and stencil size are considered to be design variables. In general radial symmetry is considered, except in corner points of the plate, with coordinates (𝑥, 𝑦) = {(0, 0), (0, 1), (1, 0), (1, 1)}. Error is computed for each variable w, 𝜙x and 𝜙y as: √ √ 𝑁 √1 ∑ | analytical − numerical | | | (43) error(𝑤,𝜙𝑥 ,𝜙𝑦 ) = √ | 𝑁 1 || max(numerical) |

5.1.3. Optimization problem 3. Optimization of stencil size and shape parameter As a final example, both distance defining stencil size and shape parameter for each grid point are optimized. Given the dimension of the optimization problem, the shape parameter interval was reduced to integer. The best solution has an error of 0.0023. The vertical deformation w, rotations 𝜙x and 𝜙y are plotted in Fig. 8 and compared to analytical solutions. Fig. 6 shows the optimized shape parameter distribution, ranging from 1 to 10, 𝑐 = 3 the most representative value. The optimal distance for each stencil is represented in Fig. 7. Optimized distances vary from 0.1 to 0.8, being 𝑑 = 0.3 for the majority of stencils. Although some punctual variations exists, both shape parameter and stencil size (radial distance) remains roughly constant through the domain.

The final error used in the cost function is given by: 𝑓 (𝐱) ≡ cost = max(error). The analytical solution used to compute the error in (43) is given by the Navier procedure. The analytical solution for displacements for simply supported boundary conditions are satisfied by the following expansions of the displacements and applied transverse load [52]: 𝑤0 (𝑥, 𝑦) =

∞ ∑ ∞ ∑ 𝑛=1 𝑚=1

𝜙𝑥 (𝑥, 𝑦) =

𝜙𝑦 (𝑥, 𝑦) =

𝑞(𝑥, 𝑦) =

∞ ∑ ∞ ∑ 𝑛=1 𝑚=1 ∞ ∑ ∞ ∑ 𝑛=1 𝑚=1

∞ ∑ ∞ ∑ 𝑛=1 𝑚=1

𝑊𝑛𝑚 sin(𝜂𝑥) sin(𝜁 𝑦)

(44)

Φ𝑥𝑛𝑚 cos(𝜂𝑥) sin(𝜁 𝑦)

(45)

Φ𝑦𝑛𝑚 sin(𝜂𝑥) cos(𝜁 𝑦)

(46)

𝑄𝑛𝑚 sin(𝜂𝑥) sin(𝜁 𝑦)

5.2. Poisson problem on an irregular geometry A Poisson problem in a amoeba-like geometry is considered. Parametric equations for the boundary are given by 𝑥 = 𝑟(𝑡) cos(𝑡), 𝑦 = 𝑟(𝑡) sin(𝑡) where 𝑟(𝑡) = 𝑒sin(𝑡) sin2 (2𝑡) + 𝑒cos(𝑡) cos2 (2𝑡), t ∈ [0, 2𝜋] [53]. 46 points are distributed on the boundary 𝜕Ω and a Halton distribution is used to distribute points in the interior domain Ω. Interior points are distributed on a square domain and only points inside the boundary are retained, resulting in 21 interior points, as shown in Fig. 9. Both shape parameter and distance for each stencil are optimized. The proposed Poisson problem is chosen in order to have a known solution, u:

(47)

∇2 𝑢 = 13𝑒(2𝑥+3𝑦) on Ω

with 𝑚𝜋 𝑛𝜋 𝜂= ;𝜁 = ; 𝑙1 𝑙2

𝑄𝑛𝑚

[m5GeSdc;November 24, 2017;9:40]

4 = 𝑙1 𝑙2 ∫0

𝑙1

𝑙2

∫0

𝑛𝜋𝑦 𝑚𝜋𝑥 𝑞(𝑥, 𝑦) sin sin 𝑑𝑥 𝑑𝑦 𝑙1 𝑙2

𝑢 = 𝑒(2𝑥+3𝑦) on 𝜕Ω

(49)

denoting u as the exact solution and uh as the numerical approximate solution, the cost function is defined as:

(48) Substituting Eqs. (44)–(47) in (30)–(34) and solving the resulting equations for Δ = (𝑊𝑚𝑛 , Φ𝑥𝑚𝑛 , Φ𝑦𝑚𝑛 ) it is possible to compute analytical solutions for the problem.

error = max(||𝑢 − 𝑢ℎ ||);

(50)

Stencil distances are defined as 0.5 < d < 4.5 with 0.5 step size and shape parameter 1 × 10−3 < 𝑐 < 10 with step size = 1. After 2000 iterations, the optimization found 12 solutions with cost= 3.08. Figs. 10 and 11 show results for optimized shape parameter and distances, respectively. As previously, circles at each grid point represent the obtained optimized values, being the diameter of the circle proportional to the corresponding value of shape parameter in Fig. 10 and stencil diameter in Fig. 11. Optimized shape parameter range form 𝑐 = 1 × 10−3 to 𝑐 = 10 and stencil distance range form 𝑑 = 0.5 to 𝑑 = 4.5.

5.1.1. Optimization problem 1. Optimization of shape parameter, c In the first optimization problem the shape parameter c is optimized for all central points of stencils. Possible values for shape parameter are 𝑐 = [0.1 ∶ 0.1 ∶ 10]. Stencil size is here defined by the radial distance to the central point of the stencil and is kept constant, being 𝑑 = 0.2. By choosing a stencil with a constant radius, the number of points can vary for each stencil. For that reason, a common choice for the shape √ parameter and that works for this type of problems could be 𝑐 = 2∕ 𝑛𝑖 , being ni the number of points in each stencil. The error for this shape parameter is 0.0989. The best solution found by optimization algorithm GLODS produced an error of 0.0699. The corresponding optimized shape parameters are represented in Fig. 2. The smaller and larger circles correspond to 𝑐 = 1.1 and 𝑐 = 9.1, respectively. As can be easily observed, the best solution uses a variable shape parameter along the domain.

6. Comments A compilation of errors for the best solutions for the composite plate in bending problem is presented in Table 1. Errors are smallest when both stencil size and shape parameter values were both optimized (optimization problem 3). In addition, optimization problem 3 was the fastest 7

ARTICLE IN PRESS

JID: EABE C.M.C. Roque, J.F.A. Madeira

[m5GeSdc;November 24, 2017;9:40] Engineering Analysis with Boundary Elements 000 (2017) 1–10

RBF numerical solution Analytical solution

8

20

10

7

15

6

10

5

5

w

φ 4

φ 5

x

y

0

0 3

-5

2

-10

-5

1

-15

0 -1

-10

-20

1

1

1

1

1

0.8 0.5

0.5

0.6 0.4

y

0.2 0

0

1

0.8 0.5

0.6 0.4

y x

0.8

0.2 0

x

0

0.6 0.4

y

0.2 0

x

0

Fig. 8. Optimization problem 3. Best solution.

2

2

1.5

1.5

1

1

0.5

0.5

0

0

-0.5

-0.5

-1

-1

-1.5 -2

-1

0

1

2

-1.5 -2

3

-1

0

1

2

3

Fig. 9. Amoeba-like irregular geometry with a total of 67 grid points.

Fig. 10. Optimization problem 3. Optimized shape parameter for each stencil. 𝑐𝑚𝑖𝑛 = 1 × 10−3 ; 𝑐𝑚𝑎𝑥 = 10.

to produce acceptable solutions. Both optimization problems 1 and 2 converged to solutions that use variable shape parameters whereas optimization problem 3 seems to converge to solutions using a more homogeneous shape parameter and stencil size. In general, the optimization technique GLODS found values for shape parameter and stencil size that produced very small errors when comparing RBF-FD numerical solutions to analytical solutions. When using

an irregular geometry with irregular stencils the optimization algorithm GLODS is able to choose an adequate shape parameter and size for each stencil. We highlight that GLODS is a global optimization method and therefore CPU times can be quit high, depending on the problem and stopping criteria. In the present case the stopping criteria was a predefined number of iterations and therefore each problem took a few hours to complete. 8

ARTICLE IN PRESS

JID: EABE C.M.C. Roque, J.F.A. Madeira

Engineering Analysis with Boundary Elements 000 (2017) 1–10

2

[10] Ferreira AJM, Roque CMC, Jorge RMN. Modelling cross-ply laminated elastic shells by a higher-order theory and multiquadrics. Comput Struct 2006;84(19–20):1288–99. [11] Shan YY, Shu C, Qin N. Multiquadric finite difference (MQ-FD) methods and its application. Adv Appl Math Mech 2009;1(5):615–38. [12] ul Islam S, Ahmad I. Local meshless method for PDEs arising from models of wound healing. Appl Math Model 2017;48(Suppl C):688–710. doi:10.1016/j.apm.2017.04.015. [13] Golbabai A, Mohebianfar E. A new method for evaluating options based on multiquadric RBF-FD method. Appl Math Comput 2017;308:130–41. doi:10.1016/j.amc.2017.03.019. [14] Shankar V. The overlapped radial basis function-finite difference (RBFFD) method: a generalization of RBF-FD. J Comput Phys 2017;342:211–28. doi:10.1016/j.jcp.2017.04.037. [15] Dehghan M, Mohammadi V. A numerical scheme based on radial basis function finite difference (RBF-FD) technique for solving the high-dimensional nonlinear schrodinger equations using an explicit time discretization: Runge kutta method. Comput Phys Commun 2017;217:23–34. doi:10.1016/j.cpc.2017.03.012. [16] Martin B, Fornberg B. Using radial basis function-generated finite differences (RBFFD) to solve heat transfer equilibrium problems in domains with interfaces. Eng Anal Bound Elem 2017;79:38–48. doi:10.1016/j.enganabound.2017.03.005. [17] Bayona V, Flyer N, Lucas G, Baumgaertner A. A 3-D RBF-FD solver for modeling the atmospheric global electric circuit with topography (gec-rbffd v1.0). Geosci Model Dev 2015;8(10):3007–20. doi:10.5194/gmd-8-3007-2015. [18] Kadalbajoo M, Kumar A, Tripathi L. Application of the local radial basis functionbased finite difference method for pricing american options. Int J Comput Math 2015;92(8):1608–24. doi:10.1080/00207160.2014.950571. [19] Martin B, Fornberg B, St-Cyr A. Seismic modeling with radial-basisfunction-generated finite differences. Geophysics 2015;80(4):T137–46. doi:10.1190/GEO2014-0492.1. [20] Tillenius M, Larsson E, Lehto E, Flyer N. A scalable RBF-FD method for atmospheric flow. J Comput Phys 2015;298:406–22. doi:10.1016/j.jcp.2015.06.003. [21] Ahmad I, ul Islam S, Khaliq A. Local RBF method for multi-dimensional partial differential equations. Comput Math Appl 2017;74(2):292–324. doi:10.1016/j.camwa.2017.04.026. [22] Khan W, ul Islam S, Ullah B. Analysis of meshless weak and strong formulations for boundary value problems. Eng Anal Bound Elem 2017;80(Suppl C):1–17. doi:10.1016/j.enganabound.2017.03.010. [23] Hardy RL. Research results in the application of multiquadric equations to surveying and mapping problems. Surv Mapp 1975;35(4):321–32. [24] Franke R. Scattered data interpolation: tests of some methods. Math Comput 1982;38(157):181–200. doi:10.1090/S0025-5718-1982-0637296-4. [25] Ballestra L, Pacelli G. A radial basis function approach to compute the firstpassage probability density function in two-dimensional jump-diffusion models for financial and other applications. Eng Anal Bound Elem 2012;36(11):1546–54. doi:10.1016/j.enganabound.2012.04.011. [26] Roque C, Martins P. Maximization of fundamental frequency of layered composites using differential evolution optimization. Compos Struct 2018;183:77–83. doi:10.1016/j.compstruct.2017.01.037. [27] Golberg M, Chen C, Karur S. Improved multiquadric approximation for partial differential equations. Eng Anal Bound Elem 1996;18(1):9–17. doi:10.1016/S0955-7997(96)00033-1. [28] Rippa S. An algorithm for selecting a good value for the parameter c in radial basis function interpolation. Adv Comput Math 1999;11(2–3):193–210. [29] Fasshauer G, Zhang J. On choosing “optimal” shape parameters for RBF approximation. Numer Algorithms 2007;45(1–4):345–68. doi:10.1007/s11075-007-9072-8. [30] Roque C, Ferreira A. Numerical experiments on optimal shape parameters for radial basis functions. Numer Methods Partial Differ Equ 2010;26(3):675–89. doi:10.1002/num.20453. [31] Roque C, Madeira J, Ferreira A. Multiobjective optimization for node adaptation in the analysis of composite plates using a meshless collocation method. Eng Anal Bound Elem 2015;50:109–16. doi:10.1016/j.enganabound.2014.08.006. [32] Roque C, Madeira J, Ferreira A. Node adaptation for global collocation with radial basis functions using direct multisearch for multiobjective optimization. Eng Anal Bound Elem 2014;39(1):5–14. doi:10.1016/j.enganabound.2013.10.012. [33] Luh L-T. The mystery of the shape parameter III. Appl Comput Harmon Anal 2016;40(1):186–99. doi:10.1016/j.acha.2015.05.001. [34] Flyer N, Barnett GA, Wicker LJ. Enhancing finite differences with radial basis functions: experiments on the Navier–Stokes equations. J Comput Phys 2016;316:39–62. doi:10.1016/j.jcp.2016.02.078. [35] Bayona V, Moscoso M, Kindelan M. Optimal constant shape parameter for multiquadric based RBF-FD method. J Comput Phys 2011;230(19):7384–99. doi:10.1016/j.jcp.2011.06.005. [36] Bayona V, Moscoso M, Kindelan M. Optimal variable shape parameter for multiquadric based RBF-FD method. J Comput Phys 2012;231(6):2466–81. doi:10.1016/j.jcp.2011.11.036. [37] Kansa EJ. Multiquadrics. a scattered data approximation scheme with applications to computational fluid-dynamics. i. surface approximations and partial derivative estimates. Comput Math Appl 1990a;19(8–9):127–45. [38] Kansa EJ. Multiquadrics. a scattered data approximation scheme with applications to computational fluid-dynamics. ii. solutions to parabolic, hyperbolic and elliptic partial differential equations. Comput Math Appl 1990b;19(8–9):147–61. [39] Sarra S. A local radial basis function method for advection-diffusion-reaction equations on complexly shaped domains. Appl Math Comput 2012;218(19):9853–65. doi:10.1016/j.amc.2012.03.062.

1.5

1

0.5

0

-0.5

-1

-1.5 -2

-1

0

1

2

3

Fig. 11. Optimization problem 3. Optimized distance for each stencil. 𝑑𝑚𝑖𝑛 = 0.5; 𝑑𝑚𝑎𝑥 = 4.5. Table 1 Errors for the best solutions, composite plate in bending problem. Optimization problem Problem 1 Problem 2 Problem 3

[m5GeSdc;November 24, 2017;9:40]

Shape parameter, c √ 2∕ 𝑛𝑖 GLODS √ 2∕ 𝑛𝑖 GLODS

Distance, d

Error

0.2 0.2 GLODS GLODS

0.0989 0.0699 0.0160 0.0023

Acknowledgment This work was supported by FCT, under LAETA, project UID/EMS/50022/2013. The support of LAETA to project “Aplicação de optimização global ao método numérico sem malha RBFFD para estudo de materiais compósitos” is also gratefully acknowledged. The support of Ministerio da Ciencia Tecnologia e do Ensino Superior and Fundo Social Europeu (MCTES and FSE) under programs POPH-QREN and Investigador FCT is also acknowledged. References [1] Medina D, Chen J. Three-dimensional simulations of impact induced damage in composite structures using the parallelized SPH method. Compos Part A Appl Sci Manuf 2000;31(8):853–60. doi:10.1016/S1359-835X(00)00031-2. [2] Liszka T, Orkisz J. The finite difference method at arbitrary irregular grids and its application in applied mechanics. Comput Struct 1980;11(1–2):83–95. doi:10.1016/0045-7949(80)90149-2. [3] Belytschko T, Lu Y, Gu L. Element free galerkin methods. Int J Numer Methods Eng 1994;37(2):229–56. doi:10.1002/nme.1620370205. [4] Mendonca P, De Barcellos C, Duarte A. Investigations on the hp-cloud method by solving timoshenko beam problems. Comput Mech 2000;25(2):286–95. [5] Garcia O, Fancello E, De Barcellos C, Duarte C. hp-clouds in mindlin’s thick plate model. Int J Numer Methods Eng 2000;47(8):1381–400. [6] Bui TQ, Nguyen TN, Nguyen-Dang H. A moving kriging interpolation-based meshless method for numerical simulation of kirchhoff plate problems. Int J Numer Methods Eng 2009;77(10):1371–95. doi:10.1002/nme.2462. [7] Sadamoto S, Ozdemir M, Tanaka S, Taniguchi K, Yu TT, Bui TQ. An effective meshfree reproducing kernel method for buckling analysis of cylindrical shells with and without cutouts. Comput Mech 2017;59(6):919–32. doi:10.1007/s00466-017-1384-5. [8] Roque CMC, Ferreira AJM, Jorge RMN. A radial basis function approach for the free vibration analysis of functionally graded plates using a refined theory. J Sound Vib 2007;300(3–5):1048–70. [9] Roque CMC, Ferreira AJM, Jorge RMN. Free vibration analysis of composite and sandwich plates by a trigonometric layerwise deformation theory and radial basis functions. J Sandw Struct Mater 2006;8(6):497–515. 9

JID: EABE

ARTICLE IN PRESS

C.M.C. Roque, J.F.A. Madeira

[m5GeSdc;November 24, 2017;9:40] Engineering Analysis with Boundary Elements 000 (2017) 1–10

[40] Custodio A, Madeira J. Glods: global and local optimization using direct search. J Global Optim 2015;62(1):1–28. doi:10.1007/s10898-014-0224-9. [41] Micchelli CA. Interpolation of scattered data distance matrices and conditionally positive definite functions. Constr Approx 1986;2(1):11–22. [42] Tolstykh A, Lipavskii M, Shirobokov D. High-accuracy discretization methods for solid mechanics. Arch Mech 2003;55(5–6):531–53. [43] Cecil T, Qian J, Osher S. Numerical methods for high dimensional Hamilton— Jacobi equations using radial basis functions. J Comput Phys 2004;196(1):327– 347. [44] Wright G, Fornberg B. Scattered node compact finite difference-type formulas generated from radial basis functions. J Comput Phys 2006;212(1):99–123. [45] Fornberg B, Wright G, Larsson E. Some observations regarding interpolants in the limit of flat radial basis functions. Comput Math Appl 2004;47:37–55. [46] Yu T, Yin S, Bui TQ, Xia S, Tanaka S, Hirose S. NURBS-based isogeometric analysis of buckling and free vibration problems for laminated composites plates with complicated cutouts using a new simple FSDT theory and level set method. Thin-Walled Struct 2016;101(Suppl C):141–56. doi:10.1016/j.tws.2015.12.008.

[47] Vu T-V, Nguyen N-H, Khosravifard A, Hematiyan M, Tanaka S, Bui TQ. A simple FSDT-based meshfree method for analysis of functionally graded plates. Eng Anal Bound Elem 2017;79(Suppl C):1–12. doi:10.1016/j.enganabound.2017.03.002. [48] Liu S, Yu T, Bui TQ, Yin S, Thai D-K, Tanaka S. Analysis of functionally graded plates by a simple locking-free quasi-3D hyperbolic plate isogeometric method. Compos Part B Eng 2017;120(Suppl C):182–96. doi:10.1016/j.compositesb.2017.03.061. [49] Bui TQ, Do TV, Ton LHT, Doan DH, Tanaka S, Pham DT, et al. On the high temperature mechanical behaviors analysis of heated functionally graded plates using fem and a new third-order shear deformation plate theory. Compos Part B Eng 2016;92(Suppl C):218–41. doi:10.1016/j.compositesb.2016.02.048. [50] Reddy JN. Mechanics of laminated composite plates and shells. CRC Press; 2004. [51] Reddy JN. Simple higher-order theory for laminated composite plates. J Appl Mech Trans ASME 1984;51(4):745–52. [52] Reddy JN. Mechanics of laminated composite plates. New York: CRC Press; 1997. [53] Zitnan P. An overdetermined b-spline collocation method for poisson problems on complex domains. Eng Anal Bound Elem 2013;37(5):860–7. doi:10.1016/j.enganabound.2013.03.002.

10