**3.1 Modal decomposition ROM**

Modes (or resonances) are inherent properties of a structure. Resonances are determined by the material properties (mass, stiffness, and damping properties) and boundary conditions of the structure. Each mode is defined by a *natural (modal or resonant) frequency*, *modal damping*, and a *mode shape*.

Fig. 3. Simple Plate Frequency Response Function

If either the material properties or the boundary conditions of a structure change, its modes will change as well. The overall response of a structure at any frequency is a *summation of responses due to each of its modes*. It is also evident that close to the frequency of one of the resonance peaks, the response of *one mode will dominate the frequency response* [6], fig.3.

Fig. 4. Overlay of Frequency and Time Response Functions

Now if we overlay the time trace with the frequency trace it is possible to notice that maximum oscillations at the time are corresponding to the frequencies of maximum frequency response function (fig.4).

Fig. 5. Flexible Body Modes [6].

158 Microelectromechanical Systems and Devices

Modes (or resonances) are inherent properties of a structure. Resonances are determined by the material properties (mass, stiffness, and damping properties) and boundary conditions of the structure. Each mode is defined by a *natural (modal or resonant) frequency*, *modal* 

If either the material properties or the boundary conditions of a structure change, its modes will change as well. The overall response of a structure at any frequency is a *summation of responses due to each of its modes*. It is also evident that close to the frequency of one of the resonance peaks, the response of *one mode will dominate the frequency response* [6], fig.3.

Now if we overlay the time trace with the frequency trace it is possible to notice that maximum oscillations at the time are corresponding to the frequencies of maximum

**3.1 Modal decomposition ROM** 

Fig. 3. Simple Plate Frequency Response Function

Fig. 4. Overlay of Frequency and Time Response Functions

frequency response function (fig.4).

Fig. 5. Flexible Body Modes [6].

*damping*, and a *mode shape*.

Let's see what happens to the deformation pattern on the structure at each one of these natural frequencies. Modes are further characterized as either *rigid body* or *flexible body*  modes. All structures can have up to *six rigid body* modes, three translational modes and three rotational modes. Many deformation problems are caused, or at least amplified by the excitation of one or more flexible body modes. Fig.5 shows some of the common fundamental (low frequency) modes of a plate. The fundamental modes are given names like those shown in Fig.5. The higher frequency mode shapes are usually more complex in appearance, and therefore don't have common names.

Modal decomposition method uses a weighted sum of *n* mode shapes (modal amplitudes or eigenvalues) *qi*, basis function (an eigenvector) *φi (x, y, z))* of the mechanical structure to represent its deflection *u*:

$$\mu(\mathbf{x}, \mathbf{y}, z, t) = \mu\_{\mathbf{e}\eta} + \sum\_{i=1}^{n} q\_i(t) \varphi\_i(\mathbf{x}, \mathbf{y}, z) \tag{5}$$

where *ueq* is the initial displacement produced by the initial load. For MEMS it is usually sufficient that few modes accurately describe dynamical response of the system. This approach is equivalent to the projection of the original PDE, describing the MEMS behavior, on the subspace defined by the basis functions.

By inserting (5) into the equation of motion (2), taking **x** = φ *eλt* and multiplying by φ*Ti* from the left, and using the orthogonality of eigenvectors, the equation (2) can be reduced to

$$\ddot{q}\_i - 2\sigma\_i q\_i + \left(\alpha\_i^2 + \sigma\_i^2\right) = q\_i^T f \quad , \quad i = 1, 2, \dots \text{n} \tag{6}$$

The eigenvectors satisfy the orthogonality conditions *φTk Mφi= φTk Dφi = φTk Kφi = 0*  for *λk \_*≠ *λi*. With the normalization of the eigenvectors *φTk Mφ<sup>k</sup> ≡ 1*, it can be shown that

$$\mathbf{q}\mathbf{p}^{\mathsf{T}}\_{\mathbf{k}}\mathbf{D}\mathbf{q}\_{\mathbf{k}} = -2\mathbf{o}\_{\mathbf{k}}\text{ and }\mathbf{q}^{\mathsf{T}}{}\_{\mathbf{k}}\mathbf{K}\mathbf{q}\_{\mathbf{k}} = \mathbf{o}\mathbf{o}\_{\mathbf{k}}^{2} + \mathbf{o}^{2}{}\_{\mathbf{k}}\text{ if }\lambda\_{\mathbf{k}} = \mathbf{o}\_{\mathbf{k}}\text{ +j}\mathbf{o}\_{\mathbf{k}}.$$

Note that the orthogonality condition does not apply in the case of multiple eigenvalues. In such case, special modal analysis techniques are needed for decoupling of modes [7]. The decoupling of modes yields that transfer functions can be written as a sum of *modal transfer functions*. Using Fourier transforming expansion and equations (6), the transfer matrix can be derived as

$$H(oo) = \sum\_{i=1}^{n} H\_k(oo) = \sum\_{i=1}^{n} \frac{\rho\_i \rho\_i^T}{(jo - \sigma\_i - jo\_i)(jo - \sigma\_i + jo\_i)}\tag{7}$$

This relation is the basis of modal analysis. It relates the *measurable* transfer functions to the modal properties ωi, σi, and *φi.* Each *i*- th mode contributes with a modal transfer matrix *Hi* to the complete transfer matrix. Fig.6 shows an example of a theoretical transfer function (7) with three *modal peaks* corresponding to modes at 1, 2, and 4 Hz which are all damped with a logarithmic decrement of 1% (*−σi/fi* = 0*.*01). Also the individual modal transfer functions (7) are plotted from which the complete transfer function is computed.

By the way the eigenvalues *qk* and the corresponding eigenvectors *φ<sup>k</sup>* for *k* = 1*,* 2*, . . ., n* can be found as an solution of the modified eigenvalue problem

$$(\mathbf{M}\lambda^2 + \mathbf{D}\lambda + \mathbf{K})\mathbf{q} = \mathbf{0} \,. \tag{8}$$

Macromodels of Micro-Electro-Mechanical Systems (MEMS) 161

1 1 1 1 [ ( , { ,..., })]

where *proj v S* (,) is the projection of the vector *v* into subspace *S*. In other words, we minimize a least squares measure of the "error" distances between the observed states and

*Singular value decomposition (SVD)* takes a rectangular matrix of modal experimental data (defined as *U*, where *U* is a *N* x *n* matrix) in which *N* is a number spatial points of each

UNxn= WNxN SNxn VTnxn , where the columns of *W* are the left singular vectors (*spatial points vectors*); *S* (the same dimensions as *U*) has singular values and is diagonal (*mode amplitudes*); and *VT* has rows that

The *SVD* represents an expansion of the original data in a coordinate system where the covariance matrix is diagonal. Calculating the *SVD* consists of finding the eigenvalues and eigenvectors of matrices *UUT* [*NxN*] and *UTU* [*n x n*]. The eigenvectors of *UTU* make up the columns of *V*, the eigenvectors of *UUT* make up the columns of *W*. Also, the singular values in *S* are square roots of eigenvalues from *UUT* or *UTU*. The singular values are the diagonal entries of the *S* matrix and are arranged in descending order: 1 2 ... 0

singular values are always real numbers. If the matrix *U* is a real matrix, then *W* and *V* are

It was shown [11] that the proper orthogonal basis functions **{***φ***1,** *φ***2,…,** *φ***n}** minimizing (9)

, {1,2,..., } *i i*

After finding basis functions *φk* the linear equations for modal amplitudes *qk* estimation are formed for each temporal snapshot on the base of selected before "snapshot" spatial points

*Karhunen-Loève decomposition (KL)* can be viewed as a statistical procedure [9] .One initially supposes that the observed system dynamics can be modeled as a second-order ergodic stochastic process. The method consists then in constructing a spatial autocorrelation tensor from data obtained through numerical or physical experiments and performing its spectral

*KL* differs from *SVD* in way of finding basis functions *φi* through the temporal snapshot

*ij i mj k*

*a u xt u xt*

1 () () *n i k i k k*

where *n* is number of temporal snapshots and *bk,i* are eigenvectors of the matrix *A*, or the

 *x bux* ,

<sup>1</sup> [ ( , ) ( , )]

1 1

*k m*

*n* 

*n*

*vi n*

 

*i i j*

*u proj u span*

snapshot, *n* is a number of snapshots. The *SVD* theorem states:

are the right singular vectors (*snapshots numbers vectors*).

also real.

can be chosen by setting

and equation (5).

decomposition.

correlation matrix entries

solutions of the equation *Ab=λb*.

as ,

where *vi* is the columns of *V*.

the basis function representation of those states.

*n n <sup>n</sup> <sup>T</sup> ii n i j i*

2

 *u*

 

 

*<sup>n</sup>* .The

, (9)

The relationships between natural frequencies *fk*, logarithmic decrements *δ<sup>k</sup>* and the eigenvalues are *fk* = 2*πω<sup>k</sup>* and *δ<sup>k</sup>* = *−σk/fk .* 

Fig. 6. Example of the modal decoupling in a transfer function [7].

#### **3.1.1 Modal ROM based on proper orthogonal decomposition methods**

Reduced order modeling using modal basis functions was originally developed by [11-13] and has been continuously improved by several authors. The basis functions (,,) *<sup>i</sup> xyz* in (5) may be chosen in two ways by:


The choice of orthogonal basis functions *φ<sup>k</sup>* can be done by the following way [8]. First the MEMS dynamics are simulated using a slow but accurate technique such as FEM or FDM. Sets of runs may be used to suitably characterize the operating range of the device. The spatial distributions of each state variable *u(x, y, z, t)* are then sampled at a series of *Ns* different times during these simulations, and the sampled distributions are stored as a series of vectors, {ui}, where each of them corresponds to a particular "snapshot" in time. Now suppose we would like to pick orthogonal basic *n* functions **{***φ***1,** *φ***2,…,** *φ***n}** in order to represent the observed state distributions as closely as possible. One way to do this is to attempt to minimize a least squares measure of the "error" distances between the observed states and the basis function representation of those states:

The relationships between natural frequencies *fk*, logarithmic decrements *δ<sup>k</sup>* and the

eigenvalues are *fk* = 2*πω<sup>k</sup>* and *δ<sup>k</sup>* = *−σk/fk .* 

may be chosen in two ways by:

Fig. 6. Example of the modal decoupling in a transfer function [7].

particularly when irregular geometries are involved.

mode shapes of the device structural elements.

states and the basis function representation of those states:

**3.1.1 Modal ROM based on proper orthogonal decomposition methods** 

Reduced order modeling using modal basis functions was originally developed by [11-13] and has been continuously improved by several authors. The basis functions (,,) *<sup>i</sup>*

 Using the undamped linear mode shapes of the undeflected microstructure as basis functions. For simple structures with simple boundary conditions, the mode shapes can be found analytically. For complex structures or complex boundary conditions, the linear mode shapes are obtained numerically using the finite element method. However, it is usually difficult to determine, *a priori*, an optimum set of eigenvectors *φk*,

 Using snapshots obtained from experiments under a training signal or full FEM model runs. Then the proper orthogonal decomposition method (Singular Value Decomposition –*SVD* , Karhunen -Loève decomposition -*KL* and neural networks-based generalized Hebbian algorithm -*GHA*) is applied to the time series for extracting the

The choice of orthogonal basis functions *φ<sup>k</sup>* can be done by the following way [8]. First the MEMS dynamics are simulated using a slow but accurate technique such as FEM or FDM. Sets of runs may be used to suitably characterize the operating range of the device. The spatial distributions of each state variable *u(x, y, z, t)* are then sampled at a series of *Ns* different times during these simulations, and the sampled distributions are stored as a series of vectors, {ui}, where each of them corresponds to a particular "snapshot" in time. Now suppose we would like to pick orthogonal basic *n* functions **{***φ***1,** *φ***2,…,** *φ***n}** in order to represent the observed state distributions as closely as possible. One way to do this is to attempt to minimize a least squares measure of the "error" distances between the observed

*xyz* in (5)

$$\sum\_{i=1}^{n} \left[ \boldsymbol{\mu}\_{i} - \operatorname{proj} \{ \boldsymbol{\mu}\_{i}, \operatorname{span} \{ \boldsymbol{\varphi}\_{1}, \dots, \boldsymbol{\varphi}\_{n} \} \} \right]^{2} = \sum\_{i=1}^{n} \sum\_{j=1}^{n} \boldsymbol{\varphi}\_{i} \boldsymbol{\varphi}\_{j}^{T} \boldsymbol{\mu}\_{i} \tag{9}$$

