CN101551835B - 微机电系统的角度参数化宏建模方法 - Google Patents

Abstract

本发明公开了一种微机电系统的角度参数化宏建模方法，属于微机电系统设计与模型降阶领域。该方法首先从相同几何形状的微结构中任意选取一个作为参考微结构，采用数值降阶方法提取其宏模型；然后根据微结构方位之间的关系，采用矩阵坐标变换方法建立其它方位微结构相对参考微结构的力与位移之间的转换方程；最后，将力与位移之间的转换方程代入宏模型的行为方程中，实现MEMS的角度参数化宏建模。本发明采用矩阵坐标变换方法实现角度参数化的宏建模，通过设置相应的角度参数，即可避免结构完全相同，初始方位不同的MEMS微结构的重复宏建模过程，从而使基于宏模型的MEMS系统级建模与仿真速度加快，进而提高整个MEMS的设计效率。

Description

微机电系统的角度参数化宏建模方法

技术领域

本发明涉及一种微机电系统(MEMS)的角度参数化宏建模方法，属于微机电系统设计与模型降阶领域。

背景技术

随着MEMS设计技术及工艺水平的迅速发展，MEMS器件结构也变得越来越复杂，大量新型器件相继产生。然而，结构对称性依然是发展中MEMS器件的一大特点，即：在MEMS的设计中，大量存在几何形状完全相同，而初始方位不同的功能部件或微结构。对这些微结构的有限元模型进行自由度缩减，建立能够表征微结构输入输出行为的低阶模型的过程称为宏建模过程，所提取的低阶模型也称为宏模型。现有的宏建模方法由于无法实现角度参数化，必须对这些结构分别处理，各自独立地进行宏建模，而宏建模本身需要大量耗时的数值计算，从而导致这种方法计算代价高，效率低。因此，实现这些结构的角度参数化宏模型对于避免重复的宏建模过程，提高整个MEMS的设计效率具有重要意义。

Carnegie Mellon大学的Qi Jing在其博士论文“Modeling and Simulation for Design ofSuspended MEMS”中，针对MEMS集总参数模型，采用坐标变换方法实现了角度参数化，但其集总参数模型采用矩阵结构分析和静电场等经典理论方法，经过解析推导而提取，与本发明将角度参数化应用于数值降阶的宏建模方法不同。此外，近年来，参数化宏建模方法取得较大进展，德国IMTEK的E.B.Rudnyi教授，美国MIT的Luca Daniel博士以及UC Davis的Bai Zhaojun教授等人对参数化宏建模方法进行了广泛而深入的研究，然而其主要针对的是MEMS结构的几何尺寸及材料属性参数化，对于角度参数化宏建模方法方面的研究尚未涉及。

发明内容

本发明的目的在于提出一种微机电系统的角度参数化宏建模方法，对于MEMS设计中几何形状完全相同，初始方位不同的微结构，只需完成其中一个微结构的宏建模，通过设置相应的角度参数，即可立即实现其它任意方位微结构的宏建模。

本发明提出的微机电系统的角度参数化宏建模方法首先从相同几何形状的微结构中任意选取一个作为参考微结构，采用数值降阶方法提取其宏模型；然后根据微结构方位之间的关系，采用矩阵坐标变换方法建立其它方位微结构相对参考微结构的力与位移之间的转换方程；最后，将力与位移之间的转换方程代入宏模型的行为方程中，实现MEMS的角度参数化宏建模。该方法主要步骤如下：

步骤一：进行参考微结构的宏建模，对于几何形状完全相同，初始方位不同的微结构，任意选择其中一个部件作为参考微结构，采用三维造型工具建立微结构的三维实体模型，并使用有限元方法进行模态分析，从分析结果中提取参考微结构的质量、刚度矩阵，建立微结构的二阶动力学行为方程：

