董金新 林謝昭 應(yīng) 濟(jì) 李 俊

1.溫州職業(yè)技術(shù)學(xué)院，溫州，325035 2.浙江大學(xué)，杭州，310027

0 引言

微流體系統(tǒng)中大都要求有微泵。其中，薄膜往復(fù)振動(dòng)式無(wú)閥微泵是以擴(kuò)散／收縮單元為閥門(mén)，通過(guò)泵膜的振動(dòng)驅(qū)動(dòng)流體，其制造簡(jiǎn)單，可以驅(qū)動(dòng)一些非均相的流體，是眾多研究的焦點(diǎn)。根據(jù)驅(qū)動(dòng)原理不同，無(wú)閥微泵可劃分為電磁、靜電、壓電、形狀記憶合金等多種類型。其中，靜電無(wú)閥微泵具有較低功耗，易于與IC工藝集成兼容等優(yōu)點(diǎn)，而得到廣泛關(guān)注。

由于微泵的尺度很小，薄膜往復(fù)振動(dòng)式無(wú)閥微泵的絕大部分物理量難以通過(guò)實(shí)驗(yàn)測(cè)量，特別是瞬態(tài)量。因此，在微泵的研制階段，利用數(shù)值方法對(duì)動(dòng)態(tài)性能進(jìn)行仿真、預(yù)測(cè)顯得很重要。它不僅可降低費(fèi)用，而且能更好地了解微泵的工作原理以及可能的潛在缺陷。但因微泵中柔性泵膜、電驅(qū)動(dòng)和流體相互耦合，加之流體在擴(kuò)散／收縮單元中的不同流動(dòng)方向表現(xiàn)出的不同壓力特性，增加了對(duì)整個(gè)泵的仿真難度和復(fù)雜性。早期的研究采用低階集總參數(shù)模型［1－2］和等效電網(wǎng)絡(luò)模型［3］等簡(jiǎn)化方法，忽略了參數(shù)的空間分布特性。目前大部分的研究沒(méi)有考慮泵的流固耦合效應(yīng)，也無(wú)法建立電驅(qū)動(dòng)信號(hào)量與流體動(dòng)力學(xué)量（比如流量）之間的關(guān)系?；蚴茄芯框?qū)動(dòng)器的動(dòng)態(tài)特性［4－7］；或是簡(jiǎn)化流體場(chǎng)，根據(jù)擴(kuò)散／收縮單元的壓力損失系數(shù)，將泵膜的運(yùn)動(dòng)與泵腔內(nèi)外壓力以及擴(kuò)散／收縮單元流量聯(lián)系起來(lái)，研究泵膜的動(dòng)態(tài)特性［8－11］；或是將驅(qū)動(dòng)器的位移輸出，作為流體場(chǎng)的移動(dòng)壁面問(wèn)題來(lái)研究泵流體的動(dòng)力學(xué)特性［12－14］。實(shí)際上，泵膜的動(dòng)態(tài)特性與流體黏滯損失的非線性和不穩(wěn)定性等是相互影響的。為了理解整個(gè)微泵的動(dòng)態(tài)特性，模型必須能有效地描述這些耦合效應(yīng)［15－16］。本文以靜電無(wú)閥微泵為對(duì)象，建立靜電－結(jié)構(gòu)－流體全耦合的3維模型，利用數(shù)值方法，仿真并分析各個(gè)變量的非線性動(dòng)態(tài)特性。

1 模型的建立

1.1 控制方程、邊界條件以及初始條件

圖1所示為靜電無(wú)閥微泵的結(jié)構(gòu)示意圖。對(duì)于泵送流體，采用任意拉格朗日－歐拉（ALE）描述比較方便。某一時(shí)刻，在參考系中＾x位置流動(dòng)的流體質(zhì)點(diǎn)速度u的時(shí)間導(dǎo)數(shù)可以表示為

圖1 靜電無(wú)閥微泵的結(jié)構(gòu)示意圖

流體域Ωf的不可壓縮流動(dòng)動(dòng)量守恒定律以及連續(xù)方程可以表示為

式中，ρf、f、σ分別為流體密度、體積力向量和Cauchy應(yīng)力張量；［0，T］為所考察的時(shí)間間隔。

牛頓流體的本構(gòu)方程為

在泵膜流固耦合界面Γf－s處，邊界條件如下：

式（5）為流固界面Γf－s處的無(wú)滑移條件；在流固界面Γf－s上，參考系的速度＾v及其位置＾x需要滿足連續(xù)條件式（6），以保證參考系的邊界準(zhǔn)確地表達(dá)結(jié)構(gòu)的當(dāng)前構(gòu)型。沿著界面的應(yīng)力平衡由式（7）表達(dá)。