where *proj v S* (,) is the projection of the vector *v* into subspace *S*. In other words, we minimize a least squares measure of the "error" distances between the observed states and the basis function representation of those states.

*Singular value decomposition (SVD)* takes a rectangular matrix of modal experimental data (defined as *U*, where *U* is a *N* x *n* matrix) in which *N* is a number spatial points of each snapshot, *n* is a number of snapshots. The *SVD* theorem states:

$$\mathbf{U\_{Nxn}} = \mathbf{W\_{NxN}} \, \mathbf{S\_{Nxn}} \, \mathbf{V^{T}\_{nxn}} \, \iota$$

where the columns of *W* are the left singular vectors (*spatial points vectors*); *S* (the same dimensions as *U*) has singular values and is diagonal (*mode amplitudes*); and *VT* has rows that are the right singular vectors (*snapshots numbers vectors*).

The *SVD* represents an expansion of the original data in a coordinate system where the covariance matrix is diagonal. Calculating the *SVD* consists of finding the eigenvalues and eigenvectors of matrices *UUT* [*NxN*] and *UTU* [*n x n*]. The eigenvectors of *UTU* make up the columns of *V*, the eigenvectors of *UUT* make up the columns of *W*. Also, the singular values in *S* are square roots of eigenvalues from *UUT* or *UTU*. The singular values are the diagonal entries of the *S* matrix and are arranged in descending order: 1 2 ... 0 *<sup>n</sup>* .The singular values are always real numbers. If the matrix *U* is a real matrix, then *W* and *V* are also real.

It was shown [11] that the proper orthogonal basis functions **{***φ***1,** *φ***2,…,** *φ***n}** minimizing (9) can be chosen by setting

$$\phi\_i = \upsilon\_{i\prime} \text{ i } = \{1, 2, \ldots, n\}$$

where *vi* is the columns of *V*.

After finding basis functions *φk* the linear equations for modal amplitudes *qk* estimation are formed for each temporal snapshot on the base of selected before "snapshot" spatial points and equation (5).

*Karhunen-Loève decomposition (KL)* can be viewed as a statistical procedure [9] .One initially supposes that the observed system dynamics can be modeled as a second-order ergodic stochastic process. The method consists then in constructing a spatial autocorrelation tensor from data obtained through numerical or physical experiments and performing its spectral decomposition.

*KL* differs from *SVD* in way of finding basis functions *φi* through the temporal snapshot

$$\text{Correlation matrix entries} \qquad \qquad a\_{\vec{\boldsymbol{\eta}}} = \frac{1}{n} \sum\_{k=1 \atop m=1}^{n} [\boldsymbol{\mu}\_{\boldsymbol{i}}(\mathbf{x}, t\_{m}) \boldsymbol{\mu}\_{\boldsymbol{j}}(\mathbf{x}, t\_{k})]^{\top}$$

as ,

where *n* is number of temporal snapshots and *bk,i* are eigenvectors of the matrix *A*, or the solutions of the equation *Ab=λb*.

1 () () *n i k i k k*

 *x bux* ,

Macromodels of Micro-Electro-Mechanical Systems (MEMS) 163

set of full finite element simulations [13]. A similar procedure is proposed in [17]. In this case, the force-displacement function and the modal strain energy are still derived from a series of FEM simulations, but a polynomial multi-variable fit is used. In order to reduce the complexity of the fitting step, dominant and relevant modes are first characterized. Both the procedures in [13] and [17] can be partially automated and extended to include other conservative energy domains. Other algorithms have also been proposed for the approximation of dissipative energy terms [18, 19]. Dedicated methods have been demonstrated for actuated microbeams, still using the linear undamped mode shapes of the device as basic functions in the Galerkin procedure [20], where two expressions of the nonlinear electrostatic term were proposed as a function of modal coordinates, each including all nonlinearities up to the fifth order, obtained via mathematical manipulation [21]. A new fuzzy-logic model (FLM) for MEMS is presented in [22] in which for reducing the number of data needed for macromodel identification, cluster estimation of a model structure and back propagation method of structure parameters adaptation are chosen to fit the data. As a result the dynamic coupled simulation of a magnetic microactuator takes only several minutes and the force macromodel yielded errors is less than 1.5% for a 5-*μ*m displacement. These FLMs combine fuzzy sets with fuzzy rules that have the capability to

(a)

(b)

Fig. 8. Basis functions for (a) displacement *z(x, t)* and (b) pressure *p(x, y, t)* [11]

model the complex nonlinear behavior.

*Generalized Hebbian algorithm (GHA)* is an Artificial Neural Network (ANN) approach of performing Principal Component Analysis (PCA) on a set of data and can be also used as a learning procedure for the approximation of *PDE* solutions by the expression (5) [10] .

The modal approach to MEMS macromodeling is illustrated by sensor device (fig.7), which is described by coupling a 1-D elastic beam equation with electrostatic force and 2-D compressible isothermal squeeze-film Reynold's equation [11]. This device consists of a deformable elastic beam microstructure that is electrostatically pulled in by an applied voltage waveform. The dynamics of the beam are first simulated using a finite-difference analysis. A quarter of the beam is initially meshed using a 20 x10 node 2-D grid.

Fig. 7. Fixed-fixed beam pull-in time pressure sensor device [11]

The state at each node consists of three quantities: *z, dz/dt* and *p*. Since *z* and *dz/dt* are simulated in 1-D, this results in coupled nonlinear ODE which must be integrated in time. Basis functions are generated for pressure *p* and displacement *z* based on runs of the finitedifference code for an ensemble of four different step voltages: 9V, 10V, 12V, 16V. One hundred samples of pressure and displacement are taken during these four runs at fixed time intervals. These samples are used to generate the basic functions. The resulting basis functions for displacement and pressure are shown in Fig.8 [11].

#### **3.1.2 Nonlinear modal ROM**

The device governing equations are generally derived from Lagrange equations, after expressing the internal (elastic and kinetic) and external (electrostatic) energy of the system in terms of modal amplitudes and symbolical calculation of the gradients. Assuming that the device undergoes small displacement, the basis chosen results in diagonal mass and stiffness matrices, which can be pre-computed. Linear elastic undamped normal modes of the undeflected device have been often chosen as basic functions to approximate the solution of an electromechanical problem discretized using finite element methods. Modal representation is very efficient since it *requires only one equation per mode* and involved conductor to describe the entire system.

The nonlinear energy terms, instead, are generally expressed as analytical functions of the modal coordinate. In [13] a single static full finite element simulation is used to determine the number of modal functions needed to capture the device behavior and their expected amplitude. Then this information is used to construct the electrostatic energy term. A 3D full model electrostatic simulation is run for values of the modal amplitude that span the operating range of the device, in order to compute the capacitance/deformation curve by using ANSYS' transducer element TRANS126*.* The results are fitted with a rational fraction of multivariate polynomials using a nonlinear function fitting scheme. In order to model large-displacement behavior of the device, the strain energy is derived by fitting data from a

*Generalized Hebbian algorithm (GHA)* is an Artificial Neural Network (ANN) approach of performing Principal Component Analysis (PCA) on a set of data and can be also used as a learning procedure for the approximation of *PDE* solutions by the expression (5) [10] . The modal approach to MEMS macromodeling is illustrated by sensor device (fig.7), which is described by coupling a 1-D elastic beam equation with electrostatic force and 2-D compressible isothermal squeeze-film Reynold's equation [11]. This device consists of a deformable elastic beam microstructure that is electrostatically pulled in by an applied voltage waveform. The dynamics of the beam are first simulated using a finite-difference

The state at each node consists of three quantities: *z, dz/dt* and *p*. Since *z* and *dz/dt* are simulated in 1-D, this results in coupled nonlinear ODE which must be integrated in time. Basis functions are generated for pressure *p* and displacement *z* based on runs of the finitedifference code for an ensemble of four different step voltages: 9V, 10V, 12V, 16V. One hundred samples of pressure and displacement are taken during these four runs at fixed time intervals. These samples are used to generate the basic functions. The resulting basis

The device governing equations are generally derived from Lagrange equations, after expressing the internal (elastic and kinetic) and external (electrostatic) energy of the system in terms of modal amplitudes and symbolical calculation of the gradients. Assuming that the device undergoes small displacement, the basis chosen results in diagonal mass and stiffness matrices, which can be pre-computed. Linear elastic undamped normal modes of the undeflected device have been often chosen as basic functions to approximate the solution of an electromechanical problem discretized using finite element methods. Modal representation is very efficient since it *requires only one equation per mode* and involved

The nonlinear energy terms, instead, are generally expressed as analytical functions of the modal coordinate. In [13] a single static full finite element simulation is used to determine the number of modal functions needed to capture the device behavior and their expected amplitude. Then this information is used to construct the electrostatic energy term. A 3D full model electrostatic simulation is run for values of the modal amplitude that span the operating range of the device, in order to compute the capacitance/deformation curve by using ANSYS' transducer element TRANS126*.* The results are fitted with a rational fraction of multivariate polynomials using a nonlinear function fitting scheme. In order to model large-displacement behavior of the device, the strain energy is derived by fitting data from a

analysis. A quarter of the beam is initially meshed using a 20 x10 node 2-D grid.

Fig. 7. Fixed-fixed beam pull-in time pressure sensor device [11]

functions for displacement and pressure are shown in Fig.8 [11].

**3.1.2 Nonlinear modal ROM** 

conductor to describe the entire system.

set of full finite element simulations [13]. A similar procedure is proposed in [17]. In this case, the force-displacement function and the modal strain energy are still derived from a series of FEM simulations, but a polynomial multi-variable fit is used. In order to reduce the complexity of the fitting step, dominant and relevant modes are first characterized. Both the procedures in [13] and [17] can be partially automated and extended to include other conservative energy domains. Other algorithms have also been proposed for the approximation of dissipative energy terms [18, 19]. Dedicated methods have been demonstrated for actuated microbeams, still using the linear undamped mode shapes of the device as basic functions in the Galerkin procedure [20], where two expressions of the nonlinear electrostatic term were proposed as a function of modal coordinates, each including all nonlinearities up to the fifth order, obtained via mathematical manipulation [21]. A new fuzzy-logic model (FLM) for MEMS is presented in [22] in which for reducing the number of data needed for macromodel identification, cluster estimation of a model structure and back propagation method of structure parameters adaptation are chosen to fit the data. As a result the dynamic coupled simulation of a magnetic microactuator takes only several minutes and the force macromodel yielded errors is less than 1.5% for a 5-*μ*m displacement. These FLMs combine fuzzy sets with fuzzy rules that have the capability to model the complex nonlinear behavior.

Fig. 8. Basis functions for (a) displacement *z(x, t)* and (b) pressure *p(x, y, t)* [11]

Macromodels of Micro-Electro-Mechanical Systems (MEMS) 165

stroke is understood as modal amplitude. Damping parameters are assigned to each

The first step of the ROM generation is to determine which modes are really significant, and to estimate a proper amplitude range for each mode. Several criteria can be applied, for instance, the lowest eigenmodes of a modal analysis, modes in operating direction, or modes, which contribute to the deflection state at a typical test load. Next the dependencies of the strain energy *Wst* and capacities are described by polynomial functions being fitted. The necessary data points are obtained by imposing each eigenmode with varying amplitude on the mechanical mode1l for the non-linear strain energy and on an electrostatic space model for capacitance. This process is computationally expensive but has to be done just once. The result is a black-box model that can be applied to any load situation. In the concept of the modal superposition method, each eigenmode represents a single

The export of the ROM to VHDL-AMS is performed in two steps [16]. At first, an initialization file containing all necessary information of the macromodel, such as the fitted polynomial coefficients and orders, is generated. Then, the source code in VHDL-AMS is automatically generated. The main problem of exporting the ROM in VHDL-AMS is to express the fitted functions of the non-linear strain energy and of the capacities which are part of coefficients of the differential algebraic equations (DAE), which can be mapped to the simultaneous statements of VHDL-AMS. If simulators support description in matrix notation properly (as MATLAB), the exported VHDL-Models will become more compact

The Modal ROM approach was implemented as the available ROM-Tool in ANSYS/Multiphysics since Release 7 **(ROM144).** It contains some terminals (fig. 9) and




independent resonator with modal mass *mi* and modal damping *ξi*.

eigenmode.

and clear.

provides necessary functions:

*FN,i* at these nodes;

modes;

of the system.

Fig. 9. ROM144 functional block

Macromodels obtained via modal basis functions methods have been demonstrated to reproduce results obtained with full physical level simulation with an accuracy of some percentage points and a reduction of the computational complexity simulation. It was shown that even when the problem is mechanically nonlinear, the linear normal modes can serve as basic functions. The approach can be outlined as follows [13]:


#### **3.1.3 VHDL-AMS export of modal ROM**

In general, the equation (5) describes a coordinate transformation of finite element displacement coordinates (mesh element coordinate) to modal coordinates of the macromodel (basic functions or a degree of freedom-d.o.f.):

$$
\mu(\mathbf{x}, \mathbf{y}, z, t) = \mu\_{e\eta} + \sum\_{i=1}^{n} q\_i(t) \phi\_i(\mathbf{x}, \mathbf{y}, z) \ .
$$

The deformation state of the structure given by *N* nodal displacements *ui* (*i=1,2,*…*,N* ) is now represented by a linear combination of *n* modes weighted by their amplitudes *qj* (*j=1,2,*…*,n* ) where *n* << *N*.

The governing equation of motion describing the ROM of electrostatic actuated MEMS structures in modal coordinates:

$$\rho m\_j \ddot{\boldsymbol{\eta}}\_j + 2 \boldsymbol{\xi}\_j \alpha\_j \boldsymbol{m}\_j \dot{\boldsymbol{\eta}}\_j + \frac{\partial}{\partial \boldsymbol{\eta}\_j} \boldsymbol{W}\_{\text{st}} (\boldsymbol{q}\_1, \dots, \boldsymbol{q}\_n) = \frac{1}{2} \sum\_r \frac{\partial}{\partial \boldsymbol{\eta}\_j} \boldsymbol{C}\_{\text{ks}} (\boldsymbol{q}\_1, \dots, \boldsymbol{q}\_n) \cdot (\boldsymbol{V}\_k - \boldsymbol{V}\_s)^2 + \sum\_{i=1}^n \boldsymbol{\varrho}\_i \boldsymbol{f}\_i \tag{10}$$

where *mj* is the modal mass, ω*j* is the eigenfrequency, ξ*<sup>j</sup>* the linear modal damping ratio, *Wst* is the modal strain energy function, *Cks* is the modal capacity-stroke function, *r* is the number of capacities involved for microsystems with multiple electrodes, *V* is the electrode voltage applied and *fi* is a local force acting at the *i*-th node. The current *Ik* at each electrode *k*  is defined by:

$$I\_k = \frac{\partial Q\_k}{\partial t} = \sum\_r \left( \mathbf{C}\_{ks} \cdot (\frac{\partial V\_k}{\partial t} - \frac{\partial V\_s}{\partial t}) + \frac{\partial \mathbf{C}\_{ks}}{\partial t} \cdot (V\_k - V\_s) \right) \tag{11}$$

An essential prerequisite to establish (10) and (11) are proper modal strain energy and capacity-stroke functions. Both are derived from a series of FEM runs at various deflection states in the operating range. The received data are used for polynomial functions fitting in order to compute the local derivatives, which describe force and stiffness terms. As a matter of fact, shape function methods can be applied to nonlinear systems, too [13]. Geometric nonlinearities, as, for instance, stress stiffening, can be regarded if the modal stiffness is computed from the first derivative of the strain energy function with respect to the modal amplitudes. Capacitance-stroke functions provide non-linear coupling between each eigenmode and the electrical quantities (i.e. electrostatic modal forces, electrical current) if

Macromodels obtained via modal basis functions methods have been demonstrated to reproduce results obtained with full physical level simulation with an accuracy of some percentage points and a reduction of the computational complexity simulation. It was shown that even when the problem is mechanically nonlinear, the linear normal modes can

2. Substitute *u (x ,y, z ,t)* in the governing equation for the deflection (e.g., the Euler–

3. Obtain a system of *n*-coupled second-order ordinary differential equations for the *qi(t).* 4. Solve the equations to compute the dynamic response either numerically or

In general, the equation (5) describes a coordinate transformation of finite element displacement coordinates (mesh element coordinate) to modal coordinates of the

> 1 ( , , ,) () ( , , ) *n eq i i i uxyzt u q t xyz*

The deformation state of the structure given by *N* nodal displacements *ui* (*i=1,2,*…*,N* ) is now represented by a linear combination of *n* modes weighted by their amplitudes *qj*

The governing equation of motion describing the ROM of electrostatic actuated MEMS

1 1

*jj j j jj st n ks n k s i i*

where *mj* is the modal mass, ω*j* is the eigenfrequency, ξ*<sup>j</sup>* the linear modal damping ratio, *Wst* is the modal strain energy function, *Cks* is the modal capacity-stroke function, *r* is the number of capacities involved for microsystems with multiple electrodes, *V* is the electrode voltage applied and *fi* is a local force acting at the *i*-th node. The current *Ik* at each electrode *k* 

( ( ) ( )) *k k s ks k ks <sup>k</sup> <sup>s</sup>*

An essential prerequisite to establish (10) and (11) are proper modal strain energy and capacity-stroke functions. Both are derived from a series of FEM runs at various deflection states in the operating range. The received data are used for polynomial functions fitting in order to compute the local derivatives, which describe force and stiffness terms. As a matter of fact, shape function methods can be applied to nonlinear systems, too [13]. Geometric nonlinearities, as, for instance, stress stiffening, can be regarded if the modal stiffness is computed from the first derivative of the strain energy function with respect to the modal amplitudes. Capacitance-stroke functions provide non-linear coupling between each eigenmode and the electrical quantities (i.e. electrostatic modal forces, electrical current) if

*Q VV C I C V V t tt t* 

<sup>1</sup> <sup>2</sup> ( ,..., ) ( ,.., ) ( ) <sup>2</sup>

*q q*

*r*

*m q m q W q q C q q V V f*

*j j r i*

(11)

, (10)

2

1

*n*

.

serve as basic functions. The approach can be outlined as follows [13]:

1. Compute the linear modes *φi* of the elastic problem.

Bernoulli or Timoshenko equations for beams).

macromodel (basic functions or a degree of freedom-d.o.f.):

**3.1.3 VHDL-AMS export of modal ROM** 

analytically.

(*j=1,2,*…*,n* ) where *n* << *N*.

is defined by:

structures in modal coordinates:

  stroke is understood as modal amplitude. Damping parameters are assigned to each eigenmode.

The first step of the ROM generation is to determine which modes are really significant, and to estimate a proper amplitude range for each mode. Several criteria can be applied, for instance, the lowest eigenmodes of a modal analysis, modes in operating direction, or modes, which contribute to the deflection state at a typical test load. Next the dependencies of the strain energy *Wst* and capacities are described by polynomial functions being fitted. The necessary data points are obtained by imposing each eigenmode with varying amplitude on the mechanical mode1l for the non-linear strain energy and on an electrostatic space model for capacitance. This process is computationally expensive but has to be done just once. The result is a black-box model that can be applied to any load situation. In the concept of the modal superposition method, each eigenmode represents a single independent resonator with modal mass *mi* and modal damping *ξi*.

The export of the ROM to VHDL-AMS is performed in two steps [16]. At first, an initialization file containing all necessary information of the macromodel, such as the fitted polynomial coefficients and orders, is generated. Then, the source code in VHDL-AMS is automatically generated. The main problem of exporting the ROM in VHDL-AMS is to express the fitted functions of the non-linear strain energy and of the capacities which are part of coefficients of the differential algebraic equations (DAE), which can be mapped to the simultaneous statements of VHDL-AMS. If simulators support description in matrix notation properly (as MATLAB), the exported VHDL-Models will become more compact and clear.

The Modal ROM approach was implemented as the available ROM-Tool in ANSYS/Multiphysics since Release 7 **(ROM144).** It contains some terminals (fig. 9) and provides necessary functions:


Fig. 9. ROM144 functional block

Macromodels of Micro-Electro-Mechanical Systems (MEMS) 167

