Numerical Approaches for Recovering the Deformable Membrane Profile of Electrostatic Microdevices for Biomedical Applications

As known, microelectromechanical systems (MEMS) allow for the creation of both sensors and actuators capable of using non-electrical devices in microchips [1,2]. Fluidic devices such as micropumps were historically the first to be made and, with the discovery of new pharmacological therapies [3,4,5,6], today they represent a highly developed research area. Currently, microfluidics is based on sophisticated physical-mathematical models which, unlike macroscopic applications, are able to take into account electrokinetic [7], magnetohydrodynamic [8], electrochemical [9], electrostatic [1] and other effects [10,11]. Concerning the design of a micropump, the actuator is one of the most important devices since it converts the applied external energy (electrical, thermal, …) into movement of the deformable element (usually, a membrane). The mechanical micropumps (piezoelectric [12], electrostatic [13], thermopneumatic [14], electromagnetic [15], bimetallic [16], ion conductive [17], phase change [18], shape memory alloy polymer films [19]) require an actuator for pumping; usually, they consist of a pumping chamber separated by a deformable diaphragm so the fluid flows by means of oscillations which generate pressures depending on the volume of the stroke inside the chamber produced by the actuator. Of course, performance is severely limited by the mechanical components. Non-mechanical micropumps, having no moving parts, do not require mechanisms for converting the energy administered into kinetic momentum. However, while easy to design and build, they can pump low conductivity fluids [20,21,22,23,24,25,26]. Research efforts are currently focused on the efficiency of drug administration to be combined with the efficiency of the fluidic device and of the micropumps (implantable and/or removable) based on MEMS technology since positive results in various clinical conditions depend on them [27,28]. Since the release and absorption of drugs strongly depends on pharmacokinetic and metabolic mechanisms, in order to reduce adverse phenomena, research is oriented towards the development of devices capable of controlling the dosage rate from the delivery system [29,30]. MEMS-based micropumps are designed to deliver drugs at a specific rate in accordance with the required pharmacokinetics also allowing prolonged infusion volume control (useful for certain pathologies) and improving the patient’s quality of life [31,32,33,34]. Electrostatic membrane MEMS micropumps currently have strategic importance as they offer a high level of performance with almost a total absence of maintenance interventions [28,35,36]. Furthermore, their easy design as well as their low cost in the manufacturing step allow for large-scale industrial production [37,38].

Although the production of electrostatic membrane MEMS micropumps has reached considerable levels, there is still a need to study new and accurate physical-mathematical models to test, in the design phase, the performance of the micropumps stressed by the different types of loads to which they will be subjected during their use. Recently, new models for electrostatic membrane micropumps based on integro-differential equations have been presented, whose solutions, also depending on the charge densities, can be obtained by analytical procedures [39,40]. Even if such models realistically simulate the multi-physics present in the devices, the high computational loads they entail do not facilitate hardware prototyping. Therefore, the development of models (especially those with a high degree of symmetry) with a reduced computational load will allow hardware to operate at a lower cost [41,42,43,44,45,46,47].

Concerning circular steady-state electrostatic membrane MEMS micropumps, recently, a dimensionless second-order semilinear elliptic model with singularity (whose solution was the profile of the deformable membrane, $u\left(r\right)$, $r\in [0,R]$, with r being the radial coordinate and R being the radius of the micropump membrane) was studied in [48,49,50], without considering effects due to the fringing field, where important results of existence [48] and uniqueness were obtained [49] without, however, guaranteeing any stability [50]. In particular, where $d\in {\mathbb{R}}^{+}$ is the distance between the membrane at rest and the counter-electrode. $\lambda$, as defined in (4), is a positive dimensionless parameter that depends on both the applied voltage, V, and the mechanical tension, T, of the membrane at rest. $d^{*}$ is the critical security distance that ensures that the membrane does not touch the counter-electrode. $\theta \in {\mathbb{R}}^{+}$ (physical parameter without limitation) considers the electrical properties of the membrane [51,52], so that $\theta \lambda$, which is subjected to any limitation if numerical recovery of the membrane is performed, expresses the electro-mechanical properties of the membrane.

