Capacitance of a doubly clamped beam mechanical resonator.

Final Project in
Introduction to Computational Physics course 118094

Stav Zaitsev
Electrical Engineering Department,
Technion - Israel Institute of Technology,
Haifa, Israel

October 26, 2009


1 Introduction
2 Doubly-clamped beam oscillator
 2.1 Device geometry
 2.2 Numerical extraction of capacitances
 2.3 Triangular mesh construction
 2.4 Results
3 Summary
4 Acknowledgment

1 Introduction

Electrostatic excitation is an important method of mechanical actuation of microelectromechanical systems (MEMS). The exciting mechanical force results from the attraction of two conductive surfaces/bodies to which a non-zero voltage difference is applied. This force strongly depends on the mutual capacitance between the two conductors.

Here, we numerically investigate the effective capacitance between two parallel doubly clamped micromechanical beams or plates. Such MEMS devices are widely used as mechanical resonators, mass detectors, parts of optical grids, mechanical switches, etc. In general, analysis of the behavior of such devices requires using the full continuous elasticity model. However, description of doubly clamped beams (as well as cantilevers) as systems with a single degree of freedom is adequate in many engineering applications.

The effective one-dimensional excitation force that exists between two charged beams can be shown to be proportional to the total capacitance between the beams and to the voltage difference squared. If the voltage is constant, the beams’ bend is static. Otherwise, if the voltage is time dependent, the beams may oscillate. In order to take the static or dynamic bending of the beams into account, one may use an appropriate form function describing the geometrical form of each bending beam. Such form function has to satisfy several border conditions. For example, a sine form function assumption is widely used for doubly clamped beams [1].

In this project, we have created a computer program for numerical calculation of the total capacitance between two beams with arbitrary dimensions and form functions. The algorithm uses the method of moments to find the electrical charge distribution on the beams, and calculate the resulting total capacitance. The results for several typical conductor configurations were fitted to approximate analytical models, such as the inverse distance squared law of parallel plate capacitor. We also compare the results for not bent (straight) beams with exact analytical solutions found in literature.

The script and graphical user interface (GUI) were written using Matlab, and run on Windows and Linux platforms. The analytical fits and results presentation were also done in Matlab.

2 Doubly-clamped beam oscillator

2.1 Device geometry

An example of MEMS topology including doubly clamped mechanical beam resonator can be seen on Fig. 1. Such a configuration can be viewed as a coplanar capacitor with the mechanical beam and the stationary electrode acting as two conductive plates [2]. The behavior of the cross-capacitance between these two conductors as function of the beam’s bending is the main interest of this project.


Figure 1: A typical device consists of a suspended doubly clamped narrow beam (length 100 m, width 1-0.25 m, and thickness 0.2 m) and a wide electrode. The excitation force is applied as voltage between the beam and the electrode. (a) Experimental setup and typical sample’s dimensions. The direction of bending of the micromechanical beam is denoted by dotted arrow. (b) SEM picture of the device with one wide electrode and two narrow doubly clamped beams.

Another widely used MEMS design includes a mechanical beam which is clamped on one end, and free on the other. This configuration is known as cantilever. The scripts developed in this project can be readily employed to investigate the electrostatic behavior of a cantilever as well.

2.2 Numerical extraction of capacitances

Self capacitance of a conductor is usually defined as the total charge on this conductor when a unit potential is applied to it, and zero potential is applied to all other conductors [3]. Self capacitance of the conductor consists of the ground capacitance, i.e., capacitance between the conductor and the ground, and cross-capacitances with all other conductors. It is assumed that the potential at infinity is zero. Therefore,

Cii = Ci0 +  Cij,

where Cii is the self capacitance of the i-th conductor, Ci0 is the ground capacitance of the same conductor, and Cij are the respective cross-capacitances. From the definition of self-capacitance it follows that

Qi = CiiVi when  Vj⁄=i = 0,

where V j denotes the potential of j-th conductor. In general, the relation between the conductors’ potentials and charges can be written in matrix form,

