硅流体芯片的仿真建模及几何参数影响分析

周铖杰1，吴央芳，夏春林，王玉翰，陆倩倩

(1.常州大学机械工程学院，江苏常州 213164； 2.浙江大学城市学院机械电子工程系，浙江杭州 310015)

引言

微型阀作为微流体控制系统中的关键部件之一，在生物医学、精密制造和制药行业有着重要的应用。它的种类繁多[1-4]，多适用于低压小流量流体的控制，难以适应工业生产上高压大流量的应用环境，因此在工业上的应用与生物医疗领域相比还较少[1,5]。硅流体芯片[6]是一种基于电热致动技术的微型阀，最大运行压力达3.45 MPa，在0.69 MPa的工作压力下最大流量可达5 L/min，是工业供热通风与空气调节领域的研究热点[7]。李勇俊等[8]对芯片的致动器进行了优化设计，提高了芯片的综合性能，以实现更大流量的输出；SKIMA等[9]提出了一种微机电系统的混合预测方法，并对芯片进行了疲劳测试；课题组已对芯片进行了实验及应用研究，实验结果表明基于硅流体芯片设计的两种不同结构(单/双芯片结构) 的气动比例压力阀都有良好的性能。硅流体芯片控制精度高，响应快，但是在气动控制领域的应用研究还相对较少。通过对硅流体芯片特性的仿真研究，可以降低芯片应用研究开发的成本，提高应用研究的工作效率。

电热微致动器是硅流体芯片的重要组成部分，它的响应时间和输出位移情况决定了硅流体芯片的性能。HUSSEIN等[10]提出了电热微致动器的电热和热机械的混合偏微分模型，为致动器的建模、设计和优化提拱了重要依据；SHAN等[11]通过建立带有致动器温度相关参数的有限差分模型研究了热致动器的温度分布、几何参数对其输出力和位移的影响以及脉冲输入对致动器工作效率的影响；LOTT等[12]考虑了自由流动、环境热辐射的影响，对热致动器进行有限元分析，比较了热致动器分别在真空和空气中的工作效率。学者们对电热微致动器的静态模型进行了大量的研究分析，但是对于动态响应的分析研究还相对较少。

本研究介绍了硅流体芯片的工作原理，利用多物理场仿真软件COMSOL Multiphysics对芯片进行了瞬态仿真，并通过AMESim软件对芯片进行了气路仿真，将仿真结果与实验结果进行对比，验证了仿真的有效性。最后探讨了芯片中致动器的几何结构参数对致动器的影响。

1 硅流体芯片的工作原理

硅流体芯片的结构如图1所示，几何规格：长10.8 mm，宽4.83 mm，高2.225 mm，更多的细节结构可参考专利文献[13]。芯片采用的是3层结构设计，中间层如图2所示，带有一个V型的电热微致动器和杠杆机构。ps为输入端口，pc为控制输出端口，po为反馈端口，在不通电的情况下，ps为常闭端口，po为常开端口。

图1 硅流体芯片结构图

微型阀的主流结构有两种：隔膜式结构[14]和滑阀式结构。硅流体芯片可以看成是一个滑阀式结构的微型阀，中间层的移动部分看成是阀芯。当芯片通入控制电压时，电流经过V型电热微致动器会使筋的温度升高，导致热膨胀，产生沿A方向的位移，如图2所示，B点则作为杠杆机构的支点将位移放大，以改变ps，po口的大小，达到比例控制输出压力或流量的目的。

图2 中间层结构图

芯片的3个工作状态如图3所示，状态1为初始状态(没有电压输入)，ps口关闭，po，pc口相通，都处于全开状态；状态2为中间的过渡状态；状态3为输入电压达到额定最大值时的状态，ps口与pc口连通处于全开状态，po口关闭。

图3 芯片的3个工作状态图

2 仿真分析

硅流体芯片的仿真分析路线如图4所示:首先对中间层建模，利用COMSOL Multiphysics软件对其进行多物理场仿真，得出其输入电压U与杠杆机构输出位移x之间的关系曲线，并利用MATLAB软件中的系统辨识工具箱对电压位移曲线进行系统辨识，得出电压与位移之间的传递函数；接着将此传递函数作为一个数学模型放入到AMESim软件中，结合芯片杠杆机构移动部分对应的端口过流面积关系进行气路仿真，并与实验结果进行比较来验证仿真的有效性；最后利用该仿真模型来分析芯片中致动器的结构参数的影响。

图4 仿真分析路线图

2.1 多物理场仿真

