Macromodels of Micro-Electro-Mechanical Systems ( М EMS )

Micro-electro-mechanical Systems (MEMS) are components with micron-scale moving parts based on materials and processes of microelectronics fabrication. This is a good example of on-chip integration of electronics, microstructures, microsensors and microactuators. Accurate simulation of MEMS requires precise modeling of all effects of mechanical and damping forces, electrostatic forces and inner stresses, heat transfer, thermal expansion, piezoelectric stresses etc. Modern methodology of MEMS design implies that the entire MEMS can be investigated only at higher abstraction levels such as schematic and system ones, where accurate macromodels can be used [1]. On the other hand, at component or device levels the physical behavior of three-dimensional continuums is described by partial differential equations (PDE) easily solvable by Finite Element or Finite Difference Element Methods (FEM or FDM) [2,3], available in ANSYS –like software. Component level simulations are classified in single domain and coupled domain simulations, both being very computer timeconsuming. The goal of this chapter is to consider methods of automatically obtaining macromodels of MEMS and their mechanical or non-electric components from ANSYS models as equivalent electric circuits or low order differential ordinary equations for further use in circuit design software. This can be done by using different model order reduction techniques developed in recent years. When dealing with the modern MEMS, the possibility for using a single environment to simulate objects, where different physical processes such as electrical, mechanical, optical, thermal etc. take place, plays an important role. Here we have to represent different subsystems of the initial MEMS as equivalent models of the same physical nature permitting to combine them for solution in a single computational process. After that, the complete behavioral model of the entire MEMS and its subsystems can be compiled either in VHDLAMS language (as sets of ODE) or in SPICE-like language (as equivalent electric circuits). The Microsystems design exploits various analytical and numerical methods for virtual prototyping of MEMS. It also demands for libraries of electromechanical, optical models and microfluid components, including springs, bulks, buffers, capacitors, inductances, operational amplifiers, transistors etc. Three basic possible approaches of MEMS design procedure are illustrated below: FEM/FDM Model, Reduced Order Model (ROM), Coupled system-level model.


Introduction
Micro-electro-mechanical Systems (MEMS) are components with micron-scale moving parts based on materials and processes of microelectronics fabrication.This is a good example of on-chip integration of electronics, microstructures, microsensors and microactuators.Accurate simulation of MEMS requires precise modeling of all effects of mechanical and damping forces, electrostatic forces and inner stresses, heat transfer, thermal expansion, piezoelectric stresses etc. Modern methodology of MEMS design implies that the entire MEMS can be investigated only at higher abstraction levels such as schematic and system ones, where accurate macromodels can be used [1].On the other hand, at component or device levels the physical behavior of three-dimensional continuums is described by partial differential equations (PDE) easily solvable by Finite Element or Finite Difference Element Methods (FEM or FDM) [2,3], available in ANSYS -like software.Component level simulations are classified in single -domain and coupled -domain simulations, both being very computer time-consuming.The goal of this chapter is to consider methods of automatically obtaining macromodels of MEMS and their mechanical or non-electric components from ANSYS models as equivalent electric circuits or low order differential ordinary equations for further use in circuit design software.This can be done by using different model order reduction techniques developed in recent years.When dealing with the modern MEMS, the possibility for using a single environment to simulate objects, where different physical processes such as electrical, mechanical, optical, thermal etc. take place, plays an important role.Here we have to represent different subsystems of the initial MEMS as equivalent models of the same physical nature permitting to combine them for solution in a single computational process.After that, the complete behavioral model of the entire MEMS and its subsystems can be compiled either in VHDL-AMS language (as sets of ODE) or in SPICE-like language (as equivalent electric circuits).The Microsystems design exploits various analytical and numerical methods for virtual prototyping of MEMS.It also demands for libraries of electromechanical, optical models and microfluid components, including springs, bulks, buffers, capacitors, inductances, operational amplifiers, transistors etc.Three basic possible approaches of MEMS design procedure are illustrated below: FEM/FDM Model, Reduced Order Model (ROM), Coupled system-level model.