對(duì)于泵膜，可以利用標(biāo)準(zhǔn)的Lagrangian描述建立其運(yùn)動(dòng)學(xué)方程，即

式中，ρs為固體的密度。

泵膜流固耦合邊界Γf－s滿足：

邊界條件式（9）、式（10）分別對(duì)應(yīng)于式（5）和式（7）。

如果沒(méi)有包含自由電荷，則描述靜電場(chǎng)的Poisson方程可以表示為

式中，φ為靜電勢(shì)；ε為介電常數(shù)張量；E為靜電場(chǎng)電場(chǎng)強(qiáng)度。

Maxwell靜電應(yīng)力張量TE由下式計(jì)算：

式中，F(xiàn)E為靜電載荷。

各個(gè)場(chǎng)變量的初始條件都設(shè)置為零。為了計(jì)算易于收斂，將一個(gè)很小的初值輸入到模型中。

1.2 模型參數(shù)及材料屬性

圖1所示靜電微泵的泵膜厚度為10μm，泵腔深度為100μm，泵膜半徑為2000μm；擴(kuò)散／收縮單元的幾何結(jié)構(gòu)相同，長(zhǎng)度為1000μm，最小端面寬度尺寸為80μm，擴(kuò)張角為10°；連接的進(jìn)出口管道長(zhǎng)度皆為2000μm，寬度為200μm。流體介質(zhì)為去離子水，物性參數(shù)如表1所示。

表1 與微泵仿真相關(guān)的材料屬性參數(shù)

假設(shè)層流流動(dòng)，考慮對(duì)稱性，在CFD－ACE中建立微泵的一半模型。利用結(jié)構(gòu)化網(wǎng)格劃分流體，單元數(shù)為52 880個(gè)，泵膜劃分為4974個(gè)殼單元。使用基于壓力修正的SIMPLE－C 算法對(duì)流動(dòng)以及能量進(jìn)行順序積分，獲得流場(chǎng)的解。流體域的對(duì)流和擴(kuò)散采用一階迎風(fēng)格式，時(shí)間積分采用了Crank－Nicolson格式。使用基于有限體積法的計(jì)算格式求解靜電場(chǎng)，這樣可以處理不同介質(zhì)介電常數(shù)問(wèn)題。利用一階殼單元來(lái)表達(dá)泵膜，并進(jìn)行大變形和接觸等幾何非線性分析?？紤]到靜電、結(jié)構(gòu)以及流動(dòng)耦合求解的需要，利用網(wǎng)格重劃分工具，使用預(yù)測(cè)－校正方法對(duì)局部變形進(jìn)行連續(xù)修正，將流動(dòng)速度和結(jié)構(gòu)速度聯(lián)系起來(lái)，實(shí)現(xiàn)流固耦合分析。流體和結(jié)構(gòu)求解器的耦合迭代頻率設(shè)置為fc＝1，使得流體場(chǎng)和結(jié)構(gòu)場(chǎng)信息能及時(shí)得到交換和更新。當(dāng)量綱一殘差小于10－4時(shí)，認(rèn)為計(jì)算是收斂的。施加的周期性驅(qū)動(dòng)電壓V＝200（1－cos（20πt））（泵腔底電壓為零），在Pentium－3.0GHz的PC機(jī)中，完成3.5個(gè)周期的求解，共350個(gè)子步數(shù)，求解時(shí)間約為59h。

2 結(jié)果及討論

2.1 泵膜的動(dòng)力學(xué)特性

圖2所示為計(jì)算得到的泵膜中心點(diǎn)垂直方向上的位移隨時(shí)間變化曲線，圖中最大變形量為28.433μm?？紤]流固耦合效應(yīng)后，在流體阻尼作用下，泵膜無(wú)法回復(fù)至初始零位置，0.653μm為新的平衡位置。穩(wěn)態(tài)情況下，泵膜位移的變化滯后于輸入電壓信號(hào)的變化，一個(gè)周期的滯后時(shí)間約為0.01s。

圖2 泵膜中心垂直方向上的位移隨時(shí)間變化曲線

2.2 流場(chǎng)的瞬態(tài)特性

