Recovering of the Membrane Proﬁle of an Electrostatic Circular MEMS by a Three-Stage Lobatto Procedure: A Convergence Analysis in the Absence of Ghost Solutions

: In this paper, a stable numerical approach for recovering the membrane proﬁle of a 2D Micro-Electric-Mechanical-Systems (MEMS) is presented. Starting from a well-known 2D nonlinear second-order differential model for electrostatic circular membrane MEMS, where the amplitude of the electrostatic ﬁeld is considered proportional to the mean curvature of the membrane, a collocation procedure, based on the three-stage Lobatto formula, is derived. The convergence is studied, thus obtaining the parameters operative ranges determining the areas of applicability of the device under analysis.


Introduction and Problem Statement
At present, micro-components are essential in the embedded engineering applications. In particular, micro-transducers' and micro-actuators' results are of paramount importance due to their role of micro-devices' interfaces [1,2]. In recent years, circular membrane MEMS technology has been exploited to this aim in a wide variety of technological fields such as thermo-elasticity [2][3][4], microfluidics [2,[5][6][7], electroelasticity [8][9][10][11][12], and biomedical engineering [13][14][15]. The focus on this geometry can mainly be attributed to the availability of theoretical models in both steady-state and dynamical conditions (see [16][17][18] and references within). However, it is worth noting that these models, even though they are very detailed, are often characterized by implicit structures so that the explicit solutions are not often obtainable. As a consequence, analytical and/or algebraic conditions ensuring both existence, uniqueness, and regularity of the solution should be achieved [19], in such a way that computational results provided by numerical procedures can be checked for validating their correctness and thus avoiding possible ghost solutions. We remember that the ghost solutions are numerical solutions in convergence conditions that do not respect both analytical existence and uniqueness conditions [20,21]. On these bases, in this study, we have focused our attention on a simplified 2D circular membrane MEMS actuator in which the membrane, anchored to the edges of a fixed plate, deforms towards another fixed plate because of an applied external electrical voltage [22,23]. During the deformation of the membrane, the electrostatic field E generated inside the device can be considered locally orthogonal to the membrane tangent line, so that the modulus |E| depends (locally) Remark 1. It is worth nothing that u 0 is a very important quantity in electrostatic. In fact, the electrostatic capacitance (C el ), co-energy function (W ), electrostatic charge on the membrane (q), electrostatic force ( f el ), and |E(r)| are dependent on u 0 . In particular, we can write (for details, see [2,22,25]): 2 (1) 2 (2) 2 (3) Figure 1. Representation of an electrostatic circular membrane MEMS device: the membrane is deformed towards the upper disk while the lower disk is at potential V = 0.

Electrostatic Membrane MEMS Actuator: Some Mathematical Backgrounds
The starting point of the model studied in [22] is the well-known steady-state model largely studied in [19]. Specifically, a MEMS device consisting of two parallel metallic plates, where the lower one (at u(x) = 0, ∀x ∈ Ω = [0.5, 0.5]) is fixed and the other one (at u(x) = 1, ∀x ∈ Ω = [0.5, 0.5]), subjected to a voltage V, is deformable but anchored at their edges, is considered. The model is: in which the existence of at least a solution has been studied exploiting Steklov boundary condition obtaining Dirichlet and Navier boundary ones. In (4), f 1 represents a bounded function that depends on the dielectric properties of the material constituting the deformable plate, and λ 1 is a function depending on V. Moreover, u ν denotes the outer normal derivative of u on ∂Ω. It is worth nothing that, in (4)d = 0, one obtains the Navier boundary conditions; ifd = ∞, one obtains the Dirichlet boundary conditions. Model (4) has been exploited in several works [3,26] in the following dimensionless version: where L is the semi-length of the dimensionless device in order to propose a version of the model where λ 2 is a parameter related to the external voltage V as follows: in which 0 represents the permittivity of the free space, T is the mechanical tension of the membrane, and d is ten distance between the two parallel plates. It is worth nothing that in (6) represents the parameter that takes into account the electro-mechanical properties of the material constituting the membrane. In (5), the quantity λ 2 (1−u(x)) 2 is physically proportional to |E| 2 , so that it makes sense to write: with θ being a parameter of proportionality. Then, in [19], since E is locally orthogonal to the tangent straight line of the membrane, |E| has been considered proportional to the curvature K of the membrane itself. Thus, taking into account that , and model (4), in [19] assumed the following form: obtaining [19]: in which the singularity 1 (1−u(x)) 2 present in (5) is not here explicitly evident. In this model, the membrane has replaced the deformable plate, and |E| has been considered as proportional to the curvature of the membrane itself. The achieved results highlighted that the analytical uniqueness condition depended on the electro-mechanical properties of the material constituting the membrane, while the numerical solutions, obtained by a shooting procedure, highlighted the range of θλ 2 ensuring convergence of the procedure (without ghost solutions). Moreover, the shooting approach has been compared with the well-known Keller-box scheme getting an optimal range of θλ 2 (in the absence of ghost solutions).