The scaling function *φ*(*x*) is supported in the interval [0*, L* − 1] while the corresponding wavelet *ψ*(*x*) is supported in the interval [1 − *L/*2*, L/*2]. The parameter *L* will be referred to as the degree of the scaling function *φ*(*x*). The coefficients *pi* are called the wavelet filter coefficients and they satisfy the same conditions. The constructed scaling function *φ*(*x*) and

The autocorrelation functions *θ*(*x*), which are used for generating basic functions, can be

() ()( ) *x* 

, ( ) (2 ) *<sup>J</sup>*

Wavelets have proven to be an efficient tool of analysis in many fields including the solution of PDE. In [23] a new wavelet interpolation Galerkin method is used for the numerical simulation of MEMS devices under the effect of squeeze film damping. The air film pressure is expressed as a linear combination of a class of basic functions generated by autocorrelation of the usual compactly supported Daubechies scaling functions, which are the first- generation wavelets. The wavelet interpolation Galerkin method was used to predict the frequency response of the accelerometer with a spring mass-damper model with a parallel-plate electrostatic force. Various numerical tests have been conducted by changing the degree of the Daubechies wavelet *L* and the number *J* of the scale. Better accuracy can be achieved by increasing *L* and *J*. The higher *L* is, the smoother the scaling function becomes. The price for the high smoothness is that its supporting domain gets larger. The higher *J* is, the more accurate the solution becomes. The number of differential equations and the CPU time increase significantly as *J* increases. The solutions for *L* = 6 and *J* = 4 have results higher

The present wavelet interpolation Galerkin method is not suitable to solve problems defined on nonrectangular domains, since higher-dimensional wavelets are constructed by employing the tensor product of the one-dimensional wavelets, so their application is restricted to rectangular domains. But usage of the second-generation wavelets which are constructed in the spatial domain can expend in future PDE solutions to complex domains.

Moment matching model order reduction is based on the approximation of the original *n*dimensional system transfer function *F(s)* of the original *n*-dimensional system with a rational function with a lower degree *q<<n* [24]. This is done by matching some terms of the

() () ()

() () *<sup>T</sup> X t AX t Bv t*

*Yt C Xt*

let's perform a Laplace transform to obtain its frequency domain transfer function

 

*x d*

 

(15)

*x xk* . (16)

(17)

<sup>1</sup> () ( ) *<sup>T</sup> F s C sI A B* (18)

*J k* 

wavelet *ψ*(*x*) have also the prescribed properties.

defined as follows:

and act as the scaling function

than the finite difference method.

**3.2 Moment matching based ROM** 

For a state space equation

Taylor expansion of *F(s)* around a certain expansion point*.*

The ROM144 takes input equations (5),(10),(11) with characteristic constants ( modal masses, modal damping ratios and eigenvectors of the master nodes) and provides the special functions of calculating the strain energy *Wmech (qi, qj , qk)* and capacitances *Cop(qi, qj , qk )*as well as their first derivatives with respect to the modal amplitudes *qi, qj* and *qk* using the information of the polynomials degrees defined in other packages.

Modal ROM144 approach speeds up computations in 40 times in comparison with the FEM model while pull-in time errors is less than 2%.

Modal ROM approach was implemented also in *INTEGRATOR system of CoventorWare*  (http:// www.coventor.com)*, MEMS Pro* (http://www.memspo.com) and *MEMSCAP*  (http:// www.memsscap.com)*.*

#### **3.1.4 Galerkin's approximated ROM**

As in (5) the desired PDE solution *u(x,y,z,t)* can be approximated by spatially varying arbitral basis functions (,,) *<sup>i</sup> xyz* with time varying coefficients : ( ) *<sup>i</sup> t*

$$\mu(\mathbf{x}, y, z, t) = \sum\_{i=1}^{n} \alpha\_i(t) \beta\_i(\mathbf{x}, y, z) \tag{12}$$

For the Galerkin's method the PDE residual (() ) *Lu f* is orthogonal to each *α<sup>i</sup>* of the basic functions in the operating range *H*:

$$\mathbb{E}\left(\alpha\_i, \mathbb{L}(\mu) - f\right) = \int \alpha\_i^T \left(\mathbb{L}(\mu) - f\right) dt = 0, \quad i = 1, n. \tag{13}$$

where *L* is a differential operator (possibly nonlinear), and *f* is an input vector.

The basic functions (,,) *<sup>i</sup> xyz* can be chosen arbitrarily, as long as their elements satisfy all of the boundary conditions and are sufficiently differentiable. It means that they can be not only eigenmodes as it was shown before, but they may be Tchebychev, Legendre, Hermit polynomials or even wavelets functions, which were introduced in the past two decades and are gaining increasing popularity. Indeed wavelets have many excellent properties: such as orthogonality, compact support, exact representation of polynomials to a certain degree, and flexibility to represent functions at different resolution levels. The wavelet-Galerkin method is a Galerkin scheme using wavelet functions as the basic functions. However, wavelet functions do not satisfy the boundary conditions. Thus the treatment of general boundary conditions is a major difficulty for the application of the wavelet- Galerkin method. The wavelet interpolation Galerkin method is a Galerkin scheme where basic functions are a class of interpolating functions generated by autocorrelation of the usual compactly supported Daubechies scaling functions [23]. Daubechies' functions are easy to construct. For an even integer *L*, we have the Daubechies' scaling function *φ*(*x*) and wavelet *ψ*(*x*) satisfying

$$\begin{aligned} \varphi(\mathbf{x}) &= \sum\_{i=1}^{L-1} p\_i \varphi(2\mathbf{x} - i) \\ \varphi(\mathbf{x}) &= \sum\_{i=2-L}^{1} (-1)^i p\_{1-i} \varphi(2\mathbf{x} - i) \end{aligned} \tag{14}$$

The scaling function *φ*(*x*) is supported in the interval [0*, L* − 1] while the corresponding wavelet *ψ*(*x*) is supported in the interval [1 − *L/*2*, L/*2]. The parameter *L* will be referred to as the degree of the scaling function *φ*(*x*). The coefficients *pi* are called the wavelet filter coefficients and they satisfy the same conditions. The constructed scaling function *φ*(*x*) and wavelet *ψ*(*x*) have also the prescribed properties.

The autocorrelation functions *θ*(*x*), which are used for generating basic functions, can be defined as follows:

$$
\partial \theta(\mathbf{x}) = \int\_{-\infty}^{\infty} \rho(\mathbf{r}) \rho(\mathbf{r} - \mathbf{x}) d\mathbf{r} \tag{15}
$$

and act as the scaling function

166 Microelectromechanical Systems and Devices

The ROM144 takes input equations (5),(10),(11) with characteristic constants ( modal masses, modal damping ratios and eigenvectors of the master nodes) and provides the special functions of calculating the strain energy *Wmech (qi, qj , qk)* and capacitances *Cop(qi, qj , qk )*as well as their first derivatives with respect to the modal amplitudes *qi, qj* and *qk* using the

Modal ROM144 approach speeds up computations in 40 times in comparison with the FEM

Modal ROM approach was implemented also in *INTEGRATOR system of CoventorWare*  (http:// www.coventor.com)*, MEMS Pro* (http://www.memspo.com) and *MEMSCAP* 

As in (5) the desired PDE solution *u(x,y,z,t)* can be approximated by spatially varying

1 ( , , ,) () ( , , ) *n*

For the Galerkin's method the PDE residual (() ) *Lu f* is orthogonal to each *α<sup>i</sup>* of the basic

( , ( ) ) ( ( ) ) 0, 1, . *<sup>T</sup>*

of the boundary conditions and are sufficiently differentiable. It means that they can be not only eigenmodes as it was shown before, but they may be Tchebychev, Legendre, Hermit polynomials or even wavelets functions, which were introduced in the past two decades and are gaining increasing popularity. Indeed wavelets have many excellent properties: such as orthogonality, compact support, exact representation of polynomials to a certain degree, and flexibility to represent functions at different resolution levels. The wavelet-Galerkin method is a Galerkin scheme using wavelet functions as the basic functions. However, wavelet functions do not satisfy the boundary conditions. Thus the treatment of general boundary conditions is a major difficulty for the application of the wavelet- Galerkin method. The wavelet interpolation Galerkin method is a Galerkin scheme where basic functions are a class of interpolating functions generated by autocorrelation of the usual compactly supported Daubechies scaling functions [23]. Daubechies' functions are easy to construct. For an even integer *L*, we have the Daubechies' scaling function *φ*(*x*) and wavelet

 

1

*L*

( ) (2 )

*x p xi*

*i i*

 

1 1

2

*i L*

1

 

( ) ( 1) (2 )

*x p x i*

*i i*

where *L* is a differential operator (possibly nonlinear), and *f* is an input vector.

*i uxyzt t xyz* 

*xyz* with time varying coefficients :

*i i*

(12)

( ) *<sup>i</sup> t*

(14)

*i i L u f L u f dt i n* (13)

*xyz* can be chosen arbitrarily, as long as their elements satisfy all

information of the polynomials degrees defined in other packages.

model while pull-in time errors is less than 2%.

(http:// www.memsscap.com)*.*

arbitral basis functions (,,) *<sup>i</sup>*

functions in the operating range *H*:

The basic functions (,,) *<sup>i</sup>*

*ψ*(*x*) satisfying

**3.1.4 Galerkin's approximated ROM** 

$$
\theta\_{l,k}(\mathbf{x}) = \theta(2^l \mathbf{x} - k) \,. \tag{16}
$$

Wavelets have proven to be an efficient tool of analysis in many fields including the solution of PDE. In [23] a new wavelet interpolation Galerkin method is used for the numerical simulation of MEMS devices under the effect of squeeze film damping. The air film pressure is expressed as a linear combination of a class of basic functions generated by autocorrelation of the usual compactly supported Daubechies scaling functions, which are the first- generation wavelets. The wavelet interpolation Galerkin method was used to predict the frequency response of the accelerometer with a spring mass-damper model with a parallel-plate electrostatic force. Various numerical tests have been conducted by changing the degree of the Daubechies wavelet *L* and the number *J* of the scale. Better accuracy can be achieved by increasing *L* and *J*. The higher *L* is, the smoother the scaling function becomes. The price for the high smoothness is that its supporting domain gets larger. The higher *J* is, the more accurate the solution becomes. The number of differential equations and the CPU time increase significantly as *J* increases. The solutions for *L* = 6 and *J* = 4 have results higher than the finite difference method.

The present wavelet interpolation Galerkin method is not suitable to solve problems defined on nonrectangular domains, since higher-dimensional wavelets are constructed by employing the tensor product of the one-dimensional wavelets, so their application is restricted to rectangular domains. But usage of the second-generation wavelets which are constructed in the spatial domain can expend in future PDE solutions to complex domains.

#### **3.2 Moment matching based ROM**

Moment matching model order reduction is based on the approximation of the original *n*dimensional system transfer function *F(s)* of the original *n*-dimensional system with a rational function with a lower degree *q<<n* [24]. This is done by matching some terms of the Taylor expansion of *F(s)* around a certain expansion point*.* For a state space equation

$$\begin{aligned} \dot{X}(t) &= AX(t) + Bv(t) \\ Y(t) &= \mathbf{C}^T X(t) \end{aligned} \tag{17}$$

let's perform a Laplace transform to obtain its frequency domain transfer function

$$F(\mathbf{s}) = \mathbf{C}^{T} \left(\mathbf{s}I - A\right)^{-1} B \tag{18}$$

Macromodels of Micro-Electro-Mechanical Systems (MEMS) 169

Note that the Arnoldi algorithm terminates when *hjj* = 0, which means that the subsequent vectors belong to the subspace already generated. The Arnoldi algorithm is basically a Gram–Schmidt procedure for orthogonalizing the Krylov vectors. The variant of the Arnoldi algorithm (also known as *PRIMA*) with some extra computational effort preserves the

The Moment matching method of model order reduction can be extended on nonlinear systems. In these cases the original nonlinear systems has to be changed previously (linearized or piecewise-linearized, approximated by a quadratic systems, divided into

A size of a second order equations system (2) can be also reduced by transforming it to the first order system (1), and then applying the methods described before. However, the reduction of second order systems by such transformation ignores the physical meaning of the original matrices and gives a reduced order model in a first order form. It is desirable for the reduced system to preserve the form of the original system (2). Approaches, that deal directly with the system (2) reduction have been proposed in the framework of Krylov

2 1 () ( ) *<sup>T</sup> H s C s M sD K B* (24)

 If the system is **undamped**, i.e. *D* = 0, the Arnoldi process can be applied for the computation of a basis for the Krylov subspace *Kq (K−1M,K−1B)* which is used for the projection matrix *V* building. The transfer function can be expanded into Taylor series as

> 0 ( ) *<sup>k</sup> k*

 ,

*k Hs ms* 

Then the reduced system matrices with variable change *x t Vx t* () () are obtained by the

Mq = VTMV, Kq = VTKV, Bq = VTB, Cq = CTV (25)

It is worth noting that, starting from a second order system in the form (2), Krylov subspace methods require the knowledge of *K−1*. For high dimensional systems, the explicit calculation of the inverse of *K* is computationally not affordable. Its computation is therefore

2

The transfer function for the system (2), with zero initial conditions, is given by:

iv. For j = 2, . . . , q:

1

*j j j i ij i w v vh* 

1

*h w v wh jj j j j jj* , / .

passivity of the original system [27].

several linear systems, etc.) [28].

**3.2.2 Second order systems** 

subspaces methods [29, 30].

where 1 1 ( )( ) *T k m C KM AB <sup>k</sup>*

orthogonal projection:

.

and the matching properties of the method are conserved [31].

 

Solve *<sup>j</sup> v* from: *Av Av j j* <sup>1</sup> . For *i* = 1*, . . . , j −* 1: *<sup>T</sup>*

*ij <sup>i</sup> <sup>j</sup> h vv* .

and expand it into Taylor series as

$$F(s) = -\sum\_{k=0}^{\infty} m\_k s^k \tag{19}$$

where <sup>1</sup> ( ) *T k m CA A B <sup>k</sup>* .

Moment matching is directly connected to the Krylov subspace formed by the pair of matrices (*A-1*, *A-1B*) [24]. The Krylov subspace is spanned by the column vectors in the following collection of matrices:

$$\{\mathbf{A}^{\text{-1}}\mathbf{B}, \{\mathbf{A}\}^{\text{-1}}\mathbf{A}^{\text{-1}}\mathbf{B}, \dots, \{\mathbf{A}\}^{\text{i}}\mathbf{A}^{\text{-1}}\mathbf{B}, \dots\}\tag{20}$$

where the column vectors are called the Krylov vectors. The *q-*th order Krylov subspace is denoted by

$$\mathbf{K}\_{\mathbf{q}}(\mathbf{A}^{\cdot 1}, \mathbf{A}^{\cdot 1}\mathbf{B})\tag{21}$$

which is spanned by the leading *q* linearly independent Krylov vectors in (20). Let *V є* R*n×q*  be any matrix which columns span the Krylov subspace *Kq*(*A-1, A*-*1B*)

$$\mathsf{Span}\,\mathsf{codim}\,\mathsf{w}\mathsf{V}\mathsf{J} = \mathsf{span}\{\mathsf{A}\cdot\mathsf{1}\mathsf{B};\mathsf{A}\cdot\mathsf{2}\mathsf{B};\ldots\mathsf{...}\mathsf{A}\mathsf{-}\mathsf{q}\mathsf{B}\}.$$

Here *V* is the orthogonal projection matrix that maps the *n*-dimensional state-space into a *q*  dimensional state-space and satisfies *VTV = I*. If the columns of *V* are orthogonal and *B* is a column vector, it can be shown that the following identities hold [24]:

$$(A) \forall A \cdot 1B = V(A\_q) \forall A\_q \cdot 1B\_q \tag{22}$$

for *i* = 0*,* 1*, . . . , q −* 1. These identities can be used to verify that at least the *q* leading moments of the full-order and reduced-order transfer functions are matched. Finally, we get the reduced order system of much smaller order (or state-space dimension) by performing variable change *x t Vx t* () () and multiplying on *VT* both sides of the equations (17):

$$\begin{aligned} \tilde{\mathbf{x}}'(t) &= A\_q \tilde{\mathbf{x}}(t) + B\_q \mathbf{v}(t) \\ \mathbf{y}(t) &= \mathbf{C}\_q^T \tilde{\mathbf{x}}(t) \end{aligned} \tag{23}$$
 
$$\text{where} \qquad \begin{aligned} \mathbf{A}\_q &= \{\mathbf{V} \mathbf{r} \mathbf{A} \mathbf{V} \}, \quad \mathbf{B}\_q = \mathbf{V} \mathbf{r} \mathbf{B}, \ \mathbf{C}\_q = \mathbf{C} \mathbf{r} \mathbf{V}. \end{aligned}$$

As a common practice, the block vectors forming the Krylov subspace are orthogonalized by using the Arnoldi algorithm for a numerical stability. Lets describe the block Arnoldi algorithm for a single column input matrix, i.e., B = b, where b є Rn, so that the algebraic operations involved can be seen.

#### **3.2.1 Arnoldi Algorithm**


$$\begin{aligned} \text{iv.} \quad \text{For } \mathbf{j} &= \mathbf{2}\_{\prime} \dots \mathbf{2}\_{\prime} \neq \\ \text{Solve } \tilde{\boldsymbol{\upsilon}}\_{j} \text{ from: } A\tilde{\boldsymbol{\upsilon}}\_{j} &= A\boldsymbol{\upsilon}\_{j-1} \dots \\ \text{For } \boldsymbol{i} &= \mathbf{1}\_{\prime} \dots \boldsymbol{i} \ \boldsymbol{j} - \mathbf{1} \colon \boldsymbol{h}\_{ij} = \boldsymbol{\upsilon}\_{i}^{T} \tilde{\boldsymbol{\upsilon}}\_{j} \ . \end{aligned}$$

$$\begin{aligned} \boldsymbol{\upsilon}\_{j} &= \tilde{\boldsymbol{\upsilon}}\_{j} - \sum\_{i=1}^{j-1} \boldsymbol{\upsilon}\_{i} \boldsymbol{h}\_{ij} \\ \boldsymbol{h}\_{jj} &= \left\| \boldsymbol{\upsilon}\_{j} \right\| \quad \boldsymbol{\upsilon}\_{j} = \boldsymbol{\upsilon}\_{j} / \left\| \boldsymbol{h}\_{jj} \right\| \ . \end{aligned}$$

0 ( ) *<sup>k</sup> k*

(19)


(23)

*k Fs ms* 

Moment matching is directly connected to the Krylov subspace formed by the pair of matrices (*A-1*, *A-1B*) [24]. The Krylov subspace is spanned by the column vectors in the

 {A-1B, (A)-1A-1B, . . . , (A)-i A-1B, . . .} (20) where the column vectors are called the Krylov vectors. The *q-*th order Krylov subspace is

 Kq(A-1, A-1B) (21) which is spanned by the leading *q* linearly independent Krylov vectors in (20). Let *V є* R*n×q* 

Spancolumn{V} =span{A-1B;A-2B; . . . ;A-qB}. Here *V* is the orthogonal projection matrix that maps the *n*-dimensional state-space into a *q*  dimensional state-space and satisfies *VTV = I*. If the columns of *V* are orthogonal and *B* is a

*A*-*1B*= *V*(*Aq*)*<sup>i</sup>*

variable change *x t Vx t* () () and multiplying on *VT* both sides of the equations (17):

() ()

*yt C xt*

 

where Aq = (VTAV ), Bq = VTB, Cq = CTV.

for *i* = 0*,* 1*, . . . , q −* 1. These identities can be used to verify that at least the *q* leading moments of the full-order and reduced-order transfer functions are matched. Finally, we get the reduced order system of much smaller order (or state-space dimension) by performing

() () ()

As a common practice, the block vectors forming the Krylov subspace are orthogonalized by using the Arnoldi algorithm for a numerical stability. Lets describe the block Arnoldi algorithm for a single column input matrix, i.e., B = b, where b є Rn, so that the algebraic

*x t Axt Bvt*

*q q T q*

*Aq*

be any matrix which columns span the Krylov subspace *Kq*(*A-1, A*-*1B*)

column vector, it can be shown that the following identities hold [24]:

and expand it into Taylor series as

.

(*A*)*<sup>i</sup>*

operations involved can be seen.

i. LU factorize matrix *A*: *A* = *LU*.

iii. Compute 11 1 *h v* and *v vh* 1 1 11 / .

**3.2.1 Arnoldi Algorithm** 

ii. Solve 1 *v* from: A <sup>1</sup> *v* = b.

following collection of matrices:

where <sup>1</sup> ( ) *T k m CA A B <sup>k</sup>*

denoted by

Note that the Arnoldi algorithm terminates when *hjj* = 0, which means that the subsequent vectors belong to the subspace already generated. The Arnoldi algorithm is basically a Gram–Schmidt procedure for orthogonalizing the Krylov vectors. The variant of the Arnoldi algorithm (also known as *PRIMA*) with some extra computational effort preserves the passivity of the original system [27].

The Moment matching method of model order reduction can be extended on nonlinear systems. In these cases the original nonlinear systems has to be changed previously (linearized or piecewise-linearized, approximated by a quadratic systems, divided into several linear systems, etc.) [28].

#### **3.2.2 Second order systems**

A size of a second order equations system (2) can be also reduced by transforming it to the first order system (1), and then applying the methods described before. However, the reduction of second order systems by such transformation ignores the physical meaning of the original matrices and gives a reduced order model in a first order form. It is desirable for the reduced system to preserve the form of the original system (2). Approaches, that deal directly with the system (2) reduction have been proposed in the framework of Krylov subspaces methods [29, 30].

The transfer function for the system (2), with zero initial conditions, is given by:

$$H(\mathbf{s}) = \mathbf{C}^{T} \left(\mathbf{s}^{2}M + \mathbf{s}D + \mathbf{K}\right)^{-1}B \tag{24}$$

 If the system is **undamped**, i.e. *D* = 0, the Arnoldi process can be applied for the computation of a basis for the Krylov subspace *Kq (K−1M,K−1B)* which is used for the projection matrix *V* building. The transfer function can be expanded into Taylor series as

$$H(\mathbf{s}) = -\sum\_{k=0}^{\infty} m\_k \mathbf{s}^{2k} \mathbf{s}^{2k}$$

where 1 1 ( )( ) *T k m C KM AB <sup>k</sup>* .

Then the reduced system matrices with variable change *x t Vx t* () () are obtained by the orthogonal projection:

$$\mathbf{M}\_{\mathbf{q}} = \mathbf{V}^{\mathrm{T}} \mathbf{M} \mathbf{V}, \ \mathbf{K}\_{\mathbf{q}} = \mathbf{V}^{\mathrm{T}} \mathbf{K} \mathbf{V}, \ \mathbf{B}\_{\mathbf{q}} = \mathbf{V}^{\mathrm{T}} \mathbf{B}, \ \mathbf{C}\_{\mathbf{q}} = \mathbf{C}^{\mathrm{T}} \mathbf{V} \tag{25}$$

and the matching properties of the method are conserved [31].

It is worth noting that, starting from a second order system in the form (2), Krylov subspace methods require the knowledge of *K−1*. For high dimensional systems, the explicit calculation of the inverse of *K* is computationally not affordable. Its computation is therefore

Macromodels of Micro-Electro-Mechanical Systems (MEMS) 171

considerable deviations at the high frequency range (fig.10,b). The model with order of 40 shows a perfect match for the lower eigenfrequencies and it is quite closer for higher

Using Krylov/Arnoldi approach, only a postprocessor to ANSYS is necessary to generate a macromodel in one of the well established model description languages: pure C code, HDL-A, MAST, Modelica and the new standardized VHDL-AMS which are supported by powerful system simulators. Such approach was implemented in *mor4ansys* (pronounced "more for ANSYS") that was developed by IMTEK [33] and which provides the reduced model of order 20 - 30 with an accuracy of a few percents when the dimension of an original

Taking in to account relations between displacements *x*, velocities *v* and accelerations *a*: , it is possible to present the equation (2) in the form [34,35]:

or , (26)

*Cv Gv Lv F t* ( )

where are equivalent matrices of capacitances, conductance and

*CGL* , ,

The elements of matrices , , *CGL* are formed from the elements of the mass, damping and

*C m ij N i j C m i N L k ij N i j*

, , 1(1) , ; , 1(1) ; 1/ , , 1(1) , ;

In this approach a capacitive-inductive-resistive model of the circuit is built which correctly reflects mass, damping and stiffness matrices. Nodal potentials in this model correspond to

Connected inductors, capacitors and conductors are placed in parallel between each two nodes and each node and ground (for the case i = j), fig.11. Values of these elements are defined by

Displacements *x* are defined by the current of inductances *Lii,* which are connected between a node and ground. Such approach suggests significant advantage if the coefficients of the

The task of reduction of MEMS model order turns now into reduction of the equivalent RLC circuit size. There are two ways of performing such circuit size decreasing with the minimal

equations (27). Similar approach was used for solving thermo-structural analysis [35].

1

1 1 1/ , 1(1) ; , , 1(1) , ; , 1(1) .

*N N ii ij ij ij ii ij j j L k i N G d ij N i j G d i N* 

*N*

*ij ij ii ij ij ij j*

(27)

frequencies, though this is not so important for the gyroscope.

( ) ( ) *<sup>d</sup> Mv Dv Kvdt F t*

FEM is up to 100 000.

*a dv dt x vdt* / ,

inductances.

**3.3 Equivalent-circuit ROM** 

stiffness matrices in the following ways:

*dt*

*C MG DL K* , ,

the displacement velocities *v*.

accuracy loss:

where N is a number of equations or nodes of the MEMS structure.

mass and stiffness matrices become time-dependant.

replaced by the solution of linear systems of equations through a *LU-* decomposition and defying *K−1*= *U−1 L−1*.

If the system is **damped**, i.e. *D* ≠ 0, in case of Rayleigh damping, it demonstrates that the damping matrix can be neglected during the reduction process and it can be computed afterwards, as a linear combination of the reduced stiffness and mass matrices [29]. That is why it is possible to recalculate the matrix *D* also similar to matrices *M* and *K* as in (25):

$$\mathbf{D}\_{\mathbf{q}} = \mathbf{V}^{\mathrm{T}} \mathbf{D} \mathbf{V}.$$

The reduction of a nonlinear system (2) can be done by using linear model order reduction techniques and considering nonlinearities as inputs. In the general case, the stiffness matrix is nonlinear since its entries dependent on the nodal displacements *x*. So damping matrix entries also are nonlinear functions of the nodal displacements *x* and the applied voltages *V*. If the nonlinearities are confined in the input function *f*, the system can be reduced using linear model order reduction technique. The only complication is that, after the reduction, the argument of *x* has to be recovered by the projection *x = Vxq* [31,32].

As a typical example Fig. 10 displays the transient simulation and frequency response of the original and reduced models of the microgyroscope being received via usage of the Arnoldi process [33].

Fig. 10. Comparison of transient behavior (a) and transfer functions (b) for the full and reduced microgyroscope models [33].

We see that the solution obtained by reduced model of order 10 is already very close to the true ANSYS excitation (fig.10,a) while the reduced models of order 15 up to 20 show considerable deviations at the high frequency range (fig.10,b). The model with order of 40 shows a perfect match for the lower eigenfrequencies and it is quite closer for higher frequencies, though this is not so important for the gyroscope.

Using Krylov/Arnoldi approach, only a postprocessor to ANSYS is necessary to generate a macromodel in one of the well established model description languages: pure C code, HDL-A, MAST, Modelica and the new standardized VHDL-AMS which are supported by powerful system simulators. Such approach was implemented in *mor4ansys* (pronounced "more for ANSYS") that was developed by IMTEK [33] and which provides the reduced model of order 20 - 30 with an accuracy of a few percents when the dimension of an original FEM is up to 100 000.

#### **3.3 Equivalent-circuit ROM**

170 Microelectromechanical Systems and Devices

replaced by the solution of linear systems of equations through a *LU-* decomposition and

If the system is **damped**, i.e. *D* ≠ 0, in case of Rayleigh damping, it demonstrates that the damping matrix can be neglected during the reduction process and it can be computed afterwards, as a linear combination of the reduced stiffness and mass matrices [29]. That is why it is possible to recalculate the matrix *D* also similar to matrices *M* and *K* as in (25):

Dq = VTDV. The reduction of a nonlinear system (2) can be done by using linear model order reduction techniques and considering nonlinearities as inputs. In the general case, the stiffness matrix is nonlinear since its entries dependent on the nodal displacements *x*. So damping matrix entries also are nonlinear functions of the nodal displacements *x* and the applied voltages *V*. If the nonlinearities are confined in the input function *f*, the system can be reduced using linear model order reduction technique. The only complication is that, after the reduction,

As a typical example Fig. 10 displays the transient simulation and frequency response of the original and reduced models of the microgyroscope being received via usage of the Arnoldi

(a)

(b)

We see that the solution obtained by reduced model of order 10 is already very close to the true ANSYS excitation (fig.10,a) while the reduced models of order 15 up to 20 show

Fig. 10. Comparison of transient behavior (a) and transfer functions (b) for the full and

reduced microgyroscope models [33].

the argument of *x* has to be recovered by the projection *x = Vxq* [31,32].

defying *K−1*= *U−1 L−1*.

process [33].

Taking in to account relations between displacements *x*, velocities *v* and accelerations *a*: , it is possible to present the equation (2) in the form [34,35]: *a dv dt x vdt* / ,

$$\frac{d}{dt}(M\mathbf{v}) + D\mathbf{v} + \int K\mathbf{v}dt = F(t) \quad \text{or} \quad \tilde{C}\dot{\mathbf{v}} + \tilde{G}\mathbf{v} + \tilde{L}\mathbf{v} = F(t) \tag{26}$$

where are equivalent matrices of capacitances, conductance and inductances. *C MG DL K* , ,

*CGL* , ,

The elements of matrices , , *CGL* are formed from the elements of the mass, damping and stiffness matrices in the following ways:

$$C\_{\vec{y}} = -m\_{\vec{y}}, \quad i, j = \text{l(l)} \, N, \quad i \neq j; \quad C\_{\vec{u}} = \sum\_{j=1}^{N} m\_{\vec{y}}, \quad i = \text{l(l)} \, N; \quad L\_{\vec{y}} = -1/\, k\_{\vec{y}}, \quad i, j = \text{l(l)} \, N, \quad i \neq j; \tag{27}$$

$$L\_{\vec{u}} = 1/\sum\_{j=1}^{N} k\_{\vec{y}}, \quad i = \text{l(l)} \, N; \quad G\_{\vec{y}} = -d\_{\vec{y}}, \quad i, j = \text{l(l)} \, N, \quad i \neq j; \quad G\_{\vec{u}} = \sum\_{j=1}^{N} d\_{\vec{y}}, i = \text{l(l)} \, N. \tag{27}$$

where N is a number of equations or nodes of the MEMS structure.

In this approach a capacitive-inductive-resistive model of the circuit is built which correctly reflects mass, damping and stiffness matrices. Nodal potentials in this model correspond to the displacement velocities *v*.

Connected inductors, capacitors and conductors are placed in parallel between each two nodes and each node and ground (for the case i = j), fig.11. Values of these elements are defined by equations (27). Similar approach was used for solving thermo-structural analysis [35].

Displacements *x* are defined by the current of inductances *Lii,* which are connected between a node and ground. Such approach suggests significant advantage if the coefficients of the mass and stiffness matrices become time-dependant.

The task of reduction of MEMS model order turns now into reduction of the equivalent RLC circuit size. There are two ways of performing such circuit size decreasing with the minimal accuracy loss:

Macromodels of Micro-Electro-Mechanical Systems (MEMS) 173

2

y2

Fig. 12. A working node of the RLC-circuit: conductance are added between node 1 and all

k

2 21 ( / ) /0 *k kk*

Note that it is equivalent to adding k*-1* new elements between first node and *k-1* former

For any two *a*-th and *b*-th nodes, which are neighbors to *i*- th node, elimination of *i*- th node

or in the p-polynomial form, taking into account existence of R, L and C elements, as shown

*ab i ab a a b b i i bb B y g pc g pc G pC pp p* 

introduced for any node in the circuit. The time constant of *i*-th node is defined as


і

In order to eliminate a fast node from *RLC* circuit, let us consider the following two cases.

 

are sums of all the capacitances, reciprocal

max*C B i i* /

*j jr*

*Y y y Y V y yV Y yV* 

*j i jj i rr*

11 11

will add a new element between these two nodes, which is equal to:

111 , , *kkk i j i j i j jjj C CB BG G* 

1. min max 2/ , *<sup>i</sup>*

inductances and conductance connected to *i*-th node. In order to simplify (33) two constant time values

and it is considered to be fast if

defined by user. So, a fast node will satisfy the following conditions:

*C GG B i ii i* , /

max max

  1

*r i*

(31)

1

y1 y2

yk

*yab a b i* ( )/ *yy Y* (32)

*RCi i i C G*/ and

max . (33)

*LCi i i C B*/ are

max being

2

neighbors nodes

on fig. 11:

where

max( , ) *i RC LC*

where min 

 

The equation (30) can be simplified:

k

yk

i

1

y1

neighbors of *i*-th node on Fig. 12.


Fig. 11. Circuit element for composing the equivalent-circuit model

#### **3.3.1 Y-∆ transformation based circuit size reduction methods**

The essence of the methods based on the Y-Δ transformation consists in the following. Let's *i*-th node and *k* its neighbors are located as shown on fig.11. Then the component equation of *i*-th row will look like

$$Y\_i V\_i - \mathbf{y}\_1 V\_1 - \mathbf{y}\_2 V\_2 - \dots - \mathbf{y}\_n V\_n = \mathbf{0} \tag{28}$$

where Let's define *Vi* as 1 . *k i j j Y y* 

$$\boldsymbol{V}\_{i} = \left(\sum\_{j=1}^{k} \mathbf{y}\_{j} \mathbf{V}\_{j}\right) \Bigg/ \mathbf{Y}\_{i} \tag{29}$$

and replace *Vi* in other *k* equations , that is equivalent to excluding *i*- th node from a circuit. Then, for example, the equation of the first node transfers to:

$$\left( (\overline{Y}\_1 + y\_1 - y\_1^2 \; / \; Y\_i) V\_1 - \left( \sum\_{j=2}^k y\_1 y\_j V\_j \right) \right) \prime Y\_i \quad - \sum\_{r=1 \atop r \neq i}^{k1} y\_r V\_r = 0 \tag{30}$$

where 1 1 1 *k r r r i Y y* is a sum of all node conductance excluding *i*- th node, *k1* is a number of

nodes being connected to node 1.

Fig. 12. A working node of the RLC-circuit: conductance are added between node 1 and all neighbors nodes

The equation (30) can be simplified:

$$(\overline{Y}\_1 + \left(\sum\_{j=2}^k y\_1 y\_j\right) / Y\_i) V\_1 - \left(\sum\_{j=2}^k y\_1 y\_j V\_j\right) / Y\_i - \sum\_{\substack{r=1 \\ r \neq i}}^{k1} y\_r V\_r = 0 \tag{31}$$

Note that it is equivalent to adding k*-1* new elements between first node and *k-1* former neighbors of *i*-th node on Fig. 12.

For any two *a*-th and *b*-th nodes, which are neighbors to *i*- th node, elimination of *i*- th node will add a new element between these two nodes, which is equal to:

$$(y\_{ab} = (y\_a y\_b) / \,\text{Y}\_i\tag{52}$$

or in the p-polynomial form, taking into account existence of R, L and C elements, as shown on fig. 11:

$$\mathbf{y}\_{ab} = \left(\mathbf{g}\_a + \frac{\mathbf{b}\_a}{p} + p\mathbf{c}\_a\right) \left(\mathbf{g}\_b + \frac{\mathbf{b}\_b}{p} + p\mathbf{c}\_b\right) \bigg/ \left(\mathbf{G}\_i + \frac{\mathbf{B}\_i}{p} + p\mathbf{C}\_i\right) \tag{33}$$

172 Microelectromechanical Systems and Devices


The essence of the methods based on the Y-Δ transformation consists in the following. Let's *i*-th node and *k* its neighbors are located as shown on fig.11. Then the component equation

, (28)

11 2 2 ... 0 *YV yV y V y V i i n n*

and replace *Vi* in other *k* equations , that is equivalent to excluding *i*- th node from a circuit.

( /) / 0

 

*Y y y Y V y yV Y yV*

(29)

1

(30)

*r i*

2 1

is a sum of all node conductance excluding *i*- th node, *k1* is a number of

*k k i jj i rr j r*

1 *k i jj i j V yV Y* 

transformation [36, 39 ];

of *i*-th row will look like

1 . *k i j j Y y* 

where Let's define *Vi* as

where

1

*k r r r i Y y* 

1

nodes being connected to node 1.

1


Fig. 11. Circuit element for composing the equivalent-circuit model

**3.3.1 Y-∆ transformation based circuit size reduction methods** 

Then, for example, the equation of the first node transfers to:

2 111 1 1 where 111 , , *kkk i j i j i j jjj C CB BG G* are sums of all the capacitances, reciprocal

inductances and conductance connected to *i*-th node.

In order to simplify (33) two constant time values *RCi i i C G*/ and *LCi i i C B*/ are introduced for any node in the circuit. The time constant of *i*-th node is defined as max( , ) *i RC LC* and it is considered to be fast if

$$1. \tag{1.1}$$

$$\tau\_i < \tau\_{\min} = 2\pi \not\subset \phi\_{\max}.$$

where min - a time constant which depends on maximal circuit frequencymax being defined by user. So, a fast node will satisfy the following conditions:

 $\left| \left| \boldsymbol{\rho}\_{\max} \mathbf{C}\_{i} < \mathbf{G}\_{i} \right. \right. \left. \mathbf{G}\_{i} < \mathbf{B}\_{i} \right. \right. \left. \left. \left. \left. \boldsymbol{\rho}\_{\max} \right. \right. \right. \right. \left. \left. \left. \left. \boldsymbol{\rho}\_{\max} \right\} \right. \right. \right. \left. \left. \left. \left. \boldsymbol{\rho}\_{\max} \right. \right. \right. \right. \right. \left. \left. \left. \left. \left. \left. \boldsymbol{\rho}\_{\max} \right\} \right. \right. \right. \right. \right. \left. \left. \left. \left. \left. \left. \boldsymbol{\rho}\_{\max} \right\} \right. \right. \right. \right. \right. \left. \left. \left. \left. \left. \left. \boldsymbol{\rho}\_{\max} \right\} \right. \right. \right. \right. \right. \right. \left. \left. \left. \left. \left. \left. \left. \boldsymbol{\rho}\_{\max} \right\} \right. \right. \right. \right. \right. \right. \right. \left. \right. \right. \right. \left. \right. \right. \left. \right. \right. \left. \right. \right. \left. \right. \right. \left. \right. \right. \left. \right. \right. \left. \right. \right. \left. \right. \left. \right. \right. \left. \right. \right. \left. \right. \right. \left. \right. \left. \right. \right. \left. \right. \left. \right. \right. \left. \right. \right. \left. \right. \left. \right. \right. \left. \right. \left. \right. \right. \left. \right. \left. \right. \right. \left. \right. \right. \left. \right. \right$ 

