**About the Editors**

**Francesco Tornabene** is a researcher at the School of Engineering, Department of Innovation Engineering, University of Salento. He was born on 13 January 1978 in Bologna, where he received the high school degree at Liceo Classico San Luigi, in 1997. In 2001 he achieved a National Patent Bologna (Italy) for the Industrial Invention: Friction Clutch for High Performance Vehicles Question BO2001A00442. He received from the University of Bologna—Alma Mater Studiorum, a M.Sc. degree in Mechanical Engineering (Curriculum in Structural Mechanics) on 23/07/2003, discussing a thesis entitled: "Dynamic Behavior of Cylindrical Shells: Formulation and Solution". In December 2003, he was admitted at the PhD course in Structural Mechanics, at the University of Bologna; in 2004, he received from the University of Bologna a Thesis prize in memory of Carlo Felice Jodi; in 2007, he received the Ph.D. degree in Structural Mechanics at the University of Bologna, discussing the Thesis entitled "Modeling and Solution of Shell Structures Made of Anisotropic Materials". From 2007 to 2009 he received a research fellowship by the University of Bologna, working on the "Unified Formulation of Shell Structures Made of Anisotropic Materials. Numerical Analysis Using the Generalized Differential Quadrature Method and the Finite Element Method". From 2011 to 2012, he became a junior researcher within the research program entitled "Advanced Numerical Schemes for Anisotropic Materials"; from 2012 to 2018, he was an Assistant Professor and Lecturer at the Alma Mater Studiorum—University of Bologna; from 2018 to date he is an Aggregate Professor in Structural Mechanics and Lecturer at the University of Salento, Department of Innovation Engineering (Lecce). For a long time, his scientific interests include structural mechanics, solid mechanics, innovative and smart materials, computational mechanics and numerical techniques, and damage and fracture mechanics. He is author of more than 280 scientific publications, and collaborates with many national or international researchers and professors all over the world, as visible from his scientific production. He is author of 11 books; see, for example, "Meccanica delle Strutture a Guscio in Materiale Composito. Il metodo Generalizzato di Quadratura Differenziale" (2012); "Mechanics of Laminated Composite Doubly-Curved Shell Structures. The Generalized Differential Quadrature Method and the Strong Formulation Finite Element Method" (2014); "Laminated Composite Doubly-Curved Shell Structures I. Differential Geometry. Higher-Order Structural Theories" (2016); "Laminated Composite Doubly-Curved Shell Structures II. Differential and Integral Quadrature. Strong Formulation Finite Element Method" (2016); "Anisotropic Doubly-Curved Shells. Higher-Order Strong and Weak Formulations for Arbitrarily Shaped Shell Structures" (2018); among many others. He is member of the Editorial Board for 44 international journals. He is also Editor-in-Chief for two international journals: *Curved and Layered Structures and Journal of Composites Science*; furthermore, he is Associate Editor for seven international journals. In the last few years, he has received different important awards, for example, "Highly Cited Researcher by Clarivate Analytics" (years 2018, 2019, 2020), "Ambassador of Bologna Award for the organization of 21st International Conference on Composite Structures ICCS21, 4-7 September 2018, Bologna, Italy" (2019), and "Member of the European Academy of Sciences" (since 2018). He collaborates as reviewer with more than 260 prestigious international journals in the structural mechanics field. From 2012, his teaching activity has included dynamics of structures; computational mechanics; plates and shells; theory of structures; structural mechanics. He is habilitated as Associate Professor and Full Professor in the area 08/B2 (Structural Mechanics) and as Associate Professor in Area 09/A1 (aeronautical and aerospace engineering and naval architecture).

**Rossana Dimitri** is an Associate Professor at the School of Engineering, Department of Innovation Engineering, University of Salento, Lecce, Italy (since 2019). She received from the University of Salento, a M. Sc. degree in "Materials Engineering" in 2004, a Ph.D. degree in "Materials and Structural Engineering" in 2009, and a Ph.D. degree in "Industrial and Mechanical Engineering" in 2013. In 2005, she received from the University of Salento the "Best M. Sc. Thesis Prize 2003–2004" in memory of Eng. Gabriele De Angelis; in 2013, she was awarded by the Italian Group for Computational Mechanics (GIMC) for the Italian selection of the 2013 ECCOMAS PhD Award. Her current interests include structural mechanics, solid mechanics, damage and fracture mechanics, contact mechanics, isogeometric analysis, high-performing computational methods, and consulting in applied technologies and technology transfer. During 2010 and 2011, she received a research fellowship by ENEA Research Centre of Brindisi (UTTMATB-COMP) for the development and the characterization of some thermoplastic composites for thermal solar panels and adhesively bonded turbine blades under severe environmental conditions. During 2011 and 2012, she was a visiting scientist with a fellowship at the Institut fur Kontinuumsmechanik Gottfried Wilhelm ¨ Leibniz Universitat Hannover to study interfacial problems with isogeometric approaches. From ¨ 2013 to 2016, she was a researcher at the University of Salento, within the ERC starting research grant "INTERFACES" on "Computational mechanical modelling of structural interfaces based on isogeometric approaches". From 2013 to 2019, she has collaborated as a researcher RTD-B with the University of Bologna and the Texas A&M University for a comparative assessment of some advanced numerical collocation methods with lower computational cost for fracturing problems and structural modelling of composite plates and shells, made by isotropic, orthotropic and anisotropic materials. She is author of 112 scientific publications, and she collaborates with many national or international researchers and professors worldwide, as visible from her scientific production. She also collaborates with different prestigious international journals in the structural mechanics field, as reviewer, member of the editorial board, and guest editor for different Special Issues.

### *Editorial* **Special Issue on Advanced Theoretical and Computational Methods for Complex Materials and Structures**

**Francesco Tornabene \* and Rossana Dimitri \***

Department of Innovation Engineering, Università del Salento, Via per Monteroni, 73100 Lecce, Italy

**\*** Correspondence: francesco.tornabene@unisalento.it (F.T.); rossana.dimitri@unisalento.it (R.D.)

#### **1. Introduction**

The large use of composite materials and shell structural members with complex geometries in technologies related to various branches of engineering, has gained increased attention from scientists and engineers for the development of even more refined approaches, to investigate their mechanical behavior. It is well known that composite materials are able to provide higher values of strength stiffness, and thermal properties, together with conferring reduced weight, which can affect the mechanical behavior of beams, plates, and shells, in terms of static response, vibrations, and buckling loads. At the same time, enhanced structures made of composite materials can feature internal length scales and non-local behaviors, with great sensitivity to different staking sequences, ply orientations, agglomeration of nanoparticles, volume fractions of constituents, and porosity levels, among others. In addition to fiber-reinforced composites and laminates, increased attention has been paid in literature to the study of innovative components such as functionally graded materials (FGMs), carbon nanotubes (CNTs), graphene nanoplatelets, and smart constituents. Some examples of smart applications involve large stroke smart actuators, piezoelectric sensors, shape memory alloys, magnetostrictive and electrostrictive materials, as well as auxetic components and angle-tow laminates. These constituents can be included in the lamination schemes of smart structures to control and monitor the vibrational behavior or the static deflection of several composites. The development of advanced theoretical and computational models for composite materials and structures is a subject of active research and this is explored here for different complex systems, including their static, dynamic, and buckling responses; fracture mechanics at different scales; as well as the adhesion, cohesion, and delamination of materials and interfaces.

#### **2. Enhanced Theoretical and Computational Formulations**

In a context where an increased theoretical/computational demand is usually required to solve solid mechanics problems, this Special Issue has collected 11 papers on the application of high-performing computational strategies and enhanced theoretical formulations to solve different static and/or dynamic problems, for different engineering applications also in coupled conditions. A wide variety of examples and topics is considered, with highly useful insights both from a scientific and design perspective. More specifically, the first paper, authored by M.H. Jalaei, R. Dimitri and F. Tornabene, combines the Hamilton's variational principle and the Eringen's constitutive model to study the dynamic stability of orthotropic temperature-dependent single-layered graphene sheets embedded in a temperature-dependent elastomeric medium and subjected to a biaxial oscillating loading in thermal environment [1]. These nanostructures are largely adopted as important components in various highly technological industries, such as nanoactuators, nanoresonators, nanosensors, and nanocomposites, such that non-classical continuum approaches are usually required for an accurate analysis of results, especially for complicated coupled problems.

**Citation:** Tornabene, F.; Dimitri, R. Special Issue on Advanced Theoretical and Computational Methods for Complex Materials and Structures. *Appl. Sci.* **2021**, *11*, 2532. https://doi.org/10.3390/ app11062532

Received: 6 March 2021 Accepted: 10 March 2021 Published: 12 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

In the second work by K.J. Huang et al. [2], the authors develop an efficient theoretical approach to investigate the nonlocal and size-dependent (NLSD) effects on the dielectric response of plasmonic nanostructures that incorporate the spatially local, nonlocal, size and even analogous quantum-size responses of the material. On the other hand, unlike traditional engineering structural problems, the design of micro-electromechanical systems (MEMS) usually involves microstructures, novel materials, and extreme operating conditions, where multi-source uncertainties usually exist. In such a context, the work by Z. Huang et al. [3] develops an evidence–theory-based robustness optimization method for a robust design of MEMS, including a micro-force sensor, an image sensor, and a capacitive accelerometer of practical interest. Another application is represented by microfluidic devices, for which an analytical and practical method is efficiently proposed in the work by U. Tuzun [4] to study the Brownian motion in nanoparticle suspensions, as used in a variety of applications involving the transport of reagents and products to and from structured material surfaces. Among geotechnical engineering applications, the work authored by H. Sun and W. Sun [5] provides a useful finite element study for an effective measurement of safety and serviceability of existing metro tunnels during adjacent excavation, in lieu of centrifuge model tests, analytical or semi-analytical methods, or a more expensive in-situ monitoring. In the further work by Y. Zhu et al. [6], a quantitative evaluation of the solid displacement induced by a ground loss and shield machine mechanical effect is provided in metro tunnel constructions, based on an elastic half-space Mindlin model. The proportion of stratum displacement caused by each factor was quantitatively analyzed in this work, as a valid guidance for the stratum displacement calculation of shield tunneling in the future. A novel method is also proposed by X. Tan, H. Zhang et al. [7] to define the three-dimensional Greenfield stratum movements caused by a shield tunnel construction, where an elastic half-space model of mirror-sink method is combined with a modified analytical approach to give a useful design tool for underground engineering constructions. An efficient numerical model is, instead, established by H. Sun et al. [8] to compute the tunnel deformation caused by a circular excavation, while adopting a hypoplasticity nonlinear constitutive model, able to account for path-and strain-dependent soil stiffnesses even at small strains. This method is useful for practical engineering soil-structure problems. The extensive use of composite materials and structures in many engineering applications with complex microstructure and manufacturing processes, requires a thorough attention to their mechanical performances, such as the structural deflection damage and load capacity [9–12], as well as the buckling and dynamic behavior [13–18], along with possible related uncertainties and stochastic variations. In this setting, the work authored by S. Zhang and X. Chen [19] provides a stochastic natural frequency analysis of typical composite structures with micro-scale (constituent-scale) and meso-scale (ply-scale) uncertainties, based on Monte-Carlo simulations and the response surface method. Another numerical model is proposed by S. Shahbazi et al. [20] to analyze the non-linear time-history of different structural members under the simultaneous effects of far- and near-field earthquakes, as useful for design computational tools. The last work, authored by A. Pavlovic, et al. [21] finally proposes a finite element study of the mechanical behavior of palletized products and polyethylene terephthalate bottles, partially filled with liquid under compressive stress conditions, as commonly occurs during a transportation process. In this last case, an accurate computational prediction of the transport-related integrity risks of bottles is desirable in order to avoid any kind of instability and damage phenomena within commercial products, with a clear economical repercussion.

#### **3. Future Developments**

Although this Special Issue has been closed, further developments on the theoretical and computational modeling of enhanced structures and composite materials are expected, including their static, dynamic, and buckling responses; fracture mechanics at different scales; as well as the adhesion, cohesion, and delamination of materials and interfaces, as useful for many industrial applications.

#### **References**


## *Article* **Dynamic Stability of Temperature-Dependent Graphene Sheet Embedded in an Elastomeric Medium**

#### **Mohammad Hossein Jalaei 1, Rossana Dimitri <sup>2</sup> and Francesco Tornabene 2,\***


Received: 30 January 2019; Accepted: 25 February 2019; Published: 1 March 2019

**Abstract:** This work applies the first-order shear deformation theory (FSDT) to study the dynamic stability of orthotropic temperature-dependent single-layered graphene sheet (SLGS) embedded in a temperature-dependent elastomeric medium and subjected to a biaxial oscillating loading in a thermal environment. Possible thermal effects are considered in the size-dependent governing equations of the problem. These last ones are derived by means of the Hamilton's variational principle combined with the Eringen's differential constitutive model. Navier's solution as well as Bolotin's approach are applied to obtain the dynamic instability region (DIR) of the graphene sheet. Thus, a parametric study is carried out to explore the sensitivity of the DIR of the graphene sheet to the temperature variation, the static load factor, the aspect ratio, the foundation type, and the nonlocal parameter (NP). Results indicate that the dimensionless pulsation frequency reduces for increasing values of temperature and NP, whereas the size effect becomes even more pronounced for increasing temperatures. In addition, the adoption of temperature-dependent mechanical properties, rather than independent ones, yields a global shift of the DIR to smaller pulsating frequencies. This proves the relevance of the temperature-dependent mechanical properties to obtain reliable results, in a physical sense.

**Keywords:** dynamic stability; elastomeric foundation; Eringen's differential constitutive model; graphene sheet; temperature-dependent properties

#### **1. Introduction**

Among nanostructures, graphene sheets have increased the interest of the scientific community, owing to their astonishing thermal, chemical, electrical, and mechanical properties. Due to the distinguishing features, theses nanostructures have been adopted as important components in various high technology industries such as nanoactuators, nanoresonators, nanosensors, and nanocomposites.

In order to investigate the mechanical characteristics of micro/nano-scale structures accurately, the small-scale effects should be considered. It has been recognized that the main properties of materials and structures at micro/nano-scale are size-dependent and differ significantly from their behavior at larger scales. Hence, non-classical continuum mechanics such as couple stress theory (CST), nonlocal strain gradient theory (NSGT), and nonlocal elasticity theory (NET) have been applied by the scientific community to analyze the mechanical behavior of micro/nanomaterials due to the difficulties in experimental investigations at the nanoscale and the high computational costs of molecular dynamics (MD) simulations.

Among these theories, the NET introduced by Eringen [1] has been broadly employed in order to include the size effects on the static and dynamic analyses of nanomaterials. For example, Murmu et al. [2] investigated the in-plane magnetic field effect on the vibration of a magnetically sensitive single-layered graphene sheet (SLGS) resting on an elastic foundation via the NET in conjunction with Kirchhoff plate theory. Wang et al. [3] presented an analytical solution for the nonlinear vibration behavior of viscoelastic double layered nanoplates. The problem was tackled through the Kirchhoff plate theory and the multiple scales method. Karliˇci´c et al. [4] performed the nonlinear vibration and dynamic instability of single-walled carbon nanotubes (SWCNTs) resting on a Kelvin–Voigt viscoelastic medium under magnetic field using NET and the Euler–Bernoulli beam theory (EBT). Arani and Jalaei [5] employed the Fourier–Laplace transformation method as well as the sinusoidal shear deformation theory (SSDT) to examine the longitudinal magnetic field effect on the transient behavior of a viscoelastic simply-supported orthotropic nanoplate embedded on a visco-Pasternak medium. Kolahchi et al. [6] studied the dynamic buckling of embedded viscoelastic laminated nanoplates on the basis of NET in conjunction with refined zigzag theory. Jalaei et al. [7] studied the dynamic stability of a viscoelastic graphene sheet embedded on viscose matrix with various boundary conditions according to the FSDT and NET by means of the Ritz technique and the Bolotin's approach. Most recently, Huang et al. [8] applied an analytical solution to perform the dynamic stability of nanobeams exposed to an axially oscillating loading on the basis of the NET and EBT. Additionally, several works in literature have utilized the NET for the study of buckling, bending, vibration, and wave propagation problems of nanobeams/rods [9–12], nanoshells [13–15], and nanoplates [16–22].

It is recognized that in practical applications, nanostructures or nano reinforced composites often work in environmental conditions. Hence, the mechanical response of these structures at different environmental conditions has been a great subject of investigation among the scientific community. In this regard, Malekzadeh et al. [23] investigated the thermal buckling behavior of arbitrary straight-sided quadrilateral orthotropic nanoplates resting on Pasternak foundation by means of the nonlocal Kirchhoff plate theory and differential quadrature method (DQM). Ansari et al. [24] applied the NET within the framework of the EBT and the Timoshenko beam theory (TBT) to examine the dynamic stability of embedded SWCNTs in a thermal environment. In line with the previous works, Karliˇci´c et al. [25] analyzed the thermal vibration and buckling behavior of the embedded multi-layered graphene sheets (MLGSs) according to the NET and Kirchhoff–Love plate theory. Sobhy [26] applied state-space concept to carry out the static bending of Levy-type nanoplates subjected to hygrothermal and mechanical loads based on the NET in conjunction with two-variable plate theory. A DQM numerical approach was applied by Ebrahimi and Hosseini [27] to investigate the nonlinear vibration behavior of bilayer viscoelastic nanoplates in a thermal environment. Ansari and Gholami [28] studied the dynamic instability of multi-walled carbon nanotubes (MWCNTs) embedded on a Winkler medium including thermal effects, while applying a nonlocal TBT. The governing equations were solved via generalized DQM in conjunction with Bolotin's approach. Wu et al. [29] employed DQM and Bolotin's approach to explore the dynamic instability behavior of functionally graded graphene-reinforced nanocomposite plates under an axially oscillating loading and uniform temperature variation according to the FSDT. Jouneghani et al. [30] employed NET to study the bending characteristics of porous functionally-graded (FG) Timoshenko nanobeams exposed to hygro-thermo-mechanical loads. Further valuable research of nanostructures in a thermal environment can be found in [31–34].