(     )     (    )
|  Q1 |     | V1 |
||  Q2 || = Ξ || V2 || ,
(  ...  )     (  ... )
  Qn          Vn


    (                         )
        C11  - C12  ⋅⋅⋅ - C1n
    ||  - C21  C22   ⋅⋅⋅ - C2n ||
Ξ = |(    ..     ..    ..    ..   |)
         .     .      .   .
       - Cn1 - Cn2  ⋅⋅⋅  Cnn

is the capacitance matrix.

In order to calculate Ξ we employ the Laplace equation in its integral form:

            ∑  ∮  σ (x)
4πϵ0ϕi(x0) =     --j----dSj,
             j   |x- x0|

where the integrals are taken over the surfaces Sj of all conductors, and only surface charge distributions σj(x) are assumed to exist, i.e., there are no volume charges. It is also assumed here that the relative electrical permittivity of the medium is unity, i.e. the system is surrounded by vacuum. The potential on the surface of a given conductor is constant, namely, ϕi(x) = V i for all points x on conductor i. It should also be noted that

Qi =   σi(x)dSi.

In order to solve (5), we use the standard method of moments (MoM) [456]. We divide the surfaces of all conductors into triangular facets and assume the charge density to be constant over each facet. The exact method employed for constructing such a triangular mesh is discussed below in Sec. 2.3.

The equation (5) can be approximated using the constructed mesh as

        ∑    ∮   dS
4πϵ0ϕi =    σj  ----j--,
         j     |x- xi|

where index j denotes all facets on all conductors, σj is the constant charge density of facet j, and ϕi is the potential at the center point xi of facet i. The last result can be rewritten in vector form as

    (     )     (    )
       ϕ1         σ1
4πϵ0|(  ...  |) = P |(  ... |) ,
       ϕn         σn


     ∮   dS
Pij =  ----j--
       |x -xi|

is an operator matrix to be calculated analytically and/or numerically.

The integral over a triangular facet which appears on the right side of (9) can be always calculated analytically following equation (19) in Ref. [7]. However, in order to save computation time, we use centroid approximation for Pij whenever the distance between the facets i and j is large enough. In this approximation,

∮   dS        S
   ---j---≈ ----j--,
   |x - xi|   |xj - xi|

where Sj is the area of facet j, and xj is the center point of facet j.

Obviously, the centroid approximation fails for adjacent facets, and the integrand in (9) diverges when i = j. Therefore, we use exact analytical solution for matrix P elements corresponding to close facets and "self-terms", i.e., diagonal terms. The maximal distance between two facet centers at which the analytical result for the corresponding matrix element is used is called close field radius.

In order to increase accuracy, we use seven point Gaussian quadrature as well [8]. The numerical quadrature is used to calculate matrix elements of P for which centroid approximation is not accurate enough, but the numerical quadrature error is insignificant when compared to exact analytical solution. In other words, numerical quadrature is employed for facets whose centers are separated further than defined by close field radius, but closer than the so called far field radius. For matrix elements corresponding to facets separated by a distance larger than the far field radius, the centroid approximation is used. The impact of close and far field radii on the accuracy of the solution will be further investigated in Sec. 2.4.

It is important to bear in mind that if the total number of facets is N, then the total number of elements in P is N2. Therefore, refining the mesh in order to increase accuracy can be very costly in terms of both computer memory and run time.

After the operator matrix P is calculated, the charge density distribution (σ1,n) can be found for any arbitrary distribution of potential over the facets, (ϕ1,n), namely,

(     )           (    )
|  σ1 |           | ϕ1 |
(  ...  ) = 4πϵ0P-1 (  ... ) ,
   σn               ϕn

Note that this operation involves solution of a large linear system of equations which includes inversion of a large matrix P. The solution can be effectively achieved by method of Gaussian reduction, for example.

After calculation of P is completed, we can calculate the self- and cross-capacitances of all conductor bodies by using the following algorithm. First, the potentials of facets belonging to the same conductor body i are set to be 1V while the potentials of all other facets are set to zero. Afterwards, (11) is solved with this potential vector. It follows that

Cii =          Sjσj,


Cij = -          Skσk.

