MEMS系统级设计中的一种变截面梁建模方法

摘要

本发明公开了一种微机电系统(MEMS)系统级设计中的变截面梁建模方法，属于微机电系统设计领域。该方法的主要过程为：首先建立变截面梁的局部坐标系，然后根据结构力学和结构矩阵分析理论分别推导出局部坐标系下变截面梁刚度矩阵、阻尼矩阵和质量矩阵，并利用这些矩阵建立变截面梁的二阶动力学方程，最后利用硬件描述语言对二阶动力学方程进行编码，实现变截面梁的参数化系统级组件建模。本发明提出的针对宽度随长度方向线性变化的变截面梁的系统级建模方法，解决了当前对具有此种结构的MEMS器件无法进行系统级仿真设计的问题。

说明书

技术领域

本发明涉及微机电系统(MEMS)系统级设计中的一种变截面梁建模方法，属于微机电系统设计领域。

背景技术

在MEMS结构中，为了实现梁的某些设计目标，如应力集中小、体积小、强度高、重量轻、稳定性佳等，并且考虑梁所受弯矩的大小沿轴线变化时，常采用宽度随长度线性变化的变截面梁代替常规的等截面梁来优化结构。如在微加速度计、微机械陀螺等惯性器件中，绝大多数悬臂梁所承受的弯矩都是沿其轴向变化的，因此必须保证悬臂梁结构自由端具有较小的刚度来产生较大的形变，以使敏感质量等惯性结构具有较大的运动幅度；同时还要考虑到悬臂梁与锚点相连的固定端处必须满足强度的要求。综合以上梁的特点，在某些特定的情况下，设计此种变截面梁结构对于MEMS器件的设计是必须的。

针对MEMS结构和器件的整体行为的快速建模与仿真要求，西北工业大学的霍鹏飞在其论文“基于组件网络方法的微加速度计建模与仿真”中，提出了基于多端口组件网络MuPEN(multi-port-element network)的MEMS系统级建模和设计方法，并建立了参数化系统级组件库。该方法把整个MEMS器件分解为多个功能结构部件，各个功能结构部件建立对应的多端口参数化系统级组件，通过这些参数化系统级组件的相互连接形成网络表述整个MEMS器件，并进行仿真分析。另外，根据相似的方法创建的组件库还有Carnegie Mellon大学的NODAS，Coventor ware^TM的Architect等。

当前，国内外MEMS设计软件的系统级参数化组件库中，梁组件均为等截面梁，都没有宽度随长度线性变化的变截面梁系统级组件，因此无法对具有此种结构的MEMS器件进行系统级仿真设计，所以MEMS系统级组件库迫切需要此类参数化的组件模型。

发明内容

为了克服上述技术问题，本发明提出一种针对宽度随长度方向线性变化的变截面梁的系统级建模方法，解决了当前对具有此种结构的MEMS器件无法进行系统级仿真设计的问题。

本发明提出的MEMS系统级设计中的一种变截面梁建模方法，包括如下步骤：

步骤一：建立变截面梁的三维局部坐标系坐标系xyz。

已知变截面梁弹性模量E，泊松比为υ，密度为ρ，阻尼系数为B′，以变截面梁的长度方向为x轴，宽度方向为y轴，厚度方向为z轴建立局部坐标系xyz，参阅图1。图中l为梁长度，t为厚度，w1为端面1的宽度，w2为端面2的宽度且w1≤w2，x_Cy_Cz_C为全局坐标系。

步骤二：根据结构力学得到变截面梁在xoy平面内的弯曲刚度矩阵[k_yij](i＝1，2、j＝1，2)和在xoz平面内的弯曲刚度矩阵[k_zij](i＝1，2、j＝1，2)。

将端面2约束建立悬臂梁力学模型，对于宽度随长度线性变化的变截面梁，任意截面上的宽度w(x)表达式为

$w\left(x\right)=w1+\frac{(w2-w1)x}{l}---\left(1\right)$

在xoy平面内，当变截面梁横截面为矩形时，梁任意截面处的转动惯量I_z(x)为

$I_z\left(x\right)={(w1+\frac{(w2-w1)x}{l})}^3\frac{t}{12}---\left(2\right)$

参阅图2，当梁端面1受到弯矩M_xoy1时，根据纯弯曲方程