Despite many works in literature have focused on the dynamic stability of nanobeams subjected to a thermal loading, there is a general lack of works focusing on the dynamic stability of nanoplates such as graphene sheets in a thermal environment. In most cases, the study of the thermal sensitivity of mechanical behavior of a SLGS is associated with temperature-independent material properties, whereas the experimental outcomes and MD-based simulations have proved the dependence of nanomaterial properties to temperature.

Hence, we propose a numerical study on the dynamic stability of a temperature-dependent SLGS embedded in a temperature-dependent elastomeric medium and subjected to a thermo-mechanical loading. This represents the key point of the present investigation, where the governing equations of motion are derived according to the FSDT in conjunction with the NET. These equations are then converted into a Mathieu–Hill type equation via the Navier's method, while obtaining the principle unstable region through the application of the Bolotin's method. A large systematic investigation evaluates the influences of the nonlocal parameter (NP), aspect ratio, temperature variation, foundation type, and static load factor on the dynamic instability region (DIR) of an orthotropic SLGS, as useful in design of nano electro-mechanical systems and micro electro-mechanical systems (NEMS/MEMS).

#### **2. Basic Formulation**

A temperature-dependent SLGS is here studied, with length *a*, width *b*, and thickness *h* subjected to a biaxial periodic compressive loading in a thermal environment (Figure 1). The graphene sheet rests on an elastomeric temperature-dependent foundation, which is described by means of the spring constant *kw* and the shear layer parameter *kp*. Among many possibilities of selecting the plate theory and solution methodology, in this work we select a FSDT, while considering possible shear effects on the structural response of the SLGS despite their atomic dimensions, as already demonstrated in [35]. The selected FSDT is here combined to the NET, whose basic notions are provided in the following overview.

**Figure 1.** Configuration of the single-layered graphene sheet (SLGS) resting on an elastomeric foundation and subjected to a biaxial periodic compressive loading in a thermal environment.

#### *2.1. Overview of NET*

The NET introduced by Eringen [1], unlike the classical elasticity theory, states the stress tensor as a nonlocal function of the strain for the whole points in the domain. Hence, for a linear homogenous body, the nonlocal stress tensor at a point *x* is given as:

$$\sigma\_{\mathbf{i}\rangle}{}^{nl}(\mathbf{x}) = \int\_{\mathcal{v}} \mathfrak{a}\left(|\mathbf{x} - \mathbf{x}'| , \mathbf{\tau}\right) \sigma\_{\mathbf{i}\rangle}{}^{l}dV(\mathbf{x}'), \qquad \forall \mathbf{x} \in \ V \tag{1}$$