In order to eliminate a fast node from *RLC* circuit, let us consider the following two cases.

Macromodels of Micro-Electro-Mechanical Systems (MEMS) 175

Table 1. Transformations for *LCi RCi*

 case

Original branch Substitution Formula

*bab=babb/Bi* 

*bab = babb/Bi Cab = baCb/Bi*

*bab=babb/Bi*

*Cab = baCb/Bi*

*Cab=ba\*Cb /Bi*

*Cab=CaCb/C* 

*Cab=Cagb/Gi*

*gab=gagb/Gi*

*gab=bagb/Gi*

*gab=bagb/Gi*

*Cab=(baCb++bbCa)/Bi*

2. If *RCi LCi* , the equation (33) can be transformed into:

$$y\_{ab} = (\mathcal{g}\_a + p\mathbf{c}\_a)(\mathcal{g}\_b + p\mathbf{c}\_b) / (\mathbf{G}\_l + p\mathbf{C}\_i) \tag{34}$$

and its decomposition into Taylor series will look like:

$$\mathbf{y}\_{ab} = \frac{\mathbf{g}\_a \mathbf{g}\_b}{G\_i} + p \frac{\mathbf{c}\_a \mathbf{g}\_b + \mathbf{c}\_b \mathbf{g}\_a}{G\_i} + p^2 \frac{4 \mathbf{c}\_a \mathbf{c}\_b}{G\_i} \,. \tag{35}$$

