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

Authors: Bian Xiangjuan; Youping Gong; Chen Guojin; Lv Yunpeng
Journal: Journal of Sensors
 (2015)

---

## Abstract

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 uncertain material 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; the method makes 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 with MOR method used in popular finite element solvers. The simulation results verify the method is effective in large scale MEMS design process.


---

## Body

## 1. 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–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–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–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; Section3 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.


## 2. Representing Input Variable Stochastic of MEMS Model

### 2.1. MEMS Stochastic Variables

The probability of the random events can be described as a positive number less than 1; the number which describes the degree of random eventA is called the sample space Ω, denoted as P(A). Probability measure events may occur only giving a measure of the standard and, for the sample point of the sample space, a similar standard is still not given. It has two basic futures: (1) MEMS stochastic variable is a real single-valued function of sample points; (2) to any real number x, ω:ξω<x is a random event. When we get the definition of random variable, then a basic event can be represented by a determined value, denoted as X=x; and any event can use domain value of random variables; we can note the following: x1≤X≤x2. The likelihood of the event occurring with probability can be described as Px1≤X≤x2.


### 2.2. 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 asXt,t∈T, to describe its probabilistic properties, the most important aspect is the distribution function of random variables at time T; we denote the function as(1)Fx,t=PXt<x,t∈T.If we consider the relationship of the stochastic process between any two random variables, then we can get the following expression:(2)Fx1,t;x2,t2=PXt1<x1,Xt2<x2,t∈T.And if there are n variables in MEMS system, then the stochastic process can be described as(3)Fx1,t1;x2,t2,…,xn,tn=PXt1<x1,Xt2<x2,…,Xtn<xn,hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhht∈T.


### 2.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 variableu=x,y,z, so we can define the stochastic field as series of random variables on a field set. And there are corresponding variables at every point ui of the field set. In fact, we just take into consideration special variables as basic parameter’s stochastic field denoted as Bu;u∈D, where D 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:(4)Bu=B0u+Bσu,where B0u is mean field function and Bσu is zero-mean stochastic field. Obviously, the covariance function of B0u is the same as the correlation function of Bσu. And we can use limited distribution function to represent the probability of random structure, like the following expression:(5)Fβ1,u1;β2,u2;…;βn,un=PBu1<β1,Bu2<β2,…,Bun<βn,where βi is value limit of variable.


### 2.4. Stochastic Input Variables of MEMS Model

The stochastic problem of MEMS can be described as the following equation:(6)Ku,ωx,θ=Fθ,where u represent the solution of response-problem; Fθ is random force; ωx,θ 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(7)ωx,β=ω-x+∑i=1∞λiζiθfix,where ω-x is the mean of the random process, λi and fix are correlation function, C(x1,x2) are the eigenvalues and eigenfunctions, respectively, where x1 and x2 are the spatial coordinates. ζiθ forms a set of uncorrelated random variables.fi and λi form the eigenvector-eigenvalue pair of the covariance kernel that satisfy(8)∫ΩCx1,x2fix2dx2=λifix1.At the same time, there is a relationship between the cofeature value function of the variance function:(9)∫Ωfixfjxdx=δij=1i=j0i≠j,(10)ζiθ=1λi∫Ωω-xθ-ω-θfixdx.In practice the expansion in (7) is truncated after a finite number of terms M, which leads to a truncation error εM. As compared to other expansion methods, which use some orthonormal functions fi, the KL expansion is optimal in the sense that the mean-square error ∫ΩεM2dx is minimized.


## 2.1. MEMS Stochastic Variables

The probability of the random events can be described as a positive number less than 1; the number which describes the degree of random eventA is called the sample space Ω, denoted as P(A). Probability measure events may occur only giving a measure of the standard and, for the sample point of the sample space, a similar standard is still not given. It has two basic futures: (1) MEMS stochastic variable is a real single-valued function of sample points; (2) to any real number x, ω:ξω<x is a random event. When we get the definition of random variable, then a basic event can be represented by a determined value, denoted as X=x; and any event can use domain value of random variables; we can note the following: x1≤X≤x2. The likelihood of the event occurring with probability can be described as Px1≤X≤x2.