where *σnl ij* is the nonlocal stress, *<sup>σ</sup><sup>l</sup> ij* is the local stress tensors, *α*(|*x* − *x*- |, *τ*) stands for the nonlocal kernel function. Additionally, *τ* is defined as *τ* = *e*0*a*/*l*, where *a* and *l* represent the internal and external characteristic lengths, respectively, and the constant *e*<sup>0</sup> is introduced for adjusting the model according to experimental data or MD simulations. The differential constitutive form of Equation (1) can be written as follows [1,25], including the thermal effects, namely:

$$\left(1 - \mu \nabla^2\right) \sigma\_{ij}^{nl} = \sigma\_{ij}^l = \mathbb{C}\_{ijkl} \left(\varepsilon\_{kl} - a\_{kl}\Delta T\right), \qquad \qquad \nabla^2 = \partial^2/\partial \mathbf{x}^2 + \partial^2/\partial y^2. \tag{2}$$

*Appl. Sci.* **2019**, *9*, 887

In Equation (2), *μ* = (*e*0*l*) <sup>2</sup> denotes the NP, *<sup>α</sup>ij* is the thermal expansion coefficient, <sup>Δ</sup>*<sup>T</sup>* <sup>=</sup> *<sup>T</sup>* <sup>−</sup> *<sup>T</sup>*<sup>0</sup> is the temperature variation with respect to the reference temperature *T*0. In addition, *εkl* indicates the strain tensor and *Cijkl* describes the fourth-order stiffness tensor.

Based on Equation (2), the constitutive relations based on the NET for the orthotropic nanoplate in thermal environment reads as follows:

$$\begin{Bmatrix} \sigma\_{xx}^{nl} \\ \sigma\_{yy}^{nl} \\ \sigma\_{yz}^{nl} \\ \sigma\_{zx}^{nl} \\ \sigma\_{xy}^{nl} \end{Bmatrix} - \mu \, \nabla^2 \begin{Bmatrix} \sigma\_{xx}^{nl} \\ \sigma\_{yy}^{nl} \\ \sigma\_{yz}^{nl} \\ \sigma\_{zx}^{nl} \\ \sigma\_{xy}^{nl} \end{Bmatrix} = \begin{bmatrix} \mathbb{C}\_{11} & \mathbb{C}\_{12} & 0 & 0 & 0 \\ \mathbb{C}\_{21} & \mathbb{C}\_{22} & 0 & 0 & 0 \\ 0 & 0 & \mathbb{C}\_{44} & 0 & 0 \\ 0 & 0 & 0 & \mathbb{C}\_{55} & 0 \\ 0 & 0 & 0 & 0 & \mathbb{C}\_{66} \end{bmatrix} \begin{Bmatrix} \varepsilon\_{xx} - a\_{xx}\Delta T \\ \varepsilon\_{yy} - a\_{yy}\Delta T \\ \gamma\_{yz} \\ \gamma\_{xz} \\ \gamma\_{xy} \end{Bmatrix} . \tag{3}$$

in which the elastic constants of orthotropic materials can be determined by

$$\mathbb{C}\_{11} = \frac{\mathbb{E}\_1(T)}{1 - \mathbb{V}\_{12}\mathbb{V}\_{21}}, \quad \mathbb{C}\_{12} = \frac{\mathbb{V}\_{12}\mathbb{E}\_2(T)}{1 - \mathbb{V}\_{12}\mathbb{V}\_{21}}, \quad \mathbb{C}\_{22} = \frac{\mathbb{E}\_2(T)}{1 - \mathbb{V}\_{12}\mathbb{V}\_{21}}, \quad \mathbb{C}\_{44} = \mathbb{G}\_{23}(T), \quad \mathbb{C}\_{55} = \mathbb{G}\_{13}(T), \quad \mathbb{C}\_{66} = \mathbb{G}\_{12}(T). \tag{4}$$

#### *2.2. Strain Displacement Relationships*

In what follows we apply the FSDT to formulate the governing equations of the problem, while including the thickness shear deformations and rotary effects. Thus, the displacement components associated to mid-surface (*u*0, *v*0, *w*0), transverse normal rotations (*φx*, *φy*), and the displacement of an arbitrary point (*u*, *v*, *w*) are related as [36]

$$\begin{aligned} u(x, y, z, t) &= u\_0(x, y, t) + z \phi\_x(x, y, t), \\ v(x, y, z, t) &= v\_0(x, y, t) + z \phi\_y(x, y, t), \\ w(x, y, z, t) &= w\_0(x, y, t), \end{aligned} \tag{5}$$

in which *t* was the time variable. The nonzero normal and transverse shear strains of the FSDT can be derived by

$$\left\{ \begin{array}{c} \varepsilon\_{xx} \\ \varepsilon\_{yy} \\ \gamma\_{xy} \end{array} \right\} = \left\{ \begin{array}{c} \varepsilon\_{xx}^{(0)} \\ \varepsilon\_{yy}^{(0)} \\ \gamma\_{xy}^{(0)} \end{array} \right\} + z \left\{ \begin{array}{c} \varepsilon\_{xx}^{(1)} \\ \varepsilon\_{yy}^{(1)} \\ \gamma\_{xy}^{(1)} \end{array} \right\}, \qquad \left\{ \begin{array}{c} \gamma\_{yz} \\ \gamma\_{xz} \end{array} \right\} = \left\{ \begin{array}{c} \gamma\_{yz}^{(0)} \\ \gamma\_{xz}^{(0)} \end{array} \right\}, \tag{6}$$

where

$$\begin{Bmatrix} \varepsilon\_{xx}^{(0)} \\ \varepsilon\_{yy}^{(0)} \\ \gamma\_{xy}^{(0)} \end{Bmatrix} = \left\{ \begin{array}{c} \frac{\partial u\_{0}}{\partial x} \\ \frac{\partial v\_{0}}{\partial y} \\ \frac{\partial w\_{0}}{\partial y} + \frac{\partial w\_{0}}{\partial x} \end{array} \right\}, \quad \left\{ \begin{array}{c} \varepsilon\_{xx}^{(1)} \\ \varepsilon\_{yy}^{(1)} \\ \gamma\_{xy}^{(1)} \end{array} \right\} = \left\{ \begin{array}{c} \frac{\partial \phi\_{x}}{\partial x} \\ \frac{\partial \phi\_{y}}{\partial y} \\ \frac{\partial \phi\_{z}}{\partial y} + \frac{\partial \phi\_{y}}{\partial x} \end{array} \right\}, \quad \left\{ \begin{array}{c} \gamma\_{yz}^{(0)} \\ \gamma\_{xz}^{(0)} \end{array} \right\} = \left\{ \begin{array}{c} \Phi\_{y} + \frac{\partial w\_{0}}{\partial y} \\ \Phi\_{x} + \frac{\partial w\_{0}}{\partial x} \end{array} \right\}. \tag{7}$$

#### **3. Energy Method**

#### *3.1. Kinetic Energy*

The kinetic energy for the nanoplate can be computed as

$$\begin{split} \mathcal{K} &= \,\_2^0 \int\_A \int\_{-h/2}^{h/2} \left( \left(\frac{\partial u}{\partial l}\right)^2 + \left(\frac{\partial v}{\partial l}\right)^2 + \left(\frac{\partial w}{\partial l}\right)^2 \right) \mathrm{d}A \mathrm{d}z \\ &= \,\_2^1 \int\_A \left\{ \log \left( \left(\frac{\partial u\_0}{\partial l}\right)^2 + \left(\frac{\partial v\_0}{\partial l}\right)^2 + \left(\frac{\partial v\_0}{\partial l}\right)^2 \right) + 2I\_1 \left( \left(\frac{\partial u\_0}{\partial l} \frac{\partial p\_1}{\partial l}\right) + \left(\frac{\partial v\_0}{\partial l} \frac{\partial p\_2}{\partial l}\right) \right) + I\_2 \left( \left(\frac{\partial q\_1}{\partial l}\right)^2 + \left(\frac{\partial p\_3}{\partial l}\right)^2 \right) \right\} \mathrm{d}A \end{split} \tag{8}$$

where *ρ* denotes the density of the orthotropic graphene sheet and *Ii* refers to the mass inertias defined by

$$I\_i = \rho \int\_{-h/2}^{h/2} z^i \mathbf{d}z \quad (i = 0, 1, 2) \,\text{.}\tag{9}$$

#### *3.2. Strain Energy*

The strain energy for the nanoplate could be obtained as

$$\mathrm{d}I = \frac{1}{2} \int\_{A} \int\_{-h/2}^{h/2} (\sigma\_{\mathrm{xx}}\varepsilon\_{\mathrm{xx}} + \sigma\_{yy}\varepsilon\_{yy} + K\_{\mathrm{s}}\sigma\_{\mathrm{xz}}\gamma\_{\mathrm{xz}} + K\_{\mathrm{s}}\sigma\_{yz}\gamma\_{yz} + \sigma\_{\mathrm{xy}}\gamma\_{\mathrm{xy}}) \mathrm{d}A \,\mathrm{d}z. \tag{10}$$

In the above equation, *Ks* stands for the shear correction factor that is here set to 5/6, as typically assumed in the literature for nanoplates and single-layer graphene sheets [20].

The combination of Equations (6), (7), and (10) yields

$$\begin{split} \mathcal{U}I &= \frac{1}{2} \int\_{A} \Bigg( N\_{\text{xx}} \frac{\partial \mathbf{u}\_{0}}{\partial \mathbf{x}} + M\_{\text{xx}} \frac{\partial \boldsymbol{\phi}\_{\text{x}}}{\partial \mathbf{x}} + N\_{\text{yy}} \frac{\partial \mathbf{v}\_{0}}{\partial \mathbf{y}} + M\_{\text{yy}} \frac{\partial \boldsymbol{\phi}\_{\text{y}}}{\partial \mathbf{y}} + N\_{\text{xy}} \Big( \frac{\partial \mathbf{u}\_{0}}{\partial \mathbf{y}} + \frac{\partial \mathbf{v}\_{0}}{\partial \mathbf{x}} \Big) \\ &+ M\_{\text{xy}} \Big( \frac{\partial \boldsymbol{\phi}\_{\text{x}}}{\partial \mathbf{y}} + \frac{\partial \boldsymbol{\phi}\_{\text{y}}}{\partial \mathbf{x}} \Big) + Q\_{\text{xz}} \Big( \frac{\partial \mathbf{w}\_{0}}{\partial \mathbf{x}} + \boldsymbol{\phi}\_{\text{x}} \Big) + Q\_{\text{yz}} \Big( \frac{\partial \mathbf{w}\_{0}}{\partial \mathbf{y}} + \boldsymbol{\phi}\_{\text{y}} \Big) \Big) dA \end{split} \tag{11}$$

where

$$\left\{ \begin{array}{c} N\_{a\beta} \\ M\_{a\beta} \end{array} \right\} = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \sigma\_{a\beta}^{nl} \left\{ \begin{array}{c} 1 \\ z \end{array} \right\} \,\mathrm{d}z \qquad a = \text{x, y} \quad \beta = \text{x, y} \tag{12}$$

$$Q\_{az} = K\_s \int\_{-\frac{b}{2}}^{\frac{b}{2}} \sigma\_{az}^{nl} \, \mathrm{d}z \qquad a = x, \, y \, . \tag{13}$$

*Nαβ* denotes the resultant force, *Mαβ* describes the resultant moment, and *Qα<sup>z</sup>* indicates the resultant shear force.

By substitution of Equations (3), (6), and (7) into Equations (12) and (13), the stress resultantdisplacement relations are determined as

$$\begin{array}{l} N\_{\text{xx}} - \mu \nabla^2 (N\_{\text{xx}}) = A\_{11} \frac{\partial u\_0}{\partial x} + A\_{12} \frac{\partial v\_0}{\partial y} - N\_{\text{xx}}^T, \\ N\_{\text{yy}} - \mu \nabla^2 (N\_{\text{yy}}) = A\_{12} \frac{\partial u\_0}{\partial x} + A\_{22} \frac{\partial v\_0}{\partial y} - N\_{\text{yy}}^T, \\ N\_{\text{xy}} - \mu \nabla^2 (N\_{\text{xy}}) = A\_{66} \left( \frac{\partial u\_0}{\partial y} + \frac{\partial v\_0}{\partial x} \right), \end{array} \tag{14}$$

$$\begin{cases} M\_{xx} - \mu \nabla^2 (M\_{xx}) = D\_{11} \frac{\partial \phi\_x}{\partial x} + D\_{12} \frac{\partial \phi\_y}{\partial y} - M\_{xx}^T, \\ M\_{yy} - \mu \nabla^2 (M\_{yy}) = D\_{12} \frac{\partial \phi\_x}{\partial x} + D\_{22} \frac{\partial \phi\_y}{\partial y} - M\_{yy}^T, \\ M\_{xy} - \mu \nabla^2 (M\_{xy}) = D\_{66} \left( \frac{\partial \phi\_x}{\partial y} + \frac{\partial \phi\_y}{\partial x} \right), \end{cases} \tag{15}$$

$$\begin{cases} Q\_{yz} - \mu \nabla^2 \left( Q\_{yz} \right) = J\_{44} \left( \phi\_y + \frac{\partial w\_0}{\partial y} \right), \\ Q\_{xz} - \mu \nabla^2 \left( Q\_{xz} \right) = J\_{55} \left< \phi\_x + \frac{\partial w\_0}{\partial x} \right>. \end{cases} \tag{16}$$

More specifically, the extensional stiffness *Aij*, the bending stiffness *Dij*, and the shear stiffness *Jii* of the SLGS, are given as

$$\left(A\_{\vec{\text{ij}}\prime}, D\_{\vec{\text{ij}}\prime}\right) = \int\_{-h/2}^{h/2} \mathbb{C}\_{\vec{\text{ij}}} \left(\mathbf{1}, z^2\right) \mathrm{d}z \qquad \qquad \qquad \left(i, j = 1, 2, 6\right), \tag{17}$$

$$\mathbf{J}\_{\rm ii} = \mathbf{K}\_{\rm s} \int\_{-h/2}^{h/2} \mathbf{C}\_{\rm ii} \, \mathbf{d}z \, \, = \, \mathbf{K}\_{\rm s} h \, \mathbf{C}\_{\rm ii} \qquad \left( i = \mathbf{4}, \mathbf{5} \right) \, . \tag{18}$$

Additionally, *N<sup>T</sup> xx*, *N<sup>T</sup> yy* and *M<sup>T</sup> xx*, *M<sup>T</sup> yy* denote the thermally induced forces and moments due to a uniform temperature variation, respectively, determined by

$$\left\{ \begin{array}{c} \text{N}\_{xx}^{T} \\ \text{N}\_{yy}^{T} \end{array} \right\} = \int\_{-\frac{1}{2}}^{\frac{1}{2}} \left\{ \begin{array}{c} A\_{11}a\_{11} + A\_{12}a\_{22} \\ A\_{12}a\_{11} + A\_{22}a\_{22} \end{array} \right\} \Delta T \text{ d}z,\tag{19}$$

*Appl. Sci.* **2019**, *9*, 887

$$\left\{ \begin{array}{c} M\_{xx}^{T} \\ M\_{yy}^{T} \end{array} \right\} = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \left\{ \begin{array}{c} A\_{11}a\_{11} + A\_{12}a\_{22} \\ A\_{12}a\_{11} + A\_{22}a\_{22} \end{array} \right\} \Delta T \, z \, \mathrm{d}z. \tag{20}$$

*3.3. External Work*

The external work done due to the temperature-dependent elastomeric medium and axially oscillating loading could be expressed as

$$V = \frac{1}{2} \int\_{A} \left( -P\_{\text{xx}} \frac{\partial^2 w\_0}{\partial \mathbf{x}^2} - P\_{yy} \frac{\partial^2 w\_0}{\partial y^2} - k\_{\text{uv}} w\_0 + k\_p \left( \frac{\partial^2 w\_0}{\partial \mathbf{x}^2} + \frac{\partial^2 w\_0}{\partial y^2} \right) \right) w\_0 \, \text{d}A \,, \tag{21}$$

where *Pxx* = *Pyy* = *P* denotes the axial dynamic forces. Additionally, the elastomeric foundation stiffness *kw* could be obtained as [37]

$$k\_{\mathcal{W}} = \frac{E\_0}{4a\left(1 - \nu\_0^2\right)\left(2 - c\_1\right)^2} \left[5 - \left(2\gamma\_1^2 + 6\gamma\_1 + 5\right)\exp\left(-2\gamma\_1\right)\right],\tag{22}$$

where

$$\alpha\_1 = (\gamma\_1 + 2) \exp(-\gamma\_1) \,, \quad \gamma\_1 = \frac{H\_\mathfrak{s}}{a} \,, \quad E\_0 = \frac{E\_\mathfrak{s}}{(1 - \nu\_\mathfrak{s}^2)} \,, \quad \nu\_0 = \frac{\nu\_\mathfrak{s}}{1 - \nu\_\mathfrak{s}} \, , \tag{23}$$

while *Es*, *ν<sup>s</sup>* and *Hs* stand for the Young's modulus, Poisson's ratio, and depth of the foundation, respectively. In this research, *Es* is supposed to be temperature-dependent, whereas *ν<sup>s</sup>* is kept constant.

#### **4. Motion Equations**

Hamilton's variational principle has been utilized in order to determine the equilibrium equations of nanoplate. Accordingly, one may write

$$\int\_{0}^{t} (\delta \mathcal{U} + \delta V - \delta \mathcal{K}) dt = 0 \,. \tag{24}$$

Inserting Equations (8), (11) and (21) into Equation (24) yields the following equations of motion

$$
\delta\mu\_0 \quad ; \qquad \frac{\partial N\_{xx}}{\partial x} + \frac{\partial N\_{xy}}{\partial y} = I\_0 \frac{\partial^2 \mu\_0}{\partial t^2} + I\_1 \frac{\partial^2 \phi\_x}{\partial t^2} \,, \tag{25}
$$

$$
\delta v\_0 \quad ; \qquad \frac{\partial N\_{xy}}{\partial \mathbf{x}} + \frac{\partial N\_{yy}}{\partial y} = I\_0 \frac{\partial^2 v\_0}{\partial t^2} + I\_1 \frac{\partial^2 \phi\_y}{\partial t^2} \,, \tag{26}
$$

$$\delta w\_{0}: \quad \frac{\partial Q\_{\rm diff}}{\partial x} + \frac{\partial Q\_{\rm gr}}{\partial y} - \left(P\_{\rm xx} + N\_{\rm xx}^{T}\right)\frac{\partial^{2}w\_{0}}{\partial x^{2}} - \left(P\_{\rm yy} + N\_{\rm yy}^{T}\right)\frac{\partial^{2}w\_{0}}{\partial y^{2}} - k\_{\rm W}w\_{0} + k\_{\rm p}\left(\frac{\partial^{2}w\_{0}}{\partial x^{2}} + \frac{\partial^{2}w\_{0}}{\partial y^{2}}\right) = I\_{0}\frac{\partial^{2}w\_{0}}{\partial t^{2}},\tag{27}$$

$$
\delta\phi\_{\mathbf{x}} : \quad \frac{\partial M\_{\mathbf{x}\mathbf{x}}}{\partial \mathbf{x}} + \frac{\partial M\_{\mathbf{x}\mathbf{y}}}{\partial y} - Q\_{\mathbf{x}\mathbf{z}} = I\_1 \frac{\partial^2 u\_0}{\partial t^2} + I\_2 \frac{\partial^2 \phi\_{\mathbf{x}}}{\partial t^2} \tag{28}
$$

$$
\delta\phi\_{\mathcal{Y}} : \quad \frac{\partial M\_{xy}}{\partial \mathbf{x}} + \frac{\partial M\_{yy}}{\partial y} - Q\_{\mathcal{Y}^2} = I\_1 \frac{\partial^2 v\_0}{\partial t^2} + I\_2 \frac{\partial^2 \phi\_{\mathcal{Y}}}{\partial t^2} \,. \tag{29}
$$

Finally, by substituting Equations (14)–(16) into Equations (25)–(29), the nonlocal governing equations of orthotropic graphene sheet resting on an elastomeric medium including the thermal effect can be determined as follows:

$$A\_{11} \left(\frac{\partial^2 u\_0}{\partial \mathbf{x}^2}\right) + (A\_{12} + A\_{66}) \left(\frac{\partial^2 v\_0}{\partial \mathbf{x} \partial y}\right) + A\_{66} \left(\frac{\partial^2 u\_0}{\partial y^2}\right) + \left(1 - \mu \nabla^2\right) \left[ -I\_0 \frac{\partial^2 u\_0}{\partial t^2} - I\_1 \frac{\partial^2 \phi\_\mathbf{x}}{\partial t^2} \right] = 0 \,, \tag{30}$$

*Appl. Sci.* **2019**, *9*, 887

$$A\left(A\_{12} + A\_{66}\right)\left(\frac{\partial^2 u\_0}{\partial x \partial y}\right) + A\_{66}\left(\frac{\partial^2 v\_0}{\partial x^2}\right) + A\_{22}\left(\frac{\partial^2 v\_0}{\partial y^2}\right) + \left(1 - \mu \nabla^2\right) \left[-I\_0 \frac{\partial^2 v\_0}{\partial t^2} - I\_1 \frac{\partial^2 \phi\_Y}{\partial t^2 \partial y}\right] = 0$$

$$\begin{split} \mathcal{J}\_{55} \begin{pmatrix} \frac{\partial^2 w\_0}{\partial x^2} + \frac{\partial \phi\_x}{\partial x} \\\\ + (1 - \mu \nabla^2) \end{pmatrix} &+ \left(1 - \mu \nabla^2 \right) \begin{pmatrix} -k\_w w\_0 + k\_p \left(\frac{\partial^2 w\_0}{\partial x^2} + \frac{\partial^2 w\_0}{\partial y^2} \right) - \left(P\_{xx} + N\_{xx}^T \right) \frac{\partial^2 w\_0}{\partial x^2} \\\ - \left(P\_{yy} + N\_{yy}^T \right) \frac{\partial^2 w\_0}{\partial y^2} - I\_0 \frac{\partial^2 w\_0}{\partial t^2} \end{pmatrix} = 0 \,,\tag{32}$$

$$\begin{array}{l} D\_{11} \left( \frac{\partial^2 \phi\_x}{\partial x^2} \right) + D\_{12} \left( \frac{\partial^2 \phi\_y}{\partial x \partial y} \right) + D\_{66} \left( \frac{\partial^2 \phi\_x}{\partial y^2} + \frac{\partial^2 \phi\_y}{\partial x \partial y} \right) - f\_{\mathbb{55}} \left( \frac{\partial \mathbf{m}\_{\mathbb{C}}}{\partial x} + \phi\_x \right) \\ + \left( 1 - \mu \nabla^2 \right) \left( -I\_1 \frac{\partial^2 \boldsymbol{u}\_0}{\partial t^2} - I\_2 \frac{\partial^2 \phi\_x}{\partial t^2} \right) = \mathbf{0} \end{array} \tag{33}$$

$$\begin{split} D\_{12} \left( \frac{\partial^2 \phi\_x}{\partial x \partial y} \right) + D\_{22} \left( \frac{\partial^2 \phi\_y}{\partial y^2} \right) &+ D\_{66} \left( \frac{\partial^2 \phi\_x}{\partial x \partial y} + \frac{\partial^2 \phi\_y}{\partial x^2} \right) - f\_{44} \left( \frac{\partial \mathbf{u}\_0}{\partial y} + \phi\_y \right) \\ &+ \left( 1 - \mu \nabla^2 \right) \left( -I\_1 \frac{\partial^2 \mathbf{v}\_0}{\partial t^2} - I\_2 \frac{\partial^2 \phi\_y}{\partial t^2} \right) = 0 \end{split} \tag{34}$$

#### **5. Solution Procedure**

The governing differential equations are solved analytically via Navier's approach for simply-supported graphene sheet in this research. According to the selected approach, the mid-plane displacement and rotation field are solved by means of the Fourier series as

$$\mu\_0(x, y, t) = \sum\_{m=1}^{\infty} \sum\_{n=1}^{\infty} l L\_{mn}(t) \cos(a\omega t) \sin(\beta y) \,\,\,\,\tag{35}$$

$$w\_0(x, y, t) = \sum\_{m=1}^{\infty} \sum\_{n=1}^{\infty} V\_{mn}(t) \sin(a\omega) \cos(\beta y),\tag{36}$$

$$w\_0(x, y, t) = \sum\_{m=1}^{\infty} \sum\_{n=1}^{\infty} \mathcal{W}\_{mn}(t) \sin(ax) \sin(\beta y),\tag{37}$$

$$\Phi\_{\mathbf{x}}(\mathbf{x}, y, t) = \sum\_{m=1}^{\infty} \sum\_{n=1}^{\infty} X\_{mn}(t) \cos(a\mathbf{x}) \sin(\beta y),\tag{38}$$

$$\phi\_y(x, y, t) = \sum\_{m=1}^{\infty} \sum\_{n=1}^{\infty} Y\_{mn}(t) \sin(ax) \cos(\beta y),\tag{39}$$

where *α* = *<sup>m</sup> <sup>π</sup> <sup>a</sup>* , *<sup>β</sup>* <sup>=</sup> *<sup>n</sup> <sup>π</sup> <sup>b</sup>* , and *m* is the half wave number in the *x* direction and *n* is the half wave number in the *y* direction.

Finally, substituting Equations (35)–(39) into Equations (30)–(34), leads to the following matrix form, namely

$$\mathbf{\dot{M}d} + (\mathbf{K}\_L - P(t)\mathbf{K}\_G)\mathbf{d} = 0,\tag{40}$$

in which **<sup>d</sup>** <sup>=</sup> {*UVWXY*}*<sup>T</sup>* denotes the displacement vector. Additionally, **<sup>M</sup>** refers to the mass, **K***<sup>L</sup>* describes the stiffness, and **K***<sup>G</sup>* stands for the geometric stiffness matrices.

We consider a periodic compressive load which includes a static and dynamic part, as follows:

$$P(t) = P\_s + P\_d \cos(\omega t) = aP\_{cr} + \beta P\_{cr} \cos(\omega t). \tag{41}$$

Here, *Ps* refers to the static force component and *Pd* indicates the dynamic force component, *Pcr* shows the static buckling load and *ω* describes the pulsation frequency. Additionally, *α* is the static load factor and *β* is the dynamic load factor. The substitution of Equation (41) into Equation (40) read

$$\mathbf{M}\ddot{\mathbf{d}} + (\mathbf{K}\_L - (\alpha + \beta \cos(\omega t))P\_{cr}\mathbf{K}\_G)\mathbf{d} = 0 \,. \tag{42}$$

Equation (42) is a Mathieu–Hill type equation discussing the dynamic instability behavior of graphene sheet under a periodic compressive load. The Bolotin's method [38] is employed to determine the boundaries of the DIR. Based on this method, the displacement vector of *d* can be written in the following series with period 2*T* as

$$\mathbf{d} = \sum\_{k=1,3,\dots}^{\infty} \left[ a\_k \sin \frac{k \,\omega \,\, t}{2} + b\_k \cos \frac{k \,\, \omega \,\, t}{2} \right] \,' \,. \tag{43}$$

in which *ak* and *bk* denote arbitrary constant vectors. On the basis of Bolotin's approach [38], the first approximation with *k* = 1 could be used to compute the instability boundary. In this case, by setting Equation (43) into Equation (42) and considering a null value for the coefficients of the harmonic functions and the sum of the constant terms to zero, we get to the following relation:

$$\left| \left( \mathbf{K}\_L - P\_{cr} \left( \mathbf{a} \pm \frac{\beta}{2} \right) \mathbf{K}\_G - \mathbf{M} \frac{\omega^2}{4} \right) \right| = 0 \,\,\,\,\tag{44}$$

At the end, we can compute the critical static buckling load, assuming *β* = *ω* = 0. On the basis of the eigenvalue problem, the variation of *ω* versus *β* can be plotted as instability regions of the nanoplate to the periodic load including thermal effect.

#### **6. Numerical Examples**

This section is devoted to the numerical study on the dynamic instability of orthotropic temperature-dependent SLGS resting on temperature-dependent elastomeric foundation undergoing biaxial oscillating loading in a thermal environment. Here, two types of zigzag SLGSs are considered with two different values of aspect ratio. The material properties of the SLGS are highly dependent to the temperature. Shen et al. [39] performed a MD simulation to evaluate these properties at three different temperatures, as listed in Table 1. In addition, a polydimethylsiloxane (PDMS) is here selected for the elastomeric foundation, while assuming the following material properties: a Poisson's ratio *ν<sup>s</sup>* = 0.48 and a Young's modulus *Es* = (3.22 − 0.0034*T*) *GPa*, where *T* = *T*<sup>0</sup> + Δ*T* and *T*<sup>0</sup> = 300 *K* refers to the room temperature [37].


**Table 1.** Thermo-mechanical properties of the material [39].

#### *6.1. Validation of Results*

A preliminary validation of the proposed mathematical formulation is performed against some results available in literature, see Tables 2 and 3. More specifically, Table 2 summarizes the results in terms of dimensionless natural frequencies for a simply-supported square graphene sheet, as computed by the present method, against predictions by Sobhy [40] based on SSDT, and or those ones by Ebrahimi and Barati [41] according to the higher-order refined plate theory and the DQM approach. A parametric investigation analyzes the variation of the NPs and foundation constants. The outcomes in Table 2 clearly show that the natural frequencies gradually decrease for an increased value of *μ* and increased for increasing values of KW and KP. The very good agreement between our results and predictions from literature, clearly confirms the reliability of the present approach to capture the

response of the problem, where few negligible differences are only related to the different selected solution methodology and plate theory.


**Table 2.** Dimensionless natural frequencies of the square graphene sheet for different nonlocal and foundation parameters (*a*/*h* = 10).

**Table 3.** Comparative evlauation of the critical buckling load (nN/nm) for isotropic square graphene sheet under biaxial compression loads.


As another attempt for validation, we compute the critical buckling loads for isotropic square graphene sheet under biaxial load with different lengths, whose results are summarized and compared in Table 3 with the ones based on an MD approach [42], a DQM approach [43], or a refined strip method (RFSM) [19]. The gradual increase of the geometrical dimension *a* yields a meaningful reduction of the critical loading, whereas the very good agreement between our results and those reported in literature, confirms once again the accuracy of the proposed method.

#### *6.2. Parametric Studies and Discussions*

In the following examples, the DIR is traced as a dimensionless excitation frequency Ω = *ωb*<sup>2</sup> &*ρh*/*D*<sup>11</sup> in which *<sup>D</sup>*<sup>11</sup> = *<sup>E</sup>*1*h*3/12(<sup>1</sup> − *<sup>ν</sup>*12*ν*21) against the dynamic load factor *<sup>β</sup>*. It was supposed, unless otherwise explained, that *a*/*b* = 1 (zigzag sheet I) and *μ* = 0.5 nm2.

Figure 2 investigates the influence of the temperature change on the dynamic instability of orthotropic graphene sheet without elastomeric foundation. As clearly visible in Figure 2, an increase in temperature change reduces the excitation frequencies and shifts the origins of the DIR to the left side and also decreases gradually its width. The temperature rise leads to an increased compressive prestress in the graphene sheet, thus weakening the nanoplate stiffness.

**Figure 2.** Temperature variation effect on the dynamic instability region (DIR) of the graphene sheet without elastomeric foundation.

As mentioned before, in this study we analyze the effect of the temperature on the mechanical properties of both the nanostructure and elastic foundation. Thus, Figure 3a,b plot the DIR of the SLGS resting on the elastomeric foundation under the twofold assumption of temperature-independent and temperature-dependent material properties. A comparative evaluation of the results is illustrated in Figure 3a,b, for two different values of temperature, i.e., *T* = 500 K and *T* = 700 K, respectively. As visible in Figure 3a,b, when temperature-dependent properties are used, the origins of the instability regions occur at smaller excitation frequencies. The difference between the two cases becomes even more pronounced as the temperature is enhanced. This means that it is significant to account for the effect of temperature on the material properties in order to obtain more reliable results.

