A Hermitian C Differential Reproducing Kernel Interpolation Meshless Method for the 3D Microstructure-Dependent Static Flexural Analysis of Simply Supported and Functionally Graded Microplates

1  Introduction

With the increasing demand for microstructures in industry and the rapid progress in material manufacturing technology, functionally graded (FG) structures have gradually shrunk from the macro scale to the micron scale. FG microstructures are gradually being used in cutting-edge technology fields, including thin films [1,2], micro-electro-mechanical systems [3,4], and atomic force microscopes [5,6]. Thus, developing an effective computational method to investigate the mechanical behavior of these microstructures has attracted considerable attention.

It is well-known that the mechanical behavior of FG macrostructures will be changed as their dimensions shrink from the macro-scale to the micro-scale [7]. The existing shell, plate, and beam theories based on the classical continuum mechanics (CCM) are inappropriate for use to analyze the dynamic and static responses of FG microshells, microplates, and microbeams due to the microstructure-dependent effect becoming significant. As a result, some non-CCM-based theoretical methods accounting for the microstructure-dependent impact have been proposed to investigate the mechanical behavior of microstructures. These theoretical methods include the couple stress theory (CST) [8,9], the strain gradient theory (SGT) [10,11], the doublet mechanics theory [12], the micropolar elasticity theory [13,14], and the nonlocal elasticity theory [15].

Hadjesfandiari et al. [16,17] and Yang et al. [18] established the consistent CST (CCST) and the modified CST (MCST) by assuming the couple-stress tensor is skew-symmetric and symmetric, respectively. As a result, instead of two material length-scale coefficients, which are required to study an elastic isotropic body’s mechanical behavior when the original CST is employed, only one material length-scale coefficient is needed when the MCST/CCST is employed. This facilitates their future application.

Within the CCST/MCST framework, some two-dimensional (2D) shear deformation theories for investigating the microstructure-dependent mechanical behavior of FG microplates/microshells have been developed by assuming particular kinematics models a priori. Beni et al. [19] presented a microstructure-dependent classical shell theory on the basis of the MCST to determine an FG circular cylindrical microshell’s smallest natural frequency and its corresponding wave number pair. Incorporating Mindlin’s kinematics model into the MCST, Ma et al. [20] established a microstructure-dependent first-order shear deformation theory (FOSDT) to analyze a homogeneous isotropic microplate’s flexural and free vibration behaviors. Arefi et al. [21] developed a novel shear deformation theory on the basis of the MCST to examine a three-layered microplate’s stress and displacement, for which the microplate of interest consists of an exponentially graded (EG) core and two piezomagnetic face sheets. Based on Hamilton’s principle combined with Reddy’s kinematics model, Lei et al. [22] and Thai et al. [23] presented a microstructure-dependent refined shear deformation theory (RSDT) on the basis of the MCST to conduct an FG microplate’s microstructure-dependent deformation and natural frequency behavior analyses. Kim et al. [24] developed a microstructure-dependent and MCSD-based third-order shear deformation theory (TOSDT) to investigate an FG microplate’s static buckling, static flexural, and free vibration behaviors. Thai et al. [25] developed a microstructure-dependent sinusoidal shear deformation theory (SSDT) on the basis of the MCST to examine an FG microplate’s static flexural and free vibration behaviors. Sobhy et al. [26] presented a microstructure-dependent and MCST-based trigonometric shear deformation theory (TSDT) with four primary variables for modeling an EG microplate’s static buckling, static flexural, and free vibration characteristics resting on Pasternak’s foundation.

Unlike these 2D microstructure-dependent and MCST-based shear deformation theories mentioned above, Wu et al. [27] established the unified microstructure-dependent shear deformation theories based on the CCST to study an FG/EG elastic microplate’s mechanical behavior. Their results showed that the CCST and MCST solutions of deformation and natural frequency associated with out-of-plane vibration modes are almost identical when setting the value of MCST’s material length-scale parameter at twice that of CCST’s material length-scale parameter. However, their solutions of natural frequency associated with the in-plane vibration modes are slightly different.