$\frac{d^2{\omega }_{\mathrm{xoy}}}{{\mathrm{dx}}^2}=\frac{M_{\mathrm{xoy}1}}{{\mathrm{EI}}_z\left(x\right)}---\left(3\right)$

式中ω_xoy为轴线上任一点在xoy平面内的挠度。结合边界条件当x＝1时，端面2处的挠角θ_xoy2＝0、挠度ω_xoy2＝0。对式(3)积分求出在端面1(x＝0)处的挠角θ_xoy1和挠度ω_xoy1表达式分别为

${\theta }_{\mathrm{xoy}1}=\frac{6l{M}_{\mathrm{xoy}1}}{\mathrm{Et}}\frac{(w1+w2)}{{w2}^2{w1}^2}---\left(4\right)$

${\omega }_{\mathrm{xoy}1}=-\frac{6l^2{M}_{\mathrm{xoy}1}}{\mathrm{Et}}\frac{1}{{w2}^{2}w1}---\left(5\right)$

参阅图3，当梁端面1处受到集中力F_xoy1时，纯弯曲方程变为

$\frac{d^2{\omega }_{\mathrm{xoy}}}{{\mathrm{dx}}^2}=\frac{x{F}_{\mathrm{xoy}1}}{{\mathrm{EI}}_z\left(x\right)}---\left(6\right)$

同理，端面1(x＝0)处的挠角和挠度表达式分别为

${\theta }_{\mathrm{xoy}1}=-\frac{6l^2{F}_{\mathrm{xoy}1}}{\mathrm{Et}}\frac{1}{{w2}^{2}w1}---\left(7\right)$

${\omega }_{\mathrm{xoy}1}=\frac{6l^3{F}_{\mathrm{xoy}1}}{\mathrm{Et}}\frac{(3w{2}^2+2\ln \left(w1\right){w2}^2-4w1w2-2\ln \left(w2\right){w2}^2+{w1}^2)}{{(w1-w2)}^3{w2}^2}---\left(8\right)$

令：$A_1=\frac{6l}{\mathrm{Et}}\frac{(w1+w2)}{{w2}^2{w1}^2}$

$B_1=-\frac{6l^2}{\mathrm{Et}}\frac{1}{{w2}^{2}w1}$

$C_1=-\frac{6l^2}{\mathrm{Et}}\frac{1}{{w2}^{2}w1}$

$D_1=\frac{6l^3}{\mathrm{Et}}\frac{(3w{2}^2+2\ln \left(w1\right){w2}^2-4w1w2-2\ln \left(w2\right){w2}^2+{w1}^2)}{{(w1-w2)}^3{w2}^2}$

可从式(4)、(5)、(7)、(8)中解出

$F_{\mathrm{xoy}1}=\frac{A_1}{(A_1{D}_1-B_1{C}_1)}{\omega }_{\mathrm{xoy}1}+\frac{B_1}{(B_1{C}_1-A_1{D}_1)}{\theta }_{\mathrm{xoy}1}---\left(9\right)$

$M_{\mathrm{xoy}1}=\frac{C_1}{(B_1{C}_1-A_1{D}_1)}{\omega }_{\mathrm{xoy}1}+\frac{D_1}{(A_1{D}_1-B_1{C}_1)}{\theta }_{\mathrm{xoy}1}---\left(10\right)$

对于xoy平面内变截面梁弯曲的刚度方程

$\left(\begin{array}{c}{F}_{\mathrm{xoy}1}\\ M_{\mathrm{xoy}1}\\ F_{\mathrm{xoy}2}\\ M_{\mathrm{xoy}2}\end{array}\right)=\left(\begin{array}{cc}{k}_{y11}& k_{y12}\\ k_{y21}& k_{y22}\end{array}\right)\left(\begin{array}{c}{\omega }_{\mathrm{xoy}1}\\ {\theta }_{\mathrm{xoy}1}\\ {\omega }_{\mathrm{xoy}2}\\ {\theta }_{\mathrm{xoy}2}\end{array}\right)---\left(11\right)$

式中[k_yij](i＝1，2、j＝1，2)为变截面梁在xoy平面内的弯曲刚度矩阵，因此其中