For typical *R,L,C* values (*R*=0,1÷1000 ohm, *L*=0,01÷10 nH, *C*=0,001÷10 pF, ω=0,1÷10 GHz) the contribution of the last term in (35) will be significantly smaller that previous ones. The constant term in (35) gives the value of the conductance, which appears between *a*- th and *b*th nodes during the elimination of *i*- th node, and the second term in (35) defines the value of capacitance.

3. If *LCi RCi* , the equation (33) is transformed into:

$$\mathbf{y}\_{ab} = \frac{1}{p} \frac{\mathbf{b}\_a \mathbf{b}\_b}{B\_i} + p \frac{\mathbf{c}\_a \mathbf{b}\_b + \mathbf{c}\_a \mathbf{b}\_b}{B\_i} \tag{36}$$

The first term with in (36) defines the value of the reactive conductance (the inductances' reciprocal value), which appears between *a*- th and *b*- th nodes during the elimination of *i*- th node, and the second term in (36) gives the value of the capacitance. 1 *p*

Final formulas for additional elements between *a*-th and *b*-th nodes during the elimination of *i*- th node for all possible cases are given in tables 1 and 2. All these formulas can be deducted from (35) and (36). The only exception is the case when *a*- th and *b*- th nodes are connected to *i*- th node through a capacitor when *C CC C ab a b i* / .

In order to eliminate a node from equivalent MEMS circuit with the smallest inaccuracy, the following two criteria should be taken into account:


In practice, usually lower eigenfrequencies for mechanical systems are most interesting. Therefore, a compromise between accuracy and size of the received circuit models can be reached by proper selection of τ min value. So the algorithm of equivalent-circuit ROM can be described as following:


For typical *R,L,C* values (*R*=0,1÷1000 ohm, *L*=0,01÷10 nH, *C*=0,001÷10 pF, ω=0,1÷10 GHz) the contribution of the last term in (35) will be significantly smaller that previous ones. The constant term in (35) gives the value of the conductance, which appears between *a*- th and *b*th nodes during the elimination of *i*- th node, and the second term in (35) defines the value

*<sup>g</sup> g cg cg cc yp p GG G* 

<sup>2</sup> 4 *a b a b b a ab*

*ii i*

 . (36) The first term with in (36) defines the value of the reactive conductance (the inductances' reciprocal value), which appears between *a*- th and *b*- th nodes during the elimination of *i*- th node, and the second term in (36) gives the value of the capacitance. Final formulas for additional elements between *a*-th and *b*-th nodes during the elimination of *i*- th node for all possible cases are given in tables 1 and 2. All these formulas can be deducted from (35) and (36). The only exception is the case when *a*- th and *b*- th nodes are

*y p pB B*

1 *ab ab ab*

*bb cb cb*

*i i*

In order to eliminate a node from equivalent MEMS circuit with the smallest inaccuracy, the

In practice, usually lower eigenfrequencies for mechanical systems are most interesting. Therefore, a compromise between accuracy and size of the received circuit models can be reached by proper selection of τ min value. So the algorithm of equivalent-circuit ROM can be

2. All the nodes should be put into a priority queue, sorted by time constants, apart from

5. Eliminated *i*- th node according to the rules described in tables 1, 2 (depending on

6. Time constants in the set *k* (which consists of the neighbor of i-th node) should be

3. Take the first i- th node in the queue (the one with the smallest time constant *<sup>i</sup>* min

should not be compared (have almost the same values).

( )( ) /( ) *ab a a b b I i y gp c g pc G pC* (34)

. (35)

 ).

, the equation (33) can be transformed into:

, the equation (33) is transformed into:

*ab*

connected to *i*- th node through a capacitor when *C CC C ab a b i* / .

1. Time constants should be calculated for all the nodes of the circuit.

7. Take i- th node from the head of the priority queue and jump to step 3.

following two criteria should be taken into account:

1 *p*

input, output and ground nodes.

4. Find the neighbors *k* of *i*- th node.

values).


described as following:

, *LCi RCi* 

recalculated

and its decomposition into Taylor series will look like:

*ab*

2. If *RCi LCi* 

of capacitance. 3. If *LCi RCi* 


Table 1. Transformations for *LCi RCi* case

Macromodels of Micro-Electro-Mechanical Systems (MEMS) 177

(37)

*bbba abaa*

*YY YY*

 

where *Ia, Ua* are current and voltage at the macromodel input, *Ib, Ub* are current and voltage at the output. It is suggested to obtain equation (37) directly from the general matrix of the

1 *bb ba*

 

, 

*b a*

*U U*

(38)

*aa bb ab aa*

,

where *Δij –* an algebraic complement to *aij* element of the initial circuit admittance matrix Y, which is equal to a determinant of the matrix being obtained from the matrix *Y* after eliminating *i* – th row and *j*- th column on the crossing of which the given element *aij* is

*Δaa,bb –* redoubled algebraic complement which equals to the determinant of the matrix obtained by eliminating *a*-th,*b*-th rows and *a*-th ,*b*-th columns and its sign is defined by the

The necessary algebraic complements can be calculated by the initial matrix inversion procedure, if matrix elements are defined for the selected frequency ω0*,* which is fixed by an

> 11 21 1 12 22 2

.... ...

*n n*

,

.

(39)

(40)

*aa bb*

...

*n n nn*

1 2

11 21 1 1 12 22 2

*gg g*

.... ...

*n n*

...

*n n nn*

*gg g*

,

/ ;

*aa bb aa bb*

*Y Y gg*

*Y gg*

*Y gg*

*bb aa aa bb*

Since macromodel parameters are determined by sum of the real and imaginary parts:

,

/ .

*ab ba ba aa bb*

,

/ ;

1 2

it is possible to find parameters of the equivalent four- terminal (38). By choosing elements ,,, *aa bb ab ba gggg* of the inverted matrix (39) and computing *aa bb aa bb ab ba* , *g gg gg* it is

*gg g <sup>Y</sup>*

circuit (e.g. admittance matrix Y) according to the expression [40]:

situated. In addition the sign of *Δij* is defined by the factor ( - 1 ) *i + j*;

Y-1 =

where ∆ *-* is the determinant of the initial matrix Y and , *bb aa ba ab*

possible to find parameters of the reduced model (37) from the relations:

So if numerical values of the inverted Y-1 matrix elements are computed:

1 

factor (-1) *a + a + b +b*.

user:

*tp*

*b a*

*I I*

*Y*


Table 2. Transformations for *RCi LCi* case

#### **3.3.2 Building MEMS macromodel as a four- terminal**

The main idea of this method is to develop a MEMS macromodel as a four-terminal circuit (n-terminal in a general case) as:

176 Microelectromechanical Systems and Devices

Table 2. Transformations for *RCi LCi*

(n-terminal in a general case) as:

 case

**3.3.2 Building MEMS macromodel as a four- terminal** 

The main idea of this method is to develop a MEMS macromodel as a four-terminal circuit

Original branch Substitution Formula

*gab=gagb/Gi Cab=(gaCb+bCa)/Gi*

*gab=gagb/Gi Cab=gaCb/Gi*

*Cab=gbCa/Gi* 

*Cab=gbCa/Gi* 

*Cab=CaCb/C* 

*gab=gagb/Gi* 

*gab=bagb /Bi*

*gab=bagb /Bi*

*gab=baCb /Bi*

*gab=babb /Bi*

$$
\begin{bmatrix} I\_a \\ I\_b \end{bmatrix} = \begin{bmatrix} Y\_{aa} & Y\_{ab} \\ Y\_{ba} & Y\_{bb} \end{bmatrix} \begin{bmatrix} U\_a \\ U\_b \end{bmatrix}, \tag{37}
$$

where *Ia, Ua* are current and voltage at the macromodel input, *Ib, Ub* are current and voltage at the output. It is suggested to obtain equation (37) directly from the general matrix of the circuit (e.g. admittance matrix Y) according to the expression [40]:

$$Y\_{tp} = \frac{1}{\Delta\_{aa,bb}} \begin{bmatrix} \Delta\_{bb} & & -\Delta\_{bu} \\\\ \Delta\_{ab} & & \Delta\_{aa} \end{bmatrix} \tag{38}$$

where *Δij –* an algebraic complement to *aij* element of the initial circuit admittance matrix Y, which is equal to a determinant of the matrix being obtained from the matrix *Y* after eliminating *i* – th row and *j*- th column on the crossing of which the given element *aij* is situated. In addition the sign of *Δij* is defined by the factor ( - 1 ) *i + j*;

*Δaa,bb –* redoubled algebraic complement which equals to the determinant of the matrix obtained by eliminating *a*-th,*b*-th rows and *a*-th ,*b*-th columns and its sign is defined by the factor (-1) *a + a + b +b*.

The necessary algebraic complements can be calculated by the initial matrix inversion procedure, if matrix elements are defined for the selected frequency ω0*,* which is fixed by an user:

$$\mathbf{Y}^{\cdot 1} = \left. \frac{1}{\Delta} \begin{bmatrix} \Delta\_{11} \,\Delta\_{21} \dots \,\Delta\_{n1} \\\\ \Delta\_{12} \,\Delta\_{22} \dots \,\Delta\_{n2} \\\\ \Delta\_{1n} \,\Delta\_{2n} \dots \,\Delta\_{nn} \end{bmatrix} \right| $$

where ∆ *-* is the determinant of the initial matrix Y and , *bb aa ba ab aa bb* . So if numerical values of the inverted Y-1 matrix elements are computed:

$$Y^{-1} = \begin{bmatrix} \mathcal{S}\_{11} \ \mathcal{S}\_{21} \cdots \ \mathcal{S}\_{n1} \\\\ \mathcal{S}\_{12} \ \mathcal{S}\_{22} \cdots \ \mathcal{S}\_{n2} \\\\ \mathcal{S}\_{1n} \ \mathcal{S}\_{2n} \cdots \ \mathcal{S}\_{nn} \end{bmatrix} \tag{39}$$

it is possible to find parameters of the equivalent four- terminal (38). By choosing elements ,,, *aa bb ab ba gggg* of the inverted matrix (39) and computing *aa bb aa bb ab ba* , *g gg gg* it is possible to find parameters of the reduced model (37) from the relations:

$$\begin{aligned} Y\_{aa} &= \mathcal{g}\_{bb} / \mathcal{g}\_{aa,bb}; \\ Y\_{ab} &= -Y\_{ba} = \mathcal{g}\_{hu} / \mathcal{g}\_{aa,bb}; \\ Y\_{bb} &= \mathcal{g}\_{aa} / \mathcal{g}\_{aa,bb}. \end{aligned} \tag{40}$$

Since macromodel parameters are determined by sum of the real and imaginary parts:

Macromodels of Micro-Electro-Mechanical Systems (MEMS) 179

Then to combine individual four-terminal circuits into one equivalent MEMS four-terminal circuit taking in account ways of connectivity of different subcircuits ( parallel, sequence or mixed) and proper recalculating four-terminal circuits parameters systems ( *y, z, h, a*) [40].

Equivalent-circuit ROM for MEMS was implemented in the circuit simulation package NetALLTED (**ALL TE**chnologies **D**esinger) which was developed not only for simulation and analysis, but for processing project procedures such as parametric optimization tasks; optimal tolerance assignments; centering availability regions; yield maximization [42]. NetALLTED is widely used for design of Nonlinear Dynamic Systems composed of either/and electronic, hydraulic, pneumatic, mechanical, electromagnetic, and other elements and it is available through the Internet (http://allted.kpi.ua/). ROM developing approach in hand provides more than 99% reduction of elements and node numbers of equivalent- circuit ROM. For example, for the accelerometer only 3 nodes and 6 elements are

Let's consider the example of the equivalent- circuit ROM for the beam working on bending and find its eigenfrequencies [34]. The left end of the beam is fixed motionlessly, right one is

The initial FEM model was constructed using ANSYS Multiphysics v.10.0 and the beam ANSYS library's BEAM3 finite element with the length of 0.5 um for following beam parameters: *L* = 25 um; beam cross-section is a square one with height of 3 um and width of 2 um. Beam material properties: coefficient of elasticity E=2·1011, Pa = 0.2 N/mkm2; Poisson coefficient μ = 0.3; material density ρ = 6·103 kg/m3 = 6·10-9 mg/mkm3. The initial

The developed equivalent-circuit ROM with 5 nodes and 14 elements is presented on fig.14.

free. Force *f* is applied to the right end perpendicularly to the beam axe.

equivalent beam circuit contains 101 nodes and 314 elements.

**3.3.3 NetALLTED equivalent-circuit ROM subsystem** 

left from initial 1,883 nodes and 62,826 elements.

Object Circuit Beam; J1(100,0)=-100;

0.116667;

&&

C\_1(82,100) = -

L\_10(100,0) = -34161.8; L\_11(23,50) = 1.35e-10; L\_12(23,0) = 1.15e-10; L\_13(50,0) = -16435.9; L\_14(50,82) = 1.6e-10; C\_2(100,0) = 6.3; C\_3(23,50) = -0.116667; C\_4(23,0) = 17.3833; C\_5(50,82) = -0.116667; C\_6(0,50) = 20.65; C\_7(0,82) = 17.5; L\_8(82,100) = 9e-11; L\_9(82,0) = -580750;

Fig. 14. Equivalent beam circuit

$$\mathbf{Y}\_{\overline{\eta}} \equiv \mathbf{a}\_0 \mathbf{+} \mathbf{i} \text{ a} \mathbf{1} \tag{41}$$

it is convenient to represent the circuit macromodel also as a sum of real and imaginary terms:

$$\boldsymbol{Y}\_{tp} = \begin{bmatrix} \boldsymbol{Y}\_{aa0}\boldsymbol{Y}\_{ab0} \\ \boldsymbol{Y}\_{ba0}\boldsymbol{Y}\_{bb0} \end{bmatrix} + \mathrm{i} \begin{bmatrix} \boldsymbol{Y}\_{aai}\boldsymbol{Y}\_{abi} \\ \boldsymbol{Y}\_{bai}\boldsymbol{Y}\_{bbi} \end{bmatrix} . \tag{42}$$

Investigation the expression (42) for different frequencies confirms that for linear circuits at least for frequency range <sup>0</sup> *,* where ω0 *–* frequency at which elements of the inverted matrix (40) are initially defined, real terms of the computed parameters (42) preserve their values and imaginary terms change proportionally to the selected frequency. This confirms the possibility to use macromodel parameters being computed at once frequency in a broad frequency range.

As it was noted earlier, the macromodel form (42) is inconvenient for direct implementation in the circuit simulation software. A real part of Yij parameters for the stable passive macromodels (a0>0) has to be positive while an imaginary part could be negative. If in (41) 0, *<sup>o</sup> a* <sup>1</sup> *a* 0 the Yij component of the four- terminal may be presented by parallel connection of conductance, capacitance and inductance [41]:

$$\mathbf{G}\_{\overline{\eta}} = \text{Re } (\mathbf{Y}\_{i,\overline{\eta}}),\\\mathbf{C}\_{\overline{\eta}} = \text{k Im } (\mathbf{Y}\_{i\overline{\eta}}) / \ \alpha\_{0\nu} \text{ L}\_{\overline{\eta}} = 1 / [(\text{k-1})\text{Im } (\mathbf{Y}\_{i\overline{\eta}}) \ \alpha\_0]. \tag{43}$$

If 0, *<sup>o</sup> a* <sup>1</sup> *a* 0 the respective component *Yij* may be presented in the same way as parallel connection of conductance, capacitance and inductance but:

$$\mathbf{G}\_{\overline{\mathbf{i}}} = \text{Re } (\mathbf{Y}\_{\overline{\mathbf{i}}\overline{\mathbf{j}}}) \text{ , } \mathbf{C}\_{\overline{\mathbf{i}}} = \text{(k-1) Im } (\mathbf{Y}\_{\overline{\mathbf{i}}}) \text{ } \alpha \mathbf{o} \text{ , and } \quad \mathbf{L}\_{\overline{\mathbf{i}}} = \mathbf{1} / (\mathbf{k} \text{ Im } (\mathbf{Y}\_{\overline{\mathbf{i}}}) \text{ } \alpha \mathbf{o}) . \tag{44}$$

where *k* is defined by reactive components values ratio.

So it is possible to select priori the equivalent-circuit macromodel, for example, shown on Fig.13, and to define parameters of its components using expressions (43) and (44).

Fig. 13. Equivalent-circuit macromodel of a four-terminal

The important feature of the equivalent circuits MEMS macromodel is a possibility to use optimization procedures for accurate adjustment of the macromodel components values to meet the device characteristics being obtained by FEM model. If a whole MEMS equivalent circuit is rather large it is possible to divide it into some subcircuits and apply transformation (37) to each subcircuit.

 Yij= a0+i a1 , (41) it is convenient to represent the circuit macromodel also as a sum of real and imaginary

0 0

Investigation the expression (42) for different frequencies confirms that for linear circuits at

matrix (40) are initially defined, real terms of the computed parameters (42) preserve their values and imaginary terms change proportionally to the selected frequency. This confirms the possibility to use macromodel parameters being computed at once frequency in a broad

As it was noted earlier, the macromodel form (42) is inconvenient for direct implementation in the circuit simulation software. A real part of Yij parameters for the stable passive macromodels (a0>0) has to be positive while an imaginary part could be negative. If in (41) 0, *<sup>o</sup> a* <sup>1</sup> *a* 0 the Yij component of the four- terminal may be presented by parallel

 Gij = Re (Yi,j), Cij = k Im (Yij)/ ω0, Lij = 1/[(k-1)Im (Yij) ω0]. (43) If 0, *<sup>o</sup> a* <sup>1</sup> *a* 0 the respective component *Yij* may be presented in the same way as parallel

Gij = Re (Yi,j) , Cij = (k-1) Im (Yij)/ ω0 and Lij = 1/(k Im (Yij) ω0), (44)

So it is possible to select priori the equivalent-circuit macromodel, for example, shown on

The important feature of the equivalent circuits MEMS macromodel is a possibility to use optimization procedures for accurate adjustment of the macromodel components values to meet the device characteristics being obtained by FEM model. If a whole MEMS equivalent circuit is rather large it is possible to divide it into some subcircuits and apply

Fig.13, and to define parameters of its components using expressions (43) and (44).

*Y i*

*aa ab aai abi*

. (42)

*ba bb bai bbi Y Y YY*

<sup>0</sup> *,* where ω0 *–* frequency at which elements of the inverted

*Y Y YY* 

0 0

*tp*

 

connection of conductance, capacitance and inductance [41]:

connection of conductance, capacitance and inductance but:

where *k* is defined by reactive components values ratio.

Fig. 13. Equivalent-circuit macromodel of a four-terminal

transformation (37) to each subcircuit.

terms:

least for frequency range

frequency range.

Then to combine individual four-terminal circuits into one equivalent MEMS four-terminal circuit taking in account ways of connectivity of different subcircuits ( parallel, sequence or mixed) and proper recalculating four-terminal circuits parameters systems ( *y, z, h, a*) [40].

#### **3.3.3 NetALLTED equivalent-circuit ROM subsystem**

Equivalent-circuit ROM for MEMS was implemented in the circuit simulation package NetALLTED (**ALL TE**chnologies **D**esinger) which was developed not only for simulation and analysis, but for processing project procedures such as parametric optimization tasks; optimal tolerance assignments; centering availability regions; yield maximization [42]. NetALLTED is widely used for design of Nonlinear Dynamic Systems composed of either/and electronic, hydraulic, pneumatic, mechanical, electromagnetic, and other elements and it is available through the Internet (http://allted.kpi.ua/). ROM developing approach in hand provides more than 99% reduction of elements and node numbers of equivalent- circuit ROM. For example, for the accelerometer only 3 nodes and 6 elements are left from initial 1,883 nodes and 62,826 elements.

Let's consider the example of the equivalent- circuit ROM for the beam working on bending and find its eigenfrequencies [34]. The left end of the beam is fixed motionlessly, right one is free. Force *f* is applied to the right end perpendicularly to the beam axe.

The initial FEM model was constructed using ANSYS Multiphysics v.10.0 and the beam ANSYS library's BEAM3 finite element with the length of 0.5 um for following beam parameters: *L* = 25 um; beam cross-section is a square one with height of 3 um and width of 2 um. Beam material properties: coefficient of elasticity E=2·1011, Pa = 0.2 N/mkm2; Poisson coefficient μ = 0.3; material density ρ = 6·103 kg/m3 = 6·10-9 mg/mkm3. The initial equivalent beam circuit contains 101 nodes and 314 elements.

The developed equivalent-circuit ROM with 5 nodes and 14 elements is presented on fig.14.

Fig. 14. Equivalent beam circuit

Macromodels of Micro-Electro-Mechanical Systems (MEMS) 181

The advantage of the equivalent circuit approach is obtaining small size of the equivalent reduced circuit as well as the possibility to get required frequencies with a high accuracy, using NetALLTED optimization possibilities (fig.16). The disadvantage is a necessity to make the additional analysis of reduced circuit in order to find the most sensitive elements

For example, for reduced circuit, shown at fig.14, four variable elements have been chosen (*L8, L11, L12, L14*) and the objective function *OF ERROR1= F8(1336.3,4009.3/T1,T2)* was constructed which contains the requirement to obtain resonance peaks at frequencies *T1=1336.3* Hz and *T2=4009.3* Hz in according to the ANSYS analysis results (table 3). The objective function contains also current values of frequency response extremes *T1=MAXA(db.K1,100,1600)* and *T2=MAXA(db.K1,1700,4100),* which are calculated with help of MAXA directive for defining time or frequency when the output (*db.K1* in our case) reaches its maximum value in the specified time or frequency range (*100-1600* Hz). Among available 12 optimization methods being incorporated in NetALLED the Random Search Method (METHOD=40) with search interval reducing has been used to optimize the beam

The modeling of MEMS provides a very challenging task in modern engineering. This field of research is inherently multiphysics of nature, since different physical phenomena are

**Node number** 101 3 Element number 314 6 1st peak, Hz 1336.3 1337.0(0.05%) 2nd peak, Hz 4009.3 4008.9(0.01%) Table 4. Results of frequency analysis of the four- terminal equivalent circuit for a beam

and take them as variable parameters during optimization procedure.

Fig. 16. ROM transfer function before (1) and after (2) optimization.

macromodel parameters.

**4. MEMS coupled system-level model** 

Source circuit

Reduced circuit (as n-ports)

The two lower eigenfrequencies of the beam are defined by computation. The results of ANSYS Multiphysics frequency analysis as well as the results of the equivalent- circuit ROM simulation for different τmin by NetALLTED [36,43] are given in table 3.

It is obvious that the accuracy of a macromodel obtained with τ min =3\*10-5 is rather high and there are only 5 nodes and 14 elements in the reduced circuit. For more accurate simulation it is possible either to use a reduced circuit obtained with smaller values of τ min (when a size of equivalent circuit ROM increases), or to adapt the received ROM with help of the optimization methods.


Table 3. Frequency beam analysis results

The macromodel of the four-terminal for the same beam is being developed in according with (38)-(44) with only 3 nodes and 6 elements that presented on fig.15.

Fig. 15. Equivalent beam circuit being considered as a four- terminal

The results of simulations of the four- terminal equivalent circuit by NetALLTED [36,43] are given in table 4.

The two lower eigenfrequencies of the beam are defined by computation. The results of ANSYS Multiphysics frequency analysis as well as the results of the equivalent- circuit ROM

It is obvious that the accuracy of a macromodel obtained with τ min =3\*10-5 is rather high and there are only 5 nodes and 14 elements in the reduced circuit. For more accurate simulation it is possible either to use a reduced circuit obtained with smaller values of τ min (when a size of equivalent circuit ROM increases), or to adapt the received ROM with help of the

Source

 , s - - 5\*10-6 10-5 3\*10-5 3\*10-5 **Node number** - 101 24 12 5 5 Element number - 314 76 38 14 14 Reduction by nodes, % - - 76,2376 88,1188 95,0495 95,0495 Reduction by elements, % - - 75,7962 87,8981 95,5414 95,5414 1st peak, Hz 1336,2 1336,3 1336,1 1334,9 1327 1336,2 2nd peak, Hz 4009,3 4009,3 4009,4 3993,3 3612,1 4009,4 Maximal error, % - - 0,01 0,3 9,9 0,003

The macromodel of the four-terminal for the same beam is being developed in according

The results of simulations of the four- terminal equivalent circuit by NetALLTED [36,43] are

with (38)-(44) with only 3 nodes and 6 elements that presented on fig.15.

Fig. 15. Equivalent beam circuit being considered as a four- terminal

ALLTED results

circuit Reduced circuit Optim.

circuit

simulation for different τmin by NetALLTED [36,43] are given in table 3.

ANSYS results

optimization methods.

Table 3. Frequency beam analysis results

L2(1,2) = 2.997781327698e-011; C2(1,2) = 3.155049077722e+002; L1(1,0) = 1.125969906329e-009; C1(1,0) = 1.890000000000e+001; L3(2,0) = 7.750537495099e-010; C3(2,0) = 1.220321457595e+001;

min 

Object Circuit Beam; J1(1,0)=-100;

&&

given in table 4.


Table 4. Results of frequency analysis of the four- terminal equivalent circuit for a beam

The advantage of the equivalent circuit approach is obtaining small size of the equivalent reduced circuit as well as the possibility to get required frequencies with a high accuracy, using NetALLTED optimization possibilities (fig.16). The disadvantage is a necessity to make the additional analysis of reduced circuit in order to find the most sensitive elements and take them as variable parameters during optimization procedure.

Fig. 16. ROM transfer function before (1) and after (2) optimization.

For example, for reduced circuit, shown at fig.14, four variable elements have been chosen (*L8, L11, L12, L14*) and the objective function *OF ERROR1= F8(1336.3,4009.3/T1,T2)* was constructed which contains the requirement to obtain resonance peaks at frequencies *T1=1336.3* Hz and *T2=4009.3* Hz in according to the ANSYS analysis results (table 3). The objective function contains also current values of frequency response extremes *T1=MAXA(db.K1,100,1600)* and *T2=MAXA(db.K1,1700,4100),* which are calculated with help of MAXA directive for defining time or frequency when the output (*db.K1* in our case) reaches its maximum value in the specified time or frequency range (*100-1600* Hz). Among available 12 optimization methods being incorporated in NetALLED the Random Search Method (METHOD=40) with search interval reducing has been used to optimize the beam macromodel parameters.
