Fluid Mechanics, 29 (1988) 337-346 Elsevier Science Publishers B.V.. Amsterdam - Printed in The Netherlands
Journal of Non-Newtonian
SINGULAR BEHAVIOR OF POWER-LAW IN HELE SHAW FLOW
337
FLUIDS
0. HASSAGER Instituttet for Kemiteknik,
Danmarks
Tekniske Hnjskole, DK-2800 Lyngby (Denmark)
and T.L. LAURIDSEN Nordiske Kabel- og Traaafabiker
A/S,
Vibeholms AU& 22, DK-2605 Brendby (Denmark)
(Received May 20, 1988)
Summary
The governing equation for Hele Shaw flow of power-law fluids is considered. The equation, which is elliptic for power-law index n > 0, is shown to have an equivalent variational integral. The solution exhibits singular behavior at comers. The singularity may be analyzed exactly, and is shown to be algebraic. Conditions under which the singularity is integrable are given.
1. Introduction
Hele Shaw flow refers to the flow between two parallel plates very close together. The boundary may take any closed shape and there may also be internal obstacles in the form of cylinders with generators perpendicular to the plates. At the boundaries either a pressure or zero normal velocity is imposed. The fluid is made to flow in the small space by differences in the pressure imposed at the boundaries. The motion was analyzed for Newtonian fluids by Hele Shaw [1,2] who noted that the streamlines become identical with those in a hypothetical two-dimensional flow with zero vorticity in a space of the same form. A particularly beautiful experimental illustration of Hele Shaw flow for Newtonian fluids, including the corner singularities has been given by Peregrine [3]. 0377”0257/88,‘$03.50
8 1988 Elsevier Science Publishers B.V.
338 The extension of the analysis of hele Shaw flow for Newtonian fluids to generalized Newtonian fluids was made by Hieber and Shen [4]. For power-law fluids the pressure is determined by a nonlinear Laplace-type equation. We here wish to demonstrate the existence of a variational principle for that equation and to investigate the local behavior of its solution near a sharp comer. 2. Analysis for power-law fluids We consider the flow in the space between two parallel planes separated by a distance S. The flow is described in a rectangular coordinate system with z-axis perpendicular to the planes, and with origin at the mid-plane. Thus the two bounding parallel planes are described by z = - 6/2 and z = 6/2. To derive the governing equation for Hele Shaw flow of power-law fluids we assume that over distances of order 6 in the x and y directions,
(1) (2) (3) (4)
u, = u,(z),
uy= u,(z), u, = 0, P’Pb,
Y>.
That is, it is assumed that the flow is locally that between infinite parallel planes. Then the equations of motion for a power-law fluid become
-=-
(5) (6)
where the power-law parameters are m and n, and the shear-rate is
(7) the boundary conditions to be used with these equations are
dux -z=Oatz=O d”u dz
and u,=u,=Oat
z=6/2.
(9)
Equations (5) and (6) are two coupled ordinary differential equations for u,
339 and uv. The equations may be uncoupled and solved for u, and u,,. We give here the results for the fluxes in the x and y directions defined by 2+1/n qx =
2
6’2v,
dz
=
_
(2 + l/n)m””
1 r=O
qu =
2
“2uy J z=o
(6/2)
dz
=
&!
ax ’
_
(10) (11)
where
This completes the derivation of the flux expressions for a given pressure gradient. Recall that eqns. (lo)-(12) were derived under the assumption that u, and u,, do not depend on x and y over distances of order 6. We continue to assume that over distances large compared to S the U, and ur may in fact depend on x and y but that they are given locally by the above analysis. A mass balance then shows that the fluxes in eqns. (10) and (11) must satisfy
&lx
x+-=0.
a% aY
(13)
Substitution of eqns. (10) and (11) in eqn. (13) gives a closed equation for the pressure in agreement with Hieber and Shen [4]. For simplicity we give the result for constant thickness independent of x and y:
Equations (12) and (14) give a closed equation for calculation of the pressure in flow between two parallel planes separated by a distance 8. The corresponding fluxes may be found subsequently by eqns. (10) and (11). Analysis shows that eqn. (14) not surprisingly is elliptic for y1> 0. Furthermore it may be shown [5] that the differential equation has the associated functional:
J=j$i$2+ (g)21”‘“;2n dV.
(15)
That is for given essential boundary conditions, the solution to eqn. (14) has the property that it minimizes the functional J defined in eqn. (15). Furthermore it may be seen that if the essential boundary conditions are specified only on a part of the boundary surrounding the region of interest,
340
then the natural boundary condition corresponding boundary is the zero normal component of the flux.
to the remaining
3. Finite element analysis In this section finite element approximation for p and (q,, q,,) are constructed. A Co approximation for p is assumed in the form N
(16) i=l
where the pi are nodal values and the &(x, JJ) are shape functions. We use here the bilinear quadrilateral Taig elements [6]. Introduction of eqn. (16) into eqn. (15) and differentiation with respect to a nodal value pk gives
1
(l-n)/2n
Aijpipj
where
A,,
dV, pm =o,
k= 1, 2 ,...,
N,
07)
With the inclusion of essential boundary conditions, eqns. (17) and (18) form a set of nonlinear algebraic equations for the pi that may be solved by iteration. The next problem to be considered is that of constructing approximations for the flux components q, and q,,. One might think of constructing these by introducing the approximation p into eqns. (10) and (11). That produces C-l approximations for the qx and qy, but we would like Co approximations. We therefore proceed by assuming Co approximations for the gradient of p: N
i=l N
gy=
C S:Si*
i=l
(20)
We then form residuals as follows:
iv
R,=gx-
ax,
(21)
R,=g,-
aF -a ay
(22)
341
co 0
d i
& d
1
0” d
-___-
i
I
,
---A-
Z d
---A--
---AC
0 d
8
,-
>‘/
I
-+-_^ _
4.00
-
,
0.01
_* 1
0.02
/ I
, c
_.
A
I
I
I
I
0.03
0.05
0.05
0.06
I 0.01
X
Fig. 1. Flow around a rounded corner (l/4 circle) for n =l.O. The magnitude and direction of the vectors indicate the total flow given by (q,, qy). The figure shows the increase in ( qx, qv) near the comer.
We now integrate Rz + Rt over the volume and minimize the result with respect to the parameters g: and gf to obtain the equations 5
Cijg;’ = /,$#ai
dV
j=l
and
where Cij=
+i+j dV. JV
(25)
Certainly the relations (23)-(25) may be obtained alternatively by applying Gale&ins method to eqns. (21) and (22). Once the approximate gradients g, and gr have been determined it is then straightforward to obtain Co approximations to the fluxes by substituting g, for ap/ax and g,, for ap/ay in eqns. (10) and (11). As an application of this technique we show in Fig. 1 a design for flow for n = 1.0 around a bend. The figure shows how the fluid velocity field
342
f
t
t
t
t
t
1
I
t
t
t
t
t
T
t
t
t
t
t
t
t
t
t
t
t
t
t
t
Itrrirt
Fig. 2. Illustration as in Fig. 1 for n = 0.55.
distributes unevenly near the bend. This distribution becomes more uneven as n decreases as shown in Fig. 2 for n = 0.55. 4. Analysis of comer singularities To analyse the singularity in p at a comer as shown in Fig. 3, we assume the following expansion [7] for the pressure near the corner: p=r”f(B)+r’q(8)+
....
(26)
We assume that the terms are ordered such that m < 1 etc. and look in the following just for the leading term f( 19).
(a)
(b)
Fig. 3. Flow near a sharp corner (a) symmetrical flow in the region - a < 19< antisymmetrical flow in the region - a c 0 < 0~.
a;
(b)
343 Introduction of eqn. (26) into eqn. (14) gives rise to the following nonlinear ordinary differential equation for f( 0): +
m(m
+
n
n
-
1)
x(1-n)/2n
f=O,
(27)
where
X=m2f2+
2
( 1 df
J-#
)
which is to be solved with the boundary conditions for either symmetric or antisymmetric flow, cf. Fig. 3: For symmetric flow: at (?=a:
df/de=O,
(29)
at 8=0:
df/de=O.
(30)
For antisymmetric flow: at 8=a:
df/de=O,
at e=O:
f=O.
m=-m,-
-- 1-n l+n’
(31)
(32) Note that the differential equations and boundary conditions are completely homogeneous in f, and that if f is a solution so is any multiple of f. Hence there will be a number of discrete eigenvalues m, I,. . . and corresponding eigenfunctions. Certainly the eigenfunctions must be such that the function in eqn. (15) is defined, and that imposes the following bound on m:
(33)
The problem will then be to find the smallest eigenvalue that satisfies the condition in eqn. (33). For n = 1 the problem is linear and it is easy to show that for symmetric flow
m, = r/a,
(34)
f, = cOs(ve/a),
(35)
and for antisymmetric flow
m, = ~/2a,
(36)
f, = sin( mB/2cw).
(37)
The solution of the nonlinear eigenvalue problem proceeds by a shooting technique as explained in the following. First we augment the boundary conditions for symmetric flow by at 8=0:
f=l,
(38)
and, for antisymmetric flow, at 8=0:
df/de=l.
(39)
344 TABLE 1 Eigenvalues m for antisymmetric flow around a 90 o reentrant comer ((Y= 311/4). Also listed is the exponent in the expression for qx - q,, - rs near the comer, and m, from eqn. (33) n
m
s
0.0100 0.1000 0.2000 0.4000 0.6000 0.8000 1.000 1.200 1.400 1.600 1.800 2.000
0.9904 0.9246 0.8734 0.8000 0.7460 0.7028 0.6666 0.6356 0.6084 0.5842 0.5624 0.5426
- 0.963 - 0.754 - 0.633 -0.500 - 0.423 - 0.372 - 0.333 - 0.304 - 0.280 - 0.260 - 0.243 - 0.229
m,
-
0.9802 0.8182 0.6667 0.4286 0.2500 0.1111 0.0000 0.0909 0.1667 0.2308 0.2857 0.3333
We have already noted that if f is a solution so is any multiple of f, and the above conditions merely serve to determine one particular choice for a multiplicative constant. In addition a value for m is guessed and eqn. (27) is integrated from 8 = 0 (with eqns. (30) and (38) for symmetric flow and eqns. (32) and (39) for antisymmetric flow) to 8 = (Yby a 4th order Runge-Kutta method. The value of m is then adjusted to satisfy the remaining condition (eqn. (29) for symmetric flow and eqn. (31) for antisymmetric flow). This adjustment is performed by Newton iteration which means that the sensitivity of the solution or its derivative with respect to m is also integrated. For example for antisymmetric flow the functions integrated would be f, g = df/de and h = ag/tbn, with initial conditions f(0)= 0, g(0)= 1, h(0) = 0. After integration to 8 = (Y an improved guess on m is obtained from m,+1= m, - g(a)/h(a)Values of m are shown in Table 1 for a 90” reentrant corner (ar = 37r/4) for antisymmetrical flow. The eigenvalues for symmetric flow are larger than for antisymmetric flow and any mixture of the two in the far field will therefore result in pure antisymmetric flow near the corner. Also shown in Table 1 is the exponent s = (m - 1)/n which determines the behavior of the total flow q = (q,2 + q;)l12 - rs near the corner. Evidently the flow is concentrated at the corner when n --*0. Also it appears that m + m, as n + 0. The asymptotic behavior of q near a sharp comer predicted from the local analysis may be compared with actual values found from the finite element analysis. Fig. 4 shows finite element values of q near a sharp comer as a function of the distance from the comer for n = 0.6, 1.0 and 1.6. The
345 100
I 50 -
I
I
q
10, - \\ 5- a&2
-
Fig. 4. Values of q = (qz + q;)l/’ near a sharp 90 o reentrant comer computed by the finite element method. A n -1.6, 0 n =1 and X n = 0.6.
approximate values of the exponent s in the expression 4 - T’ are called sFEM. These compare with values from the local analysis slocal (Table 1) as follows: n = 0.6, n = 1.0, n =
1.6,
sFEM= -0.44, SFEM= -0.35, SFEM = -0.25,
SPA = -0.423; Siocd= -0.333; slocal= -0.260.
The agreement between the local analysis and the finite element analysis is considered to be good, the deviation being due to the coarseness of the discretization. 5. Concluding remarks
The governing equation considered, eqn. (14) with eqn. (12), is invariant under a change of length scale. For this reason it is perhaps not surprising that a local solution near a corner may exhibit an algebraic singularity in cylindrical coordinates (i.e. the radius to a fractional power times some function of the polar angle). Similar considerations apply to the governing equations for 2-dimensional plane flow of power-law fluids, where the velocities and stress also exhibit algebraic singularities near sharp corners. It has been suggested [8] that the Oldroyd-B model also shows algebraic singularities in velocities and stresses near a sharp corner. Note, however, that in contrast to the power-law model, the Oldroyd-B model is not invariant under a change of length scale. Hence it is not possible to develop similarity solutions for the flow of Oldroyd-B fluids near corners. A satisfactory analysis of the flow of Oldroyd-B fluids near a sharp corner still remains to be developed.
346 Acknowledgement
One of the authors (TLL) wishes to acknowledge partly support by Akademiet for de Tekniske Videnskaber. The technique for the construction of a smooth gradient of a C ’ function was suggested to us by Dr. R. Keunings. References 1 H.J.S. Hele Shaw, Nature, 58 (1898) 34. 2 G.K. Batchelor, An Introduction to Fluid dynamics, Cambridge Univ. Press, 1967, p. 222. 3 D.H. Peregrine, in M. van Dyke “An album of fluid motion”, the Parabolic Press, Stanford, 1982, photo nr. 5. 4 C.A. Hieber and SF. Shen, J. Non-Newt. Fluid Mech, 7 (1980) l-32. 5 R.B. Bird, R.C. Armstrong and 0. Hassager, Dynamics of Polymeric Liquids, Vol 1, Fluid Mechanics, 2nd revised ed., Wiley, 1987, p. 775. 6 M.J. Crochet, A.R. Davies and K. Walters, Numerical Simulation of Non-Newtonian Flow, Elsevier, 1984, p. 201. 7 J.W. Hutchinson, J. Mech. Phys. Solids, 16 (1968) 13. 8 M.R. Apelian, R.C. Armstrong and R.A. Brown, J. Non-Newtonian Fluid Mech., 27 (1988) 299-321.