實(shí)際上，泵流體動(dòng)態(tài)特性更為設(shè)計(jì)者所關(guān)注，為了更好地詮釋泵膜振動(dòng)與流體特性之間的耦合關(guān)系，用圖3所示的泵膜流固耦合界面上中心點(diǎn)處的垂直方向位移與泵腔內(nèi)流體壓力來(lái)顯示其對(duì)應(yīng)關(guān)系，從圖中可以明顯地看出，流體壓力的變化與泵膜的位移存在著密切關(guān)系。在靜電力的作用下，泵膜下拉變形的速度大于回彈的速度，所以在流體壓出泵腔的過(guò)程中其壓力變化較大，最高壓力為945.1Pa；而在流體吸入泵腔過(guò)程中，其壓力變化相對(duì)較小，最大負(fù)壓為－558.0Pa；整個(gè)泵送過(guò)程中，流體壓力呈現(xiàn)非對(duì)稱性變化。圖3還表明，泵初始工作點(diǎn)從零點(diǎn)開(kāi)始，經(jīng)過(guò)第一個(gè)周期的瞬態(tài)過(guò)程后，第二個(gè)周期的軌跡與第三個(gè)周期的軌跡已經(jīng)重合，表明泵已經(jīng)進(jìn)入穩(wěn)定工作狀態(tài)。流體壓力的波動(dòng)周期同樣滯后于驅(qū)動(dòng)電壓的變化（約0.01s）。

圖3 泵膜中心點(diǎn)處垂直方向上位移與流體壓力的關(guān)系

圖4所示為進(jìn)出口單元處的瞬態(tài)體積流量變化曲線，從圖中可以看出，在泵膜下拉（往固定電極）運(yùn)動(dòng)階段（泵送階段），出口單元的瞬態(tài)流量大于進(jìn)口單元的瞬態(tài)流量；而在泵膜回彈階段（泵吸入階段），進(jìn)口單元的流量則大于出口單元的流量。最大的瞬態(tài)流量差發(fā)生在泵腔內(nèi)壓力達(dá)到極大值時(shí)刻，而在泵腔內(nèi)壓力與進(jìn)出口壓力相差不大的時(shí)間段，進(jìn)出口單元的流量相差不大。表明要使擴(kuò)散／收縮單元發(fā)揮“整流”作用，其兩端的壓力差必須達(dá)到一定值，并且壓力差要盡可能地大。對(duì)穩(wěn)態(tài)情況下的泵凈流量進(jìn)行積分計(jì)算，得到的凈輸出流量為10.764nL／s。

圖4 進(jìn)出口瞬態(tài)流量隨時(shí)間的變化

圖5所示為微泵橫截面方向上流體壓力的分布情況，由圖可知泵腔內(nèi)的壓力幾乎分布一致，說(shuō)明在平面布置的無(wú)閥微泵中，采用Reynolds方程來(lái)描述泵腔內(nèi)流體的瞬態(tài)動(dòng)力學(xué)特性具有一定的可行性。另外，微泵所連接的進(jìn)出口管道對(duì)微泵的工作性能也有影響。

圖5 不同時(shí)刻下的微泵橫截面方向上的壓力分布情況

考慮了流固耦合特性后，模型不僅能夠仿真泵流場(chǎng)的動(dòng)態(tài)特性，還能夠仿真泵膜的應(yīng)力特性變化。圖6所示為t＝28ms時(shí)刻的泵腔內(nèi)流體壓力分布與泵膜應(yīng)力分布情況，以及這一時(shí)刻的擴(kuò)散／收縮單元內(nèi)流體的速度場(chǎng)分布。此時(shí)，泵腔內(nèi)流體壓力達(dá)到最大值945.3Pa，泵膜中心撓度為13.20μm；泵 膜 邊 緣 處 VonMises 應(yīng) 力 為28.39MPa。微泵腔內(nèi)的流速較小，擴(kuò)張管／收縮管頸部的平均流速達(dá)到最大，計(jì)算得到的平均流速為0.0867m／s，雷諾數(shù)為7.71，遠(yuǎn)小于宏觀條件下通常認(rèn)定的臨界雷諾數(shù)（2000左右），這也證實(shí)了前面假設(shè)層流流動(dòng)模式的正確性。另外，圖6b和圖6c表明，收縮／擴(kuò)散管內(nèi)未出現(xiàn)流體與固體壁面分離流動(dòng)現(xiàn)象。

圖6 t＝28ms時(shí)刻的泵膜應(yīng)力和微泵的流體場(chǎng)瞬態(tài)結(jié)果

3 結(jié)論

（1）擴(kuò)散／收縮單元的最大流量差都發(fā)生在泵腔內(nèi)壓力與進(jìn)出口壓力相差最大的時(shí)刻，表明要使得擴(kuò)散單元發(fā)揮“整流”作用，其單元兩端的壓力差必須達(dá)到一定值。

（2）在仿真的每個(gè)時(shí)刻，泵腔內(nèi)流體壓力分布幾乎一致，表明用Reynolds方程來(lái)描述泵腔內(nèi)流體動(dòng)力學(xué)特性具有可行性。