**Figure 3.** Comparison between the DIR of SLGS resting on the elastomeric foundation for constant and temperature-dependent properties.

The surrounding elastic foundation effect on the DIR of SLGS is displayed in Figure 4a,b, for *T* = 300 K and *T* = 700 K, respectively. As visible in Figure 4a,b, the presence of an elastic foundation increases the excitation frequency and hence, the origin of unstable region moves to the right side. The existence of an elastomeric foundation magnifies the nanoplate stiffness which yields an increased stability of the nanoplate. In addition, by considering an elastic foundation, the instability region of the graphene sheet tends to become wider, especially at increased levels of temperature. Moreover, it can be observed that the effect of the elastomeric Pasternak medium (KW = 0 and KP = 0) is higher

than the elastomeric Winkler (KW = 0 and KP = 0) one, on the DIR of the SLGS. The elastomeric Winkler foundation is able to describe only the normal loading condition of the elastomeric foundation, whereas the Pasternak model considers normal stresses as well as transverse shear loading conditions. From these figures, it is also worth observing that the foundation effects become even more significant for increased temperatures.

**Figure 4.** Effect of the elastomeric foundation on the DIR of the SLGS.

Figure 5a,b illustrates the aspect ratio *a*/*b* effect on the unstable region of graphene sheets embedded in the elastomeric foundation for *T* = 300 K and *T* = 700 K, respectively. In this example, two different values of aspect ratio (i.e., *a*/*b* 1 (zigzag sheet I) and *a*/*b* 2 (zigzag sheet II)) are considered [39]. This figure indicates that the excitation frequency reduces and the dynamic instability moves to the left side as *a*/*b* increases. This trend is almost predictable, since the nanoplate becomes softer for increasing aspect ratios, with a global reduction of the structural stiffness. Moreover, the width of the instability region increases with a decrease in the slenderness ratio.

**Figure 5.** Aspect ratio effect on the DIR of the SLGS embedded on an elastomeric foundation.

The NP effect on the DIR of embedded SLGS is displayed in Figure 6a,b, for *T* = 300 K and *T* = 700 K, respectively. The dimensionless pulsation frequency seems to decrease as the NP becomes higher, since the DIR of graphene sheet forms at lower excitation frequencies. Such phenomenon is owing to the fact that NP indicates the softening effect on stiffness which leads to lower pulsation frequency. Furthermore, comparing these figures exhibits that an increase in temperature change increases the effect of NP on the unstable region of the nanoplate.

**Figure 6.** Effect of the NP on the DIR of the SLGS resting on an elastomeric foundation.

Finally, Figure 7a,b demonstrates the static load factor effect on the unstable region of the SLGS resting on the elastomeric foundation for *T* = 300 K and *T* = 500 K, respectively. As can be concluded, with an increase in the static load factor the excitation frequency decreases and therefore, dynamic instability shifted to the left side due to the reduction in the stiffness of the nanoplate. Furthermore, the results indicate that the DIR becomes wider for a higher static load factor. Any static force component is present, in the particular case when *α* = 0.

**Figure 7.** Static load factor effect on the DIR of the SLGS embedded on an elastomeric foundation.

#### **7. Conclusions**

This study investigated the dynamic instability of temperature-dependent SLGS embedded in an elastic medium in thermal environment. This was modelled by means of a temperature-dependent elastomeric foundation. Based on FSDT within the framework of Eringen's differential constitutive model, the equations of motion were here derived from the Hamilton's variational principle as well as the energy method. The Navier's solution in conjunction with the Bolotin's approach were here applied to calculate the unstable region for the graphene sheet. A large numerical investigation examined the influence of different parameters such as the NP, temperature change, aspect ratio, static load factor, and foundation type on the dynamic stability behavior. Based on the numerical results, it was found that the elastic foundation increases the width of the instability region and the origin of the unstable region shift to the right side which leads to higher excitation frequencies. An increased static load factor and/or temperature weakens the nanoplate stiffness and leads to a lower pulsating frequency, while moving the DIR to the left side. It seemed also that the DIR occurs at smaller values of excitation frequency, when temperature-dependent mechanical properties were assumed for both the SLGS and the foundation, compared to the assumption of constant properties.

**Author Contributions:** Conceptualization, M.H.J., R.D. and F.T.; Formal analysis, M.H.J., R.D. and F.T.; Investigation, M.H.J. and F.T.; Validation, M.H.J., R.D. and F.T.; Writing—Original Draft, M.H.J., R.D. and F.T.; Writing—Review & Editing, R.D. and F.T.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Nonlocal and Size-Dependent Dielectric Function for Plasmonic Nanoparticles**

**Kai-Jian Huang 1,2, Shui-Jie Qin 1,2,\*, Zheng-Ping Zhang 1, Zhao Ding <sup>1</sup> and Zhong-Chen Bai <sup>2</sup>**


Received: 4 July 2019; Accepted: 26 July 2019; Published: 31 July 2019

**Abstract:** We develop a theoretical approach to investigate the impact that nonlocal and finite-size effects have on the dielectric response of plasmonic nanostructures. Through simulations, comprehensive comparisons of the electron energy loss spectroscopy (EELS) and the optical performance are discussed for a gold spherical dimer system in terms of different dielectric models. Our study offers a paradigm of high efficiency compatible dielectric theoretical framework for accounting the metallic nanoparticles behavior combining local, nonlocal and size-dependent effects in broader energy and size ranges. The results of accurate analysis and simulation for these effects unveil the weight and the evolution of both surface and bulk plasmons vibrational mechanisms, which are important for further understanding the electrodynamics properties of structures at the nanoscale. Particularly, our method can be extended to other plasmonic nanostructures where quantum-size or strongly interacting effects are likely to play an important role.

**Keywords:** EELS; plasmons vibrational modes; nanoparticles; nonlocal and size-dependent dielectric

#### **1. Introduction**

Classical local electrodynamics in nanostructures under the assumption that the polarization characteristics at a given point are locally related to the external electromagnetic excitation, which has been thoroughly explored from isolated nanoparticles and dimers to other complex structures [1–5]. On the other hand, the nonlocal response of a metal nanostructure to an external disturbance stimulus determines many important plasmonic properties and has attracted much theoretical interest [6–9]. In recent years, the Mortensen group [10] have studied plasmons in metal nanoparticles by using a popular hydrodynamic model which incorporates nonlocal effects by introducing a pressure term that inhibits squashing electrons into small plasmonic nanoparticles. Similarly, the work of Christensen et al. [11] included nonlocality in the optical response of metallic nanoparticles by introducing quantum corrections obtained from first principles calculations, also using the hydrodynamic model. Ming et al. [12,13] have demonstrated a Maxwell-hydrodynamic model for studying the charge, energy and angular momentum conservation laws of high-order harmonics due to nonlocal and nonlinear interactions between electromagnetic waves and metallic metamaterials in the terahertz regime.

Meanwhile, parallel to the above remarkable progress of nonlocal properties, there has been a growing interest in exploring the quantum properties of surface plasmon and the prospects of plasma devices based on reliable operation at the quantum level. However, theoretical investigation of plasmonic resonances for nanostructure sizes ranging between 2 and 20 nm remains a challenge. This length scale is too small for classical electromagnetic theories to be valid but too large for full-fledged quantum simulations to be feasible. Therefore, some of the studies are committed to provide a unified theoretical framework that breaks wide discrepancy between classical

electromagnetic theory and quantum mechanical theory for plasmonics and nanophotonics, and ultimately control plasmonic responses in nanostructures. Crucial to this endeavor is to understand finite-size effects in nanostructures since the control of nanoparticle size represents one of the most effective approaches to harness the tremendous potentials of plasmonics and nanophotonics.

Recent studies have shown a theoretical method with a suitable modification of the bulk complex dielectric function for free and bound electrons within metallic nanoparticles [14–19], for the determination of plasma frequency, damping constants, and size distributions of noble metal nanoparticles, such as gold, silver and copper; in addition, the corresponding optical characterization deduced by Mie's theory has been tested to fit full UV-Vis experimental extinction spectra. The above theoretical framework follows the idea of Kreibig and von Fragstein [20], who considered that the augmentation of the damping constant in the Drude model results from additional collisions of free electrons with the particle boundary, as well as the idea of Inouye et al. [21], who took into account the interband transitions of bound electrons for the complex dielectric function. Consequently, the metal dielectric function can be regarded as the sum of the free (intraband transitions) and bound (interband transitions) electron contributions; then, a size-dependent (SD) dielectric function (Equation (4)) based on a stepwise modification of the Drude model is put forward. Under this theoretical framework, the significant finite-size effects of optical properties for metallic nanoparticles can be expected to be observed below a certain limited size range (e.g., valid restriction of 20 nm drawn from Refs. [14,15]).

A developed microscopic approach based on accurate inclusion of the Lorentz friction for the plasmon damping within the random phase approximation (RPA) model [22–24] has also been performed in recent years. Unlike the conventional Mie's theory that only gives reasonably intrinsic size effects for nanospheres of radii lower than 5 nm in vacuum due to neglecting the Lorentz friction losses that dominate over all attenuation channels, the specific size dependence for plasmons due to the Lorentz friction induced corrections derived upon this theoretical approach can be applied to account for the observed irregular size-effect of plasmon resonance for relatively large metallic nanoparticles (e.g., for Au nanospheres of radii 20–60 nm in vacuum drawn from Ref. [24]). This theoretical approach has also been verified to improve experimental data fitting of metallic nanospheres measured in water colloidal solutions, especially pronounced for a large size limit (e.g., of radii 10–75 nm for Au nanospheres in water colloidal solutions drawn from Ref. [22]). One should note that, although accurate accounting for the Lorentz friction is important for nanospheres of radii above 20 nm (for Au in vacuum), the major drawback of our work is that adequate corrections of the Lorentz friction are not included, since integrating them with the nonlocal and size-dependent effects within a macroscopic material dielectric function is not straightforward.

Extending beyond the nonlocal and size-dependent effects we have focused on, in the following, there is a huge amount of research concerning about other "non-classical electromagnetic effects" by incorporating some "quantum effects" within classical electrodynamics. Quantum-size effects due to the confinement of electron motion in small dimensions can lead to a discrete electronic energy levels and a size-dependent conductivity [25]. Accessing plasmonic modes in the quantum tunneling regime also remains as a challenge, quantum tunneling effects become significant for subnanometer interparticle gaps and can be reflected through a quantum corrected model (QCM), which models the junction between adjacent nanoparticles by means of a local dielectric response that includes electron tunneling and tunneling resistivity at gap [26,27]. However, such "non-classical" effects are of great interest and significance in the field of nano-optics, and can be incorporated into classical electrodynamical descriptions for some simple systems; directly integrating them within a general framework is not straightforward. In particular, quantum tunneling effects between metal nanoparticles will not be considered in this work.

Nonlocal effects are a direct consequence of the inhomogeneity of the electron gas, and can hardly be ignored when the mean free path of an electron in the simple systems is larger than the wavelength of light [28,29]. To a certain extent, nonlocal electrodynamics partially account for finite-size and quantum mechanical effects inherent in a nanostructure in a more rigorous way. Such as in the case of

applying the hydrodynamic dielectric function to calculate the absorption coefficient from metallic nanoparticles, its size-dependent effects will show up by an effective dielectric function due to the additional boundary conditions for solving polarizabilities [30]. Although it has been argued that the hydrodynamic model yields only semiquantitative results in the linear regime, recent improvements regarding the incorporation of interband transitions and Landau damping have demonstrated that fully quantitative results may be obtained [31,32]. In principle, as most of the modified hydrodynamic models shown in the aforementioned applications, a spatially nonlocal relationship between the material polarization and the field of plasmonics can typically be mirrored through the use of a spatially nonlocal dielectric function. Some of them take finite-size effects of the nanostructure into consideration by using a modified collision frequency. However, the corresponding changes in optical, plasmonic and dielectric responses are almost negligible, and the effective dimension of nanostructure for nonlocal electrodynamics is also being limited or obscured consequently.

Inspired by the studies of the size-dependent dielectric function as mentioned above, in this paper, we move a step forward in this direction by developing nonlocal and size-dependent (NLSD) dielectric approach that incorporates the spatially local, nonlocal, size and even analogous quantum-size responses of the material. As an example, we investigate how the electron energy loss spectroscopy (EELS) and optical responses of gold spherical dimers of various radii and inter-particle gaps are modified once the nonlocality and finite-size effects are considered. These examples allow for a straightforward interpretation of our results. Unlike those complex configurational and resource consuming methods, such as quantum mechanical simulations based on the time-dependent Kohn–Sham density functional theory (TD-KSDFT), our approach provides a powerful and simple tool to calculate and analyze the EELS and optical response ranging from classical to quantum regimes.

#### **2. Methods and Theoretical Framework**

Figure 1 displays the fundamental system under study; a dimer consists of two closely spaced freestanding metallic nanospheres of radius *R*, and a certain distance *d* apart. The system is excited by an electron beam with high kinetic energy up to 200 keV passing by or penetrating through at the positions as indicated in the figure, and illuminated from the top by a plane wave polarized along the dimer axis, respectively. For this study, the metal is assumed to be gold, and the experimental and theoretical values of permittivity for contrast taken from Johnson and Christy [33] and Werner et al. [34], respectively. We chose the data models above for comparison for the following reasons. First, of all, in general, ab initio density functional theory (DFT) approaches need to deal with pseudopotentials, and include quantum effects that completely change the spectral distribution and the field enhancements of the plasmon resonant modes supported by a given nanostructure; they are computationally demanding and are currently limited to small-sized nanoparticles with a few thousand electrons (e.g., valid restriction of 12.3 nm drawn from Ref. [35]). Second, the application of classic local experimental dielectric data for small-sized nanoparticles ranging between 2 and 20 nm becomes questionable as described in the Introduction. These factors will be discussed in detail below combined with calculation results. All of the numerical simulations are implemented via the general public license toolbox MNPBEM [36,37], which is based on a boundary element method (BEM) approach.

The dielectric function of noble metal nanoparticles (e.g., Cu, Au, and Ag) is well described by using the hydrodynamical Drude model limit by three components [38–40]:

$$
\varepsilon\_{\text{Normal}}(\mathbf{k}, \omega) = \varepsilon\_{\text{\textquotedblleft}} + \varepsilon\_{\text{inter}}(\omega) + \varepsilon\_{\text{intra}}(\mathbf{k}, \omega), \tag{1}
$$

where the **k** dependence in *ε*Nonlocal(**k**,*ω*) is necessary to describe the spatially nonlocal response of structures with features less than 10 nm [41,42]. *ε*<sup>∞</sup> is the value as *ω* → ∞. *ε*inter(*ω*) represents the contribution from interband electron transitions and can be physically described in the form of a Lorentz oscillator model,

$$
\varepsilon\_{\rm inter}(\omega) = \sum\_{j} \frac{\Delta \varepsilon\_{\rm Lj} \omega\_{\rm Lj}^2}{\omega\_{\rm Lj}^2 - \omega(\omega + i2\delta\_{\rm Lj})}. \tag{2}
$$

In this work, we take index *j* = 2 labeling the individual d-band to sp-band electron transitions since there are two interband transitions in Au at optical frequencies [33]. *ε*intra(**k**,*ω*) is responsible for the intraband effects which are usually referred to as free-electron effects, where both the plasmonic optical response of metals and nonlocal effects are contained and can be described by

$$\varepsilon\_{\rm intra}(\mathbf{k},\omega) = -\frac{\omega\_{\rm D}^2}{\omega\left(\omega + i\gamma\right) - \beta^2 \mathbf{k}^2},\tag{3}$$

where *ω*<sup>D</sup> is the plasma frequency, *γ* is the collision frequency, and *β* is the hydrodynamical parameter measuring the degree of nonlocality. For the three-dimensional situation *β*<sup>2</sup> = 3*υ*F/5, where *<sup>υ</sup>*<sup>F</sup> = 1.39 × <sup>10</sup><sup>6</sup> m/s is the Fermi velocity for the Au nanoparticles. Then, all parameters in Equations (2) and (3) can be estimated by fitting Equation (1) to the experimentally determined dielectric data of bulk Au in Ref. [33]. Here, we use all the parameters taken from Ref. [39], *ω*<sup>D</sup> = 8.812 eV, *γ* = 0.0752 eV, *ε*<sup>∞</sup> = 3.559, Δ*ε*L1 = 2.912, *ω*L*<sup>j</sup>* = 4.693 eV, *δ*L1 = 1.541 eV, Δ*ε*L2 = 1.272, *ω*L*<sup>j</sup>* = 3.112 eV, *δ*L1 = 0.525 eV.

**Figure 1.** Schematic illustration of the system under study: the interaction of the electron beam passes by or penetrates through a dimer consisted with two Au nanospheres with radius R and gap d, as well as the plane wave polarized along the dimer axis. The inset shows the exact located positions (or impact parameters) of electron beams: the red ring represents the middle of the gap and the blue cross represents the centre in one of the nanospheres.

