Bending analysis of magnetoelectroelastic nanoplates resting on Pasternak elastic foundation based on nonlocal theory

Based on the nonlocal theory and Mindlin plate theory, the governing equations (i.e., a system of partial differential equations (PDEs) for bending problem) of magnetoelectroelastic (MEE) nanoplates resting on the Pasternak elastic foundation are first derived by the variational principle. The polynomial particular solutions corresponding to the established model are then obtained and further employed as basis functions with the method of particular solutions (MPS) to solve the governing equations numerically. It is confirmed that for the present bending model, the new solution strategy possesses more general applicability and superior flexibility in the selection of collocation points. Finally, the effects of different boundary conditions, applied loads, and geometrical shapes on the bending properties of MEE nanoplates are evaluated by using the developed method. Some important conclusions are drawn, which should be helpful for the design and applications of electromagnetic nanoplate structures.


Introduction
By virtue of the mechanical-electric-magnetic coupling effect, magnetoelectroelastic (MEE) solids have been widely used in smart devices [1] . In the past two decades, much attention has been paid to the bending analysis of the MEE plates [2][3][4][5] , and some achievements on the bending problems of MEE nanoplates have been made [6][7] . Considering the size effects of nanostructures, the nonlocal theory initiated by Eringen [8][9] has been successfully used to analyze magnetoelectric nanostructures in recent years [10][11][12][13] . Among them, Li et al. [10] and Guo et al. [12] investigated the buckling problems of the MEE nanoplates and the multilayered MEE nanoplates under the nonlocal theory, respectively. Li et al. [13] also studied the bending analysis of piezoelectric nanobeams. However, to the best of the authors' knowledge, up till now, the bending model of the MEE nanoplates based on the nonlocal theory has not been proposed yet, let alone the corresponding problem of the MEE nanoplates resting on Pasternak elastic foundation.
On the other hand, the plate's bending problems can be actually reduced to solve different types of partial differential equations (PDEs) under the corresponding boundary conditions. Due to the limitations of analytical methods, numerical methods are necessary to be utilized. As one of them, the meshless methods including the wave-based method [14][15] , the boundary node method [16][17] , and the method of fundamental solutions [18][19] have been all developed for plate bending problems to a certain extent. The above methods seek the approximations with the expansion functions to satisfy the governing equations and then obtain the unknown weight parameters by requiring the approximations satisfying the given boundary conditions [20] . It is difficult to solve the problems with complex boundary conditions, such as Ref. [19], which utilizes the polynomial expansions approximations and simulates the static deformation of pyramids. In view of this, the method of particular solutions (MPS) [21][22][23][24] is further proposed to relax such strong requirement. Among them, the radial basis functions (RBFs) are more successful to find the particular solutions for some PDEs [22][23] ; however, the RBFs cannot solve the optimization problem of shape parameters well [22,25] . Additionally, although the particular solutions based on Chebyshey polynomials can avoid such difficulties and simultaneously possess relatively high accuracy [26][27] , there is an obvious disadvantage in using Chebyshey polynomials as basis functions, i.e., the forcing term of the corresponding differential equation needs to be smoothly extendable to the exterior of the solution domain in the case of irregular domains [20,25] , and the selection of the corresponding collocation points at this time must be treated specially. Recently, a novel method to derive the particular solutions to a PDE by using the standard polynomial basis functions was presented [25] , where the collocation points could be selected arbitrarily inside the solution domain. The method was recently used to investigate the bending problems of elastic plates, where the restriction that the constant term in the differential operator must be non-zero was relaxed [20] .
It is known that the polynomial basis functions are not ideal for a global approach since the resultant matrix is extremely ill-conditioned with the higher order of the polynomial basis [25] . Hence, the method of polynomial particular solutions are useless without proper treatment of the resultant matrix. In this paper, we adopt the multiple-scale technique [28][29] , which is a convenient and effective pre-conditioning technique to reduce the condition number of the resultant matrix of the MPS. Certainly, there also exist other superior methods to deal with the above ill-conditioned problem, such as the finite point method [30] and the localized method of fundamental solutions [31] . Besides, due to the utilization of the concept of localization [31] , the condition number can be markedly decreased, and the error can be analyzed theoretically. Therefore, introducing the concept of localization into the MPS will become an emphasis in our future works.
Finally, in this paper, we tentatively investigate the bending problems of the MEE nanoplates resting on the Pasternak elastic foundation based on the nonlocal theory by the method of polynomial particular solutions. Quite different from Refs. [20] and [25], the governing equations obtained here should be a system of PDEs.