Instead of the CCST and MCST, other non-CCM-based analytical and numerical methods, including the SGT, the differential quadrature method (DQM), the iso-geometric analysis technique (IGAT), etc., have also been employed to study an FG microplate’s mechanical behavior. Incorporating Kirchhoff-Love’s kinematics model into the SGT, Deng et al. [28] established a non-CCM-based theory to determine an FG microplate’s smallest natural frequency with variable thickness. Within the SGT framework, Balobanov et al. [29] presented a microstructure-dependent classical thin shell theory to investigate a circular cylindrical microscale shell’s static flexural behavior. Integrating the advantages of the finite element method (FEM) and the DQM, Zhang et al. [30] developed a Hermitian C1 four-node quadrilateral element for conducting a moderately thick microplate’s mechanical behavior analysis. Integrating the MCST and the IGAT, Thanh et al. [31] established a seventh-order shear deformation theory for analyzing a porous FG microplate’s microstructure-dependent nonlinear thermal stability behavior. Nguyen et al. [32] developed a computational approach for analyzing an FG microplate’s geometrically nonlinear behavior on the basis of the IGAT and the RSDT. In conjunction with a modified nonlocal CST and the IGAT, Pham et al. [33] conducted an FG microplate’s static flexural and free vibration characteristics analyses, where the microplate rested on an elastic foundation. Based on the CCST, Wu and his colleagues [34,35] established a semi-analytical Hermitian Cn FEM to conduct elastic and piezoelectric microscale plates’/shells’ microstructure-dependent static and dynamic behavior analyses.

Meshless methods have also been employed to investigate a microscale structure’s mechanical behavior. Incorporating Mindlin’s kinematics model and radial basis functions into the MCST, Roque et al. [36] proposed a point collocation method for analyzing a homogeneous isotropic microplate’s static flexural behavior. Incorporating HSDT’s kinematics model into the MCST, Tran et al. [37] and Thai et al. [38] presented a moving Kriging interpolation meshless method to investigate an FG sandwich microplate’s static buckling, static flexural, and free vibration characteristics. Nguyen et al. [39] incorporated a four-variable kinematics model into the MCST to develop a non-uniform rational B-splines (NURBS) meshless method, which was used to investigate an FG microplate’s microstructure-dependent static buckling, static flexural, and free vibration behaviors. Finally, Thai et al. [40] employed the NURBS meshless method to conduct an FG microplate’s static buckling and free vibration behavior analyses.

In their series of papers, Li et al. [41], Simkins et al. [42], Liu et al. [43], and Lu et al. [44] established the reproducing kernel element method to solve Galerkin weak forms of a system of higher order partial differential equations which are associated with Dirichlet boundary conditions.

Chen et al. [45] and Wang et al. [46] established the Hermitian C1 and Lagrange C0 differential reproducing kernel interpolation meshless (DRKIM) methods, respectively, for investigating laminated composite and FG macroscale structures’ mechanical behavior. The novelty of these DRKIM methods is that the shape functions of the differential reproducing kernel (DRK) interpolant’s derivatives are obtained by setting up the differential reproducing conditions (DRCs) without using direct differentiation, as is necessary for the conventional reproducing kernel interpolation and approximation methods [47]. It has been shown that the solutions obtained using these DRKIM methods closely agree with the available 3D solutions of macroscale plates, rather than those of microplates.

As we see in the strong form of the 3D CCST, the primary variables’ highest order derivative is the third order for microplates, which differs from the first order designation for macroscale plates. This situation will reduce the accuracy of the early proposed Hermitian C1 and Lagrangian C0 DRKIM methods and slow down their convergence rate. To enhance these DRKIM methods’ accuracy and speed up their convergence rate, in this paper, we aim to establish a Hermitian C2 DRKIM method by making some modifications: the Hermitian C2 DRK interpolant should satisfy the nodal interpolation properties and the continuity conditions up to primary variables’ second-order derivatives at each sampling node. Moreover, we also aim to establish the Hermitian C2 DRKIM method, which is a point collocation method, by incorporating our Hermitian C2 DRK interpolant into the strong form of the 3D CCST to carry out an FG microplate’s 3D microstructure-dependent static flexural analysis. After validating the Hermitian C2 DRKIM’s accuracy using the relevant 3D solutions reported in the literature, we will carry out a parametric study for the FG microplate’s 3D microstructure-dependent static flexural behavior to examine how the impact of essential factors affects the induced deformations, in-plane stresses, transverse stresses, and couple stresses, including the length-to-thickness ratio, the material length-scale parameter, and the inhomogeneity index.