（3）仿真得到的最大雷諾數(shù)遠(yuǎn)小于通常認(rèn)定的臨界雷諾數(shù)，表明微泵流體具有層流特性。

（4）泵腔內(nèi)流體動(dòng)態(tài)特性與泵膜的運(yùn)動(dòng)存在著密切關(guān)系，如果忽略了流固耦合效應(yīng)，簡(jiǎn)單地將微泵的驅(qū)動(dòng)器與腔內(nèi)的流場(chǎng)分離處理，必然歪曲了泵流體場(chǎng)的動(dòng)態(tài)特性。泵膜變形、流體壓力和進(jìn)出口流速等響應(yīng)與輸入電壓信號(hào)之間存在著滯后現(xiàn)象。對(duì)微泵的3D全耦合仿真，能夠獲得驅(qū)動(dòng)電信號(hào)與流體流量、泵送壓力等輸出變量之間的直接關(guān)系，有利于從整體意義上實(shí)現(xiàn)微泵的優(yōu)化設(shè)計(jì)。

［1］Morris C J，F(xiàn)orster F K.Low－order Modeling of Resonance for Fixed－valve Micropumps Based on First Principles［J］.Journal Microelectromech.Syst.，2003，12：325－334.

［2］Olsson A，Stemme G，Stemme E.A Numerical Design Study of the Valveless Diffuser Pump Using a Lumped－mass Model［J］.Journal Micromech.Microeng.，1999，9：34－44.

［3］Bourouina T，Grandchamp J P.Modeling Micropumps with Electrical Equivalent Networks［J］.Journal Micromech.Microeng.，1996，6：398－404.

［4］Lin X Z，Ying J，Chen Zichen.Research on Terminal Behavior of Electrostatically Actuated Micropump Membrane Based on Modal Analysis［C］／／Wei Y L，Chong K T，Takahashi T，et al.ICMIT 2005：Mechatronics，MEMS and Smart Materials，SPIE.Bellingham，2005，6040：P60400U.

［5］Faris W F，Abdel－Rahman E M，Nayfeh A H.Mechanical Behavior of an Electrostatically Actuated Micropump［C］／／Agnes G，Chamis C，Noor A，et al.43rd AIAA／ASME／ASCE／AHS／ASC Structures，Structural Dynamics，and Materials Conference.Denver，2002：AIAA2002－1303.

［6］Morris C J，F(xiàn)orster F K .Optimization of a Circular Piezoelectric Bimorph for a Micropump Driver［J］.Journal Micromech.Microeng.，2000，10：65－459.

［7］王立文，高殿榮，楊林杰，等.壓電驅(qū)動(dòng)微泵泵膜振動(dòng)有限元分析［J］.機(jī)械工程學(xué)報(bào)，2006，42（4）：230－235.

［8］Ullman A，F(xiàn)ono I.The Piezoelectric Valve－less Pump—improved Dynamic Model［J］.Journal Microelectromech.Syst.，2002，11：655－664.

［9］Pan L S，Ng T Y，Wu X H，et al.Analysis of Valveless Micropumps with Inertial Effects［J］.Journal Micromech.Microeng.，2003，13：390－399.

［10］MacHauf A，Nemirovsky Y，Dinnar U.A Membrane Micropump Electrostatically Actuated Across the Working Fluid［J］.Journal Micromech.Microeng.，2005，15（12）：2309.

［11］白蘭，馮志慶，吳一輝.無(wú)閥微泵動(dòng)態(tài)特性的固液耦合分析［J］.機(jī)械工程學(xué)報(bào)，2008，44（7）：69－74.

［12］謝海波，傅新，楊華勇.微型無(wú)閥泵流動(dòng)特征仿真與試驗(yàn)研究［J］.機(jī)械工程學(xué)報(bào)，2002，38（7）：54－57.

［13］楊愷祥.壓電無(wú)閥式微泵浦制造與測(cè)量分析［D］.臺(tái)灣云林縣：臺(tái)灣云林科技大學(xué)，2004.

［14］Tsui Y Y，Lu S L.Evaluation of the Performance of a Valveless Micropump by CFD and Lumpedsystem Analyses［J］.Sensors and Actuators A：Physical，2008，148（1）：138.

［15］Yao Q，Xu D，Pan L S，et al.CFD Simulation of Flows in Valveless Micropumps［J］.Eng.Appl.Comput.Fluid Mech.，2007，1（3）：181－188.

［16］林謝昭，應(yīng)濟(jì)，陳子辰.靜電無(wú)閥微泵兩種仿真模型的比較研究［J］.哈爾濱工業(yè)大學(xué)學(xué)報(bào)，2010，42（12）：1943－1947.