SLOW VISCOUS PLOW IN A SYRINGE LAYNE T. WATSON~, S.C. BILLUPS1 C.Y. WANG, Member ASME2 E.A. EVERETT3 Technical Report 85-31 October 1, 1985 Abstract. The slow viscous flow in a syringe is modelled by the quasi-steady axisymmetric Stokes equation with a point sink for the needle hole. The governing equations are approximated using nonstandard finite difference formulas optimized for the boundary conditions, and solved numerically using a SOR technique. Streamlines and pressure profiles are computed for a variety of syringe configurations. 1. Introduction. The syringe is probably the most important and most indispensable instrument in medicine. The syringe consists of three parts, the needle, the glass cylinder and the plunger. As the plunger is squeezed, the fluid inside the glass cylinder is ejected through the hollow needle. The present paper is concerned with the low Reynolds number flow inside a syringe. This situation occurs in the slow controlled intravenous infusion of fluids and drugs in the hospital or laboratory. For example, typical vlaues for the Reynolds number in indicator-dilution experiments are Re= 0.07(50cc syringe at 1cc/min or Re= 0.015(10cc syringe at 0.1 cc/min). The infusion rates are kept constant by a calibrated infusion pump. In Figure 1 the plunger is moving to the left in a circular cylinder. Fluid is forced through a hollow needle of diameter usually less than 3% of the diameter of the cylinder. Due to the geometry the fluid dynamics can be separated into two parts: the flow in the needle and the flow in 0 Departments of Electrical Engineering and Computer Science, Industrial and Operations Engineering, and Mathematics, University of Michigan, Ann Arbor, MI 48109. 1 Department of Computer Science, Virginia Polytechnic Institute & State University, Blacksburg, VA 24061. 2 Departments of Mechanical Engineering and Physiology, Michigan State University, East Lansing, MI 48824. 3 Department of Computer Science, University of Wisconsin, Madison, WI 53706. 1

the cylinder. The flow in the needle (orifice) has been solved by Dagan, Weinbaum and Pfeffer [1], who concluded the velocity profile is dominantly Poiseuille, with end effects limited to a distance of only 1/4 the inner needle diameter. In what follows we shall study the cylinder region, with the needle hole approximated by a sink. 2. Formulation Let the origin of cylindrical coordinates (r, z) be situated at the sink. The cylinder is at z' = 0, r' < a and z' > 0, r' = a. At the time of investigation, the plunger is at z' = 1, r' < a traveling with velocity U. Let v' and w' be the velocity components in the r' and z' directions, respectively. Since the Reynolds number is much less than one, the flow is governed by the quasi-steady axisymmetric Stokes equation: ( 2 1 a9 2 2 (^2I- OrW) 2 01, (1) where bt' is a stream function defined by, 1 9b', 1 a' v - I ' w, (2) v = - ^'' O r'' The boundary conditions for the domain 0 < zt < 1, 0 < r' < a are: on z = O, vt = w' = O, except at the sink at the origin, z'= 1, v' = O w' =-U r = a, v = w = 0 (3) r' = 0, = = o rl We normalize all lengths by a, the velocities by U, the stream function by - a2 U and drop primes. The governing equations become ( a2 I 82 2 t —+ a- =0, Or2 rdr O z2r on z = 0, ' = 1, = 0 except at the origin, I az84 18 a - =,a = 0, 9 ~= 2, r=0, 0,= 8O3r k4r)rI= a r or r r=1,,= 1, = 0 dr is the instantaneous ratio o the axial length of tf he cylinder to its radius. There exists a weak singularity on the circle z =,, r = 1 where the velocity w is discontinuous and a strong singularity 2

at the origin where b is discontinuous. These singularities pose immediate problems in solving for the stream function directly. To remove the strong singularity analytically we set )(r, z) = )o(r, Z) + x(r, z) (5) where b0 represents the Stokes flow due to a sink on an infinite solid plane [2]: o =1- z3(r2 + Z2)-3/2 (6) Thus 8( 2 18 )2X = 2 The boundary conditions are ax on z= 0, X = 0 = 0, (8) z on z=,=, X =r2 - 1 +P3(r2 +p2)-3/2 aX = 32 r2(r2 + 2)-5/2 (9) az onr=0, X=0 lim Or(rr) = 0, (10) on r= 1, X = 3(1 + 2)-/2, X = -3z3(1 + z2)-5/2 (11) Although Equations (7-11) are regular at the origin, a discontinuity in the derivative of X still exists on z = f3, r = 1. The boundary conditions are complicated enough such that an analytical solution is not worthwhile. In what follows we shall use numerical means to solve the problem. 3. Description of the numerical method The partial differential equation for X is solved numerically using finite difference methods [3], which produce an approximation for X at each grid point. In the middle away from all boundaries, standard difference equations are used to approximate the partial derivatives. Near the boundaries, different nonstandard formulas are needed. These new difference formulas, using both values and partials of X on the boundaries, were derived using the symbolic manipulation package SMP [4]. All the finite difference approximations used are given in Appendix A. The finite difference approximation results in a linear system of equations Ax = b, with dimension 2500 for a modest 50 x 50 grid. Since the boundary conditions are not symmetric and difference formulas using known derivatives on the boundaries are employed, the coefficient matrix A is not symmetric. The matrix A is nearly diagonally dominant, however, suggesting that successive over-relaxation (SOR) may work. SOR does not, in fact, converge for all choices of the relaxation parameter w anywhere between 0 and 2 (A is not symmetric and positive definite). Nevertheless, excellent convergence for SOR was achieved with w = 1.6. 4. The streamlines. After X (r, z) is computed, it is added to 4o to obtain the stream function t1. Figure (2) shows the streamlines for i/ = 5 or when the plunger is 5 radius' distance from the sink. For such large p values, the fluid flow shows 3 distinct regions: the plunger region at right, the middle Poiseuille 3

region and the sink region at left. The flow in the plunger region is similar to the entry flow into a circular tube [5] where the uniform velocity profile is changed into the parabolic profile. Our results agree with Lew and Fung's [5] analytical result and Wagner's [6] numerical result for the infinite tube at zero Reynolds number. The entrance length is about one diameter. The middle regions shows almost Poiseuille parallel flow. For comparison, the dashed lines are the Poiseuille streamlines = 2r2 -r4 (12) In the sink region the fluid turns and flows into the origin. The streamlines are dominated by the Stokes sink represented by 4'o. There are no regions of recirculation although the velocity is quite low in the corner at r = 1, z = 0. Figures 3-7 show the changes in streamline patterns as the plunger moves in closer. The three distinct regions disappear due to mutual interaction. At /3 = 0.25 the velocity becomes almost radially inwards in contrast to the mainly axial Poiseuille flow at /3 = 5. 5. The pressure distribution. The pressure is numerically computed by finite differences from the normalized Stokes formula: ap_ 1 ( a3 3 1 I a2 a —r = r - z\ar2- +-a-) Z3 r2 azar ' l'aS~ 031)z^T2+12) I r(13) ap _1 (a3 a3 1 a2b 1 a az r r3 + Qrz2" r2 dr2 r3 ar The constant pressure lines are shown in Figures 8-13. We see a high pressure region exists near the corner at z = 3, r = 1 while near the center of the plunger (z = P/, r = 0) the pressure shows a local high for large /3 and a local low for P/ < 1. In the vicinity of the sink, the pressure is high on the surface z = 0 and is low on the axis r = 0. The pressure pattern near the sink persists for all P/ and agrees with the Stokes sink pressure distribution [2]. Notice that for /3 = 5 there exists a middle Poiseuille region where pressure is independent of radial distance and the pressure gradient is uniform axially. 6. Discussion The Stokes flow in a syringe is now solved. Although the streamlines do not show recirculation, the flow in a syringe is not uniform. The pressure distribution is more complicated. Regions of local high and local low pressure may cause particulate matter in the fluid to migrate, resulting in inhomogeneity of the injectate. Stokes flow is not valid in the immediate neighborhood of the sink where velocities increase indefinitely. To eliminate the singularity, one must consider a hole of finite radius at the origin. Dagan, Weinbaum and Pfeffer [1] considered a finite cylindrical hole in an infinite plate and computed the streamlines and pressure distribution. They found the effect of the hole is limited to one hole radius in the vicinity of the junction. We conclude that the shape and size of the hole has little effect on the flow in the syringe presented in this paper. What is the force on the plunger? Due to the difference in size, almost all hydrodynamic resistance comes from the flow through the needle. Reference [1] showed the pressure difference 4

through a finite sized hole of radius ah and length Lh is a linear sum of Poiseuille resistance and Sampson resistance [7]. Thus the force on the plunger is a constant: F-= ra2 8Lh + 3) (14) w' ah 7. References [1] Dagan, Z., Weinbaum, S. and Pfeffer, R., "An infinite-series solution for the creeping motion through an orifice of finite length," J. Fluid Mech., Vol. 115, 1982, pp. 505-523. [2] Happel, J. and Brenner, H., Low Reynolds Number Hydrodynamics, Prentice Hall, N.J., 1965, Chap. 4. [3] Smith, G.D., Numerical Solution of Partial Differential Equations, Oxford U. Press, N.Y., 1965. [4] Wolfram, S., "Symbolic Mathematical Computation," Comm. ACM, Vol. 28, 1985, pp. 390 -395. [5] Lew, H.S., and Fung, Y.C., "On the low-Reynolds number entry flow into a circular cylindrical tube," J. Biomechanics, Vol. 2, 1969, pp. 105-119. [6] Wagner, M.H., "Developing flow in circular conduits: transition from plug flow to tube flow," J. Fluid Mech., Vol. 72, 1975, pp. 257-268. [7] Sampson, R.A., "On Stoke's Current Function," Phil. Trans. Roy. Soc. London, Vol. A182, 1891, pp. 449-518. FIGURE CAPTIONS Figure 1. Syringe geometry and coordinate system. Figure 2.,b streamlines for P = 5.0 (solid) and Poiseuille streamlines (dashed). Figure 3. t/ streamlines for / = 2.0. Figure 4. b streamlines for 3 = 1.0. Figure 5. 0b streamlines for / =.75. Figure 6. ~b streamlines for # =.50. Figure 7. b streamlines for / =.25. Figure 8. Constant pressure lines (equispaced values) for / = 5.0. Figure 9. Constant pressure lines (equispaced values) for / = 2.0. Figure 10. Constant pressure lines (equispaced values) for / = 1.0. Figure 11. Constant pressure lines (equispaced values) for / =.75. Figure 12. Constant pressure lines (equispaced values) for / =.50. Figure 13. Constant pressure lines (equispaced values) for / =.25. 5