## 2.2. 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 asXt,t∈T, to describe its probabilistic properties, the most important aspect is the distribution function of random variables at time T; we denote the function as(1)Fx,t=PXt<x,t∈T.If we consider the relationship of the stochastic process between any two random variables, then we can get the following expression:(2)Fx1,t;x2,t2=PXt1<x1,Xt2<x2,t∈T.And if there are n variables in MEMS system, then the stochastic process can be described as(3)Fx1,t1;x2,t2,…,xn,tn=PXt1<x1,Xt2<x2,…,Xtn<xn,hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhht∈T.


## 2.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 variableu=x,y,z, so we can define the stochastic field as series of random variables on a field set. And there are corresponding variables at every point ui of the field set. In fact, we just take into consideration special variables as basic parameter’s stochastic field denoted as Bu;u∈D, where D 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:(4)Bu=B0u+Bσu,where B0u is mean field function and Bσu is zero-mean stochastic field. Obviously, the covariance function of B0u is the same as the correlation function of Bσu. And we can use limited distribution function to represent the probability of random structure, like the following expression:(5)Fβ1,u1;β2,u2;…;βn,un=PBu1<β1,Bu2<β2,…,Bun<βn,where βi is value limit of variable.


## 2.4. Stochastic Input Variables of MEMS Model

The stochastic problem of MEMS can be described as the following equation:(6)Ku,ωx,θ=Fθ,where u represent the solution of response-problem; Fθ is random force; ωx,θ 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(7)ωx,β=ω-x+∑i=1∞λiζiθfix,where ω-x is the mean of the random process, λi and fix are correlation function, C(x1,x2) are the eigenvalues and eigenfunctions, respectively, where x1 and x2 are the spatial coordinates. ζiθ forms a set of uncorrelated random variables.fi and λi form the eigenvector-eigenvalue pair of the covariance kernel that satisfy(8)∫ΩCx1,x2fix2dx2=λifix1.At the same time, there is a relationship between the cofeature value function of the variance function:(9)∫Ωfixfjxdx=δij=1i=j0i≠j,(10)ζiθ=1λi∫Ωω-xθ-ω-θfixdx.In practice the expansion in (7) is truncated after a finite number of terms M, which leads to a truncation error εM. As compared to other expansion methods, which use some orthonormal functions fi, the KL expansion is optimal in the sense that the mean-square error ∫ΩεM2dx is minimized.


## 3. 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ζiθ be a set of orthonormal Gaussian random variables. Using this, the polynomial chaos expansion of a second-order output variables random process u(x,θ) can be described as (11)ux,θ=a0xΓ0+∑i1=1∞ai1xΓ1ζi1θ+∑i1=1∞∑i2=1∞ai1ai2xΓ2ζi1θ,ζi2θ+⋯,where Γn(ζi1,ζi2,…,ζin) denotes the polynomial chaos of order n in terms of the multidimensional Gaussian random variables. For convenience, (11) can be written as(12)ux,θ=∑i=0∞a^ixΨζθ,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 p. Thus, the expansion in (12) can now be written as(13)ux,θ=∑i=0Na^ixΨζθ.The total number of terms included in the polynomial chaos expansion (N+1) depends on both the dimensionality n and the highest order p of the multidimensional polynomials used, which can be given as(14)N+1=n+p!n!p!.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:(15)λ~kζ-≈λk·∑b=1mabk·ϕbζ-∀k∈1,…,nU~kζ-≈∑r=1nβ~rkζ-·Urwithβ~rkζ-≈∑l=1mβlr·ϕlζ-∀k∈1,…,n,where βi∈1,…,mj∈1,…,n are coefficients of the truncated expended, ϕiζ-i∈n is family of Hermite polynomial chaos. Introducing a classical normalization condition, UkTInU-kζ-=1∀k∈1,…,n, we can get(16)U~kζ-≈Uk+∑r=1n∑l=1mβlk·ϕζ-·Ur.We also can express the other input random element of the problem with the same approximation on Ψ.