A New 2D Model for Electrostatic Circular Membrane MEMS: |E| in Terms of Mean Curvature
In [22], the attention was focused on a 2D circular membrane MEMS actuator, exploited in several biomedical and industrial applications, characterized by an axial symmetry in the geometry of the membrane. In particular, considering the z-axis as a rotation one, the profile u of the membrane has been thought as a surface generated by a vertical rotation of a curve C around z located in the first quadrant, when R ≥ r ≥ 0. In such conditions, u was only dependent on r transforming the 2D problem under study into a 1D one replacing x by r. In such a way, ∆ 2 u(x) in (4) reduced into its radial part, so that (4) became: in which the singularity 1/r appeared. As shown in [22], λ 2 /(1 − u(r)) 2 ∝ |E| 2 so that model (13) in [22] has been rewritten as follows: In particular, θ has been supposed, for physical reasons, to be a continuous function depending on r ∈ Ω = [−R, R]. Moreover, K(r, u(r)) represented the curvature of the deformed membrane, so that |E| = µ(r, u(r))K(r, u(r)), with represented the proportionality function. Finally, model (14) became (see [22]): where K(r, u(r)), being in 2D formulation, becomes the mean curvature. As shown in [22], the mean curvature of the membrane, H(r), exploiting some procedures based on differential geometry concepts, can be easily written as so that model (17) became: (20) in which d * represents the critical security distance, ensuring that the deflection of the membrane does not touch the fixed upper plate. In [22], the following result of existence of at least one solution for (20) has been presented and proved as well.

Remark 2.
It is worth nothing that the mathematical models that adhere to the physical reality of MEMS are extremely complex and do not allow analytical studies on them. Therefore, it is necessary to implement some simplifications in the geometry of the MEMS device so that the simplified model obtained can be studied mathematically. It follows that the qualitative assessments obtained from the mathematical study of the model (20) (which is a simplified model) will hardly agree with any experimental data but will give (qualitative) indications of the behavior of the MEMS device characterized by simplified geometry.

Theorem 1.
Let us consider the model (20). Moreover, let 1 (r) and 2 (r) be two functions, both defined in [0, R] and twice continuously differentiable, in order that 1 (r) < 2 (r). In addition: be a continuous function obviously, (except for r = 0) satisfying the Generalized Lipschitz condition where 0 is the permittivity of the free space and k is the constant of proportionality between the electrostatic pressure p el and the displacement at the center of the membrane u 0 ( u 0 = kp el ), there exists at least one solution for the problem (20).

Remark 3.
It is worth remembering that, as explained in [22], the greater k is, the lower the value of θλ 2 will be, so that the concavity of the membrane will rise (see (24)). However, even if (20) admits at least one solution, as specified in Theorem 1, its uniqueness has not been ensured (see Theorem 2 in [22]).

On the Applicability of the Numerical Procedure
As well known, the BVPs are much harder to solve than IVPs [27]. Unlike IVPs that have a unique solution, a BVP may not have a solution, or may have a finite number, or may have infinitely many. Moreover, BVPs defined on an infinite or semi-infinite interval are usual. The classical numerical approaches for solving BVPs are the shooting method and the collocation one. The first one combines a numerical method based on the solution of IVPs for ordinary differential equations (ODEs) and one for the solution of nonlinear algebraic equations. The main difficulty with the shooting method is that a perfectly well-behaved BVP can require the integration of IVPs that are not stable. This means that the solution of a BVP can be insensible to the changes in the boundary values, yet the solutions of the IVPs of shooting are sensible to the changes in the initial values. On the other side, the collocation methods are efficient and reliable tools, though, often, the underlying method may not be appropriate for high accuracies and for problems with hardly sharp changes in their solutions. The authors have used the collocation with piecewise polynomial functions and found it an appropriate method for solving the singular BVP under study because only a coefficient is singular and the solutions are smooth [27]. Then, they propose the collocation method implemented in the bvp4c routine in Matlab. The bvp4c solver is a code that implements the three-stage Lobatto IIIa formula. This is a collocation formula, and the collocation polynomial provides a C 1 continuous solution that is fourth-order accurate uniformly in the whole computational domain. Finite difference methods for BVPs are proposed in [28].

