Finite Differences for Recovering the Plate Proﬁle in Electrostatic MEMS with Fringing Field

: Global existence and uniqueness conditions for a dimensionless fourth-order integro-differential model for an electrostatic-elastic MEMS device with parallel plates and fringing ﬁeld contribution were recently achieved by the Authors. Moving from this work, once the dielectric proﬁle of the deformable plate according with experimental setups has been assigned, some technical conditions of applicability for the intended use of the device as well as the mechanical tension of the deformable plate are presented and discussed. Then, highlighting the link between the fringing ﬁeld and the electrostatic force, ﬁnite differences were exploited for recovering the deformable plate proﬁle according both global existence and uniqueness conditions. Moreover, the inﬂuence of the electro-mechanical properties of the deformable plate on both the numerical approach and on the intended uses of the device is discussed, comparing the results with experimental setups regarding pull-in voltage and electrostatic pressure.


Introduction
In the last decade, a high synergy between the skills of scientific research and the world of industry has emerged for the development of technologically transferable physicalmathematical models capable of formalizing the increasingly performing behaviors of micro-electro-mechanical systems (MEMS) [1][2][3]. Today, these devices are considered "intelligent" because they combine electrical, electronic, mechanical, optical and other behaviors, managing highly complex industrial processes [4][5][6][7][8][9]. Among these, electrostatic MEMS with parallel metal plates are widely used devices in the industry, as they are easy to construct and exhibit high versatility [1,10,11]. Today, theoretical modeling of these devices is so advanced as to evaluate any correspondences between mechanical stresses of the deformable elements and the intended use of the device without neglecting the possibility of carrying out functionality tests in operation, which would normally require the destruction of the device itself [12,13]. Whatever the intended use of the device, it is necessary that there are no electrostatic discharge phenomena inside it, caused by contact between deformable and fixed elements, which would damage the device itself [14][15][16]. Therefore, it appears necessary to reduce, as much as possible, the physical causes, such as the fringing field, to produce an excessive approach between deformable and fixed elements [17][18][19]. The fringing field, which strongly depends on the length/width ratio of the device, produces important effects on the bending of the lines of force of the electric field, E, inside it, manifesting this influence near the edges; however, in the center, this effect is almost nil [20,21]. Among these models is an important dimensionless integrodifferential model of the fourth order, studied with a high level of attention to detail in [22], which, being endowed with global existence results, opened interesting scenarios regarding future developments. However, due to its formulation, the dimensionless model studied in [22] is not very suitable for technology transfer, since some of its analytical elements do not correctly model the behavior of some elements present in the device. This impasse was overcome in [17] where, starting from [22], a new formulation of the dimensionless model was presented and studied, more in keeping with the industrial reality of the MEMS devices produced, which also considers the effects due at the fringing field. This dimensionless model takes the form: (1) where • Ω represents the smooth bounded domain (MEMS device); • u : Ω → R defines the profile of the deforming plate; • f : Ω → R + is a bounded function which considers the dielectric properties of the material constituting the deformable plate; • λ is a positive dimensionless parameter depending on the applied voltage, V, defined as in (3) which represents the ratio of a reference electrostatic force to a reference elastic force; • β is a positive dimensionless parameter which takes into account the stiffness of the deformable plate; • γ is a positive dimensionless parameter which takes into account the stretching effect in the deformable plate; • χ, defined as in (8), is a positive dimensionless parameter which takes into account the non-local dependence of V on the solution due to a possible non-uniform electric charge distribution; • δ is a positive dimensionless parameter (usually constant) that weighs the effects due to the fringing field. • λδ|∇u(x)| 2 takes into account, according to the Pelesko-Driscoll theory, the effects due to the fringing field [17,23].
In the past, simplified physical-mathematical models, containing terms due to the fringing field but poorly adhering to industrial realities, have been proposed, studied, and validated [23,24]. Model (1) appears interesting for many reasons. On one hand, λδ|∇u(x)| 2 allows easy software/hardware implementations, while on the other hand, it is voltage controllable (presence of λ that, as specified in (3), depends on V), even if the dependence of δ on the applied external voltage has not yet been highlighted (satisfactory limitations of δ have been obtained for membrane MEMS devices [17,23,24]). However, even if the stability of (1) has not yet been studied (for which any unstable positions of the deformable plate would risk causing unwanted electrostatic discharges), the device is controlled in voltage (presence of λ). In such cases, appropriate destinations of use (such as biomedical ones, at reduced voltage) reduce the risk of electrostatic discharges. We also observe that, in (1), there are terms due to stiffness and self-stretching (identifiable by the parameters β and γ) of the deformable element to take into account the mechanical phenomena of fatigue due to the continuous rising-lowering of the deformable element. Furthermore, with the deformation of the plate, significant variations of electrostatic capacitance occur inside the MEMS, making it necessary to express its dependence on the geometric parameters of the device and on an auxiliary electrostatic capacitance (C f ) capable of opposing sudden voltage variations applied (presence of χ in (1)) [1,17]. Finally, f (x) in (1) (present in the "capacitive" term of the model) models the dielectric properties of the plates.
Equation (1) does not allow to obtain u(x) explicitly, thus in [17], important conditions of global existence and uniqueness for the solution have been obtained by opening possible scenarios for the numerical reconstruction of u(x) profiles without representing ghost solutions (i.e., approximate solutions not satisfying the conditions of existence and uniqueness for (1)). Therefore, the main aims of this paper can be summarized in the following items: • Exploiting the conditions of global existence and uniqueness for (1), fundamental results are discussed highlighting that, starting from the reference energy state, the solution is unique ensuring that the profile of the deformable plate is uniquely determined without highlighting anomalies in the curvature of the deformable plate (especially in the start-up phase). This item is very important because, compared to what is already known in the literature, the existence and uniqueness of the solution for (1) is studied starting from a reference energy state, thus linking the important properties of existence and uniqueness with the minimum energy level necessary to move the deformable plate. • After selecting the dielectric profile, f (x) (in order to satisfy binding physical requirements allowing comparison with experimental setups known in the literature), we discuss important industrial implications in terms of limitations both for the applied voltage, V, and for the mechanical tension of the deformable plate at rest, T, providing both a graphical differentiation of areas allowed in the plane (V, T) and reference energy states. So, unlike what is already known in the literature (where in the dimensionless mathematical models f (x) is set equal to the unit), in this paper we provide a useful criterion for selecting the intended use of the device, based on (V, T), starting at f (x). • We study how the fringing field affects the electrostatic force in the device, obtaining an increase for it depending on both V (which affects the intended use of the device) and T (which influences the choice of material of the deformable plate). In the literature, there are no increases for either V or T, thus failing to provide maximum admissible values for both. This item of the present work fills this gap. • Following a simplification in the model that does not affect the goodness of the results obtained, the numerical recovering of the deformable plate profile was obtained by the finite difference method (implemented in MatLab® R2019a running on an Intel Core 2 CPU at 1.45 GHz) "gold standard" for model as (1), under different operating conditions, giving results compatible with the conditions of global existence and uniqueness of the solution (thus not representing ghost solutions). • We also provide an effective criterion for choosing the intended use of the device (strongly linked to V) starting from the choice of the material constituting the deformable plate (strongly dependent on T) and vice versa very useful for any industrial applications (unlike the scientific works known in the literature where such a criterion has never been elaborated). • Furthermore, how the electro-mechanical properties of the deformable plate affect the numerical profile recovering is studied and discussed also selecting the most important intended uses of the studied device. • Finally, important comparative results with experimental setups concerning pull-in voltage and electrostatic pressure are presented.
The reminder of the paper is organized as follows: Beginning with the description of the device and analyzing both behavior and analytical models (Section 2), the most important results of global existence and uniqueness for the solution for (1) are summarized in Section 3. Next, after selecting the dielectric profile of the deformable plate as specified in Section 4 (according to several industrial applications), Section 5, starting the abovementioned inequality governing the global existence and uniqueness for (1), discussions and important limitations for both the intended use of the device and the mechanical properties of the deformable plate are deduced, respectively. Once Section 6 formalizes the link between the fringing field and electrostatic force in the device, a numerical approach is proposed and applied in Section 7 to recover the profile of the deformable plate, the results of which are presented and discussed in Section 8. Section 9 discusses an interesting limitation for the mechanical tension of the deformable plate, while Sections 10 and 11, respectively, discuss some details concerning the possible influence of the properties of the deformable plate on the numerical procedure and the possible intended uses of the device. Finally, Section 12 presents the results concerning the pull-in voltage and the electrostatic pressure. Further conclusions and perspectives on these issues complete this work.