$k_{y11}=\left(\begin{array}{cc}\frac{A_1}{A_1{D}_1-B_1{C}_1}& \frac{B_1}{B_1{C}_1-A_1{D}_1}\\ \frac{B_1}{B_1{C}_1-A_1{D}_1}& \frac{D_1}{A_1{D}_1-B_1{C}_1}\end{array}\right)---\left(12\right)$

利用受力平衡可求出端面2所受反力F_xoy2＝-F_xoy1和弯矩M_xoy2＝F_xoy1l-M_xoy1，可知

$k_{y21}=\left(\begin{array}{cc}\frac{-A_1}{A_1{D}_1-B_1{C}_1}& \frac{-B_1}{B_1{C}_1-A_1{D}_1}\\ \frac{A_{1}l-B_1}{B_1{C}_1-A_1{D}_1}& \frac{B_{1}l-D_1}{A_1{D}_1-B_1{C}_1}\end{array}\right)---\left(13\right)$

同理，当约束端面1，在自由端端面2处分别施加弯矩M_xoy2和集中力F_xoy2，可分别解出[k_yij]中的k_y12、k_y22，其表达式为

$k_{y22}=\left(\begin{array}{cc}\frac{A_2}{A_2{D}_2-B_2{C}_2}& \frac{B_2}{B_2{C}_2-A_2{D}_2}\\ \frac{B_2}{B_2{C}_2-A_2{D}_2}& \frac{D_2}{A_2{D}_2-B_2{C}_2}\end{array}\right)$$k_{y22}=\left(\begin{array}{cc}\frac{-A_2}{A_2{D}_2-B_2{C}_2}& \frac{-B_2}{B_2{C}_2-A_2{D}_2}\\ \frac{A_{2}l-B_2}{B_2{C}_2-A_2{D}_2}& \frac{B_{2}l-D_2}{A_2{D}_2-B_2{C}_2}\end{array}\right)---\left(14\right)$

式中$A_2=\frac{6l}{\mathrm{Et}}\frac{(w1+w2)}{{w2}^2{w1}^2}$

$B_2=\frac{6l^2}{\mathrm{Et}}\frac{1}{{w1}^{2}w2}$

C₂＝B₂

$D_2=\frac{6l^3}{\mathrm{Et}}\frac{(-3{w1}^2-2\ln \left(w2\right){w1}^2+4w1w2-{w2}^2+2\ln \left(w1\right){w1}^2)}{{(w1-w2)}^3{w1}^2}$

在xoz平面内，梁任意截面处的转动惯量I_y(x)为

$I_y\left(x\right)=(w1+\frac{(w2-w1)x}{l})\frac{t^3}{12}---\left(15\right)$

与xoy面内弯曲同理，可求出在xoz平面内的弯曲刚度矩阵[k_zij](i＝1，2、j＝1，2)，其表达式为

$k_{z11}=\left(\begin{array}{cc}\frac{A_3}{A_3{D}_3-B_3{C}_3}& \frac{B_3}{B_3{C}_3-A_3{D}_3}\\ \frac{B_3}{B_3{C}_3-A_3{D}_3}& \frac{D_3}{A_3{D}_3-B_3{C}_3}\end{array}\right)$$k_{z21}=\left(\begin{array}{cc}\frac{-A_3}{A_3{D}_3-B_3{C}_3}& \frac{{-B}_3}{B_3{C}_3-A_3{D}_3}\\ \frac{A_{3}l-B_3}{B_3{C}_3-A_3{D}_3}& \frac{B_{3}l-D_3}{A_3{D}_3-B_3{C}_3}\end{array}\right)$

$k_{z22}=\left(\begin{array}{cc}\frac{A_4}{A_4{D}_4-B_4{C}_4}& \frac{B_4}{B_4{C}_4-A_4{D}_4}\\ \frac{B_4}{B_4{C}_4-A_4{D}_4}& \frac{D_4}{A_4{D}_4-B_4{C}_4}\end{array}\right)$$k_{z12}=\left(\begin{array}{cc}\frac{-A_4}{A_4{D}_4-B_4{C}_4}& \frac{-B_4}{B_4{C}_4-A_4{D}_4}\\ \frac{A_{4}l-B_4}{B_4{C}_4-A_4{D}_4}& \frac{{B_{4}l-D}_4}{A_4{D}_4-B_4{C}_4}\end{array}\right)---\left(16\right)$