$\left\{\begin{array}{c}M\stackrel{\·\·}{q}\left(t\right)+D\stackrel{\·}{q}\left(t\right)+\mathrm{Kq}\left(t\right)={\mathrm{PF}}_c\\ u_c=E^{T}q\left(t\right)\end{array}\right.---\left(1\right)$

式中：t为时间，q(t)为节点位移向量，P为载荷矩阵，F_c为输入载荷向量，u_c为输出位移向量，E为输出矩阵，M，D，K∈R^N×N分别为参考微结构的质量、阻尼和刚度矩阵，MEMS结构一般采用Rayleigh结构阻尼模型，即D＝αM+βK，其中α、β是和结构模态阻尼比相关的参数。

上述二阶动力学行为方程经等价变换后，可转换为一阶状态方程：

$\left\{\begin{array}{c}C\stackrel{\·}{x}\left(t\right)+\mathrm{Gx}\left(t\right)={\mathrm{BF}}_c\\ u_c=L^{T}x\left(t\right)\end{array}\right.---\left(2\right)$

其中： $x\left(t\right)=\left[\begin{array}{c}q\left(t\right)\\ \stackrel{\·}{q}\left(t\right)\end{array}\right],$ $C=\left[\begin{array}{cc}D& M\\ W& 0\end{array}\right],$ $G=\left[\begin{array}{cc}K& 0\\ 0& -W\end{array}\right],$ $B=\left[\begin{array}{c}P\\ 0\end{array}\right],$ $L=\left[\begin{array}{c}E\\ 0\end{array}\right].$

W是一个N×N的非奇异矩阵，在结构力学中，通常取W＝M，这样上述状态方程的系统矩阵C、G保持了原二阶状态方程的系数矩阵的对称特性。从上述变换过程可知，二阶动力学行为方程变换为一阶状态方程后，其方程规模也由N增加到2N。

以A＝-G^-1C、R＝G^-1B作为输入，应用Arnoldi算法对状态方程(2)进行降阶，得到参考微结构的宏模型：

$\left\{\begin{array}{c}{C}_r{\stackrel{\·}{x}}_r\left(t\right)+G_r{x}_r\left(t\right)=B_r{F}_c\\ u_c=L_r^T{x}_r\left(t\right)\end{array}\right.---\left(3\right)$

其中：t是时间变量，x_r(t)是宏模型的状态变量，C_r、G_r为宏模型的系统矩阵，B_r是宏模型的输入矩阵，L_r是宏模型的输出矩阵，r为宏模型的规模，由于r＜＜N，所以宏模型的求解效率相对原始状态方程大幅提高。

步骤二：根据微结构方位之间的关系，采用矩阵坐标变换方法建立其它方位微结构相对参考微结构的力与位移转换方程。参阅图1，本发明首先定义两个坐标系以描述微结构在三维空间内的运动。把与参考微结构初始方位重合的坐标系称为总体坐标系，用c-x_cy_cz_c表示，简称c坐标系。为了便于描述微结构的转动行为，把固连于参考微结构某点p的坐标系称为局部坐标系，用p-x_py_pz_p表示，简称p坐标系。根据定义可知，c坐标系是静止坐标系，而p坐标系随部件运动。这样，通过p坐标系的转动可描述微结构的转动行为。本发明采用以下转动顺序引入确定p坐标系在c坐标系的方位角，即p坐标系绕c坐标系x、y、z坐标轴的相对转角，也即其它微结构相对参考微结构的转角，用(α，β，γ)表示：最初两坐标系完全重合，而后顺序通过三次转动达到p坐标系所在方位，这三次转动为：绕z_c转一个γ角，使p-x_py_pz_p由最初与c-x_cy_cz_c重合的位置转到p-x_p1y_p1z_c，第二次绕y_p1轴转β角，使p-x_p1y_p1z_c到达新的p-x_p2y_p1z_p2位置，第三次绕x_p2转α角，使p-x_p2y_p1z_p2到达p-x_py_pz_p的最终位置，把这种转动方式记为 $c\underset{z_c}{\overset{\mathrm{\γ}}{\→}}1\underset{y_1}{\overset{\mathrm{\β}}{\→}}2\underset{x_2}{\overset{\mathrm{\α}}{\→}}p,$ 其转动过程参阅图2。

设T_c ^p为变换过程中c坐标系到p坐标系的转换矩阵，则根据矩阵坐标变换原理，可知：

$T_c^p=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \mathrm{\α}& \sin \mathrm{\α}\\ 0& -\sin \mathrm{\α}& \cos \mathrm{\α}\end{array}\right]\left[\begin{array}{ccc}\cos \mathrm{\β}& 0& -\sin \mathrm{\β}\\ 0& 1& 0\\ \sin \mathrm{\β}& 0& \cos \mathrm{\β}\end{array}\right]\left[\begin{array}{ccc}\cos \mathrm{\γ}& \sin \mathrm{\γ}& 0\\ -\sin \mathrm{\γ}& \cos \mathrm{\γ}& 0\\ 0& 0& 1\end{array}\right]=$

$\left[\begin{array}{ccc}\cos \mathrm{\γ}\cos \mathrm{\β}& \sin \mathrm{\γ}\cos \mathrm{\β}& -\sin \mathrm{\β}\\ -\sin \mathrm{\γ}\cos \mathrm{\α}+\cos \mathrm{\γ}\sin \mathrm{\β}\sin \mathrm{\α}& \cos \mathrm{\γ}\cos \mathrm{\α}+\sin \mathrm{\γ}\sin \mathrm{\β}\sin \mathrm{\α}& \cos \mathrm{\β}\sin \mathrm{\α}\\ \sin \mathrm{\γ}\sin \mathrm{\α}+\cos \mathrm{\γ}\sin \mathrm{\β}\cos \mathrm{\α}& -\cos \mathrm{\γ}\sin \mathrm{\α}+\sin \mathrm{\γ}\sin \mathrm{\β}\cos \mathrm{\α}& \cos \mathrm{\β}\cos \mathrm{\α}\end{array}\right]$

设与p坐标系固连的微结构在c坐标系的位移与转角为u_cp＝[u_x，u_y，u_z，a_x，a_y，a_z]^T，微结构所受力及转矩为F_cp＝[F_x，F_y，F_z，T_x，T_y，T_z]^T。则可得其它方位微结构相对参考微结构的力与位移转换方程：

$u_{\mathrm{cp}}=\left[\begin{array}{c}{T}_c^p\\ T_c^p\end{array}\right]u_c,F_{\mathrm{cp}}=\left[\begin{array}{c}{T}_c^p\\ T_c^p\end{array}\right]F_c---\left(4\right)$

其中u_x，u_y，u_z，a_x，a_y，a_z分别为沿x、y、z方向的平动位移和绕x、y、z轴的转角，F_x，F_y，F_z，T_x，T_y，T_z分别为沿x、y、z方向的力和绕x、y、z轴的转矩。u_c为参考微结构在c坐标系中的位移与转角，F_c为参考微结构在c坐标系中的力与转矩，式(4)表明，微结构在c坐标系与p坐标系中的动力学关系可通过矩阵坐标变换来表征。

步骤三：将力与位移之间的转换方程代入宏模型的行为方程中，实现宏模型的角度参数化。将步骤二所得力与位移之间的转换方程(4)代入到步骤一所得的参考微结构的宏模型(3)中，同时根据其它微结构相对参考微结构的方位，设置相对转角参数值(α，β，γ)，迅速提取其它方位微结构的宏模型，实现角度参数化的宏建模。

本发明的有益效果是：采用矩阵坐标变换方法实现角度参数化的宏建模，通过设置相应的角度参数，即可避免结构完全相同，初始方位不同的MEMS微结构的重复宏建模过程，从而使基于宏模型的MEMS系统级建模与仿真速度加快，进而提高整个MEMS的设计效率。

附图说明

图1c坐标系与p坐标系的定义示意图

图2c坐标系到p坐标系的转换过程示意图

图3本发明具体实施例中平面加速度计的结构示意图

图中，1-左上角折叠梁，2-右上角折叠梁，3-左下角折叠梁，4-右下角折叠梁，5-质量块

具体实施例

下面结合一种平面加速度计的折叠梁角度参数化的宏建模过程对本发明进行进一步说明。参照附图3，加速度计的框架结构由左上角折叠梁1、右上角折叠梁2、左下角折叠梁3、右下角折叠梁4和质量块5构成，它关于坐标系oxyz的原点o中心对称以及x、y坐标轴的轴对称性，在加速度计的左上、右上、左下、右下方均有一个结构复杂的折叠梁。折叠梁的角度参数化宏建模过程具体步骤如下：

步骤一：选择左上角折叠梁1作为参考结构，采用有限元软件ANSYS，选择beam188单元进行模态分析，其有限元模型含有356个节点，其中包含两个输入输出节点，由于beam188的每个节点包含6个自由度，即沿x、y、z坐标轴的平动和绕x、y、z坐标轴的转动，因此左上角折叠梁1共包含12个输入输出端口，模型的质量刚度矩阵的规模为N＝356×6＝2136阶。提取有限元模型的质量、刚度矩阵，并令α＝0，β＝6×10^-6，可得左上角折叠梁1的二阶动力学行为方程为：

$\left\{\begin{array}{c}M\stackrel{\·\·}{q}\left(t\right)+6\×{10}^{-6}K\stackrel{\·}{q}\left(t\right)+\mathrm{Kq}\left(t\right)={\mathrm{PF}}_{\mathrm{cp}1}\left(t\right)\\ u_{\mathrm{cp}1}=E^{T}q\left(t\right)\end{array}\right.---\left(5\right)$

其中，t为时间，q(t)为节点位移向量，P∈R^2136×12是输入矩阵，M，K∈R^2136×2136分别是质量、刚度阵，F_cp1＝[F_x1，F_y1，F_z1，T_x1，T_y1，T_z1，F_x2，F_y2，F_z2，T_x2，T_y2，T_z2]^T∈R¹²是输入载荷向量，表示作用于左上角折叠梁1的全部12个端口上的力与转矩，E∈R^2136×12为输出矩阵，u_cp1＝[u_x1，u_y1，u_z1，a_x1，a_y1，a_z1，u_x2，u_y2，u_z2，a_x2，a_y2，a_z2]^T∈R¹²是输出位移向量，表示左上角折叠梁1的全部12个端口上的输出位移及转角。某些情况下还可能包括部件的速度向量，取决于所关注的物理量。将方程(6)转换为一阶状态方程如下：

$\left\{\begin{array}{c}C\stackrel{\·}{x}\left(t\right)+\mathrm{Gx}\left(t\right)={\mathrm{BF}}_{\mathrm{cp}1}\left(t\right)\\ u_{\mathrm{cp}1}\left(t\right)=L^{T}x\left(t\right)\end{array}\right.---\left(6\right)$

以A＝-G^-1C、R＝G^-1B作为输入，应用Arnoldi算法对状态方程(6)进行降阶，并取宏模型的阶数为4，得到左上角折叠梁1的宏模型：

$\left\{\begin{array}{c}{C}_r{\stackrel{\·}{x}}_r\left(t\right)+G_r{x}_r\left(t\right)=B_r{F}_{\mathrm{cp}1}\left(t\right)\\ u_{\mathrm{cp}1}\left(t\right)=L_r^T{x}_r\left(t\right)\end{array}\right.---\left(7\right)$

其中：t是时间变量，x_r(t)∈R⁵是宏模型的状态变量，C_r，G_r∈R^5×5是宏模型的系统矩阵，B_r∈R^4×12是输入矩阵，L_r∈R^4×12是输出矩阵。

步骤二：定义oxyz为c坐标系，参阅图3，它是一个静止坐标系。由于参考微结构所在p坐标系和c坐标系重合，也即左上角折叠梁1相对其本身的方位角为(0，0，0)。而右上角折叠梁2、左下角折叠梁3和右下角折叠梁4所在方位相对左上角折叠梁1的转角分别下表所示：

折叠梁方位	α	β	γ
				右上角	0	180	0
左下角	180	0	0
				右下角	0	0	180

因此它们所在的p坐标系相对oxyz坐标系的转换矩阵分别如下：

$T_c^{p2}=\left[\begin{array}{ccc}-1& 0& 0\\ 0& 1& 0\\ 0& 0& -1\end{array}\right],$ $T_c^{p3}=\left[\begin{array}{ccc}1& 0& 0\\ 0& -1& 0\\ 0& 0& -1\end{array}\right],$ $T_c^{p4}=\left[\begin{array}{ccc}-1& 0& 0\\ 0& -1& 0\\ 0& 0& 1\end{array}\right]$

其中T_c ^p2、T_c ^p3、T_c ^p4分别表示右上角折叠梁2、左下角折叠梁3和右下角折叠梁4所在的p坐标系相对oxyz坐标系的转换矩阵，将它们代入方程(4)中，可得右上角折叠梁2、左下角折叠梁3和右下角折叠梁4相对左上角折叠梁1的力与位移之间的转换方程分别为：

$u_{\mathrm{cp}2}=\left[\begin{array}{cccccc}-1& & & & & \\ & 1& & & 0& \\ & & -1& & & \\ & & & -1& & \\ & 0& & & 1& \\ & & & & & -1\end{array}\right]u_{\mathrm{cp}1},$ $F_{\mathrm{cp}2}=\left[\begin{array}{cccccc}-1& & & & & \\ & 1& & & 0& \\ & & -1& & & \\ & & & -1& & \\ & 0& & & 1& \\ & & & & & -1\end{array}\right]F_{\mathrm{cp}1}---\left(8\right)$

$u_{\mathrm{cp}3}=\left[\begin{array}{cccccc}1& & & & & \\ & -1& & & 0& \\ & & -1& & & \\ & & & 1& & \\ & 0& & & -1& \\ & & & & & -1\end{array}\right]u_{\mathrm{cp}1},$ $F_{\mathrm{cp}3}=\left[\begin{array}{cccccc}1& & & & & \\ & -1& & & 0& \\ & & -1& & & \\ & & & 1& & \\ & 0& & & -1& \\ & & & & & -1\end{array}\right]F_{\mathrm{cp}1}---\left(9\right)$

$u_{\mathrm{cp}4}=\left[\begin{array}{cccccc}-1& & & & & \\ & -1& & & 0& \\ & & 1& & & \\ & & & -1& & \\ & 0& & & -1& \\ & & & & & 1\end{array}\right]u_{\mathrm{cp}1},$ $F_{\mathrm{cp}4}=\left[\begin{array}{cccccc}-1& & & & & \\ & -1& & & 0& \\ & & 1& & & \\ & & & -1& & \\ & 0& & & -1& \\ & & & & & 1\end{array}\right]F_{\mathrm{cp}1}---\left(10\right)$

步骤三：将步骤二所得力与位移的转换方程(8)、(9)、(10)依次代入方程(7)输入载荷向量和输出位移向量中，即可获取右上角折叠梁2、左下角折叠梁3和右下角折叠梁4的宏模型行为方程，进而实现角度参数化的宏建模。

这样，本实例首先选取左上角折叠梁为参考部件并进行宏建模，然后通过设置相应的角度参数，以提取其它方位折叠梁的宏模型，从而实现了角度参数化的宏建模。

Claims (1)

1.一种微机电系统的角度参数化宏建模方法，其特征在于，包括如下步骤：

步骤一：进行参考微结构的宏建模，对于几何形状完全相同，初始方位不同的微结构，任意选择其中一个部件作为参考微结构，建立其宏模型：

$\left\{\begin{array}{c}{C}_r{\stackrel{\·}{x}}_r\left(t\right)+G_r{x}_r\left(t\right)=B_r{F}_c\\ u_c=L_r^T{x}_r\left(t\right)\end{array}\right.$

其中：t是时间变量，x_r(t)是宏模型的状态变量，C_r、G_r为宏模型的系统矩阵，B_r是宏模型的输入矩阵，F_c为输入载荷向量，u_c为输出位移向量，L_r是宏模型的输出矩阵；

步骤二：根据微结构方位之间的关系，采用矩阵坐标变换方法建立其它方位微结构相对参考微结构的力与位移转换方程，首先定义两个坐标系以描述微结构在三维空间内的运动，把与参考微结构初始方位重合的坐标系称为总体坐标系，用c-x_cy_cz_c表示，简称c坐标系，把固连于参考微结构某点p的坐标系称为局部坐标系，用p-x_py_pz_p表示，简称p坐标系，采用以下转动顺序引入确定p坐标系在c坐标系的方位角，即p坐标系绕c坐标系x、y、z坐标轴的相对转角，也即其它微结构相对参考微结构的转角，用(α，β，γ)表示：最初两坐标系完全重合，而后顺序通过三次转动达到p坐标系所在方位，这三次转动为：绕z_c转一个γ角，使p-x_py_pz_p由最初与c-x_cy_cz_c重合的位置转到p-x_p1y_p1z_c，第二次绕y_p1轴转β角，使p-x_p1y_p1z_c到达新的p-x_p2y_p1z_p2位置，第三次绕x_p2转α角，使p-x_p2y_p1z_p2到达p-x_py_pz_p的最终位置，把这种转动方式记为 $c\underset{z_c}{\overset{\mathrm{\γ}}{\→}}1\underset{y_1}{\overset{\mathrm{\β}}{\→}}2\underset{x_2}{\overset{\mathrm{\α}}{\→}}p;$

变换过程中c坐标系到p坐标系的转换矩阵T_c ^p为：

$T_c^p=\left[\begin{array}{ccc}\cos \mathrm{\γ}\cos \mathrm{\β}& \sin \mathrm{\γ}\cos \mathrm{\β}& -\sin \mathrm{\β}\\ -\sin \mathrm{\γ}\cos \mathrm{\α}+\cos \mathrm{\γ}\sin \mathrm{\β}\sin \mathrm{\α}& \cos \mathrm{\γ}\cos \mathrm{\α}+\sin \mathrm{\γ}\sin \mathrm{\β}\sin \mathrm{\α}& \cos \mathrm{\β}\sin \mathrm{\α}\\ \sin \mathrm{\γ}\sin \mathrm{\α}+\cos \mathrm{\γ}\sin \mathrm{\β}\cos \mathrm{\α}& -\cos \mathrm{\γ}\sin \mathrm{\α}+\sin \mathrm{\γ}\sin \mathrm{\β}\cos \mathrm{\α}& \cos \mathrm{\β}\cos \mathrm{\α}\end{array}\right]$

其它方位微结构相对参考微结构的力与位移转换方程为：

$u_{\mathrm{cp}}=\left[\begin{array}{c}{T}_c^p\\ T_c^p\end{array}\right]u_c,$ $F_{\mathrm{cp}}=\left[\begin{array}{c}{T}_c^p\\ T_c^p\end{array}\right]F_c---\left(1\right)$

式中u_cp＝[u_x，u_y，u_z，a_x，a_y，a_z]^T，表示与p坐标系固连的微结构在c坐标系的位移与转角，F_cp＝[F_x，F_y，F_z，T_x，T_y，T_z]^T，表示微结构所受力及转矩，u_x，u_y，u_z，a_x，a_y，a_z分别为沿x、y、z方向的平动位移和绕x、y、z轴的转角，F_x，F_y，F_z，T_x，T_y，T_z分别为沿x、y、z方向的力和绕x、y、z轴的转矩，u_c为参考微结构在c坐标系中的位移与转角，F_c为参考微结构在c坐标系中的力与转矩；

步骤三：将步骤二所得力与位移之间的转换方程(1)代入到步骤一所得的参考微结构的宏模型中，同时根据其它微结构相对参考微结构的方位，设置相对转角参数值(α，β，γ)，代入方程(1)中，提取其它方位微结构的宏模型，实现角度参数化的宏建模。

Images