Model (1) is interesting because it models a device that is controllable in terms of voltage because, as shown in (4), $\lambda$ depends on V (so that $\lambda$, like V, is a bounded parameter). Moreover, it does not consider the effects due to fringing [53,54], and in parallel, it is not particularly adherent to the main requirements of micropumps for drug delivery systems because the dielectric properties of the membrane are modelled by a function, $f\left(r\right)=1$, $\forall r\in [0,R]$, as if the membrane, dielectrically exhibits the same behavior at each of its points. In the recent past, important results have been obtained concerning analogous models of MEMS devices, even in the presence of fringing field, without defining a suitable dielectric profile of the membrane [51], while the dielectric profile was formulated according to some experimental evidences, where, however, the deformable element was a rectangular metal plate [52]. Finally, any instability of the membrane could cause electrostatic discharge. It is worth noting that, when the membrane deforms, strong variation in the electrostatic capacitance occurs in the device, for which an auxiliary electrostatic capacitance, $C_f$, is needed to oppose the variation of V.

With these premises, the main purpose of this work can be summarized by the following points:

The remainder of the work is organized as follows: Section 2 details the studied micropump, framing the problem from the points of view of both the actuator and the transducer. Then, Section 3 describes the proposed physical-mathematical model, presenting important results related to the existence, uniqueness, and stability of the solutions (in Section 4, Section 5, Section 6 and Section 7). After the dielectric properties of the membrane have been modeled, an important limitation for the electrostatic force is presented and discussed (Section 8). Next, Section 9 is dedicated to the numerical recovery of the membrane profile for both different loading conditions and different effects due to the fringing field. Then, after quantifying the V needed to overcome the mechanical inertia of the membrane (Section 10), a link between V and T is proposed in Section 11. Finally, before concluding the paper, possible future developments are provided (Section 13), and some considerations regarding both the pull-in voltage and electrostatic pressure are discussed in Section 12.

