Analysis of discretized parabolic problems modeling electrostatic micro-electromechanical systems

Our aim in this paper is to study discretized parabolic problems modeling electrostatic micro-electromechanical systems (MEMS). In particular, we prove, both for semi-implicit and implicit semi-discrete schemes, that, under proper assumptions, the solutions are monotonically and pointwise convergent to the minimal solution to the corresponding elliptic partial differential equation. We also study the fully discretized semi-implicit scheme in one space dimension. We finally give numerical simulations which illustrate the behavior of the solutions both in one and two space dimensions.


1.
Introduction. MEMS (micro-electromechanical systems) combine electronics and micro-size mechanical devices in order to decrease the scaling of electromechanical systems to micro-scale; this is similar to NEMS (nano-electromechanical systems) which go to nano-scale (see [27]). The idea of micro-machineries was presented by Feynman in his famous lecture (see [10]) at the end of the 1950's. Several years later, the earliest micro-machinery, named resonant gate transistor and which served as a tuner for micro-electronic radios (see [24]), was created by Nathanson and his coworkers (see [25]) at Westinghouse research labs. From then on, MEMS devices have been extensively applied to many commercial systems, including inkjet printers, MEMS microphones in portable devices, accelerometers for airbag deployment and electronic stability control in modern cars, biosensors, silicon pressure sensors, such as disposable blood pressure sensors, and so on (see more examples in [6] and [27]).
From a mathematical point of view and with the fundamental works by Pelesko and Bernstein (see [27]), we consider an idealized MEMS device which is described in the following sketch (Fig. 1). The device mainly contains a thin and deformable elastic membrane with supported boundary and a parallel rigid ground electric plate.
Elas%c membrane with conduc%ng film and supported boundary The upper surface of the membrane, which is normally dielectric, is coated with a metallic conducting film and the thickness of the film is considered to be negligible. When applying a voltage to the conducting film, the elastic membrane deforms towards the ground plate. Considering both the dynamics and electrostatic processes (see [9] and [27] for details) and applying dimensionless analysis, we obtain the following idealized parabolic MEMS problem: 2 in Ω, u(t, x) = 0 on ∂Ω; u(0, x) = 0 in Ω, (1) and the corresponding elliptic problem: 2 in Ω, u(x) = 0 on ∂Ω; 0 ≤ u < 1 in Ω, where u = 1 − d and d corresponds to the dimensionless distance between the membrane and the plate. Furthermore, f describes the dielectric profile of the elastic membrane and λ > 0 characterizes the applied voltage.
There are several central problems in the study of problem (1) and its corresponding elliptic problem (2). For instance, when the applied voltage λ increases to some threshold, the device cannot remain stable and a touchdown phenomenon appears, which means that u goes to 1 in finite time. The threshold value is called the pull-in voltage, denoted by λ * , and is defined by λ * (Ω, f ) = sup{λ > 0 | problem (2) possesses at least one solution}.
The questions of the evaluation of λ * and of how λ, as well as λ * , influence the solutions to MEMS problems are among the central questions in the mathematical study of MEMS. Moreover, we list below the definitions of quenching time and quenching set (see [15], [16] and [21]) which are also important problems in the study of MEMS parabolic problem. Definition 1.1. We call T the quenching time of problem (1) if T satisfies Definition 1.2. We call Σ the quenching set of problem (1) if Σ satisfies Besides, the diversification of materials for the membrane leads to different dielectric profiles f and corresponding solutions. We refer the readers to [9], [11], [12], [17], [19], [20], [23], [27] and references therein for more details, including the evaluation of λ * , discussions on the touchdown phenomenon and the existence and properties of the global solutions to both the elliptic and parabolic problems, as well as the study of the equations with varying dielectric properties from a theoretical perspective. Furthermore, we refer the readers to [22] for modified MEMS problems with a non-local term, to [29] for an advection term and to [8] and references therein for fourth-order problems. We are interested in this paper in the numerical approximation of problem (1). In particular, we prove that, both for semi-implicit and implicit semi-discrete schemes, the solutions are monotonically and pointwise convergent to the minimal solution to the corresponding elliptic partial differential equation, under proper assumptions. We also study the fully discretized semi-implicit scheme in one space dimension. We finally give numerical simulations which illustrate the behavior of the solutions, as well as the touchdown phenomenon, with different schemes and different initial conditions.

