A New MEMS Stochastic Model Order Reduction Method : Research and Application

Modeling and simulation of MEMS devices is a very complex tasks which involve the electrical, mechanical, fluidic, and thermal domains, and there are still some uncertainties that need to be accounted for during the robust design of MEMS actuators caused by uncertainmaterial and/or geometric parameters. According to these problems, we put forward stochastic model order reduction method under random input conditions to facilitate fast time and frequency domain analyses; themethodmakes use of polynomial chaos expansions in terms of the random input variables for the matrices of a finite element model of the system and then uses its transformation matrix to reduce the model; the method is independent of the MOR algorithm, so it is seamlessly compatible withMORmethod used in popular finite element solvers.The simulation results verify the method is effective in large scale MEMS design process.


Introduction
MEMS are attractive for many applications because of their small size and weight, which allow systems to be miniaturized [1].Modeling and simulation is at the basis of the prediction of the device behavior and optimization of its performance.The complexity of the simulation analysis adds to the complexity of the model.MEMS characterization often requires computationally expensive analysis such as transient analysis.Moreover, the derived models are often nonlinear.The primary source of nonlinearity is the electrostatic force, which couples the electrical and mechanical energy domains.Simulation of nonlinear models makes use of iterative time consuming algorithms.To improve the simulation computation efficiency, MOR (model order reduction) method was used to extract the lower order ODE system that reproduces the input/output behavior with good accuracy.The method is thus solely based on the mathematical properties of the original system and is, therefore, formal, robust, and in great part automatic [2,3].These properties render the use of mathematical model order reduction more and more popular in the study of MEMS devices.Variations during fabrication lead to uncertain material and/or geometric parameters causing a significant impact on MEMS device performance [4].While there have been significant advances in numerical simulation methods that allow better understanding of the underlying multiphysics [5][6][7], they, however, assume that the geometrical and physical properties of the device are known in a deterministic sense.Recently there have been efforts towards developing computational methods that can handle input uncertainties.A standard way to deal with these uncertainties is brute-force Monte Carlo simulations; these basically involve considering a large number (typically greater than 10000) of realizations (or samples) of the geometry and solving the deterministic problem for each one of these realizations [8].Several techniques have been developed, for improving convergence, for example, Latin hypercube sampling, [9,10], the quasi-Monte Carlo (QMC) method, and the Markov chain Monte Carlo method [11].The statistics such as the mean and standard deviation of the required output quantity such as the electrostatic force are then generated.This is a very time consuming process.Also, there is some uncertainty research work about MEMS: Kong et al. [12] studied the performance variability of a ceramic MEMS actuator under random variations in the shape of the actuator and the air gap in the condenser.Han and Kwak [13] presented the use of robust optimization during the design of a microgyroscope using MC simulations to compare predicted yields.Liu et al. [14] presented a robust design method to minimize the sensitivity of a laterally vibrating resonator against width variations due to fabrication errors.Stochastic model order reduction can offer an efficient way for optimization of dynamic problems [15][16][17].This has led to a lot of efforts towards developing MOR algorithms for variational and parametric uncertainty analysis.In [18], a variational balanced truncation method was introduced for model reduction of variable geometry interconnects.In [19], a method for model order reduction of RCL interconnects was described.Various algorithms for parametric model order reduction have been proposed in [20][21][22][23][24][25].
In this paper, we input parameter uncertainty.We put forward a new method for MEMS stochastic MOR computing process in the presence of input parameter uncertainty.Firstly, we expressed the random input variables as standard random variables by KL series expansion method; then, we expressed the random output variables as standard random variable polynomials by PC method; while obtaining the polynomials of output variables, performing deterministic MOR for each system to generate the corresponding transformation matrix at each point on the Smolyak sparse grid, at last, calculate the mean, standard deviation and other statistics of the system response.The remainder of this paper is organized as follows.In the next section, we represent random input variables matrix of MEMS models; Section 3 is dedicated to the PC representation of the output random variables; polynomial chaos expansion coefficient solution by sparse grid method was described in Section 4; stochastic reduced order model is described in Section 5.The computational gain provided by the sparse PC expansions simulation is illustrated in Section 6 by numerical examples.

Representing Input Variable
Stochastic of MEMS Model