FEM/FDM model
MEMS typically involve multiple energy domains such as kinetic energy, elastic deformation, electrostatic or magneto static stored energy and fluidic interactions.The difficulty in the modeling of MEMS devices is mainly caused by the tight coupling between the multiple energy domains.Individual physical effects are governed by partial differential equations (PDE), typically nonlinear.When these equations become coupled, the computational challenges of highly meshed numerical simulation become formidable.FEM relies on highly localized interpolation functions (or mesh element functions) for approximation of the solution of PDE.These mesh element functions are generated by meshing the domain of interest and parameterize the desired solution locally on each mesh element.This parameterized solution converts a continuous (PDE) problem to a coupled system of ordinary differential equations (ODE) that can be integrated in time.The resulting ODE system usually has many degrees of freedom (perhaps several variables per mesh element).If a fine mesh is required, the problem size grows rapidly, with a corresponding rapid growth in computational cost for explicit dynamic simulation.Consequently, it is very expensive to use FEM model in system-level simulations during MEMS iterative design.As a result, FEM models are mostly used to analyze the performance of MEMS components and to couple their multiphysics effects.By reading ANSYS binary FULL file it is possible to assemble a MEMS component statespace model in the form of first order systems or second order ordinary differential equations (ODE) Mx'' + Dx' + Kx =Bf , Y=Q T x +R T x' , where A r , E r ,C r , B r , M ,D, K, B, C-are the system matrices, B r , B are the input and the C r , Coutput matrices, f is input force.In mechanics matrices M, D and K are known as the mass, damping and stiffness matrices correspondingly.Usually damping is included in the model as Rayleigh damping.The damping matrix D is computed as a linear combination of the stiffness K and the mass M matrices: where , are constant coefficients.
In (1) the state space vector z is defined through the unknowns deflections u(x,t) and pressures p(x,y,t) into the node points being automatically generated in MEMS structure: By defining (4) second equations ( 2) can be transfer to the first (1).
The FULL file contains all the information about the system: the system element matrices, Dirichlet boundary conditions, equation constrain and the load vector.It is generated using ANSYS partial solver, which enables to assemble system element matrices for the desired analysis without solving them and it therefore computationally fast.The speed of the reading operation has been optimized taking into account that the element matrices are sparse.The load vector directly gives the matrix-vector product Bf and thus describes the distribution of all loads being applied.In order to obtain the B matrix, and thus being able to modify the inputs singularly, it is necessary to repeat the partial solution for each input of interest.

ROM (Reduced Order Model)
It would be easier and more intuitive for the designer to explore the design space if the MEMS model had only a few variables with a clear relationship between them and the overall device performance.Reduced-order models (ROM), also called macromodels, lend themselves very well to these purposes.The main idea behind the reduced order model is that the number of ordinary differential equations (ODE) needed to simulate the system has been reduced from perhaps many thousands in the case of the full FEM simulation, to just a few basis function coordinates (fig.1).Fig. 2. Compact reduced order model in MEMS design [33] www.intechopen.com

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. www.intechopen.com 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) q i , basis function (an eigenvector) φ i (x, y, z)) of the mechanical structure to represent its deflection u: where u eq 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 φ T i from the left, and using the orthogonality of eigenvectors, the equation ( 2) can be reduced to 22

2( )
The eigenvectors satisfy the orthogonality conditions With the normalization of the eigenvectors φ T k Mφ k ≡ 1, it can be shown that 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 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 H i 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 /f i = 0.01).Also the individual modal transfer functions (7) are plotted from which the complete transfer function is computed.By the way the eigenvalues q k and the corresponding eigenvectors φ k for k = 1, 2, . .., n can be found as an solution of the modified eigenvalue problem The relationships between natural frequencies f k , logarithmic decrements k and the eigenvalues are Fig. 6.Example of the modal decoupling in a transfer function [7].