The Exploited Numerical Approach: A Scheme with the Three-Stage Lobatto IIIa Formula
This procedure solves a system of ODEs in the following form [24]: where g[u(a), u(b)] = 0 represents the boundary conditions.

Remark 4.
Generally, (25) can be written as du dr = F(r, u(r), p) in which p represents a vector of unknown parameters. However, for our purposes, the system can be written as (25) (for more details, see [24]).
In this work, the authors have exploited the MatLab R R2017a bvp4c solver (Natlick, MA, USA), running on an Intel Core 2 CPU 1.45 GHz machine (Santa Clara, CA, USA) because: it implements a collocation technique exploiting a piecewise cubic polynomial, S(r), whose coefficients are determined by the requirement that S(r) be continuous on (a, b); 2.
the mesh and the estimation of the error are both based on the evaluation of the residual of S(r); 3.
the control of the residual is particularly useful to handle poor guesses for the mesh and the solution as well; 4.
the computational complexity to obtain the Jacobian is reduced; 5.
being a vectorized solver, it reduces the run-time vectorizing F(r, u(r)).
Let us start with the following two definitions.

Definition 1 (step-size definition).
Let h m be the step-size computed as r m+1 − r m in the mesh-grid 0 = a = r 0 < r 1 < ... < r n = b = R.

Remark 5. The cubic polynomial S(r) satisfies the boundary conditions in (25) that is g[S(a), S(b)] = 0.
Moreover, ∀(r m , r m+1 ) of the mesh, the subdivision 0 = a = r 0 < r 1 < ... < r n = b = R is considered, and S(r) collocates at the ends of the sub-interval and the midpoint. Obviously, S(r) is continuous at the endpoints of each sub-interval.
This collocation approach is totally equivalent to the three-stage Lobatto IIIa implicit Runge-Kutta procedure [24].
Then, after pressing the previous definitions, we are able to show the fourth-order Lobatto IIIa formula as follows: Remark 6. Obviously, when the procedure is applied to a quadrature problem, it reduces to the well-known Simpson formula:

Remark 7.
It is worth noting that the global discretization error of the procedure only contains even power of the step-size. This occurrence makes the method particularly attractive for both extrapolation and deferred correction.

Remark 8.
We observe that S(r), together its derivatives, satisfies the following condition [24]: In addition, at each intermediate point, S(r) satisfies the Equations (25) at the midpoint of each interval justifying the term "collocation polynomial". Moreover, MatLab R chooses the form of S(r) by determining any unknown parameters. Finally: Conditions (30)-(32) are nonlinear equations iteratively solved by a linearization procedure; then, the linear equation solver of MatLab R is used. Moreover, if necessary, the cubic polynomial can be evaluated ∀r ∈ (a, b) by the bvpval function.

Remark 9.
It is worth noting that a BVP could have more than one solution so that the users should supply a guess for both the solution and the initial mesh. Then, the solver is able to adapt the mesh to achieve a solution with a reduced number of mesh points.
Obviously, a good guess is often tricky so that the solvers control the residual defined as follows.
Definition 3 (the residual). The residual in a ODE is defined as follows:

In addition, in the boundary conditions, it can be written as g[S(a), S(b)].
Remark 10. We observe that, if res(r) is small, then S(r) can be considered as a good solution. In addition, if the problem is well-conditioned, small res(r) implies that S(r) is next to u(r).

A More Suitable Writing of the Analytical Model
The two-point BVP under study (20) is singular with respect to the independent variable at the initial point r = 0, due to the division by zero of the first term on the right-hand side of (20). We must deal with the singularity in the coefficient at r = 0 because the bvp4c solver always evaluates the ODEs at the endpoints. For this problem, we show how providing the correct value at the singular point is enough for computing the numerical solution. By the theoretical results presented in the previous sections, a smooth solution is expected. Indeed, symmetry implies that this smooth solution has a derivative that vanishes at the origin r = 0. This physical condition allows for demonstrating that so that the singular term in the ODE is well defined at the origin for a well-behaved solution. Then, if we let r → 0 in the equation, we find that and, solving for dr 2 , we obtain the value that must be used for evaluating the ODEs at the point To solve the problem under study (20) by means of the collocation method, we rewrite it as a system of first order ODEs of the form by setting u 1 (r) = u(r) and u 2 (r) = du 1 (r) dr = du(r) dr with the supplementary condition given by Equation (36) at initial point r = 0. As specified above, when an external V is applied, the membrane inside the device moves towards the upper disk. The higher is the applied V, the closer the membrane to the upper disk will be. In addition, instability phenomena of the membrane can arise when V grows too much; so that it is important to know the range of V generating instability. However, V is linked to θλ 2 by inequality (24). Then, it follows that knowing the behavior of the membrane when θλ 2 increases gives us, once both d * and k are fixed, the range of V producing instability of the membrane in absence/presence of ghost solutions.