MEMS Stochastic
Process.MEMS stochastic process refers to a cluster of random variables defined in the parameter set; every point at the parameter set is corresponding to a random variable.One-dimensional stochastic process can be viewed as a natural extension of random vectors.If a stochastic process was denotes as {(),  ∈ }, to describe its probabilistic properties, the most important aspect is the distribution function of random variables at time ; we denote the function as  (, ) =  [ () < ] ,  ∈ . ( If we consider the relationship of the stochastic process between any two random variables, then we can get the following expression: And if there are  variables in MEMS system, then the stochastic process can be described as ∈ . (3)

MEMS Stochastic Field
. MEMS stochastic field is the natural expansion of stochastic process in special field; to a stochastic field, its basic parameters is special variable  = {, , }, so we can define the stochastic field as series of random variables on a field set.And there are corresponding variables at every point   of the field set.In fact, we just take into consideration special variables as basic parameter's stochastic field denoted as {();  ∈ }, where  is defined area.To a uniform stochastic field, because of its mean function is a constant, so it can be turned into the following form: where  0 () is mean field function and   () is zero-mean stochastic field.Obviously, the covariance function of  0 () is the same as the correlation function of   ().And we can use limited distribution function to represent the probability of random structure, like the following expression: where   is value limit of variable.

Stochastic Input Variables of MEMS Model.
The stochastic problem of MEMS can be described as the following equation: where  represent the solution of response-problem; () is random force; (, ) is the MEMS material stochastic field; and  is the stochastic event.From expression (5), we discrete the stochastic field and realize the conversion from stochastic field to discrete random variable set.Then we will use KL expansion to make these variables independent from each other.The KL expansion can be written as where () is the mean of the random process,   and   () are correlation function, ( 1 ,  2 ) are the eigenvalues and eigenfunctions, respectively, where  1 and  2 are the spatial coordinates.{  ()} forms a set of uncorrelated random variables. and   form the eigenvector-eigenvalue pair of the covariance kernel that satisfy At the same time, there is a relationship between the cofeature value function of the variance function: In practice the expansion in ( 7) is truncated after a finite number of terms , which leads to a truncation error   .
As compared to other expansion methods, which use some orthonormal functions {  }, the KL expansion is optimal in the sense that the mean-square error ∫ Ω  2   is minimized.

MEMS Output Variables by Polynomial Chaos Expansion
The polynomial chaos expansion is a spectral expansion of the random process in terms of the orthogonal polynomials in multidimensional random variables.Let {  ()} be a set of orthonormal Gaussian random variables.Using this, the polynomial chaos expansion of a second-order output variables random process (, ) can be described as where Γ  (  1 ,   2 , . . .,    ) denotes the polynomial chaos of order  in terms of the multidimensional Gaussian random variables.For convenience, (11) can be written as where there is a one-to-one correspondence between the functions Γ[⋅] and Ψ[⋅] and between the coefficients.Finite number of random variables is used in the expansion to represent finite number of random parameters in the system.Also, the order of the polynomial used in the expansion is restricted to .Thus, the expansion in ( 12) can now be written as The total number of terms included in the polynomial chaos expansion ( + 1) depends on both the dimensionality  and the highest order  of the multidimensional polynomials used, which can be given as The eigenvalues and coefficients of eigenvectors are random processes; they belong to ; we write the random eigenvalues and eigenvectors using the Karhunen-Loeve truncated expansion with Hermite polynomial chaos: where { ∈{1,...,} ∈{1,...,} } are coefficients of the truncated expended,   () ∈ is family of Hermite polynomial chaos.Introducing a classical normalization condition, {  }  [  ]{  ()} = 1 ∀ ∈ {1, . . ., }, we can get We also can express the other input random element of the problem with the same approximation on Ψ.

Polynomial Chaos Expansion Coefficient Solution by Sparse Grid Method
The most important aspect of polynomial chaos method is to calculate the each coefficient of the polynomial chaos expansion.On a multidimensional model containing random variables polynomial chaos unfolds; the output at some point can be calculated with its polynomial coefficients chaos which unfolds.In a vector space expansion, if we us {  ()} as base, and every sample {  ()} has corresponding point, we call these points as configuring points.There are many methods to compute the configuring points, such as tensor product method and Stroud-2 method.We use sparse grid method which derived from Smolyak method.This approach can significantly reduce the required number of points and this algorithm provides a linear combination of the tensor product of the interpolation error which can be made very close to full tensor product method.Firstly, let: Δ  =   −  −1 ; then where (, ) is a linear function and (, )() is determined by the limited point function value.When  0 = 0 and || =  1 + ⋅ ⋅ ⋅ +   especially, there is And all polynomial chaos expansion coefficient solution is just processed by these grids points, and the point's sets are where the number of points used in each dimension vary depending on the method of polynomials and quadrature.

Stochastic Reduced Order Model
The stochastic reduced order model method's basic theory is when we get the each point on the Smolyak sparse grid; then we can calculate the full-finite element system matrices at the point, and we can perform deterministic MOR for each system to generate the corresponding transformation matrix.Using this matrix information, we can compute the augmented system using the coefficients.

Deterministic Reduced Order
Model.MEMS devices can be modeled by partial differential equations (PDEs).To simulate such models, spatial discretization via, for example, finite element discretization is necessary, which results in a system of ordinary differential equations (ODEs) or differential algebraic equations (DAEs).After spatial discretization, the number of degrees of freedom usually is very high.Therefore, it is time consuming to simulate the large-scale systems of ODEs or DAEs.Developed from well-established mathematical theories and robust numerical algorithms, model order reduction (MOR) has been recognized as being very efficient in reducing the simulation time of largescale systems.Through MOR, a small system of ODEs with reduced number of equations (reduced model) is derived.The reduced model is simulated instead, and the solution of the original PDEs or ODEs can then be recovered from the solution of the reduced model.As a result, the simulation time of the original large-scale system can be shortened by several orders of magnitude.The reduced model as a whole can also replace the original system and be reused for many times during the design process, which can save much time further.To illustrate the MOR formulation procedure and moment matching we will take the first-order system's model order reduction shown in (1) as an example and rewrite it as follows: where Here the matrices , , and  are positive semidefinite matrices.The model order reduction based method can be written as follows.
Step 2. To implement Arnoldi process, specific algorithm is as follows: (a) given matrix,  ∈   ,  ∈  × , and QR decomposition of matrix , let  =  0 ; Step 3.According to the transformation matrix, then get the original system of reduced order system: where x() ∈   , Ẽ =   , Ã =   , B =   .
PRIMA algorithm fully considered the system structure characteristics, so it keeps the passivity of the system and keeps the reduced order model system on the basis of structural characteristics, which has higher calculation accuracy.

Stochastic Reduced Order Model.
To get stochastic reduced order model of MEMS, the first thing is to represent "randomness" in the original system matrices and the state vector; this is done by making use of polynomial chaos expansion.Suppose two input variables are assumed to be random and uniformly distributed.Considering a linear expansion in both of the input random variables, and using one-dimensional polynomial chaotic maps, the Ψ[⋅] in (13) can be written as one-dimensional Hermite polynomial: Assuming the input of a linear variable, there are where  0 ,  1 ,  2 ,  0 ,  1 ,  2 ,  0 ,  1 ,  2 is reduced order model coefficient matrix of deterministic model. 1 ,  2 are two-dimensional random space orthogonal polynomials.Once the coefficient matrix was determined by sparse grid method, the stochastic mapping reduced order model can be written as (11): Let Equation ( 18) may be cast in the following "deterministic" form: Equation ( 19) is a deterministic system which can be easily used for time and frequency domain analyses.Then we can use deterministic model order reduction method to solve the stochastic reduced order model.

Computation of the Final Results: Expected Values and
Standard Deviations.Once we have solved the  systems numerically, that is to say, when we have got the matrices Ã and B , then we plug back the results into (11) and (12) to get an approximation of the probability density functions of the random eigenvalues and eigenvectors.As we have seen before this is the solution to the full problem with just one chosen mistuning parameter projected onto the Ψ subspace.Nevertheless we are not directly interested in this solution.Indeed, practically what is useful is the average value and standard deviation of each eigenvalue and eigenvector of the mistuned systems.But this is easy to get from the approximations of the probability density functions we have got.We directly use the linearity of the operator and the orthogonality of polynomial chaos on ( 16) to get From ( 27) we get Hence using the linearity and orthogonality properties as previously mentioned, we get Finally using ( 28) and ( 29), we directly get the standard deviation: We do exactly the same thing for eigenvectors and we directly get the expected value from the linearity and orthogonality properties applied on (28): . (32)

Numerical Studies: MEMS Examples
In order to illustrate the model reduction technique, we give an example of the well-known beam structure in a fluid environment.Figure 1 shows the beam structure; when a voltage is applied, the top plate of the structure bends downwards due to the resultant electrostatic force.Also when the beam bends, the pressure distribution of the ambient air under the beam increases.This pressure increase produces a backward pressure force that damps the beam motion.The beam structure has been used in many sensor applications.
The beam can be modeled by coupling the Euler beam equation with the electrostatic force and the Reynolds squeeze-film damping equation as follows: Among them (, ) is the immunity of the  direction,  is young's modulus,  = ℎ 3 /12-moment of inertia's is residual stress (Figure 5),  is the density,  0 is dielectric constant,  is the air damping of the pressure,  0 is the ambient pressure,  is the air viscosity coefficient,  is the distance from the plate under initial state, and  is air mean free path (see Table 1).
Let ỹ = ((, ) − )/, p(, ) = ((, ) −  0 )/ 0 , substitute them into (3) and in  and  0 using Taylor series expansion, and take position, speed, and pressure as state variables like we can get the micro beam for original system state equation: To realize the stochastic reduced order model simulation, we firstly established the corresponding finite element model of the beam, to get the initial mathematical model such as rigidity matrix and mass matrix, which can be used as random input parameters, do standard KL series expansion with these matrix and do polynomial chaos expansion with output matrix; after expansion finish, we use Smolyak method to decide which point should be used as reduced order model simulation; then we use SPRMI method at the Smolyak point; at last, we can get the probability distribution from order reduction model simulation results.In MEMS, because micromechanical usually driven by electrostatic force, the machine is usually a conductor in static electric field, subject to the effect of electrostatic force deformation occurring, and the structure size and the location of the components will influence the distribution of electric field in turn, so it is a kind of strong coupling.Here we adopt finite element method combined with the boundary element method to  solve the force electric coupling problem of MEMS, solving process as follows.
Step 1.Under the condition of the plate without displacement, giving an initial voltage on the plate and solving problems of electrostatic field with the boundary element, we can get the distribution of electrostatic force of plate.
Step 2. Solve the displacement of plate under electrostatic force through the finite element method; if it reaches precision displacement, skip to Step 4; otherwise, go to Step 3.
Step 3. Based on the displacement of the plate to form the boundary of the electrostatic field, use the boundary element method to solve the problem of electrostatic field, obtain the distribution plate of the electrostatic force, and let  =  + 1; go to Step 2.
Step 4. The process is over.
We use ( + 1) × ( + 1) mesh as shown in Figure 2, where  represents the number of inner grid points in the  direction and  is the number of inner grid points in the  direction, according to the two-dimensional model establishment of beam parameters as shown in Figure 2. The beam was divided with 2D quadrilateral element PLANE82, and the air medium was divided with PLANE121; displacement constraints were applied to fix the microbeam anchors in ANSYS (Figure 4); each degree of freedom beam base and left side of beam were restricted; the finite element model is shown in Figure 3. Then we applied effect of voltage on the beam structure.After the mesh is generated, we then project the unknowns (, ) and (, , ) onto the mesh points and apply the trapezoidal rule to discretize the special integral   operator and the central difference method to discretize the special derivative operators as follows: Accelerometer beam clamped at one end, so the boundary conditions are  0 = 0, −1 = where  = [  0 0  ] ∈  × ,  = [   −  0 ] ∈  × , and  = [  1 0 ] ∈  × () are the microdisplacement of the beam end and V() is the input voltage of the system.
Applying the boundary conditions, the finite element model is solved and can change the endpoint microbeam under different voltage conditions.We get the deflection of the beam's end point as shown in Figure 6.The voltage load vector along with the mass and stiffness matrices were extracted as input random parameters.
A generalized polynomial chaos expansion is now generated to represent the dependence of the maximum displacement of the membrane as a function of the input parameter space.This includes geometric parameters thickness ℎ, length , width , Young's modulus , and the initial gap  between the membrane and the pull-down electrode.The variability in the thickness and gap is estimated from known processing conditions.The parameter limits are selected based on the requirement; the device cannot experience pull-in instability over the time span of interest.The parameter limits and distribution types are presented in Table 2.
Then making use of Smolyak algorithm consisting of only 29 grid points in the two-dimensional random parameter spaces for numerical integration, the grid result is shown in Figure 7.
Then, we compute 29 times of model reduction simulation at Smolyak point; the results were shown in Figure 8. From the figure, we learn when time ≤ 0.7 s.The most simulation results are almost the same.After that, we get kinds of deflection.
We also compute that the resulting probability distribution of the beam maximum deflection is presented in Figures 9,10,and 11.In Figure 9, we just show geometric parameters thickness ℎ, length , and width  as random parameters, and Figure 10 shows Young's module as random parameter, and Figure 11 shows the initial gap distance as random parameters.

Conclusions
We have developed a new stochastic model order reduction frame to MEMS model uncertainty and application in single fixed micro-beam.In the computing process, firstly, we construct finite model, then extract the matrix of analysis, get the random input matrix, and use KL series expansion method for input parameters matrix; then, we express the random output variables as standard random variable polynomials by PC method; while obtaining the polynomials of output variables, we use Smolyak method to compute coefficients of PC equation, then take the reduced matrix to PRSIM model reduction method, and compute the statistics such as the mean and standard deviation of the desired system response and deflection of microbeam.We can get that all the computational effort is focused on the estimation of a small set of PC coefficients in the proposed methodology.This leads to a considerable gain compared to the techniques based on Monte Carlo simulation.

Figure 6 :
Figure 6: The endpoint deflection of the beam.

Table 1 :
The parameters of single-ended fixed electrostatic actuating microbeam structure.

Table 2 :
The GPC response inputs parameters.