The ground capacitance Ci0 can be now calculated using (1). The whole capacitance matrix Ξ can be found in this manner. Ideally, the matrix Ξ should be symmetrical. However, due to numerical errors, small differences exist between every Cij and Cji. In order to insure the matrix is physically correct, it is customary to use the mean value of Cij and Cji as the cross-capacitance estimate.

2.3 Triangular mesh construction

Building a mesh to cover the surfaces of all conducting bodies in the system is a very important step in capacitance extraction. The quality of the mesh has a tremendous impact on the accuracy of calculations [549]. The meshing process used in our project is described here.

In this project, a triangular mesh is employed. Triangular mesh has several advantages over rectangular and other meshes. First, straightforward algorithms such as Delaunay triangulation [9] exist that create a triangular mesh on the surface of any reasonable geometry. Second, an exact solution of integrals such as the one in (9) over triangular domains exist in literature [710].

At the beginning, two conductors are defined as rectangular prisms, and their faces are numbered as shown in Fig. 2.


Figure 2: The coordinate system, original prism dimensions, and face numbering. Each of two beams is a prism defined by 3 dimensions: li in x direction, ti in y direction, and wi in z direction, where i = 1,2 denotes the beam number.

Each face is divided into triangular facets, as shown in Fig. 3.


Figure 3: Triangular mesh of a single face, its transformation according to form function (in this case - cosine), and numbering order of the triangular facets.

The numbers of grid divisions in the x, y, and z directions for i-th beam are denoted as Nli, Nti and Nwi, respectively. For example, the total number of facets on face number 8, which belongs to the second beam and lies in the x - y plane (see Fig. 2), is 2 × Nl2 × Nt2.

This process results in two data structures - node list and facet list. Node list is an ordered set of coordinate triplets, corresponding to the coordinates of all nodes in the system. Facet list is also an ordered list, each row of which contains the numbers of three nodes that define the corresponding facet [59]. Figure 4 shows the numbers of important nodes on a given rectangular prism.


Figure 4: Numbers of important nodes in the mesh of a single rectangular prism. Numbers denote different prism faces. Nl, Nt and Nw denote the number of grid divisions of the facets in the x, y, and z directions, respectively.

After the mesh of the original two rectangular prism conductors is created, form functions are applied to it, as shown in Figures 3 and 5. Any form function can be employed, allowing one to model any bending or displacement of the mechanical beams. In this work, we apply a cosine form function to the first beam, in order to model bending of doubly clamped beam. We also apply a simple shift in y direction to the second beam, which represents the stationary electrode in our experimental setup, as shown in Fig. 1.


Figure 5: The two beams after form functions were applied. The lower beam has a cosine form function, modeling a bending of a doubly clamped mechanical beam. The form of upper beam remains unchanged, but it is shifted in y direction.

2.4 Results

Below, we present the results of capacitance extraction and their comparison to analytical model for a typical coplanar doubly clamped mechanical beam resonator shown in Fig. 1. The device dimensions for the doubly clamped beam are length l1 = 100m, and thickness t1 = m; for the stationary electrode the dimensions are length l2 = 100m, and thickness t2 = 10m; the width of both conductors is w = 0.5m. The distance between the center lines is d = 15m in y direction.

Both the accuracy of the results and the run time of the scripts strongly depend on the correct choice of close and far field radii values. Therefore, we investigate the Gaussian quadrature and centroid methods of evaluating matrix elements (9), and compare the results to exact analytical solution given in Ref. [7].


Figure 6: The absolute values of relative errors of quadratures calculated using seven point Gaussian quadrature method (black dashed line) and centroid approximation (red solid line). The approximated integral is dS∕|x0 -x| (9), and the surface of integration S is a triangle in x - y plane with nodes (-0.5,-0.5,0), (0,0.5,0), and (0.5,-0.5,0). The position of the test point x0 is swept along the z axis from (0,0,0) to (0,0,50). The horizontal axis represents the distance between x0 and the center of the triangle normalized by the length of the longest side of the triangle, i.e., by √-2-----2
 1 + 0.5 = √ ----
  1.25. Note the logarithmic scales.