Modal ROM based on proper orthogonal decomposition methods
Reduced order modeling using modal basis functions was originally developed by [11][12][13] and has been continuously improved by several authors.The basis functions (,,) i xyz  in ( 5) may be chosen in two ways by:  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 , particularly when irregular geometries are involved. 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 mode shapes of the device structural elements.The choice of orthogonal basis functions φ k 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 N s different times during these simulations, and the sampled distributions are stored as a series of vectors, {u i }, 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: 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: 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 V T 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 UU T [NxN] and U T U [n x n].The eigenvectors of U T U make up the columns of V, the eigenvectors of UU T make up the columns of W. Also, the singular values in S are square roots of eigenvalues from UU T or U T U. The singular values are the diagonal entries of the S matrix and are arranged in descending order: 12  where v i is the columns of V.
After finding basis functions φ k the linear equations for modal amplitudes q k 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 correlation matrix entries where n is number of temporal snapshots and b k,i are eigenvectors of the matrix A, or the solutions of the equation Ab= b.
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].

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 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.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]: 1. Compute the linear modes φ i of the elastic problem.
2. Substitute u (x ,y, z ,t) in the governing equation for the deflection (e.g., the Euler-Bernoulli or Timoshenko equations for beams).3. Obtain a system of n-coupled second-order ordinary differential equations for the q i (t). 4. Solve the equations to compute the dynamic response either numerically or analytically.

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.): .
The deformation state of the structure given by N nodal displacements u i (i=1,2,…,N ) is now represented by a linear combination of n modes weighted by their amplitudes q j (j=1,2,…,n ) where n << N.
The governing equation of motion describing the ROM of electrostatic actuated MEMS structures in modal coordinates: where m j is the modal mass, ω j is the eigenfrequency, ξ j the linear modal damping ratio, W st is the modal strain energy function, C ks 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 f i is a local force acting at the i-th node.The current I k at each electrode k is defined by: 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 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 W st 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 m i 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: -the master node terminals which describe the displacement u i and the inserted forces F N,i at these nodes; -the modal terminals with the modal amplitude q i and modal force F M,i for the chosen modes; -the electrical terminals which provide the voltages V i and currents I i for the electrodes of the system.

Fig. 9. ROM144 functional block
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 W mech (q i , q j , q k ) and capacitances C op (q i , q j , q k )as well as their first derivatives with respect to the modal amplitudes q i , q j and q k using the information of the polynomials degrees defined in other packages.Modal ROM144 approach speeds up computations in 40 (,,) For the Galerkin's method the PDE residual (() ) Lu f  is orthogonal to each i of the basic functions in the operating range H: where L is a differential operator (possibly nonlinear), and f is an input vector.The basic functions (,,) i 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 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 p i 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: and act as the scaling function 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.

Moment matching based ROM
Moment matching model order reduction is based on the approximation of the original ndimensional 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 let's perform a Laplace transform to obtain its frequency domain transfer function 1 () ( ) www.intechopen.comand expand it into Taylor series as 0 () where 1 () Moment matching is directly connected to the Krylov subspace formed by the pair of matrices (A -1 , A -1 B) [24].The Krylov subspace is spanned by the column vectors in the following collection of matrices: where the column vectors are called the Krylov vectors.The q-th order Krylov subspace is denoted by 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 Here V is the orthogonal projection matrix that maps the n-dimensional state-space into a q dimensional state-space and satisfies V T V = I.If the columns of V are orthogonal and B is a column vector, it can be shown that the following identities hold [24]: 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 () () xt Vxt   and multiplying on V T both sides of the equations ( 17): where 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 є R n , so that the algebraic operations involved can be seen.
www.intechopen.comiv.For j = 2, . . ., q: Solve j v  from: For i = 1, . . ., j − 1: Note that the Arnoldi algorithm terminates when h jj = 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].

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: 21 () ( ) 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 K q (K −1 M,K −1 B) which is used for the projection matrix V building.The transfer function can be expanded into Taylor series as 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 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): 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 = Vx q [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].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.

Equivalent-circuit ROM
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]: where are equivalent matrices of capacitances, conductance and inductances.

