Final Project in

Introduction to Computational Physics course 118094

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

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

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.

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.

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.

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,

| (1) |

where C_{ii} is the self capacitance of the i-th conductor, C_{i0} is the ground capacitance of the same conductor,
and C_{ij} are the respective cross-capacitances. From the definition of self-capacitance it follows
that

| (2) |

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,

| (3) |

where

| (4) |

is the capacitance matrix.

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

| (5) |

where the integrals are taken over the surfaces S_{j} of all conductors, and only surface charge distributions σ_{j}()
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}() = V _{i} for all points on conductor i. It should also be noted
that

| (6) |

In order to solve (5), we use the standard method of moments (MoM) [4, 5, 6]. 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

| (7) |

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 _{i} of facet i. The last result can be rewritten in vector form
as

| (8) |

where

| (9) |

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 P_{ij} whenever the distance between the facets i and j is large enough. In this
approximation,

| (10) |

where S_{j} is the area of facet j, and _{j} 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 N^{2}. 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,

| (11) |

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

| (12) |

and

| (13) |

The ground capacitance C_{i0} 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 C_{ij} and C_{ji}. In order to insure the matrix is physically correct, it is customary to use the
mean value of C_{ij} and C_{ji} as the cross-capacitance estimate.

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 [5, 4, 9]. 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 [7, 10].

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

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

The numbers of grid divisions in the x, y, and z directions for i-th beam are denoted as Nl_{i}, Nt_{i} and Nw_{i},
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 × Nl_{2} × Nt_{2}.

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 [5, 9]. Figure 4 shows the numbers of important nodes on a given rectangular prism.

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.

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 l_{1} = 100µm, and thickness t_{1} = µm; for the stationary electrode the dimensions are
length l_{2} = 100µm, and thickness t_{2} = 10µm; the width of both conductors is w = 0.5µm. The distance between
the center lines is d = 15µm 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].

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 _{0} 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 _{0} 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.

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

| (14) |

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.,

| (15) |

where C_{0} 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.

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%.

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.

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.

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 C_{0} 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%.

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.