式中$A_3=\frac{12l}{{\mathrm{Et}}^3}\frac{\ln (w1/w2)}{w1-w2}$

$B_3=\frac{12l^2}{E{t}^3}\frac{w1\ln \left(w1\right)-w1+w2-w1\ln \left(w2\right)}{{(w1-w2)}^2}$

C₃＝B₃

$D_3=\frac{6l^3}{E{t}^3}\frac{-3{w1}^2-2{w1}^2\ln \left(w2\right)+4w1w2-{w2}^2+2{w1}^2\ln \left(w1\right)}{{(w1-w2)}^3}$

$A_4=\frac{12l}{E{t}^3}\frac{\ln (w1/w2)}{w1-w2}$

$B_4=\frac{-12l^2}{E{t}^3}\frac{w1-w2+w2\ln \left(w2\right)-w2\ln \left(w1\right)}{{(w1-w2)}^2}$

C₄＝B₄

$D_4=\frac{6l^3}{E{t}^3}\frac{3{w2}^2+2{w2}^2\ln \left(w1\right)-4w1w2-2{w2}^2\ln \left(w2\right)+{w1}^2}{{(w1-w2)}^3}$

步骤三：计算变截面梁的轴向拉压刚度系数k_x1、k_x2和扭转刚度系数k_t1、k_t2。

参阅图4，建立变截面悬臂梁模型，当变截面悬臂梁端面2约束、端面1受到轴向力F_N和扭矩T作用时，变截面梁的轴向位移Δl和扭角的表达式分别为

$\mathrm{\Delta l}={\int }_0^l\frac{F_N\mathrm{dx}}{\mathrm{Ew}\left(x\right)t}---\left(17\right)$

将(1)代入(17)、(18)，并令

$k_{x1}=\frac{1}{{\int }_0^l\frac{\mathrm{dx}}{\mathrm{Ew}\left(x\right)t}}---\left(19\right)$

$k_{t1}=\frac{1}{{\int }_0^l\frac{7}{2G}[\frac{1}{{\mathrm{tw}}^3\left(x\right)}+\frac{1}{t^{3}w\left(x\right)}]\mathrm{dx}}---\left(20\right)$

其中，G＝E/(2+υ)为剪切模量。则k_x1、k_t1分别为端面1的轴向拉压刚度系数和扭转刚度系数。

由受力平衡可知，端面2的轴向拉压刚度系数、扭转刚度系数分别为k_x2、k_t2，即k_x2＝-k_x1、k_t2＝-k_t1。

步骤四：将上述推导出的弯曲刚度矩阵[k_yij]、[k_zij](i＝1，2、j＝1，2)、轴向拉压刚度系数k_x1、k_x2和扭转刚度系数k_t1、k_t2按照矩阵方式表出，最终得出在局部坐标系下维数为12×12的变截面梁刚度矩阵K，其具体表达式为