电热微致动器的仿真分析主要涉及传热、电场和结构力学3个物理场的耦合。传热的机制有3种：传导、对流以及辐射。在V型电热微致动器中，传热的主要形式为传导和对流，其传热方程形式如下：

式中，ρ为材料密度；Cp为恒压热容；T为温度；u为外部流体的速度场；k为导热系数；Q为热源项是电场中电流通过致动器产生的焦耳热，其计算公式如下：

Q=J·E

(2)

式中，J为电流密度矢量；E为电场强度矢量。

V型电热微致动器的结构如图5所示，焦耳热导致温度变化使致动器发生热膨胀产生致动效应，其作用于筋末端的合力F为[9]：

F=2EAh·ΔT·sinθ=αU

(3)

式中，E为材料的杨氏模量；A为筋的横截面积；h为材料的热膨胀系数； ΔT为温度差；θ为筋的倾角；U为输入电压；α为比例系数，与筋的材料属性和几何参数有关。

图5 V型电热微致动器结构图

将该V型电热微致动器看成是一个质量弹簧阻尼系统[9]，有：

ma=Fs+Ff+Fr+F

(4)

式中，Ff为筋运动时所受的摩擦力，由于筋与上下层表面之间没有接触，故摩擦力项为0；a为加速度；m为质量，因致动器的的尺寸较小，ma项可忽略不计；Fs为等效弹簧力。

Fs=-ks·x0

(5)

式中，ks为等效弹簧刚度；x0为致动器输出位移；Fr为等效粘滞力。

(6)

式中，β为等效阻尼系数；对式(4)进行化简，得：

(7)

对其进行拉式变换，得到输入电压与致动器输出位移之间的传递函数如下：

(8)

其中，静态增益ζ=α/ks；时间常数τ=β/ks；杠杆机构的输出位移与致动器位移成正比，其输出位移与输入电压之间的传递函数也可由式(8)表示。

使用COMSOL Multiphysics 5.2a软件对硅流体芯片中的V型电热微致动器进行仿真分析。仿真使用材料库中的单晶硅材料，材料属性及致动器参数如表1所示。软件中添加的物理场及设置的边界条件如图6所示。实验研究表明[9]，当输入电压为8 V时，杠杆机构产生的位移适中，故将其作为参考电压来进行仿真研究。

表1 材料参数及致动器几何参数

图6 物理场及边界条件图

输入电压为8 V的阶跃信号时，杠杆机构的输出位移x随时间变化如图7实线所示，利用MATLAB中的系统辨识工具箱进行系统辨识，得到传递函数如下：

(9)

辨识结果如图7虚线所示，辨识度为97.65%。在t为0.4 s时，V型电热微致动器的温度分布如图8所示，最高温度为559 K。

图7 仿真结果与MATLAB系统辨识对比图

图8 温度分布图

2.2 AMESim气路仿真

利用输入电压与杠杆机构输出位移的传递函数关系，通过AMESim软件对芯片进行气路仿真。输入电压为0～12 V时，对应的过流面积如图9所示，阴影部分表示端口的面积变化情况。

图9 输入电压对应过流面积图

仿真模型中假设芯片封装良好，没有泄漏，AMESim气路仿真模型如图10所示。由于芯片的尺寸太小，难以对其特性进行直接的实验测试，故将其贴在底座上如图11所示进行相应的特性研究。

气源压力为0.7 MPa，分别输入三角波信号和阶跃信号，双芯片结构的输出压力p如图12所示，近似线性可控制范围在30%～60%左右，闭环阶跃响应时间小于0.04 s，存在一定的泄漏，可能与封装工艺有关。此外，芯片的闭环滞回特性实验及仿真结果表明：在输入频率大于4 Hz后，随着输入频率的增加，不规则的循环回路开始逐渐扩大。输入为5 V，5 Hz的指数衰减信号，其实验与仿真的闭环滞回特性曲线如图13所示。

图10 AMESim仿真模型图

图11 单、双芯片结构图

图12 仿真与实验结果对比图

图13 闭环滞回特性曲线

对比以上实验与仿真结果可知仿真是有效的。接着利用仿真模型来讨论静态增益ζ和时间常数τ值的影响。由图12b可知，ζ值的改变会影响双芯片结构的上升时间，输入相同的控制电压，ζ越大阶跃响应上升时间越短。