2.
Setting of the problem. We consider the following initial and boundary value problem: in Ω, where f describes the permittivity profile of the elastic membrane and λ > 0 characterizes the applied voltage. We make the following assumptions: • Ω is a bounded and regular domain of R N , N = 1, 2 or 3; • f ∈ C α (Ω) for some α ∈ (0, 1] and f satisfies the condition 0 ≤ f ≤ 1, but is not reduced to the null function. • u 0 ∈ L 2 (Ω) and 0 ≤ u 0 < 1 a.e..
In particular, when u 0 = 0, there exists λ * > 0 such that, if 0 ≤ λ ≤ λ * , then (5) possesses a unique solution which globally converges as t → +∞, monotonically and pointwise, to its unique minimal steady state. Furthermore, when λ > λ * , the unique solution reaches the singular value 1 in finite time. We refer the interested reader to [9] for more details.
Proposition 2. There holds, for n ∈ N ∪ {0}, where u was given in (8), and c 0 > 0 is the optimal constant in the Poincaré inequality Proof. We have, Multiplying (15) by v n+1 and integrating over Ω and by parts, we obtain which yields, noting that v n and v n+1 are nonnegative and recalling that u n ≤ u λ and the assumptions on f , It thus follows from (12) and the Poincaré inequality that This completes the proof.
As a consequence of Proposition 1, we deduce the following corollary.
Corollary 1. We assume that Then, u n converges to u λ in L 2 (Ω) as n → +∞.
Remark 1. We can note that (16) holds for In particular, this estimate does not depend on the choice of τ > 0.

Remark 2.
Let us take u 0 = 0. It is clear that u 1 ≥ u 0 = 0. Let us then assume that, for a given n ∈ N, u n ≥ u n−1 . We have, and v n+1 = 0 on ∂Ω.
Therefore, for every x ∈ Ω, the sequence {u n (x)} is monotone increasing and bounded from above by u λ . Hence it converges monotonically as n → +∞.
Moreover, it is verified numerically (see [9] for details) that, if λ < λ * and is close to λ * , with given f and N , there may exist another stable solution to the elliptic problem (7) which is denoted by u + λ and satisfies u λ < u + λ < 1. We now assume that the stationary problem possesses at least 2 solutions u λ < u + λ < 1 and that the initial condition satisfies u λ ≤ u 0 < u + λ < 1. Then, we additionally have the following proposition.
Proposition 3. We consider the semi-implicit scheme (6), with 0 < λ < λ * and f satisfying the assumptions mentioned in Section 2. Then, if the initial condition Proof. We proceed as above. On the one hand, for n = 0, then u λ ≤ u 0 < 1. We further assume that, for a given n ∈ N ∪ {0}, there holds u λ ≤ u n < 1, a.e. x ∈ Ω.
where u λ is the minimal solution to (7).
Proof. It is obvious that 0 is a subsolution to problem (27). Furthermore, we note that so that u λ is a supersolution to (27). The proof is completed.

Proceeding as in the proof of Proposition 2 and setting
Multiplying (29) by v n+1 and applying the Poincaré inequality, we obtain This yields the following result.
Remark 4. We again take u 0 = 0. Then, u 1 ≥ u 0 . Let us assume that, for a given n ∈ N, u n ≥ u n−1 . Then, the function and v n+1 = 0 on ∂Ω.
Multiplying (33) by −v − n+1 and integrating over Ω and by parts, we find Therefore, for a given λ, if τ ≤ τ 0 , where τ 0 is small enough so that then v − n+1 = 0 and thus u n+1 ≥ u n . It thus follows from Proposition 4 that u n converges monotonically and pointwise.

5.
The fully discretized semi-implicit scheme. We only consider the onedimensional problem in this section. We believe that similar results hold in twodimensional space and will address this elsewhere.
Let M > 0 be an integer, h = (b − a)/(M + 1) denote the spatial mesh size, and x i = a + ih, i = 0, · · · , M + 1. We consider the fully discretized semi-implicit scheme as follows: for n ∈ N ∪ {0}, where u n i ≈ u(t n , x i ). Problem (39) can be rewritten equivalently as We then rewrite (40) in vector form as where U n = (u n 1 , u n 2 , · · · , u n M ) t and We can rewrite (41) equivalently as where We note that It is well known that A is positive definite, and thus invertible. Furthermore,