Statement of problem and basic theories 2.1 Problem statement
As shown in Fig. 1, an MEE nanoplate with length l, width b, and thickness h rests on a Pasternak elastic foundation. A Cartesian coordinate system (x, y, z) is introduced with the z-axis along the plate thickness direction and the xy-plane sitting on the midplane of the undeformed plate. The transverse mechanical load q is distributed on the surfaces of the MEE plate. The transversely isotropic MEE solid is polarized in the z-direction. As mentioned in Refs. [1], [10], [32], and [33], for the MEE thin plates, applying electric and magnetic potentials ϕ 0 and φ 0 on the surface of the nanoplate should be more reasonable. For the present bending problem of the MEE nanoplate, the constitutive equations can be found in Refs. [1] and [10]. The electromagnetic boundary conditions can be described as follows [1,10] : As an example, the mechanical boundary conditions of an edge of the nanoplate parallel to the y-axis (i.e., x = 0 or x = l) can be given as w = ψ y = M xx = 0 for the simply supported boundary, M xx = M xy = Q xx = 0 for the free boundary, and w = ψ x = ψ y = 0 for the clamped boundary. w is the deflection of the MEE nonaplate. ψ x and ψ y are the rotations of the normal to the mid-plane about the x-axis and y-axis, respectively. M xx , M xy , and Q xx are the resultant bending moments and twisting moments per unit length in terms of stresses σ xx , σ xy , and σ zx , respectively [34] . For the three kinds of mechanical boundary conditions at one edge of the nanoplate parallel to the x-axis (i.e., y = 0 or y = b), the resultant bending moments and twisting moments can be similarly given, which are omitted here.

Nonlocal theory and Mindlin plate theory
The nonlocal elastic theory assumes that the stress state at a reference point ϑ is dependent not only on the strain state at ϑ but also on the strain state at all other points in the whole domain [9,35] . Thus, for the present multi-field problems, based on the nonlocal theory, the following equations should hold [12][13] : where e 0 is a constant related to the considered material, a is an internal characteristic length of the material, and e 0 a is the corresponding constant parameter which can be obtained by molecular dynamics, experimental studies, or molecular structure mechanics [9] . (σ nl ij , D nl i , B nl i ) and (σ ij , D i , B i ) are the components of the nonlocal and local generalized stresses (i.e., stresses, electric displacements, and magnetic inductions), respectively.
In the meantime, according to the Mindlin plate theory, the displacement components u x (x, y, z), u y (x, y, z), and u z (x, y, z) of the present problem can be expressed as Accordingly, the non-zero strains can be expressed as follows:

Derivation of governing equations
Considering the electromagnetic quasi-static approximation of the Maxwell's vector equations, the electric and magnetic fields can be simplified as the gradients of φ and ϕ, i.e., In the meantime, for the considered thin MEE nanoplate, the in-plane electric and magnetic fields can be ignored [1] , i.e., In order to obtain the governing equations of the present problem, the variational principle is utilized, where δU and δW are the variations of the strain energy and the external work, respectively.
Considering the constitutive equations and Eqs. (4) and (5), we can obtain the strain energy of the MEE nanoplates as where and k = π 2 /12 is the shear correction coefficient [1] . On the other hand, the external work can be written as where f is related to the Pasternak foundation and the applied generalized loads, and it can be expressed as where k w and k g are the spring and shear coefficients of the Pasternak foundation [36] , N x = e 31 φ 0 + f 31 ϕ 0 and N y = e 31 φ 0 + f 31 ϕ 0 are the electromagnetic forces in the xand y-directions, respectively, whilst e 31 and f 31 are explained later. Then, by using the variational principle, i.e., Eq. (6), the following equations can be easily obtained: Substituting Eqs. (4) and (5) and the constitutive equations into Eq. (12) and noting the electromagnetic boundary conditions (1), we can obtain the electric and magnetic fields as [10] ⎧ ⎪ ⎪ ⎨ where . e ij , f ij , and g ij denote the piezoelectric, piezomagnetic, and magnetoelectric constants, respectively, and h ij and μ ij are the dielectric and magnetic permeability coefficients, respectively.
Then, substituting Eq. (4) and the constitutive equations into Eq. (8) and considering the nonlocal theory, i.e., Eq. (2), we get where and c ij are the elastic constants. Substituting the above equation into Eq. (11) and introducing the following dimensionless quantities: , we derive the dimensionless governing equations for the first time as follows: Thus, the bending problem of the MEE nanoplate considered here is finally reduced to solving the system of PDEs (15) under the given corresponding boundary conditions.

The strategy of the MPS for solving the present system of PDEs
We first briefly describe the governing equations (15) as where the linear partial differential operators 1 , 2 , and 3 are defined by in which χ is a chosen supplementary parameter to avoid λ 3 = 0 when k w = 0, since λ 3 is the denominator in the following equation. Meanwhile, other operators and functions are given as Based on the MPS, the solutions to Eq. (16) can be approximated by the linear combination of particular solutions, i.e., Introducing the polynomial basis functions b ξ τ (x, y), the above particular solutions should satisfy , the superscript ξ and subscript τ denote the order and sequence of the monomial basis functions, respectively, denotes the dimension of the polynomial space, and = (ξ + 1)(ξ + 2)/2 for the two-dimensional problems.
where m and n are the exponents corresponding to x and y, respectively. With the similar method in Refs. [20] and [25], the τ th particular solutions (τ = τ ) of rotations and deflection can be derived from Eq. (18) as follows: Without loss of generality, the boundary conditions are assumed as where ∂Ω 1 and ∂Ω 2 simply denote two different boundaries, and the boundary operators ij and functions G i can be determined by the concrete boundary conditions. Finally, substituting Eqs. (17) and (18) into Eqs. (16) and (20) and applying the collocation method, we can obtain the following algebraic equations: where n i is the number of the points (x k , y k ) chosen inside the computational domain, n b1 and n b2 are the number of the points (x k1 , y k1 ) chosen on the boundary ∂Ω 1 and the number of (x k2 , y k2 ) on ∂Ω 2 , respectively. If 3 × n i + n b1 + n b2 = 3 , α τ , β τ , and γ τ can be easily obtained by solving the corresponding linear system. If 3 × n i + n b1 + n b2 > 3 , all of them can be then obtained from the corresponding over-determined system by using of the least-square approach. For convenience, the above algebraic equations can be further reformulated as a matrix equation, i.e., where Bending analysis of magnetoelectroelastic nanoplates resting on Pasternak elastic foundation 1777 To solve Eq. (22), the multiple-scale technique [28][29] is further employed to reduce the condition number of the coefficient matrix. After introducing the pre-conditioning matrix ω = diag(ω 1 , ω 2 , · · · , ω 3 ), Eq. (22) can be rewritten as Kω −1 ωP = Φ, where ω τ = 3×ni+n b1 +n b2 k=1 K 2 kτ , τ = 1, 2, · · · , 3 , and K kτ are the elements of the coefficient matrix K.
Finally, the unknown coefficients can be obtained by solving K P = Φ, in which K = Kω −1 and P = ωP. The introduction of the pre-conditioning matrix ω can make all the columns of the new coefficient matrix K treated by the multiple-scale technique own the same norm, and ensure K being well conditioned [20] .

