精簡微機電元件模型之建立與動態反應之研究

行政院國家科學委員會補助專題研究計畫成果報告

※※※※※※※※※※※※※※※※※※※※※※※※※

※ ※

※

精簡微機電元件模型之建立與動態分析研究

※

※

※

※※※※※※※※※※※※※※※※※※※※※※※※※

計畫類別：█個別型計畫

□整合型計畫

計畫編號：NSC 90-2218-E-002-031

執行期間：

90 年 8 月 1 日至 91 年 7 月 31 日

計畫主持人：楊燿州

共同主持人：

本成果報告包括以下應繳交之附件：

□赴國外出差或研習心得報告一份

□赴大陸地區出差或研習心得報告一份

□出席國際學術會議心得報告及發表之論文各一份

□國際合作研究計畫國外研究報告書一份

執行單位：台大機械系

中

華

民

國

91 年

10

月

28

日

行政院國家科學委員會專題研究計畫成果報告

精簡微機電元件模型之建立與動態分析研究

Methodology of Compact Model Generation and Study of Dynamic

Response for MEMS Devices

計畫編號：NSC 90-2218-E-002-031

執行期限：90 年 8 月 1 日至 91 年 7 月 31 日

主持人：楊燿州 台大機械系

計畫參與人員：曾思遠 蔡宗勳 台大機械系

中文摘要

： 本篇文章的主旨是將動態模型，藉由有限差分 法(Finite Difference Method)計算出的解，建構 出來有效率且精確度高的精簡模型，用此精簡 模型分析動態系統，可用來分析動態方程式為 非線性且藕合的系統。精簡模型所建構的方式 是紀錄數次有限差分的結果(snapshot)，經由 Karhunen-Loeve Decomposition 求 出 精 簡 模型。此精簡模型可有效率且高精確的模擬出 動態系統暫態反應。並以微機電元件壓力感測 器(pressure sensor)作為驗證，精簡模型計算出 的結果，平均誤差<1%，計算的速度(speedup factor)可達到 5 以上，可驗證降階模型具有高 精確度與高效率。

關鍵詞： 模型降階演算法、微機電系

統、

Karhunen-Loeve Decomposition

Abstr act

:

In this work, we present a reduced order modeling technique for extracting compact macromodel from complex MEMS coupled systems. This methodology uses the snap-shots of finite-element/finite-difference simulation runs, and applies Karhunen-Loeve decomposition to extract reduced-order models. Comparison of computational efficiency as well as experimental verification are also provided.

一、 簡介

微機電系統(MEMS)通常整合數個不同 的領域能量[1][2]，如：動能(kinetic energy)、 彈 性 位 能 (elastic energy) 、 電 能 (electrostatic stored energy) 以 及 磁 能 (magnetostatic stored energy)，在多個能量維度的問題中，需要用有 限差分的方式求解聯立偏微分方程式，可是利 用有限差分法求解偏微分方程式的暫態解，需 要大量的計算，尤其在系統中含有多個元件更 是如此，如果可以不直接計算有限差分，而利 用降階的精簡模型計算，可以快速且精確的求 出系統的解，成為一個很重要的課題。 在文獻中，有利用線性化將原系統統馭 方程式轉換為線性方程式[3]，但因許多系統 本身為非線性，因此在線性化之後，不能得到 理想的結果，也有文獻是將原系統作二次展開 [4]，也因為會受到系統非線性程度的影響， 所得到的解與真實的解比較，誤差也是相當 大。 針對非線性系統無法單純用一線性系統 取 代 之 的 問 題 ， 於 是 衍 生 出 片 段 線 性 (piecewise linear)的方法[5][6]，此方法是在非 線性的系統中，針對不同狀態分別作線性化的 展 開 ， 再 將 線 性 化 的 系 統 作 Arnoldi algorithm[7]降階計算，此方法可以有效的模擬 出系統真實的解，但是在建構出片段線性的系 統 時 ， 需 要 先 輸 入 一 訓 練 輸 入 (training input)，再根據原系統所模擬的結果找出片段 線性的模型，此方法所得出之精簡模型與原訓 練輸入有很大的相關性，若改變環境的參數， 則精簡模型的輸出會有很大的誤差。 另一用來分析非線性系統的方式為 Karhunen-Loeve Decomposition [8][9][10]，此 方法是利用有限差分或有限元素法模擬數次 原系統，紀錄原系統的暫態解(snapshot)，再 將此結果使用 Karhunen-Loeve Decomposition 求出基底函數(basis function)，再利用這些基 底函數找出精簡模型，雖然在找出此精簡模型 之前，需要先使用有限差分或有限元素法模擬 原系統，但是只要經過一次這樣完整的模擬所 建構出的精簡模型是可以有高精確度、高效率 且具有很大的彈性的模型，此降階演算法是可 以用來分析單一元件，或是分析一整個系統， 都可以得到令人滿意的結果。 本 篇 論 文 是 採 用 Karhunen-Loeve