$K=\left(\begin{array}{cccccccccccc}{k}_{x1}& 0& 0& 0& 0& 0& k_{x2}& 0& 0& 0& 0& 0\\ 0& k_{y11}\left(\mathrm{1,1}\right)& 0& 0& 0& k_{y11}\left(1.2\right)& 0& k_{y12}\left(\mathrm{1,1}\right)& 0& 0& 0& k_{y12}\left(\mathrm{1,2}\right)\\ 0& 0& k_{z11}\left(\mathrm{1,1}\right)& 0& k_{z11}\left(\mathrm{1,2}\right)& 0& 0& 0& k_{z12}\left(\mathrm{1,1}\right)& 0& k_{z12}\left(\mathrm{1,2}\right)& 0\\ 0& 0& 0& k_{t1}& 0& 0& 0& 0& 0& k_{t2}& 0& 0\\ 0& 0& k_{z11}\left(\mathrm{2,1}\right)& 0& k_{z11}\left(\mathrm{2,2}\right)& 0& 0& 0& k_{z12}\left(\mathrm{2,1}\right)& 0& k_{z12}\left(\mathrm{2,2}\right)& 0\\ 0& k_{y11}\left(\mathrm{2,1}\right)& 0& 0& 0& k_{y11}\left(\mathrm{2,2}\right)& 0& k_{y12}\left(\mathrm{2,1}\right)& 0& 0& 0& k_{y12}\left(\mathrm{2,2}\right)\\ k_{x2}& 0& 0& 0& 0& 0& k_{x1}& 0& 0& 0& 0& 0\\ 0& k_{y21}\left(\mathrm{1,1}\right)& 0& 0& 0& k_{y21}\left(\mathrm{1,2}\right)& 0& k_{y22}\left(\mathrm{1,1}\right)& 0& 0& 0& k_{y22}\left(\mathrm{1,2}\right)\\ 0& 0& k_{z21}\left(\mathrm{1,1}\right)& 0& k_{z21}\left(\mathrm{1,2}\right)& 0& 0& 0& k_{z22}\left(\mathrm{1,1}\right)& 0& k_{z22}\left(\mathrm{1,2}\right)& 0\\ 0& 0& 0& k_{t2}& 0& 0& 0& 0& 0& k_{t1}& 0& 0\\ 0& 0& k_{z21}\left(\mathrm{2,1}\right)& 0& k_{z21}\left(\mathrm{2,2}\right)& 0& 0& 0& k_{z22}\left(\mathrm{2,1}\right)& 0& k_{z22}\left(\mathrm{2,2}\right)& 0\\ 0& k_{y21}\left(\mathrm{2,1}\right)& 0& 0& 0& k_{y21}\left(\mathrm{2,2}\right)& 0& k_{y22}\left(\mathrm{2,1}\right)& 0& 0& 0& k_{y22}\left(\mathrm{2,2}\right)\end{array}\right)$

步骤五：根据矩阵结构分析理论得到变截面梁在局部坐标系下的阻尼矩阵B、质量矩阵M。按照弯扭工程理论，并略去剪切变形，可得具有两节点、每节点6自由度梁在局部坐标系下的形函数矩阵

其中无因次参数的表达式为ξ＝x/l、η＝y/l、ζ＝z/l。对于变截面梁来说，任意截面上的宽度不再为常数，因此y的变化范围为，z的变化范围为(-t/2，t/2)。因此局部坐标系下阻尼矩阵B和质量矩阵M分别为

B＝B′∫_VN^TNdV                        (22)

M＝ρ∫_VN^TNdV                         (23)

其中V为变截而梁体积。

步骤六：根据结构矩阵分析理论将局部坐标系下变截面梁的刚度矩阵K、阻尼矩阵B和质量矩阵M，转换成全局坐标系下变截面梁的刚度矩阵K_C、阻尼矩阵B_C和质量矩阵M_C

$\left(\begin{array}{c}{K}_C={\lambda }^T\mathrm{K\lambda }\\ B_C={\lambda }^T\mathrm{B\lambda }\\ M_C={\lambda }^T\mathrm{M\lambda }\end{array}\right)---\left(24\right)$

式中λ为局部坐标系向总体坐标系投影的坐标变换矩阵，其表达式为

$\lambda =\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos \alpha & \sin \alpha \\ 0& -\sin \alpha & \cos \alpha \end{array}\right)\left(\begin{array}{ccc}\cos \beta & 0& -\sin \beta \\ 0& 1& 0\\ \sin \beta & 0& \cos \beta \end{array}\right)\left(\begin{array}{ccc}\cos \gamma & \sin \gamma & 0\\ -\sin \gamma & \cos \gamma & 0\\ 0& 0& 1\end{array}\right)=$(25)

$\left(\begin{array}{ccc}\cos \gamma \cos \beta & \sin \gamma \cos \beta & -\sin \beta \\ -\sin \gamma \cos \alpha +\cos \gamma \sin \beta \sin \alpha & \cos \gamma \cos \alpha +\sin \gamma \sin \beta \sin \alpha & \cos \beta \sin \alpha \\ \sin \gamma \sin \alpha +\cos \gamma \sin \beta \cos \alpha & -\cos \gamma \sin \alpha +\sin \gamma \sin \beta \cos \alpha & \cos \beta \cos \alpha \end{array}\right)$

其中α、β、γ分别为局部坐标系x、y、z轴与全局坐标系x_C、y_C、z_C轴之间的夹角。

步骤七：利用全局坐标系下的刚度矩阵K_C、阻尼矩阵B_C和质量矩阵M_C建立结构的二阶动力方程

$K_C{X}_C+B_C{\stackrel{·}{X}}_C+M_C{\stackrel{··}{X}}_C=F_C---\left(26\right)$