Numerical examples and discussion
In this section, numerical examples for the MEE nanoplate bending problems are presented. The MEE composite made of the PE material BaTiO 3 as the inclusions and the PM material CoFe 2 O 4 as the matrix is considered, and the material properties of them are given in Refs. [10] and [37].

Verification of the present method and the determination of concrete computational parameters
For comparison, we first reduce the current MEE nanoplate problem to the model for an MEE nanobeam based on the nonlocal theory and Mindlin plate theory by removing the items involving the y-dependence in Eq. (15). Furthermore, to estimate the computational accuracy appropriately, the maximum absolute error (E MA ), the root-mean-squared error (E RMS ), and the relative error (E Rel ) are, respectively, introduced as follows [20,25,38] : where ψ x (j) and ψ x (j) are, respectively, the numerical and analytical solutions at the test points (x j , y j ), and n t is the number of the test points. The errors of the polynomial particular solutions method for simulating the bending of reduced nanobeam compared with the analytical results given in Ref. [13] are presented in Table 1, where the MEE nanobeam is simply supported. The length l, thickness h, and nonlocal parameter e 0 a of the beam are chosen as 10 nm, 0.4 nm, and 4 nm, respectively; the dimensionless load is s = q 0 sin(πx) with q 0 = 0.69; either the electric or magnetic load is assumed to be zero; and the coefficients of the Pasternak foundation are not taken into account in this example. It is pointed out that for the present nanobeam, the number of the test points is taken as n t = 101, and that the number of points on the boundaries is taken as n b = 2, which are assumed to be located at both ends of the beam. Obviously, as shown in Table 1, the present numerical method possesses a high level of accuracy for both the rotation ψ x and the normalized transverse displacement w (i.e., deflection). Also, both increasing the number of the inside points n i and increasing the order of polynomial basis functions ξ suitably can effectively improve the computational accuracy. Moreover, in order to verify the validity of the method for solving the present MEE nanoplate bending problem, the analytical solution for the MEE nanoplate bending under the simply supported boundary and sinusoidal-distribution of the transverse load is additionally derived, which is, for brevity, given in Appendix A. The effects of the supplementary parameter, the number of the order of polynomial basis functions, and the number of the collocation points on the computational accuracy, stability, convergence, and efficiency are evaluated in Figs. 2-4 and Table 2 in detail, where the size and nonlocal parameter of the plate are, respectively, chosen as l = b = 40 nm, h = 0.34 nm, and e 0 a = 1 nm; the dimensionless spring and shear coefficients of the Pasternak foundation are, respectively, taken as k w = 3 and k g = 3 × 10 −2 ; the applied electric and magnetic potentials are adopted as ϕ 0 = 2 × 10 −4 V and φ 0 = 2 × 10 −6 A; and the dimensionless transverse load here is assumed as s = q 0 sin(πx) sin(πy) and with q 0 = 0.69 (see Ref. [10]). Moreover, it is pointed out here that the number of the test points is determined by our numerical tests as n t = 961, which is enough to ensure the convergence in the following numerical examples.  Figure 2 shows the effects of the supplementary parameter χ on the computational accuracy and convergence as n i = 361, n b = 76, and ξ = 14. It can be seen that, with the increase in χ, the three different errors all monotonically decrease, and the computational results converge when χ 250. Thus, for the sake of conservatism, χ = 500 is chosen in the following examples. Moreover, in the case of both the order of the polynomial basis functions and the number of the collocation points being similar to those in Ref. [20], where the bending of an elastic Winkler plate with a simply supported boundary was simulated by the polynomial particular solution method, the accuracy of the present method for simulating the bending of the MEE nanoplate can reach a similar level or even better level.
Besides, Fig. 2 implies that for the three different errors, the accuracy for ψ x and ψ y is practically the same. The reason for this is perhaps that due to the uniform distribution of the collocation points and the symmetry of the structure, boundary, and loads. The analytical solutions and numerical results have the properties of ψ x (x, y) = ψ y (y, x) and ψ x (x, y) = ψ y (y, x), respectively. Therefore, in the following discussion, the error analyses for ψ y are omitted. Figure 3 displays the effects of the order of the polynomial basis functions on the computational stability and accuracy when n i = 361 and n b = 76. As shown in Fig. 3(a), the condition number of the interpolation matrix immoderately increases with the increase in ξ as the multiple-scale technique is not used. Therefore, the proposed method cannot acquire the acceptable numerical solutions without the multiple-scale technique when ξ is large [20,25] . Also, it is found that the condition number keeps moderate as ξ 10, when the multiple-scale technique is used. In fact, the E RMS with the multiple-scale technique is obviously less than that without the multiple-scale technique (see Fig. 3(b)). It is further found from Fig. 3(b) that the E RMS is minimum with ξ = 14. Thus, ξ = 14 is then used in the following examples.  Table 2 The effects of the number of interior and boundary collocation points on the computational accuracy and efficiency  Table 2 reveals the effects of the number of the collocation points on the computational accuracy and efficiency, where both the interior points and the boundary points are always uniformly distributed. It is noted that the present method is highly accurate, and the computational errors converge fast with the increase in either n i or n b . However, when the number of the collocation points is large (such as n i > 361 and n b > 76), there is a slight increase in all the errors. Besides, the required computation time, carried out on the MATLAB (2017) platform in OS Window 10 (64 bit) with 2.80 GHz CPU and 8 GB memory, remarkably increases with the increase in n i and n b . By considering the computational accuracy and efficiency synthetically, n i = 361 and n b = 76 are used in the following examples. Figure 4 displays the distribution of the absolute errors for ψ x and w by surface plots which are based on 2 601 test points. From Fig. 4(a), the maximum and minimum of the absolute errors for ψ x are 3.00 × 10 −5 and 1.18 × 10 −10 , respectively, and the maximum appears at the four corners of the plate. On the other hand, Fig. 4(b) shows that the maximum and minimum of the absolute errors for w are 7.49 × 10 −6 and 2.39 × 10 −6 , respectively, and the location of the maximum error is in the center of the plate. It is remarked that Fig. 4 further demonstrates the high accuracy of the present method.