## 4. 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ζiθ as base, and every sample ζiθ 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: Δi=ui-ui-1; then(17)Aq,N=∑q-N+1i≤q-1q-iN-1q-iui1⊗⋯⊗uiN,where Aq,N is a linear function and Aq,Nf is determined by the limited point function value. When u0=0 and i=i1+⋯+iN especially, there is(18)Aq,N=∑i≤qΔi1⊗⋯⊗ΔiN.And all polynomial chaos expansion coefficient solution is just processed by these grids points, and the point’s sets are(19)θN≡Hq,N=Uq-N+1≤i≤qθi1×⋯×θNi1,where the number of points used in each dimension vary depending on the method of polynomials and quadrature.


## 5. 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.


### 5.1. 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 large-scale 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:(20)Edx(t)dt=Ax(t)+Bu(t)yt=BTxt,where E=M00H∈Rn×n, A=NF-FT0∈Rn×n, B=B10∈Rn×p. Here the matrices M, H, and N are positive semidefinite matrices. The model order reduction based method can be written as follows.Step 1.
We build block Krylov subspace:KrA-1E;A-1B.Step 2.
To implement Arnoldi process, specific algorithm is as follows:(a)
given matrix,A∈Rn, B∈Rn×p, and QR decomposition of matrix B, let B=V0T;(b)
computeV^1=AV0-V0H00, where H00=V0TAV0;(c)
computeV^2=AV1-V1H11-V0H01, where H11=V1TAV1, H01=V0TAV1;(d)
then compute the following parameter intrun:V^r=AVr-1-Vr-1Hr-1,r-1-V1H1,r-1-V0H0,r-1, where Hr-1,r-1=Vr-1TAVr-1,…,H1,r-1=V1TAVr-1, H0,r-1=V0TAVr-1.Step 3.
According to the transformation matrix, then get the original system of reduced order system:(21)E~dx~tdt=A~xt+B~uty~t=B~Tx~t,where x~t∈Rrp, E~=VTEV, A~=VTAV, B~=VTB.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.


### 5.2. 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:(22)Ψ0ζ=1,Ψ1ζ=ζ,Ψ2ζ=ζ2-1,Ψ3ζ=ζ3-3ζ,Ψ4ξ=ζ4-6ζ2+3,…Assuming the input of a linear variable, there are(23)E~r=E0+E1ζ1+E2ζ2,A~r=A0+A1ζ1+A2ζ2,B~r=B0+B1ζ1+B2ζ2,where E0,E1,E2,A0,A1,A2,B0,B1,B2 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):(24)E0E13E23E1E00E20E0dx~0(t)dtdx~1(t)dtdx~2(t)dt=A0A13A23A1A00A20A0x~(t)+B0B13B23B1B00B20B0u(t)y~t=B0B13B23B1B00B20B0Tx~t.Let(25)Esr=E0E13E23E1E00E20E0,Asr=A0A13A23A1A00A20A0,Bsr=B0B13B23B1B00B20B0.Equation (18) may be cast in the following “deterministic” form:(26)E~srdx~(t)dt=A~srx~(t)+B~sru(t)y~t=B~srTx~t.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.


### 5.3. Computation of the Final Results: Expected Values and Standard Deviations

Once we have solved then systems numerically, that is to say, when we have got the matrices A~sr and B~sr, 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(27)λ~kζ-≈λk·∑b=1mabk·ϕb=λk·∑b=1mabk·δ1b=λk·a1khhhhhhhh∀k∈1,…,n(28)U~k=Uk+∑r=1n∑l=1mβlr·ϕl·Ur=Uk+∑r=1nβlr·Ur.From (27) we get(29)λ~k2=λk·a1k2.Hence using the linearity and orthogonality properties as previously mentioned, we get(30)λ~k2=λk·a1k2+λk2·∑b=2mabk2·ϕb2.Finally using (28) and (29), we directly get the standard deviation:(31)δλ~k=λ~k2-λ~k2=λ·∑b=2mabk2·ϕb2.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)σU~kq=∑r=1n∑l=2mβlr2ϕl2·Urq2∑r=1r1,r2≠kn∑l=2mβlr1βlr2ϕl2·Ur1qUr2qhhhhh+2·∑r=1r1,r2≠kn∑l=2mβlr1βlr2ϕl2·Ur1qUr2q1/2.