其中X_C、F_C为全局坐标系下梁节点产生的位移和所受的合外力。

步骤八：应用硬件描述语言对方程(26)编码，即得变截面梁的系统级组件模型，其示意图参阅图5。变截面梁的系统级组件端口与其自由度对应，端口信号用跨量和通量表示。组件模型中机械平动端口的跨量是沿x、y、z轴方向的线位移，通量是力F_x、F_y、f_z；机械转动端口的跨量是绕x、y、z轴转动的角位移a、b、c，通量是转矩M_x、M_y、M_z。梁的尺寸参数包括长度l、宽度w1、w2和厚度t；位置参数包括梁的初始方位角α、β、γ；材料属性参数包括弹性模量E、泊松比υ、密度ρ；环境参数包括阻尼系数B′；以上参数均可变。

本发明的有益效果是：根据结构矩阵分析理论和结构力学理论建立变截面梁的二阶动力学方程，并用硬件描述语言对动力学方程进行编码，建立了变截面梁系统级组件模型，实现了具有变截面梁结构的MEMS器件的系统级仿真设计。

附图说明

图1变截面梁结构示意图

图2自由端受弯矩时示意图

图3自由端受集中力时示意图

图4自由端受轴向力和扭矩时示意图

图5变截面梁的系统级组件示意图

具体实施例

下面结合对一个实际变截面梁结构进行系统级组件建模，对本发明进行进一步说明。已知梁长l＝200um、w1＝30um、w2＝5um、t＝5um、E＝130GP_a、ρ＝2330kg/m³、υ＝0.22、B′＝0；初始方位角α＝0，β＝0，γ＝0；本发明实例的具体步骤如下：

步骤一：以变截面梁的长度方向为x轴，宽度方向为y轴，厚度方向为z轴建立局部坐标系.

步骤二：推导变截面梁在xoy平面内弯曲刚度矩阵[k_yij](i＝1，2、j＝1，2)和在xoz平面内的弯曲刚度矩阵[k_zij](i＝1，2、j＝1，2)。

将端面2约束建立悬臂梁力学模型，在xoy平面内，通过式(1)-(13)可求出

$k_{y11}=\left(\begin{array}{cc}291.293& 0.049936\\ 0.049936& 8.08679e-6\end{array}\right)$$k_{y21}=\left(\begin{array}{cc}-291.293& -0.049936\\ 0.00832267& 1.07853e-6\end{array}\right)$

根据受力平衡得出

$k_{y12}=\left(\begin{array}{cc}-291.293& 0.00832267\\ -0.049936& 1.07853e-6\end{array}\right)$$k_{y22}=\left(\begin{array}{cc}291.293& -0.00832267\\ -0.00832267& 5.86005e-7\end{array}\right)$

在xoz平面内通过式(16)，可得

$k_{z11}=\left(\begin{array}{cc}29.8228& -0.00382857\\ -0.00382857& 5.85975e-7\end{array}\right)$$k_{z21}=\left(\begin{array}{cc}-29.8228& 0.00382858\\ -0.00213598& 1.7974e-7\end{array}\right)$

$k_{z21}=\left(\begin{array}{cc}-29.8244& -0.00213609\\ 0.0038288& 1.79755e-7\end{array}\right)$$k_{z22}=\left(\begin{array}{cc}29.8244& 0.00213609\\ 0.00213609& 2.474634e-7\end{array}\right)$

步骤三：推导变截面梁的拉压刚度系数k_x1、k_x2和扭转刚度系数k_t1、k_t2。

通过式(17)-(20)，可得k_x1＝45346.5，k_t1＝1.0686e-7，k_x2＝-45346.5，k_t2＝-1.0686e-7。

步骤四：将上述推导出的弯曲刚度矩阵、轴向拉压和扭转刚度系数按矩阵方式表出，最终得出在局部坐标系下维数为12×12的变截面梁刚度矩阵K：