The results are shown in Fig. 6. The relative error of both approximate methods is plotted as a function of the distance between the test point x0 and the center of the triangular facet. This distance is normalized by the length of the largest side of the facet, which is used as a measure of the facet’s size. It is evident that if seven point Gaussian quadrature is used at distances which are twice the facet’s size or more, then the relative error is smaller than 10-5. In order to achieve the same amount of error with centroid approximation, the minimal distance from test point x0 to the center of the triangular facet should be 30 times larger than the size of the facet. We employ these results in setting the far and close field radii.

Next, we try to find the optimal number of grid divisions along different dimensions of the beams. The result of refining the mesh in x direction is presented in Fig. 7. It follows that in the case used in this study, the relative difference in the extracted cross-capacitance for Nl = 20 and Nl = 40 is less than 0.2%. The difference between Nl = 4 and Nl = 40 is 1.2%. These results were taken into account while estimating the inaccuracy in extracted capacitances, which are presented below.


Figure 7: Impact of mesh refining on the accuracy of capacitance extraction. Cross-capacitance between two straight parallel beams is calculated as function of the number of grid divisions along the longest dimension of the beams Nl. The dimensions of the beams are: l1 = l2 = 100m, t1 = 1m, t2 = 10m, w1 = w2 = 0.5m. The distance between the center lines is 15m along the y axis. The grid parameters are Nt1 = 2, Nt2 = 4, Nw1 = Nw2 = 2, Nl is the same for both beams and is swept from 4 to 40. For each value of Nl the cross-capacitance C12 is calculated, and normalized by the value of C12 corresponding to Nl = 40.

After establishing the extraction parameters, we calculate the self- and cross-capacitances for a typical configuration of MEMS doubly clamped beam coplanar resonator, sweeping a wide range of bending amplitudes and separation distances between the resonator beam and the stationary electrode. Our goal is to compare these results with simple analytical model of parallel plates capacitor.

The cross-capacitance of two parallel plates in vacuum, whose dimensions are much larger than the distance between them, is

C0 = 0--,

where A is the area of a plate, and d is the separation distance. This capacitance is a function of the separation distance, i.e.,

C (y) =---0-y,
       1- α d

where C0 is the capacitance at original separation distance d, y is the change in the separation distance (positive when the separation distance decreases), and α is a fitting coefficient, which is equal to unity for an ideal parallel plates capacitor.

In order to establish the general correctness of our extraction algorithm, we use it to calculate the capacitance of parallel plates capacitor with finite area. The dimensions of each square plate are set to 1m × 1m, and the plates’ separation value is swept from 1m down to 1cm. The numerical results are compared to the ideal value (14) in Fig. 8.


Figure 8: Comparison of numerically calculated capacitance of finite parallel plates capacitor (blue dots), and an ideal value (14) (green dashed line). The dimensions of the plates are 1m × 1m, and the separation distance d is swept from 1m down to 1cm. The linear size of the plates is denoted by l = 1m. The width of the plates is 1cm. The absolute relative error due to ideal approximation quickly decreases as the ratio between the linear size of the plates l and the plates’ separation distance d increases (red solid line). Note the logarithmic horizontal scale.

The results of the ideal parallel plates model are reasonably accurate when the separation between the plates is at least 100 times smaller than the linear size of the plates. In this case, the error due to fringe electrical field around the edges of the plates is less than 3.5%.

The results of fitting analytical model (15) to numerical cross-capacitance results in a typical case of doubly clamped beam and stationary electrode are shown in Fig. 9 and Fig. 10. Note that if the bending amplitude y does not exceed 10% of the total separation distance d, then the relative error of a simple parallel plates model (α = 1) is less than 3%, while for the modified parallel plates model with fitting coefficient α = 0.88, the maximal relative error is reduced to 1.5%.


Figure 9: Cross-capacitance between a doubly clamped beam and a stationary electrode as function of the normalized bending amplitude y∕d. The dimensions of the beam are length 100m, and thickness 1m; the dimensions of the stationary electrode are length 100m, and thickness 10m; the width of both conductors is 0.5m. The distance between the center lines is d = 15m in y direction. The blue dots represent the cross-capacitance extracted from simulation. Red solid line represents a parallel plates approximation (15) in which d is taken to be equal to the actual center line separation distance, i.e., α = 1. Black dashed line represents the results of a parallel plates model in which α = 0.88, which minimizes the mean square error of (15) in this case.


