热电耦合微执行器温度分布的节点分析法(1).pdf

第26卷 第3期 2005年3月 半 导 体 学 报 CHINESE JOURNAL OF SEMICONDUCTORS Vol. 26 No. 3 Mar. ,2005 3 国家杰出青年科学基金资助项目批准号50325519 黎仁刚 男,1978年出生,硕士研究生,主要从事热2电2机械宏模型的研究. 黄庆安 男,1963年出生,教授,博士生导师,主要从事微电子技术教学与MEMS研究. 2004204201收到,2004207211定稿○ c 2005中国电子学会 热电耦合微执行器温度分布的节点分析法 3 黎仁刚 黄庆安 李伟华 东南大学MEMS教育部重点实验室,南京 210096 摘要 MEMS系统设计的节点分析法已经被成功地用于机械或机电耦合器件的仿真,但其无法用于热执行器中 热2电耦合问题的分析仿真.本文提出一种通过傅里叶变换使用节点分析法动态分析仿真热执行器中热电耦合问 题的方法,并建立了热执行器中基本单元 梁单元的热电耦合模型.这种模型的分析计算结果与有限元软件 ANSYS吻合较好. 关键词 MEMS热执行器;梁单元;热电耦合;节点分析法 EEACC 2575 ; 8460 中图分类号 TN402 ; TN405 文献标识码 A 文章编号 02532417720050320562205 1 引言 随着MEMS系统复杂度的不断增加,设计过程 中分析仿真的难度也在不断增加,通常使用有限元 方法或者有限差分方法对具体的MEMS器件进行 仿真分析.有限元方法具有精度高的优点,但是仿真 速度慢,建模时间长,更重要的是无法对MEMS系 统进行系统级的仿真.目前,MEMS系统级仿真通 常使用节点分析法[1],虽然节点分析法具有仿真速 度快、 建模方便的优点,但是由于其本质上是一种集 总化的分析方法,所以到目前为止,还没有用于热分 布这一类问题的分布模型. 目前使用的热执行器,从运动方向可以分为垂 直于衬底和平行于衬底两种,垂直于衬底的热执行 器主要为双金属结构[2];而平行于衬底的热执行 器[3～7],虽然形状和结构变化较大,但总可以从中抽 象出简单的梁和锚点两个基本结构,绝大多数平行 于衬底运动的热执行器都可以由这两种基本结构组 合而成.理想锚点的分析较为简单,此处不加赘述. 本文主要以梁结构为研究对象,建立其热电耦合节 点分析模型. 本文从最基本的传热方程出发,建立了梁的热 电耦合模型,并使用MA TLAB和HSPICE软件求 解所建模型,同时利用有限元软件ANSYS对计算 结果加以检验,结果验证了此方法的合理性. 2 理论模型 2. 1 热电耦合梁单元基本方程 对于一长度为l ,宽度为w ,厚度为b的梁,梁中 通以大小为 it 的时变电流,如果lw ,且lb,则 可以忽略梁的截面上的温度分布,将梁上的温度视 为在梁的长度方向上的一维分布下文称为一维传 热梁 . 同时仅仅考虑梁的上表面的对流传热,忽略 下表面与衬底之间的热交换作为近似,也可以用热 阻表示 , 这样关于梁上各部分相对于环境的温度 T x ,t的瞬态传热方程为 ρcwb 5T x , t 5t - kwb 52T x , t 5x2 - hw - i2tρ0ξ wb T x ,t i2tρ0 wb 1 其中 ρ,c, k,ρ 0,ξ分别为梁材料的密度、 定容比热、 环境温度下的电阻率和电阻率的温度系数; h为对 流传热系数.这里考虑忽略电阻率随温度变化的简 第3期黎仁刚等 热电耦合微执行器温度分布的节点分析法 单情况,方程1可简化为 ρcwb 5T x ,t 5t - kwb 52T x , t 5x2 - hw T x , t i2tρ0 wb 2 为了推导方便,令αpk/ρc,ξh/ kb, mρ0/ kw2b2, 考虑梁两端的边界条件,则方程2表述为 1 αp 5T x ,t 5t - 52T x ,t 5x2 -ξT x , t i2t m T |x 0φ1 t , T |x lφ2t 3 此方程为非齐次偏微分方程,且方程的边界条件也 非齐次,为求解此方程,令T x , tv x , tw x , t ,其中w x , tφ 2t-φ1t l xφ1 t , 方程3 演变为 [8] 1 αp vt- vxxξv i2t m λ1t x λ2t v |x 00, v |x l0 4 其中 λ1t ξ[φ1t -φ2 t ] l φ ′ 1t -φ ′ 2t αpl λ2t -ξ φ1 t - φ′ 1t αp 5 方程4解的基本形式为 v x , t ∑ ∞ n1 Cntsin nπx l 6 代入到方程4 , 得出 ∑ ∞ n0 1 αC ′ nt ξ nπ l 2 Cntsin nπx l i2t m λ1t x λ2t7 同时,将方程的右端在梁的长度方向即x坐标方 向上作周期性奇延拓,之后再对空间坐标x作空 间傅里叶变换,并将λ1t和λ2t代入,得 1 αp C ′ nt ξ nπ l 2 Cnt - 2 nα π φ′ 1t 2ξ nπ φ1t - 1 n2 nα π φ′ 2t 2ξ nπ φ2t [1- -1 n ] 2i2 t m nπ 8 方程的初始条件为 Cn0 2 l∫ l 0 T x ,0 - φ20 - φ10 l x -φ10sin nπx l dx 9 6、8、9三式即为瞬态一维传热梁的传热方程. 通过以上的分析计算,将带有两个自变量的二阶偏 微分方程化为一系列的常微分方程,每个方程仅与 时间有关,而与空间坐标无关.由于方程8的右端 是由空间傅里叶变换而来的,其值一般情况下都随 方程的阶数n的增大而很快减小,方程的解Cn随着 阶数n的增加而不断减小,因而在一定的精度要求 之下可以忽略某一阶数以上的方程,而只保留有限 项.另外,一般情况下,梁的初始温度分布均匀且与 环境温度相等,这时9式中的边界条件为零.这样 的常微分方程用许多常微分求解工具都可以求解, 而现有节点法模型的求解也都是基于常微分方程求 解工具的,在这一点上就可以将热电耦合一维传热 梁模型和现有的节点法分析模型统一起来. 2.2 方程分析 方程8的右端是由空间奇延拓之后再作傅里 叶变换而来,流经梁上各个坐标点的电流值相等,在 经过奇延拓之后,梁的两个节点上的值都不连续,在 作傅里叶变换后会存在明显的吉布斯现象 [9] ,节点 的温度分布同样如此,所以空间分布的各次谐波分 量随阶数增大的衰减相对较慢.由于这个原因描述 一个一定精度的一维传热梁所需的阶数较多,且在 节点处的精度也相对较低,这是该方法的主要缺点. 梁上任意坐标的温度T x , t都是函数v x , t 与w x , t之和,而且v x , t与w x , t在坐标上任 意点均连续且可导.流过梁的节点的热流量可以表 示为 f t ∑ ∞ n0 fnt10 对于左右节点分别有 fleft ,0t kwbwx0, t kwbφ 2 t - φ1t l fleft ,n kwbvnx0,t kwbnπCnt l 11 fright ,0t - kwbwx0,t kwbφ 1t -φ2t l fright ,n - kwbvnx0,t - 1 n1 kwbnπCnt l 12 v x , t中的每一阶都对应于节点分析模型中的一个 子单元, w x , t也对应于一个子单元,通过对8、 11、12式的分析,可以将这些子单元分为两类 365 半 导 体 学 报第26卷 1集总子单元,即w x , t ,这个单元将梁视 为一个简单的热阻, R l kwb ,温度在梁中呈线性分 布,除节点外与梁外界无热交换; 2分布子单元,即v x , t ,根据v x , t级数 中n的奇偶性,亦可将此单元分为两类 Ⅰ n 为奇数时,梁的左右节点的热流大小相 等,方向相反即同时流出或同时流入 , 且8式中 含电流的项不为0,可将此类项视为梁中流过的电 流产生焦耳热和梁上各部分对外传热共同作用的结 果; Ⅱ n 为偶数时,梁的左右节点的热流大小相 等,方向相同即从一端流入另一端流出 , 8式中 含电流的项为0,梁中无焦耳热产生,可将其视为梁 两个节点的温度梯度作用的结果. 考虑联系梁两端所加电压和梁中流过电流的欧 姆定律,综合8 , 10 , 11 , 12式,即可建立热电 耦合一维传热梁的节点方程,每个节点的横越量有 节点电压和节点温度,流量有电流和热流量.由于热 流量为8 , 11 , 12式中各阶的热流量之和,而在 各阶方程中使用相同的节点温度,很显然,各阶之间 的连接方式应该为并联.同时因为热流项不是节点 温度的显式函数,而是有关热流的微分方程,故此节 点模型无法使用传统的节点法求解,必须使用改进 的节点法 [10] . 3 结果与讨论 3.1 单个一维传热梁的仿真结果 本文使用MATLAB编写了单个热电耦合一维 传热梁的仿真程序,并与ANSYS的仿真结果加以 比较.实验中梁的几何参数为l 200μm ,b 2μm ,w 2μm.梁的上表面与空气的对流系数为1104W m - 2 ℃ - 1 ,不考虑下表面与衬底之间以及结构 各部分之间的传热,不考虑辐射传热,密度为 2330kgm - 3 ,电阻率为11310 - 5Ω m - 1 ,热导率 为131Wm - 1 ℃ - 1 ,比 热 容 为700Jkg - 1 ℃ - 1 ;环境温度为0℃,梁的初始温度分布均匀且 与环境温度相等. 实验1 ,梁中无电流通过,梁左端温度固定于 0℃,右端在1ms的时间内从0℃ 开始均匀地上升到 300℃.分别对距梁左端50 ,100和150μm三点的温 度变化使用MATLAB和ANSYS进行仿真,其中 MATLAB程序对一维传热梁取10阶精度,AN2 SYS中的选取SOLID69单元,仿真结果如图1所 示. 图1 实验1中各坐标温度随时间变化曲线 Fig. 1 Temperature versus time in the 1st simulation 实验2 ,梁两端的温度固定于0℃,0时刻在梁 的两端施加3V的理想阶越电压.分别对距梁的左 端30 ,60和100μm三点的温度变化使用MATLAB 和ANSYS进行仿真,其中MATLAB程序对一维 传热梁取10阶精度,ANSYS中的选取SOLID69 单元,仿真结果如图2所示. 图2 实验2中各坐标温度随时间变化曲线 Fig. 2 Temperature versus time in the 2nd simulation 由以上的试验结果可以看出,使用本文所提模 型的计算结果与有限元软件的计算结果吻合较好. 但是在梁从初始温度开始上升时,梁上的温度分布 对空间坐标的二阶导数较大,因而产生的空间高次 谐波分量较大,而模型只取了有限项进行计算,因而 与有限元软件的计算结果的偏差相对较大;而在温 度逐渐趋向稳定之后,产生的空间高次谐波分量较 少,所以精度更高一些. 465 第3期黎仁刚等 热电耦合微执行器温度分布的节点分析法 312 基本热执行器仿真结果 图3中的基本热执行器几何尺寸为l 240μm ,lc 200μm ,lf 40μm ,wh 2μm ,wc 10μm ,g 2μm.上表面与空气的对流系数为104W m - 2 ℃ - 1 ,不考虑下表面与衬底之间以及结构 各部分之间的传热,不考虑辐射传热,密度为 2330kgm - 3 ,电阻率为11310 - 5Ω m - 1 ,热导率 为131Wm - 1 ℃ - 1 ,比 热 容 为700Jkg - 1 ℃ - 1 ;环境温度为0℃,梁的初始温度与环境温度 相同,两个锚点中的一个接地,另一个在0时刻施加 5V理想阶越电压.在HSPICE软件上使用节点法 对基本热执行器结构进行分析.实验中HSPICE软 件计算过程中对一维传热梁取10阶精度,ANSYS 仿真中使用SOLID69单元建模,下图是ANSYS瞬 态分析仿真结果和HSPICE瞬态仿真结果. 图3 a基本热执行器示意图 ;b 热执行器节点示意图 Fig. 3 a Basic thermal actuator ; b Node structure of basic thermal actuator 由图4可见,使用HSPICE软件计算的结果和 ANSYS的计算结果的相对误差控制在5 以内,其 中误差最大的节点为节点2 ,究其原因是因为交于 节点2的两根梁梁2和梁 3 的宽度相差较大,在 节点处实际的温度不能再视作一维分布,而必须考 虑其二维分布,所以使用节点法所产生的误差会相 对稍大一些. 此外,将图4与图2、 图3的数据相比较,图4 的精度明显比图2、 图3低,这是梁两端的电流和温 度分布值作空间傅里叶变换的吉布斯现象造成的, 在梁的两端傅里叶变换的误差最大,而在梁的内部 误差相对小一些;另一方面,多根梁之间误差的相互 叠加和积累也是造成误差增大的原因. 图4 基本热执行器温度阶越响应 Fig. 4 Temperature versus time the step response of basic thermal actuator 4 结论 本文提出了一种使用节点法计算热电耦合热执 行器上温度分布的方法,在MATLAB和HSPICE 软件上实现其算法,并将这两种软件的计算结果和 有限元软件ANSYS的计算结果加以比较,证实这 种方法有两个优点 1 精度较高,最大误差可控制 在5 以内 ;2 速度快,在同一台计算机上,完成基 本热执行器节点的瞬态温度分布计算, HSPICE使 用的时间为0194s ,而ANSYS软件的计算时间在 5min以上.因而,此方法可以为热执行器的设计提 供参考. 参考文献 [ 1 ] Mukherjee T ,Fedder G K,Blanton R D. Hierarchical design and test of integrated microsystems. IEEE Design beam element ; coupled thermomechanics; nodal analysis EEACC 2575 ; 8460 Article ID 02532417720050320562205 3Project supported by National Science Fund for Distinguished Young ScholarNo. 50325519 Li Ren’gang male ,was born in 1978 ,master candidate. His research interests are in coupled electro2thermo2mechanical macromodel. Huang Qing’an male ,was born in 1963 ,professor. His research interests include MEMS and microelectronic devices. Received 1 April 2004 ,revised manuscript received 11 July 2004○ c 2005 Chinese Institute of Electronics 665