Remark 11.
It is interesting to observe that θλ 2 , being dependent on both the electromechanical properties of the material constituting the membrane and V, once the ranges of stability/instability of the membrane are known, it is possible to know the operation parameters in convergence area respecting (24) and the engineering areas of applicability of the device.
Here, it is reasonable to set k = 1. In other words, we consider that p el and p are totally equivalent. This makes sense because, when V is externally applied, |E| and p el are generated inside the device. If we suppose that the losses are negligible, then p el = p (that is, k = 1). Moreover, we set d * = 10 −9 because lower values of d * do not make sense from a physical point of view. The numerical procedure has been applied for different values of θλ 2 and u 1 = 0.1 and u 2 = 0 as initial guesses. All the simulations have been carried out by bvp4c MatLab R solver on an Intel Core 2 CPU 1.45 GHz machine, using both the default relative and absolute error tolerances. Moreover, the optimal number of grid points, on the computational domain [0, R] = [0, 1], has results of 100. In fact, by a greater number of grid points, the performance of the numerical procedure, qualitatively did not improve.
Moreover, by numerical experiments, we obtain different solutions starting by different initial guesses. The following three cases hold: 1. ∀θλ 2 ∈ (10 −6 , +∞) the numerical procedure converges, and no problem arises, and instabilities do not occur. In other words, when θλ 2 increases starting from 10 −6 , (20) increases more and more (from negative values towards zero), so that the concavity of the deformed membrane more and more decreases avoiding instabilities next to the edge of the membrane. This is confirmed by the condition (24), which ensures the existence of the solution for (20). In fact, from (24), the higher θλ 2 is, the lower V 2 will be; that is, increasing θλ 2 , the membrane does not deform too much as if a low external V is applied. Figure 2 displays an example of recovering of the membrane when θλ 2 = 0.5, with initial guesses u 1 ≤ 2.446 and u 2 = 0. Being V reduced, the membrane moves towards the upper disk just a little so that instability phenomena do not appear. It is worth noting the symmetry of the first derivative with respect to (0, 0). The behavior of the procedure is different when the initial guess of u 1 increases (with u 2 = 0). In this case, the solution is very different with respect to the previous cases. In fact, Figures 3-5 show examples of recovering of the membrane with θλ 2 = 0.5 and with an initial guess for u 1 belonging to [2.447, 2.453], [2.454, 9.474], [9.63, 12.7], and [15.1, 19.978], [9.475, 9.62], [12.71, 15] and [19.979, ∞), respectively (initial guess for u 2 is zero). In particular, increasing the initial guess for u 1 and u 2 = 0, the recovering of the membrane results in being erratic but still symmetrical (see Figure 3) until the profile assumes a bell-shape (see Figure 4). The erratic behavior of the solution appears again when the initial guess increases (see Figure 5). However, it is worth noting that, even if in Figures 3-5 show simulations that make sense from a numerical point of view, they are not realistic because the achieved profile u 1 is greater than d.