,, CGL 
The elements of matrices , , CGL  are formed from the elements of the mass, damping and stiffness matrices in the following ways: (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 L ii , 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: , ; , 1(1) ; 1 / , , 1(1) , ; .
sequential excluding of the internal circuit nodes by application of Y-∆ (star-triangle) transformation [36,39 ]; building a macromodel as a four-terminal [40].
where Let's define V i as (29) and replace V i 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:   is a sum of all node conductance excluding i-th node, k1 is a number of nodes being connected to node 1.The equation ( 30) can be simplified: 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: and its decomposition into Taylor series will look like: 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 bth nodes during the elimination of i-th node, and the second term in (35) defines the value of capacitance.

RCi
   , the equation ( 33) is transformed into: 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 connected to i-th node through a capacitor when / 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: 1.Time constants should be calculated for all the nodes of the circuit.
2. All the nodes should be put into a priority queue, sorted by time constants, apart from input, output and ground nodes.3. Take the first i-th node in the queue (the one with the smallest time constant 4. Find the neighbors k of i-th node. 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 recalculated 7. Take i-th node from the head of the priority queue and jump to step 3.  Original branch Substitution Formula

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: (37) where I a , U a are current and voltage at the macromodel input, I b , U b 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]: , where Δ ij -an algebraic complement to a ij 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 a ij 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: , 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: it is possible to find parameters of the equivalent four-terminal (38).By choosing elements ,,, Since macromodel parameters are determined by sum of the real and imaginary parts: it is convenient to represent the circuit macromodel also as a sum of real and imaginary terms: Investigation the expression (42) for different frequencies confirms that for linear circuits at least for frequency range 0    , where ω 0frequency 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 Y ij parameters for the stable passive macromodels (a 0 >0) has to be positive while an imaginary part could be negative.If in (41) 0, the Y ij component of the four-terminal may be presented by parallel connection of conductance, capacitance and inductance [41]: 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).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.
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.Maximal error, % --0,01 0,3 9,9 0,003 Table 3. Frequency beam analysis results