2  The Hermitian C2 DRKIM Method

2.1 The Hermitian C2 DRK Interpolant

We consider np discrete sampling nodes placed at ξ = ξ1,  ξ2, ⋯,  ξnp, respectively, in a natural coordinate system ξ, the values of which are ξ1=−1, ξnp=1, and the others are randomly selected between −1 and 1. The Hermitian C2 DRK interpolant fh(ξ) and its first- and second-order derivatives, dfh(ξ)/dξ and d2fh(ξ)/dξ2 (i.e., θh(ξ) and κh(ξ)), are required to satisfy the nodal interpolation properties. Thus, it is defined as

fh(ξ)=∑l=1np  [Nl(ξ)fl+N^l(ξ)θ l+N¯l(ξ) κl]=∑l=1np { [ϕl(ξ)+ψl(ξ)]fl+[ϕ^l(ξ)+ψ^l(ξ)]θl+[ϕ¯l(ξ)+ψ¯l(ξ)]κl},  (1)

where Nl(ξ), N^l(ξ), and N¯l(ξ) (l = 1, 2,…, np) denote the Hermitian C2 DRK interpolant’s shape functions at ξ = ξl; fl, θl, and κl are the nodal values of fh(ξ), θh(ξ), and κh(ξ) at ξ=ξl, respectively. ψl(ξ), ψ^l(ξ), and ψ¯l(ξ) (l = 1, 2,…, np) are the primitive functions for fh(ξ), θh(ξ), and κh(ξ), respectively, which are selected to satisfy the Kronecker delta properties. The primitive functions chosen in this article are ψl(ξ)=wq(ξ), ψ^l(ξ)=(ξ−ξl)wq(ξ), and ψ¯l(ξ)=(ξ−ξl)2wq(ξ)/2, in which wq(ξ) is defined as a normalized eighth-degree (octic) polynomial with the support size a0=(0.99)min(|ξl−ξl+1|,|ξl−ξl−1|), such that these primitive functions and their first- and second-derivatives satisfy the Kronecker delta properties (i.e., ψl(ξk)=δlk, dψ^l(ξk)/dξ=δlk, d2ψ¯l(ξk)/dξ2=δlk, ψ^l(ξk)=ψ¯l(ξk)=0, dψl(ξk)/dξ=dψ¯l(ξk)/dξ=0, and d2ψl(ξk)/dξ2=d2ψ^l(ξk)/dξ2=0). The symbols ϕl(ξ), ϕ^l(ξ), and ϕ¯l(ξ) (l = 1, 2,…, np) are defined as the enrichment functions for fh(ξ), θh(ξ), and κh(ξ), respectively, which are determined by setting up the nth-order DRCs. In our Hermitian C2 DRK interpolant, the enrichment functions are arranged as ϕl(ξ)=wa(ξ−ξl)PT(ξ−ξl)b0c2(ξ), ϕ^l(ξ)=wa(ξ−ξl)P^T(ξ−ξl) b0c2(ξ), and ϕ¯l(ξ)=wa(ξ−ξl)P¯T(ξ−ξl) b0c2(ξ), in which b0c2(ξ) and wa(ξ−ξl) denote the undetermined function vector and a Gaussian function, respectively, and

PT(ξ−ξl)={1(ξ−ξl)(ξ−ξl)2⋯(ξ−ξl)n},(2)

P^T(ξ−ξl)=(−1)dPT(ξ−ξl)/d(ξ−ξl)=dPT(ξ−ξl)/dξl={0,−1,−2(ξ−ξl),−3(ξ−ξl)2,⋯,−n(ξ−ξl)n−1},(3)

P¯T(ξ−ξl)=(−1)2d2PT(ξ−ξ)/d(ξ−ξl)2=d2PT(ξ−ξl)/dξl2={0,0,2,6(ξ−ξl),⋯,n(n−1)(ξ−ξl)n−2}.(4)

In order to determine the undetermined function vector b0c1(z), we select the complete nth-order polynomials as the basis functions that are to be reproduced and set up (n + 1) DRCs as follows:

∑l=1np{[ϕl(ξ)+ψl(ξ)]ξlm+[ϕ^l(ξ)+ψ^l(ξ)]m ξlm−1+[ϕ¯l(ξ)+ψ¯l(ξ)]m (m−1)ξlm−2}=ξmm=0,1, ⋯,n.(5)