L-: J% k\ L { \

i-I, 2 1 0 5.0

L I~~~~~~ 1% f^ iiiiii iiiiiii > *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~h

0

(y,

IMa O, ll" o \

tl 0I I r) u'

O'E 0 V- 2

I I I 1 I1 I 1I I I I I I/h |' H \ I> I I I I I 0

0 0 C)~~~~~~~~~~~~~~~

I 0 ^i * L4 mmm Ul~~~~~~~~~~~~~~~

0 | I I I a ' I I I I I I I I a I I I I I I I I 01 Ul~~~~~~~~~~

-"I C,? t!

APPENDIX A FINITE DIFFERENCE FORMULAS For the first derivative with respect to r at all points, d(CHI(z,r)) dr CHI(z,r+h) - CHI(z,r-h) 2h For the second partial with respect to r at all points, 2 d (CHI(z,r)) 2 dr CHI(z,r+h) - 2CHI(z,r) + CHI(z,r-h) = --- —— _ --- —--- --------------— ____ 2 h For the third partial with respect to r in the middle away from the boundaries r=O and r=l, 3 d (CHI(z,r)) d dr CHI(z,r+2h) - 2CHI(z,r+h) + 2CHI(z,r-h) - CHI(zr-2h) 3 2h Near the boundary r=O (i.e. at points (z,h) ) 3 d (CHI(z,r)) 3 dr d(CHI(z,r-h)) ) -3CHI(z,r) + 2h ------------- + 8 CHI(z,r-h) ____________ CHI(z,r+2h) + -._ ----. --- dr 3 3 3 h

