- Email: [email protected]
Microfluidic system modeling using VHDL-AMS and circuit simulation
Microelectronics Journal ELSEVIER
Microelectronics Journal 29 (1998) 791-797
Microfluidic system modeling using VHDL-AMS and circuit simulation P. Voigt, G. Schrag, G. Wachutka Institute.for Physics of Electrotechnology, Technical University of Munich, Arcisstrafle 21, 80290 Munich, Germany
Abstract
It is demonstrated that the operation of electrofluidic microsystems such as micropumps can be efficiently modeled by standard electrical circuit simulators, using the analog hardware description language (VHDL-AMS) to describe the models. As a representative example, compact models were developed for all constituent parts of an electromechanical micropump fabricated by silicon-based microstmcture technology, including ele.ctrostatically actuated membrane drives, membrane valves and tubes together with the power supply and control circuitry. With the resulting macromodel of the micropump, the static characteristics as well as the dynamic behavior is correctly reproduced. Thus, with this approach it is possible to include fluidmechanical system components straightforwardly in conventional microelectronic design environments, allowing for easy full system analysis and optimization of even complex microsystems. © 1998 Elsevier Science Ltd. All rights reserved.
1. Introduction
The rapid progress in microstructure technology allows the fabrication of increasingly complex microsystems, thereby raising a high demand for powerful simulation tools for the design and optimization of such systems [1]. However, the development of special purpose software packages tends to be quite expensive and time consuming. On the other hand, the electronic industry has driven the development of efficiently performing design systems for microdevices. Since many aspects of microsystem design are closely related to the silicon-based integrated circuit standard technologies, notably with a view to the integration of electronics and transducers in 'smart systems', it seems advantageous to simuhtte the nonelectrical parts of a microsystem using conventional design packages available from the semiconductor industry. The appropriate way of modeling at system level is the method of tailored compact modeling as described in [2]. In this work, the practicability of this idea is demonstrated with reference to electrofluidic microsystems. It is typical of such microsystems that they are equipped with one or several electromechanical driving units such as electrostatically, electrothermally or piezoelectrically actuated micropumps, passive fluidic components such as static or dynamic valves and tubes, sensor elements for control and analysis, and a co-integrated or hybrid coupled electronic circuitry for power supply, signal conditioning and system control including (self-)test and (self-)calibration as well as error detection and compensation. By their nature, the constituent components link different energy and signal domains such 0026-2692/98/$19.00 © 1998 Elsevier Science Ltd. All rights reserved PII S 0 0 2 6 - 2 6 9 2 ( 9 7 ) 0 0 0 9 3 - 1
as mechanical, fluidic, thermal, electrical and other physical or chemical quantities. Typical applications include, among others, chemical microanalyzers, mass flowmeters, fuel drop generators for automobile heaters, ink jet print heads or microdosed drug infusion. The natural description of the involved transducer effects is provided by the concepts of continuum theory, which couples the relevant mechanical, electrodynamical and thermal field quantities in terms of a (usually quite complex) system of partial differential equations [2]. However, it is hardly feasible to solve the model equations in their continuous form for realistic microsystems in their full size and complexity. Therefore, to keep the computational effort for such problems within acceptable limits, the proper and adequate level of simulation is analog circuit simulation based on a network theory which, on the one hand, preserves the basic conservation laws expressed in the continuum models and, on the other hand, reduces the set of variables to a number that is still manageable with one of today's standard circuit simulation tools as they are routinely used in the microelectronic industry. In this approach, the microsystem is partitioned into functional blocks, which interact with each other as constituent parts of a flux-conserving ('Kirchhoffian') network. For each of them, a compact model with only few degrees of freedom is formulated, thus defining a macromodel for all parts of the microsystem. Then VHDL-AMS (analog hardware description language) is used to code the models. This has the advantage of providing a simulator-independent model code, since VHDLAMS is a quasi-standard, accepted as input language by many electronic simulation systems.
792
P. Voigt et al./Microelectronics Journal 29 (1998) 791-797
In the following, the method of tailored compact modeling (Section 2) is briefly introduced. Then the general principles basic to this approach are exemplified with reference to an electrofluidic micropump, the detailed description of which is given in [3],[4]. The compact models of the micropump components are then derived and discussed (Sections 3 and 4), and finally the transient operation of the micropump by means of time-dependent system simulation is studied (Section 5).
2. Physically-based compact modeling System models of complex microstructures can be constructed by partitioning the full system into simpler subsystems. The interfaces between interacting subsystems must be chosen in such a way that the exchange of energy, mass, charge and other physical quantities can be described by pairs of conjugate state variables. A systematic method of selecting such pairs of variables is provided by the basic principles of thermodynamics [2]. A pair of conjugate variables consists of the flux of a physical quantity and an associate driving force. Examples are given in Table 1. In the context of VHDL-AMS, the driving force is called 'across quantity' and the flux 'through quantity'. A subsystem is called 'block' and is defined by the number and nature of its terminals. The terminals allow for the exchange of flux quantities across subsystem boundaries, concentrated ('lumped') at single points instead of the distributed flow of the respective field quantities. The intemal behavior of each of the blocks is described by a physically-based compact model extracted from continuous field theory and/or from experimental data. In practice, the compact models are formulated using an appropriate analog hardware description language such as, for instance, VHDL-AMS. The full microsystem is synthesized by linking the individual HDL-specifled compact models of the constituent blocks through flux-conserving interface conditions for the respective pairs of conjugate state variables in adjacent components. The resulting 'Kirchhoffian network' model (= system macromodel) is governed by generalized mesh rules and node rules for each pair of conjugate state variables and, thus, constitutes a natural approach for system simulation, because the coupling between the different energy and signal domains is governed by the basic conservation laws for energy, particle numbers, mass, charge, etc.
One should note that only one of two conjugate quantities is determined inside the block, whereas the respective conjugate variable can be calculated from the conservation laws governing the network. So the Kirchhoffian network description is, in a certain sense, the natural extension of electric circuit theory, where the branch currents and node voltages are determined by the systematic application of the 'current sum rule' at the nodes and the 'voltage mesh rule' along a closed loop of the network [5], but now with the extension that also physical quantities other than electrical charge are allowed to flow through the network (cf Table 1). Since the system description is based on a meta language, practically any general system simulator with the required HDL interface such as, for instance, Spectre, ELDO or SABER, can be used for implementation.
3. Case study: electrostatically driven micropump 3.1. Operational principle As an illustrative demonstration of the capabilities of the compact modeling method, a macromodel of an electrostatically actuated micropump is derived and discussed. This micropump, as it was developed by Zengerle [3], consists of a pump chamber with an electrostatically actuated membrane as driving element, an inlet and an outlet valve, and externally located connecting tubes (Fig. 1). The flexible membrane is separated from a solid counterelectrode by an airgap and an isolating oxide, which prevents electrical shortage when the membrane touches the counterelectrode. During operation, a voltage is applied between membrane and counterelectrode. The resulting attractive force between the two electrodes causes underpressure in the pump chamber, thus forcing the inlet valve to open. When the voltage is switched off again, the membrane gets released from the counterelectrode, causing a pressure increase in the chamber and forcing the outlet valve to open and the inlet valve to close. In order to obtain maximum membrane deflection and, consequently, maximum pumped volume, the applied voltage must be high enough to make the membrane touch the counterelectrode. membrane
counterelectrode
Table 1 Driving forces and resulting fluxes in a Kirchhoffian network Domain
'Across' -quantity
'Through'-quantity
Electrical Mechanical Thermal Fluidic Chemical Any
Voltage Velocity Temperature Pressure Chemical potential Driving force
Electric current Force Heat flow Mass flow Particle flow Flux
inlet and outlet tubes
Fig. 1. Schematic view of a micropump fabricated by silicon-based microstructure technology.
793
P. Voigt et al./Microelectronics Journal 29 (1998) 791-797
3.2. System analysis For the analysis of the system operation, an appropriate set of state variables has to be determined, and the entire system has to be decomposed into its basic constituent parts, which then will be described by compact models. The natural set of conjugate state variables for fluidic systems consists of pressure and mass flow rate. For the application considered, the fluid .:an be assumed to be incompressible. Therefore, we may introduce the volume flow rate as a variable equivalent to the mass flow rate. The system is divided into the following parts: tubes; valves; and pump chamber with electrostatic membrane drive. The interfaces which separate the valves from the membrane are displayed in Fig. 2. The assumptions underlying this block segmentation are as follows: the pressure distribution within the chamber is spatially uniform, and the flow rates through the interfaces, given as the integrals of the flow density along the respective interface, are adequate to describe the operation of the system parts. This implies that all local fluid redistributions inside the pump chamber occur quasi-statically (i.e. in:~tantaneously on the time scale considered), thus maintaining uniform pressure conditions. The resulting blocks and their interconnections, as they appear in an equivalent 'generalized Kirchhoffian network' description, are shown in Fig. 3. The interfaces separating the valves from the connecting tubes are defined correspondingly.
as a function ofy and Ap. This relation may be introduced as an empirical function or extracted from a fluid-mechanical numerical simulation using, for example, finite element analysis. Fig. 4 (left part) shows the final relation between the flow rate and the pressure difference as obtained from FEM simulation, after the internal variable y has been eliminated. When Ap changes from positive (valve open) to negative values (valve closed), the membrane hits the valve seat. This causes a discontinuity in the slope of the function y(Ap). The volume Vdis displaced by the moving membrane can be expressed as a function of y. Hence, the time derivative of Vdis contributes to the flow rate through the valve by a displacement flow according to: dVdis
dt
dVdis dy d(Ap) dy d(Ap) dt
Due to the nonsmooth behavior ofy(Ap), the corresponding fluidic capacitance Cf = dVdis/d(Ap) (see Fig. 4, right part) exhibits a discontinuity, when the membrane touches the valve seat. One should note that this model describes only the quasistatic valve behavior. Transient effects due to the inertia of the membrane (as described in [7]) are not taken into account. Therefore the application of the model is limited to operation frequencies well below the resonance frequency of the membrane, which amounts to ca 1 kHz in the case considered.
4.2. Tubes 4.
Basiccompactmodels
4.1. Valves According to the system decomposition introduced above, the compact model of the valves has two terminals which represent the inlet and the outlet. The function of the compact model is to establish the relationship between the pressure values at the terminals and the resulting flow rates passing through the terminals. To this end, first the membrane center deflection y is calculated as a function of the applied pressure difference Ap between input and output terminals [6]. Here, y appears as a quasi-static internal (intermediate) state variable, which is eliminated afterwards. Next, the tlow rate through the valve is modeled membrane
counterelectrode
For a given pressure difference along a tube, the resulting flow rates at each end of the tube need to be modeled. Our approach consists in dividing the tube into small slices of equal length A1, with each of the slices being described by the following elementary subcircuit (in the sense of generalized Kirchhoffian network theory):
Here the subcircuit elements represent generalized resistors, capacitors and inductors, understood in the sense that these
membrane]
valvel Fig. 2. Block decomposition of a micropump in valves, pump chamber and membrane.
valve2
Fig. 3. Part of the equivalent generalized Kirchhoffian network describing the pump chamber with membrane drive and the valves.
P. Voigt et al./Microelectronics Journal 29 (1998) 791-797
794
elements constitute a symbolic notation of a linear relation between an across-quantity A and the conjugate throughquantity T or their respective time-derivatives: dA dT
A=RT, T=C~, A=L~.
The exterior nodes attached to 'ground' represent the environment outside the tube (i.e. the ambient pressure). In the above-displayed subcircuit, Rs stands for the static flow resistance derived under the assumption of purely laminar flow, which is a good approximation in view of the small tube diameter. L represents the inertia of the fluid, and C accounts for the volume change caused by the elastic property of the tube material. Rd models the local dynamic damping originating from the viscosity of the fluid. This approach resembles very much that of an electrical transmission line model, apart from that Rd is put in series to C and not in parallel, because there is no static leakage flow through the tube wall. According to [8], an alternative model of a tube slice could be expressed as follows:
R1
delay
Whether or not this second approach is more efficient for numerical simulation depends sensitively on the actual implementation in the simulation tool used. Concerning the approximation error introduced by the discretization of the tube into slices, a practical estimate of the maximum length of the tube slices is given by 2~1 < v/lOfup,where v is the pressure propagation velocity of the fluid in the tube, andfup is the highest frequency present in the oscillations of the coupled system of fluid and tube wall. The pressure pulse propagation in a tube as obtained from our model is shown in Fig. 5. The pulse travels through the tube and gets reflected at the end. The phaseshift by 180° results from the 'open'-condition maintained by forced constant pressure at the end of the tube. 30
.............................................
r ...................
4.3. Membrane The membrane model has three terminals, one belonging to the fluidic domain and two referring to the electrical domain. At the two electrical terminals the control voltage between membrane and counterelectrode is applied, whereas the fluidic terminal couples the membrane to the two valves (Fig. 3). The flow rate into the fluidic terminal is the time-derivative of the displaced volume below the membrane. The electrical function of the membrane is modeled as a condensor the capacitance of which is determined by the deflection of the membrane. A fairly simple and yet accurate model of the geometric shape of the displaced membrane is based on the assumption that the displacement surface has approximately cylindrical symmetry, generated by the rotation of a 'deflection curve' around an axis situated in the center of the membrane. As long as the membrane does not touch the fixed counterelectrode, the deflection curve may be approximated by a parabola. Thus the membrane deflection can be parametrized by a single state variable, the center deflection y. Since the membrane thickness is much larger than the displacement, we arrive at a linear relation between deflection and applied force. Once the membrane touches the counterelectrode, the touching segment is excluded from the force balance, resulting in a reduced effective membrane length l (see Fig. 6). So 1 is the second internal state variable introduced in the model. For deflections smaller than the airgap between membrane and counterelectrode, l is constant, equal to the geometrical membrane length. When the membrane has snapped down and hit the counterelectrode, y stays constant, but now l decreases with the applied pressure, following approximately a cubic relationship. The electrostatic attractive force between membrane and counterelectrode is approximated by an array of infinitesimal plate capacitors in parallel, using an weighted average distance between the electrodes. This is justified by the large aspect ratio between lateral and vertical dimensions. A correction for the real three-dimensional membrane geometry
i
m ~20-
i
i
i .............
i ............................
i ............
o
0 -0.4
-0.2
0.0
pressure (bar)
0.2
0.4
0.00 . . . . . . . . . ' . . . . . . . . . ' . . . . . . . . . . . . . . . . -0.4 -0.2 0.0 0.2
0.4
pressure (bar)
Fig. 4. Static characteristics o f the valves: (a) mass flow rate a n d (b) fluidic c a p a c i t a n c e vs applied pressure difference.
P. Voigt et al./Microelectronics Journal 29 (1998) 791-797 1.6-
795
14-_ ........ ................... ~................... •. . . . . . 1._~_ input,]
iiiSiiiiiiiiiiii? iiiiiiiiiiii? iiiiiiill i]
1.4-
i ,,u,l,u,i 12 ........................................................................ ! [
I
10-
1.2~
.............................................................................
~1.0 ~0.8~"0.6-
i
"~ 8- .................. ~ ............ . ............................................... ,
ii!i!i!!ili !!
.................. "~........ ~............................................ ~4o
,
0.4 .......... ~...........,;~-'*:l:............................................
0.2 .......................
'l'l II II ~
1
';, ................................................
i
',.u
0.10
i
0.12
0.14
2
]
0.16
0 -2-
, , , , i ...........................
0.10
0.12
time (s)
0.14
0.16
time (s)
Fig. 5. Transient pulse response of a tube. In the left figure, the pressure is displayed at the input side and at positions 10%, 50% and 90% of the total tube length. The right figure shows the flow rate at the input and output sides of the tube.
is obtained as follow,;: the membrane is divided into three regions Fig. 7. The displaced volume beneath the regions At and Ap is calculated from the above-discussed twodimensionally rotated displacement curve, whereas the displaced volume associated with the comer regions Ao is calculated from y and l and a fitted comer-fill factor. The static characteristics (deflected volume vs applied pressure) is shown in Fig. 8 (left part). The model reproduces accurately the hysteresis originated from the snap-down of the membrane, which occurs when the pressure is increased with sufficiently high voltage applied between the electrodes. Even though some details of this effect (for instance that the membrane suddenly snaps into another deflection shape when a volume change is forced) have been neglected in the compact model, the results are accurate enough to co~rrectly reproduce the membrane behavior in a system model (see right part of Fig. 8).
are connected to constant pressure sources, which are zero in our case. The ground nodes represent the pressure of the environment. The voltage source is a symbolic notation which actually comprises all the attached electric control circuitry, built up by a network of microelectronic compact models. The response of the pump to a 200 V voltage pulse is shown in the left part of Fig. 10. The right part of this figure shows the flow rate and the pumped volume (time integral of flowrate), with the pump operated at a driving voltage of 200 V and a frequency of 500 Hz. It is interesting to note that it takes c a 50 msec for the pump to reach the full (steady-state) pump rate. The results are in good agreement with the measured data given in [3]. In particular, there is evidence for pressure oscillations arising from the inertia of the fluid in the tubes. It is important to note that the system behavior is strongly influenced or even dominated by the 'parasitics', that is, the inertia of the fluid and the elasticity of the tube walls and the valve membrane.
5. System synthesis and simulation results Using the above compact models, a macromodel of the entire micropump can be built as shown in Fig. 9. The tubes
6. Conclusions Using the concept of tailored compact modeling, it is possible to arrive at a generalized Kirchhoffian network
i_membran
I/2
I/2 Ap
Ac .
~
.........i........i........ .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
Xt
Fig. 6. Membrane deflection curve. The segment touching the counterelecrode from - x t to xt is excluded from the force balance, reducing the effective membrane length to small values of 1.
Fig. 7. Top view of a membrane, indicating different regions: At = counterelectrode touched, Ap = parabolic approximation of the deflection curve, Ac = a corner fill fit applies.
796
P. Voigt et al./Microelectronics Journal 29 (1998) 791-797
!
i' isOV
i......
~.........
!
i
i 0v
!
_
0-
~
.o
~-100-
1: /~/ .....
1
'
i
!
°''°''
tl
::
// i
~0V measured
~
-
i
-200
I'
~J--°l~0Vf°rw~d
iL--" 150V reverse
.oo.o,
3ooo2;1o ...... o'.o...... 2;.o ..... '~.o ..... '~.o' ..... ' ~ o " -0.2
0.0
0.2
volume[nl]
0.4
pressure (bar) Fig. 8. Static characeristicsof a membrane.The left figureshows volumedeflectionvs. pressure.The right figureshowsthe resultingpressurewhen volumeis forced. A voltage of 150 V is applied to the membrane,causing the snap effect. The differencebetween the compactmodel and measurementis due to the parabolic approximationof the membranebending curve. description of electrofluidic microsystems that constitutes the proper framework for the numerical simulation and analysis of all effects relevant to the operation of the microsystem and its individual components. Following the methodology proposed for describing and implementing models on the base of VHDL-AMS it is possible to simulate any microfluidic component by the use of any system simulator which is capable of understanding VHDL-AMS. Since the underlying coupling of physical transducer effects is consistently accounted for in the electrofluidic macromodels, including the desired and undesired (parasitic) effects occurring in all the constituent system components (as, for example, the inertia and the damping of the fluid in the tubes or the discontinuous capacitive behavior of flap valves), the system behavior can be correctly reproduced by the simulation. The practicality and efficiency of this approach with reference to a micropump fabricated by the use of silicon microstructure technology has been demonstrated.
Introducing an appropriate decomposition of the micropump system in functional blocks, all parts of the micropump are modeled using a description in VHDLAMS. The resulting macromodel of the entire micropump covers all relevant functional details and can be easily passed to any conventional electronic circuit simulator, yielding a quick and accurate visualization and prediction of the pump operation. All the typical effects which were found in measured characteristics and which were previously analyzed using a special tool exclusively dedicated to micropumps [3], could be obtained and confirmed in a much more general context and approach. As a consequence and future perspective, the nonelectrical parts as well as the electronic circuitry in general microsystems may be simultaneously simulated using an uniform and consistent description language and computer-aided design environment, thus allowing the easy design and optimization of even complex microsystems.
F®q I U
-.,
I
i> valvel
I
,,.
i>-.,I valve2
Fig. 9. GeneralizedKirchhoffiannetwork macromodelof the entire mircopump.
t
P. Voigt et al./Microelectronics Journal 29 (1998) 791-797
797
0.2 ,~ o o , 1 ............. i .............. i.................. ~i ........... i ................. i ............... I
~_~0.0 • ~" "°-~4 ~i
i i ~ i ', I ............ t .................. ! .................. i .................. i .................. i .................. I
-o.o~].
=.....
i .........
0.00
0.02
i .........
i .........
0.04
0.06
i .........
0.08
i .........
I
0.I0
0.12 -0.2
o.o6~. . . . . . . . f................... i]~-....... 7 - ' ~ . p_membrane] -~ 0 0 , t ............... i ................ ~................ i I ~ - i e0o2t -o.o2
................
............... }.................. }...............
}l-'v
i
..........
i ...............
i ..............
0.00
]
r
i ..................
i ................
0.02
0 . 0 6 ~T . . . . . . . . .
~..........
0.04
r- . . . . . . .
r
0.06
1- . . . . . . . .
i
.......... i ............... :
o.og
7 .......
!
i .................
O.lO
~--
q
0.12
o m~etp
,40t'°iiiiiiiiliiiiiiiiiiii12ili1711ili ~ ~80-
0.o2~-................ i!.................. .o.o2
i ................. i~-............. i.................. i ................. 4 !
i
]
............. !.................! ................. I................. I.................. I................I
............... i.................. i.................. i.................. i................ i -0.061..~
0.00
.,,i
.........
0,02
i .........
0.04
i .........
0.06
time (s)
i .........
0.08
i .........
0,10
~ 60~,400.
i
0.12
0.00
i .
.
. I 0.02
: .
.
.Y" .
.
.
I 0.04
I 0.06
time
i .
.
. I 0.08
ii ~i I 0.10
(s)
Fig. 10. Transient pulse response of the pump. The pressure at the inlet, at the outlet and in the pump chamber near the membrane is shown in the left figure, with a single pulse of 200 V applied from 0 to 60 ms. The flow rate and pumped volume are displayed at the right, with the pump pulsed at 200 V and 500 Hz.
References [1] S'D" Senturia' CAD f°r micmelectr°mechanical systems' Digest °fTech" Papers of Transducers ' 95, Eurosensors IX, Stockholm, 1995, pp. 5-8. [2] G. Wachutka, Tailored modeling: a way to the "virtual microtransducer fab"? Sensors and Actuators A 46-47 (1995) 603-612. [3] R. Zengerle, Mikro-Membranpumpen als Komponenten for MikroFluidsysteme, Verlag Shaker, Aachen, 1994. [4] R. Zengede, J. Ulrich, S. Kluge, M. Richter, A. Richter, A bidirectional silicon micropump, Sensors and Actuators A 50 (1995) 81-86.
[5] R. Boite, J. Neirynck, Theorie des Reseaux de Kirchhoff, Editions Georgi, St-Saphorin, Switzerland, 1976. [6] S. Timoshenko, Theory of Plates and Shells. McGraw-Hill, New York, 1959. [7] J. Ulrich, R. Zengerle, Static and dynamic flow simulation of a KOHetched microvalve using the finite-element method, Sensors and Actuators A 53 (1996) 379-385. [8] R. Gupta, S. Kim, L. Pileggi, Domain characterization of transmission line models and analyses, IEEE Trans. Computer-Aided Design 15 (Feb.) (1996) 184-193•
Microelectronics Journal 29 (1998) 791-797
Microfluidic system modeling using VHDL-AMS and circuit simulation P. Voigt, G. Schrag, G. Wachutka Institute.for Physics of Electrotechnology, Technical University of Munich, Arcisstrafle 21, 80290 Munich, Germany
Abstract
It is demonstrated that the operation of electrofluidic microsystems such as micropumps can be efficiently modeled by standard electrical circuit simulators, using the analog hardware description language (VHDL-AMS) to describe the models. As a representative example, compact models were developed for all constituent parts of an electromechanical micropump fabricated by silicon-based microstmcture technology, including ele.ctrostatically actuated membrane drives, membrane valves and tubes together with the power supply and control circuitry. With the resulting macromodel of the micropump, the static characteristics as well as the dynamic behavior is correctly reproduced. Thus, with this approach it is possible to include fluidmechanical system components straightforwardly in conventional microelectronic design environments, allowing for easy full system analysis and optimization of even complex microsystems. © 1998 Elsevier Science Ltd. All rights reserved.
1. Introduction
The rapid progress in microstructure technology allows the fabrication of increasingly complex microsystems, thereby raising a high demand for powerful simulation tools for the design and optimization of such systems [1]. However, the development of special purpose software packages tends to be quite expensive and time consuming. On the other hand, the electronic industry has driven the development of efficiently performing design systems for microdevices. Since many aspects of microsystem design are closely related to the silicon-based integrated circuit standard technologies, notably with a view to the integration of electronics and transducers in 'smart systems', it seems advantageous to simuhtte the nonelectrical parts of a microsystem using conventional design packages available from the semiconductor industry. The appropriate way of modeling at system level is the method of tailored compact modeling as described in [2]. In this work, the practicability of this idea is demonstrated with reference to electrofluidic microsystems. It is typical of such microsystems that they are equipped with one or several electromechanical driving units such as electrostatically, electrothermally or piezoelectrically actuated micropumps, passive fluidic components such as static or dynamic valves and tubes, sensor elements for control and analysis, and a co-integrated or hybrid coupled electronic circuitry for power supply, signal conditioning and system control including (self-)test and (self-)calibration as well as error detection and compensation. By their nature, the constituent components link different energy and signal domains such 0026-2692/98/$19.00 © 1998 Elsevier Science Ltd. All rights reserved PII S 0 0 2 6 - 2 6 9 2 ( 9 7 ) 0 0 0 9 3 - 1
as mechanical, fluidic, thermal, electrical and other physical or chemical quantities. Typical applications include, among others, chemical microanalyzers, mass flowmeters, fuel drop generators for automobile heaters, ink jet print heads or microdosed drug infusion. The natural description of the involved transducer effects is provided by the concepts of continuum theory, which couples the relevant mechanical, electrodynamical and thermal field quantities in terms of a (usually quite complex) system of partial differential equations [2]. However, it is hardly feasible to solve the model equations in their continuous form for realistic microsystems in their full size and complexity. Therefore, to keep the computational effort for such problems within acceptable limits, the proper and adequate level of simulation is analog circuit simulation based on a network theory which, on the one hand, preserves the basic conservation laws expressed in the continuum models and, on the other hand, reduces the set of variables to a number that is still manageable with one of today's standard circuit simulation tools as they are routinely used in the microelectronic industry. In this approach, the microsystem is partitioned into functional blocks, which interact with each other as constituent parts of a flux-conserving ('Kirchhoffian') network. For each of them, a compact model with only few degrees of freedom is formulated, thus defining a macromodel for all parts of the microsystem. Then VHDL-AMS (analog hardware description language) is used to code the models. This has the advantage of providing a simulator-independent model code, since VHDLAMS is a quasi-standard, accepted as input language by many electronic simulation systems.
792
P. Voigt et al./Microelectronics Journal 29 (1998) 791-797
In the following, the method of tailored compact modeling (Section 2) is briefly introduced. Then the general principles basic to this approach are exemplified with reference to an electrofluidic micropump, the detailed description of which is given in [3],[4]. The compact models of the micropump components are then derived and discussed (Sections 3 and 4), and finally the transient operation of the micropump by means of time-dependent system simulation is studied (Section 5).
2. Physically-based compact modeling System models of complex microstructures can be constructed by partitioning the full system into simpler subsystems. The interfaces between interacting subsystems must be chosen in such a way that the exchange of energy, mass, charge and other physical quantities can be described by pairs of conjugate state variables. A systematic method of selecting such pairs of variables is provided by the basic principles of thermodynamics [2]. A pair of conjugate variables consists of the flux of a physical quantity and an associate driving force. Examples are given in Table 1. In the context of VHDL-AMS, the driving force is called 'across quantity' and the flux 'through quantity'. A subsystem is called 'block' and is defined by the number and nature of its terminals. The terminals allow for the exchange of flux quantities across subsystem boundaries, concentrated ('lumped') at single points instead of the distributed flow of the respective field quantities. The intemal behavior of each of the blocks is described by a physically-based compact model extracted from continuous field theory and/or from experimental data. In practice, the compact models are formulated using an appropriate analog hardware description language such as, for instance, VHDL-AMS. The full microsystem is synthesized by linking the individual HDL-specifled compact models of the constituent blocks through flux-conserving interface conditions for the respective pairs of conjugate state variables in adjacent components. The resulting 'Kirchhoffian network' model (= system macromodel) is governed by generalized mesh rules and node rules for each pair of conjugate state variables and, thus, constitutes a natural approach for system simulation, because the coupling between the different energy and signal domains is governed by the basic conservation laws for energy, particle numbers, mass, charge, etc.
One should note that only one of two conjugate quantities is determined inside the block, whereas the respective conjugate variable can be calculated from the conservation laws governing the network. So the Kirchhoffian network description is, in a certain sense, the natural extension of electric circuit theory, where the branch currents and node voltages are determined by the systematic application of the 'current sum rule' at the nodes and the 'voltage mesh rule' along a closed loop of the network [5], but now with the extension that also physical quantities other than electrical charge are allowed to flow through the network (cf Table 1). Since the system description is based on a meta language, practically any general system simulator with the required HDL interface such as, for instance, Spectre, ELDO or SABER, can be used for implementation.
3. Case study: electrostatically driven micropump 3.1. Operational principle As an illustrative demonstration of the capabilities of the compact modeling method, a macromodel of an electrostatically actuated micropump is derived and discussed. This micropump, as it was developed by Zengerle [3], consists of a pump chamber with an electrostatically actuated membrane as driving element, an inlet and an outlet valve, and externally located connecting tubes (Fig. 1). The flexible membrane is separated from a solid counterelectrode by an airgap and an isolating oxide, which prevents electrical shortage when the membrane touches the counterelectrode. During operation, a voltage is applied between membrane and counterelectrode. The resulting attractive force between the two electrodes causes underpressure in the pump chamber, thus forcing the inlet valve to open. When the voltage is switched off again, the membrane gets released from the counterelectrode, causing a pressure increase in the chamber and forcing the outlet valve to open and the inlet valve to close. In order to obtain maximum membrane deflection and, consequently, maximum pumped volume, the applied voltage must be high enough to make the membrane touch the counterelectrode. membrane
counterelectrode
Table 1 Driving forces and resulting fluxes in a Kirchhoffian network Domain
'Across' -quantity
'Through'-quantity
Electrical Mechanical Thermal Fluidic Chemical Any
Voltage Velocity Temperature Pressure Chemical potential Driving force
Electric current Force Heat flow Mass flow Particle flow Flux
inlet and outlet tubes
Fig. 1. Schematic view of a micropump fabricated by silicon-based microstructure technology.
793
P. Voigt et al./Microelectronics Journal 29 (1998) 791-797
3.2. System analysis For the analysis of the system operation, an appropriate set of state variables has to be determined, and the entire system has to be decomposed into its basic constituent parts, which then will be described by compact models. The natural set of conjugate state variables for fluidic systems consists of pressure and mass flow rate. For the application considered, the fluid .:an be assumed to be incompressible. Therefore, we may introduce the volume flow rate as a variable equivalent to the mass flow rate. The system is divided into the following parts: tubes; valves; and pump chamber with electrostatic membrane drive. The interfaces which separate the valves from the membrane are displayed in Fig. 2. The assumptions underlying this block segmentation are as follows: the pressure distribution within the chamber is spatially uniform, and the flow rates through the interfaces, given as the integrals of the flow density along the respective interface, are adequate to describe the operation of the system parts. This implies that all local fluid redistributions inside the pump chamber occur quasi-statically (i.e. in:~tantaneously on the time scale considered), thus maintaining uniform pressure conditions. The resulting blocks and their interconnections, as they appear in an equivalent 'generalized Kirchhoffian network' description, are shown in Fig. 3. The interfaces separating the valves from the connecting tubes are defined correspondingly.
as a function ofy and Ap. This relation may be introduced as an empirical function or extracted from a fluid-mechanical numerical simulation using, for example, finite element analysis. Fig. 4 (left part) shows the final relation between the flow rate and the pressure difference as obtained from FEM simulation, after the internal variable y has been eliminated. When Ap changes from positive (valve open) to negative values (valve closed), the membrane hits the valve seat. This causes a discontinuity in the slope of the function y(Ap). The volume Vdis displaced by the moving membrane can be expressed as a function of y. Hence, the time derivative of Vdis contributes to the flow rate through the valve by a displacement flow according to: dVdis
dt
dVdis dy d(Ap) dy d(Ap) dt
Due to the nonsmooth behavior ofy(Ap), the corresponding fluidic capacitance Cf = dVdis/d(Ap) (see Fig. 4, right part) exhibits a discontinuity, when the membrane touches the valve seat. One should note that this model describes only the quasistatic valve behavior. Transient effects due to the inertia of the membrane (as described in [7]) are not taken into account. Therefore the application of the model is limited to operation frequencies well below the resonance frequency of the membrane, which amounts to ca 1 kHz in the case considered.
4.2. Tubes 4.
Basiccompactmodels
4.1. Valves According to the system decomposition introduced above, the compact model of the valves has two terminals which represent the inlet and the outlet. The function of the compact model is to establish the relationship between the pressure values at the terminals and the resulting flow rates passing through the terminals. To this end, first the membrane center deflection y is calculated as a function of the applied pressure difference Ap between input and output terminals [6]. Here, y appears as a quasi-static internal (intermediate) state variable, which is eliminated afterwards. Next, the tlow rate through the valve is modeled membrane
counterelectrode
For a given pressure difference along a tube, the resulting flow rates at each end of the tube need to be modeled. Our approach consists in dividing the tube into small slices of equal length A1, with each of the slices being described by the following elementary subcircuit (in the sense of generalized Kirchhoffian network theory):
Here the subcircuit elements represent generalized resistors, capacitors and inductors, understood in the sense that these
membrane]
valvel Fig. 2. Block decomposition of a micropump in valves, pump chamber and membrane.
valve2
Fig. 3. Part of the equivalent generalized Kirchhoffian network describing the pump chamber with membrane drive and the valves.
P. Voigt et al./Microelectronics Journal 29 (1998) 791-797
794
elements constitute a symbolic notation of a linear relation between an across-quantity A and the conjugate throughquantity T or their respective time-derivatives: dA dT
A=RT, T=C~, A=L~.
The exterior nodes attached to 'ground' represent the environment outside the tube (i.e. the ambient pressure). In the above-displayed subcircuit, Rs stands for the static flow resistance derived under the assumption of purely laminar flow, which is a good approximation in view of the small tube diameter. L represents the inertia of the fluid, and C accounts for the volume change caused by the elastic property of the tube material. Rd models the local dynamic damping originating from the viscosity of the fluid. This approach resembles very much that of an electrical transmission line model, apart from that Rd is put in series to C and not in parallel, because there is no static leakage flow through the tube wall. According to [8], an alternative model of a tube slice could be expressed as follows:
R1
delay
Whether or not this second approach is more efficient for numerical simulation depends sensitively on the actual implementation in the simulation tool used. Concerning the approximation error introduced by the discretization of the tube into slices, a practical estimate of the maximum length of the tube slices is given by 2~1 < v/lOfup,where v is the pressure propagation velocity of the fluid in the tube, andfup is the highest frequency present in the oscillations of the coupled system of fluid and tube wall. The pressure pulse propagation in a tube as obtained from our model is shown in Fig. 5. The pulse travels through the tube and gets reflected at the end. The phaseshift by 180° results from the 'open'-condition maintained by forced constant pressure at the end of the tube. 30
.............................................
r ...................
4.3. Membrane The membrane model has three terminals, one belonging to the fluidic domain and two referring to the electrical domain. At the two electrical terminals the control voltage between membrane and counterelectrode is applied, whereas the fluidic terminal couples the membrane to the two valves (Fig. 3). The flow rate into the fluidic terminal is the time-derivative of the displaced volume below the membrane. The electrical function of the membrane is modeled as a condensor the capacitance of which is determined by the deflection of the membrane. A fairly simple and yet accurate model of the geometric shape of the displaced membrane is based on the assumption that the displacement surface has approximately cylindrical symmetry, generated by the rotation of a 'deflection curve' around an axis situated in the center of the membrane. As long as the membrane does not touch the fixed counterelectrode, the deflection curve may be approximated by a parabola. Thus the membrane deflection can be parametrized by a single state variable, the center deflection y. Since the membrane thickness is much larger than the displacement, we arrive at a linear relation between deflection and applied force. Once the membrane touches the counterelectrode, the touching segment is excluded from the force balance, resulting in a reduced effective membrane length l (see Fig. 6). So 1 is the second internal state variable introduced in the model. For deflections smaller than the airgap between membrane and counterelectrode, l is constant, equal to the geometrical membrane length. When the membrane has snapped down and hit the counterelectrode, y stays constant, but now l decreases with the applied pressure, following approximately a cubic relationship. The electrostatic attractive force between membrane and counterelectrode is approximated by an array of infinitesimal plate capacitors in parallel, using an weighted average distance between the electrodes. This is justified by the large aspect ratio between lateral and vertical dimensions. A correction for the real three-dimensional membrane geometry
i
m ~20-
i
i
i .............
i ............................
i ............
o
0 -0.4
-0.2
0.0
pressure (bar)
0.2
0.4
0.00 . . . . . . . . . ' . . . . . . . . . ' . . . . . . . . . . . . . . . . -0.4 -0.2 0.0 0.2
0.4
pressure (bar)
Fig. 4. Static characteristics o f the valves: (a) mass flow rate a n d (b) fluidic c a p a c i t a n c e vs applied pressure difference.
P. Voigt et al./Microelectronics Journal 29 (1998) 791-797 1.6-
795
14-_ ........ ................... ~................... •. . . . . . 1._~_ input,]
iiiSiiiiiiiiiiii? iiiiiiiiiiii? iiiiiiill i]
1.4-
i ,,u,l,u,i 12 ........................................................................ ! [
I
10-
1.2~
.............................................................................
~1.0 ~0.8~"0.6-
i
"~ 8- .................. ~ ............ . ............................................... ,
ii!i!i!!ili !!
.................. "~........ ~............................................ ~4o
,
0.4 .......... ~...........,;~-'*:l:............................................
0.2 .......................
'l'l II II ~
1
';, ................................................
i
',.u
0.10
i
0.12
0.14
2
]
0.16
0 -2-
, , , , i ...........................
0.10
0.12
time (s)
0.14
0.16
time (s)
Fig. 5. Transient pulse response of a tube. In the left figure, the pressure is displayed at the input side and at positions 10%, 50% and 90% of the total tube length. The right figure shows the flow rate at the input and output sides of the tube.
is obtained as follow,;: the membrane is divided into three regions Fig. 7. The displaced volume beneath the regions At and Ap is calculated from the above-discussed twodimensionally rotated displacement curve, whereas the displaced volume associated with the comer regions Ao is calculated from y and l and a fitted comer-fill factor. The static characteristics (deflected volume vs applied pressure) is shown in Fig. 8 (left part). The model reproduces accurately the hysteresis originated from the snap-down of the membrane, which occurs when the pressure is increased with sufficiently high voltage applied between the electrodes. Even though some details of this effect (for instance that the membrane suddenly snaps into another deflection shape when a volume change is forced) have been neglected in the compact model, the results are accurate enough to co~rrectly reproduce the membrane behavior in a system model (see right part of Fig. 8).
are connected to constant pressure sources, which are zero in our case. The ground nodes represent the pressure of the environment. The voltage source is a symbolic notation which actually comprises all the attached electric control circuitry, built up by a network of microelectronic compact models. The response of the pump to a 200 V voltage pulse is shown in the left part of Fig. 10. The right part of this figure shows the flow rate and the pumped volume (time integral of flowrate), with the pump operated at a driving voltage of 200 V and a frequency of 500 Hz. It is interesting to note that it takes c a 50 msec for the pump to reach the full (steady-state) pump rate. The results are in good agreement with the measured data given in [3]. In particular, there is evidence for pressure oscillations arising from the inertia of the fluid in the tubes. It is important to note that the system behavior is strongly influenced or even dominated by the 'parasitics', that is, the inertia of the fluid and the elasticity of the tube walls and the valve membrane.
5. System synthesis and simulation results Using the above compact models, a macromodel of the entire micropump can be built as shown in Fig. 9. The tubes
6. Conclusions Using the concept of tailored compact modeling, it is possible to arrive at a generalized Kirchhoffian network
i_membran
I/2
I/2 Ap
Ac .
~
.........i........i........ .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
Xt
Fig. 6. Membrane deflection curve. The segment touching the counterelecrode from - x t to xt is excluded from the force balance, reducing the effective membrane length to small values of 1.
Fig. 7. Top view of a membrane, indicating different regions: At = counterelectrode touched, Ap = parabolic approximation of the deflection curve, Ac = a corner fill fit applies.
796
P. Voigt et al./Microelectronics Journal 29 (1998) 791-797
!
i' isOV
i......
~.........
!
i
i 0v
!
_
0-
~
.o
~-100-
1: /~/ .....
1
'
i
!
°''°''
tl
::
// i
~0V measured
~
-
i
-200
I'
~J--°l~0Vf°rw~d
iL--" 150V reverse
.oo.o,
3ooo2;1o ...... o'.o...... 2;.o ..... '~.o ..... '~.o' ..... ' ~ o " -0.2
0.0
0.2
volume[nl]
0.4
pressure (bar) Fig. 8. Static characeristicsof a membrane.The left figureshows volumedeflectionvs. pressure.The right figureshowsthe resultingpressurewhen volumeis forced. A voltage of 150 V is applied to the membrane,causing the snap effect. The differencebetween the compactmodel and measurementis due to the parabolic approximationof the membranebending curve. description of electrofluidic microsystems that constitutes the proper framework for the numerical simulation and analysis of all effects relevant to the operation of the microsystem and its individual components. Following the methodology proposed for describing and implementing models on the base of VHDL-AMS it is possible to simulate any microfluidic component by the use of any system simulator which is capable of understanding VHDL-AMS. Since the underlying coupling of physical transducer effects is consistently accounted for in the electrofluidic macromodels, including the desired and undesired (parasitic) effects occurring in all the constituent system components (as, for example, the inertia and the damping of the fluid in the tubes or the discontinuous capacitive behavior of flap valves), the system behavior can be correctly reproduced by the simulation. The practicality and efficiency of this approach with reference to a micropump fabricated by the use of silicon microstructure technology has been demonstrated.
Introducing an appropriate decomposition of the micropump system in functional blocks, all parts of the micropump are modeled using a description in VHDLAMS. The resulting macromodel of the entire micropump covers all relevant functional details and can be easily passed to any conventional electronic circuit simulator, yielding a quick and accurate visualization and prediction of the pump operation. All the typical effects which were found in measured characteristics and which were previously analyzed using a special tool exclusively dedicated to micropumps [3], could be obtained and confirmed in a much more general context and approach. As a consequence and future perspective, the nonelectrical parts as well as the electronic circuitry in general microsystems may be simultaneously simulated using an uniform and consistent description language and computer-aided design environment, thus allowing the easy design and optimization of even complex microsystems.
F®q I U
-.,
I
i> valvel
I
,,.
i>-.,I valve2
Fig. 9. GeneralizedKirchhoffiannetwork macromodelof the entire mircopump.
t
P. Voigt et al./Microelectronics Journal 29 (1998) 791-797
797
0.2 ,~ o o , 1 ............. i .............. i.................. ~i ........... i ................. i ............... I
~_~0.0 • ~" "°-~4 ~i
i i ~ i ', I ............ t .................. ! .................. i .................. i .................. i .................. I
-o.o~].
=.....
i .........
0.00
0.02
i .........
i .........
0.04
0.06
i .........
0.08
i .........
I
0.I0
0.12 -0.2
o.o6~. . . . . . . . f................... i]~-....... 7 - ' ~ . p_membrane] -~ 0 0 , t ............... i ................ ~................ i I ~ - i e0o2t -o.o2
................
............... }.................. }...............
}l-'v
i
..........
i ...............
i ..............
0.00
]
r
i ..................
i ................
0.02
0 . 0 6 ~T . . . . . . . . .
~..........
0.04
r- . . . . . . .
r
0.06
1- . . . . . . . .
i
.......... i ............... :
o.og
7 .......
!
i .................
O.lO
~--
q
0.12
o m~etp
,40t'°iiiiiiiiliiiiiiiiiiii12ili1711ili ~ ~80-
0.o2~-................ i!.................. .o.o2
i ................. i~-............. i.................. i ................. 4 !
i
]
............. !.................! ................. I................. I.................. I................I
............... i.................. i.................. i.................. i................ i -0.061..~
0.00
.,,i
.........
0,02
i .........
0.04
i .........
0.06
time (s)
i .........
0.08
i .........
0,10
~ 60~,400.
i
0.12
0.00
i .
.
. I 0.02
: .
.
.Y" .
.
.
I 0.04
I 0.06
time
i .
.
. I 0.08
ii ~i I 0.10
(s)
Fig. 10. Transient pulse response of the pump. The pressure at the inlet, at the outlet and in the pump chamber near the membrane is shown in the left figure, with a single pulse of 200 V applied from 0 to 60 ms. The flow rate and pumped volume are displayed at the right, with the pump pulsed at 200 V and 500 Hz.
References [1] S'D" Senturia' CAD f°r micmelectr°mechanical systems' Digest °fTech" Papers of Transducers ' 95, Eurosensors IX, Stockholm, 1995, pp. 5-8. [2] G. Wachutka, Tailored modeling: a way to the "virtual microtransducer fab"? Sensors and Actuators A 46-47 (1995) 603-612. [3] R. Zengerle, Mikro-Membranpumpen als Komponenten for MikroFluidsysteme, Verlag Shaker, Aachen, 1994. [4] R. Zengede, J. Ulrich, S. Kluge, M. Richter, A. Richter, A bidirectional silicon micropump, Sensors and Actuators A 50 (1995) 81-86.
[5] R. Boite, J. Neirynck, Theorie des Reseaux de Kirchhoff, Editions Georgi, St-Saphorin, Switzerland, 1976. [6] S. Timoshenko, Theory of Plates and Shells. McGraw-Hill, New York, 1959. [7] J. Ulrich, R. Zengerle, Static and dynamic flow simulation of a KOHetched microvalve using the finite-element method, Sensors and Actuators A 53 (1996) 379-385. [8] R. Gupta, S. Kim, L. Pileggi, Domain characterization of transmission line models and analyses, IEEE Trans. Computer-Aided Design 15 (Feb.) (1996) 184-193•