Eq. (5) can be rearranged in the following explicit forms:

m=0:∑l=1np[ϕl(ξ)+ψl(ξ)]=1⇒∑l=1np ϕl(ξ)=1−∑l=1np ψl(ξ),(6)

m=1: ∑l=1np[ϕl(ξ)+ψl(ξ)]ξl+∑l=1np[ϕ^l(ξ)+ψ^l(ξ)]=ξ⇒∑l=1np (ξ−ξl) ϕl(ξ)+∑l=1np(−1)ϕ^l(ξ)=−∑l=1np (ξ−ξl) ψl(ξ)−∑l=1np(−1)ψ^l(ξ),(7)

m=2∑l=1np[ϕl(ξ)+ψl(ξ)]ξl2+∑l=1np[ϕ^l(ξ)+ψ^l(ξ)](2ξl) +∑l=1np[ϕ¯l(ξ)+ψ¯l(ξ)](2)=ξ2⇒∑l=1np (ξ−ξl)2ϕl(ξ)+∑l=1np(−2)(ξ−ξl)  ϕ^l(ξ)+∑l=1np2ϕ¯l(ξ)=−∑l=1np (ξ−ξl)2ψl(ξ)−∑l=1np(−2)(ξ−ξl)  ψ^l(ξ)−∑l=1np2ψ¯l(ξ),(8)

⋮

m=n:∑l=1np[ϕl(ξ)+ψl(ξ)]ξln+∑l=1np[ϕ^l(ξ)+ψ^l(ξ)](n ξln−1)+∑l=1np [ϕ¯l(ξ)+ψ¯l(ξ)][n(n−1) ξ ln−2]=ξn⇒∑l=1np (ξ−ξl)n ϕl(ξ)+∑l=1np(−n) (ξ−ξl)n−1ϕ^l(ξ)+∑l=1np(n) (n−1)(ξ−ξl)n−2ϕ¯l(ξ)=−∑l=1np (ξ−ξl)n ψl(ξ)−∑l=1np(−n) (ξ−ξl)n−1ψ^l(ξ)−∑l=1np(n) (n−1)(ξ−ξl)n−2ψ¯l(ξ).(9)

We rewrite Eqs. (6)–(9) in matrix form as follows:

∑l=1np P(ξ−ξl) ϕl (ξ)+∑l=1np P^(ξ−ξl)  ϕ^l (ξ)+∑l=1np P¯(ξ−ξl)  ϕ¯l (ξ)=P(0)−∑l=1np P(ξ−ξl)  ψl (ξ)−∑l=1np P^(ξ−ξl)  ψ^l (ξ)−∑l=1np P¯(ξ−ξl)  ψ¯l (ξ),(10)

where P(0)=[100⋯0]T.

We substitute the enrichment functions into the DRCs to yield the following expression for the undetermined function vector b0c2(ξ):

b0c2(ξ)=Ac2−1(ξ)[P(0)−∑l=1np P(ξ−ξl)  ψl (ξ)−∑l=1np P^(ξ−ξl)  ψ^l (ξ)−∑l=1np P¯(ξ−ξl)  ψ¯l (ξ)],(11)

where Ac2(ξ)=∑l=1np [P(ξ−ξl) wa(ξ−ξl) PT(ξ−ξl)+P^(ξ−ξl) wa(ξ−ξl) P^T(ξ−ξl)+P¯(ξ−ξl)wa(ξ−ξl) P¯T(ξ−ξl)].

We substitute Eq. (11) into Eq. (1) to yield the shape functions of the Hermitrian C2 DRK interpolant as follows:

Nl(ξ)=ϕl(ξ)+ψl (ξ)(l=1,2,…,np),(12)

N^l(ξ)=ϕ^l(ξ)+ψ^l(ξ)(l=1,2,…,np),(13)

N¯l(ξ)=ϕ¯l(ξ)+ψ¯l(ξ)(l=1,2,…,np),(14)

where ϕl(ξ)=wa(ξ−ξl)PT(ξ−ξl)b0c2(ξ),

ϕ^l(ξ)=wa(ξ−ξl)P^T(ξ−ξl)b0c2(ξ),