The Behavior of the MEMS Device: The Analytical Model with Fringing Field
The electrostatic MEMS studied in this paper is an elastic system consisting of two parallel plates (of the same material and equal thickness) of which the upper one (nondeformable) is at V > 0 electric potential, while the lower plate (bound to the edges and at electric potential reference, V = 0) deforms towards the top plate without touching it (to avoid unwanted electrostatic discharges) (for further details, ref. [17]). Figure 1 displays a 3D schematic of the device (the deformable plate is in its rest position). As already detailed in [17], if V is the drop voltage, the load function F(x) (which establishes how the deformable plate is stressed when V is applied, determining E which locally activates an electrostatic pressure deforming the plate) can be written as [17] with where 0 is the permittivity of the free space, d is distance between the plates, and L is the length of the MEMS. (2), according the Pelesko's procedure [1,25], is necessary to formulate where K 1 (x) and K 2 (x) are specific weight functions as below defined. It should be noted that (2) cannot electrically control the MEMS because the applied V must be controlled in order to avoid sudden deformation of the deformable plate. Then, a capacitive control de-vice (such as the one hypothesized in [17]) can be used to avoid this drawback. Particularly, considering the series of a source voltage, V s , and a capacitor, C f , it follows that achieving Therefore, setting from (7), we can write V V s 2 control action Finally, (2) becomes (10) so that (4) assumes the final form We also observe that χ mainly depends on the capacity [1] and, furthermore, χ ∈ [0, 1) [17], thus avoiding a dangerous bifurcation phenomena which can make the device highly unstable (particularly if χ → 1 − the rupture of the device can takes place [17]). Finally, ∆u(x) being an operator indicative of the curvature assumed by the plate during deformation [26], it can be proved that K 1 (x)∆ 2 u(x), in elastic regimes, takes the form [17]: and being K 1 (x) = D L 2 T ≈ 1 [1,17,26] (D, flexural rigidity of the material constituting the deformable plate) achieves the following model: studied in [1,25], highlighting important results regarding the bifurcation. Finally, to take into account the effects due to the fringing field (physical phenomenon, heavily dependent on the ratio L/d, that manifests itself with the curvature of the E lines of force, more evident at the edges of the device and negligible in its center), we use the Pelesko-Driscoll formulation [23] that quantifies these effects through the addend λδ|∇u(x)| 2 (for details, [17,22,23]), achieving model (1). Table 1 lists all symbols present in (1).  (13).