Decomposition 分析壓力感測器通電壓時產生 的 pulling-in 現象(如圖一)，壓力感測器的動 態方程式是整合桿的方程式(固力部分)與雷 諾方程式(流體部分)，為一非線性且藕合的系 統，分析的方式是先用有限差分法求解此聯立 偏微分方程式，並紀錄完整有限差分法不同時 間 的 高 度 與 壓 力 (snapshot) ， 進 而 利 用 Karhunen-Loeve Decomposition 求出基底函數 (basis function)進而建構出精簡模型，將原本 非線性且藕合的偏微分方程式，化簡為二個常 微分方程式的精簡模型，且精簡模型狀態的數 目遠少於原有限差分狀態的數目，此精簡模型 與有限差分法模擬結果顯示，精簡模型有著高 計算效率與高精確度的優點。 圖一 壓力感測器示意圖

二、 Kar hunen-Loeve Decomposition

介紹 Karhunen-Loeve Decomposition 之 前，先定義數個符號及函數，N 個任意的函 數 ， 我 們 稱 之 為 snapshots ， 符 號 為 {vn} ， n=1,2,… ,N ， 要 探 討 的 重 點 是 如 何 從 這 些 {vn}，找出最具代表性的的函數

φ

(

x

)

，幾個數 學的符號如下：

)

,

(

x

y

v

_n ：

a function defined in a function space(1)

{ }

v

n ： ensemble of snapshots (2)

∫

Ω

Ω

=

f

x

g

x

d

g

f

,

)

(

)

(

)

(

：

inner product in the function space (3)

∑

=

≡

N n n n

v

x

N

v

1

)

(

1

：

ensemble average of snapshots (4)

我們的目標為最大化

)

,

(

)

,

(

2

φ

φ

φ

λ

v

n

=

(5)， 利用(3)內積的定義，可將(5)式改為如下的型 式

∫

Ω

∫

Ω

=

(

)

(

)

(

'

)

(

'

)

'

)

,

(

φ

v

_n 2

φ

x

v

_n

x

dx

φ

x

v

_n

x

dx

{

v

n

(

x

)

v

n

(

x

'

)

φ

(

x

)

dx

}

φ

(

x

'

)

dx

'

∫ ∫

Ω Ω

=

(6) 接著定義函數

k

(

x

,

x

'

)

如下：

∑

=

=

=

N n T n n n n

v

x

v

x

N

x

v

x

v

x

x

k

1

)

(

)

(

1

)

'

(

)

(

)

'

,

(

(7) 以及線性運算子 R：

∫

Ω

•

≡

•

k

(

x

,

x

'

)

dx

'

R

(8) 因此可將(6)表示為

{ }{ }

∫

Ω

=

=

(

,

)

)

,