ϕ¯l(ξ)=wa(ξ−ξl)P¯T(ξ−ξl)b0c2(ξ).

Eq. (11) shows that the enrichment functions vanish at all the sampling points (i.e., ϕl (ξk)=ϕ^l (ξk)=ϕ¯l (ξk)=0, for all l and k = 1, 2, …, np). When we select the primitive functions mentioned above for fh(ξ), such that ψl (ξk)=δlk, ψ^l (ξk)=0, and ψ¯l (ξk)=0 a priori, the shape functions will satisfy the Kronecher delta properties, which are Nl (ξk)=δlk, N^l (ξk)=0, and N¯l (ξk)=0

2.2 The Hermitian C2 DRK Interpolant’s Derivatives

The Hermitian C2 DRK interpolant fh(ξ) in Eq. (1) has the first-order derivative with respect to ξ:

dfh(ξ)dξ=∑l=1np  [Nl(1)(ξ)fl+N^l(1)(ξ)θl+N¯l(1)(ξ)κl]=∑l=1np [ (ϕl(1)(ξ)+ψl(1)(ξ))  fl+(ϕ^l(1)(ξ)+ψ^l(1)(ξ))  θl+ (ϕ¯l(1)(ξ)+ψ¯l(1)(ξ))  κl] ,(15)

where Nl(1)(ξ), N^l(1)(ξ), and N¯l(1)(ξ) (l = 1, 2,…, np) are the shape functions of the Hermitian C2 DRK interpolant’s first-order derivative at the node ξ = ξl, which satisfy the Kronecker delta properties; ψl(1)(ξ), ψ^l(1)(ξ), and ψ¯l(1)(ξ) (l = 1, 2,…, np) are primitive functions’ first-order derivatives (i.e., ψl(1)(ξ)=d ψl(ξ)/dξ, ψ^l(1)(ξ)=d ψ^l(ξ)/dξ, and ψ¯l(1)(ξ)=d ψ¯l(ξ)/dξ); and ϕl(1)(ξ), ϕ^l(1)(ξ), and ϕ¯l(1)(ξ) (l = 1, 2,…, np) denote enrichment functions’ first-order derivatives, which are obtained by imposing the nth-order DRCs, and are expressed as ϕl(1)(ξ)=wa(ξ−ξl)PT(ξ−ξl) b1c2(ξ), ϕ^l(1)(ξ)=wa(ξ−ξl)P^T(ξ−ξl) b1c2(ξ), and ϕ¯l(1)(ξ)=wa(ξ−ξl)P¯T(ξ−ξl) b1c2(ξ), for which b1c2(ξ) is the undetermined function vector.

In order to determine the undetermined functions b1c2(ξ) in Eq. (15), again, we select the complete nth-order polynomials as the basis functions to be reproduced and set up (n + 1) DRCs as follows:

∑l=1np {[ϕl(1)(ξ)+ψl(1)(ξ)]ξlm+[ϕ^l (1)(ξ)+ψ^l(1)(ξ)]m ξlm−1+[ϕ¯l (1)(ξ)+ψ¯l(1)(ξ)]m (m−1)ξlm−2}=m ξm−1,(16)

where m=0,1,2,  ⋯,  n.

We rearrange Eq. (16) in the explicit forms as follows:

m=0:∑l=1np [ϕl(1)(ξ)+ψl(1)(ξ)]=0⇒∑l=1np ϕl(1)(ξ)=−∑l=1np ψl(1)(ξ),(17)

m=1:∑l=1np [ϕl(1)(ξ)+ψl(1)(ξ)]  ξl+∑l=1np [ϕ^l(1)(ξ)+ψ^l(1)(ξ)]=1⇒∑l=1np (ξ−ξl)ϕl(1)(ξ)+∑l=1np (−1)ϕ^l(1)(ξ)=−1−∑l=1np (ξ−ξl)ψl(1)(ξ)−∑l=1np (−1)ψ^l(1)(ξ),(18)

m=2:∑l=1np [ϕl(1)(ξ)+ψl(1)(ξ)]  ξl2+∑l=1np [ϕ^l(1)(ξ)+ψ^l(1)(ξ)]  (2ξl)+∑l=1np [ϕ¯l(1)(ξ)+ψ¯l(1)(ξ)]  (2)=2ξ⇒∑l=1np (ξ−ξl)2ϕl(1)(ξ)+∑l=1np (−2)(ξ−ξl)  ϕ^l(1)(ξ)+∑l=1np (2)ϕ¯l(1)(ξ)=−∑l=1np (ξ−ξl)2ψl(1)(ξ)−∑l=1np (−2)(ξ−ξl) ψ^l(1)(ξ)−∑l=1np (2)ψ¯l(1)(ξ),(19)