In (1), $\lambda /{(1-u\left(r\right)-d^{*})}^2$ represents a term proportional to $\left|\mathbf{E}\right|$ [50,51,53]. As anticipated above, we introduce a positive bounded function, $f\left(r\right)$, to define the dielectric properties of the membrane (usually referred to in the literature as the dielectric profile of the membrane). Furthermore, to take into account the effects due to the fringing field, we consider an additional term that is equal, according to the Pelesko–Driscoll approach, to ${\lambda f\left(r\right)\delta |\nabla u\left(r\right)|}^2$ [51], where $\delta \in {\mathbb{R}}^{+}$ weighs the effects due to the fringing field. Then, (1) becomes Furthermore, by exploiting the radial symmetry with respect to the vertical axis and considering that $\nabla u\left(r\right)=\frac{du\left(r\right)}{dr}$, (17) can be written as Moreover, by exploiting the idea that ${\left|\mathbf{E}\right|}^2\propto \lambda f\left(r\right){(1-u\left(r\right)-d^{*})}^{-2}$, it makes sense to write In the past, it was proven that $\mathbf{E}$ is locally orthogonally to the straight-line tangent to the membrane [50,51,53]. Therefore, considering $K(r,u(r\left)\right)$ as the curvature of the membrane (see (3)), we can write $\left|\mathbf{E}\right|\propto K(r,u(r\left)\right)$, whose function of proportionality $\mu (r,u\left(r\right),\lambda ,f\left(r\right))=\lambda f\left(r\right){(1-u\left(r\right)-d^{*})}^{-1}$ belongs to $C^0([-0.5,0.5]\times [0,1)\times [0,1]\times [0,1])$ which in [51] has been shown to have no limitation. Therefore, (18), after some calculations, becomes where represents the electromechanical load on the membrane when V is applied (V generates $\mathbf{E}$ between the electrodes which determines $f_{el}$ which, for each surface unit, is transformed into $p_{el}$. Obviously, model (20) is difficult to solve analytically; therefore, in the first instance, we verified the existence, uniqueness, and stability of the solution in this new formulation where $f\left(r\right)$ takes place. It is worth noting that, in (20), local mechanical stresses were not considered in order not to make the model depend on tensor variables which would further reduce $u_0$. However, not considering this contribution allows for a possible overestimated numerical recovery “with a safety advantage”, since there is the certainty that the real value of $u_0$ is lower than the numerically obtained one. Thus, the membrane will certainly not touch the counter-electrode of the micropump. We also observed that $f\left(r\right)$, which influences the electrostatic charge in the micropump, is decisive for the functioning of the micropump, especially when the deformation of the membrane is the maximum allowed, as already experimented in similar devices [56,57]. This is due to the fact that $u\left(r\right)<1-d^{*}$ because, physically, the membrane must not touch the top plate (mathematically imposed in (20)). Then, if $u\left(r\right)=1-d^{*}$, from (20), it follows that $f\left(r\right)=0$, so (21) will vanish.

Furthermore, (20) models circular electrostatic micropumps, even for small displacements of the membrane. Indeed, if $\pi R^2\gg d^2$, $\delta$ will be negligible. However, from $\pi R^2\gg d^2>u^2\left(r\right)$, it follows that $u\left(r\right)\ll \sqrt{\pi R^2}\approx {10}^{-6}$ m, which provides the possibility of performing small displacements of the membrane. Obviously, (20) is formulated with an extremely simplified micropump geometry, allowing a preliminary analytical study to be carried out, of which the results, still unconformed by experimentation, provide an important theoretical contribution.

In this section, we ask ourselves whether the presence of $f\left(r\right)$ in (20) substantially modifies the algebraic condition of existence already known in the literature. In this regard, the following result holds.

In Figure 2, one can see the trend of (30) as $\delta$ increases, which shows that the area under each curve (forbidden area) decreases as the effects due to the fringing field increase. There is a further decrease of the forbidden area when moving towards the edges of the device. However, unlike what was obtained in [51], the modeling of the dielectric properties of the membrane using $f\left(r\right)$ modifies, by increasing the values of $\theta \lambda$, the extension of these areas by reducing the allowed areas on the $r-\theta \lambda$ plane. It follows that, qualitatively, Figure 2 is a valid criterion for choosing the material constituting the membrane when the intended use of the micropump has been identified.

In addition, regarding the uniqueness of the solution for (20), with respect to what was obtained in [51], a new algebraic condition is obtained. Formally, the following result holds.

As in [51], the dependence of (31) on H has an important practical implication. In fact, since the membrane deforms symmetrically with respect to the vertical axis passing through $r=0$, the greatest slope of the membrane profile is found at the edges, and this is confirmed by the numerical recovering performed below. However, unlike what was obtained in [51], here the uniqueness of the solution (governed by (31)) is ensured not only by the presence of $\theta \lambda$ and by the geometry of the micropump but also by the presence of $f\left(r\right)$, confirming that the dielectric properties strongly influence the deformation of the membrane. Therefore, the greater the $f\left(r\right)$ value together with R and $\delta$, the greater V must be in order to overcome the mechanical inertia of the membrane.

To achieve this condition, it is sufficient to solve the algebraic system constituted by both (30) and (31). It is easy to verify that because, on the contrary, one would obtain which represents an implausible result. In fact, taking into account that $H=146$ (as already proven in [51]) when substituting in (33) a plausible value for each parameter, it would follow that $H>{10}^{12}$. Therefore, condition (30) ensures both the existence and uniqueness of the solution for (20) so that any numerical recovery of the membrane that does not satisfy it will give a ghost solution.

For this purpose, we will use the first Lyapunov’s criterion. Therefore, by exploiting two functions, $g_1\left(r\right)$ and $g_2\left(r\right)$, such that model (20) can be rewritten as Therefore, according the fist Lyapunov’s criterion, by imposing ${\displaystyle \frac{d{g}_1\left(r\right)}{dr}}={\displaystyle \frac{d{g}_2\left(r\right)}{dr}}=0$, we achieve the following critical point: Point (36) identifies the position of $u_0$ of the membrane next to the counter-electrode because $d^{*}$ is very small [51]. However, it is easy to see that so that, once V is selected, $p_{el}$ (and consequently p) is not affected by appreciable variation. Thus, by linearizing (35) around (36), since both $g_1\left(r\right)$ and $g_2\left(r\right)$ are analytical functions, and considering the following new variables (with a small enough $ϵ$) the Taylor series made of both $\overline{f}(g_1\left(r\right),g_2\left(r\right))$ and $\overline{g}(g_1\left(r\right),g_2\left(r\right))$ (neglecting the terms of an order higher than the linear one) becomes where $\tau =\sqrt{{\xi }^2+{\eta }^2}$. Furthermore, ${\displaystyle \frac{\partial \overline{f}(g_1^0,g_2^0)}{\partial g_1}}=0$, ${\displaystyle \frac{\partial \overline{f}(g_1^0,g_2^0)}{\partial u_2}}=1$, ${\displaystyle \frac{\partial \overline{g}(g_1^0,g_2^0)}{\partial g_1}}=-r^{-1}$ and ${\displaystyle \frac{\partial \overline{g}(g_1^0,g_2^0)}{\partial u_2}}=0$; thus (39) can be written as whose matrix has no eigenvalues with a positive real part; therefore, (36) represents a stable equilibrium position for (20), whatever $f\left(r\right)$.

As highlighted in [56,57,58,59,60,61], pull-in instability could occur in this kind of device, regardless of $f\left(r\right)$. Therefore, here, we question where there are formulations for $f\left(r\right)$ that affect $u\left(r\right)$, which is usually symmetrical and concave with $u\left(r\right)<1$. Experimentally, the area of the membrane around $u_0$ could present instability because, according to (19), $\left|\mathbf{E}\right|$ is very strong ($1-u\left(r\right)-d^{*}\ll 1$) while on the edge it is very weak ($u\left(r\right)=u\left(R\right)=0$). Thus, we need to formulate $f\left(r\right)$ to reduce $\left|\mathbf{E}\right|$ at $r=0$ by allowing it to be stronger near the edges. A good formulation for $f\left(r\right)$, as a large experimental application highlight, seems to be [58,59,60,61] so that which represents the electrostatic load function for the micropump. However, we question whether specific formulations of $f\left(r\right)$ can affect the numerical solutions (and any multiplicity) without neglecting the pull-in voltage and the stable operating range of the micropump.

Here, we exploit numerical techniques considered the “gold standard” for this type of problem with a singularity of $1/r$ [51]. Particularly, the Keller–Box scheme was implemented in MatLab R2022 together with shooting techniques and the III/IV Lobatto IIIa formulas, which were respectively implemented by exploiting the subroutines ode23, ode45, bpv4c, and bpv5c (mathematical details can be easily found in [51,52]). The Keller–Box procedure, although computationally expensive, offers a better performance. This computational load could create problems in special real-time applications. However, there are very few applications where near-instantaneous membrane recovery is required.

We highlight that, by retracing the procedure used in [51], it is easy to achieve from which, for each numerical approach j, even in the presence of a fringing field and in convergence conditions, the corresponding profile of the membrane, $u_{j,\delta }\left(r\right)$ satisfying (49) does not represent a ghost solution. Therefore, the range of values that ensure convergence without ghost solutions is $\forall j$ and $\forall \delta$, $[{({\left(\theta \lambda \right)}_{conv-noghostsolutions})}_{j,\delta },+\infty )$, where Obviously, if $\theta \lambda >inf{\left(\theta \lambda \right)}_{\mathrm{conv}}$, all of the numerical procedures will converge, but one must pay attention to possible ghost solutions. However, unlike what was obtained in [51], given the presence of $f\left(r\right)$, there is an increase in the values of $\theta \lambda$ which determine the convergence intervals (even as the of $\delta$). Table 1 highlights the range of possible values for $\theta \lambda$ when $\alpha =0.2$, which ensures the convergence of all numerical procedures without ghost solutions (i.e., the numerical solutions obtained satisfy the condition (30)). Indeed, all maximum values of u obtained numerically (denoted by $u_i$) satisfy (30).

From the analysis of the numerical data obtained, agreement with the qualitative results highlighted in Figure 2 is evident. According to this, there is an increase in the values of $\theta \lambda$ compared to what was obtained in [51]. More specifically, the more intense the effects caused by the fringing field, the lower are the minimum values of the intervals of $\theta \lambda$ capable of guaranteeing the convergence of numerical procedures both in the presence and absence of ghost solutions. Our study highlighted how the numerical procedure based on a Keller–Box scheme provides a minimum value of $\theta \lambda$ which guarantees that all numerical procedures converge without the profiles of the recovered membrane representing ghost solutions (see Table 1). Finally, we underline the fact that both Keller–Box and Lobatto’s formulas were set with a number of nodes equal to 45 (below which Keller–Box did not converge, while both Lobatto’s formulas provided inaccurate results) to allow the results to be compared with those set automatically by ode23 and ode45. In Figure 3, Figure 4, Figure 5, Figure 6 and Figure 7, the profiles of the membrane obtained numerically are displayed, increasing $\delta$ in compliance with the range of possible values for $\theta \lambda$ and ensuring the convergence of all the numerical procedures. From these figures, it is clear that, as $\delta$ increases, $u_0$ significantly decreases, as required in (20). Furthermore, the approach used here, reducing $\theta \lambda$ (i.e., reducing externally applied V) establishes regimes of small membrane displacements allowing regular and symmetrical recovery of the membrane with limited deformations, avoiding tearing of the membrane especially at the edges (in [51], the recovery of the membrane showed high entity trapezoidal deformations) favouring the administration of small amount of drugs, as required in specific therapies.

As anticipated above, from (30), with $r<R$, by exploiting (4), it is easy to achieve which represents the externally applied V that is necessary to overcome the mechanical inertia of the membrane ((51) is analogous to the one obtained in [51] but with the presence of $f\left(r\right)$). From (51), it can be seen that the presence of fringing field, by curving the lines of force of $\mathbf{E}$ towards the outside of the micropump, favors the deformation of the membrane, so the device is also recommended for the administration of low-dose drugs [5,6].

From (51) it can be deduced that the choice of the material constituting the membrane (i.e., T) strongly conditions the minimum value of volts to be attributed to the external source in order to overcome the inertia of the membrane itself. We ask ourselves what condition T and V must satisfy to ensure the correct functioning of the device. Therefore, by retracing the same approach used in [51], it is possible to achieve from which the link between T (which represents the material constituting the membrane) and V (which represents the intended use of the micropump) depends on both the electromechanical and dielectric properties of the membrane, $\theta$ and $\lambda$, which influence both the concavity of $u\left(r\right)$ (the greater the T, the lower the concavity) and the convergence of the numerical procedures. This confirms what has already been highlighted in [50], where a similar micropump had been studied in a one-dimensional regime. The results obtained here, with the necessary variations due to $f\left(r\right)$, are completely superimposable since the model has evident radial symmetry. Moreover, from (52), according to [50,51], it can be seen that the highly stressed membranes have a greatly reduced risk of ghost solutions as they allocate the device to applications with low electrical potential values. Particularly, by setting the possible value for each parameter, a graphical representation of (52) is shown in Figure 8, in which both T and V are no longer dimensionless and both allowed and forbidden areas are depicted. Furthermore, from Figure 8, it can be seen that the modeled micropump allows a wide range of uses, since the allowable T of the membrane has an extremely low minimum allowed value. The trend displayed in Figure 8 is similar to that obtained in [52] where the link between T and V assumes the same trend but with more high numerical values because it was considered a MEMS device whose deformable element was a plate.

High values of V (beyond the “pull-in” value) could create instability due to high electrostatic force values that exceed the elastic ones [57,58]. Specifically, if $\lambda >{\lambda }^{*}$, (20) does not provide solutions (for details, see [55,57,58,62,63]). Starting from $\lambda =0$, we traced the trends of $u_0$ to highlight the values of ${\lambda }^{*}$ as a function of $\delta$. The values above ${\lambda }^{*}$ show that (20) has no solution. Figure 9, starting from the absence of a fringing field, highlights the trend of the bifurcation diagram, which has characteristics similar to the experimental evidence [53,55].

To confirm this, Figure 10 displays the trend of the pull-in voltage, $V_{pull-in}$, as a function of the space between the two electrodes. It appears qualitatively superimposable to what is reported in [53]. It should be noted that the increase of d produces a deviation of the pull-in voltage, as $\delta$ varies. Then, as the effects due to the fringing field increase, the pull-in voltage values rise, removing the risk of bifurcations. Moreover, from (4), it is easy to achieve which is very similar to the analogous formulation presented in [53].

Furthermore, from (47) together with (4), it is very easy to achieve $p_{el}$ due to the fringing field effects: With (53), this gives us the link between ${\left(p_{el}\right)}_{fringing}$ and $V_{pull-in}$, as displayed in Figure 11. It is worth noting that (54) is completely analogous to the one obtained in [52,64] where the deformable element was a metallic plate.

Biomedical micropumps for intravenous drug delivery represent an important facility when treatment is reserved for uncooperative patients. Recently, devices have been developed with physics that fall into various fields such as electrostatic micropumps that are of significant importance, since they are easy to make and do not require particular maintenance operations. Physical-mathematical modeling currently exhibits models that realistically describe the behavior of the device in various operating phases. However, given their complexity, they are not suitable for any real-time applications. Circular electrostatic membrane micropumps are no exception, and their analytical models present singularities. In this work, a dimensionless 2D model with second-order differential circular symmetry was presented and studied to model an electrostatic membrane micropump for biomedical purposes, and the effects due to the fringing field (according to Pelesko–Driscoll theory which allows the addition of just one a term via software/hardware) were presented. The modeling of the device, which is voltage controlled by a suitable capacitive electric circuit, expresses $\left|\mathbf{E}\right|$ in terms of the mean curvature of the membrane to allow the use of important results for the study of the existence, uniqueness, and stability of the solution. Therefore, results concerning the existence, uniqueness, and stability of the solution were proved, and the dielectric properties of the material constituting the membrane (here modeled to simulate the capacitive effects during the deformation of the membrane) were exploited. These are strongly linked to the intended use of the device (links to be improved using more sophisticated formulations of curvatures). By using appropriate numerical techniques, the recovery of the membrane was carried out for different electrostatic load conditions, highlighting limitations for some electrical and mechanical quantities, including the electric voltage and the electrostatic force, and highlighting good overlap with experimental evidence in the literature. This allows the recovery of membrane profiles that are advantageous in terms of the safety of the micropump (absence of contact of the membrane with the counter-electrode).

Numerical recovery in our work is a means of recognizing any ghost solutions and understand the link between the properties of the constituent material the membrane and any intended use of the device in the various operating conditions. Thus, it is necessary in the near future to proceed with a FEM analysis of the analytical model under study in order to obtain a software tool that can be easily translated into hardware for any industrial prototypes. In this work, we have limited ourselves to using the numerical procedures since, currently, they are considered a “gold standard” for this type of analytical problems. Furthermore, a criterion for choosing the intended use of the micropump was proposed on the basis of the mechanical tension of the membrane at rest being an important limitation associated with using both V and T to reconstruct $u\left(r\right)$ (model solution) while providing complete safety. The results obtained are comparable with experimental evidence in the literature, especially for cases where small displacements of the membrane are required. The work represents a step of a research project that has the main objective of investigating physical-mathematical models of both membrane and deformable plate microdevices. The cornerstone of this project is the physical-mathematical idea according to which the amplitude of the electric field inside the electrostatic device is locally proportional to the average curvature of the deformable element at the point considered. Specifically, the submitted paper developed the step concerning the 2D modeling of a membrane device in the presence of a modeled fringing field, through one of the most accredited theories developed by Pelesko and Driscoll, with relative numerical recovery of the deformable element and analysis of possible intended uses of the device starting from some properties of the membrane itself. At present, we do not have sufficient elements to define with certainty any practical and specific applications for this type of device. However, we are certain that the proposed device, at least in theory, can be used in biomedical applications in which the external voltage applied is reduced. In the near future, once the dependence of $\delta$ on V has been defined and formalized, it will be desirable to develop more complete models (including the dynamic aspect of the problem and the dependence on the electrical conductivity and on the temperature) so that the adherence with the experimental evidence is more marked.