专利名称：：微机电系统的温度参数化降阶建模方法
技术领域：
：本发明涉及一种微机电系统(MEMS)的温度参数化降阶建模方法，属微机电系统设计与模型降阶领域。
背景技术：
：MEMS器件结构通常采用硅为材料，通过半导体工艺加工而成。由于硅为一种热敏材料，它的物理性质受温度的影响较大，当温度变化时，硅的杨氏模量、热扩散系数和热应力会发生较大变化，而这些物理性质的变化会影响MEMS器件的谐振频率，从而导致输出信号有较大误差。为了縮短研发周期，降低产品成本，需要提高设计和仿真效率。一般可以对原始模型进行处理，采用一定的降阶方法，在维持输入输出特性基本不变的情况下使其自由度大幅縮减，建立能反应原始模型特性的低阶模型来达到高效快速分析和评价的目的。但是传统的降阶方法往往不能保留原始模型中的任何参数，这样如果原始模型的参数发生改变时就需要多次降阶，导致较大的计算量。为了在降阶过程中保护原参数系统的参数，就需要采用参数化的降阶算法。复旦大学冯丽红在申请号为200510025270.9的专利"一种参数系统的模型降阶方法"中，针对热学模型发明了一种参数系统的模型降阶方法，但是其主要针对的是一阶单参数系统，不能解决二阶系统的温度参数化问题。
发明内容本发明的目的在于提供一种微机电系统降阶建模方法，能够实现快速、精确地分析温度对微机械结构性能的影响，并且可以实现在系统级求解器中与接口电路进行联合仿真，进而可以分析出温度对整个微机电系统性能的影响。本发明解决其技术问题所采用的技术方案为一种微机电系统的温度参数化降阶建模方法，该方法包括以下步骤第一步提取质量、刚度、温度应力刚度矩阵，建立与温度相关的二阶动力学方程+z)雄)+[五0(1+rc五Ar)ii:+￡0(i+rc五Ar)rc;腐];c(o<=A+￡rQA7^(1)洲=其中，^为时间变量,x^状态空间变量，M,D，《和S分别代表质量、阻尼、刚度、温度应力刚度矩阵；A为驱动力载荷分布，A为热载荷分布矩阵；X^为输出向量，丄为输出矩阵；参数rc/z,rcs，zir分别为材料密度、常温下的杨氏模量、热扩散系数、杨氏模量的热系数、温度变化量；阻尼可用瑞利阻尼来表示，即阻尼矩阵/)是质量和刚度矩阵的线性组合，即"="^1/+-《，其中a，-为瑞利阻尼系数，其与结构模态阻尼比相关；第二步构造正交映射矩阵K首先构造映射矩阵P^和F^co/印朋^}=(〖—'M,〖—'S)(2)c咖—K2}=&(〖-、,〖-'5)(3)其中，矩阵S为输入向量的组合，即3=^+6,;,、1，2指的是々阶Krylov子空间，"阶Krylov子空间定义为《(爿，B)二co/^a"(S,爿25,…，爿"—接下来应用Amoldi算法对映射矩阵K、^进行正交化，然后计算"、^的并集f=^UK2(4)最后对P中所有列向量经正交化处理后得到映射矩阵K。第三步将方程(1)投影到正交映射矩阵F所在的子空间上，从而得到温度参数化的低阶模型"《+￡rc,,r(5)其中，Afr=rrA/T，JDr=FrDr=",,+/^/r〖K=aMr+y^:r，《=^'1「，S,)Z，&=，6,，丄,.="丄。本发明的有益效果如下51.高的计算效率和降阶精度通过降阶算法可以使原始模型的自由度规模大幅縮减，从而使基于降阶模型的系统级建模与仿真速度加快。同时，采用矩匹配原理，可以使降阶模型与原始模型的传递函数较好近似，从而有较高的精度。2.降阶模型保留了温度参数，实现了温度参数化的降阶建模在降阶过程中，可以生成与参数无关的投影矩阵，最终使原系统中温度参数在低阶模型中得以保留。下面结合附图和实施例对本发明进一步说明。图l是实施例中的Z轴微机械陀螺；图2是常温下微机械陀螺x方向和y方向的谐响应曲线图3是常温下微机械陀螺降阶模型x方向和y方向的相对误差图4是10(TC时微机械陀螺x方向和y方向的谐响应曲线图5是10(TC时微机械陀螺降阶模型x方向和y方向的相对误差图6是-100°C时微机械陀螺x方向和y方向的谐响应曲线图7是-10(TC时微机械陀螺降阶模型x方向和y方向的相对误差图中标号，l为锚点，2为平板质量块，3为支撑梁，4为驱动电极，5为检测电极。具体实施例方式下面结合一个Z轴微机械陀螺对本发明进行进一步说明。如图l，该Z轴微机械陀螺结构由一个平板质量块、16根支撑梁、两个驱动电极、两个检测电极和四个锚点组成。该结构通过施加于梳齿上的直流偏交流电压产生的静电力来驱动陀螺在x方向做往复运动。当外界有绕Z轴的角速度输入时，就会科氏力产生，从而使质量块在y方向振动。于是，通过检测敏感方向梳齿电容的变化来输出检测信号。该实施例中，Z轴微机械陀螺的温度参数化降阶建模方法包括以下步骤第一步利用有限元软件Ansys将该Z轴陀螺的三维几何结构模型进行网格划分，得到的网格模型由45092个四面体结构单元(Solidl87)禾Q93139个节点组成。在锚点处施加狄利克雷边界条件后，共有1656个自由度被约束。为了观察陀螺结构的动态响应，选取驱动方向和检测方向两个悬浮梁的中心节点作为输出点。采用有限元法提取微机械陀螺的质量、刚度、温度应力刚度矩阵，建立与温度相关的微机械陀螺二阶动力学方程+w)+[￡0(1+rc五Ar)ii:+￡0(1+rc五Ar)rc力A73>(o<=&+￡7^A7S,+(6)另外，微机械陀螺通过科氏效应来检测外界角速度信号，因此还需要提取科氏矩阵。方程(6)中w为角速度，矩阵C为阻尼矩阵D和科氏矩阵G组成，即C:D+z^G。另外，陀螺结构旋转将会产生离心力，A为离心载荷分布矩阵。第二步构造正交映射矩阵F。首先在计算投影矩阵R和^。根据式(2)和(3)所构成的Krylov子空间，这里S二&+5,+^，利用块Arnoldi算法分别得到映射矩阵R和F2的列向量数分别为35和40。对映射矩阵)^和^按照式(4)取并集得到的p，然后对组合形成的矩阵列向量进行正交得到投影矩阵F，此时F中有70列向量。第三步将方程(6)投影到r所在的子空间得到与温度相关的微机械陀螺的低阶模型，其阶数为70。++[￡0(i+rc五Ar)^+￡0(1+7r^Ar)rc;A7^];0)肌，，2《(7)这样，本实例首先通过提取系统矩阵建立与温度相关的陀螺二阶系统方程，然后根据系统矩阵构造映射矩阵F，最后将陀螺二阶系统方程投影到矩阵F所在的子空间，得到了与温度相关的陀螺低阶模型，从而实现了温度参数化的降阶建模。图2-图8为仿真结果。在常温常压下，当外界有100rad/s的角速度输入时，陀螺在驱动方向(x方向)和检测方向(y方向)的谐响应如图2所示。从图2可以看出降阶模型与原始模型的仿真结果即使在高频处也能较好的重合。图3为陀螺降阶模型在x方向和y方向的相对误差，可以看出两者相对误差在0.2%以内。为了研究温度变化对陀螺结构的影响，在10(TC时陀螺的谐响应曲线如图4所示。从图4可以看出当温度为IO(TC时，降阶模型与原始模型的仿真结果较好地重合。从图5可以看出两者相对误差在0.17%以内。在-10(TC时陀螺的谐响应曲线如图6所示。从图6可以看出当，降阶模型与原始模型的仿真结果也较好地重合，从图7可以看出两者相对误差在0.42%以内。从图2、图4和图6可以看出，当温度从-10(TC变化到IO(TC的过程中，由于Z轴陀螺结构上的对称性，其驱动模态和敏感模态仍较好的匹配。在100-100kHz的范围内对微机械陀螺进行谐响应分析，当计算190个频率点时，有限元模型的计算时间为33377s，降阶模型的计算时间为1220s。由此可见，降阶模型的计算效率与有限元模型计算效率相比提高了27倍，但最大相对误差保持在0.5%以内，具有较高的精度。权利要求1.一种微机电系统的温度参数化降阶建模方法，其特征在于，该方法包括以下步骤第一步提取质量、刚度、温度应力刚度矩阵，建立与温度相关的二阶动力学方程<mathsid="math0001"num="0001"><math><![CDATA[<mfencedopen='{'close=''separators=''><mtable><mtr><mtd><mi>ρM</mi><mover><mi>x</mi><mrow><mo>·</mo><mo>·</mo></mrow></mover><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>+</mo><mi>D</mi><mover><mi>x</mi><mo>·</mo></mover><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>+</mo><mo>[</mo><msub><mi>E</mi><mn>0</mn></msub><mrow><mo>(</mo><mn>1</mn><mo>+</mo><msub><mi>TC</mi><mi>E</mi></msub><mi>ΔT</mi><mo>)</mo></mrow><mi>K</mi><mo>+</mo><msub><mi>E</mi><mn>0</mn></msub><mrow><mo>(</mo><mn>1</mn><mo>+</mo><msub><mi>TC</mi><mi>E</mi></msub><mi>ΔT</mi><mo>)</mo></mrow><msub><mi>TC</mi><mi>h</mi></msub><mi>ΔTS</mi><mo>]</mo><mi>x</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mo>=</mo><msub><mi>B</mi><mi>d</mi></msub><mo>+</mo><msub><mi>ETC</mi><mi>h</mi></msub><mi>Δ</mi><msub><mi>TB</mi><mi>t</mi></msub></mtd></mtr><mtr><mtd><mtext>y</mtext><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>=</mo><msup><mi>L</mi><mi>T</mi></msup><mi>x</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mtd></mtr></mtable><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow></mfenced>]]></math></maths>其中，t为时间变量，x(t)状态空间变量，M，D，K和S分别代表质量、阻尼、刚度、温度应力刚度矩阵；Bd为驱动力载荷分布，Bt为热载荷分布矩阵；y(t)为输出向量，L为输出矩阵；参数ρ，E0，TCh，TCE，ΔT分别为材料密度、常温下的杨氏模量、热扩散系数、杨氏模量的热系数、温度变化量；第二步构造正交映射矩阵V首先构造映射矩阵V1和V2<mathsid="math0002"num="0002"><math><![CDATA[<mrow><mi>colspan</mi><mo>{</mo><msub><mi>V</mi><mn>1</mn></msub><mo>}</mo><mo>=</mo><msub><mi>K</mi><msub><mi>r</mi><mn>1</mn></msub></msub><mo>(</mo><msup><mi>K</mi><mrow><mo>-</mo><mn>1</mn></mrow></msup><mi>M</mi><mo>,</mo><msup><mi>K</mi><mrow><mo>-</mo><mn>1</mn></mrow></msup><mi>B</mi><mo>)</mo><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>2</mn><mo>)</mo></mrow></mrow>]]></math></maths><mathsid="math0003"num="0003"><math><![CDATA[<mrow><mi>colspan</mi><mo>{</mo><msub><mi>V</mi><mn>2</mn></msub><mo>}</mo><mo>=</mo><msub><mi>K</mi><msub><mi>r</mi><mn>2</mn></msub></msub><mrow><mo>(</mo><msup><mi>K</mi><mrow><mo>-</mo><mn>1</mn></mrow></msup><mi>S</mi><mo>,</mo><msup><mi>K</mi><mrow><mo>-</mo><mn>1</mn></mrow></msup><mi>B</mi><mo>)</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>3</mn><mo>)</mo></mrow></mrow>]]></math></maths>其中，矩阵B为输入向量的组合，即B＝Bd+Bt；id="icf0004"file="A2009100220270002C4.tif"wi="7"he="5"top="158"left="142"img-content="drawing"img-format="tif"orientation="portrait"inline="yes"/>i＝1，2指的是ri阶Krylov子空间，n阶Krylov子空间Kn(A，B)定义为Kn(A，B)＝colspan{B，AB，A2B，...，An-1B}接下来应用Arnoldi算法对映射矩阵V1、V2进行正交化，然后计算V1、V2的并集<mathsid="math0004"num="0004"><math><![CDATA[<mrow><mover><mi>V</mi><mo>^</mo></mover><mo>=</mo><msub><mi>V</mi><mn>1</mn></msub><mo>∪</mo><msub><mi>V</mi><mn>2</mn></msub><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>4</mn><mo>)</mo></mrow></mrow>]]></math></maths>最后对id="icf0006"file="A2009100220270002C6.tif"wi="2"he="4"top="218"left="43"img-content="drawing"img-format="tif"orientation="portrait"inline="yes"/>中所有列向量经正交化处理后得到映射矩阵V。第三步将方程(1)投影到正交映射矩阵V所在的子空间上，从而得到温度参数化的低阶模型<mathsid="math0005"num="0005"><math><![CDATA[<mrow><mfencedopen='{'close=''><mtable><mtr><mtd><mi>ρ</mi><msub><mi>M</mi><mi>r</mi></msub><msub><mover><mi>x</mi><mrow><mo>·</mo><mo>·</mo></mrow></mover><mi>r</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>+</mo><msub><mi>D</mi><mi>r</mi></msub><msub><mover><mi>x</mi><mo>·</mo></mover><mi>r</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>+</mo><mo>[</mo><msub><mi>E</mi><mn>0</mn></msub><mrow><mo>(</mo><mn>1</mn><mo>+</mo><msub><mi>TC</mi><mi>E</mi></msub><mi>ΔT</mi><mo>)</mo></mrow><msub><mi>K</mi><mi>r</mi></msub><mo>+</mo><msub><mi>E</mi><mn>0</mn></msub><mrow><mo>(</mo><mn>1</mn><mo>+</mo><msub><mi>TC</mi><mi>E</mi></msub><mi>ΔT</mi><mo>)</mo></mrow><msub><mi>TC</mi><mi>h</mi></msub><mi>Δ</mi><msub><mi>TS</mi><mi>r</mi></msub><mo>]</mo><msub><mi>x</mi><mi>r</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mo>=</mo><msub><mi>B</mi><msub><mi>d</mi><mi>r</mi></msub></msub><mo>+</mo><msub><mi>ETC</mi><mi>h</mi></msub><mi>Δ</mi><msub><mi>TB</mi><msub><mi>t</mi><mi>r</mi></msub></msub></mtd></mtr><mtr><mtd><mi>y</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>=</mo><msubsup><mi>L</mi><mi>r</mi><mi>T</mi></msubsup><msub><mi>x</mi><mi>r</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>5</mn><mo>)</mo></mrow></mrow>]]></math></maths>其中，Mr＝VTMV，Dr＝VTDV＝αVTMV+βVTKV＝αMr+βKr，Kr＝VTKV，Sr＝VTSV，Bdr＝VTBd，Btr＝VTBt，Lr＝VTL。2.—种如权利要求1所述的微机电系统的温度参数化降阶建模方法，其特征在于所述第一步中的阻尼使用瑞利阻尼来表示，即阻尼矩阵Z)是质量和刚度矩阵的线性组合Z)二a^f+戶ii:，其中Ot,yS为瑞利阻尼系数。全文摘要本发明公开了一种微机电系统的温度参数化降阶建模方法，属微机电系统设计与模型降阶领域。该方法包括以下步骤提取质量、刚度、温度应力刚度矩阵，建立与温度相关的二阶动力学方程；构造正交映射矩阵V；将原二阶动力学方程投影到正交映射矩阵V所在的子空间上，从而得到温度参数化的低阶模型。该降阶方法通过降阶算法可以使原始模型的自由度规模大幅缩减，从而使基于降阶模型的系统级建模与仿真速度加快。同时，采用矩匹配原理，可以使降阶模型与原始模型的传递函数较好近似，从而有较高的精度。在降阶过程中，可以生成与参数无关的投影矩阵，最终使原系统中温度参数在低阶模型中得以保留。文档编号B81C99/00GK101567018SQ20091002202公开日2009年10月28日申请日期2009年4月15日优先权日2009年4月15日发明者常洪龙,张亚飞,徐景辉,苑伟政申请人:西北工业大学