Near the boundary r=l (i.e. at points (z,l-h) ) 3 d (CHI(z, r)) 3 dr d(CHI(z,r+h)) 3CHI(z,r) + 2h ----------- - dr 8 CHI(z,r+h) CHI(z,r-2h) -— _ —_ ---_ --- - — _ --- — 3 3 3 h For the fourth partial with respect to r away from the boundaries r=O and r=l. 4 d (CHI(z,r)) 4 dr CHI(z,r-2h) - 4CHI(z,r-h) + 6CHI(z,r) - 4CHI(z,r+h) + CHI(z,r+2h) 4 h Near r=O 4 d (CHI(z,r)) -------— _ --- = 4 dr CHI(z,r+3h) 16CHI(zr) - -------- 4 8CHI(z,r+2h) + ------------ -9CHI(z,r+h) 3 113CHI(z,r-h) d (CHI(z,r-h)) - ------------- - 5h -------------- 12 dr 4 h

Near r=l, 4 d (CHI(z,r)) 4 dr CHI(z,r-3h) 8CHI(z,r-2h) 16CHI(z,r) - ---- ---- - ---------- - 9CHI(z,r-h) 4 3 113CHI(z,r+h) d (CHI(z,r+h)) - ------------ + 5h -------------- 12 dr h For the fourth partial with respect to z away from the boundaries z=O and z=beta. d (CHI(z,r)) -------------- = 4 dz CHI(z-2k,r) - 4CHI(z-k,r) + 6CHI(z,r) - 4CHI(z+kr) + CHI(z+2k,r) 4 k