Symbol
Meaning parameter that weighs the deformation energy of the deformable plate λ ratio of a reference electrostatic force to a reference elastic force γ stretching parameter We observe that, under the effect of the fringing field, the electrostatic capacity of MEMS varies considerably and can also be formulated through empirical approaches [1,18,25]. Recently, many MEMS models with fringing field have been studied according to the Pelesko-Driscoll theory [24] (which, however, propose studies on highly simplified models). Obviously, even if (1) models a detailed behavior of the MEMS under study, it does not consider further non-linearities (local mechanical stresses, bifurcations), and thus introduces a degree of uncertainty. However, to take into account local mechanical stresses would mean to introduce in the model a terms depending, for example, on the Piola-Kirchhoff tensor which, usually, exhibits discontinuities along the surface of separation between different zones of the deformable plate. This additional term would reduce the amplitude of the deformation of the membrane profile by approximately 5%. Therefore, not considering the contribution due to this term allows us to affirm that the numerical recovering we will obtain will be overestimated but for the sake of safety since we are sure that the real deformation of the deformable plate is below that reconstructed numerically (certainty that the deformable plate does not touch the upper wall of the device). Furthermore, considering λ < λ * (λ * , pull-in voltage) assures us that bifurcation phenomena cannot take place.

Remark 1.
It is worth noting that model (1) represents of a generalization of the model extensively studied in [27]. In fact, when β = γ = χ = 0 (i.e., in the presence of a deformable membrane), it is easy to achieve which represents the model studied in [27]. (3), it is easy to deduce that by selecting the type of device (that is, by choosing T), its intended use is determined (λ is proportional to V 2 ). The converse is also valid: once the use of MEMS has been identified, the material constituting the deformable plate can be selected (in other words, the intended use of the device sets an important limitation for λ). Therefore, as already verified in the past for other devices [28], the link between the intended use of the device and the mechanical properties of the deformable plate is confirmed. (2), as already mentioned, represents the type of load on the deformable plate applying V externally (V determines E between the two plates by generating the electrostatic force which, per surface unit, is transformed into electrostatic pressure acting on the deformable plate). Then, the need for (2) to depend on V is confirmed.

Remark 4.
We observe that the dielectric properties of the material constituting the deformable plate are decisive for the operation of the device when the deformation is the maximum allowed [1,25].
In fact, indicating with d * the minimum distance allowed between the two plates, it follows that u(x) < 1 − d * (this is due to the fact that, physically, the deformable plate must not touch the upper plate and, mathematically, in (1) u(x) < 1.). If u(x) = 1 − d * , from the equation of model (1), it follows that f (x) = 0, so that (2) is canceled. Industrially, this is highly important, because f (x) contributes to determine the electrostatic load in the device. Therefore, (1) represents a good formulation for the device, since it takes into account f (x) (as already experimentally tested for similar devices [29]).

Remark 5. Equation
(1) models electrostatic MEMS for industrial applications that require small displacements of the deformable plate (u(x) ≈ 10 −6 ). In fact, if A d 2 (d the distance between the two plates at rest, A the surface of the plates) the effects due to the fringing field would be negligible which does not fit well with the hypothesis of small displacements. This confirms that (1) is a good model for industrial developments.

Remark 6.
Unlike other analytical models [29,30], Equation (1) was built on very simple device geometry; this was necessary to carry out the analytical study [17] obtaining results which, even if they have not yet achieved an experimental confirmation, have provided a significant theoretical contribution.
We note that Ω (deformable plate), industrially, has dimensions of the order of 10 −6 . Therefore, the diameter of the domain can be considered 1 (as predicted by Theorem 1). Furthermore, the constraint N < 4, even if more stringent than N < 8 [31], does not affect the industrial validity of the result (which requires 3D formulations).
The fact that the result requires Then, for each point of the deformable plate, | f (x)| must be bounded. So, (15) makes sense because f (x) represents the dielectric profile of the deformable plate which must be measurable (and in any case positive and bounded). Furthermore, Theorem 1 dictates that ∃λ * ∀λ ∈ (0, λ * ) : λ < λ * because λ * , is pull-in voltage (i.e., λ such that there are no solutions for λ ≥ λ * ). Finally, the continuity of higher order curvatures (u ∈ H 4 (Ω)) is required by imposing that the deformable element, during its movement, must not undergo substantial deformations such as to make the device unusable (especially when fatigue phenomena caused by prolonged and continuous exploitation take place [26]). Obviously, this condition is more restrictive with respect to the conditions required when membrane MEMS are considered. It is worth nothing that, even though λ is controllable by V, δ is not; this dependence is expressed by the product λδ in (1). In the future, this would require an incisive effort on the possible generalization of the Pelesko-Driskoll theory.
The proof of Theorem 1, exploiting a well-known topological result presented in [32], is divided into five points preliminarly defining the following two suitable sets: and where M, M 1 and M 2 are suitable positive constants. Moreover, if Z = L 2 (Ω), B = ∆ 2 and the following important inequalities governing both global existence and uniqueness for (1) have been achieved:

Analytical Modeling of f (x)
As proved in [33], regardless of how the dielectric constant profile is chose, the pull-in instability cannot be avoided [1,25,29]. Thus, it seems legitimate to ask whether a physically valid formulation of f (x) can influence the solutions. It is worth nothing that in [17] it has been proved that the smooth solution u(x) are symmetric, concave and, moreover, u(x) < 1. From the heuristic point of view, the most unstable area of the device is its center, where the influence of E is very strong and the influence in the supporting boundaries appears rather weak [33]. This is due to the fact that |E| is proportional to (23) so that at the center of the device, by virtue of the symmetry and concavity of the solutions, 1 − (x) assumes a minimum value (i.e., maximum value of |E|). Therefore, it seems reasonable to try to reduce |E| in the center of the device as such possible by allowing it to be stronger near the boundaries (more stable parts of the device). So, by choosing it follows that which represents the load function (see (2)). However, we speculate whether specific formulations of f (x) can affect the numerical solutions (and any multiplicity) without neglecting the pull-in voltage and the stable operating range of the device. In this work, f (x) is chosen to satisfy a symmetric power law such as In addition, a large number of experimental and industrial applications are based on (26) [1,[33][34][35][36]. Finally, for α = 0, f (x) = 1.

On the Reference Energy State
Let us first observe that the satisfaction of (18) by u 0 means that, in (x 0 , u 0 ), the deformable plate behaves like a membrane. In fact, (18) is a typical semi-elliptical nonlinear model describing an electrostatic membrane MEMS device. This is confirmed by the fact that, by definition, Therefore, β = 0 means that the energy accumulated by the deformable plate is zero (reference energy state). Moreover, γ = 0 implies that the effects due to stretching are nil, while χ = 0 imposes that C f → ∞ with immediate consequence that V = V s (i.e., the drop-voltage equivalent to the external voltage applied). This is in accordance with the experimental and industrial experience according to which the deformable plate, in the start-up phase, is similar to a membrane [1,25,26].
Furthermore, in [17], it has been shown that B(Y) is a neighborhood of z 0 = B(y 0 ). Therefore, there exists a ball S(z 0 , r) ⊂ B(Y), with radius r, and a neighborhood of has a single solution y : Then, there exists a neighborhood of and, starting from (31), the solution is unique in the neighborhood of x 0 . This highlights the fact that, starting from the reference energy state, it is guaranteed that the solution is unique (due also to the fact that the deformable plate behaves like a membrane [1,25,26]).

5.2.
On the Continuity of x → F(x, y 0 ) Considering the elastic deformation regime of the plate valid, the elastic force will be proportional to the increase of the surface, so that the stretching energy is writable as so that (19) (deriving, according to [32], from the study of the continuity of the mapping However, being y 0 ∈ Y, Therefore, by (26), from which || f || 2 ∞,Ω < M 4 , so that, also considering (3), (33) becomes where (once the geometry of the device and the material constituting the deformable plate have been chosen) C 1 and C 2 are constants. This means that in these cases, the continuity of x → F(x, y 0 ), is voltage controlled, highlighting that, in the reference energy state, there are no appreciable anomalous phenomena of curvature of the deformable plate. This is in line with the experimental and industrial experience according to which the deformable plate does not show evident deformations (and/or curvatures) in the start-up phase [37,38].

On the Injectibility of B On Y
The deformable plate structurally responds to the theory of Sophie-Germain Lagrange [26,37,38] according to which, qualitatively, where p(x) is the mechanical pressure, ∆ 2 = B, u(x) are admissible profiles of the deformable plate (u(x) ∈ Y) and D flexural stiffness of the plate, is in which t is the thickness of the deformable plate and E and ν are the Young and Poisson modules, respectively. Then, due to both the uniqueness of ∆ 2 and its injectivity on Y, it is legitimate to state that, given the load p(x), the distribution of the deformations of the plate is completely defined. Then, by means of constitutive laws [26], it is possible to obtain T(x) (mechanical tension of the deformable plate which, however, is a bounded function, and therefore increased by a constant) particularly useful for the choice of the material constituting the deformable plate.

An Interesting Limitation for Applied Voltage V
It seems only right to specify that in (20) In addition, being d Ω 1, it follows that d is a sufficiently wide range of values capable of verifying the (20) for a variety of devices.
Remark 7. It is worth noting that the decimal numbers obtained in both (45) and (46) are the result of numerical approximations. Figure 2 offers a graphical representation of (46) (where both T and V are no longer dimensionless), differentiating the areas allowed from those forbidden. Particularly, it is evident that the studied device allows a high number of intended uses (V ∈ [0, 100] Volt), allowing the employment of a wide range of materials for the deformable plate (T ∈ [1.968 × 10 −20 , +∞) N/m 2 ).

Remark 9.
It is worth noting that the maximum allowed voltage, as obtained in (44), is strongly dependent on the electromechanical properties of the material constituting the deformable plate.

Remark 10.
To achieve (44), χ → 1 to take into account possible incipient breakage of the device (bifurcation). Moreover, as specified in [17], all constants (M, M 1 and M 2 ) were set to 1. (1) does not allow explicit obtaining of solutions. Therefore, to obtain approximate solutions satisfying the analytical model (no ghost solutions), we rely on a numerical procedure which is considered the "gold standard" for this type of model.

How Fringing Field Affects Electrostatic Force
It is known that [17], so that the electrostatic force, f el , becomes: Moreover, if E f ringing f ield is the fringing field, one achieves Therefore, the related electrostatic force becomes (1 − u(x)) 2 (50) so that, by (3), one easily achieves Therefore, the total electrostatic force becomes: confirming that the ( f el ) TOT inside the device depends not only on δ but also on L/d (algebraic ratio indicative of the possible presence of the fringing field).
Here, ( f el ) TOT also depends on both x and T so that, increasing T, ( f el ) TOT decreases as required physically. To further confirm the applicability of (52) in industrial applications [39,40], when increasing V, ( f el ) TOT also increases. We finally note that, since 1 (1−u(x)) 2 < 1 d * 2 , and in our case (no longer dimensionless) d * ≈ 10 −9 and L ≈ 10 −6 , from (52), we achieve a limitation for ( f el ) TOT ( f el ) TOT < (24.12 · 10 38 + 21.14 · 10 −24 δ) highlighting, once again, that the increase of T causes a decrease in the effects due to the fringing field.

Remark 12.
According to Section 5.3, T(x) is a bounded function such that T = sup x {T(x)}, so that (53) is valid.

Remark 13
. While it appears that (44) does not depend on δ, (53) ensures that T limits ( f el ) TOT where δ exists.

Remark 14.
Both (52) and (53) confirm the important experimental issue according to which the effect due to E f ringing f ield depends on L/d. Particularly, if L/d 1, the effects due to the E f ringing f ield will be significantly reduced (as experimentally verified in [41,42]).

Deformable Plate Profile Recovering: A Numerical Approach
Equation (1), due to the way it is structured, is not suitable for numerical processing. Therefore, some simplifications in compliance with the conditions of global existence and uniqueness must be made.

From (1), taking into account (8) and considering that [17]
we can easily achieve so that (1) becomes Assuming that γ takes precedence over β, (57) takes its final form: which represents the simplified version of (1) numerically implementable by the finite difference approach (the "gold standard" procedure for this type of problem which enables us to avoid the risk of approximate solutions representing ghost solutions).

Remark 15.
We observe that (58), even if it is a simplified version of (1), does not invalidate the verification of all the properties tested in the previous sections (we omit the proof for reasons of space, as it retraces the five steps of the proof of Theorem 1 [17]).

Derivation of the Numerical Model: Finite Difference Approach
In this section, we derive the finite difference method for the non linear boundary values problems (58) in 2D space. In particular: where we assume v(x) = u x (x) and w(x) = u y (x). We consider the computational domain −1 ≤ x ≤ 1, −1 ≤ y ≤ 1 and use a uniform Cartesian grid consisting of grid points (x i , y j ), where x i = i∆x and y j = j∆y, for i = 0, · · · , I and j = 0, · · · , J, with ∆x and ∆y increments in x−direction and y−direction, respectively. In this context, we consider the special case where ∆x = ∆y = h, so that I = J = N number of intervals in both directions.
Let U ij be the approximation to the exact solution u(x i , y j ) at the grid points (x i , y j ). To discretize the model (59) at each grid point (x i , y j ) of the computational domain, we use the following second order centered finite differences Let V ij and W ij be the approximations to the x− and y−derivatives u x (x i , y j ) and u y (x i , y j ) of the exact solution, respectively; thus, we obtain the finite difference method for i, j = 1, · · · , N − 1, with s 1 = 1 − 0.5hγ and s 2 = 1 + 0.5hγ and where The finite difference Equations (60) at points near the boundary involve the known boundary values, which, are generally moved to the right-hand side. However, the zero boundary conditions do not contribute to the model under study. We thus obtain a system of 3(N − 1) 2 nonlinear equations in 3(N − 1) 2 unknowns, U ij ,V ij and W ij , that can be expressed in the vector form where U is the vector of the unknowns with U, V, W assigned by the natural row-wise ordering Z [2] . . .
with a block structure as shown in Figure 3. Since each finite difference equation involves at most five unknowns, each row of the matrix A has at most five non-zeros elements. with U, V, W assigned by the natural row-wise ordering Z [2] . . .
with a block structure as shown in Fig. 5. Since each finite difference equation involves 271 at most five unknowns, each row of the matrix A has at most five non-zeros elements.

273
In particular for the sub-matrix A 1 , T has a tridiagonal form Each sub-matrix A k , k = 1, · · · , 9, has a block form where each block, T, S 1 and S 2 , is a (N − 1) × (N − 1) matrix. In particular, for the sub-matrix A 1 , T has a tridiagonal form with s 1 = 1 − 0.5hγ and s 2 = 1 + 0.5hγ. S 1 and S 2 are diagonal matrices with elements equal to s 1 and s 2 , respectively. The non-zeros values of the matrices S 1 and S 2 are separated from the diagonal by N − 1 zeros, since these coefficients correspond to grid points lying above or below the central point in the stencil, and are hence in the next or previous row of unknowns. For the sub-matrix A 2 , it has with s 1 = −1 and s 2 = 1. The matrices S 1 and S 2 are null. For the sub-matrix A 3 , it has T null matrix and S 1 and S 2 diagonal matrices with element equal to s 1 = 1 and s 2 = −1.
The sub-matrices A 5 and A 9 are diagonal, with principal elements equal to −2h; the other sub-matrices are null. Next, the fsolve.m Matlab®R2019a routine, running on an Intel Core 2 CPU at 1.45 GHz , is used to solve the nonlinear algebraic system (62), with initial guess U 0 = 0.

Remark 16.
We observe that the obtained matrix is sparse and (at times) large, requiring the use of iterative techniques for solving linear systems based on Krylov sub-spaces ( [43]). However, in this paper, the dimensions of the array are not such as to require the imperative use of these methods. Furthermore, since the application under study is of an off-line type, the calculation times performed by fsolve.m do not affect the validity of the proposed procedure.

Remark 17.
The finite difference method approximates the value of the derivative of a function in a point (for which it would be necessary to know all the values of the function (therefore infinite) in a neighborhood of the same point), with an expression that takes them into account only a finite number (often very small). This is passed from the limit operation to the incremental ratio operation. This allows to transform a differential equation into an algebraic problem whose sparse matrix strongly depends on the number of values used in the approximation of the derivatives. Starting from the considerations on the orders of accuracy of the derivative approximation formulas, it is easy to prove that the error made by approximating the second derivative is a O(h 2 ). Moreover, the approximation error between the true solution u(x) of the problem (1) and the approximate one is of the second order. In fact, introducing the truncation error τ j , and exploiting the Taylor series developments, it is easy to achieve: where g must be at least C 2 . Moreover, introducing the vector of the global error, H, we can easily achieve: ||H|| ≤ ||A −1 || ||τ|| (69) and again ||τ|| → 0 (which ensure the consistency of the numerical approach) and ||A −1 || ≤ C, with C independent on h (which ensure the stability of the numerical approach).

Approximate u(x) and Ghost Solutions
Equations (19), (20) and (21) depict the conditions which numerical solutions must satisfy to avoid representing ghost solutions. Once implemented, the numerical procedure was tested for several experimental cases intended for industrial applications. Particularly, assuming f (x) = |x| 0.2 , h = 0.05, K 1 = 0.9 [33], max U ij have been achieved when δ increases (see Table 2 where it is also highlighted that max U ij satisfy (19), (20) and (21)).
The max U ij obtained shows that the hypothesis of small displacements is also confirmed in the numerical recovering of the deformable plate profile, with a slight increase due to the growth in the effects caused by the fringing field. The results shown in Table 2 refer to a deformable plate with γ = 0.1; however, different values of γ did not produce significant changes in max U ij with respect to the case γ = 0.1. This is because γ, in (1), has little impact, compared to the other addend, in weighing ∆u(x) (which governs the curvature of the deformable plate). Moreover, as δ increases, the contribution due to V increases, thus even high values of V (19), (20) and (21) are verified. In addition, the effects due to δ affect ∆ 2 u(x): there is a prevalence of the effect due to δ as sanctioned by (1), because u(x) increases its concavity when δ increases (E fringing field favors the deformation of the deformable plate). However, |∇u(x)| 2 in (1), has little effect numerically (for implementation reasons), therefore the underestimation of max U ij could be compensated for by amplifying factors to be researched experimentally.
We observe that the numerically constructed regularity of u(x) suggests that the performance of the approach is satisfactory even if the high computational complexity (with respect to other numerical approaches, [44]) would make the approach inconvenient for any real-time industrial applications (i.e., when short plate recovery times are required). However, these applications are currently not widespread on a large scale.   Figures 5-9 show some recovering of deformable plate profile as the amplitude of λ increases. From these profiles it can be deduced once again how, as this amplitude increases (and therefore as the externally applied V increases) the profile of the deformable plate rises more and more. Obviously, the values of max U ij are not such as to fear the risk of the deformable plate touching the top plate of the device.

Remark 18.
It is worth noting that in (1) λ is a parameter which, apparently, would turn out to be positive. However, the numerical recovering of the deformable plate profile required negative values of λ. This is due to the fact that λ > 0 would derive from an external V inversely polarizing the device.

A Further Limitation for T
The following results yield: Proposition 1. For (1), the following limitation for T holds: Proof. From (57), we achieve from which, by (3), (47) and (49), we achieve From (72), considering that |∇u(x)| 2 is bounded [17], therefore |∇u(x)| 2 < H (H constant). Moreover, being u(x) < 1, it is easy to achieve with H 1 = 1 is a bounded function). Therefore, by (46) and (73), we achieve (70), which represents a good limitation for T, depending on both x (due to f (x)) and δ. Moreover, (73) also depends on |E| (by θ). In addition, fixed L, plates with higher T can be destined in MEMS with plates very close to each other (reduced values of d), because high values of T force the deformable plate not to touch the plate superior. On the contrary, if fixed d, deformable plates with reduced T can be exploited in MEMS with reduced value of L.

Remark 19.
We observe that d 9 (with d = 10 −9 ) in (70) allows us to validate Remark 8. Figure 10 displays the link between T and V as governed by (73). Therefore, as qualitatively deduced in [27], T = KV 2 , so that once the intended use of the MEMS has been established, the straight line parallel to the ordinate axis (passing through (V, 0)) intercepts the curve at a point, in this case, all the points below it select T eligible. As in [27], Figure 10 can be employed to choose the intended use, beginning with the selected material.

On the Influence of the Properties of the Deformable Plate on the Proposed Procedure
As verified, obtaining the limitations for T and V depend on the electro-mechanical properties of the material constituting the deformable plate (presence of both θ and λ) influencing the concavity of u(x), highlighting that the increments of T reduce the amplitudes of u(x). On the other hand, the increase of T produces smaller deformations at the edges. Furthermore, the numerical approach depends on the above parameters which, in a certain sense, govern the convergence. Finally, as in [27], these electro-mechanical properties influence the choice of T starting from V, and vice versa. Particularly, the presence of T in the denominator in the limitation of V strongly influences ( f el ) TOT , and therefore |E| total which, through (1), affects ∆ 2 u(x) (higher order curvatures of the deformable plate); thus, the greater the T, the lower the concavity of the deformable plate (also confirming that plates that are more mechanically resistant show reduced deformations at the edges). Furthermore, λ in (19), (20), and (21) ensure that T also influences the algebraic conditions of existence and uniqueness of the solution for (1).
We highlight that the numeric procedure depends on λ, γ, and δ, which govern the convergence of the method obtaining u(x) which satisfy (19), (20), and (21). Finally, by (53), it becomes clear that both V and T act effectively on ( f el ) TOT . In other words, both the intended use of the device and the mechanical properties of the deformable plate affect the behavior of the device. To summarize, the model studied is more versatile than other simplified models where the intended use of the device imposed T and vice versa.

Possible Intended Uses of the MEMS Under Study
Industrially, the device presented here does not require particular construction requirements, and therefore the production costs would be low. Furthermore, both (45) and (46) ensure that the device is usable for a wide range of industrial applications (i.e., biomedical applications, micro-pumps for intravenous drug administration, and surgical micro-systems). Obviously, the exact value of V strictly determines the intended use of the device.
Nowadays, parallel plate devices are industrially widespread, especially when their use is integrated in further electro-mechanical devices. However, even if the electronic tests on MEMS have reached high levels of reliability, the mechanical tests still do not present similarly high levels of performance. Despite this, in recent years, the joint use of electronic and mechanical components mounted on a single chip has allowed very high performances of micro-sensors used in robotics and electronics for use, for example, in large metropolitan areas/Smart Cities.

Concerning the Pull-In Voltage and Electrostatic Pressure
In MEMS, applying V, is its deformable element. However, high values of V (beyond the "pull-in" value) generates instability, because the electrostatic forces exceed the elastic ones (a potentially destructive phenomenon). In other words, from (3), if λ > λ * (with λ * , pull-in value), (1) does not admit solutions. Conversely, if λ < λ * the solution for (1) exists.
As known in the literature, λ * has been proved both mathematically [1] and experimentally [45][46][47]. Particularly, following the approach proposed in [27], we plot the trends of u 0 (maximum deflection of the deformable plate) versus λ (starting from λ = 0, there will be a value of λ , indicated with λ * below which there are two solutions while above there is no solution, thus λ * represents the bifurcation point). Figure 11 depicts the rends of u 0 − λ, both without fringing field (δ = 0) and with it (δ = 1), respectively. Both trends highlight the superimposition with experimentally obtained bifurcation diagrams [45]. In fact, in non dimensionless conditions, the pull-in voltage, V pull-in , depending on the gap between the two plates, is displayed in Figure 12 whose trend is completely similar to those displayed in Figure 6 in [45] and Figure 6 in [48] where an experimental setup, consisting of two plates separated by a distance d (each plate was acrylic sheet coated with conductive aluminum) highlighting that (1), even if referred to a simplified geometry, in terms of pull-in voltage, has a well-established experimental confirmation in the literature.  This important experimental finding is confirmed by the fact that, from (3), it is easy to achieve which is completely analogous to (11) achieved in [45]. In addition, (74) is still analogous to (3) in [49] where a MEMS switch with perforated serpentine (Au) was designed, sim-ulated, fabricated and characterized. Not last, (74) is quite similar to (1) in [48] where an electrostatic MEMS with parallel plates was considered. We also note that, from (45), we can write: V max admissible = 23.7 · 10 7 10 21 T − 19.68.
Therefore, considering the usual values for each parameter, we achieve V pull-in ≈ V max admissible , with the consequence that the obtained V max admissible is compatible with the experimental setups known in the literature. Particularly, Table 3, after fixing the geometry of the device and varying the material constituting the deformable plate as in [49], reports the values of V max admissible achieved which result similar to those obtained in [49]. As already highlighted, the effects due to the fringing field depend on the aspect ratio L/d. Then, the pull-in voltage is also affected by this relationship. Figures 13 and 14 highlight, for aspect ratio L/d = 1, 2, the trends of V pull-in using both the numerical procedure used in [45] and with the procedure proposed here. It should be noted that the d increase produces a deviation from the trend obtained using [45], while agreeing with the experimental trend (as also highlighted in [45]). So, once again, we highlight the adherence of the theoretical results discussed in this paper (Deriving from (1)) with important experiments of electrostatic deflections.  Finally, from (51), the electrostatic pressure due to E f ringing f ield can be written as: which, combined with (74), provide the link between [p el ] E f ringing f ield and V pull-in that produce the following trend (as shown in Figure 15): Quite similar to that obtained in [50], concluding once again that (1) is in line with many important experimental results already known in the literature.

Remark 21.
We observe that, as discussed in Section 12, the results obtained are mainly comparable with experimental setups obtained on electrostatic membrane MEMS devices. This is due to the fact that the global existence and uniqueness conditions for (1) have been obtained starting from (31) (reference energy state), which models an electrostatic membrane MEMS device with fringing field. Furthermore, the hypothesis of small displacements, together with the fact that β and γ in (1), do not vary during the deformation, the results for (1) agree with the experimental setups for membrane devices. In the near future, we hope to use β and γ variables to extend the comparison to further experimental setups.

Remark 22.
Regarding the possible experimental confirmation of the results achieved for T, Figure 12 depicts the link obtained between the pull-in voltage and the distance between the plates of the device. This link was found to be superimposable to experimental links known in the literature. Therefore, by (74), also T, indirectly, is verified experimentally.

Conclusions and Perspectives
In this paper, a dimensionless fourth-order integro-differential model known in the literature was considered while modeling the behavior of an electrostatic MEMS device with metal parallel plates in which the effects due to the fringing field are modeled, according to the Pelesko-Driscoll theory, by means of an easily implementable addend via software/hardware whose weight can be controlled in voltage (parameter λ). This model was formalized so that the device can be controlled in voltage by a suitable capacitive control circuit. Starting from the conditions of global existence and uniqueness, the analytical results were discussed, highlighting the most important and captivating implications for our industrial realities. In particular, the choice of a particular dielectric profile of the deformable plate (satisfying mandatory physical requirements) led to obtaining important limitations for the external electrical voltage, for the mechanical tension of the deformable plate, and for the total electrostatic force (including the contribution due to the fringing field) highlighting relevant adherence with experimental evidence known in the literature. It is undoubtedly worth noting that both the limitations for V and for T allow you to reconstruct u(x) profiles (model solutions) in complete safety. In other words, once a particular device is fixed (i.e., by fixing T of the material constituting the deformable plate), the external V applied satisfying the above limitations, it will allow u(x) profiles of the deformable plate complying with the safety standards according to which the deformable element must not touch the upper wall of the device. Conversely, by selecting the intended use of the device (i.e., V satisfying the above limitations), it is possible to choose a particular material for the deformable plate (whose T satisfies the relative limitation) capable of guaranteeing solutions u(x) of high security. Since the model does not allow the explicit recovering of the deformable plate, a numerical technique based on finite differences (the "gold standard" for this particular type of analytical models) was used to obtain profiles that are fully compatible with the conditions of global existence and uniqueness. Furthermore, a simple criterion for choosing the intended use has been provided, beginning with the mechanical tension of the deformable plate and vice versa, which can be particularly useful for industrial applications. In particular, the high value of V max admissible allows the device studied to be used also in intended uses where V would be high (industrial applications that require up to Class 3B micro devices (characterized by high fault electrical voltage) according to the ESD/CEI classification). Furthermore, the direct dependence of V max admissible on T allows us to confirm what is required by the industry according to which the device characterized by a deformable plate with high rigidity allows uses with an even higher V max admissible (with profiles u(x) such as not to compromise the functionality of the device). The quality of the results obtained was also confirmed by comparison with experimental setups known in the literature, highlighting that the hypothesis of small displacements, together with the non-variability of some mechanical parameters during the deformation of the plate, allows an exhaustive comparison of performing models of electrostatic devices membrane, which are famously easier to implement. Thus, the possibility emerges of opening a new line of research into comparisons between complete (but difficult to implement) models and simplified models representing devices that are well-suited to experimental performances. We underline that the model, being an ordinary derivative one, does not lend itself to being solved using FEM techniques (particularly suitable for dynamic partial derivative models). Therefore, it is necessary to reformulate the model in the near future so that its most important dynamic aspects are considered (including the dependence on electrical conductivity and temperature), on one hand, to obtain the more general conditions of global existence and uniqueness, and, on the other, to proceed with the recovering of the deformable profile using FEM techniques (performing for software/hardware applications).