Further analyses for the rectangle MEE nanoplate bending
Adopting the computational parameters determined above, the effects of the nonlocal parameter, Pasternak foundation, applied electromagnetic loads, and boundary conditions on the bending for the rectangle MEE nanoplate with same size are further investigated in this part, where the applied transverse mechanical load on the plate is the same as before as well.
In fact, for the case of the simply supported boundary condition on all four edges, there is symmetry with respect to ψ x and ψ y . Therefore, the maximum absolute values of ψ x and ψ y are always equal, and max(|ψ y |) is then left out in Figs. 5 and 6. As shown in these two figures, both max(|ψ x |) and max(|w|) monotonously increase with the increase in the nonlocal parameter, which implies that the absolute values of the deflection and rotation predicted by the nonlocal theory are always larger than those predicted by the classical theory (i.e., when e 0 a = 0 nm). The similar conclusions can be found in the studies about buckling problems of the multilayered MEE nanoplates [12] and composite nanoplates with coated one-dimensional quasicrystal [39] , i.e., the increase in nonlocal parameter will always weaken the flexural rigidity of nanoplates. It perhaps can verify the validity of the present results to some extent.
From the first three rows of Table 3, where the structures are symmetric with respect to both x = 0.5 and y = 0.5, it is further shown that the variation of ψ x is symmetric about y = 0.5 and antisymmetric about x = 0.5, whilst the variation of ψ y is respectively symmetric and antisymmetric with respect to x = 0.5 and y = 0.5. On the other hand, the variations of w are symmetric about both x = 0.5 and y = 0.5, and the minimum of w is always located in the center of the plate. It is also noted that for the S-S-S-S boundary condition, ψ x and ψ y are symmetric about y = x (i.e., ψ x (x, y) = ψ y (y, x)). Moreover, when the boundary condition orderly varies to F-S-F-S, S-S-S-S, and C-S-C-S, the maximum absolute values of any one of ψ y and w decrease in turn, which should result from the enhancement of the constraints of the boundary conditions. In fact, the similar phenomena have been found in the related bending analysis of the elastic Mindlin plate [34] . Finally, from the fourth row of Table 3, it is simultaneously observed that all of the antisymmetry of ψ x about x = 0.5, and the symmetries of ψ y and w about x = 0.5 do not appear. Additionally, the maximum absolute values of ψ x , ψ y , and w under the C-S-F-S boundary condition are all larger than the corresponding ones under the C-S-C-S boundary condition, and their locations generally are close to the boundary x = 1.

