Tian-Shiang Yang
email: tsyang@mail.ncku.edu.tw

Chi-Chung Wang

EFFECTS OF ACTUATOR IMPACT ON THE NONLINEAR DYNAMICS OF A VALVELESS PUMPING SYSTEM

Keywords: Valveless pumping, Liebau effect, actuator impact, nonlinear dynamics

HAL is

Introduction

Many fluid-dynamical systems in nature and engineering more or less rely upon valveless pumping to transport the working fluids therein. For instance, blood circulation in the cardiovascular system is maintained to some extent even when the heart's valves fail. [START_REF] Liebau | Über ein ventilloses Pumpprinzip[END_REF][START_REF] Moser | Impedance defined flow: Generalisation of William Harvey's concept of the circulation -370 years later[END_REF] Meanwhile, the embryonic vertebrate heart begins pumping blood long before the development of discernable chambers and valves. [START_REF] Forouhar | The embryonic vertebrate heart tube is a dynamic suction pump[END_REF][START_REF] Manopoulos | Modelling of the blood flow circulation in the human foetus by the end of the third week of gestation[END_REF] (The reader is also referred to Ref. 2 for a discussion of several other biomedical systems that involve valveless pumping.) In microfluidics, valveless "impedance pumps" have been fabricated, and are expected to be particularly suitable for handling sensitive biofluids. [START_REF] Rinderknecht | A valveless micro impedance pump driven by electromagnetic actuation[END_REF] A rather effective mechanism of valveless pumping, namely that of impedance pumps, can be explained as follows. Consider a pliant tube periodically compressed at an asymmetric site from the tube's interfaces to different tubing or reservoirs. Waves then are excited and propagate on the tube. At the tube ends, the waves would partly reflect due to impedance mismatch, thereby building up a mean pressure head that drives the fluid flow in the tube. [START_REF] Hickerson | Experimental study of the behavior of a valveless impedance pump[END_REF][START_REF] Hickerson | On the resonance of a pliant tube as a mechanism for valveless pumping[END_REF] The study of valveless pumping can be traced back to the work of Liebau in 1954, [START_REF] Liebau | Über ein ventilloses Pumpprinzip[END_REF] in which he demonstrated the generation of a mean flow in a valveless openloop mechanical apparatus by periodically compressing a plaint part of the apparatus. The aim of such a demonstration was to support the view that the heart is assisted by other forces to sustain the blood circulation (see Ref. 2 for a brief discussion of the evolution of cardiovascular circulation theories). Liebau's work then inspired a number of experiments and mathematical models for various valveless pumping systems. As these works have been very nicely reviewed by Manopoulos et al., [START_REF] Manopoulos | One-dimensional model of valveless pumping in a closed loop and a numerical solution[END_REF] a detailed account of them will not be presented here. Instead, to give some quick ideas about such works, a few of them are sketched below. For example, in 1963, Mahrenholtz [START_REF] Mahrenholtz | Ein Beitrag zum Förderprinzip periodisch arbeitender, ventilloser Pumpen[END_REF] carried out experiments using a valveless closed-loop apparatus as a mechanical model for the blood circulation system with a heart of valvular disease. A lumped-parameter mathematical model for the system was also derived in that work. In 1978, using the method of characteristics in conjunction with numerical computations, Thomann [START_REF] Thomann | A simple pumping mechanism in a valveless tube[END_REF] studied a simple torus model consisting of an elastic tube and a rigid tube joined together (as described earlier by Liebau [START_REF] Liebau | Aus welchem Grunde bleibt die Blutförderung durch das Herz bei valvulärem Versagen erhalten?[END_REF] and physically reconstructed in a more recent work by Moser et al. [START_REF] Moser | Impedance defined flow: Generalisation of William Harvey's concept of the circulation -370 years later[END_REF] ). A one-dimensional (1-D) model also was constructed and solved numerically by Rath et al. [START_REF] Rath | Der Fördereffekt in ventillosen, elastischen Leitungen[END_REF] to study an open-loop valveless pumping system. Meanwhile, in 1983, Takagi and Saijo [START_REF] Takagi | Study of a piston pump without valves (1st report, on a pipecapacity-sysem with a T-junction)[END_REF] studied, both theoretically and experimentally, the performance of a valveless piston pump under non-resonant conditions. In a follow-up work, Takagi and Takahashi [START_REF] Takagi | Study of a piston pump without valves (2nd report, pumping effect and resonance in a pipe-capacity-system with a T-junction)[END_REF] then modified the piston pump configuration and examined the effects of resonance phenomenon on the pumping characteristics.

In recent years, experimental studies and theoretical works continue to be devoted to understanding the many interesting characteristics of valveless pumping. For example, in 2003, Ottesen [START_REF] Ottesen | Valveless pumping in a fluid-filled closed elastic tube-system: onedimensional theory with experimental validation[END_REF] derived by averaging the Navier-Stokes equations a 1-D nonlinear mathematical model for a fluid-filled torus (like that originally described by Liebau 11 ). The mathematical model was then analyzed partly analytically and partly numerically, and the theoretical findings were validated by some experiments carried out on a physical realization of the torus system. In such experiments, a remarkable unidirectional mean flow in the torus was produced by compressing the soft tube at an asymmetric site with respect to the system's configuration. It was also observed that the size and direction of the mean flow depend on the frequency and specific temporal form of the compression, the elasticity of the tubes, and the compression ratio in a complicated manner.

Meanwhile, aiming to answer a number of important questions regarding the fundamental mechanisms of valveless pumping, Hickerson et al. [START_REF] Hickerson | Experimental study of the behavior of a valveless impedance pump[END_REF] carried out a compre-Effects of Actuator Impact on the Nonlinear Dynamics of a Valveless Pumping System 3 hensive experimental study using an impedance pump. The pump was constructed by connecting an elastic tube to two water reservoirs; in some experiments, the two reservoirs were also connected by an additional rigid tube to form a closed-loop system. The elastic tube was compressed by a set of pinchers moving sinusoidally at a controllable frequency and fixed nominal duty cycle (i.e., the fraction of the pinching period during which the pinchers are separated by a distance smaller than the resting diameter of the tube). Several distinguishing features of impedance-driven flows were then precisely identified, including that the net flow is highly sensitive to the pinching frequency. In particular, it was observed that an open-loop system can also create and sustain a pressure head, and that an elastic material is not a necessary condition for impedance-driven flows. [START_REF] Hickerson | Experimental study of the behavior of a valveless impedance pump[END_REF] Hickerson and Gharib 7 then proposed a simple 1-D "wave-pulse model" accounting for wave amplitude attenuation and reflection at impedance mismatched interfaces in an impedance pump; the model predicted many of the characteristics exhibited by the experiments of Hickerson et al. [START_REF] Hickerson | Experimental study of the behavior of a valveless impedance pump[END_REF] A more sophisticated numerical model involving the use of commercial softwares to solve the Navier-Stokes and structural dynamics equations was also adopted by Avrahami and Gharib, [START_REF] Avrahami | Computational studies of resonance wave pumping in compliant tubes[END_REF] to carry out a more detailed and comprehensive study of the wave dynamics in impedance pumps. It is also noteworthy that Forouhar et al. [START_REF] Forouhar | The embryonic vertebrate heart tube is a dynamic suction pump[END_REF] used the concept of impedance pumps to interpret the images of zebrafish heart tube contractions before valve formation, and proposed that in early cardiac dynamics the pumping action results from suction due to elastic wave propagation in the heart tube. Moreover, inspired by the embryonic heart structure, Loumes et al. [START_REF] Loumes | Resonance pumping in a multilayer impedance pump[END_REF] have introduced a novel concept of "multilayer impedance pump."

Despite the existence of a sizable research literature on valveless pumping, a full understanding of its mechanisms has not yet been obtained. This, to a certain extent, is because valveless systems that differ in configuration, or merely in operation conditions (e.g., frequency of compression), may indeed rely upon different mechanisms to pump fluids. Moreover, in a particular system setup, several different mechanisms of valveless pumping may interact synergistically. It is intuitively clear, however, that for a valveless pumping system to produce any mean flow, a necessary-but not sufficient (to be explained below)-condition is to drive (i.e., compress) the system at an asymmetric site with respect to its configuration.

It is also arguable that appropriate theoretical modeling has been as important as experimental studies in furnishing a better understanding of valveless pumping mechanisms. Broadly speaking, existing mathematical models of valveless pumping systems can be divided into the following three categories: namely lumpedparameter models where the state of the system is represented by a small number of lumped parameters (which, in turn, are governed by a set of coupled ordinary differential equations; see, for example, Refs. 2, 18-20), computational models that solve relevant fluid and structural dynamics equations (e.g., Refs. 3, 8, 15, 21, 22), and the wave-pulse model of Hickerson and Gharib 7 mentioned above. Basically, lumped-parameter models are suitable for systems having a relatively low driving frequency so that the effects of wave propagation are subdominant, whereas the wave-pulse model presumably would work well for systems driven at a relatively high frequency. Also, if all the relevant physical factors are accounted for, computational models supposedly can be used to simulate the system dynamics for a broad range of driving frequencies (but usually at higher computational costs).

Here we would like to point out that, in existing models of valveless pumping, the interaction between the compression actuator and the pliant part of the system in question has often been modeled in rather simple ways. As a matter of fact, in some of the models cited above, the instantaneous deformation of the pliant part of the system is prescribed in a manner not inconsistent (but not absolutely consistent either, as we shall explain in Sect. 4) with the actuator motion (see, for example, Ref. 15), while in some other models the temporal variation of the external pressure at the compression site is specified (e.g., Ref. 22). Meanwhile, in the wave-pulse model of Hickerson and Gharib, 7 a pair of pressure waves of a prescribed waveform are emitted from the compression site at the driving frequency, and then travel in opposite directions. Under suitable conditions, the aforementioned approaches have produced theoretical predictions that can be satisfactorily verified by experiments. This suggests that, in such cases, the actuator and the pliant part of the system may indeed interact in a rather simple way. However, in the present work we shall demonstrate that, in some cases, the collision of the compression actuator and the system's pliant part would be the only cause that produces a mean flow in the system, and therefore deserves more careful modeling.

To that end, a fluid-filled closed-loop system is considered here, which consists of two distensible reservoirs connected by two rigid tubes, with one of the reservoirs (hereafter referred to as the V 0 reservoir) compressed by an actuator at a prescribed frequency; see Fig. 1. It is assumed that the actuator motion specifies the upper bound of the instantaneous volume of the V 0 reservoir. Moreover, during the whole compression cycle, that upper bound is taken to be smaller than the free (i.e., unstressed) volume of the V 0 reservoir, and so the nominal duty cycle of the compression actuator is unity. To suppress wave propagation effects, the driving frequency is assumed to be sufficiently low, and a lumped-parameter model (with constant coefficients) accounting for mass and momentum balance in the system is constructed in Sect. 2.

Based upon such a model, it will be shown in Sect. 3 that the V 0 reservoir may separate from the compression actuator during some part of the compression cycle if the driving frequency of the actuator exceeds a certain threshold value. The actual duty cycle (i.e., the fraction of the compression cycle during which the actuator actually compresses the system's plaint part) would then become less than unity. In fact, only when that happens would a nonzero mean flow be produced in the fluid loop by the synergetic interaction between configurational asymmetry and the nonlinear effects associated with actuator impact. A physically based model for the nonlinear impact effects will then be proposed in Sect. 4 to complete the mathe-Effects of Actuator Impact on the Nonlinear Dynamics of a Valveless Pumping System 5 Fig. 1. Schematic of the model system studied here. matical formulation. In Sect. 5, on the basis of the impact model proposed here, we shall apply some asymptotic techniques to examine the existence and stability of steady periodic system responses in which the V 0 reservoir never stays in contact with the actuator. Such responses are interesting because they have an actual duty cycle of zero, and are sustained solely by periodic collisions of the V 0 reservoir with the actuator. Then, through numerical solutions of the complete model, a systematic parameter study will be carried out. As it turns out, the system dynamics is rather rich, and strongly depends upon the driving frequency of the actuator and other geometrical and material properties of the system (see Sect. 6). In addition to clarifying the different types of actuator-reservoir interactions, the corresponding mean flowrate and reservoir pressures will also be calculated and discussed in Sect. 6. Finally, some remarks will be given in Sect. 7 to conclude this paper.

Basic Formulation

Model system setup and constitutive relations

Figure 1 depicts the closed-loop valveless pumping system studied here; it consists of two distensible reservoirs-the V 0 and V 1 reservoirs-connected by two rigid tubes. The system is completely filled with an incompressible liquid (say, water) of density ρ and viscosity µ. The instantaneous volumes of the two reservoirs are V 0 (t) and V 1 (t), respectively, where t is time. The rigid tube on the left-hand side of the system (hereafter referred to as the L-tube) has radius a L and length l L , while that on the right-hand side of the system (the R-tube) has radius a R and length l R . Since the tubes are assumed to be rigid, a L,R and l L,R have constant values.

Suppose also that, around the V 0 reservoir, there is a compression actuator contracting and expanding periodically. In particular, the physical presence of the actuator sets the maximum admissible instantaneous volume of the V 0 reservoir. We shall denote such a maximum admissible volume by V a (t). Also, as a specific example, it is assumed that V a (t) varies sinusoidally, with an angular frequency ω, between V min (> 0) and V max (> V min ):

V a (t) = 1 2 (V max + V min ) + 1 2 (V max -V min ) cos ωt. (1)
When V 0 (t) < V a (t), the V 0 reservoir is not in contact with the actuator, and so its volume would vary freely (but in a way that conserves mass and satisfies momentum balance; see Sect. 2.2). For convenience, in the sequel we shall say that the system is in the "free mode" when V 0 (t) < V a (t). On the other hand, whenever the V 0 reservoir is in contact with the actuator-and the system is said to be in the "contact mode," we shall simply set V 0 (t) = V a (t). In view of this requirement, the present model system of valveless pumping essentially is "volume-controlled" (as an anonymous reviewer of this paper kindly pointed it out), although V 0 (t) is not affected by the actuator motion during the free mode. Further details for various interactions between the actuator and the V 0 reservoir will be discussed in Sect. 4.

To construct a model where nonlinearities arise only from interactions between the V 0 reservoir and the actuator, the reservoir volume is assumed to vary linearly with the pressure difference across its elastic wall. Specifically, with p 0 (t) being the fluid pressure inside the V 0 reservoir, and p e (t) its external pressure, one has

V 0 (t) = V 0f + C 0 • {p 0 (t) -p e (t)}, (2)
where V 0f is the free (i.e., unstressed) volume of the V 0 reservoir, and C 0 is a constant compliance coefficient. Meanwhile, it is assumed that the actuator can only exercise a positive pressure on the V 0 reservoir, so that p e > 0 when the system is in the contact mode, and p e = 0 in the free mode. A similar constitutive relation for the V 1 reservoir having an internal fluid pressure p 1 (t) is

V 1 (t) = V 1f + C 1 • p 1 (t). (3)
Here, in the absence of an additional actuator, the V 1 reservoir has no external pressure. Like C 0 , the compliance coefficients C 1 is also taken to be constant.

Mass conservation and momentum balance

Denoting the volumetric flowrates (from the V 0 reservoir to the V 1 reservoir) inside the L-and R-tubes by Q L and Q R , respectively, mass conservation requires that

V0 = -Q L -Q R , (4)
V1 = Q L + Q R , (5)
where the overdots represent derivatives of the state variables with respect to time t. Recall that, when the system is in the contact mode, V 0 (t) = V a (t) is prescribed by Eq. (1); Eq. ( 4) then requires the flowrates Q L and Q R to vary consistently. On the other hand, in the free mode Eq. ( 4) should be used in conjunction with other governing equations to determine V 0 (t). Again, in order not to introduce additional nonlinearities, the resistance to the fluid flow in each tube is assumed to be linearly proportional to the volumetric flowrate therein. Likewise, the inertia of the fluid in each tube is taken to be constant. The momentum equations for the L-and R-tubes therefore are

I L QL + R L Q L = p 0 -p 1 , (6)
I R QR + R R Q R = p 0 -p 1 , (7)
where R L,R and I L,R are the constant resistance and inertia coefficients, respectively. It can be readily deduced from Eq. ( 6) that, if the L-tube stands alone (i.e., not connected to the system), the intrinsic time scale for the variation of

Q L is T L = I L /R L .
Similarly, the intrinsic time scale for Q R variation when the R-tube stands alone is

T R = I R /R R .
So far we have introduced 7 state variables (namely V 0,1 , p 0,1 , p e , and Q L,R ) and 11 input parameters (ω,

V max , V max , V 0f , V 1f , C 0,1 , I L,R , R L,R
) of the model system. As non-dimensionalization of the above mathematical formulation does not significantly reduce the number of variables and parameters (by three at most, since three independent physical dimensions-namely mass, time, and length-appear in the formulation), here we shall proceed with the dimensional formulation. However, in the sequel, the main analytical and numerical results will be put in suitable dimensionless form, so as to bring out some similarities of the system dynamics.

Initial conditions

Here we shall consider the special case that the V 0 and V 1 reservoirs are first filled to their free volumes before the compression actuator is installed. As can be readily deduced from Eqs. (4) and (5), mass conservation implies that the sum of the fluid volumes in the two distensible reservoirs remains invariant with time: V 0 (t)+V 1 (t) = V 0f +V 1f . Suppose then that an actuator is installed and compresses the V 0 reservoir to a certain extent before it starts moving. The initial volume of the V 0 reservoir, V max , therefore is smaller than its free volume V 0f , and we may define the "precompression ratio" as γ = (V 0f -V max )/V 0f (> 0). In terms of γ, the initial volumes of the two reservoirs are given by V 0 (0) = V 0f (1 -γ) and V 1 (0) = V 1f + γV 0f , respectively. Also, with ∆V a = (V max -V min )/2, Eq. (1) can be rewritten as

V a (t) = V 0f (1 -γ) -∆V a (1 -cos ωt). (8)
Meanwhile, since the fluid in the system initially is motionless, one has Q L (0) = Q R (0) = 0. It can also be readily calculated from the constitutive relations, Eqs. ( 2) and ( 3), that initially the uniform internal pressure of the system is p 0 (0) = p 1 (0) = γV 0f /C 1 , and the external pressure exerted on the V 0 reservoir by the stationary actuator is p e (0) = γV 0f /C 01 , with C 01 given by 1/C 01 = 1/C 0 + 1/C 1 . Note that the parameter C 01 may be viewed as a "characteristic compliance" of the system.

Sample parameter values

Before we proceed, let us envision a sample system having realistic dimensions and material properties. A dimensionless compliance parameter will then be identified to be relatively small under realistic conditions (see Sect. 3.2), and some asymptotic analysis of the system responses will be carried out accordingly. Specifically, suppose that the V 0 and V 1 reservoirs are made of latex rubber, so that a realistic estimate for the elastic modulus is E = 2 MPa. Let us also think of the two reservoirs as elastic balls, with free (unstressed) radii of a 0 = 25 mm and a 1 = 50 mm, respectively, and the same skin thickness of h = 5 mm. The free volumes of the two reservoirs therefore are V 0f = 65.4 ml and V 1f = 523.6 ml, respectively. Meanwhile, despite that relatively large volume changes will be allowed in our numerical calculations, the compliance coefficients of the two reservoirs are estimated by that of thin-walled spherical shells having mean radii a 0,1 and undergoing small volume changes, namely C 0,1 = 2πa 4 0,1 /Eh. This somewhat crude approximation yields C 0 = 2.45 × 10 -4 ml/Pa and C 1 = 3.93 × 10 -3 ml/Pa. (We have also used expressions of elasticity for a thick-walled hollow sphere to estimate the compliance coefficients, and the results are only about 20% smaller than that listed above. So, the values of C 0,1 specified above indeed are reasonable estimates of the compliance coefficients. And of course, more accurate estimates of the parameter values may be needed if one wishes to compare the model predictions with relevant experimental data in future works.)

Suppose also that the L-and R-tubes have the same length l L = l R = 200 mm, but differing inner radii, a L = 5 mm and a R = 7.5 mm. Water at room temperature, with density ρ = 1000 kg/m 3 and viscosity µ = 0.00112 Pa•sec, is taken to be the working fluid in the system. Using the Poiseulle's law, the flow resistance coefficients R L,R appearing in Eqs. ( 6) and ( 7) are related to the corresponding tube radii a L,R by R L,R = 8µl L,R /πa 4 L,R , and so R L = 0.913 Pa•sec/ml and R R = 0.180 Pa•sec/ml, respectively. Meanwhile, the inertia coefficients I L,R in Eqs. ( 6) and ( 7) are taken to be proportional to the mass of the fluid in the L-and R-tubes:

I L,R = ρl L,R /πa 2 L,R . It is then calculated that I L = 2.55 Pa•sec 2 /ml and I R = 1.13 Pa•sec 2 /ml.

Threshold Frequency for the Onset of Free Mode

Note that choosing a positive pre-compression ratio γ implies that the nominal duty cycle is unity, regardless of the driving frequency of the system. However, as we shall show here, when the driving frequency of the system exceeds a certain welldefined threshold value, the V 0 reservoir would separate from the actuator. When that happens, the actual duty cycle would be less than unity, and the mathematical Effects of Actuator Impact on the Nonlinear Dynamics of a Valveless Pumping System 9 model discussed above has to be supplemented by other components (see Sect. 4).

Steady periodic system response in the contact mode

First we show that, when the system is in the contact mode, the governing equations can be reduced to a single first-order differential equation for Q L . Briefly, using Eq. ( 4) one may write

Q R = -V0 -Q L = -Va -Q L (
the last equality holds in the contact mode). Meanwhile, eliminating the pressure difference p 0 -p 1 from Eqs. ( 6) and ( 7) yields

I L QL + R L Q L = I R QR + R R Q R = -I R • ( Va + QL ) -R R • ( Va + Q L ).
For the particular choice of V a (t) given by Eq. ( 8), this result gives

(I L + I R ) QL + (R L + R R ) Q L = -R R • Va -I R • Va = ∆V a • (ωR R sin ωt + ω 2 I R cos ωt). (9)
Given a suitable initial value of Q L , one can solve Eq. ( 9) analytically for Q L (t), and then calculate all the other state variables accordingly using Eqs. ( 2)- (7). Equation ( 9) also indicates that, when the L-and R-tubes are connected with the V 0 and V 1 reservoirs, the intrinsic time scale for the flowrate variation in the closed loop is

T LR = (I L + I R )/(R L + R R ), which generally is different from the time scales T L = I L /R L and T R = I R /R R for stand-alone tubes.
The sample parameter values specified in Sect. 2.4 give T L = 2.79 sec, T R = 6.28 sec, and T LR = 3.37 sec. Assume now that the system has a steady periodic response completely in the contact mode. Accordingly, one may write

Q L (t) ω∆V a = A cos ωt + B sin ωt, (10)
where the coefficients A and B can be readily determined by substituting Eq. ( 10) into Eq. ( 9) and then comparing coefficients of the sine and cosine functions in the resulting expression. This yields

A = R LR R L • ω(T R -T LR ) 1 + ω 2 T 2 LR , B = R LR R L • 1 + ω 2 T R T LR 1 + ω 2 T 2 LR ,
with R LR , given by 1/R LR = 1/R L + 1/R R , being the overall resistance coefficient when the L-and R-tubes are connected in parallel. It then transpires that, if the system is completely in the contact mode-and hence is perfectly linear, no pumping effects are produced on average. This also points to the fact that asymmetry alone is not sufficient for producing a mean flow in a valveless pumping system; certain nonlinearities of the system also are needed (i.e., a perfectly linear valveless pumping system excited sinusoidally at an asymmetric site still produces no mean flow).

Threshold frequency

Let us proceed with the analysis anyway. With the above exact result, Eq. ( 6) gives

p 0 (t) -p 1 (t) ωR L ∆V a = (A + BωT L ) cos ωt + (B -AωT L ) sin ωt. ( 11
)
When the system is in the contact mode,

V 1 (t) = V 0f + V 1f -V a (t)
, and so Eq. ( 3)

gives p 1 (t) = [V 0f -V a (t)]/C 1 .
Moreover, using Eqs. ( 2) and ( 8), the external pressure of the V 0 reservoir can be written as p e (t) = p 0 (t) -p 1 (t) + [γV 0f + ∆V a (1cos ωt)]/C 01 . Combining this result with Eq. ( 11) yields

p e (t) ∆V a /C 01 = S 0 -S 1 cos ωt -S 2 sin ωt, (12)
where S 0 = 1 + γ/δ and

S 1 = 1 - ǫ 2 ω 2 T 2 LR 1 + ω 2 T 2 LR T L + T R -T LR T LR + ω 2 T L T R , S 2 = - ǫ 2 ωT LR 1 + ω 2 T 2 LR 1 + ω 2 T LR (T L + T R ) -ω 2 T L T R .
In the above expressions, δ = ∆V a /V 0f is a dimensionless "actuator squeeze volume" and ǫ 2 = C 01 R LR /T LR is a dimensionless compliance parameter.

Using the sample parameter values specified in Sect. 2.4, the characteristic compliance and resistance parameters of the system are calculated to be C 01 = 2.31 × 10 -4 ml/Pa and R LR = 0.150 Pa•sec/ml. Recall also that T LR = 3.37 sec as calculated above, and so we have ǫ 2 = 1.03 × 10 -5 . The extremely small value of ǫ 2 implies that the reservoirs are relatively rigid. Of course, it is possible to specify other material properties and dimensions of the system to have a greater value of ǫ.

In the ensuing analysis, however, we shall exploit the smallness of the dimensionless compliance parameter ǫ to simplify some of the results.

Recall now that the external pressure of the V 0 reservoir, p e , has to be nonnegative; hence Eq. ( 12) is an admissible solution only if S 2 0 ≥ S 2 1 + S 2 2 . For weakly compliant systems with ǫ ≪ 1, if the driving frequency of the actuator ω = O(1)

then S 2 1 + S 2 2 = 1 + O(ǫ 2 ), while S 2 0 = (1 + γ/δ) 2 > 1.
The above admissibility condition therefore is satisfied. However, it can be deduced by use of dominant balance analysis [START_REF] Johnson | Singular Perturbation Theory[END_REF] that the "threshold frequency" ω = ω th -at which

S 2 0 = S 2 1 + S 2 2 - is of O(ǫ -1
). Accordingly, ωT LR , ωT L and ωT R are O(ǫ -1 ) quantities, and the expressions for S 1,2 derived above simplify to S 1 = 1 -(ǫω) 2 T L T R + O(ǫ 2 ) and S 2 = O(ǫ). One may then approximate the aforementioned admissibility condition by

S 2 0 ≥ [1 -(ǫω) 2 T L T R ] 2
to leading order, and deduce that the system has a steady periodic response completely in the contact mode only when ω < ω th , with

ω th ∼ ǫ -1 1 + S 0 T L T R 1/2 (13)
for ǫ → 0 + . As will be shown in Sect. 5, Eq. ( 13) gives an accurate estimate for ω th , provided that the dimensionless compliance parameter ǫ is sufficiently small.

Mode Switching and Jump Conditions

In Sect. 2, we have deliberately used linear constitutive relations and conservation laws to construct the basic mathematical model. Also, the system is assumed to

Effects of Actuator Impact on the Nonlinear Dynamics of a Valveless Pumping System 11 be pre-compressed so that the nominal duty cycle is unity. Consequently, no mean flow would be produced unless the driving frequency of the system ω exceeds a well defined threshold value ω th , and the V 0 separates from the actuator for some part of the compression cycle. To complete the mathematical formulation, we now discuss a number of mode switching and jump conditions for the interactions between the V 0 reservoir and the actuator. Nonlinearities resulting from such conditions are responsible for the mean flow generation in the system.

Separation of the V 0 reservoir from the actuator

As shown in Sect. 3.2, while the system is pre-compressed, as time goes on the external pressure of the V 0 reservoir, p e given by Eq. ( 12), may decrease to zero, say, at t = T s . The V 0 reservoir then separates from the actuator, and the system switches from the contact mode to the free mode. After t = T s , the volume of the V 0 reservoir varies freely (i.e., without forcing), and we set p e = 0. Accordingly, Eqs. ( 2) and ( 3)

give p 0 (t) = -[V 0f -V 0 (t)]/C 0 , p 1 (t) = [V 0f -V 0 (t)]/C 1 ,
and so in the free mode

p 0 (t) -p 1 (t) = -[V 0f -V 0 (t)]/C 01 < 0. ( 14
)
This implies that, when the system is in the free mode, the difference in the reservoir pressures tends to drive the working fluid from the V 1 reservoir to the V 0 reservoir. (Due to inertial effects, however, the flowrates Q L,R still may be positive.) Note also that all of the state variables vary continuously at t = T s , and the only discontinuity occurs in ṗe , which "jumps" from a finite negative value at t = T - s to zero at t = T + s .

Blocking of the fluid loop

Now, when the system is in the free mode, V 0 (t) may decrease to zero. (Note that this would not happen for V 1 (t), because then we would have

V 0 = V 0f + V 1f > V max ; but V max is the maximum value of V 0 (t) for ∀t > 0.
) When the V 0 reservoir is completely depleted, say, at t = T 0 , the fluid loop would be blocked, and so we set both Q L and Q R to zero, regardless of what values they may have at t = T - 0 . Other state variables, however, are continuous at t = T 0 . At t = T + 0 , because p 0 -p 1 = -V 0f /C 01 as indicated by Eq. ( 14), the working fluid would immediately start refiling the V 0 reservoir, thereby increasing its volume again.

Collision of the V 0 reservoir with the actuator

With or without blocking of the fluid loop happening first, during the free mode of the system, the volume of the V 0 reservoir would eventually become equal to V a (t), say, at t = T c . Of course, for such a "collision" to occur, it is necessary that V0 (T - c ) > Va (T c ). However, upon collision the physical existence of the actuator would render V0 (T + c ) = Va (T c ), or, using Eq. ( 4),

Q L (T + c ) + Q R (T + c ) = -Va (T c ). ( 15
)
Let us also assume that the collision at t = T c produces a delta-function-like variation of the external pressure p e , and a similar variation of p 0 consistent with the constitutive relation, Eq. ( 2). This means that, while the collision lasts for practically zero time, the "impulse" produced by the collision

I = T + c T - c p 0 (t) dt = 0.
Integrating then Eqs. ( 6) and ( 7) yields

I L • Q L (T + c ) -Q L (T - c ) = I = I R • Q R (T + c ) -Q R (T - c ) . (16)
Using Eqs. ( 15) and ( 16), we may express the flowrates after the collision by

Q L (T + c ) = I L • Q L (T - c ) -I R • [Q R (T - c ) + Va (T c )] I L + I R , (17)
Q R (T + c ) = I R • Q R (T - c ) -I L • [Q L (T - c ) + Va (T c )] I L + I R . (18)
Also, upon substitution of Eqs. ( 17) and ( 18) into Eq. ( 16), a simple expression for the impulse produced by the collision is obtained:

I = I LR • V0 (T - c ) -Va (T c ) (> 0), (19)
where the "characteristic inertia" I LR of the system is defined by 1/I LR = 1/I L + 1/I R . It then transpires that the magnitude of the impulse and the resulting flowrate jumps are directly proportional to the decrease in V0 ; see Eq. ( 16).

Mode switching at the reservoir-actuator collision

After the collision at t = T c , the V 0 reservoir may stay in contact with the actuator (and the system switches to the contact mode), or simply rebound (so that the system stays in the free mode). To decide what would happen, first we assume that the V 0 reservoir rebounds, so that p e (T + c ) = 0 and hence p 0 (T

+ c ) = -[V 0f - V a (T c )]/C 0 , p 1 (T + c ) = [V 0f -V a (T c )]/C 1 (see Sect. 4.1)
. Of course, one would also have Q L (T + c ) and Q R (T + c ) given by Eqs. ( 17) and (18), in addition to

V 0 (T + c ) = V a (T c ) and V 1 (T + c ) = V 0f + V 1f -V a (T c
). If these state conditions render V 0 < V a right after t = T c , it is interpreted that the V 0 reservoir just has to rebound. On the other hand, if V 0 > V a subsequently (which is inadmissible), the V 0 reservoir stays in contact with the actuator, and the system is switched to the contact mode.

Figure 2 summarizes the decision rules discussed above. Such rules can be readily integrated with the numerical procedures for solving Eqs. ( 2)- (7). Briefly, with the initial conditions specified in Sect. 2.3, the fourth-order Runge-Kutta method is used to integrate Eqs. ( 2)- (7). At each time step, the decision rules are applied to check if any mode switching occurs, and, if necessary, reset the values of the system's state variables. In general, the computation is terminated when a steady state of the system is reached. In our computations, when the relative variations Effects of Actuator Impact on the Nonlinear Dynamics of a Valveless Pumping System 13 Fig. 2. Summary of the mode decision rules for the model system studied here. of all the system state variables after a complete compression period are less than 10 -13 (or 10 -11 %), which in fact is a rather stringent criterion, it is decided that a steady state is established. Generally speaking, it takes a few dozen compression cycles for the system to reach steady periodic responses with complete depletion of the V 0 reservoir, while a few hundred compression cycles may be needed for steady periodic responses without complete depletion of the V 0 reservoir. Care is also taken to ensure that the time step is sufficiently small to adequately resolve the system dynamics. Typically, the number of time steps in a complete compression cycle is on the order of 10 3 -10 4 . The numerical results will be discussed in due course.

Steady Periodic Purely-Free System Response

As shown in Sect. 3.1, when the system is driven at a frequency ω lower than the threshold frequency ω th , it would have a steady periodic response completely in the contact mode. A related interesting question then is, would it be possible for the system to have a steady periodic response that is completely in the free mode? If so, the actuator would collide with the V 0 reservoir at a particular phase of the compression cycle, but the collision would send the V 0 reservoir right back to the free mode. The actual duty cycle of the system then is zero, while the nominal duty cycle quite misleadingly is unity. Here we examine the existence and stability of such responses-hereafter referred to as the steady periodic "purely-free" responses for convenience.

Asymptotic solution of the free-mode system response

To begin with, recall that when the system is in the free mode, p e (t) = 0 and

p 0 -p 1 = -(V 0f -V 0 )/C 01 = I L QL + R L Q L = I R QR + R R Q R . ( 20
)
Substituting Eq. ( 4) into the last equality in Eq. ( 20) yields (

I L + I R ) QL + (R L + R R )Q L = -I R V0 -R R V0 . Also, Eq. ( 20
) gives V 0 -V 0f = C 01 (I L QL + R L Q L ).
Combining these two results, it is found that in the free mode Q L (t) satisfies

ǫ 2 d 3 Q L dt 3 + T L + T R T L T R d 2 Q L dt 2 + 1 + ǫ 2 T L T R • dQ L dt + Q L T L T R T LR = 0. ( 21
)
It can be readily shown that Q R (t) and V 0 (t) -V 0f also are governed by Eq. ( 21).

Proceed with the assumption that the dimensionless compliance parameter ǫ 2 ≪ 1 (see Sect. 3.2). The presence of the small parameter ǫ 2 multiplying the highestderivative term d 3 Q L /dt 3 in Eq. ( 21) implies that Q L varies on a fast time scale. Specifically, a dominant-balance analysis reveals that the time scale is of O(1/ǫ). [START_REF] Johnson | Singular Perturbation Theory[END_REF] The same conclusion applies for Q R (t) and [V 0 (t) -V 0f ] as well. One may then introduce the fast temporal variable τ = t/ǫ and transform Eq. ( 21) into

Q ′′′ L + ǫ(T L + T R ) T L T R Q ′′ L + 1 + ǫ 2 T L T R Q ′ L + ǫQ L T L T R T LR = 0, ( 22
)
where the primes denote differentiations with respect to τ , i.e., (•) ′ = d(•)/dτ . The leading-order approximation of Eq. ( 22) for ǫ → 0 + is

Q ′′′ L + Q ′ L /T L T R = 0, (23)
which has the general solution

Q L (τ ) = C 1 + C 2 cos(τ / √ T L T R )+ C 3 sin(τ / √ T L T R
). However, Eq. ( 23) is a valid approximation of Eq. ( 22) only for fast varying Q L , and hence the constant solution must be discarded. Accordingly, for ǫ → 0 + , we have

Q L (t) ∼ C L ǫI L √ T L T R sin t ǫ √ T L T R + φ L (24)
to leading order, where the amplitude factor C L and phase constant φ L are to be determined. Note also that some extra factors are introduced in the coefficient of Eq. ( 24) for some convenience in subsequent algebra. Similarly, we have for

ǫ → 0 + Q R (t) ∼ C R ǫI R √ T L T R sin t ǫ √ T L T R + φ R (25)
and

V 0 (t) -V 0f ∼ C V I LR cos t ǫ √ T L T R + φ V , (26)
Effects of Actuator Impact on the Nonlinear Dynamics of a Valveless Pumping System 15 with undetermined (as yet) amplitude factors C R,V and phase constants φ R,V . Now, with Eqs. ( 14) and ( 24)-( 26), balancing the leading-order, namely O(1/ǫ 2 ), terms in Eqs. ( 6) and ( 7) requires that

C L,R ǫ 2 T L T R cos t ǫ √ T L T R + φ L,R = C V ǫ 2 T L T R cos t ǫ √ T L T R + φ V .
It is then deduced that

C L = C R = C V ≡ C, φ L = φ R = φ V ≡ φ.
Without any loss of generality, it is assumed that C > 0.

Periodicity and separation conditions

As we are seeking steady periodic system responses, the V 0 reservoir must have the same volume every time when it collides with the actuator. In view of Eq. ( 26), these observations imply that, in a steady periodic purely-free system response, V 0 (t) must first decrease and then increase. In addition, if at one collision, say at t = t 1 , the phase of the cosine function in Eq. ( 26)

t 1 ǫ √ T L T R + φ = π -α (27)
for a certain α ∈ (0, π/2), then at the next collision, t = t 2 , one must have

t 2 ǫ √ T L T R + φ = π + α. ( 28
)
Note also that between the two consecutive collisions, the motion of the compression actuator must complete an integral number of cycles, so that ω(t 2 -t 1 ) = 2nπ where n = 1, 2, 3, • • •. However, here we are concerned with the case that the actuator motion completes exactly one cycle between two consecutive collisions, and write ωt 1 = θ, ωt 2 = θ + 2π. Recall also that the driving frequency ω must exceed the threshold frequency ω th , which is of O(1/ǫ) as indicated by Eq. ( 13). Let us then take ω = O(1/ǫ) as well, and explicitly write ω = ω/ǫ, with ω = O(1) being a rescaled driving frequency parameter. Equations ( 27) and (28) then yield

α = π ω √ T L T R , φ = π - θ + π ω √ T L T R . (29)
Formally, for a specified driving frequency (i.e., ω), we have obtained an explicit expression for the phase parameter α appearing in Eqs. ( 27) and (28). However, the other phase parameter φ has only been related to the collision phase θ of the actuator motion, which remains to be determined. Meanwhile, the amplitude parameter C appearing in Eqs. ( 24)-(26) also is undetermined as yet.

To determine θ, φ, and C, note first that V 0 (t 1,2 ) = V a (t 1,2 ). With Eqs. ( 8) and (26), this implies that -C cos α/I LR = -γV 0f -∆V a (1 -cos θ); hence we have

cos θ = S 0 - C cos α I LR ∆V a . ( 30
) T.-S. Yang & C.-C. Wang
Next, using Eqs. ( 24) and ( 25), the jumps in the flowrates Q L and Q R upon a collision are calculated as follows:

Q L (t + 1 ) -Q L (t - 1 ) = Q L (t 1 ) -Q L (t 2 ) ∼ 2C sin α/ǫI L √ T L T R ; similarly, Q R (t + 1 ) -Q R (t - 1 ) ∼ 2C sin α/ǫI R √ T L T R .
These results are consistent with Eq. ( 16), which also gives the impulse of the collision at t = t 1 : I ∼ 2C sin α/ǫ √ T L T R (> 0). Meanwhile, Eq. ( 19) gives

I = I LR [ V0 (t - 1 )-Va (t 1 )] = I LR [ V0 (t 2 )-Va (t 1 )] ∼ ǫ -1 (C sin α/ √ T L T R + ωI LR ∆V a sin θ).
Equating the above two equivalent results then yields

sin θ = C sin α ωI LR ∆V a √ T L T R = Cα sin α πI LR ∆V a , (31)
where the last equality is obtained by use of Eq. ( 29). Moreover, as explained in Sect. 4.4, to ascertain that the collision at t = t 1 does not switch the system to the contact mode, V 0 (t) must decrease faster than V a (t) does at t = t + 1 . Since V0 = Va after collision (see Sect. 4.3), this implies that V0 (t + 1 ) -Va (t 1 ) < 0. Using Eqs. ( 26), ( 29) and (30), it is calculated that V0

(t + 1 ) -Va (t 1 ) ∼ ǫ -2 ω 2 ∆V a [S 0 -(1 - α 2 /π 2 )(C/I LR ∆V a ) cos α]; hence it is required that C I LR ∆V a ≥ S 0 (1 -α 2 /π 2 ) cos α . (32)
Now, making use of the basic trigonometric identity that sin 2 θ + cos 2 θ = 1, it is derived from Eqs. ( 30) and (31) that

α 2 π 2 + 1 - α 2 π 2 cos 2 α C I LR ∆V a 2 -2S 0 cos α C I LR ∆V a + S 2 0 -1 = 0.
One root of the above algebraic equation clearly violates the separation condition given by Eq. (32), and the only possibly admissible root is written as follows:

C I LR ∆V a = S 0 cos α + cos 2 α -(α 2 /π 2 )(S 2 0 -1) sin 2 α α 2 /π 2 + (1 -α 2 /π 2 ) cos 2 α . (33)
After some straightforward yet tedious algebra, it can be shown that, for the above root to satisfy Eq. (32), one must have

cos α ≥ S 0 • (α/π) (1 -α 2 /π 2 ) 1/2 [1 + (α 2 /π 2 )(S 2 0 -1)] 1/2 . ( 34
)
The right-hand side of Eq. ( 34) clearly vanishes at α = 0 and increases monotonically with α for α ∈ (0, π/2). Meanwhile, in the same interval, cos α decreases from unity to zero with increasing α. It then transpires that steady periodic purely-free responses exist only when the value of α is less than a critical value, denoted by α cr below, which can be determined numerically by equating the two sides of Eq. (34).

In view of the first part of Eq. ( 29), this means that such responses are possible only when the driving frequency of the system exceeds a corresponding critical value:

ω ≥ π/α cr ǫ √ T L T R ≡ ω cr . ( 35
)
Effects of Actuator Impact on the Nonlinear Dynamics of a Valveless Pumping System 17

Loop blocking condition

In the above calculations, it is implicitly assumed that V 0 (t) > 0 (so that the fluid loop is never blocked). Here let us examine the parametric condition under which such an assumption is valid. To do so, note first that the minimum of V 0 (t) is V 0f -C/I LR during the steady periodic purely-free system response; see Eq. ( 26). Denoting this minimum by γ m V 0f , one then has C/(I LR ∆V a ) = (1-γ m )V 0f /∆V a = (1 -γ m )/δ, which can be combined with Eq. (33) to yield the following expression:

δ = [(1 -γ m ) cos α -γ] 2 + (1 -γ m ) 2 (α/π) 2 sin 2 α 2[(1 -γ m ) cos α -γ] . (36)
For given values of γ and γ m , Eq. ( 36) can be used to plot δ against α (or ω).

As a specific example, let us use the parameter values specified in Sect. 2.4, and take the maximum actuator volume V max = 0.6V 0f (so that the pre-compression ratio γ = 0.4). The threshold frequency ω th predicted by the asymptotic result, Eq. ( 13), then is plotted as the solid curve in Fig. 3; the corresponding numerical results (obtained by use of the numerical procedure summarized in Fig. 2) also are as the overlapping symbols. For convenience, both the dimensional values (for the present parameter setting) and normalized values (by the characteristic frequency ǫ -1 (T L T R ) -1/2 ) of ω th are shown in Fig. 3. Also, while Eq. ( 13) predicts that the leading-order approximation of the normalized threshold frequency, ǫω th √ T L T R , depends on the ratio of γ/δ only, here the more straightforward (from a practical standpoint) volume ratio V min /V max = (1 -γ -2δ)/(1 -γ) is used as the abscissa in Fig. 3. It is clearly seen that, given the extremely small value of ǫ 2 = 1.03 × 10 -5 here, the numerical and asymptotic results are in excellent agreement.

Moreover, the critical frequency beyond which the purely-free system responses arise, i.e., ω cr defined by Eq. ( 35), and the parameter boundary for such responses to have a non-negative V 0 (t), specified by Eq. ( 36) with γ m = 0, are shown as the short-dashed and dot-dashed curves in Fig. 3, respectively. It is important to note that the ω cr curve intersects the γ m = 0 curve at about V min /V max = 0.30 (with ǫω th √ T L T R = 3.7) and terminates there, since Eq. ( 35) is valid only when V 0 (t) is non-negative. Again, the asymptotic and numerical results of ω cr agree extremely well. Other curves in Fig. 3 will be explained later.

Linear stability analysis

Since the steady periodic purely-free system responses (having an actual duty cycle of zero) are sustained solely by collision of the actuator with the V 0 reservoir, they may become unstable. To analyze the stability of such responses, suppose that one collision occurs at a slightly deviated phase (rather than the correct phase θ determined in Sect. 5.2), so that

ωt 1 = θ + δ a1 , t 1 ǫ √ T L T R = π -α + δ v1 , T.-S. Yang & C.-C. Wang
ωt 2 = θ + 2π + δ a2 , t 2 ǫ √ T L T R = π + α + δ v2 , with |δ a2 |, |δ v2 | ≪ 1.
The steady periodic purely-free responses would be unstable if the system parameters render |δ a2 /δ a1 | > 1.

From the above expressions one derives that ω(t

2 -t 1 ) = 2π + δ a2 -δ a1 = ω √ T L T R (2α + δ v2 -δ v1
), and so

δ v2 = δ v1 + δ a2 -δ a1 ω √ T L T R .
Note also that, right after the first collision, V 0 (t + 1 ) = V a (t 1 ) and V0 (t + 1 ) = Va (t 1 ). Using Eqs. ( 8) and ( 26), with a small error δ c added to the amplitude factor C in the latter (|δ c | ≪ C), these conditions are written as

-γV 0f -∆V a [1 -cos(θ + δ a1 )] = C + δ c I LR cos (π -α + δ v1 ) , ( 37
) ω T L T R ∆V a sin(θ + δ a1 ) = C + δ c I LR sin (π -α + δ v1 ) . ( 38
)
Effects of Actuator Impact on the Nonlinear Dynamics of a Valveless Pumping System 19

Equations ( 37) and (38) then are linearized with respect to the correct response (having δ a1 , δ v1 , δ c = 0), and the resulting linear equations are solved to yield

δ c /δ a1 I LR ∆V a = sin θ cos α + π α cos θ sin α, ( 39
)
Cδ v1 /δ a1 I LR ∆V a = sin θ sin α - π α cos θ cos α. ( 40
)
Since only small collision phase deviations are considered here, there is no need to impose an additional condition for the separation of the V 0 reservoir from the compression actuator. Similarly, at the next collision, a volume matching condition analogous to Eq. ( 37) can be written down and then linearized. Using also Eqs. ( 39) and ( 40), the result is reduced to

δ a2 δ a1 sin θ + Cα/π I LR ∆V a sin α = sin θ cos 2α + π α cos θ sin 2α + Cα/π I LR ∆V a sin α.
Substituting Eq. (33) into the above result, one then deduces (with some algebra) that the stability condition

|δ a2 /δ a1 | ≤ 1 implies that 1 + π 2 S 2 0 α 2 1 - α 2 π 2 1 - α 2 π 2 cos 2 α - α 2 π 2 ≤ (1 -α 2 /π 2 ) cos 2 α + α 2 /π 2 (1 -α 2 /π 2 ) cos 2 α -α 2 /π 2 .( 41
)
With the separation condition, Eq. ( 34), it can be shown that the left-hand side of Eq. ( 41) is positive and decreases monotonically with α (from an infinitely large value at α = 0), while the right-hand side increases monotonically with α (from unity at α = 0). So, for the steady periodic purely-free response to be stable, the phase parameter α must exceed a certain critical value, denoted below by α st , which can be determined numerically from the equality in Eq. ( 41). In terms of the driving frequency ω of the system, this implies that such responses are stable when

ω ≤ π/α st ǫ √ T L T R ≡ ω st . (42)
Again, with the same parameter values specified in Sect. 2.4 and the precompression ratio γ = 0.4, the stability-margin frequency ω st calculated from Eq. ( 41) is shown as the long-dashed curve in Fig. 3. The ω st curve intersects the γ m = 0 curve at about V min /V max = 0.31 (with ǫω st √ T L T R = 4.2) and terminates there, since in deriving Eq. (41) a non-negative V 0 (t) is assumed. As shown in Fig. 3, when the volume ratio V min /V max exceeds about 0.34 (at which ǫω st √ T L T R = 4.2 and γ m ≈ 0.02), the numerical results (indicated by the overlapping symbols) agree well with the asymptotic results of ω st . However, when the V min /V max ratio is decreased below 0.34 (but still greater than the limiting value 0.31), a sudden increase in ω st is observed in the numerical results. The nature and implications of this will be explained in the next section through a detailed discussion of numerical results.

Further Numerical Results and Discussion

6.1. System dynamics for various volume ratios V min /V max The above asymptotic calculations have revealed the following features of the system dynamics: For a given pre-compression ratio γ and a sufficiently large volume ratio V min /V max (specifically, for γ = 0.4 and V min /V max > 0.34 here), when the driving frequency ω of the system is below the threshold frequency ω th given by Eq. ( 13), the V 0 reservoir would never separate from the compression actuator. In such cases, the system does not produce any net pumping effects. When ω exceeds ω th , but is lower than the critical frequency ω cr defined by Eq. ( 35), in the steady state the V 0 reservoir would separate from the actuator during a certain part of the compression cycle, thereby inducing some nonlinear effects that generally produce a non-zero mean flowrate in the fluid loop. Then, for ω > ω cr , the system response becomes purely-free, in which case the V 0 reservoir rebounds from the compression actuator when they collide, so that the actual duty cycle is zero. However, the purely-free system responses will lose their stability when ω exceeds a well-defined stabilitymargin frequency ω st given by Eq. ( 42).

As a specific example, let us take γ = 0.4 and V min /V max = 0.6 (corresponding to the thin vertical line marked with the Roman numeral "I" in Fig. 3). For such a parameter combination, the asymptotic results give ǫω th √ T L T R = 2.31 (or

f th = ω th /2π = 27.4 Hz), ǫω cr √ T L T R = 4
.28 (f cr = ω th /2π = 50.7 Hz), and ǫω st √ T L T R = 4.73 (f st = ω st /2π = 55.9 Hz). Also, the temporal variations of V 0 (t) are calculated for a number of driving frequencies, and the results are shown in Fig. 4. It is clearly seen that when the driving frequency f = ω/2π = 35 Hz, which is between f th and f cr , in the steady state the V 0 reservoir is in contact with the actuator for some part of the compression cycle (and so V 0 (t) = V a (t)), but moves freely for the rest of the cycle. When the system is driven at 52 Hz, which is between f cr and f st , the response indeed becomes purely-free. Note also that if the asymptotic results were included in Fig. 4(a), they would be almost indistinguishable from the numerical results. Now, for f = 56 Hz, which slightly exceeds the stability-margin frequency f st , it is clearly demonstrated in Fig. 4(b) that the precise collision phase varies from cycle to cycle, even after the initial transient has essentially died out (which is discarded in Fig. 4(b)). This, of course, is precisely a manifest of the instability explained in Sect. 5.4. It is also interesting to note that the response shown in Fig. 4(b) appears to repeat itself every two compression cycles. However, here we shall not attempt to trace the detailed transition from the steady periodic responses repeating themselves in one compression cycle to other more complex responses. (Yet some ongoing and future works in this direction will be remarked upon in Sect. 7.)

For smaller values of V min /V max , the results presented in Fig. 3 suggest that somewhat more complicated system responses may arise. While asymptotic analysis can be carried out (in principle at least) to clarify the system dynamics in such cases, (a) (b) Fig. 4. The temporal variations of V 0 (t) for γ = 0.4 and V min /Vmax = 0.6 when the system is driven at frequencies (a) 35 Hz, 52 Hz, and (b) 56 Hz, respectively. For comparison, the volume constraint Va(t) set by the actuator motion is also shown.

the algebra involved would be quite formidable. So, here we shall discuss the system dynamics by use of numerical results only. First, as shown in Fig. 5, for γ = 0.4 and V min /V max < 0.34, there exists a critical driving frequency ω v0 beyond which the V 0 reservoir would be completely depleted at a certain phase of the compression cycle (and start being refilled afterwards). It also appears that the ω v0 curve in Fig. 5 passes through the intersection of the ω cr and γ m = 0 curves (at V min /V max = 0.31 and ǫω st √ T L T R = 4.2), and terminates when it reaches the ω st curve. This has an important implication. Specifically, for 0.31 < V min /V max < 0.34 (with γ = 0.4), the asymptotic results still are applicable, and predict the existence of purely-free system responses with positive definite V 0 (t) for ω ∈ (ω cr , ω st ). The additional T.-S. Yang & C.-C. Wang Fig. 5. For γ = 0.4 and V min /Vmax smaller than about 0.34, when the driving frequency exceeds a critical value ω v0 -plotted here for various V min /Vmax ratios-the V 0 reservoir would be completely depleted at a certain phase of the compression cycle. Some results in Fig. 3 are duplicated here for reference.

presence of the ω v0 curve in this parameter region then implies the coexistence of system responses in which V 0 (t) ↓ 0 at a certain phase of the compression cycle.

As a demonstration, let us take V min /V max = 0.33 (corresponding to line II in Figs. 3 and5). For this case, the asymptotic results give ǫω th √ T L T R = 2.00 (or f th = 23.7 Hz), ǫω cr √ T L T R = 3.72 (f cr = 44.1 Hz), and ǫω st √ T L T R = 4.18 (f st = 49.4 Hz). Meanwhile, it is determined numerically that ǫω v0 √ T L T R = 3.97 (f v0 = ω v0 /2π = 47.0 Hz). When the system is driven at the frequency f = 48 Hz, which is higher than both f cr and f v0 and lower than f st , two different system responses indeed are found in the steady state, and are both shown in Fig. 6(a). In one of them, the V 0 reservoir has a positive volume throughout the compression cycle, and the numerically calculated response is visually indistinguishable from the asymptotic solution given by Eq. ( 26) (not shown in Fig. 6(a)). In the other response, the V 0 reservoir is completely depleted at a particular compression phase. Since in such a response V0 (t) has a non-zero (negative) value when V 0 (t) vanishes, it will be referred to as the "hard-landing" response in the sequel. Note also that both of the coexisting responses shown in Fig. 6(a) are purely-free, and are sustained by a collision of the V 0 reservoir with the actuator in each compression cycle. However, to find the purely-free steady periodic system responses having a positive definite V 0 (t) (or, the "never-landing" responses), the initial conditions of the system need to be carefully tuned. When the arbitrary initial conditions specified in Sect. 2.3 are used, we generally find the hard-landing purely-free responses instead. This seems to suggest that the hard-landing responses are more stable than the never-landing ones, in the sense that certain accidental depletions of the V 0 reservoir accompanying the initial transient would cause the system to respond in the hard-landing way.

Note also that, for driving frequencies f < f v0 (but exceeding f cr ), only the purely-free never-landing responses are found. Meanwhile, for values of f that are slightly higher than the stability-margin frequency f st predicted by the asymptotic theory, only the hard-landing purely-free responses are found. However, as indicated in Fig. 5, if f is further increased, the hard-landing purely-free responses also will become unstable. In particular, for V min /V max = 0.33, this higher stability-margin frequency is about 65 Hz (see Fig. 5). Taking then f = 70 Hz as an example, the resulting system response is shown in Fig. 6(b). It is interesting to note that the response appears to repeat itself every six compression cycles. However, as mentioned above, tracing the detailed transition of the simple periodic responses to other more complex responses is beyond the scope of this paper. Now, for even smaller volume ratios, i.e., V min /V max < 0.31, Fig. 5 indicates that the never-landing purely-free system responses would cease to exist. In such cases, the only useful asymptotic result is the still accurate estimate of the threshold frequency ω th for the onset of free-mode system responses. Here we shall take V min /V max = 0.15 as a particular example, corresponding to line III in Figs. 3 and5. For this parameter combination, it is determined numerically that f v0 is 33.1 Hz (or ǫω v0 √ T L T R = 2.80) and f st = 128 Hz (ǫω v0 √ T L T R = 10.8); see Fig. 3. Also, the asymptotic theory predicts that f th = 22.4 Hz (ǫω th √ T L T R = 1.89). Moreover, it is found that, for f ∈ (f th , f v0 ) (say, f = 28 Hz; see Fig. 7(a) for the resulting system response), the V 0 reservoir would separate from the compression actuator for some part of the compression cycle, and its volume during that part would be positive. Meanwhile, when f slightly exceeds f v0 , the V 0 reservoir still would be in contact with the compression actuator for some part of the compression cycle, but would be completely depleted during the free mode of the system response; see Fig. 7(a) for an example with f = 50 Hz. When the driving frequency is further increased (beyond a critical frequency f cr of about 71 Hz), however, the system response becomes purely-free, and the V 0 reservoir would be completely depleted at a certain phase of the compression cycle. An example with f = 100 Hz is also shown in Fig. 7(a). Finally, when f > f st (say, f = 130 Hz; see Fig. 7(b)), steady periodic responses repeating themselves in one compression cycle no longer can be found. While the response shown in Fig. 7(b) also appears to repeat itself every six compression cycles, more complex system responses would result from rather small variations of the driving frequency, which will not be discussed here.

Frequency dependence of actual duty cycle and mean flowrate

We now make some quantitative comparisons of the steady periodic system responses discussed above. First, the actual duty cycles for cases I, II, and III are plotted in Fig. 8. Recall that the actual duty cycle is the fraction of the compression cycle in which the V 0 reservoir is in contact with the actuator. In all our computations, it is found that the actual duty cycle becomes zero once the driving frequency ω exceeds the critical frequency ω cr for the onset of purely-free responses. And of course, the system responses would have a unity duty cycle if ω is lower than the threshold frequency ω th . So, only the results for ω ∈ (ω th , ω cr ) are shown in Fig. 8. It is interesting to note that, when the actual duty cycle is plotted against the nor-Effects of Actuator Impact on the Nonlinear Dynamics of a Valveless Pumping System 27 malized (dimensionless) driving frequency parameter ω = (ω -ω th )/(ω cr -ω th ), the curves for cases I and II (having V min /V max ratios of 0.6 and 0.33, respectively) appear to be quite similar, while that for case III (with V min /V max = 0.15) looks rather different. In particular, the curve for case III has a kink at the critical frequency ω v0 (beyond which the V 0 reservoir would be completely depleted at a certain phase of the compression cycle). Intuitively, when the V 0 reservoir is completely depleted, the system is restarted (by resetting both Q L and Q R to zero), so that the V 0 reservoir would collide with the actuator earlier than it would when negative values of V 0 (t) were allowed. This results in the small increase of the actual duty cycle when ω slightly exceeds ω v0 (followed by further decrease of the duty cycle for even higher driving frequencies). Meanwhile, the curves in Fig. 8 suggest that the frequency dependence of the actual duty cycle for different V min /V max ratios is somewhat similarity-with ω being the similarity parameter-as long as the V 0 reservoir has a positive volume throughout the compression cycle.

Next we discuss the frequency dependence of the mean volumetric flowrate in the L-tube, which is calculated by averaging Q L (t) over one compression period:

Q L = T 0 Q L (t) dt/T , with T = 2π/ω. (By mass conservation, Q R = -Q L in
the steady state; there therefore is no need to discuss Q R separately.) After being normalized by the characteristic scale ∆V a /T LR , the calculated results of Q L for cases I, II, and III are plotted as functions of the normalized driving-frequency parameter ω * = (ω -ω th )/(ω st -ω th ) in Fig. 9(a). Unlike in Fig. 8, here we present the results for driving frequencies ω up to the stability-margin frequencies ω st in the three case. Also, since the mean flowrate Q L may change sign as ω varies, shown in Fig. 9(a) is the absolute value of Q L .

Let us discuss the frequency dependence of Q L in case I (with V min /V max = 0.6) first; see the solid curve Fig. 9(a). Note that in case I the mean flowrate Q L is positive (i.e., there is a mean flow in the L-tube from the V 0 reservoir to the V 1 reservoir) when the driving frequency ω slightly exceeds the threshold frequency ω th . As ω further increases, Q L reaches a maximum near ω * = 0.10, and changes its sign near ω * = 0.62. (The dependence of the mean-flow direction on the driving frequency is a common feature shared by many valveless pumping systems. For example, in Ref. 15 such a dependence was demonstrated both experimentally and numerically.) The reversed mean flow reaches its maximum magnitude at about ω * = 0.73, and vanishes at about ω * = 0.81, which corresponds to the onset frequency ω cr for the purely-free responses in case I (marked with the abbreviation "ω I = ω cr " in Fig. 9(a)). When the driving frequency ω exceeds ω cr (but still is within the stability margin), the numerically calculated values of Q L become exceedingly small and cannot be shown in Fig. 9(a). While this result may seem surprising at first, it is actually predicted by the asymptotic calculations. Specifically, according to the leading-order results, Eqs. ( 24) and ( 25), the purely-free system responses would produce zero average flowrate in the fluid loop, which is now confirmed numerically. In case II (with V min /V max = 0.33), the frequency dependence of Q L depicted by the (blue) dashed curve in Fig. 9(a) is qualitatively similar to that in case I for driving frequencies up to ω cr (marked with "ω II = ω cr "). The only quantitative differences are that Q L now changes sign (from positive to negative) at about ω * = 0.38, and vanishes at about ω * = 0.50. However, while the purely-free responses having positive definite V 0 (t) again have a vanishingly small mean flowrate, a relatively large Q L is produced when the V 0 reservoir is completely depleted. Specifically, this occurs when the driving frequency ω exceeds the critical frequency ω v0 , which corresponds to ω * = 0.56 and is marked with "ω II = ω v0 " in Fig. 9(a). Note in particular that the resulting mean flowrate is about 100 times greater than the maximum mean flowrate that would be produced by system responses having a nonzero actual duty cycle. This makes good intuitive sense, because the system restarting following the complete depletion of the V 0 reservoir generally would produce a much stronger nonlinearity than that produced by the separation of the V 0 reservoir and the actuator. (It was also pointed out by Kenner et al. [START_REF] Kenner | The Liebau-effect or on the optimal use of energy for the circulation of blood[END_REF] that complete compression of a fluid-filled valveless tube, like the complete depletion of the V 0 reservoir here, is very effective in generating a mean flow in the tube. In the present model system, however, the V 0 reservoir has to separate from the actuator first, before it can possibly be completely depleted. Careful modeling of the interactions between the V 0 reservoir and the compression actuator therefore still is important.) For the same reason, the mean flowrate Q L in case III, with V min /V max = 0.15 and depicted by the (red) short-dashed curve in Fig. 9(a), also experiences a sudden increase of about 500 times as the driving frequency crosses ω v0 , corresponding to ω * = 0.10 and marked with "ω III = ω v0 ." Note also that since the depletion of the V 0 reservoir occurs first (i.e., at a lower critical driving frequency) before purely-free responses arise in case III, the calculated (positive) values of Q L are never vanishingly small. Some interesting observations also can be made by examining the temporal variations of the volumetric flowrates Q L,R in various types of system responses (as an anonymous reviewer kindly pointed it out). In particular, Q L,R (t) in three representative responses with V min /V max = 0.33 (i.e., case II discussed above) and driving frequencies f = 40, 45, and 50 Hz (corresponding to ω * = (ω -ω th )/(ω stω th ) = 0.39, 0.52, and 0.64), respectively, are are plotted in Fig. 10. It is clearly seen that, while Q L and Q R vary with time in similar manners, the instantaneous flowrate in the R-tube generally has a larger magnitude than that in the L-tube, because both the inertia and resistance coefficients of the R-tube are smaller than that of the L-tube; see Sect. 2.4. Note also that the sudden increases in Q L,R observed in Fig. 10 (around ωt/2π = 0.25) are caused by a collision of the V 0 reservoir with the actuator. In the case of f = 40 Hz, the impact of the collision relatively is smaller, and so the V 0 reservoir stays in contact with the actuator after the collision. As a result, Q L,R both continue to increase for a while after the collision. By contrast, in the cases of f = 45 and 50 Hz, the impact of the collision is somewhat stronger, and the V 0 reservoir separates from the actuator right afterwards. Also, the subsequent refilling of the V 0 reservoir causes Q L,R to decrease immediately after the collision. There is a qualitative difference between the cases of f = 45 and 50 Hz, however. Specifically, for the case of f = 50 Hz, the V 0 reservoir is completely depleted at ωt/2π ≈ 0.69, causing both Q L and Q R to be reset to zero. As already discussed above, this reset produces a rather substantial mean flowrate in the loop. It is also worth noting that the variation amplitudes of the normalized flowrates, Q L,R T LR /∆V a , are of O( 102 ) and so are much higher Effects of Actuator Impact on the Nonlinear Dynamics of a Valveless Pumping System 31 that accounts for mass and momentum balance in a closed fluid loop. As demonstrated in Sect. 3, when the system is pre-compressed and the driving frequency ω of the system is below a well defined threshold frequency ω th , the pliant part of the model system (i.e., the V 0 reservoir in Fig. 1) would never separate from the compression actuator, and the perfectly linear governing equations discussed in Sect. 2 would produce zero mean flowrate in the fluid loop. In the present model, a nonzero mean flowrate therefore can only be produced by the somewhat complex interactions of the compression actuator with the V 0 reservoir arising when ω > ω th .

To account for such interactions, the mathematical model is supplemented with a number of decision rules (discussed in Sect. 4) that are derived from fundamental mechanical considerations. The decision rules provide some nonlinearities for the overall mathematical model, while an additional nonlinearity comes with the restarting of the system when the V 0 reservoir is completed depleted at a certain phase of the compression cycle under suitable conditions (also discussed Sect. 4).

As the overall mathematical model involves rather intricate decision rules (summarized in Fig. 2), its solution generally has to be obtained numerically. However, by using a set of sample parameter values specified in Sect. 2.4, a dimensionless compliance parameter is identified to be relatively small under realistic conditions. We have then taken advantage of this small parameter, and carried out a series of asymptotic calculations to obtain accurate estimates for some critical driving frequencies across which the system response undergoes qualitative changes. In particular, the threshold frequency ω th for the separation of the V 0 reservoir from the compression actuator, the onset frequency ω cr and stability-margin frequency ω st of the purely-free responses (which have an actual duty cycle of zero) are obtained in Sect. 5. Furthermore, numerical solutions of the model for three representative cases are carefully discussed in Sect. 6 to elucidate the complex dynamics of the model system studied here, and to clarify the frequency dependence of the mean flowrate and mean pressures in the system.

While we have also studied the model system using other pre-compression ratios (than the particular choice of γ = 0.4 here), the dependence of the numerical results upon the driving frequency and other system parameters appear to be qualitatively similar to that discussed in this article. Such results therefore are not included here for brevity of the presentation. However, it should be noted that when the system is not pre-compressed, i.e., when the compression actuator initially is not in contact with the V 0 reservoir, they would collide with each other at any driving frequency, and are expected to produce a significant mean flowrate thereby. Although some numerical results have been obtained, a detailed discussion of them is beyond the scope of this article, and will be presented elsewhere.

In fact, even if we limit our scope to pre-compressed systems only, the system dynamics still is extremely complicated. For example, as suggested by the responses depicted in Figs. 6(b) and 7(b), there may exist steady periodic system responses having a period that is a multiple of the compression cycle period for higher driving frequencies. In fact, chaotic system responses may also arise under suitable conditions. While we have clarified some transition routes of such more complicated responses, including them here would render this article overly lengthy and poorly focused. Such transition routes therefore will be discussed in a follow-up work. Furthermore, even if we focus only on steady periodic system responses repeating themselves in one compression cycle, it will also be interesting to locate the precise boundaries between system responses of different types and completely map the parameter space shown in Fig. 3. These, again, will be presented in another work.

In closing, we would like to point out while the mathematical model proposed here may seem to be simple, it is capable of providing valuable insights into the actuator impact effects that have often been overlooked in previous studies of valveless pumping systems. Moreover, the results discussed in previous sections suggest that such effects very likely would be important in other more complicated valveless pumping systems as well. The actuator impact effects therefore deserve to be carefully considered in other modeling approaches (such as the computational and wave-pulse models mentioned in Sect. 1). Also, from a dynamical systems viewpoint, the system dynamics of such a seemingly simple model system is rich and complicated. Although some interesting results have been obtained and discussed in this article, there clearly still are more to be explored. And of course, it would be interesting as well to verify such results experimentally. Finally, from the viewpoint of biomedical applications, we expect that a deep and complete understanding of valveless pumping mechanisms would motivate certain innovative designs of blood pumps that can be used in cardiac assist devices. (An interesting axial blood pump, which does not rely upon valveless pumping though, is studied in Ref. 25.)

Fig. 3 .

 3 Fig.3. The threshold frequency ω th for the onset of free-mode responses, and the critical frequencies ωcr and ωst between which steady periodic purely-free responses stably exist. The curves are the asymptotic results, and the overlapping symbols are the corresponding numerical results. (The four numerical data of ωst for V min /Vmax < 0.34 also are connected by a curve for clarity.)

Fig. 6 .

 6 Fig. 6. The temporal variations of V 0 (t) for γ = 0.4 and V min /Vmax = 0.33 when the system is driven at frequencies (a) 48 Hz and (b) 70 Hz, respectively.

Fig. 7 .

 7 Fig. 7. The temporal variations of V 0 (t) for γ = 0.4 and V min /Vmax = 0.15 when the system is driven at frequencies (a) 28 Hz, 50 Hz, 100 Hz and (b) 130 Hz, respectively.

25 Fig. 8 .

 258 Fig. 8. The frequency dependence of the actual duty cycle for γ = 0.4 and V min /Vmax = 0.15, 0.33, and 0.60.

Fig. 9 .

 9 Fig. 9. The frequency dependence of various temporally averaged system characteristics for γ = 0.4 and V min /Vmax = 0.15, 0.33, and 0.60: (a) the mean volumetric flowrate in the L-tube (normalized by ∆Va/T LR ), and (b) the mean pressures p 0,1 in the V 0,1 reservoirs (compared with the contact-mode mean pressures p 0,1 0 , and normalized by ∆Va/C 0 and ∆Va/C 1 , respectively).

  T.-S. Yang & C.-C. Wang

Fig. 10 .

 10 Fig. 10. Temporal variations of the volumetric flowrates (a) Q L and (b) Q R for driving frequencies f = 40, 45, and 50 Hz. Here V min /Vmax = 0.33, and so the three cases correspond to case II in Fig. 9(a), with ω * = (ω -ω th )/(ωst -ω th ) = 0.39, 0.52, and 0.64, respectively.

Acknowledgments

The authors gratefully acknowledge the Taiwan National Science Council for supporting this work through grant No. NSC96-2628-E-006-241-MY3. They would also like to thank Professors M.-S. Ju, S.-Y. Lee, K.-S. Chen, and C.-Y. Wen of NCKU for a number of fruitful discussions on this work and other related topics.

than their mean values (see Fig. 9a). This is another common feature shared by many valveless pumping systems; see Ref. 15 for example.

We have also calculated the mean pressures in the V 0 and V 1 reservoirs for the three cases discussed above, and the results are presented in Fig. 9(b) for the same driving frequency ranges shown in Fig. 9(a). Specifically, the mean pressures are determined by evaluating p 0,1 = T 0 p 0,1 (t) dt/T , with T = 2π/ω. Also, to bring out some similarity of the results, the mean pressures in the two reservoirs when the system is in the contact mode (i.e., having an actual duty cycle of unity) are subtracted from p 0 and p 1 , and the resulting mean pressure differences are normalized by the characteristic pressure scales ∆V a /C 0 and ∆V a /C 1 , respectively. Using the exact results presented in Sect. 3, it can be readily calculated that the contact-mode mean pressures, denoted by p 0 0 and p 1 0 in Fig. 9(b), are equal and independent of the driving frequency:

Note that, when the driving frequency ω exceeds the threshold frequency ω th , the V 0 reservoir would separate from the actuator for some portion of the compression cycle, so that the temporal mean of its volume would be smaller than that in the contact mode. Also, by mass conservation, the mean volume of the V 1 reservoir would be greater than that in the contact mode. Accordingly, the V 1 reservoir is expanded to a greater extent in the free mode than in the contact mode, so that the mean pressure difference p 1 -p 1 0 is positive irrespective of the driving frequency, as shown in Fig. 9(b). Meanwhile, for ω > ω th , since the V 0 reservoir has a smaller mean volume and is in less complete contact with the actuator, its mean pressure is lower than that in the contact mode; see also Fig. 9(b).

Unlike on the mean flowrates, however, the driving frequency appears to have a milder effect upon the mean reservoir pressures. This, of course, results from the fact that the volumes of the distensible reservoirs in our model (and in reality as well) have a well defined lower limit of zero, and less apparent upper limits set by the magnitude of the pressure variations in the fluid loop. On the other hand, the mean flowrate in the loop is produced by various nonlinear effects having substantially differing strengths. Nevertheless, when the driving frequency ω crosses a particular critical frequency associated with a qualitative change of the system responses, the mean pressure variations also are affected to some extent. In particular, in cases II and III (see Fig. 9(b)), the effect is more apparent when ω crosses ω v0 (and the V 0 reservoir starts to be completely depleted at a certain phase of the compression cycle) than when it crosses ω st (and the system response becomes purely-free). Consistently, since the V 0 reservoir would never be completely depleted in case I, the mean pressures vary rather smoothly with ω.

Summary and Concluding Remarks

In order to focus upon the often-overlooked actuator impact effects on the dynamics and performance characteristics-such as the mean flowrate and mean pressure-of valveless pumping systems, here we have constructed a lumped-parameter model