$K=\left(\begin{array}{cccccccccccc}45346.5& 0& 0& 0& 0& 0& -45346.5& 0& 0& 0& 0& 0\\ 0& 291.293& 0& 0& 0& 0.049936& 0& -291.293& 0& 0& 0& 0.00832267\\ 0& 0& 29.8228& 0& -0.00382857& 0& 0& 0& -29.82244& 0& -0.00213609& 0\\ 0& 0& 0& 1.0686e-7& 0& 0& 0& 0& 0& -1.0686e-7& 0& 0\\ 0& 0& -0.00382857& 0& 5.85975e-7& 0& 0& 0& 0.0038288& 0& 1.79755e-7& 0\\ 0& 0.0049936& 0& 0& 0& 8.08679e-6& 0& -0.049936& 0& 0& 0& 1.07853e-6\\ -45346.5& 0& 0& 0& 0& 0& 45346.5& 0& 0& 0& 0& 0\\ 0& -291.293& 0& 0& 0& -0.049936& 0& 291.293& 0& 0& 0& -0.00832267\\ 0& 0& -29.8228& 0& 0.00382858& 0& 0& 0& 29.8244& 0& 0.00213609& 0\\ 0& 0& 0& -1.0686e-7& 0& 0& 0& 0& 0& 1.0686e-7& 0& 0\\ 0& 0& -0.00213598& 0& 1.7974e-7& 0& 0& 0& 0.00213609& 0& 2.474634e-7& 0\\ 0& 0.00832267& 0& 0& 0& 1.07583e-6& 0& -0.00832267& 0& 0& 0& 5.86005e-7\end{array}\right)$

步骤五：通过(21)、(22)、(23)求出局部坐标系下的阳尼矩阵B和质量矩阵M，由于本例B′＝0故B＝[0]，质量矩阵M为

$M=\left(\begin{array}{cccccccccccc}1.84458e-11& 0& 0& 0& 0& 0& 6.79583e-12& 0& 0& 0& 0& 0\\ 0& 2.097e-11& 0& 0& 0& 5.38119e-16& 0& 5.24939e-12& 0& 0& 0& -2.64683e-16\\ 0& 0& 2.09725e-11& 0& -5.38131e-16& 0& 0& 0& 5.21192e-12& 0& 2.66213e-16& 0\\ 0& 0& 0& 1.00724e-21& 0& 0& 0& 0& 0& 2.40686e-22& 0& 0\\ 0& 0& -5.38131e-16& 0& 1.83225e-20& 0& 0& 0& -2.38535e-16& 0& -1.16528e-20& 0\\ 0& 5.38119e-16& 0& 0& 0& 1.75177e-20& 0& 4.33504e-16& 0& 0& 0& -1.17226e-20\\ 6.79583e-12& 0& 0& 0& 0& 0& 8.7375e-12& 0& 0& 0& 0& 0\\ 0& 5.24939e-12& 0& 0& 0& 4.33504e-16& 0& 9.35804e-12& 0& 0& 0& -3.18597e-16\\ 0& 0& 5.21192e-12& 0& -2.38535e-16& 0& 0& 0& 9.32255e-12& 0& 3.16287e-16& 0\\ 0& 0& 0& 2.40686e-22& 0& 0& 0& 0& 0& -6.71493e-22& 0& 0\\ 0& 0& 2.66213e-16& 0& -1.16528e-20& 0& 0& 0& 3.16287e-16& 0& 1.28493e-20& 0\\ 0& -2.64683e-16& 0& 0& 0& -1.17226e-20& 0& -3.18597e-16& 0& 0& 0& 1.28493e-20\end{array}\right)$

步骤六：已知初始方位角为α＝0，β＝0，γ＝0，所以转换矩阵λ为单位矩阵。因此全局坐标下的刚度K_C、阻尼B_C和质量M_C矩阵等于局部坐标系下的刚度K、阻尼B和质量M矩阵。

步骤七：将矩阵K_C、B_C和M_C代入式(26)，建立变截面梁的二阶动力学方程。

步骤八：利用硬件描述语言对式(26)进行编码，即得变截面梁的系统级组件模型，参阅图5，图中x1、y1、z1表示端面1处的线位移端口；x2、y2、z2表示端面2处的线位移端口；a1、b1、c1表示端面1处的角位移端口；a2、b2、c2表示端面2处的角位移端口；F_x1、F_y1、F_z1表示端面1处的力端口；F_x2、F_y2、F_z2表示端面2处的力端口；M_a1、M_b1、M_c1表示端面1处的转矩端口；M_a2、M_b2、M_c2表示端面2处的转矩端口。