⋮

m=n:∑l=1np [ϕl(1)(ξ)+ψl(1)(ξ)]  ξln+∑l=1np [ϕ^l(1)(ξ)+ψ^l(1)(ξ)]  n ξln−1+∑l=1np [ϕ¯l(1)(ξ)+ψ¯l(1)(ξ)]  n (n−1)ξln−2=n ξn−1⇒∑l=1np (ξ−ξl)nϕl(1)(ξ)+∑l=1np (−n) (ξ−ξl)n−1ϕ^l(1)(ξ)+∑l=1np (n) (n−1)(ξ−ξl)n−2ϕ¯l(1)(ξ)=−∑l=1np (ξ−ξl)nψl(1)(ξ)−∑l=1np (−n)(ξ−ξl)n−1ψ^l(1)(ξ)−∑l=1np n (n−1)(ξ−ξl)n−2ψ¯l(1)(ξ).(20)

We rewrite the above Eqs. (17)–(20) in matrix form as follows:

∑l=1np P(ξ−ξl)  ϕ l(1)(ξ)+∑l=1np P^(ξ−ξl)  ϕ^l(1)(ξ)+∑l=1np P¯(ξ−ξl)  ϕ¯l(1)(ξ)=P^(0)−∑l=1np P(ξ−ξl)  ψl(1)(ξ)−∑l=1np P^(ξ−ξl)  ψ^l(1)(ξ)−∑l=1np P¯(ξ−ξl)  ψ¯l(1)(ξ),(21)

where P^(0)=dP(0)/dξl=[0−10⋯0]T.

We substitute the enrichment functions into the DRCs to yield the undetermined function vector b1c2(ξ) as follows:

b1c2(ξ)=Ac2−1(ξ)[P^(0)−∑l=1np P(ξ−ξl)  ψl(1)(ξ)−∑l=1np P^(ξ−ξl)  ψ^l(1)(ξ)−∑l=1np P¯(ξ−ξl)  ψ¯l(1)(ξ)].(22)

We substitute Eq. (22) into Eq. (15) to obtain the shape functions of the Hermitian C2 DRK interpolant’s first-order derivatives as follows:

Nl(1)(ξ)=ϕ l (1)(ξ)+ψ l (1)(ξ)(l=1,2,…,np),(23)

N^l(1)(ξ)=ϕ^l (1)(ξ)+ψ^l (1)(ξ)(l=1,2,…,np),(24)

N¯l(1)(ξ)=ϕ¯l (1)(ξ)+ψ¯l (1)(ξ)(l=1,2,…,np),(25)

where ϕ l(1)(ξ)=wa(ξ−ξl)PT(ξ−ξl)b1c2(ξ), ϕ^l(1)(ξ)=wa(ξ−ξl)P^T(ξ−ξl)b1c2(ξ), ϕ¯l(1)(ξ)=wa(ξ−ξl)P¯T(ξ−ξl)b1c2(ξ).

From Eqs. (23)–(25), it can be seen that the values of the enrichment functions’ first-order derivatives at all sampling nodes are zero (i.e., ϕ l(1)(ξk)=ϕ^l(1)(ξk)=ϕ¯l(1)(ξk)=0). Subsequently, suppose we select the first-order primitive functions for dfh(ξ)/dξ such that ψl(1)(ξk)=0, ψ^l(1)(ξk)=δlk, and ψ¯l(1)(ξk)=0, a priori. Finally, the above shape functions satisfy the Kronecker delta properties (i.e., Nl(1)(ξk)=0, N^l(1)(ξk)=δlk, and N¯l(1)(ξk)=0).

Similarly, the above derivation procedure can proceed to the rth-order derivative of the Hermitian C2 DRK interpolant fh (ξ), which is thus expressed in Appendix A.

2.3 Weight Functions and Primitive Functions

In implementing our Hermitian C2 DRKIM method, we must select the weight function and the primitive function in advance. This work uses the normalized Gaussian function as the weight function, which is expressed as follows [47]:

Normalized Gaussian function: wa(s)={e−(s/α)2−e−(1/α)21−e−(1/α)2fors≤10fors>1,(26)

where s=|ξ−ξl|/al, in which al denotes the support size at the reference sampling point l, and the value of α is set at α=0.3.

As mentioned above, we define the primitive functions for the Hermitian C2 DRK interpolant as ψl(ξ)=wq(ξ), ψ^l(ξ)=(ξ−ξl)wq(ξ), and ψ¯l(ξ)=(ξ−ξl)2wq(ξ)/2, respectively, for which wq(ξ) is a normalized eighth-degree (octic) polynomial, which is given as follows [47]:

wq(s)={−3s8+8s6−6s4+1fors≤10fors>1,(27)

where s=|ξ−ξl|/a0, in which a0 is defined as a0=(0.99)min(|ξl−ξl+1|,|ξl−ξl−1|) to ensure the Kronecker delta properties are satisfied (i.e., ψl(ξ=ξk)=δlk, dψ^l(ξk)/dξ=δlk, and d2ψ¯l(ξk)/dξ2=δlk).

It is noticed that for a meshless method, the support size al for the selected weight function wa(ξ) will not be a very small value, often resulting in numerical errors; whereas, it also has to be small enough to preserve the meshless method’s local character due to an increase in the support size also resulting in numerical errors. Chen et al. [45] and Wang et al. [46] thus recommended a compromise range of the value of al to ensure the Hermitian C1 and Lagrange C0 DRKIM methods’ accuracy and convergence rate. It has been recommended as follows: In the case of a uniform sampling node distribution, the appropriate value of al is al=(n+0.1)Δξ, where  Δξ denotes the spacing between the adjacent nodes, and the value of al is constant for each node. In the case of a randomly scattered node distribution, the appropriate value of al is selected to include (2n + 1) nodes and the value of al is variable for each node. This guidance is adopted in this paper.

To have a clear picture related to how the values of these shape functions vary in the natural coordinate, in Fig. 1a–c, we consider a case of 11 sampling nodes with uniform spacing and present the distributions of the enrichment function (ψ6(ξ)), the primitive function (ϕ6(ξ)), and the shape function (N6(ξ)) of node 6 along the natural coordinate axis, respectively, for which N6(ξ)=ϕ6(ξ)+ψ6(ξ). It can be seen in Fig. 1a–c that the Kronecker delta properties, ψ6(ξi)=δi6  and  N6(ξi)=δi6(i=1−11), are satisfied, and ϕ6(ξi)=0(i=1−11). Furthermore, in Fig. 2, we present the distribution of each sampling node’s shape function along the natural coordinate axis, i.e., Ni(ξ)(i=1−11). Again, each shape function is shown to satisfy the Kronecker delta properties and localize in a region of the support size.

Figure 1: Distributions of (a) the enrichment function, (b) the primitive function, and (c) the shape function of node 6 in the natural coordinate in the case of np = 11 with uniform spacing

Figure 2: Distributions of the shape functions (a) (i = 1–3), (b) (i = 4–6), (c) (i = 7–9), and (d) (i = 10 and 11), in the natural coordinate in the case of np = 11 with uniform spacing

3  3D Microstructure-Dependent Static Flexural Analysis of FG/EG Microplates

3.1 The Quasi-State Space Equations of the CCST

This work considers the 3D microstructure-dependent static flexural problem of a simply-supported FG microplate under either a sinusoidally distributed load or a uniform load, and the former loading case is shown in Fig. 3. The symbols h, Lx, and Ly represent the microplate’s height, length, and width, respectively. A Cartesian coordinate system (x, y and z) is oriented so that the xy-plane is the microplate’s mid-plane.

Figure 3: An FG microplate of interest that is subjected to a sinusoidally distributed load

The displacement vector u of the deformed microplate is expressed as u=ux i+uy j+uz k, where i, j, and k represent the unit basis vectors in the x, y, and z directions, respectively.

The strain tensor ε is symmetric, and its relationships with the displacement tensor in Cartesian coordinates are expressed as

εxx=ux,x,(28a)

εyy=uy,y,(28b)

εzz=uz,z,(28c)

γxz=2εxz=ux,z+uz,x,(28d)