(

φ

v

_n 2

R

φ

φ

dx

R

φ

φ

(9) 因此最大化的問題(5)，可以轉變為特徵值 (eigenvalue)的問題：

λφ

φ

=

R

(10) 將最大化(5)的

λ

，變為求出(10)中最大的特徵 值，利用 Schmidt-Hilbert technique，將

φ

(

x

)

假 設如下：

)

(

)

(

x

v

_k

x

k k

∑

=

α

φ

(11) 將(11)帶入(10)可得：

∫ ∑

Ω

∑

= N n k k k T n n

x

v

x

v

x

dx

v

N

1

'

)

'

(

)

'

(

)

(

1

α

∑

=

n n n

v

(

x

)

α

λ

(12) 可將上式化簡為 n k nk

C

α

=

λα

(13)

∫

Ω

≡

1

v

(

x

'

)

v

(

x

'

)

dx

'

N

C

_nk T_{n k} (14) nk

C

是對稱(symmetric)且正定(positive definite) 的矩陣，由(13)特徵值問題所求出的特徵向量 (eigenvector)，可代入(11)中找出

φ

_k，特徵值 依 大 小 排列

λ

₁

>

λ

₂

>

...

>

λ

_n，分別對應了 N

φ

φ

φ

₁

,

₂

,...,

，

φ

₁對應了最大的特徵值，也代 表了{vn}中最主要的結構，

φ

₂對應第二大的特 徵值，也代表了{vn}次主要的結構，依此類 推，因此可以將非線性系統作幾次完整的計 算 ， 紀 錄 不同 時 間的 結 果(snapshot)，利用 Karhunen-Loeve Decompositon 轉換為較精簡 的模型，以便我們可以快速且精確的計算出系 統的解。

三、壓力感測器分析

壓力感測器(圖一)，是一微結構之彈性 桿(elastic beam)，當施以一電壓，此彈性樑會 因為靜電力的作用而吸附(pull in)，在分析吸 附過程中，此桿的暫態反應非常敏感於桿下方 的壓力[11][12]，此吸附的現象可用一維的桿

方程式與二維可壓縮絕熱雷諾方程式聯立分 析[13]， 2 2 2 2 4 4

t

z

F

F

x

z

S

x

z

EI

_{elec air}

∂

∂

−

+

=

∂

∂

−

∂

∂

_ρ

(15)

t

pz

p

p

z

K

∂

∂

=

∇

+

⋅

∇

((

1

6

)

3

)

12

η

(

)

(16)

)

2

(

2 2 0

z

V

F

_elec

=

−

ε

ω

此為靜電力，

∫

−

=

w air

p

p

dy

F

0

(

0

)

此為桿下方與上方壓力 差所產生的負載，z(x,t)為桿的高度，p(x,y,,t) 為 桿 下 方 的 空 氣 壓 力 ， Knudsen’s number

z

t

x

K

(

,

)

=

λ

/

，空氣的平均自由路徑(mean free path)

λ

=

0

.

064

µ

m

， 楊 氏 彈 性 模 數

Gpa

E

=

149

，桿的轉動慣量 3

/

12

wh

I

=

， 桿 之 寬 度

w

=

40

µ

m

， 桿 之 長 度

m

l

=

610

µ

，桿之厚度

m

h

=

2

.

2

µ

，桿之原始高度

z

₀

=

2

.

3

µ

m

，內 應 力

S

/

hw

=

−

3

.

7

Mpa

， 密 度 為 3

/

2330

]

/

[

ρ

hw

=

kg

m

， 空 氣 黏 滯 係 數 5

10

82

.

1

×

−

=

η

。此桿是假設為兩端固定，壓 力沿著桿的方向是開放的，在桿的兩個尾端是 封閉的(沒有氣流)。

四、精簡模型的推導

首先我們假設偏微分方程式為以下的型 式：

f

u

L

(

)

=

(17) L為一微分運算子，且L允許為非線性的運算 子，u為此偏微分方程式的解，接著假設解是 u(x,t)，外力為f(x,t)，根據 KL Decomposition 可 得知解可以寫成以下的型式：

∑

=

=

N i i i

t

a

x

t

x

u

1

)