Note that the above dielectric function calculation model neglects quantum-mechanical exchange and correlation effects, but the significant nonlocal effects from the reduced mean free path of the conduction electrons due to electron–interface scattering in small metal structures can be taken into account by using a modified collision frequency *γ*- = *γ* + 4*Aυ*F*V*/*S* in Equation (3) for the three-dimensional problems [43]. In principle, conclusions from the previews studies that adopt such modified collision frequency have shown that the Nonlocal model of Equation (1) can model the spatially nonlocal dielectric response of arbitrarily shaped and sized structures [39,40]. However, although the volume of the structure *V* and surface area *S* in the modified collision frequency implies the influence of size factor, only the impact of free electrons was considered to be modified by size. In addition, correctly selecting the value of dimensionless parameter *A* remains a challenging since its complex details are of discrepancy in different theoretical dielectric models (consult Section 1 in Ref. [44] and its references for a detailed discussion). Therefore, in order to overcome this difficulty in a classical mean free path approach, a more reasonable alternative way for describing size-dependent dielectric functions is needed by considering the dominated conditions and a scope of the contributions of intraband transitions (free electrons) and interband transitions (bound electrons). Herein, according to

previous studies about suitable modification of the bulk complex dielectric function for free and bound electrons within metallic nanoparticles [14–16], the full size-corrected dielectric function can be written as the sum of three terms:

$$
\varepsilon\_{\rm SD} \left( \omega, \boldsymbol{\mathcal{R}} \right) = \varepsilon\_{\rm bulk} \left( \omega \right) + \Delta \varepsilon\_{\rm free} \left( \omega, \boldsymbol{\mathcal{R}} \right) + \Delta \varepsilon\_{\rm bound} \left( \omega, \boldsymbol{\mathcal{R}} \right), \tag{4}
$$

where *ε*bulk (*ω*) corresponding to the experimentally measured bulk values from UV-visible to NIR-FIR, and the last two terms are size-corrective contributions for free and bound electrons, respectively:

$$\Delta \varepsilon\_{\text{free}}\left(\omega, \mathbb{R}\right) = -\omega\_{\text{D}}^{2} \sum\_{n=1}^{\infty} \left(-1\right)^{n} \frac{\left(i\omega \text{C} \upsilon\_{\text{F}}/\text{R}\right)^{n}}{\left(\omega^{2} + i\omega \gamma\_{\text{free}}\right)^{n+1}}\tag{5}$$

$$\begin{split} \Delta \varepsilon\_{\text{bound}} \left( \omega, R \right) &= -K\_b e^{R/R\_0} \int\_{\omega\_\S}^{\infty} \frac{\sqrt{\mathbf{x} - \omega\_\S}}{\mathbf{x}} \left[ 1 - F(\mathbf{x}, T) \right] \\ &\times \frac{d\mathbf{x}}{\mathbf{x}^2 - \omega^2 + \gamma\_{\text{bound}} - 2i\omega\gamma\_{\text{bound}}}, \end{split} \tag{6}$$

where *h*¯ *ω<sup>g</sup>* is the gap energy; *F*(*x*, *T*) is the Fermi distribution function of conduction electron of energy *hx*¯ and temperature *T* with Fermi energy *E*F, *γ*free = *γ*, *γ*bound stands for the damping constant in the band to band transition; *Kb* is a proportionality factor with units of *s*<sup>−</sup>3/2; *C* is a constant that depends on the material. For noble metals, it has a value of about 0.8 that has been theoretically justified from first principles calculations [45]; *R*<sup>0</sup> is a reference radius value that represents the range for which the density of states can be considered to reach the value of the bulk. For Au, the following values determined in Refs. [14–16] were used: *<sup>ω</sup><sup>g</sup>* = 2.1062 eV, *<sup>γ</sup>*bound = 0.158 eV, *Kb* = 2.3 × <sup>10</sup><sup>24</sup> <sup>s</sup><sup>−</sup>3/2, *E*<sup>F</sup> = 2.5 eV.

Note that the interband transitions from the d-band to the conduction sp-band near the L point in the Brillouin zone can be taken into account by the size-corrective contribution for bound electrons in the form of Equation (6) [46]. Because of the nonlocal dielectric function, *ε*Nonlocal(**k**,*ω*) is estimated by fitting Equation (1) to the complex bulk dielectric function, which is the local, frequency-dependent part of the response, taken from optical measurements given in Ref. [33] and exactly the same as the first term *ε*bulk (*ω*) using the size-dependent dielectric function *ε*SD (*ω*, *R*) of Equation (4); we adopt here the former as a substitute for the latter. The resulting nonlocal and size-dependent (NLSD) dielectric function of the metal is thus approximated by

$$
\varepsilon\_{\text{NLSD}}(\mathbf{k},\omega,\mathbf{R}) = \varepsilon\_{\text{Noxlocal}}(\mathbf{k},\omega) + \Delta\varepsilon\_{\text{free}}\left(\omega,\mathbf{R}\right) + \Delta\varepsilon\_{\text{bound}}\left(\omega,\mathbf{R}\right). \tag{7}
$$

#### **3. Results and Discussion**

The first example of the compatible effectiveness produced by our NLSD approach is shown in Figure 2 for comparisons of EELS calculated with different dielectric models and gold spherical dimers of different sizes, where electron beam penetrates the centre in one of the nanospheres (Figure 2a) and passes by the middle of the gap (Figure 2b) as shown in Figure 1, and the gap *d* fixed on 1 nm. Figure 2a shows that the agreement of the total energy loss probabilities between the NLSD approach(black solid line) and the theoretical dielectric data set (green dashed line), which is obtained from density functional theory (DFT) calculations by Werner et al. [34], is remarkable for radii up to 20 nm. Moreover, for the NLSD approach, the evolves of the size-dependent response as the radius increases are observably more outstanding than the full size-dependent dielectric model (blue dotted line, Equation (4)), and the total energy loss probabilities become almost consistent with those calculated by using the local, frequency-dependent and experimental complex bulk dielectric values (red dash-dotted line) as radii up to 100 nm.

In addition, blueshifts of the main energy loss peaks become more apparent as the radius increases and it is accompanied by plasmon broadening. These nonlocal effects are known as the availability of the additional loss channels [7], and can have an effect in the higher energy region when compared to classical local electrodynamics because of the interplay between *ω* and **k** in Equation (3). All of the above agreements between the model predictions support that our approach provides a straightforward theoretical framework to describe metal behavior combining local, nonlocal and size-dependent effects in broader energy and size ranges simultaneously.

**Figure 2.** EELS (electron energy loss spectroscopy) spectrums for Au spherical dimers of different radii *R*, and for different dielectric models (black solid curves, NLSD (nonlocal and size-dependent) model; blue dotted curves, SD model; red dash-dotted curves, experimental dielectric data from Johnson et al. [33]; green dashed curves, DFT (density functional theory) calculated data from Werner et al. [34]; thinner gray solid curves, Nonlocal model), where the gap *d* fixed on 1 nm. (**a**) electron beam penetrating the centre in one of the nanospheres, and (**b**) passing by the middle of the gap, respectively, as indicated in Figure 1.

Moreover, according to the various developed theories of the SD model as mentioned in Section 1, the significant finite-size effects on the optical properties of nanoparticles can only be expected to be observed below a certain limited size range of about 20 nm. Our model captures the main aspects of the finite-size effects for small-sized nanoparticles. For radii below 20 nm, Figure 2b shows a significant difference of the total (and surface) energy loss probabilities between the NLSD approach and nonlocal dielectric model (thinner gray solid curves) of Equation (1). *A* = 0.1 was used in the modified collision frequency *γ*- = *γ* + 4*Aυ*F*V*/*S* when using the Nonlocal model [39,40]. In addition, unlike the outcomes between the SD model and the experimental dielectric data, which show remarkable differences in all size range, the total (and surface) energy loss probabilities calculated by using the NLSD model are almost identical with the Nonlocal model for radii larger than 20 nm this is why the forecasted results of Nonlocal model are not shown in Figure 2b where the radii larger than 20 nm. Similarly, for clarity, Figure 2a does not show the outcomes from the Nonlocal model because they are identical with those from the NLSD model in all size ranges.

In consideration of the size-corrective terms from SD and NLSD model being the same, these results might appear unexpected at the first view. The nonlocal effects themselves are also size-dependent. This can be seen from Figure 2b where radii are below 20 nm; besides capturing the blueshifts, the main energy loss peaks of the Nonlocal model are obviously less than those by using experimental dielectric data, but in a low-energy region (e.g., 0.65–2.4 eV), the outcomes of the Nonlocal model are higher than those by using experimental dielectric data, which have the same changing tendency as the results of SD and NLSD models. However, the NLSD model is a more suitable implementation for revealing the contribution of the surface vibrational mode with nanoparticles of the small radius. When compared with the Nonlocal model, the outcomes of the NLSD show that the relative reductions in the main energy loss peaks, and the increases of the outcomes in the low-energy region are more remarkable, as shown in Figure 2b. All of these are because the size-dependent effects due to the size-corrective terms are contained in the NLSD model.

For further confirmation of the present findings, EELS calculated with different dielectric models for isolated Au spherical nanoparticles of different sizes are presented in Figure 3, where electron beam passing near the nanoparticles and the contribution of the surface vibrational mode is dominant. The difference of the total (and surface) energy loss probabilities between the NLSD model and Nonlocal model becomes more pronounced with the reduction of particle size, as shown in Figure 3 where the radius *R* = 1 nm, and the prediction of our NLSD model are much more close to the outcomes from the SD model and the DFT model. The discrepancy in predicted energy loss probabilities between the NLSD model and the Nonlocal model would gradually narrow and ultimately disappear when the radius is greater than 26 nm, as shown in Figure 3b, in addition to the discrepancy in predicted energy loss probabilities between the SD model and experimental dielectric data. These results for isolated nanoparticles are more compatible with the conclusions that have been verified both experimentally and theoretically in the previous studies of the SD model in which the significant size-dependent effects can only be expected to be observed below a certain limited size range of approximately 20 nm.

**Figure 3.** EELS spectrums for isolated Au spherical nanoparticles of different radii *R*, and for different dielectric models as defined in the text, when electron beam passing near the nanoparticles as depicted in the inset. Distances *d*between electron beam and nanoparticles were fixed at 2 nm.

Although, under what specific physical mechanism will the differences due to the size-dependent effects between the SD model and experimental dielectric data disappear in the outcomes between the NLSD and Nonlocal model for dimers of a large radius merit further investigation. However, according to all of the above results of EELS, one can see that the NLSD model presented in this paper has the advantage of clearly separating the size regions where local, nonlocal or size-dependent effects are dominant, for both simulation cases (isolated particle or dimer).

First, the size-dependent effects somehow can be revealed in small isolated nanoparticles by using the SD model, whereas the outcomes of this model do not contain any appreciable impact of the nonlocal effects. In addition, at least, for dimers of a large radius, the outcomes of the SD model should be much more consistent with those by using the classic local experimental dielectric data—particularly for an electron beam penetrating the centre in one of the nanospheres and the contribution of the bulk vibrational mode being dominant, as shown in Figure 2a. By contrast, the size-dependent effects due to the same size-corrective terms in the NLSD model have been

demonstrated that, while they can be completely neglected for Au dimers of all size when the contribution of the bulk vibrational mode is dominant (Figure 2a), and, for Au dimers of a large radius, when the contribution of the surface vibrational mode is dominant (Figure 2b), they still have a great influence on dimers of a small radius, just as the results have been demonstrated for isolated particles (Figure 3).

Second, we should notice that the classic local experimental dielectric data used in this paper was obtained from reflection and transmission measurements on vacuum-evaporated thin films at room temperature [33]; the application of such local continuum electrodynamics for small-sized nanoparticles becomes questionable. This can be attributed to two effects. First of all, their imaginary part of the interband contribution due to size-dependent effects is not negligible in the spectral region of surface plasmon polariton resonance [14,15]. Second, nonlocal effects due to the scattering of the conduction electrons arise as the nanoparticle size is decreased, relative to metal films [39,40]. By contrast, either the EELS results or the optical properties that are going to be demonstrated in the following have shown that outcomes from our NLSD model can fit the DFT model better than other models for small-sized nanoparticles. These agreements benefit from the first term in the NLSD model, namely, the Nonlocal model of Equation (1), which have the advantage of describing the dynamical optical response of structures that are too large to treat using quantum mechanics but too small for classical electromagnetic theories to be valid [39,40].

It should be noted that a consistent energy loss with peak near 1.25 eV when the radius *R* = 100 nm (as indicated by the inset gray arrow in Figures 2a and 4(upper)), which is similar to the main plasmon resonances, will be presented below in Figure 6a for the extinction cross sections excited by optical plane wave, from approximately 955 nm. In order to gain physical insight into the origin of these common resonances in both situations of excitations, the surface and bulk configurations that form the total loss probabilities where *R* = 100 nm are plotted in Figure 4, respectively. The results shown in Figure 4 indicate that the total loss probabilities near 1.25 eV for any dielectric model are all dominated by surface loss probabilities. In addition, it is interesting that we found that the total energy loss probabilities shown in Figure 2b are also dominated by the surface loss probabilities, since all the bulk loss probabilities are almost zero when the electron beam located in the middle of the gap (not shown). We attribute this surface energy loss dominating effect to the contribution of the surface vibrational mode.

**Figure 4.** (Upper panel) Surface and (lower panel) bulk loss probabilities for gold spherical dimers of radius *R* = 100 nm, where the other simulation configurations remain the same as defined in Figure 2a.

One should note that, at least for now, the main energy loss peaks as depicted in Figure 2a are the coupling contribution of the surface and the bulk vibrational modes. However, it can be seen that these surface energy loss peaks near 1.25 eV as shown in Figure 4(upper) do not show up in Figure 2b where the contribution of the surface vibrational mode is dominant. In addition, according to our surface loss probability calculation results (not shown, similar to the form of Figure 4(upper)), they also do not show up when *R* = 50 nm where the other simulation configurations remain the same as defined in Figure 2a. However, the reducing contribution of the surface vibrational mode with nanoparticles of the small radius can be distinctly revealed by the NLSD model. For example, Figure 5 plots the simulated EELS probability maps of our dimer system with four different dielectric models for the 1.25 eV resonances. Our model captures the main aspects of the nonlocal electrodynamics for small-sized nanoparticles.

**Figure 5.** Total (top panels), bulk (middle panels) and surface (bottom panels) EELS probability maps for the gold dimers with particles radii *R* = 20 nm (left column) and *R* = 50 nm (right column), and for different dielectric models as defined in the text. The inter-particle gap *d* = 1 nm.

For radius *R* = 20 nm, the surface and bulk vibrational modes by using the NLSD model are higher than the others as shown in the middle and bottom panel of Figure 5 (left column). Thus, the energy loss extended throughout the whole surface of the gold spherical dimer, whereas those by using experimental dielectric data are confined in narrow rings formed by the surface and the horizontal cross-section of the dimer, as the total EELS probabilities shown in the left-top panel in Figure 5. We found similar results for dimers with radii less than 20 nm, except that the intensities of the energy loss calculated by using the NLSD approach the DFT model as the radius decreases.

For radius *R* = 50 nm, the total energy loss confined more strongly in narrow rings formed by the surface and the horizontal cross-section of the dimer for all dielectric models, and they are dominated by surface vibrational mode, as indicated by the EELS probability maps shown in Figure 5 (right column). We also found similar results for dimers with radii greater than 50 nm, except that the patterns of the total (and surface) energy loss changing from the narrow rings to crescents as the radius increases. These results mean that the nonlocal impact on the mode evolution of EELS and responses from the quantum mechanical mechanism can also be depicted by using our NLSD model. Although there are similar outcomes of mode evolution by using the SD and DFT models, they are obviously incompatible with the classic local experimental data where the size of the dimer becomes larger, such as shown in Figure 2 (e.g., *R* = 100 nm).

Finally, we use our NLSD approach to explore another interesting compatible effectiveness for the optical properties of the dimer system illuminated by a plane wave as defined in Figure 1. Although the effects of quantum natural and nonlocality have different origins, their impact on the optical response of the nanostructure can be measured by its physical parameters (in our case, which means the radius and gap). Figure 6a renders the maximum plasmonic resonances of extinction cross-sections that correspond to different radius for gold spherical dimers with the constant gap *d* = 1 nm, and for different dielectric models. The comparison between predictions from our NLSD model and the

other dielectric models clarifies how nonlocality and quantum mechanical effects modify the optical properties of the system.

**Figure 6.** (**a**) maximum plasmonic resonances of extinction cross-sections that correspond to different radius for gold spherical dimers with constant inter-particle gap *d* = 1 nm, and for different dielectric models as defined in text; (**b**) maximum field enhancement at the middle point of the inter-particle gap (the red circle depicted in Figure 1) as a function of *R* for gold spherical dimers and their corresponding plasmonic resonances as shown in (**a**).

For small radii (approximately less than 35 nm), there are significant blueshifts for the prediction of our NLSD model when compared to the predictions of SD model and experimental dielectric data, which are similar to the comparison results by using another entirely different nonlocal dielectric model of the metal from previous studies (see Figure 2 or Figure 4a in Ref. [7]), and can be accounted for the inclusion of nonlocal effects in the NLSD model, which also means that the size-dependent effects due to the size-corrective terms in the NLSD model can not affect the maximum plasmonic resonances of extinction cross-sections, since the predictions of our NLSD model and the SD model are approximately consistent with those of the Nonlocal model and experimental dielectric data, respectively, as shown in Figure 6a. One should note that these maximum plasmonic resonances of extinction cross-sections that correspond to different radius for gold spherical dimers of radii ranging between 20 and 60 nm may be more accurate if our model can be integrated with the adequate corrections of the Lorentz friction. However, nevertheless, either our NLSD model or the RPA approach included with the Lorentz friction induced corrections as discussed in Section 1 has observed anomalous shift of plasmonic resonance for nanoparticles with radii significantly exceeding the previously suggested limiting 5 nm by Mie's theory.