时间常数τ对开环滞回特性的影响如图14所示，随着τ值减小，双芯片结构的开环滞回特性有所改善，当减小至0.001时，τ值再减小对开环滞回特性的影响不再明显。此外，当τ值为0.01时，对应的闭环滞回特性仿真结果表明：在输入频率大于8 Hz后，随着输入频率的增加，不规则的循环回路开始逐渐扩大，τ等于0.01的双芯片结构滞回特性优于τ等于0.032452。

图14 时间常数对开环滞回特性的影响

2.3 芯片致动器的仿真分析

下面利用仿真模型来分析筋的几何参数对杠杆机构位移x、致动器输出位移x0、静态增益ζ值和时间常数τ值的影响。

(1) 致动器截面及筋对数的影响。输入8 V的阶跃信号，致动器倾角不变，分别改变筋的对数及横截面积，芯片中杠杆机构的位移x变化情况如图15所示，位移从稳态值的10%上升到90%的时间均小于0.07 s。当电压变至0时，致动器没有完全回到初始位置，四对筋结构的返回情况比单对筋的要好，这可能与致动器的机械结构(如筋的对数)有关。

图15 筋截面变化及对数对输出位移影响图

由表2可知，杠杆机构的位移放大倍数约为11.3倍。4对筋截面变化的致动器输出位移为11.7 μm，与单对筋的输出位移(7.7 μm)相比增加了51.9%，与截面恒定的致动器(10.7 μm)相比输出位移增加了9%，筋的对数对致动器的输出位移影响较大。

表2 筋对数及截面变化对应的输出位移

致动器筋上温度T随跨长L的分布情况如图16所示，筋的中心位置温度最高，截面变化的筋上最高温度为559 K，截面恒定的为569.8 K，截面变化的筋的平均温升(491 K)要小于截面恒定的筋(493 K)，在相同的输入电压下，截面变化的筋与截面恒定的筋相比，能以更低的温度产生更大的位移。

(2) 致动器筋的倾角θ、宽度D、厚度H、跨长L的影响。以截面恒定采用4对筋结构的致动器为仿真对象，分别做如下改变：改变倾角1°～10°(编号1～10)；改变宽度为原来0.12 mm的3/4，7/8，1，9/8，5/4倍(编号11～15)；改变厚度为原来0.75 mm的1/4，1/2，3/2，2倍(编号16～19)；改变跨长为原来5 mm的1/2，2/3，1倍(编号20～22)；输入8 V的电压信号，改变倾角后对应致动器的最大输出位移x0max如图17所示，倾角为4 °时致动器输出的位移最大为10.86 μm，输出位移最小为5.83 μm对应的筋倾角为1 °。

图16 温度分布图

图17 改变倾角对输出位移影响图

改变筋宽度、厚度和跨长对应的致动器最大输出位移如图18所示 ，筋跨长为原来1倍的致动器的输出位移(10.7 μm)约为原来1/2倍的致动器(3.055 μm)的3.5倍。可见筋的倾角和跨长是影响致动器输出位移的两个关键因素，筋宽度和厚度的影响相对较小。此外，跨长也是影响致动器响应时间的重要因素，筋跨长分别为原来1/2倍、2/3倍、1倍致动器的输出位移从稳态值10%上升到90%的时间如图19所示，分别为0.026，0.032，0.056 s，跨长增加一倍相应的响应时间增加了约1.15倍。

输入为8 V，改变致动器几何参数对应的静态增益ζ值及时间常数τ值如表3所示，致动器筋的对数、倾角、宽度以及跨长对静态增益ζ值的影响明显，时间常数τ值主要与筋的跨长有关，随着跨长的减小，时间常数τ值减小。

图18 筋宽度、厚度、跨长对输出位移的影响图

图19 跨长对响应时间的影响图

表3 几何参数对应的静态增益ζ值及时间常数τ值

3 结论

(1) 建立的气路仿真模型和多物理场仿真模型是有效的；双芯片结构的闭环阶跃响应上升时间与静态增益ζ值有关，ζ值越大上升时间越短，降低时间常数τ值能改善结构的滞回特性；

(2) 电热微致动器作为芯片的重要组成部分，其多物理场仿真结果表明：筋的对数、倾角对致动器的输出位移影响较大；变化截面的筋与截面恒定的筋相比能以更低的温升产生更大的位移；与筋宽度和厚度的影响相比，跨长改变对输出位移带来的影响更为明显；

(3) 致动器几何参数的改变对静态增益ζ值的影响明显，时间常数τ值的大小主要与跨长有关，具体的影响函数还有待进一步讨论。

本研究研究成果，对硅流体芯片的性能优化及在气动领域的拓展应用具有一定意义。