Figure 10: The relative error of two parallel plates models shown in Fig. 9. Red solid line represents the error of parallel plates approximation (15) in which d is taken to be equal to the actual center line separation distance, i.e., α = 1. Black dashed line represents the error of a parallel plates model in which α = 0.88, which minimizes the mean square error of (15).

Finally, we numerically calculate the cross-capacitance between the beam and the stationary electrode for a wide range of separation distances d and bending amplitudes y. The results are presented in Fig. 11. The extracted values of cross-capacitance are compared with the analytical parallel plates model in Fig. 12, where the relative error between the numerically extracted values and the analytical model with optimal fitting coefficient α = 0.8 is plotted.


Figure 11: Cross-capacitance between a coplanar doubly clamped beam and a stationary electrode for a range of bending amplitudes y = -3 3m, and separation distances d = 12 30m.


Figure 12: The relative error of parallel plates model (15) with fitting coefficient α = 0.8.

According to the results for the particular device geometry studied, if the relative bending amplitude y∕d does not exceed 10%, then the relative error of parallel plates approximation (15) with fitting coefficient α = 0.8 is smaller than 2%. The largest error occurs when the bending amplitude is large enough to be comparable to the distance between the beam and the electrode, in which case the beam and the electrode almost touch, and parallel plates approximation fails.

3 Summary

In this study, we have investigated whether the capacitance of a simple MEMS device consisting of a coplanar doubly clamped beam and a stationary electrode can be adequately described using the model of parallel plates capacitor.

A Matlab script for capacitance extraction was written, including a GUI interface. The charge distribution is calculated using the method of moments. Main sources of numerical errors were investigated, including the finite mesh facet size and numerical quadratures, and the error was shown to be less than 0.5% for the case discussed here. It was also shown that the script correctly predicts the behavior of the capacitance of a finite parallel plates capacitor, and the numerical cross-capacitance value converges to the ideal one when the size of the plates greatly exceeds the distance between them.

Finally, we have shown that the cross-capacitance between a doubly clamped mechanical beam with cosine form function and a stationary electrode can be approximated by the ideal parallel plates capacitance (15), where C0 and d are respectively the cross-capacitance and the distance between the center lines of the electrode and the relaxed (not bent) beam, α = 0.8 0.9, and the bending amplitude y does not exceed 10% of the separation distance d. In this case, the relative error of paralllel plates approximation is less than 2%.

4 Acknowledgment

We would like to thank Dr. Joan Addler for her invaluable help and guidance.


[1]   S. P. Timoshenko. A Course in Strength of Materials. Kiev: Kiev Polytechnic Institute, 1911. In Russian.

[2]   K. C. Gupta, R. Garg, I. Bahl, and P. Bhartia. Microstrip lines and Slotlines. Artech House, 2nd edition, 1996.

[3]   J. D. Jackson. Classical Electrodynamics. Wiley, 1999.

[4]   R. F. Harrington. Field Computation by Moment Methods. Wiley-IEEE Press, 1993.

[5]   W. C. Gibson. The method of moments in electromagnetics. Chapman & Hall/CRC, 2008.

[6]   K. Nabors and J. White. FastCap: a multipole accelerated 3-D capacitance extraction program. IEEE Trans. Computer-Aided Design, 10(11):1447–1459, Nov 1991.

[7]   R. D. Graglia. On the numerical integration of the linear shape functions timesthe 3-D Green’s function or its gradient on a plane triangle. IEEE Trans. Antennas and Propagation, 41:1448–1455, Oct 1993.

[8]   D. A. Dunavant. High degree efficient symmetrical Gaussian quadrature rules for the triangle. Internat. J. Numer. Methods Engrg., 21(6):1129–1148, 1985.

[9]   J. F. Thompson, B. K. Soni, and N. P. Weatherill. Handbook of grid generation. CRC Press, 1999.

[10]   T. F. Eibert and V. Hansen. On the calculation of potential integrals for linear sourcedistributions on triangular domains. IEEE Trans. Antennas and Propagation, 43:1499–1502, Dec 1995.