Figure 6b renders the field enhancement <sup>|</sup>E<sup>|</sup> <sup>|</sup>E0<sup>|</sup> at the middle point of the gap as a function of *R* for gold spherical dimers and their corresponding plasmonic resonances as shown in Figure 6a. Quite different from the extinction spectral profiles predicted, Figure 6b shows approximate tendency between the NLSD model and the DFT model. This means that the impacts of nonlocality and quantum confinement can be compatible, and both reduce the enhancement of fields. The quantum confinement effects are important when the mean free path of valence electrons is comparable to the particle size, and the mean free path of electrons in gold at room temperature is below 50 nm. These play a large role in interpreting the discrepancy in predicted field enhancements between the NLSD model and

experimental dielectric data would gradually narrow when the radius is greater than 50 nm, as shown in Figure 6b. Similar behavior occurs for the prediction of the DFT model, except that the radius which stands out for such change is approximately 65 nm. We attribute this difference to the inclusion of nonlocal and finite-size effects in the NLSD model, which lead to an increase in the contribution of the bulk vibrational mode and to reduce the contribution of the surface vibrational mode when the radius becomes less than 50 nm and gradually decreases, as the EELS results (Figure 5) discussed above. Both effects are not exclusive of small particles as the inset shown in Figure 6b.

In addition, they are also important for small inter-particle gaps since the strongly interacting between metallic geometries encountered in plasmonics. For example, as one can see from the maximum field enhancement as a function of separation for gold dimer with *R* = 50 nm shown in Figure 7, the predictions of the NLSD model and the DFT model are closely similar, and the gaps in predicted field enhancements between them and the SD model and the experimental dielectric data would gradually narrow or even vanish as the inter-particle gap *d* increases. As we can see both in Figures 6b and 7, between NLSD and nonlocal models, there is also a slight difference in their outcomes with a radius of less than approximately 65 nm (Figure 6b), or with a gap of less than approximately 0.5 nm (Figure 7). These are entirely caused by the size-dependent effects due to the size-corrective terms of the NLSD model, which also lead to a decreasing trend in the outcomes of the NLSD model toward the DFT model, particularly for dimers with small radii as the inset shown in Figure 6b. Furthermore, we also should note that the maximum field enhancement shown in Figure 7 may be not very accurate when the inter-particle gap *d* decreases below 0.5 nm, since the NLSD model may not be capable of calculating the quantum tunneling effects under such subnanometer gaps, as mentioned in Section 1.

**Figure 7.** Maximum field enhancement as a function of separation for gold dimer with *R* = 50 nm, and for different dielectric models as defined in text (and depicted in Figure 2).

#### **4. Conclusions**

In conclusion, we presented a simplified theoretical model to describe the nonlocal and size effects of nanoparticles. Our implementation is based on the hydrodynamic Drude model and two size-corrective terms derived from the size-dependent dielectric function. We have demonstrated the compatible effectiveness and the accuracy of our theoretical model through comparisons to analytical results of the EELS and the optical responses for gold spherical dimers. These analytical results demonstrated a number of effects that result from nonlocality and finite-size effects in the dielectric response, including analogous quantum confinement effects that decrease the electromagnetic field enhancement and their corresponding weight relations of the contributions of surface and bulk vibrational modes reflected by EELS, blueshifts of surface plasmon resonances and main energy loss peaks. From these results, we can state that the nonlocal effects are of great significance as they obviously determine the blueshifts of surface plasmon resonances for the gold dimers with radii smaller than approximately 35 nm and the electromagnetic field enhancement for the gold dimers with radii ranging from 20 to 50 nm. In addition, blending the adequate size-dependent effects and nonlocal effects can more accurately reflect reductions in the main energy loss peaks and the electromagnetic

field enhancement due to quantum confinement effects for the gold dimers with radii smaller than approximately 20 nm, as well as the weight relations of the contributions of surface and bulk vibrational modes for the gold dimers with radii ranging from 1 to 50 nm; for the gold dimers with radii larger than 50 nm, the nonlocal effects in the EELS and optical responses decrease as the radii increase further, and the local effects almost dominate the EELS and optical responses as radii up to 100 nm.

Our NLSD model used here accounts for the plasmonic nanoparticles behavior combining local, nonlocal and size-dependent effects, and reconciles the discrepancy between classical electromagnetic theory and quantum mechanical theory for plasmonics and nanophotonics in broader energy and size ranges. Importantly, although we have focused on gold spherical dimers, our approach can be extended to motivate new, and more precise experimental studies on other metallic nanostructures with complex configurations, particularly those where quantum-size effects and strongly interacting effects are likely to play a large role. In the future, we expect that our approach will make possible a deep exploration of the dissipative quantum electromagnetics, where couple loss or dissipation in a quantum system to a bath of oscillators maintained a strong connection with the classical electromagnetic system [47].

**Author Contributions:** Conceptualization, K.-J.H.; Formal analysis, K.-J.H.; Funding acquisition, Z.-C.B.; Methodology, K.-J.H.; Supervision, S.-J.Q.; Writing—original draft, K.-J.H.; Writing—review and editing, S.-J.Q., Z.-P.Z., Z.D. and Z.-C.B.

**Funding:** This research was funded by the National Natural Science Foundation of China (Grant No. 61741505) and the Guizhou Province Science and Technology Projects (No. QKHZ [2017]2887).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Evidence-Theory-Based Robust Optimization and Its Application in Micro-Electromechanical Systems**

#### **Zhiliang Huang 1,2, Jiaqi Xu 1, Tongguang Yang 1, Fangyi Li 3,\* and Shuguang Deng <sup>1</sup>**


Received: 21 February 2019; Accepted: 1 April 2019; Published: 7 April 2019

**Featured Application: This paper develops an evidence-theory-based robustness optimization (EBRO) method, which aims to provide a potential computational tool for engineering problems with epistemic uncertainty. This method is especially suitable for robust designing of micro-electromechanical systems (MEMS). On one hand, unlike traditional engineering structural problems, the design of MEMS usually involves micro structure, novel materials, and extreme operating conditions, where multi-source uncertainties inevitably exist. Evidence theory is well suited to deal with such uncertainties. On the other hand, high performance and insensitivity to uncertainties are the fundamental requirements for MEMS design. The robust optimization can improve performance by minimizing the effects of uncertainties without eliminating these causes.**

**Abstract:** The conventional engineering robustness optimization approach considering uncertainties is generally based on a probabilistic model. However, a probabilistic model faces obstacles when handling problems with epistemic uncertainty. This paper presents an evidence-theory-based robustness optimization (EBRO) model and a corresponding algorithm, which provide a potential computational tool for engineering problems with multi-source uncertainty. An EBRO model with the twin objectives of performance and robustness is formulated by introducing the performance threshold. After providing multiple target belief measures (Bel), the original model is transformed into a series of sub-problems, which are solved by the proposed iterative strategy driving the robustness analysis and the deterministic optimization alternately. The proposed method is applied to three problems of micro-electromechanical systems (MEMS), including a micro-force sensor, an image sensor, and a capacitive accelerometer. In the applications, finite element simulation models and surrogate models are both given. Numerical results show that the proposed method has good engineering practicality due to comprehensive performance in terms of efficiency, accuracy, and convergence.

**Keywords:** epistemic uncertainty; evidence theory; robust optimization; sensor design

#### **1. Introduction**

In practical engineering problems, various uncertainties exist in terms of the operating environment, manufacturing process, material properties, etc. Under the combined action of these uncertainties, the performance of engineering structures or products may fluctuate greatly. Robust optimization [1,2] is a methodology, and its fundamental principle is to improve the performance of a product by minimizing the effects of uncertainties without eliminating these causes. The concept of robustness optimization has long been embedded in engineering design. In recent years, thanks to the rapid development of computer technology, it has been widely applied to many engineering fields, such as electronic [3], vehicle [4], aerospace [5], and civil engineering [6]. The core of robustness optimization lies in understanding, measuring, and controlling the uncertainty in the product design process. In mechanical engineering disciplines, uncertainty is usually differentiated into objective and subjective from an epistemological perspective [7]. The former, also called aleatory uncertainty, comes from an inherently irreducible physical nature, e.g., material properties (elasticity modulus, thermal conductivity, expansion coefficient) and operating conditions (temperature, humidity, wind load). A probabilistic model [8–10] is an appropriate way to describe such uncertain parameters, provided that sufficient samples are obtained for the construction of accurate random distribution. Conventional robustness optimization methods [11–13] are based on probabilistic models, in which the statistical moments (e.g., mean, variance) are employed to formulate the robustness function for the performance assessment under uncertainties. On the other hand, designers may lack knowledge about the issues of concern in practice, which leads to subjective uncertainty, also known as epistemic uncertainty. The uncertainty is caused by cognitive limitation or a lack of information, which could be reduced theoretically as effort is increased. At present, the methods of dealing with epistemic uncertainty mainly include possibility theory [14,15], the fuzzy set [16,17], convex model [18,19], and evidence theory [20,21]. Among them, evidence theory is an extension of probability theory, which can properly model the information of incompleteness, uncertainty, unreliability and even conflict [22]. When evidence theory treats a general structural problem, all possible values of an uncertain variable are assigned to several sub-intervals, and the corresponding probability is assigned to each sub-interval according to existing statistics and expert experience. After synthesizing the probability of all the sub-intervals, the belief measure and plausibility measure are obtained, which constitute the confidence interval of the proposition and show that the structural performance satisfies a given requirement. Compared with other uncertainty analysis theories, evidence theory may be more general. For example, when the sub-interval of each uncertain variable is infinitely small, evidence theory is equivalent to probability theory; when the sub-interval is unique, it is equivalent to convex model theory; when no conflict occurs to the information from different sources, it is equivalent to possibility theory [23].

In the past decade, some progress has been made in evidence-theory-based robust optimization (EBRO). For instance, Vasile [24] employed evidence theory to model the uncertainties of spacecraft subsystems and trajectory parameters in the robust design of space trajectory and presented a hybrid co-evolutionary algorithm to obtain the optimal results. For the preliminary design of a space mission, Croisard et al. [25] formulated the robust optimization model using evidence theory and proposed three practical solving technologies. Their features of efficiency and accuracy were discussed through the application of a space mission. Zuiani et al. [26] presented a multi-objective robust optimization approach for the deflection action design of near-Earth objects, and the uncertainties involved in the orbital and system were qualified by evidence theory. A deflection design application of a spacecraft swarm with Apophis verified the effectiveness of this approach. Hou et al. [27] introduced evidence-theory-based robust optimization (EBRO) into multidisciplinary aerospace design, and the strategy of an artificial neural network was used to establish surrogate models for the balance of efficiency and accuracy during the optimization. This method was applied to two preliminary designs of the micro entry probe and orbital debris removal system. The above studies employed evidence theory to measure the epistemic uncertainties involved in engineering design, and expanded robustness optimization into the design of complex systems. However, the studies of EBRO are still in a preliminary stage. The existing research has mainly aimed at the preliminary design of engineering systems. Most of them have been simplified and assumed to a great extent. In other words, the performance functions are based on surrogate models and even empirical formulas. So far, EBRO applications in actual product design, with a time-consuming simulation model being created

for performance function, are actually quite few. After all, computational cost is a major technical bottleneck limiting EBRO applications. First, evidence theory describes uncertainty through a series of discontinuous sets, rather than a continuous function similar to a probability density function. This usually leads to a combination explosion in a multidimensional robustness analysis, and finally results in a heavy computational burden. Secondly, EBRO is essentially a nested optimization problem with performance optimization in the outer layer and robustness analysis in the inner layer. The direct solving strategy means a large number of robustness evaluations using evidence theory. As a result, the issue of EBRO efficiency is further exacerbated. Therefore, there is a great engineering significance in developing an efficient EBRO method in view of actual product design problems.

In this paper, a general EBRO model and an efficient algorithm are proposed, which provide a computational tool for robust product optimization with epistemic uncertainty. The proposed method is applied to three design problems of MEMS, in which its engineering practicability is discussed. The remainder of this paper is organized as follows. Section 1 briefly introduces the basic concepts and principles of robustness analysis using evidence theory. The EBRO model is formulated Section 2. The corresponding algorithm is proposed in Section 3. In Section 4, this method is validated through the three applications of MEMS—a micro-force sensor, a low-noise image sensor and a capacitive accelerometer. Conclusions are drawn in Section 5.

#### **2. Robustness Analysis Using Evidence Theory**

Consider that uncertainty problem is given as *f*(*Z*), where *Z* represents the *nZ*-dimensional uncertain vector, *f* is the performance function which is uncertain due to *Z*. Conventional methods [16–18] of robust optimization employ probability theory to deal with the uncertainties. The typical strategy is to consider the uncertain parameters of a problem as random variables, and thereby the performance value is also a random variable. The mean and variance are used to formulate the robustness model. In practical engineering, it is sometimes hard to construct accurate probability models due to the limited information. Thus, evidence theory [20,21] is adopted to model the robustness. In evidence theory, the frame of discernment (FD) needs to be established first, which contains several independent basic propositions. It is similar to the sample space of a random parameter in probability theory. Here, 2<sup>Θ</sup> denotes the power set of the FD (namely Θ), and 2<sup>Θ</sup> consists of all possible propositions contained in Θ. For example, for a FD with the two basic propositions of <sup>Θ</sup><sup>1</sup> and <sup>Θ</sup>2, the corresponding power set is 2<sup>Θ</sup> <sup>=</sup> {∅, {Θ1}, {Θ2}, {Θ1, <sup>Θ</sup>2}}. Evidence theory adopts a basic probability assignment (BPA) to measure the confidence level of each proposition. For a certain proposition *A*, the BPA is a mapping function that satisfies the following axioms:

$$\begin{cases} 0 \le m(A) \le 1, \forall A \in 2^{\Theta} \\ m(\Phi) = 0 \\ \sum\_{A \in 2^{\Theta}} m(A) = 1 \end{cases} \tag{1}$$

where if *m*(*A*) ≥ 0, *A* is called a focal element of *m*. The BPA of *m*(*A*) denotes the extent to which the evidence supports Proposition A. When the information comes from multiple sources, *m*(*A*) can be obtained by evidence combination rules [28]. Evidence theory uses an interval consisting of the belief measure (*Bel*) and the plausibility measure (*Pl*) to describe the true extent of the proposition. The two measures are defined as:

$$\begin{array}{l}Bel(A) = \sum\_{\mathbb{C}\subseteq A} m(\mathbb{C})\\Pl(A) = \sum\_{\mathbb{C}\cap A\neq \emptyset} m(\mathbb{C}) \end{array} \tag{2}$$

As can be seen from Equation (2), *Bel*(*A*) is the summary of all the BPA that totally support Proposition *A*, while *Pl*(*A*) is the summary of the BPA that support Proposition *A* totally or partially.