(

)

(

)

,

(

ˆ

α

(18) 引入 Galerkin condition： 0 ) ) ˆ ( ( ) ) ˆ ( , (

a

i

L

u

−

f

=

∫

a

iT

L

u

−

f

dx

= (19) 對於所有的

i

∈

{

1

,

2

,

3

,...,

N

}

。 接著寫出壓力感測器 pull-in 解的型式：

∑

=

+

=

M i i i

t

b

x

z

t

x

z

1 0

(

)

(

)

)

,

(

ˆ

β

(20)

∑

=

+

=

N i i i

t

a

x

p

t

y

x

p

1 0

(

)

(

)

)

,

,

(

ˆ

α

(21) 0

z

為桿未變形的電壓，

p

0為一大氣壓， 利用 Galerkin condition 可將原偏微分方程式 化為兩聯立的常微分方程式， 桿的偏微分方程式可推導為如下的常微分方 程式：

0

=

+

+

K

f

M

β&&

β

(22)

∫

=

L i j ij

b

b

dx

M

ρ

dx

x

b

x

b

S

x

b

x

b

EI

K

L j i j i ij

∫

_∂

∂

∂

∂

+

∂

∂

∂

∂

=

(

₂

)

2 2 2

∫

+

−

=

L i elec air i

b

F

F

dx

f

(

)

L是桿的長度。 雷諾方程式可以推導如下： 0 = + +

B

C

A

α

&

α

(23)

dxdy

z

a

a

A

_ij

∫

_{i j} Ω

=

12

η

ˆ

∫

Ω









∂

∂

∂

∂

+

∂

∂

∂

∂

+

=

y

a

y

a

x

a

x

a

p

z

K

B

i j i j ij

{(

1

6

)

ˆ

ˆ

3

x

z

a

a

_{i j}

∂

∂

+

12

η

ˆ

}

dxdy

dxdy

x

z

a

p

c

_i

∫

_i Ω

∂

∂

=

12

η

₀

ˆ

Ω

為桿的面積。 而

a

_i

(

x

)

與

b

_i

(

x

)

(basis function)所選取 的方式是利用有限差分法模擬原系統，分別模 擬不同電壓 9v, 10v, 12v, 16v，以固定時間間 隔紀錄桿的高度與壓力(snapshot)，再將這些 snapshots 作 SVD(singular value

decomposition)[14]， T

W

V

U

=

Σ

(24)

V

的行向量[v1,v2,v3,… ]，分別對應了基底函 數，由最具代表性的基底、次代表性的基底， 依此類推。

五、結果與討論

分析壓力感測器的方式是先將桿作網格 化(mesh)，將桿面分為41×20個格點，再利 用有限差分法，將(15)(16)對 x 與 y 偏微分的 部分改為差分，對時間積分採用的是 Runge-Kutta 的方法來計算，桿的方程式是一 維，壓力的方程式是二維，用有限差分法共需 要解 ODE 數目為 41

×

2+820=902 個。 精簡模型模擬的結果與有限差分法模擬 的結果比較，如圖二，M 代表高度的基底函 數個數，N 代表壓力的基底函數個數，模擬的 電壓為 10v，精簡模型相對於有限差分法的誤 差，如圖三。

有限差分與KL Decomposition比較 0 0.5 1 1.5 2 2.5

0.0E+00 4.0E-05 8.0E-05 1.2E-04 1.6E-04 2.0E-04 時間 (sec) 桿之高度(um) 有限差分法 KL M=1, N=4 KL M=1, N=2 KL M=1, N=3 圖二 桿的最低點對時間的變化圖，輸入為 10v。 誤差 vs 時間 0.00001 0.0001 0.001 0.01 0.1 1

0.0E+00 5.0E-05 1.0E-04 1.5E-04 2.0E-04

時間 (sec) 誤差 M=1, N=4 M=1, N=3 圖三 誤差對時間的變化圖，電壓 10v 由圖二可知，當基底函數達到 M=1, N=3，精簡模型所繪出曲線與有限差分的曲線 相當接近，由圖三可知，基底函數在 M=1, N=4 時，平均誤差在 1%以下。

time elapsed (sec) speedup factor Full mesh 394 1 M=1, N=4 78 5.1 M=1, N=3 22 17.9 M=1, N=2 4 98.5 表一 精簡模型與有限差分模擬時間比較 如表一，可知精簡模型計算的時間較完 整有限差分法的時間大幅減小，可驗證精簡模 型為有效率，且高精確度的模型。

六、參考文獻

[1] Lynn D. Gabbay, Jan E. Mehner, and Stephen D. Senturia, “Computer-Aided Generation of Nonlinear

Reduced-Order Dynamic Macromodels-I,”Journal of MICROELECTROMECHANICAL SYSTEMS, Vol. 9, NO. 2, June 2000

[2] Jan E. Mehner, Lynn D. Gabbay, and Stephen D. Senturia, “Computer-Aided Generation of Nonlinear Reduced-Order Dynamic Macromodels-II,” Journal of MICROELECTROMECHANICAL SYSTEMS, Vol. 9, NO. 2, June 2000

[3] Jinghong Chen and Sung-Mo (Steve) Kang, “An Algorithm for Automatic Model-Order Reduction of Nonlinear MEMS Devices,”IEEE International Symposium on Circuit and Systems, May 28-31 [4] Y. Chen, J. White, “A Quadratic Method for Nonlinear Model Order Reduction,”in proceedings of the International Conference on Modeling and Simulation of Microsystems, pp. 477-480, 2000. [5] M. Rewienski, J. White, “A Trajectory Piecewise-Linear Approach to Model Order Reduction and Fast Simulation of Nonlinear Circuit and Micromachined Devices,”in the proceedings of the International Conference on Computer-Aided Design, pp.252-257, 2001.

[6] M. Rewienski, J. White, “Improving Trajectory Piecewise-Linear Approach to Nonlinear Model Order Reduction for Micromachined Devices Using an Aggregated Projection Basis,”Modeling and Simulation of Microsystems, 2002.

[7] F. Wang, J. White, “Automatic model order reduction of a microdevice using the Arnoldi approach,” pp. 527-530, ASME, 1998.

[8] H. M. Park, M. W. Lee, “An Efficient Method of Solving the Navier-Stokes Equations for Flow Control,” International Journal forNumerical Methods in Engineering, Vol. 41, 1133-1151, 1998. [9] Elmer S. Hung, Stephen D. Senturia, “Generation Efficient Dynamical Models for

Microelectromechanical Systems from a Few Finite-Element Simulation Runs,”Journal of

Microelectromechanical Systims, Vol. 8, No. 3, 1999. [10] Jinghong Chen, Sung-Mo Kang, “Model-Order Reduction of Nonlinear MEMS Devices through Arclength-Based Karhunen-Loeve Decomposition,”

Circuit and Systems, 2001.

[11] C. –L. Chen, J. J. Yao, “Damping control of MEMS Devices Using structural design approach,”

IEEE Solid-State Sensor and Actuator Workshop, pp.72-75. 1992.

[12] Y. –J. Yang, M. –A. Gretillat, and S. D. Senturia, “Effect of air damping on the dynamics of

nonuniform deformations of microstructures,” in Transducers, pp. 1093-1096, ’97.

[13] B. J. Hamrock, Fundamentals of Fluid Film Lubrication. New York: McGraw-Hill, 1994. [14] Gene H. Golub, Charles F. Van Loan, Matrix Computations. Th