Bending analysis of the MEE nanoplates with other shapes
In order to demonstrate the general applicability of the present method, the bending problems of the MEE nanoplates with other shapes under the clamped boundary condition and uniform distribution of the transverse load (s = q 0 = 0.69) are investigated in this section. Meanwhile, the nonlocal parameter, coefficients of the Pasternak medium, and applied electric and magnetic potentials are the same as in Subsection 5.1.
First, we consider a circular MEE nanoplate as shown in Fig. 7(a). In the numerical computation, χ = 500, ξ = 16, n i = 1 000, n b = 100, and n t = 2 200 are chosen. The variations of w, ψ x , and ψ y are shown in Figs. 7(b)-7(d), respectively. Obviously, w is symmetric about x = 0.5 and y = 0.5. Also, it is noted that the minimum of w is −0.008 2 and appears at (0.5, 0.5). ψ x is symmetric about y = 0.5 and antisymmetric about x = 0.5. The maximum of ψ x is 0.025 1, which appears at (0.2, 0.5). Besides, it is found that ψ x (x, y) = ψ y (y, x) holds true. Then, another irregular MEE nanoplate as shown in Fig. 8(a) is considered. Numerical results are plotted in Figs. 8(b)-8(d). In order to ensure the convergence of the numerical solutions, χ = 1 000, ξ = 20, n i = 1 331, n b = 140, and n t = 2 276 are adopted. Figure 8(b) shows that w is symmetric about y = x, and its minimum is −0.007 6 and appears at (0.45, 0.45). Figure 8(c) indicates that the maximum and minimum of ψ x are 0.026 6 and −0.021 1, respectively, which appear at (0.16, 0.48) and (0.72, 0.52), respectively. Moreover, the relation of ψ x (x, y) = ψ y (y, x) holds true as well.

Conclusions
In this work, based on the nonlocal theory and Mindlin plate theory, the bending model of MEE nanoplates resting on the Pasternak elastic foundation is established. The method of polynomial particular solutions is developed to solve the driven governing equation system. From the numerical examples, some conclusions can be drawn. With the suitable computational parameters and introducing the multiple-scale technique, the present numerical method possesses superior accuracy, efficiency, and stability to simulate the bending problem of the MEE nanoplates. The present method is also flexible and can be generally applied to carry out the bending analysis of the MEE nanoplates with different structural shapes, boundary conditions, and mechanical loads. Compared with the classical continuum mechanical theory, the presence of nonlocal effect can increase the absolute values of the deflection and rotation. In addition, for the considered model, the increase in either electric potential or the coefficients of the Pasternak medium and the decrease in magnetic potential all lead to increasing in the absolute values of the deflection and rotation.