A two-dimensional design problem is taken as the example to illustrate the process of robustness analysis using evidence theory. The performance function contains two uncertain parameters (*a*, *b*), which are both considered as evidence variables. The FDs of *a*, *b* are the two closed intervals, i.e., *A* = ' *AL*, *AR*( and *B* = ' *BL*, *BR*( . *A* contains *nA* number of focal elements, and the subinterval of *Ai* = ' *BL <sup>i</sup>* , *<sup>B</sup><sup>R</sup> i* ( represents the *i*-th focal element of *A*. The definitions of *nB* and *Bj* are similar. Thus, a Cartesian product can be constructed:

$$D = A \times B = \left\{ D\_k = \left[ A\_i, B\_j \right], A\_i \in A, B\_j \in B \right\} \tag{3}$$

where *Dk* is the *k*-th focal element of *D,* and the total number of focal elements is *nA* × *nB*. For ease of presentation, assuming that *a*, *b* are independent, a two-dimensional joint BPA is obtained:

$$m(D\_k) = m(A\_i) \cdot m\left(B\_j\right) \tag{4}$$

More general problems with parametric correlation can be handled using the mathematical tool of copula functions [29].

As analyzed above, the performance function of *f* is uncertain. The performance threshold of *v* is given to evaluate its robustness. Given that the design objective is to minimize the value of *f*, the higher the trueness of Proposition *f* ≤ *v*, the higher the robustness of *f* relative to *v*. Proposition *f* ≤ *v* is defined as the feasible domain:

$$F = \{ f : f(a, b) \le v \}, (a, b) \in D\_k,\\ D\_k = \left[ A\_{i\prime} B\_{\bar{j}} \right] \subset D \tag{5}$$

Substituting *A*, *C* with *F*, *Dk* in Equation (2), the belief measure and plausibility measure of Proposition *f* ≤ *v* are expressed as follows:

$$\begin{aligned} Bel(F) &= \sum\_{D\_k \subseteq F} m(D\_k) \\ PI(F) &= \sum\_{D\_k \cap F \neq \Phi} m(D\_k) \end{aligned} \tag{6}$$

In evidence theory, the probabilistic interval composed by the two measures can describe the trueness of *f* ≤ *v*, written as *R*(*F*) ∈ [*Bel*(*F*), *Pl*(*F*)]. The accumulation of *Bel*, *Pl* needs to determine the positional relationship between each focal element and the *F* domain. As a result, the performance function extrema of each focal element must be searched. For this example, the *nA* × *nB* pairs of extremum problems are established as:

$$\begin{aligned} f\_k^{\min} &= \min\_{(a,b)\in D\_k} f(a,b) \\ f\_k^{\max} &= \max\_{(a,b)\in D\_k} f(a,b) \end{aligned} \quad \begin{cases} k = 1, 2, \dots, n\_A \times n\_B \end{cases} \tag{7}$$

where *f* min *<sup>k</sup>* , *<sup>f</sup>* max *<sup>k</sup>* are the minimum and maximum of the *k*-th focal element. The vertex method [30] can efficiently solve the problems in Equation (7) one by one. If *f* max *<sup>k</sup>* ≤ *v*, *Dk* ⊆ *F*, and *m*(*Dk*) is simultaneously accounted into *Bel*(*F*) and *Pl*(*F*); If *f* min *<sup>k</sup>* < *v*, *Dk* ∩ *F* = Φ, and *m*(*Dk*) is only accounted into *Pl*(*G*). After calculating the extrema for all focal elements, the *Bel* and *Pl* can be totaled.

#### **3. Formulation of the EBRO Model**

As mentioned above, evidence theory uses a pair of probabilistic values [*Bel*, *Pl*] to measure the robustness of the performance value related to the given threshold. However, engineers generally tend to adopt conservative strategies to deal with uncertainties in the product design process. Thus, the robustness objective of EBRO can be established as max *Bel*(*f* ≤ *v*). Meanwhile, in order to improve

product performance, the performance threshold is minimized. The EBRO model is formulated as a double-objective optimization problem:

$$\begin{aligned} \min v, \; & \max Bel(f(d, \mathbf{X}, \mathbf{P}) \le v) \\ \text{s.t. } d^l \le d \le d^u, \; & \mathbf{X}^l \le \mathbf{X} \le \mathbf{X}^u \end{aligned} \tag{8}$$

where *d* is the *nd*-dimensional deterministic design vector; *X* is the *nX*-dimensional uncertain design vector; *P* is the *nP*-dimensional uncertain parameter vector; the superscripts of *l*, *u* represent the value range of a design variable; and *X* represents the nominal value of *X*. Note that the threshold of *v* is usually difficult to give a fixed value to, while it should be treated as a deterministic design variable.

The proposed model is an improvement on the existing model [24] because it can handle more types of uncertainty, such as the perturbations of design variables resulting from production tolerances, and the variations of parameters due to changing operating conditions. As for the solving process, the EBDO involves the nested optimization of the double-objective optimization in the outer layer and the robustness assessment in the inner layer. Due to the discreteness introduced by the evidence variables, each of the robustness analyses need to calculate the performance extrema of all focal elements. Essentially, extremum evaluation is an optimization problem involving the performance function based on time-consuming simulation models, and therefore the robustness analysis bears a high computational cost. More seriously, the double-objective optimization in the outer layer requires a large number of robustness evaluations in the inner layer. Eventually, the EBDO solving becomes extremely inefficient.

#### **4. The Proposed Algorithm**

To improve efficiency, this paper proposes a decoupling algorithm of EBRO, and its basic idea is to convert the nested optimization into the sequence iteration process. Firstly, the original problem is decomposed into a series of sub-problems. Secondly, the uncertainty analysis and the deterministic optimization are driven alternately until convergence. The framework of the proposed method is detailed below.

#### *4.1. Decomposition into Sub-Problems*

Robust optimization is essentially a multi-objective problem that increases product performance at the expense of its robustness. Therefore, robust optimization generally does not have a unique solution, but a set of solutions called the Pareto optimal set [2]. It is a family of solutions that is optimal in the sense that no improvement can be achieved in any objective without degradation in others for a multi-objective problem. The Pareto-optimal solutions can be obtained by solving appropriately formulated single objective optimization problems on a one-at-a-time basis. At present, a number of multi-objective genetic algorithms have been suggested. The primary reason for this is their ability to find multiple Pareto-optimal solutions in parallel. From the viewpoint of mathematical optimization, genetic algorithms are a kind of suitable method for solving a general multi-objective optimization. However, the efficiency of a genetic algorithm is usually much lower than the gradient-based optimization algorithms, which has become the main technical bottleneck limiting its practical application [31,32]. Although a priori information is not required when using genetic algorithms, most designers have some engineering experience in practice. Therefore, for the specific problem shown in Equation(8), the robustness objective of *Bel*(*f*(*d*, *<sup>X</sup>*,*P*) <sup>≤</sup> *<sup>v</sup>*) is often handled as a reliability constraint [2,33]. In this paper, the EBRO problem is transformed into a series of sub-problems under the given target belief measures:

$$\begin{array}{c} \min v\\ \text{s.t. } Bel(f(d, \mathbf{X}, \mathbf{P}) \le v) \ge Bel\_j^T\\ d^l \le d \le d^\mu, \overline{\mathbf{X}}^l \le \mathbf{X} \le \overline{\mathbf{X}}^\mu \end{array} \tag{9}$$

where *Bel<sup>T</sup> <sup>j</sup>* represents the *<sup>j</sup>*-th target belief measure; and *Bel*(*<sup>f</sup>* ≤ *<sup>v</sup>*) ≥ *Bel<sup>T</sup> <sup>j</sup>* is the reliability constraint derived from the robust objective. In many cases, the designer may focus on the performance values under some given conditions based on the experience or quality standard. This condition is usually a certain probability of *<sup>f</sup>* ≤ *<sup>v</sup>*, namely *Bel<sup>T</sup> j* .

#### *4.2. Iteration Framework*

Theoretically, the EBDO problems in Equation (9) can be solved by existing methods [34]. However, the resulting computational burden will be extremely heavy. To address this issue, a novel iteration framework is developed, in which the uncertainty analysis and design optimization alternate until convergence.

In the *k*-th iteration, each optimization problem in Equation (9) requires the performance of an uncertainty analysis at the previous design point:

$$\text{Bel}\left(f(\mathbf{Z}) \le v\_j^{(k-1)}\right), \mathbf{Z} = (\mathbf{X}, \mathbf{P}), j = 1, 2, \dots, n\_{\text{T}} \tag{10}$$

This mainly consists of two steps, illustrated by the example in Figure 1. Step 1 is to search for the most probable focal element (MPFE) along the limit-state boundary of *<sup>f</sup>*(*Z*) <sup>≤</sup> *<sup>v</sup>* (*k*−1) *<sup>j</sup>* . The MPFE [35] is similar to the most probable point (MPP) in probability theory, which is the point with the most probability density on the limit-state boundary. Compared to other points on the boundary, the minimal error of reliability analysis can be achieved by establishing the linear approximation for the performance function at the MPP [36]. Similarly, the MPFE contains the maximal BPA among the focal elements that are crossed by the limit-state boundary. The searching process of MPFE is formulated as:

$$\begin{aligned} \max\_{D\_k} m(D\_k) \\ \text{s.t. } f(\mathbf{Z}) = v\_j^{(k-1)} \end{aligned} \qquad \begin{aligned} \max\_{\mathbf{Z}} m(D\_k) \\ \vdots \end{aligned} \qquad \begin{aligned} \max\_{\mathbf{Z}} m(D\_k) \end{aligned} $$

where *m*(*Dk*) represents the BPA of the focal element where the *Z* point is located. Note that there is a difference between *v* (*k*−1) *<sup>j</sup>* , *<sup>j</sup>* = 1, 2, ... , *nT* at each iteration step due to the minor difference of *Bel<sup>T</sup> j* . Consequently, different MPFEs may be obtained for Equation (11). However, the difference between the MPFEs is minor relative to the entire design domain. To ensure efficiency, the unique MPFE is investigated at each iteration. Equation (11) can be rewritten as:

$$\begin{aligned} \max\_{\begin{subarray}{c} D\_k \\ \mathbf{s} \in \mathcal{L} \end{subarray}} m(D\_k) \\ \text{s.t. } f(\mathbf{Z}) = v^{(k-1)} \end{aligned} \tag{12}$$

where *v*(*k*−1) represents the performance threshold that has not yet converged.

Step 2 is to establish linear approximation for the performance function at the central point <sup>+</sup> ¯ *Z M*, of MPFE:

$$L^{(k)}(\mathbf{Z}) = f\left(\mathbf{Z}^{M(k)}\right) + \left(\mathbf{Z} - \mathbf{Z}^{M(k)}\right)^{\mathrm{T}} \cdot \nabla f\left(\mathbf{Z}^{M(k)}\right) \tag{13}$$

The *L*-function is used to replace the *f*-function to calculate *Bel*, and thereby the optimization processes in Equation (7) no longer requires the calculation of any performance function.

The efficient calculation of *Bel* has been achieved in the iterative process, but the overall process of EBDO still requires dozens or even hundreds of *Bel* evaluations due to the nested optimization. To eliminate the nested optimization, a decoupling strategy is proposed similar to that in the probabilistic method [37]. At each iteration step, the reliability constraint is transformed into a deterministic constraint by constructing the shifting vector of *S*(*k*) *<sup>j</sup>* ; and then a deterministic optimization is updated and solved to obtain the current solution. In the *k*-th iteration, the deterministic optimization can be written as:

$$\begin{array}{c} \min v\\ \text{s.t.} \, f\left(d, \bar{\mathbf{Z}} - \mathbf{S}\_{\bar{\boldsymbol{\gamma}}}^{(k)}\right) \le v\\ d^{\boldsymbol{l}} \le d \le d^{\boldsymbol{u}}, \; \overline{\mathbf{X}}^{\boldsymbol{l}} \le \mathbf{X} \le \overline{\mathbf{X}}^{\boldsymbol{u}} \end{array} \right\} = 1, 2, \ldots, n\_{\boldsymbol{T}} \tag{14}$$

The shifting vector determines the deviation between the original reliability boundary and the deterministic boundary at the *k*-th iteration step. For the *j*-th problem in Equation (14), the formulation of the shifting vector is explained as in Figure 2. For convenience of presentation, the constraint contains only two evidence variables *<sup>Z</sup>* <sup>=</sup> (*a*, *<sup>b</sup>*). *<sup>F</sup>* represents the domain of *<sup>f</sup>* <sup>≤</sup> *<sup>v</sup>* (*k*−1) ¯ *Z* (*k*−1) *<sup>j</sup>* is the

*<sup>j</sup>* . previous design point, which is based on the previous equivalent boundary of *f <sup>Z</sup>* <sup>−</sup> *<sup>S</sup>*(*k*−1) *j* = *v* (*k*−1) *<sup>j</sup>* . The rectangular domain represents the FD at the previous design point. *F* represents the domain of (*k*−1)

*f* ≤ *v* (*k*−1) *<sup>j</sup>* . If the FD is entirely in the *<sup>F</sup>* domain, *Bel f* ≤ *v* (*k*−1) *j* <sup>=</sup> 100%. In Figure 2, the FD of ¯ *Z* is partially in the *<sup>F</sup>* domain, and *Bel f* ≤ *v* (*k*−1) *j* is still less than *Bel<sup>T</sup> <sup>j</sup>* . To satisfy *Bel <sup>f</sup>*(*Z*) <sup>≥</sup> *<sup>v</sup>* (*k*−1) *j* ≥ ¯ (*k*−1)

*Bel<sup>T</sup> j* , *Z* needs to move further into the *F* domain. Therefore, the equivalent boundary needs to move further toward the *F* domain. The updated equivalent boundary is constructed as follows:

$$f\left(\mathbf{Z} - \mathbf{S}\_{\dot{j}}^{(k)}\right) = v\_{\dot{j}}^{(k-1)},\ \mathbf{S}\_{\dot{j}}^{(k)} = \mathbf{S}\_{\dot{j}}^{(k-1)} + \Delta \mathbf{S}\_{\dot{j}}^{(k)}\tag{15}$$

where Δ*S*(*k*) *<sup>j</sup>* denotes the increment of the previous shifting vector. The principle for calculating <sup>Δ</sup>*S*(*k*) *j* is set as *Bel <sup>f</sup>*(*Z*) <sup>≥</sup> *<sup>v</sup>* (*k*−1) *j* ≥ *Bel<sup>T</sup> <sup>j</sup>* and is just satisfied. Thus, the mathematical model of <sup>Δ</sup>*S*(*k*) *<sup>j</sup>* is created as:

$$\min\_{\mathbf{s}} \|\|\mathbf{s}\|\|$$

$$\text{s.t. } Rel\left(f(\mathbf{Z} + \mathbf{s}) \le v\_{\mathbf{j}}^{(k-1)}\right) = Bel\_{\mathbf{j}}^T$$

Equation (16) can be solved by multivariable optimization methods [38]. To further improve efficiency, the *f* -function is replaced by the *L*-function formulated in Equation (13).

**Figure 1.** Uncertainty analysis for the performance function. FD: frame of discernment; MPFE: most probable focal element.

**Figure 2.** Formulation of the shifting vector. FD: frame of discernment.

Uncertain analysis and design optimization are carried out alternatively until they meet the following convergence criteria:

$$\left| \frac{Bel\_j^{(k)} \ge Bel\_j^T}{\frac{\nu\_j^{(k)} - \nu\_j^{(k-1)}}{\nu\_j^{(k)}}} \right| \le \varepsilon\_r \quad \begin{cases} j = 1, 2, \dots, \nu\_T \\\\ \end{cases} \tag{17}$$

where *<sup>ε</sup><sup>r</sup>* is the minimal error limit. The solutions of <sup>+</sup> *d*∗ *j* , ¯ *X* ∗ *j* , , *j* = 1, 2, ... , *nT* form the final optimal set.

The flowchart of the EBRO algorithm is summarized as Figure 3.

#### **5. Application Discussion**

In the previous sections, an EBDO method is developed for engineering problems with epistemic uncertainty. This method is especially suitable for the robust design of micro-electromechanical systems (MEMS). On one hand, unlike traditional engineering structural problems, the design of MEMS usually involves micro structure, novel materials, and extreme operating conditions, where epistemic uncertainties inevitably exist. Evidence theory is well suited to deal with such uncertainties. On the other hand, high performance and insensitivity to uncertainties are the fundamental requirements for MEMS design. Over the past two decades, robust optimization for MEMS has gradually attracted the attention of both academics and engineering practice [39–41]. In this section, this method is applied to three applications of MEMS: a micro-force sensor, a low-noise image sensor, and a capacitive accelerometer. The features of the proposed approach are investigated in terms of efficiency and accuracy. Performance function evaluations are accounted to indicate efficiency, and the reference solution is compared to verify accuracy. The reference solution is obtained by the double-loop method, where sequential quadratic programming [38] is employed for performance optimization, and a Monte-Carlo simulation [42] is used for robust assessment.

**Figure 3.** The flowchart of the proposed method. EBRO: evidence-theory-based robust optimization; BPA: basic probability assignment.

#### *5.1. A Micro-Force Sensor*

A piezoelectric micro-force sensor [43] has several advantages, including a reliable structure, fast response, and simple driving circuits. It has been extensively applied in the fields of precision positioning, ultrasonic devices, micro-force measurement, etc. Given that uncertainties are inevitable in structural sizes and material parameters, robust optimization is essential to ensure the performance of the sensor.

As shown in Figure 4, the core part of the micro-force sensor is a piezoelectric cantilever beam, which consists of a piezoelectric film, a silicon-based layer, and two electrodes. The force at the free end causes bending deformation on the beam, which drives the piezoelectric film to output polarization charges through the piezoelectric effect. The charge is transmitted to the circuit by the electrodes and converted into a voltage signal. According to the theoretical model proposed by Smits et al. [43], this voltage can be formulated as:

$$\mathcal{U} = \frac{3 \cdot d\_{31}^P \cdot S\_{11}^{Si} \cdot S\_{11}^P \cdot h \cdot h^P \cdot \left(h + h^P\right) \cdot L \cdot F}{K \cdot \varepsilon\_{33}^P \cdot w} \tag{18}$$

where

$$\begin{array}{l} \mathbb{K} = \mathbf{4} \cdot \mathbb{S}\_{11}^{\mathrm{Si}} \cdot \mathbb{S}\_{11}^{P} \cdot h \cdot \left(h^{P}\right)^{\mathfrak{T}} + \mathbf{4} \cdot \mathbb{S}\_{11}^{\mathrm{Si}} \cdot \mathbb{S}\_{11}^{P} \cdot h^{\mathfrak{T}} \cdot h^{P} \\\ + \left(\mathbb{S}\_{11}^{P}\right)^{2} \cdot h^{\mathfrak{A}} + \left(\mathbb{S}\_{11}^{\mathrm{Si}}\right)^{2} \cdot \left(h^{P}\right)^{\mathfrak{A}} + \mathbf{4} \cdot \mathbb{S}\_{11}^{\mathrm{Si}} \cdot \mathbb{S}\_{11}^{P} \cdot \left(\mathbb{S}\_{11}^{\mathrm{Si}}\right)^{2} \cdot \left(h^{P}\right)^{2} \end{array} \tag{19}$$

where *F* is the concentration force; *L*, *w* represent the length and width of the beam; *h*, *h<sup>P</sup>* denote the thickness of the silicon-base layer and piezoelectric film; *SSi* <sup>11</sup>, *<sup>S</sup><sup>P</sup>* <sup>11</sup> are the compliance coefficient of the silicon-based layer and piezoelectric film; and *d<sup>P</sup>* <sup>31</sup>, *<sup>ε</sup><sup>P</sup>* <sup>33</sup> is the piezoelectric coefficient and dielectric constant of the piezoelectric film. The constants in Equation (18) include *<sup>h</sup><sup>P</sup>* = <sup>5</sup> × <sup>10</sup>−<sup>4</sup> mm, *SP* <sup>11</sup> = 18.97 × <sup>10</sup><sup>−</sup>12m2/N, and *<sup>S</sup>Si* <sup>11</sup> = 7.70 × <sup>10</sup>−12m2/N, The structural sizes of *<sup>L</sup>*, *<sup>w</sup>*, *<sup>h</sup>* and the material parameter of *d<sup>P</sup>* <sup>31</sup>, *<sup>ε</sup><sup>P</sup>* <sup>33</sup> are viewed as evidence variables. The marginal BPAs of the variables are shown in Figure 5, and the nominal values of *d<sup>P</sup>* <sup>31</sup>, *<sup>ε</sup><sup>P</sup>* <sup>33</sup> are, respectively, 1.8 C/N, 1.6 F/m.

**Figure 4.** A piezoelectric cantilever beam.

**Figure 5.** Marginal BPAs of variables in the micro-force sensor problem. BPA: basic probability assignment; FD: frame of discernment.