Near z=O 4 d (CHI(z,r)) 4 dz CHI(z+3k,r) 8CHI(z+2k,r) 16CHI(z,r) - -------- + --------- - 9CHI(z+k,r) 4 3 113CHI(z-k,r) d (CHI(z-kr)). ---- ----. —. 5k -----—. ---12 dr k Near z=beta, d (CHI(zr)) 4 dz CHI(z-3k,r) 16CHI(z,r) - --------- 4 8CHI(z-2k,r) + ----------- - 9CHI(z-k,r) 3 113CHI(z+k,r) d (CHI(z+k,r)) - ------------ + 5k -------------- 12 dr 4 k

For the mixed fourth partial in the middle away from all boundaries. d (CHI(z,r)) _ --- —--------- = 2 2 dr dz 2 ( 30CHI(z,r) + CHI(z,r-2h) - 16CHI(z,r-h) - 16CHI(z,r+h) + CHI(z,r+2h) ) - ( 30CHI(z-kr) + CHI(z-kr-2h) - 16CHI(z-kr-h) - 16CHI(z-k,r+h) + CHI(z-k,r+2h) ) - ( 30CHI(z+kr) + CHI(z+kr-2h) - 16CHI(z+kr-h) - 16CHI(z+kr+h) + CHI(z+k,r+2h) ) _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ __ — ---- m_ _ _ _ _ _ 2 2 12 h k

Near z=O, 4 d (CHICz,r)) 2 2 dr dz -10CHI(z,r) 257CHI(z-k,r) -2 ( ~- +. ----_ ---- _+ 7CHI(z+k,r) 4 2CHI(z+2k,r) 9 3 144 CHI(z+3k,r) + -.-.-..-48 + 48 5k d CHI(z-k,r) _ — 12 -- ) 12 dz -1OCHI(z, r-h) 257CHI(z-kr-h) + ( -..-+ --- ——.- - + 7CHI(z+k, r-h) 4 2CHI (z+2k, r-h) 9 3 144 CHI(z+3k,r-h) + -.. --- —. --- + 48 5k d CHI(z-k,r-h) 12 dz- -- ) 12 dz -10CHI(z, r+h) + ( -------- 3 257CHI(z-k,r+h) +.-+. --- —----- + 144 7CHI(z+k,r+h) 4 2CHI(z+2k, r+h) 9 CHI(z+3k, r+h) + — _ --- —--- + 5k d CHI(z-k,r+h) ~ - ----— I —. — ) 48 12 dz 2 2 h k