## 5.1. 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 large-scale 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:(20)Edx(t)dt=Ax(t)+Bu(t)yt=BTxt,where E=M00H∈Rn×n, A=NF-FT0∈Rn×n, B=B10∈Rn×p. Here the matrices M, H, and N are positive semidefinite matrices. The model order reduction based method can be written as follows.Step 1.
We build block Krylov subspace:KrA-1E;A-1B.Step 2.
To implement Arnoldi process, specific algorithm is as follows:(a)
given matrix,A∈Rn, B∈Rn×p, and QR decomposition of matrix B, let B=V0T;(b)
computeV^1=AV0-V0H00, where H00=V0TAV0;(c)
computeV^2=AV1-V1H11-V0H01, where H11=V1TAV1, H01=V0TAV1;(d)
then compute the following parameter intrun:V^r=AVr-1-Vr-1Hr-1,r-1-V1H1,r-1-V0H0,r-1, where Hr-1,r-1=Vr-1TAVr-1,…,H1,r-1=V1TAVr-1, H0,r-1=V0TAVr-1.Step 3.
According to the transformation matrix, then get the original system of reduced order system:(21)E~dx~tdt=A~xt+B~uty~t=B~Tx~t,where x~t∈Rrp, E~=VTEV, A~=VTAV, B~=VTB.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.


## 5.2. 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:(22)Ψ0ζ=1,Ψ1ζ=ζ,Ψ2ζ=ζ2-1,Ψ3ζ=ζ3-3ζ,Ψ4ξ=ζ4-6ζ2+3,…Assuming the input of a linear variable, there are(23)E~r=E0+E1ζ1+E2ζ2,A~r=A0+A1ζ1+A2ζ2,B~r=B0+B1ζ1+B2ζ2,where E0,E1,E2,A0,A1,A2,B0,B1,B2 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):(24)E0E13E23E1E00E20E0dx~0(t)dtdx~1(t)dtdx~2(t)dt=A0A13A23A1A00A20A0x~(t)+B0B13B23B1B00B20B0u(t)y~t=B0B13B23B1B00B20B0Tx~t.Let(25)Esr=E0E13E23E1E00E20E0,Asr=A0A13A23A1A00A20A0,Bsr=B0B13B23B1B00B20B0.Equation (18) may be cast in the following “deterministic” form:(26)E~srdx~(t)dt=A~srx~(t)+B~sru(t)y~t=B~srTx~t.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.


## 5.3. Computation of the Final Results: Expected Values and Standard Deviations

Once we have solved then systems numerically, that is to say, when we have got the matrices A~sr and B~sr, 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(27)λ~kζ-≈λk·∑b=1mabk·ϕb=λk·∑b=1mabk·δ1b=λk·a1khhhhhhhh∀k∈1,…,n(28)U~k=Uk+∑r=1n∑l=1mβlr·ϕl·Ur=Uk+∑r=1nβlr·Ur.From (27) we get(29)λ~k2=λk·a1k2.Hence using the linearity and orthogonality properties as previously mentioned, we get(30)λ~k2=λk·a1k2+λk2·∑b=2mabk2·ϕb2.Finally using (28) and (29), we directly get the standard deviation:(31)δλ~k=λ~k2-λ~k2=λ·∑b=2mabk2·ϕb2.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)σU~kq=∑r=1n∑l=2mβlr2ϕl2·Urq2∑r=1r1,r2≠kn∑l=2mβlr1βlr2ϕl2·Ur1qUr2qhhhhh+2·∑r=1r1,r2≠kn∑l=2mβlr1βlr2ϕl2·Ur1qUr2q1/2.