In engineering, the greater the output voltage, the higher the theoretical accuracy of the sensor. Thus, *U* is regarded as the objective function. The design variables are *L*, *w* and *h*. The constraints of shape, stiffness and strength are considered, which are expressed as *η* ≥ 0.83, *δ* ≤ 2.5 μm, and *σ* ≤ 32.0 MPa, where *η* is the ratio of *w* to *h*, *δ* denotes the displacement at the free end of the beam, and *σ* denotes the maximum stress of the beam. *δ* and *σ* can be written as [43].

$$\begin{array}{c} \sigma = \frac{\frac{6 \cdot F \cdot L \cdot S\_{11}^{Si} \cdot \left( S\_{11}^{p} \cdot h + S\_{11}^{Si} \cdot h^{P} \right) \cdot \left( h + h^{P} \right)}{K \cdot w}}{\delta = \frac{4 \cdot F \cdot L^{3} \cdot S\_{11}^{Si} \cdot S\_{11}^{P} \cdot \left( S\_{11}^{P} \cdot h + S\_{11}^{Si} \cdot h^{P} \right)}{K \cdot w}} \end{array} \tag{20}$$

Due to uncertainties in the structure, *η*, *δ* and *σ* are also uncertain. Theoretically, the three constraints should be modeled as reliability constraints. To focus on the topic of robust optimization, the constraints are considered as deterministic in this example. That is, the nominal values of the uncertain variables are used to calculate *η*, *δ* and *σ*. In summary, the EBRO problem is formulated as follows:

$$\begin{aligned} \max & \text{l.o.} \max \text{Bel}(\text{l.l}(\mathbf{X}, \mathbf{P}) \ge \text{l.l}\_0) \\ \text{s.t. } \frac{w}{h} \le 0.83, \delta \left( \overline{\mathbf{X}} \right) \le 2.5 \mu \text{m}, \sigma \left( \overline{\mathbf{X}} \right) \le 32.0 \text{MPa} \\ & 0.40 \text{mm} \le \overline{L} \le 1.20 \text{mm}, 0.06 \text{mm} \le \overline{w} \le 0.10 \text{mm}, 0.04 \text{mm} \le \overline{h} \le 0.10 \text{mm} \end{aligned} \tag{21}$$

where *X* = (*L*, *w*, *h*), *P* = *dP* <sup>31</sup> ,*ε<sup>P</sup>* 33 ; *U*<sup>0</sup> represents the performance threshold, which is set as the deterministic design variable.

The steps to solve this problem using the proposed method are detailed below. Firstly, according to the marginal BPAs of the five variables in Figure 5, the joint BPAs of the focal elements (85 = 32768) are calculated by Equation (4). Secondly, Equation (21) is converted into a series of sub-problems, which are expressed as:

$$\begin{array}{c} \max \mathcal{U}\_0\\ \text{s.t. } Bel(\mathcal{U}(\mathbf{X}, \mathbf{P}) \ge \mathcal{U}\_0) \ge Bel\_j^T\\ \frac{w}{\hbar} \le 0.83, \ \delta \begin{pmatrix} \mathbf{X} \\ \mathbf{X} \end{pmatrix} \le 2.5 \mu \text{m}, \ \sigma \begin{pmatrix} \mathbf{X} \\ \mathbf{X} \end{pmatrix} \le 32.0 \text{MPa} \end{array} \tag{22}$$
  $0.40 \text{mm} \le \overline{L} \le 1.20 \text{mm}, 0.06 \text{mm} \le \overline{w} \le 0.10 \text{mm}, 0.04 \text{mm} \le \overline{h} \le 0.10 \text{mm}$   $Bel\_j^T = (80\%, 85\%, 90\%, 95\%, 99.9\%)$ 

where *Bel<sup>T</sup> <sup>j</sup>* represent a series of target *Bel* for the proposition of *U* ≥ *U*0, which are given by the designer according to engineering experience or quality standards. Thirdly, the iteration starts from the initial point of *L* (0) , *w*(0) , *h* (0) , *U*(0) 0 = (0.60 mm, 0.08 mm, 0.06 mm, 35.6 mV), where *L*(0), *w*(0), *h*(0) are selected by the designer and *U*(0) <sup>0</sup> is calculated by Equation (18). At each iteration step, the approximate function of *U* is established as Equation (13), and then 10 numbers of <sup>Δ</sup>*S*(*k*) *<sup>j</sup>* are obtained through Equation (16). Correspondingly, the 10 optimization problems as Equation (14) are updated. By solving them, the optimal set in the current iteration is obtained. After four iteration steps, the optimal set is converged as listed in Table 1. The results show that the performance threshold decreases gradually with increase in *Bel.* In engineering, a designer can intuitively select the optimal design option from the optimal set by balancing the product performance and robustness. In term of accuracy, the solutions of the proposed method are very close to the corresponding reference solutions, and the maximal error is only 2.5% under the condition of *Bel<sup>T</sup> <sup>j</sup>* = 95%. In efficiency, the proposed method calculates performance function only 248 times, and the computational cost is much less than that of evolutionary algorithms [31]. From a mathematical point of view, it is unfair to compare the efficiency of the proposed method with the evolutionary algorithms. From the view of engineering

practicality, however, the solutions of the proposed method may help the designers create a relatively clear picture of the problem with high efficiency and acceptable accuracy.


**Table 1.** Optimal set of the micro-force sensor problem.

#### *5.2. An Ultra-Low-Noise Image Sensor*

Recently, a type of ultra-low-noise image sensor [44] was developed for applications requiring high-quality imaging under extremely low light conditions. Such a type of sensor is ideally suited for a variety of low light level cameras for surveillance, industrial, and medical applications. In application, the sensor and other components are assembled on a printed circuit board (PCB). Due to the mismatch in the thermal expansion coefficient of the various materials, thermal deformation occurs on the PCB under the combined action of self-heating and thermal environment. As a result, the imaging quality of the sensor is reduced. Moreover, to acquire more image information under low-light conditions, the sensor is designed in a large format. Thus, the imaging quality is more susceptible to deformation. This issue has become a challenging problem in this field and needs to be solved urgently.

A robust optimization problem is considered for the camera module as Figure 6, in which the image-sensor-mounted PCB is fastened with the housing. The sensor is designed as a 4/3-inch optical format and features an array of five transistor pixels on a 6.5 μm pitch with an active imaging area of 2560(H) × 2160(V) pixels. It delivers extreme low light sensitivity with a read noise of less than 2.0 electrons root mean square (RMS) and a quantum efficiency above 55%. In order to analyze the thermal deformation of the sensor under the operating temperature (20 ◦C ∼ 45 ◦C), the finite element model (FEM) is created as shown in Figure 7, in which the power dissipation *P* = (*P*1, *P*2) of the codec chip and the converter is given as 1.2 W and 0.2 W, according to the test data. It can be observed that a certain deformation appears on the sensor die, and the peak–peak value (PPV) of the displacement response achieves about 3.0 μm. Consequently, the image quality of the sensor will decrease. To address the issue, the design objective is set to minimize the PPV, and the design variables of *X* = (*X*1, *X*2, *X*3) are the normal positions of the PCB-fixed points. In engineering, manufacturing errors are unavoidable and power dissipation fluctuates with changing loads, and thereby *X* and *P* are treated as evidence variables. Their BPAs are summarized on the basis of limited samples, as listed in Table 2. This robust optimization is constructed as follows:

$$\min \upsilon\_{\prime} \max \operatorname{Bel}(PPV(\mathbf{X}\_{\prime}\mathbf{P}) \le \upsilon) \tag{23}$$

**Figure 6.** The camera module (**a**) with an ultra-low-noise image sensor (**b**). PCB: printed circuit board.

**Figure 7.** The finite element model (FEM) of the image-sensor-mounted PCB. PCB: printed circuit board.


**Table 2.** The marginal BPA of variables in the image sensor problem.

As mentioned above, the performance function of PPV is implicit and based on the time-consuming FEM, which consists of 88,289 8-node thermally coupled hexahedron elements. The computational time for solving the FEM is about 0.1 h, if using a computer with the i7-4710HQ CPU and 8 G of RAM. To realize the parameterization and reduce the computational cost of obtaining reference solutions, a second-order polynomial response surface is created for the performance function by sampling 200 times on the FEM.

$$\begin{array}{l} \text{PPV} = 2.396 - 9.924X\_1 - 6.495X\_2 + 14.178X\_3 + 0.311P\_1 - 0.226P\_2 + 8.564X\_1^2 + 16.960X\_2^2 \\ + 14.104X\_3^2 - 0.019P\_1^2 + 0.794P\_2^2 - 1.540X\_1X\_2 - 13.168X\_1X\_3 - 13.822X\_2X\_3 - 0.074P\_1P\_2 \end{array} \tag{24}$$

In order to analyze the efficiency of the proposed method for problems with different dimensional uncertainty, three cases are considered: only *P* is uncertain in Case 1; only *X* is uncertain in Case 2; and both of them are uncertain in Case 3. The initial design option is selected as ¯ *X* 0 <sup>=</sup> (0.10 mm, 0.10 mm, 0.10 mm), and *<sup>v</sup>*<sup>0</sup> <sup>=</sup> 2.92 <sup>μ</sup>mis obtained by *PPV*<sup>+</sup> ¯ *X* 0 , ¯ *P* , . After giving *Bel<sup>T</sup> <sup>j</sup>* = (85%, 95%, 99.9%), the original problem is converted into three sub-problems, as in Equation (9). They are solved by the proposed method and the double-loop method; all results are listed in Table 3. Firstly, the results of the proposed method and the reference solutions are almost identical for all cases, and thus the validity of the results is presented. Secondly, each of the cases converges into a stable optimal set after three or four iteration steps. For this problem, the convergence of the proposed method is little affected by the number of uncertain variables. Thirdly, the performance function evaluations (*NF*) increase with the increasing dimensional number of uncertainties, while overall the efficiency of the proposed method is relatively high. Case 3 is taken as an example. Even if the FEM is called directly by EBRO, *NF* = 198 means a computational time of only about 20 h.


**Table 3.** The optimal set of the micro-force sensor problem.

#### *5.3. A Capacitive Accelerometer*

The capacitive accelerometer [45] has become very attractive for high-precision applications due to its high sensitivity, low power consumption, wide dynamic range of operation, and simple structure. The capacitive accelerometer is not only the central element of inertial guidance systems, but also has applications in a wide variety of industrial and commercial problems, including crash detection for vehicles, vibration analysis for industrial machinery, and hovering control for unmanned aerial systems.

Most of the capacitive accelerometers consist of two main modules: the sensing structure and the signal processing circuit; the former plays a critical role in the overall product performance. The sensing structure in this example, as in Figure 8, mainly includes five parts: a fixed electrode, a movable electrode, a coil, a counter weight, block 1, and block 2. The material they are made of is listed in Table 4. The capacitance between the two electrodes varies with the vertical displacement of the movable plate under the excitation of acceleration, which can be clearly presented through the finite element simulation, as in Figure 9. The nodes in the effective area on the movable electrode have been offset relative to the original position under the excitation of acceleration. The increment of capacitance is expressed as [44]:

$$\begin{array}{l} \Delta \mathbb{C} = \varepsilon \cdot A\\ A = \sum\_{i=1}^{l} \left( \frac{S\_i}{h + \delta\_i} - \frac{S}{h} \right) \end{array} \tag{25}$$

where *ε* is the dielectric constant; h represents the original distance between the electrodes; *S* denotes the effective area on the movable electrode; *δ<sup>i</sup>* is the displacement response of the *i*-th node, and *Si* is the area of corresponding element. Note that the performance function of *A* is based on the FEM, which contains 114,517 8-node thermally coupled hexahedron elements in total, and it takes about 1/3 h to solve each time, when using a personal computer. The displacement of the movable electrode, in addition to the response to acceleration, may be caused by varying ambient temperature. This can be found from the simulation result in Figure 9, where the load is changed from acceleration to varying temperature. Reducing the effect of thermal deformation on accuracy has become a problem that must be faced in the design process. Therefore, the EBRO model of the capacitive accelerometer is formulated as:

$$\min v, \max \text{Bel}\left(f = \frac{A(\mathbf{X}, \mathbf{a})}{\Delta T} \le v\right) \tag{26}$$

where *f* represents the sensitivity of error to temperature at 35 ◦C; *v* denotes the performance threshold; and *X*, *P* denote the design vector and parameter vector. The components of *X*, *α* are the structural sizes, as shown in Figure 8, where *α*1, *α*2, *α*<sup>3</sup> are the mounting angle of the counter weight, block 1 and block 2, respectively. The value ranges are given as 6.0 mm ≤ *X*1, *X*<sup>2</sup> ≤ 12.0 mm, and the

nominal values are *α<sup>i</sup>* = 0 mm, *i* = 1, 2, 3. All of them are uncertain variables, respectively caused by machining errors and assembly errors. According to the existing samples, the marginal BPAs are listed in Figure 10.

**Figure 8.** The sensing structure of a capacitive accelerometer.

**Table 4.** Material of accelerometer parts.


**Figure 9.** The FEM of the capacitive accelerometer.

After being given a series of *Bel<sup>T</sup> <sup>j</sup>* = 80%, 85%, 90%, 95%, 99.9%, Equation (26) can be rewritten as Equation (27):

$$\begin{aligned} \min v\\ \text{s.t.} \operatorname{Bel} \left( f = \frac{A(\mathbf{X}, \mathbf{a})}{\Delta T} \le v \right) \ge \operatorname{Bel}\_j^T \end{aligned} \quad \begin{aligned} \text{s.t.} \quad & j = 1, 2, \dots, 6 \end{aligned} \tag{27}$$

For easy reproduction of results, the response surface of *A*(*X*, *α*) is constructed as follows:

$$\begin{array}{l} A = 1.207X\_1^2 - 0.430X\_1X\_2 - 18.06X\_1 + 1.004X\_2^2 - 13.974X\_2 \\ + 100.9a\_1^2 - 9.0a\_1(a\_1 - 1) + 89.0a\_2^2 - 7.2a\_2 + 40.9a\_3^2 - 6.7a\_3 + 144.0 \end{array} \tag{28}$$

Next, the EBRO is performed by the proposed method and the double-loop method. The initial design point is selected as ¯ *X* (0) = (9.6mm, 9.6*mm*), and *f* + ¯ *X* (0) , = 0.565μm/◦C. All results are given in Table 5. The proposed method converges the optimal set after four iteration steps. Each element of the optimal set is very close to that of the reference solution. This indicates to some extent the convergence and accuracy of the proposed method. As for efficiency, the performance function evaluations of the proposed method are done 171 times. Compared to the double-loop method (12,842 times), the efficiency of this method has a definite advantage. Given that hundreds of simulations or dozens of hours of computation are acceptable for most engineering applications, it is feasible to directly call the time-consuming simulation model when performing the EBRO in practice.


**Table 5.** The optimal set of the accelerometer problem.

On the other hand, the proposed method provides six design options under different robustness requirements in this example. The higher the robustness, the greater the performance threshold. From a designer's point of view, choosing a lower yield (i.e., *Bel* = 81.8%) means that a higher cost and less error (i.e., *v* = 0.244 μm/◦C) are introduced by temperature varying. Usually, the final design option is selected from the optimal set after balancing the cost and performance of the accelerometer. Objectively speaking, the proposed method does not provide a complete Pareto optimal set, but rather solutions under the given conditions. However, for the design of an actual product, the information of *Bel<sup>T</sup> <sup>j</sup>* can usually be obtained on the basis of the engineering experience or the quality standards. Therefore, the proposed method is suitable for most product design problems.

#### **6. Conclusions**

Due to inevitable uncertainties from various sources, the concept of robust optimization has been deeply rooted in engineering designs. Compared to traditional probability models, evidence theory may be an alternative to model uncertainties in robustness optimization, especially in the cases of limited samples or conflicting information. In this paper, an effective EBRO method is developed, which can provide a computational tool for engineering problems with epistemic uncertainty. The contribution of this study is summarized as follows. Firstly, the improved EBDO model is formulated by introducing performance threshold as a newly-added design variable, and this model can handle the uncertainties involved in design variables and parameters. Secondly, the original ERBO is transformed into a series of sub-problems to avoid double-objective optimization, and thus the difficulty of solving

is reduced greatly. Thirdly, an iterative strategy is proposed to drive the robustness analysis and the optimization solution alternately, resulting in nested optimization in the sub-problems achieving decoupling. The proposed method is applied to the three MEMS design problems, including a micro-force sensor, an image sensor, and a capacitive accelerometer. In the applications, the finite element simulation models and surrogate models are both given. Numerical results show that the proposed method has good engineering practicality due to comprehensive performance in terms of efficiency, accuracy, and convergence. Also, this work provides targeted engineering examples for peers to develop novel algorithms. In the future, the proposed method may be extended to more complex engineering problems with dynamic characteristics or coupled multiphysics.

**Author Contributions:** Conceptualization, Z.H. and S.D.; Data curation, J.X.; Formal analysis, J.X.; Funding acquisition, T.Y.; Investigation, Z.H.; Methodology, Z.H.; Project administration, T.Y.; Resources, S.D. and F.L.; Software, J.X.; Validation, S.D. and F.L.; Visualization, J.X.; Writing–original draft, Z.H.; Writing–review & editing, S.D.

**Funding:** This research was supported by the Major Program of National Natural Science Foundation of China (51490662); the Educational Commission of Hunan Province of China (18A403, 17A036, 17C0044); and the Natural Science Foundation of Hunan Province of China (2016JJ2012, 2017JJ2022, 2019JJ40296, 2019JJ40014).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