Near z=beta, 4 d (CHI(z,r)) 2 2 dr dz -10CHI(z,r) -2 ( ------ + 3 257CHI(z+k,r) 144 7CHI(z-k,r) 4 2CHI(z-2k,r) 9 CHI(z-3k,r) + -- --- - 48 - 48 5k d CHI(z+kr) 12 dz -10CHI(z,r-h) 257CHI(z+k,r-h) + ( ------- + --------- + 7CHI(z-k,r-h) emmlllmlllllll - 2CHI(z-2kr-h) 9 3 144 CHI(z-3k, r-h) + ------------- - 48 5k d CHI(z+k,r-h) 12 dz ) -10CHI(z, r+h) + ( ---- -- 3 257CHI(z+k r+h) -+.. --- —----- -+ 144 7CHI(z-k,r+h) 4 2CHI(z-2k,r+h) 9 CHI(z-3k,r+h) + ---------- - 5k d CHI(z+kr+h) -m - - - - - - - - - - - - - - -m ) 48 12 dz 2 2 h k

Near r-O, 4 d (CHI(z.r)) 2 2 dr dz -10CHI(z,r) -2 ( -----—. + 3 257CHI(z, r-h) -— 144 + 144 7CHI(z,r+h) 4 2CHI(z, r+2h) 9 CHI(z,r+3h) 5h d CHI(z,r-h) 48 12 dr. + 48 12 dr -10CHI(z-k, r) + ( -.... — - 3 257CHI(z-k, r-h) + 44 --- —---- + 144 7CHI(z-k, r+h) 4 2CHI(z-k,r+2h) 9 CHI (z-k, r+3h) +.-+ --- —---- + 48 5h d CHI(z-k,r-h) ~- ----— ~ ---- ) 12 dr -10CHI(z+k, r) + ( ------- 3 257CHI z+k,r-h) 7CHI(z+k,r+h) + --------------- + ------------- - 2CHI(z+k,r+2h) 9 144 4 CHI(z+k,r+3h) +..8 --- —-- + 48 5h d CHI(z+k,r-h) 1 dr 12 dr 2 2 h k

Near r=l, 4 d (CHI(z,r)) 2 2 dr dz -1OCHI(z,r) -2 ( _. --- —- + 3 257CHI(z, r+h) -144- 144 7CHI(z,r-h) 4 2CHI(z,r-2h) 9 + CHI(z,r-3h) 48 5h d CHI(z,r+h) _ _12 )dr 12 dr -10CHI(z-k,r) 257CHI(z-kr+h) + --— + ( --- - + 3 144 7CHI(z-k r-h) 4 2CHI(z-k,r-2h) 9 CHI(z-k,r-3h) + 48 --- —---- - 4B 5h 12 d CHI(z-k,r+h) dr- ) dr -10CHI(z+k,r) 257CHI(z+k,r+h) + ( -------- + -----— __ ---- + 7CHI(z+k, r-h) _ _ _ _ _ _ _ _ _ _ _ _ _ - 2CHI(z+k,r-2h) 9 3 144 CHI(z+k,r-3h) — 48 — 48 5h d CHI(z+k,r+h) 12 dr — ) 12 dr 2 2 h k

For the mixed third partial in the middle asway from the boundaries z=O and z=beta. 3 d (CHI(z,r)) 2 dr dz 30CHI(z,r-h) + CHI(z-2k,r-h) - 16CHI(z-k,r-h) - 16CHI(z+k,r-h) +CHI(z+2k, r-h) + 30CHI(z,r+h) + CHI(z-2k,r+h) - 16CHI(z-k,r+h) - 16CHI (z+k, r+h) +CHI z+2k, r+h) 2 24 h k Near z=O, 3 d (CHI(z,r)) -------------- = 2 dr dz -1OCHI (z,r-h) _( _ —. --- — 3 257CHI(z-k, r-h) + --------------- + 144 7CHI(z+k, r-h) 4 2CHI(z+2k,r-h) 9 CHI(z+3k,r-h) +... --- —----- + 5k d CHI(z-k,r-h) --------- ) 12 dz -10CHI(z,r+h) 257CHI(z-kr+h) + ( --------- + ---—... --- —- + 7CHI(z+k,r+h) 4 2CHI(z+2k, r+h) 9 3 144 CHI(z+3k,r+h) + ~ —. --- —---- + 5k d CHI(z-k,r+h) ~- ---------- ) 48 12 dz 2 2h k

UNIVERSITY OF MICHIGAN 3 9015 04733 Near z=beta, 3 d (CHI(z,r)) 2 dr dz -10CHI(z,r-h) ( ----- 3 257CHI(z+k,r-h) +...4 --- —----- + 144 7CHI(z-k,r-h) 4 2CHI(z-2k,r-h) 9 CHI(z-3kr-h) — 48 — 48 5k d CHI(z+k,r-h) 1- -- -- )dz 12 dz -10CHI(z,r+h) + ( --------- 3 257CHI(z+k,r+h) + —. --- —--—._ + 144 7CHI(z-k, r+h) 2CHI(z-2k, r+h) 4 9 CHI(z-3kr+h) 5k d CHI(z+kr+h) + -------— _ — - -- --—... ---- ) 48 12 dz 2 2h k