## 6. 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. Figure1 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.Figure 1
Single-ended fixed electrostatic actuating microbeam structure.The beam can be modeled by coupling the Euler beam equation with the electrostatic force and the Reynolds squeeze-film damping equation as follows:(33)EI∂4y∂x4-s∂2y∂x2=-ξ0ωV22y2+∫0∞p-p0dy-ρ∂2y∂t2∇·y3p∇p=12μ1+6k∂py∂t.Among themy(x,t) is the immunity of the z direction, E is young’s modulus, I=wh3/12—moment of inertia’s is residual stress (Figure 5), ρ is the density, ε0 is dielectric constant, p is the air damping of the pressure, p0 is the ambient pressure, μ is the air viscosity coefficient, g is the distance from the plate under initial state, and λ is air mean free path (see Table 1).Table 1
The parameters of single-ended fixed electrostatic actuating microbeam structure.
l (beam length)
610
w (beam width)
40
h (beam thickness)
2.2
g
0 (initial gap distance)
2.3
E (material Young’s modulus)
149 GPa
ρ (material density)
2330 kg/m3
μ (air viscosity)
1.82 × 10−5 kg/(m·s)
λ (mean-free path of air)
0.064Lety~=y(x,t)-g/g, p~(x,t)=p(x,t)-p0/p0, substitute them into (3) and in g and p0 using Taylor series expansion, and take position, speed, and pressure as state variables like(34)x=yy˙p=x1x2x3;we can get the micro beam for original system state equation:(35)x˙=y˙y¨p˙=x2-EIρ·∂4yi∂x4+sρ∂2yi∂x2+ε0ωV2g3ρx1+p0gρ∫0ωx3dy-3ε0ωV22g3ρx12-ε0ωV22g3ρ1+6λg+1+6λgx3+2+6λgx1∗g2p012μ∇2x3-1+x3+x1.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 leti=i+1; go to Step 2.Step 4.
The process is over.We useN+1×M+1 mesh as shown in Figure 2, where N represents the number of inner grid points in the x direction and M is the number of inner grid points in the y 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 u(x,t) and p(x,y,t) 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:(36)∂2yi∂x2=yi+1-2yi+yi-1Δx2∂4yi∂x4=yi+2-4yi+1-6yi+4yi-1+yi+2Δx4∫0ωxi,yj,tdy=∑j=1j=MΔy2pxi-1,yj-1,t+pxi,yj,t.Accelerometer beam clamped at one end, so the boundary conditions are y0=0,y-1=y1; substitute (36) into (35) to discrete; then map the n state vector X=y1⋯yN∂yi/∂t⋯∂yN/tp11⋯pMNT to the various nodes; we can convert discrete equations into n dimensional state space equation (9), where n=2N+N∗M,(37)X˙t=AX(t)+Bv(t)yt=BTXT,where E=M00H∈Rn×n, A=NF-FT0∈Rn×n, and B=B10∈Rn×py(t) are the microdisplacement of the beam end and v(t) is the input voltage of the system.Figure 2
Structure model of microbeam.Figure 3
Finite element model of microbeam.Figure 4
The displacement of microbeam.Figure 5
The stress of microbeam.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 Figure6. The voltage load vector along with the mass and stiffness matrices were extracted as input random parameters.Figure 6
The endpoint deflection of the beam.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 thicknessh, length l, width w, Young’s modulus E, and the initial gap g 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.Table 2
The GPC response inputs parameters.
Variable
Distribution
Range
Units
h
h
N
(
0.05,3.40
)
1.6–3.6
μm
l
N
(
250,650
)
549–671
μm
w
N
(
6,32
)
27–33
μm
E
N
(
190.5,6.9
)
177–216
GPa
g
N
(
4.5,0.1
)
4.0–5.0
μmThen 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 Figure7.Figure 7
Smolyak grid (total number of points 29).Then, we compute 29 times of model reduction simulation at Smolyak point; the results were shown in Figure8. 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.Figure 8
Model reduction simulation results.We also compute that the resulting probability distribution of the beam maximum deflection is presented in Figures9, 10, and 11. In Figure 9, we just show geometric parameters thickness h, length l, and width w as random parameters, and Figure 10 shows Young’s module as random parameter, and Figure 11 shows the initial gap distance as random parameters.Figure 9
The geometric(w,h,l).Figure 10
Young’s modulus.Figure 11
The initial gap distance.


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



---

*Source: 125621-2015-04-05.xml*