Source
The macromodel of the four-terminal for the same beam is being developed in according with ( 38)-( 44 The results of simulations of the four-terminal equivalent circuit by NetALLTED [36,43]  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.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.

MEMS coupled system-level model
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 tightly intertwined at microscale.Typically, up to four different physical domains are usually considered in the analysis of microsystems: mechanical, electrical, thermal and fluidic.
For each of these separate domains, well-established reduced order modeling and analysis techniques are available.However, one of the main challenges in the field of microsystems engineering is to connect models for the behavior of the device in each of these domains to equivalent lumped or reduced-order models without making unacceptably inaccurate assumptions and simplifications and to couple these domains correctly and efficiently.Micromechanical membrane devices (capacitive pressure transducers, ultrasonic transducers), surface micro machined devices (RF switches, micro optical devices) as well as bulk micro machined devices (accelerometers, inclinometers, laser scanning mirrors) are driven or sensed by nearly parallel electrode pairs in many cases.The motion of these electrodes is strictly normal to the surfaces.It means that MEMS electrical parts of these MEMS had to be combined with mechanical ones (fig.16).Typically, block-diagram descriptions and lumped-element circuit models for components are connected into a full system.Mostly, this description is used for functional analysis of a design concept.The coefficients and electrostatic nodal force are obtained from the capacitancedisplacement function C(w) of the associated electrode portions and gap space.The function can be input by one of three means [44]:  as analytical function if the electrode portions make up a plate capacitor geometry with homogeneous intermediate field;  as polynomial approximation of a function given by data points;  as data table wherein the element subroutine interpolates values during solution.The electrostatic forces acting on the movable conductor of the device are included in the model as nonlinear input forces, which are applied on nodes distributed over the conductor surface.Nodes divide the surface into N smaller portions.A lumped force is applied to kth node at the center of each portion, in its preferential direction of movement x i .The capacitance C k between the k-th portion and the fixed electrode of the device is computed, for its undeformed configuration, using an electrostatic analysis The entity of each force f k is then approximated as: where d 0 is the initial distance between the conductor and the electrode, is the relative dielectric permittivity, k i w is the displacement of the k-th node along the direction x i and (V p -V n ) is the voltage difference between the structure and the fixed electrode.The capacitance C k can be calculated from If C k varies considerably with the deformation of the structure, then a series of electrostatic computations for different device deflections in its operation range can be performed.The results can be used for extracting the dependency C k ( k i w ) and calculating the electrostatic force.Fig. 16.MEMS system-level simulation approach [33] According to the assumption that conductors are equipotential, all the nodes connected to a certain conductor are subjected to the same voltage boundary conditions.The total current flowing in the conductor is simply given by the sum of the currents at those nodes: The nodes used for electrostatic force application can be also used for monitoring the distance between the movable structure and the electrode.When this is equal to the transduction gap, the contact condition is reached and contact forces with the stiffness of the contact K n , given by equation can be applied to nodes.Computation of the electrostatic forces adds some complexity to the development of the MEMS system-level model, but this is largely compensated by the speed-up simulation of the full model.There is the special ANSYS' transducer element TRANS126 which has the possibility to calculate the capacitance of a parallel plate capacitor model ( at particular nodes or as a whole) as [47] The element has also contact capabilities (47).It is possible to specify a minimal gap and a spring stiffness K n for the repelling force.The drawback of this kind of element is that it is limited to the case where the electrodes are (almost) parallel plates, so that the stroke/capacitance function can be evaluated from a single degree of freedom.But the extensions for rotation plates and for 2D cases were developed.The last one with a triangular shape (element TRANS109) is useful for simulating structures such as comb drivers and optical MEMS, in which capacitance between the device parts is generally a function of a two-directional displacement.TRANS126 and TRANS109 elements enable a huge reduction of the complexity of the system-level simulation.Let's consider for example the system-level macromodel of an ultrasonic transducer which has two plates with the bottom electrode area A c and the plate dimension L and which can be presented by nonlinear capacitance: (50) where C 0 is the smallest capacitance in the absent of voltage V: ); 0  is the absolute dielectric permittivity of the vacuum, 1  is the relative dielectric permittivity of the poly silicon , ins  is the relative dielectric permittivity of the insulator; w is a plate deflection .
The largest value of this capacitance corresponds to the value .The electrostatic force acting on the capacitor surfaces is the Coulomb force:  where R is the interaction force.The process of interaction simulated by the force R sometimes can be highly sensitive to the values of the model parameters (rigidity and damping factor), especially if .The interaction can induce high-frequency motions and slow down the rate of convergence considerably.Therefore, a good deal of attention must be paid to accurately modeling and representing this process using experimental data.If the input voltage is increased more, there will be no equilibrium and the plate collapse takes place.In this case, a hard stop or some other arrangement must be introduced to limit the plate motion.During the plate collapse, the difference between the electrostatic force and the elastic force of the spring will continue to increase.Therefore, when the plate drops down into the hard stop, it is not enough to reduce the input voltage below to release the plate.The input voltage should be reduced more to make the electrostatic force at least equal to the elastic force.Hence, the plate capacitance exhibits the hysteretic behavior with respect to voltage change.But if the insulation layer is rather thick its restrictive effect should be taken into account.If the maximal center displacement is equal to initial thickness air gap the moving plate touches the insulator top surface of the electrode when the input V n voltage reaches the value of ( max V may be calculated from simulation).But as soon as the input voltage will be decreased under this value the plate will leave the hard stop.So, its capacitance does not demonstrate the hysteretic behavior (fig.17).Parameter in (50) defines an ultrasonic transducer frequency band and can be calculated through the plate displacement and its velocity , which are defied from ROM equations in the following way: (51 The coefficient 3 appears in (51) due to the fact that a capacitance recharge to 98% for a time value which is equal approximately 3 .It is possible to see two included procedures in according to fig.18: one for development of ROM for an ultrasonic transducer plate, where a deflection and speed of central point`s deflection of transducer plate is calculated for the value V n , and second -for determine of MEMS system-equivalent capacity value (SLM), using values of plate central point coordinates .Then the cycle of calculations recurs whereupon.Providing a single representation of a MEMS operating in multiple physical domains, the electrical circuit approach is very convenient.Moreover, powerful mathematical techniques and circuit simulation programs are available for solving design tasks.It is possible to develop a library of schematic model for different MEMS elements and then use their combinations to build a system-level macromodels for entire rather complicated MEMS constructions.

Conclusion
In this chapter, the methods and issues encountered in the of MEMS macromodels at the system level have been presented.System level modeling is the highest and most abstract level of modeling.This level requires various devices` linking of MEMS component level models -both electronic and micromechanical -into a micro-electromechanical system.System-level models of MEMS components are needed to allow a fast and sufficiently exact investigation of their behavior to simulate entire MEMS.Starting point for the extraction of a reduced order model (ROM or a macromodel) is already its description with a large ODE system, which is typically derived using physical modeling techniques based on Finite Element Method (FEM) which is rather time consuming.Macromodels application allows the extraction of lower order ODE system that reproduces the input/output behavior with good accuracy.Particular attention has been posed in the chapter on the possibility to get a macromodel circuit presentation.There are special methods for generating ROM for MEMS components and entire MEMS based on FEM descriptions.To derive macromodels of smaller sizes different approaches (Modal decomposition, Moment matching, Equivalent circuit presentation) were developed.Usage of the reduced MEMS components models allows applying successfully modern circuit simulators in workflow for MEMS design on system level.Three automatic procedures to generate device reduced order macromodels, being based on full FEM/FDM models, were demonstrated in this chapter.Two of them are suitable for simulators with possibilities to get input information in the equation forms (ODE or OAE).The third one in opposite produces macromodels in circuit presentation and so it is more suitable for circuit simulators.The Modal ROM approach is based on using natural (modal or resonant) frequencies of MEMS structure and it is spread mostly in the USA and Asia.The Moment matching ROM approach is based on using the Krylov subspace for transfer functions and it is popular in Western Europe and Asia.The Equivalent circuit ROM approach is based on using a capacitive-inductive-resistive circuit model for mass, damping and stiffness matrices and it is used mostly in Eastern Europe.It is worth to notice that the Modal ROM approach requires some full ANSYS runs to perform a proper orthogonal decomposition during basic functions determination in opposite to the Moment matching and the Equivalent circuit ROM approaches for which it is enough to use ANSYS only for FEM model matrices building.It seems to be interesting and perspective trying to combine mentioned approaches, for example, to start with Krylov/Arnoldi reduction of ODE dimension, then to build the proper equivalent circuit for obtained ODE systems and finally to apply Y/Δ transformation or n-port transformation for further reducing macromodel order.

Fig. 1 .
Fig. 1.Reduced order model illustration Such the macromodel simulation can be very efficient computationally compared to the FEM model.A designer can use the FEM model for different component geometry and materials trying and the ROM model for investigation of different input forces effect (fig.2).

Fig. 4 .Fig. 5 .
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. 10 .
Fig. 10.Comparison of transient behavior (a) and transfer functions (b) for the full and reduced microgyroscope models [33].

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

Fig. 12 .
Fig. 12.A working node of the RLC-circuit: conductance are added between node 1 and all neighbors nodes eliminate a node from equivalent MEMS circuit with the smallest inaccuracy, the following two criteria should be taken into account: -a node should be fast; -time constants , LCi RCi   should not be compared (have almost the same values).
a b b /B i C ab =(b a C b ++b b C a )/B i b ab = b a b b /B i C ab = b a C b /B i b ab =b a b b /B i C ab = b a C b /B i C ab =b a *C b /B i C ab =C a C b /C C ab =C a g b /G i g ab =g a g b /G i g ab =b a g b /G i g ab =b a g b /G i respective component Y ij may be presented in the same way as parallel connection of conductance, capacitance and inductance but:

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

Fig. 16 .
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.3Hz and T2=4009.3Hz in according to the ANSYS analysis results (table3).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.

Fig. 18 .
Fig. 18.ROM-system-level model coupled simulationsInstead of using two sequence procedures mentioned above it is possible using functional possibilities of the circuit simulator NetALLTED to built a single system-level equivalent circuit model for an ultrasonic transducer by introducing directly into the equivalent-circuit ROM of mechanical MEMS part the additional arbitrarily connected element (a Depended Source) with an informative function which is determined by equation (50)[42].Optimization procedures of NetALLTED allow getting the desirable values of this transducer capacity and through it to get a desirable value of output signal of an ultrasonic transducer system-level model by the changing ROM parameters, which, in turn, are depended upon an ultrasonic transducer construction sizes and used material properties.

Table 2 .
Transformations for RCi case

Table 4 .
are given in table 4. Results of frequency analysis of the four-terminal equivalent circuit for a beam (44)w are the initial distance and the displacement between the plate.The element has two nodes; the gap distance is calculated as the sum of the initial displacement and the difference of the nodal displacements in the direction of the element.The force is calculated by equation being similar to(44)so that in the constant voltage case