so that
This yields that A −1 ≥ 0.
We now consider the equation which is the centered difference scheme of the steady state problem (7).

Equation (44) can be rewritten as
The solvability of the above problem is given in the following theorem.
Theorem 5.1. We assume that 0 ≤ λ ≤ 8(1 − δ)δ 2 (b − a) 2 , for some δ ∈ (0, 1). Then (45) possesses a solution U * such that 0 ≤ U * < E, where E = (1, 1, · · · , 1) t 1×M . Proof. First note that B is invertible and (as above) B −1 ≥ 0. We rewrite (45) in the form We then consider the sequence Note that W 1 exists and, since B −1 ≥ 0, W 1 ≥ 0. Actually, as long as it exists, as long as this makes sense. Now, recalling the assumptions on f , we have We denote by D M the determinant of B. It is easy to see that Therefore, since D 2 = 3 and D 3 = 4, it follows that D M = M + 1. It thus follows from [3] (see also [14]) that B −1 is the factorizable matrix

Then, we have
This yields that {W k } exists for all k ∈ N ∪ {0} and Hence, each component of {W k } is bounded and monotone increasing, and thus converges. This finishes the proof.
We assume from now on that the assumptions of Theorem 5.1 hold, so that 0 ≤ U 1 ≤ U * . Let us assume that 0 ≤ U n ≤ U * . Then, setting V n = U n − U * , there holds Furtheremore, since A −1 ≥ 0, U n+1 ≥ 0. It thus follows that {U n } is well-defined and 0 ≤ U n ≤ U * , n ∈ N ∪ {0}. We can also prove, proceeding as above, that where ((·, ·)) denotes the usual Euclidean scalar product, with associated norm · . It follows from the discrete Poincaré inequality that We now have AV n+1 = V n + G(U n ) − G(U * ), so that, taking the scalar product with V n+1 , Therefore, Therefore, if Remark 8. The above results, as well as those in the previous sections, can be improved by taking an assumption on f of the form In particular, in Theorem 5.1, we then have the more general assumption which allows for a larger λ whenf is small.
6. Numerical simulations. In this section, we give several numerical simulations which show the behavior of the solutions u corresponding to different schemes, different λ's and different initial conditions (see Proposition 3 and Remark 3). In particular, it is verified that, when 0 < λ ≤ λ * , the solutions to problem (5) tend to a stable solution u λ as time grows. Furthermore, if λ > λ * , one can observe the so-called touchdown phenomenon. The numerical simulations are performed with MATLAB.
The branches and solutions for the one-and two-dimensional elliptic problems are helpful in view of the study of the parabolic problems.
6.2. The parabolic problem. We set t n = nτ , n = 0, · · · , K, and write the time steps as superscripts and the spatial nodes as subscripts.
which is larger than u λ and less than but close to u + λ show that the solution decreases and converges to u λ (see Fig. 6(a)). In Fig. 6(b), the solution for the nonsymmetric initial condition also converges to u λ . However, for λ = 4.45 > λ * and a smaller time step τ = 0.001, we observe in Fig. 6(c) that, after 2654 steps, there does not exist a stable solution to the parabolic problem. We have similar results when applying the implicit scheme; these are not displayed here.
is called the tangent vector induced by A, where (·) * denotes the Hermitian transpose.
Here and below, M = M 2 . We then use an approximate Euler predictor and Newton-type iterations as corrector steps, described in the following algorithm, see [1].

Algorithm 1 Continuation Method
Input N pas, and w such that H(w) = 0 for n = 1, N pas do In this algorithm, A + denotes the Moore-Penrose inverse of A and stands for the current step size. We refer the readers to [5] and [26] for Newton's method based on the Moore-Penrose inverse and [2] for a convergence result which ensures that the above algorithm is applicable and effective under reasonable assumptions. To compute A + and t(A), we use a QR factorization (see [1] for details), that is to say, A being a maximal rank M × (M + 1) matrix and A * representing its conjugate transpose, we have with the QR factorization, A * = Q R 0 * , where Q is an (M +1)×(M +1) orthogonal matrix and R is a nonsingular M × M upper triangular matrix. Then, we compute t(A) = σq, where q denotes the last column of Q and, in order to satisfy the orientation which has been defined in [1](2.5), we can choose σ = sign(det Q det R). The main idea of the continuation method is displayed in Fig. 10, using a rough prediction and several steps of corrections.