2.
∀θλ 2 ∈ (0, 10 −7 ), the numerical procedure does not work. It means that ( (20) increases too much, (theoretically , so that the numerical procedure is forced to stop because the Jacobian matrix becomes singular.

3.
∀θλ 2 ∈ [10 −7 , 10 −6 ] strong instabilities take place. In this case, even if the numerical procedure does not stop, instabilities next to the edge of the recovered membrane occurred. Figures 6 and  7 display two examples of recovering of the membrane when θλ 2 ∈ [10 −7 , 10 −6 ] highlighting strong instabilities of the membrane. In particular, Figure 6 has been achieved setting u 1 = 0.1 and θλ 2 = 5 · 10 −7 , while Figure 7 has been obtained setting u 1 = 1.2 and θλ 2 = 10 −6 . Obviously, depending on the considered initial data, the values of θλ 2 change for which one has no instabilities.

Remark 12.
In the resolution of BVPs, an important aspect concerns the choice of initial conditions. The difficulty therefore lies in choosing suitable initial conditions that lead to a convergence of the method. Usually, if you know the progress of the solution, they are chosen in such a way as to be as accurate as possible, especially to achieve convergence and calculation time that is not too long. It starts by assigning values and proceeding with testing until the solution we expect is found. The results shown in Figure 3 highlight a non-convergence of the method; this is due to not "good" initial values. Furthermore, the different acceptable solution (Figures 2 and 4) obtained with different initial values show how the BVPs cannot admit a single solution. (We only consider Figure 2 because the solution of the model must be lower than 1).

Order of Magnitude of θλ 2 According to the Analytical Condition of Uniqueness of the Solution
From inequality (24), setting d * = 10 −9 , k = 1 and 0 ≈ 10 −12 , we can write: from which:
If the numerical procedure converges with instability phenomena, θλ 2 ∈ [10 −7 , 10 −6 ]. In this case, we can write: for which it follows: obtaining that, ∀V 2 k ∈ (10 −12 , 10 −11 ), the numerical procedure converges even if instability phenomena can occur next to the edges of the membrane.
Then, it follows that the analytical condition (3) can be written as θλ 2 > 10 −12 . Finally, ∀θλ 2 ∈ (0, 10 −12 ], any solutions can be considered as ghost solutions. However, because ∀θλ 2 ∈ (0, 10 −7 ], the numerical procedure does not converge; then, it follows that each numerical solution achieved is not a ghost solution. The results discussed in this Section are summarized in Table 1.
No Ghost Solutions No Ghost Solutions

Electromechanical Properties of the Membrane, Applied Voltage, and Ghost Solutions: Engineering Areas of Applicability of the Device
As discussed above, θλ 2 is very important for convergence of the numerical procedure and for its physical meaning as well. In particular, θ represents the proportionality between |E| and λ 2 (1−u(r)) 2 [20]: Moreover, λ 2 depends on both V and the electromechanical properties of the material constituting the membrane. In fact, λ 2 can be expressed as follows [20]: where ρ = 0 (2R) 2 2d 3 T (ρ takes into account the electromechanical properties of the material constituting the membrane) and T is the mechanical tension of the membrane. Then, combining equality (49) with equality (50), we obtain: In addition, multiplying (51) by λ 2 and taking into account the (50), we achieve: obtaining Considering that, in dimensionless conditions, (1 − u(r)) 2 < 1, d = R = 1 and |E| 2 < sup{|E| 2 }, from (53), the following inequality makes sense: As known, θλ 2 can be expressed as (for details, see [19,20]): from which, it follows that: If the non convergence of the numerical procedure occurs, it makes sense to write the following chain of inequalities: so that we can write: Then, from (58), we can obtain: (60) Once the intended use of the device has been chosen (that is, once the pair {V, sup{|E| 2 } has been fixed), from the inequality (59), inf{T} is computable. In other words, once {V, sup{|E| 2 } has been chosen, the material of the membrane is selected. Conversely, if the material of the membrane has been chosen (that is T has been fixed), it is possible to achieve the pair {V, sup{|E| 2 } satisfying the inequality (60) (that is the intended use of the device).

Conclusions and Perspectives
In this paper, a numerical approach based on a three-stage Lobatto technique for recovering the profile of the membrane in an electrostatic circular MEMS is proposed. Mainly exploited as a numerical collocation technique for the resolution of BVPs, this procedure has been preferred to shooting methods because it usually requires the integration of IVPs being heavily unstable. The proposed numerical procedure has been applied to a well-known 2D nonlinear second-order differential model related to a circular membrane MEMS actuator in which the amplitude of the electrostatic field has been locally expressed in terms of the mean curvature of the membrane. The study of the convergence properties of the proposed method has allowed for highlighting the ranges of characteristic parameters (depending on both V and electromechanical properties of the material constituting the membrane) that ensure that the method can achieve convergence with or without instabilities. Furthermore, exploiting an algebraic inequality governing the existence of the solution for the differential model, the ranges of parameters (depending on the externally applied voltage) that ensure convergence (with or without instabilities) bringing out the ghost solutions' areas have been evaluated. Finally, an algebraic inequality has been obtained, which allows for managing the choice of the intended use of the device once the electromechanical properties of the material constituting the membrane have been selected (and vice versa). The results we have obtained can be considered as encouraging. The tests have shown that the recovering of the membrane profile (in the absence of ghost solution, and a condition of convergence without instability) reach a maximum value at the center of the membrane of not less than 0.3, thus proving that the membrane offers an excellent response to the applied external voltage. Furthermore, the tests carried out have shown that the non-convergence area is extremely limited, thus guaranteeing a wide range of applicability of the device in industrial applications. However, we point out that the presence of instabilities makes the numerical approach not yet optimal, requiring further effort to reduce the instability zone as much as possible. This will be the subject of a future work.