- Email: [email protected]
Analytic results for viscoelastic flow in a porous tube
Journal Mechanics, 28 (1988) 349-371 Else&r Science Publishers B.V., Amsterdam - Printed in The Netherlands
ANALYTIC RESULTS IN A POROUS TURE
FOR VISCOELASTIC
349
FLOW
RG. LARSON AT&T Bell Laboratories,
Muray
Hill New Jersey 07974 (U.S.A.)
(Received October 28,1987)
We consider the problem of an upper-convected Maxwell fled that is injected into or sucked from a cylindrical tube with porous walls. Many salient features of the flow are revealed by solving for the stresses assuming Newtonian kinematics. It is thereby found that the stress components have eigensolutions that in suction are eliminated by boundary conditions imposed on the centerline, which is the upstream boundary for suction. For injection, the centerline is the downstream boundary; boundary conditions applied at the centerline then do not eliminate the eigensolutions. As a result, numerical integration of the differential equations starting at the centerline is stable for suction, unstable for injection. This conclusion applies whether or not Newtonian kinematics are assumed. It is shown here that the eigensolution in injection is eliminated, and the integration stabilized, if integration is started at a position of zero radial flow outside the tube wall, in the region of “pre-history”, and carried out in the direction of flow, towards the centerline of the tube. The solution obtained in this way shows steep stress gradients that are driven by a biaxial extensional near-singularity in the region of pre-history.
1. Introduction Calculation of multi-dimensional viscoelastic processing flows, such as mold filling, cable jacketing, or film blowing, requires that one chooses both a constitutive equation that incorporates elastic effects and a numerical technique capable of solving the chosen constitutive equation coupled with the continuity and momentum balance equations. Unfortunately to date the numerical techniques used to solve these sets of nonlinear partial differential 0377-0257/88/$03.50
0 1988 Elsevier Science Publishers B.V.
350 Vr =
-I
k&-__,-___ Vr =$=O Fig. 1. The porous tube geometry with injection; from Kim-E [S].
equations have generally broken down when the elastic effects were made important. Numerical breakdown or loss of convergence is often associated with or accompanied by limit points [l], noisy or osciZZatory stresses or velocities [2,3] and steep boun&v Zuyers [3]. The origin of these effects and even whether they stem from the constitutive equation as opposed to the numerical technique is unclear. There is therefore a great need to develop analytic tools that can distinguishtypes of pathological behavior and trace them to their source in the constitutive equation or in the numerical technique. The behavior of constitutive equations in one-dimensionalflows can be tackled analytically and does yield some useful insight into possible sources of the numerical difficultiesencountered in the process geometries, but these possibilitiesseem now to be nearly exhausted. Needed are analytic results in multiple dimensions, or in axisymmetric geometries so that light can be shed on the difficulties associated with multi-dimensionality in process simulations. An important recent discovery is therefore the axisymmetric similarity structure developedby Brady and Acrivos [4] and applied by Kim-E [5] and by Menon et al. [6] to the porous tube problem with the upper-convected Maxwell constitutiveequation [5,6]. The porous tube geometry is depicted in Fig. 1. Fluid is uniformlyinjected into or sucked from an infinite cylinder of unit radius. A plane of symmetry perpendicular to the centerline at z = 0 separates fluid flowing to the right from that flowing to the left. The similaritystructure also holds in the planar analog of the axisymmetric tube, but this case was not analyzed. The upper-convected Maxwell (UCM) equation has been the most commonly used constitutive equation in numerical calculations, because of its relative simplicity and because it arises from the earliest and simplest molecular models. In dimensionlessform the UCM equation is De ?+T= -20, (1) where V-m T=T-Vv=~T-r~~v
351 is the upper-convected time derivative; i is the substantial time derivative. The Deborah number, De, can be thought of as the material relaxation time. The sign convention chosen for the stress tensor, T, is such that the momentum balance equation without body forces and inertia (the influence of inertia is considered in a forthcoming publication [7]) is vp+
v-7=0,
(2)
with p the pressure. D is the symmetric part of the velocity gradient tensor. The flow is assumed isochoric, i.e., v*u=o.
(3)
The upper convected Maxwell equation is not a particularly realistic equation. The dimensionless shear and first normal stress difference in simple shearing flow as functions of shear rate, 9, are . 712= -Y, (4) 711 -
722 =
-De
q2.
(5)
These dependencies on shear rate are usually only seen in viscoelastic materials at small p. In extensional flows, the equation is even more unrealistic; the normal stresses in uniaxial extension, where z is the direction of extension, are -2; 7ZL =l-2Dei’ where i = au,/az. In biaxial extension the normal stresses are -2i ?&= l-De<;
4i ‘ZZ= l+2Dei’
(7)
where in this case i = Elu,/gr. Therefore each extensional viscosity, (7,, ~,,)/i in uniaxial and (r,, - r,,)/i in biaxial extension, is singular at a critical extension rate. 2. SImIIarity solution Despite its inadequacies, the UCM equation is useful in that it is simple and does show pronounced elastic effects. For the porous tube flow, it is exceptional in that it admits as a solution the similarity form [5,6]: u, =fW/r,
u, = -zf’(
T,== k(r)z,
T,, = g(r) + @(r)z2,
7,, = i(r),
+,
% =_W
P = POW + iPd
(f-0
352
The dependencies of velocity and stress on z exhibited by this similarity solution are consistent with the continuity, momentum balance, and constitutive equations. The forms assumed for u, and u, satisfy the continuity equation for any f(r). 2.1 Ordinary differential equations to be solved The functions of r: f(r), k(r), g(r), h(r), i(r), j(r), p,,(r) and the constant p1 must be determined from a set of 7 equations arising from the momentum balance and constitutive equations of which 4 must be solved as a coupled set for k(r), h(r), i(r), and f(r); namely, r ( z component of momentum balance),
k’=
-PI-h-k/
fll=
fl. + kr+De(fk’--f’k+kf/r) l-Dei
(rz component of constitutive equation), h’=
(9)
i,i,:;
+4kj$)‘)
i’c~{2(~)‘[i--L]--L)
(zz component of constitutive equation),
00) (11)
( rr component of constitutive equation). (12)
Once the functions, f(r), h(r), i(r), and k(r), have been obtained, g(r), and p,(r), can be obtained by solving the remaining three equations, one at a time. For example, g(r) is obtained by solving
j(r),
2.2 Boundary conditions The boundary conditions on the kinematics for uniform injection and suction are f(1) = -1 injection, f (1) = 1 suction.
(13)
The no-slip boundary condition on the tube wall implies that f ‘(1) = 0.
(14) We shah occasionally cite results obtained by Kim-E [5], in which slip at the wall as allowed by replacing eqn. (14) with +f’=(l+)k.
05)
353 The slip coefficient /3, when set to unity, gives the no-slip condition, eqn. (14). Symmetry and continuity of both velocity components at the centerline imply that f(0) =f’(O)
= 0.
06)
If i’ and h’ are bounded at r = 0, the following boundary conditions can be derived from eqns. (9)-(12): k(0) = h(0) = 0, -f “(0) m= l_Def”(0)’
(17)
The linear dependence of u, on z and the z-independence of u, assure that the flow is a combination of uniuxial extension - for injection, biaxial extension - for suction, and shearing that increases linearly with z. From this and the linear dependence of Tagon p and quadratic dependence of 7 - 7,, on q2 in simple shearing flows of upper-convected Maxwell fluids, i;e might surmise that at large z in the porous tube, 7rZought to be linear, and 7,, - T,, quadratic in z. This is exactly the form the similarity solution takes. Other constitutive equations for which a similarity solution exists, namely the Oldroyd-B equation, the lower-convected Maxwell equation, other special cases of the Oldroyd Eight-Constant equation, and the second-order fluid, also have a linear dependence of shear stress and a quadratic dependence of the first normal stress difference on the shear rate. The simple behavior of these equations in shear therefore seems to be responsible for the existence of the similarity solution. More realistic constitutive equations are therefore not expected to have exact similarity solutions, although for other constitutive equations similarity forms that apply at asymptotically large z can be sought. The elongational character of the porous tube flow is exposed if the shear is removed from the problem by replacing the no-slip condition on the tube wall with perfect slip; i.e., zero shear stress. The injection problem is then equivalent to homogeneous uniaxial extension, and the suction problem is equivalent to biaxial extension, both at strain rates, i = 2. According to eqns. (6) and (7), the stresses for these flows exhibit simple poles at De = l/4 and l/2, respectively. 3. Review of the work of Kim-E and Menon et al. 3. I Stress singularity The resealing technique used by Kim-E and Menon et al., and discussed more detailed in Section 5.1, reduces the equations and boundary conditions
354 to an initial value problem. The scaled initial value problem was solved by an integration, using an implicit Gear integration routine, of the set of ordinary differential equations from t = 0, where the boundary conditions given by eqns. (15-(M) were applied, to t = 1, where the boundary conditions were satisfied by appropriately choosing the values of the scaling parameters. For suction, solutions were obtained for De up to 0.38. At De = 0.263, a singularity, namely a simple pole, appears in the solution for pi, and beyond De = 0.263, pi becomes aphysically negative. If perfect slip is introduced into the wall boundary condition, this simple pole shifts to De = 0.5; this is the value of De at which pure biaxial extension becomes singular. Thus the simple pole in the suction problem is closely related to the one in biaxial extension. 3.2 Numerical breakdown and limit points Based on these results, one might expect that in the injection case another. simple pole should be seen, one that shifts to De = 0.25 as slip is introduced. However, the integration from r = 0 fails to give any converged solution for De > 0. The finite difference technique, with boundary conditions applied at both r = 0 and r = 1, gave the same results for the suction problem as were obtained by the integration. For injection, finite-difference solutions were obtained by starting at De = 0 and then increasing De in small increments using arclength continuation, out to De = 0.073, where a limit point appeared. At a limit point, a family of solutions reaches a maximum Deborah number and bends back toward smaller values of De; see Fig. 2. If the constitutive equation, under conditions at which the limit point occurs,
I
!
I
I
-II-
I= -12 5 ; c
-
-13 -
Y g
-14-
# 2
-15-
B g
-16 -
I 0.02
I 0.04 DEBORAH
I 0.06 NUMBER,
I 0.06
D6
Fig. 2. The injection solution family, as computed by Runge-Kutta solid line), and by finite difference (Menon et al. [6] dashed line).
integration (this work,
355 represents real fluid behavior, then the limit point important transition in the flow.
corresponds
to an
3.3 Steep boundary layers and oscillations As the limit point was approached, large stress boundary layers developed near r = 1, especially in h(r). These boundary layers grew as the limit point was passed from below and solutions were obtained for ever lower values of De on the upper branch. On the upper branch, not only were the boundary layers steeper, but there were large oscillations in stresses on coarser meshes [5]; these were eliminated by ultra-refinement of the mesh [6]. Through all these variations in stresses, the velocity function, f(r), remained almost unchanged for the no-slip problem; it remained almost indistinguishable from its zero-Deborah-number limit. These results: limit points, steep boundary layers in stress, oscillations, and numerical breakdown, are banes familiar to practitioners of finite element analysis of similar problems. It is significant that the same pathology as is seen in finite element analysis is also exhibited by the much simpler semi-analytical techniques that were applied to porous tube problem by Kim-E and Menon et al. This suggests that a large part of the numerical difficulty lies in the constitutive equation; here the upper convected Maxwell equation. The simple pole that appears in the suction problem, for example, is clearly related to the constitutive equation. Before abandoning the UCM equation, however, it is important that its pathology be clearly understood so that symptoms of the same afflictions can more readily be spotted in other constitutive equations. In particular, we need to understand the differences between the injection and suction problems. Why can’t the injection problem be solved in resealed variables by the initial value integration? Why does the finite difference analysis of injection show a limit point? What causes the steep boundary layers? 4. Solutions using Newtonian kinematics To address these questions, we simplify the problem still more by taking advantage of the near-invariance of the computed velocity field throughout all calculations described above. Thus we shall compute the stresses from the constitutive equation using the Newtonian velocity field as the assumed kinematics. 4. I Kinematics For injection, the Newtonian solution for
f(r) =r4-2r2,
f,
namely
(18)
356 gives the velocities u, = ;f = -r(2
-2f’
ur = -
- t2),
= 4(1-
r
(1%
T2)Z.
(20)
From these kinematics one obtains explicit formulas for the coordinates, r(t; ro, ro, 0 and z(t; ro, zo, t’) in terms of time and the initial position:
(21)
z =
+z,r,” 1 +
2
- 1 exp(4(t - t’))
[i
I
1
* exp( -4(t
- t’)).
(22)
4.2 Lodge integral equation With the kinematics fixed, the stresses that satisfy the UCM equations can be obtained using the Lodge integral: Q = 1’ exp( - (t - t’)/De)C-‘(t, --cQ
t’) dt’.
(23)
The stress tensor, T, can be obtained from u by T= -(a-s)/De,
(24)
where 8 is the unit tensor. 4.2. I Finger tensor Stresses computed from the Lodge equation must satisfy the differential UCM equation. C- ’ is the Finger strain tensor, whose components are obtained from the kinematics: 2
c-l=
g
i cil=
1
($+
,
(g-t
G’=(g)(gJ
c,-,l=
r ’ . ( r0 1
(25)
357 These components can be calculated with eqns. (21) and (22). The results are Cr;‘=+[r2+(2-r2)exp(-4(t-t’))13exp(8(t-t’)), ~~,‘=16exp(-8(t-t’))[r~+(2-r’)exp(-4(t-t’))]-~ + 2.r2r2 [ r2 + (2 - r2) exp( -4( t - t’))]
(1 - exp(4( t - t’)))2,
Cr.l=+r,[r2+(2-r2)exp(-4(t-t’))12 X (exp(4( t - t’)) - exp(8(t - t’))), C;l=
f[r2 + (2 - r2) exp( -4(t
- t’))].
(26)
Since the kinematics of suction are those of injection with the flow direction reversed, the components of C-l for suction are obtained from eqns. (26) simply by changing each t - t’ to -(t - t’). 4.2.2 Pre-history To compute the stresses, the strains given in eqns. (26) must be multiplied by exp( -( t - t’)/De) and integrated over all times, t’ from - cc to the present time, t; see eqn. (23). For the injection problem, this requires that C-l be evaluated at times t’ when the particle had not yet entered the tube. Thus, for the integral equation, the particle pre-history must be provided before the stresses can be computed. To supply the needed pre-history for the integral equation, we analytically continue the kinematics given by eqns. (19) and (20) into the region outside the tube; that is, we apply these equations to the region r > 1. From eqn. (19), v, goes to zero at r = n. Figure (3) shows the Newtonian streamlines; there is a separating cylinder at r = fi. Thus all needed kinematics come from the region, r I a. With this particular pre-history, the velocity field is analytic inside and outside the tube and across the wall of the tube; the stresses computed from the Lodge integral are therefore analytic functions of De, at least near De = 0.
4.2.3 Injection: pole at De = l/8 For injection, C,S ‘, C;5 ‘, and C,T’ all grow as positive exponentials of 8( t - t’) at large (t - t’). From eqn. (23), the stresses, T,,, T,,, and rrZ therefore all become singular at De = l/8. The singularity in r,, is strongest at the wall. The singularities in T,, and TV=at large z are also strongest at the wall; these are caused by shearing components of the flow, although the elongational singularity in 7,,, which is strongest on the centerline, dominates at small z. In the numerical solution of the injection problem by finite difference with a two-point boundary condition, no simple poles were encountered; De = l/8 was not reached because of the appearance of a limit point at De - l/12.
: 0
I
I
I
1
J
1
2
2
4
5
1rx111l,
COORDINATE,
s
Fig. 3. Streamlines for Newtonian flow; the dashed and dotted lines show the positions of the tube wall and separating cylinder, respectively.
4.2.4 Suction: pole at De = I / 4
In suction, however, the simple pole found at De = 0.263 by both the two-point finite difference calculation and by the scaled initial value integration corresponds closely with the prediction of the Lodge model with Newtonian kinematics, namely a singularity at De = l/4. Exchanging (t - t ‘) for -(t - t ‘) in eqn. (26) gives the components of C-’ for suction. All four of these grow as exp(4( t - t’)) at large (t - t’); thus the corresponding stress components all show poles at De = l/4. For T,, and T,, this singularity is strongest at the wall; for ree and r,,, it is strongest at the centerline. Thus the singularities at De = l/4 and l/2 for biaxial and uniaxial extension occur at exactly twice the values of the Deborah number at which singularities occur in injection and suction, if the Newtonian kinematics are assumed. For suction, assumption of Newtonian kinematics gives almost precisely the behavior exhibited by the numerical calculations of the stresses with coupled non-Newtonian kinematics; for numerical calculation of the injection problem, a limit point intervenes at De - l/12, a limit point that is not predicted by the Lodge integral with Newtonian kinematics. 4.3 Upper-convected
Maxwell equation with Newtonian
kinematics
We again assume the Newtonian kinematics within the tube and solve the upper-convected Maxwell equation for the stress components. Any stress
)t
359
tensor obtained from the Lodge equation must be a solution the UCM equation; the latter equation can, however, have other solutions. The UCM equation can be solved sequentially: with f known,eqn. (12) can be solved for i; with f and i known, eqn. (10) can be solved for k; with f and k known, eqn. (11) can be solved for h. For each stress component, an inhomogeneous linear ordinary differential equation must be solved, for example,
(27)
i’ + P( r)i = Q(r),
(28) Q(r)=
(29)
-k--(f)‘.
4.3. I Eigensolutions The general solution to eqn. (27) is i=exp[-jPdx]/Q(x)
exp[jP(y)
dy] dx=c,
exp[-/Pdx],
(30)
where the integrals are all indefinite. The eigensolution,
=4 _r2)2-1/4Der2+1/2D (31)
comes from the solution to the homogeneous part of eqn. (27): i’ + P( r)i = 0,
(32)
where we are assumingthe Newtonian kinematics for injection given in eqn. (18). The suction case, which is related to injection by a change in the sign of f, can be obtained from eqn. (31) merely by changing De to -De. The unknownconstant, cl, multiplyingthe eigensolution,must be determined by the boundary condition. The boundary condition for i(0) was obtained from differential equation by assuming that i’(0) is bounded. But the eigensolution is always zero at r = 0 in injection, and is infinite in suction unless De is greater than l/4, the value at which the pole appears. Thus in suction, below the simple pole at De = l/4, the eigensolutionand its derivatives are unbounded; its coefficient must therefore be zero, since i’(0) is assumed finite. For injection, however, the coefficient of the eigensolution is undetermined. Furthermore, since for injection the eigensolutionand at least two of its derivatives are zero at r = 0, no reasonable boundary condition
360 applied to i at r = 0 will determine the This behavior is understandable: since in injection the flmd enters the tube at r = 1, it is there that the state of stress needs to be specified, just as with the Lodge intergral the prehistory of fluid entering the tube at r = 1 had to be specified. The eigensolutions for injection can be interpreted as the memory of the fluid as it enters the domain. As De decreases, the memory decays faster and so the eigensolution eigensolution becomes singular at the wall as De approaches zero. Therefore solution techniques that involve analytic or regular away from the Newtonian solution at De = 0 can avoid the indeterminate
(33) The result of the Lodge integration, eqn. (31), is obtained when c1 = 9/8. Added to the Lodge result is an eigensolution, r*/(2 - r’), with undetermined coefficient, which becomes large only at r = 1, and so adds a steep boundary layer there. Using the general solution for i at De = l/12, eqn. (10) can be solved for k: kc
24
(2 - r’)’
(r + cIr9 + c2r7).
(34)
k has a term, c,r9/(2
- r2)2, that comes from the eigensolution for i, as well as its own eigensolution, c2r7/(2 - r2)2. With k known, the linear ordinary
differential for h can be solved: h=
(2~1~~j3 ( r2 - clrlo - 2c2r8 + c&j),
(35)
which has terms coming from the eigensolutions for i and k as well as its own eigensolution, c//(2 - r2)3. This solution of the UCM equation with Newtonian kinematics does not exactly satisfy the momentum balance equation, eqn. (9), for any choice of the constants, cr, c2, and c3. Had we chosen a smaller De, say De = l/16, the eigensolutions would have involved even higher powers of r.
361
4.2.3 Eigensolutions in homogeneous flow Even for homogeneous flows on an unbounded domain there are eigensolutions for the UCM equation. For example, in homogeneous elongation, where the Lodge integral gives the stress components of eqn. (6), the UCM equation has these solutions with eigensolutions added to them: 7rr
-2i
1 *“=l+Dei
i %8=
+ qr
=l-2De1
l+De{
-4+2/Dei
+
c r2+2/Dei 3
+
c
4
+
c
2
=Zr2/Dei
3
,
(36b)
r2+2/Dei
These satisfy momentum balance if c,=Oand
c4=(3+&)c3.
The eigensolutions can be eliminated on a bounded domain by an appropriate choice of boundary conditions, or on an unbounded domain by requiring the stresses to be homogeneous, or by requiring the solution to be an analytic continuation of the zero-Deborah-number solution. Other solutions of the UCM equation that are not obtained from the integral are the aphysical stresses obtained in homogeneous elongation when De i > l/2. (These solutions are not on the solution branch passing through De = 0.) The Lodge integral for these kinematics is unbounded. 4.3.4 Implications for boundary conditions Thus the solution for injection to the UCM equation with Newtonian kinematics contains the results of the integral equation plus eigensolutions with undetermined coefficients. For the suction problem, the coefficients are determined by requiring bounded stress derivatives at r = 0, but for injection, this is not a strong enough condition to determine them. In general, it seems that one needs to impose conditions on the stress at the upstream boundary of the domain. For suction this is at r = 0, where the boundary conditions were applied in the numerical techniques of Kim-E, but for injection the upstream boundary is at r = 1. There is thus a history effect that obviously must be reckoned for in the integral model, but takes on a disguised form of upstream stress boundary conditions in the differential model. This phenomenon has been encountered by others, for example in the problem of drawing a filament from the end of a die [8].
362 5. In~tion
starting upstreamfor
iqjection
If our hypothesis is correct, the scaled integration technique of Kim-E and Menon et al. can be made to work by choosing an upstream starting point, and integrating in the direction of frow. The tube wall is not an acceptable starting point, because the stresses are not known there. However, as noted in the newtonian kinematics, there is a position outside the tube where f = 0. Here, as is true at the centerline, the stress boundary conditions can be obtained from the differential equations themselves, if the stresses are assumed to have bounded first derivatives. We shall assume that such a position exists for the non-Newtonian kinematics, and begin the integration there. 5.1 Resealing Kim-E and Menon et al. found that eqns. (9)-(12) can be resealed if one defines t = De’/*pi/*r, f(r) = De-*p,‘F(
p* =pI/I 0,
k(r) = De-“*p~‘*K(~),
p1 I,
h(r) =p,HW, i(r) = De-‘I(6).
(37)
Note that p* = f 1. Using these definitions of the resealed functions, H(6), K(6), and 40, eons. (9)-(12) become
F( 0,
K’=
-p*-H-K/&
(3W
F,,
F’ : -Kt+ 6
F’K-F(K’+ I-l
K/t)
,
(3W
H’=f(4K(F’,t-F//)-HE), 1’=+{2(1-l)(F’-F,<)-It}. The prime denotes differentiation with respect to 1. The resealing has eliminated De and p1 from the problem. The no-slip boundary condition at the wall transforms to F’(6,)
= 0,
or, if there is slip ,to
(39)
363 where 6 = I, corresponds to r = 1. Once eqns. (38) have been solved subject to appropriate boundary conditions, and a point, & at which I;‘(&) = 0 is found, the Deborah number and p1 can be obtained from the scaling relationships: De=lF(E,)/eZ(;
P15
(41)
r,“/mJ
5.2 Bounaluy conditions We seek a solution to eqn. (38) for which the radial velocity becomes zero for some to > I,; i.e., F(&) = 0. When this condition is inserted into eqns. (38), the boundary conditions at &, on the stresses can be obtained:
%J = to K(I
)
0
=
-
2F’(I‘,) _ 2F’([,)
(424
’
[F”(lO)EO - F'(toW(Io)
- 1)
to(F’t~o)- 60)
H(Eo) = 4~(Io)(~‘(I,)
(42b)
’
- F”(IO)IO)&.
(42c) Determination of these stresses therefore requires specification of F’(Eo) and F”(<,). The two degrees of freedom introduced by the need to specify these values are canceled by the two constraints arising from the boundary conditions at 6 = 0: F(O)=F’(O)
=o.
5.3 Transformation to a set of first or&r ODE’s This problem can be solved by Runge-Kutta formed to a set of first order ODE’s of the form
integration
if it is trans-
(43) where F is a vector function of the set of dependent variables, y, and the independent variable 6. The transformation to a problem of this kind can be carried out by augmenting the set of dependent variables with F’(t) and by using eqn. (38a) to substitute for K’ appearing in eqn. (38b), and eqn. (38b) to substitute for F” in eqn. (38~). 5.4 Perturbation expansions To proceed with the Runge-Kutta (RK) integration, we pick a value of to where we apply the boundary condition, F(6,) = 0. Note that eqns. (38~)
364 and (38d) become singular as I; + 0. To begin the RK integration, therefore expand each function about 6 = & in a Taylor’s series:
we
~=f,+fi(E-Eo)+f2(1-Eo)2+~(~-~~)3+f4(~-Eo)4, H=ho+h,(~-50)+h*(~-Eo)l, I=i,+il(E-~0)+i2(~-50)*, K=
k,
+
k,(5
-
6)
+
k2(6
-
&id*.
w
Using the differential equations, (38), we find the values of the coefficients in terms of fO, fi, and f2. From the boundary condition, F( &,) = 0, we find that f. = 0. There remain two unknown coefficients, fi and f2, corresponding to F’( &J and F”(&,). Th ese must be chosen so that the boundary conditions on the centerline, namely F(0) = 0 and F’(0) = 0, are satisfied. Our procedure, then, is to repeatedly guess values of F’(&) and F”( &,), and integrate from E = [,, towards 6 = 0, until values of F’( &,) and F”( to) are found that allow the centerline boundary conditions to be satisfied. In practice, the differential equations become infinitely stiff as t = 0 is approached, because I; must tend toward zero if a solution is to be found. Therefore, instead of integrating to E = 0, we integrate to a small value of E, and require that I;, H, I, and K match the values obtained from a perturbation expansion about t = 0: j7=
$2
+
p*(t6-
1)
54,
H = h,[* + h4t4, I = -/-y
+ i,.(* + i4i4,
K = k,( + k3t3,
(4%
where the values of the above coefficients are given by Kim-E [S] and by Larson [7].
5.5 Results inside the separating cylinder Starting the integration at the separating cylinder and integrating in the flow direction allows satisfaction of the boundary conditions and gives results that are insensitive to RK step size, A& when the step size is 0.0005 or less. Kim-E was not able to obtain such converged results when he integrated against the flow direction. Figure (4) shows the values of F’(&,) and F”(&) which, when chosen as initial conditions for the RK integration, lead to values of F, H, I, and K
365 ‘.’
0
I
I
I
I
(al
0.4
0.6
1.2
Fig. 4. F’(&,) and F”(&) to, for injection.
1.6
as functions of the resealed position of the separating cylinder,
that satisfy eqns. (45) near the centerline of the tube. A limit point at to = 1.32 is evident. The upper solution branch asymptotically approaches the line, F’(&,) = &,. On this line, I(&,) becomes singular according to eqn. (42a), and therefore K(&) and II(&) also become singular; see eqns. (42b) and (42~). This singularity is associated with the biaxial extensional singularity of the UCM equation, and is approached ever more close as the Deborah number is decreased along the upper solution branch! Note from Fig. (4) that P’(&,) is becoming larger as this singularityis approached. For each solution, there is a point, [, < &,, where we find F’(&) = 0, the no-slip wall boundary condition. From this ES, De and p1 can be obtained from eqn. (41). A plot of p1 against De in Fig. (2) shows a limit point at De - 0.08. The value of &, at this limit point in Deborah number is co - 1.32; thus the limit point in De occurs at or near the limit point shown in Fig. (4) for to. This limit point also occurs at or near the limit point in A = F”(O), which is A - -0.31. Thus, to within numerical accuracy, the limit points in lo, De, and A all occur simultaneously.Kim-E [5,6] using finite difference, found a limit point at a. somewhat different Deborah number, De - 0.073; see Fig. (2). The differencebetween the value obtained by Kim-E and that reported here apparently arises because the discretization in Kim-E’s finite difference technique implicitly set stress boundary conditions at the tube wall that are different from those that are set by the flow outside the tube [6]. Figure (5) shows f(r) and -f’(r)/ r, which is proportional to the axial velocity, both computed by RK integration for &, = 0.8 on the upper
366
I(-
N
-205
0.0
I
0.2
I 0.4
I I
0.6
I
I
0.6 RAolus,
1.0
I
12
I
1.4
6 1.1
r
Fig. 5. (a) Velocity simihuity function, - f(r), and (b) axial velocity at z = 1 for De = 0.031 on the upper solution branch. The dashed lines are the correspondingNewtonian results.
branch. As Kim-E found, the velocity inside the tube is nearly Newtonian. Just outside the tube, however, the velocity becomes highly non-Newtonian. Kim-E observed steep boundary layers in the stress similarity functions, especially h, near the tube wall [5,6]. Figure 6(a) shows that the steep boundary layer in h just inside the tube wall is dwarfed by the growth of h between the tube wall and the separating cylinder. The growth appears steep even on a plot of log h versus r, Fig. 6(b). The existence of two injection solutions for values of -A less than 0.31 raises a question concerning the perturbation expansion about t = 0, the first few terms of which are given in eqn. (45). The series have only one parameter, A. Once it is given, all coefficients of eqn. (45) can be obtained through recursion relations. But for each A, there are two solutions. To which solution, then, does the perturbation expansion converge? Close inspection of the two numerical solutions as I approaches zero shows that they approach each other, but neverthelessare not identical until 6 actually reaches zero. If both solutions are analytic near 6 = 0, there can be only one conclusion: the perturbation expansion about t = 0 must have a zero radius of convergence! In fact, a ratio test performed on successive coefficients of solution
367
I 0.0
I 0.2
I
I
0.4
0.6 RADIUS.
0.6
I
1.0
I!2
I: 4
r
Fig. 6. 7zz normal force hilarity function, h, as a function of radial position, r, on (a) a linear plot, and (b) a semi-log plot, showing the steep boundary layer near the tube wall and between the tube wall and separating cylinder.
the perturbation expansion showed a radius of convergence that tends towards zero as higher order terms are ratioed. Thus, for 6 not equal to zero, the perturbation expansion cannot exceed a limiting accuracy no matter how many terms are included. Likewise, the instability of the numerical integration starting from < = 0 should not be too surprising,since for each value of -A < 0.31 chosen, there are at least two solutions to which the integration could conceivably converge. For each solution to the resealed problem, there is not only a no-slip solution, but also a slip solution for each b < 1. One need only search the no-slip solution for a different stopping condition, the one corresponding to eqn. (40), to find a solution corresponding to slip, that is to j3 < 1. As one decreases I, from the value at which F’(&) equals zero to { = 0; the slip coefficient, j3, varies from unity to zero. Thus each no-slip solution has a set of slip solutionsembeddedin it; the smaller value of the slip coefficient, the smaller, in resealed units, is the diameter of the embedded tube with the slip wall condition. The existence, for each value of -A < 0.31, of two no-slip solutions that remain distinct for t > 0, therefore implies that for any /I < 1,
368 there are again two solutions for each value of -A less than -A,( 8). A,,(p) is the limit point for /3 < 1. In addition to these two slip solutions, Menon et al. [6] found a fully analytic solution for any A. The velocity solution of this solution is a pure homogeneous alongation and so cannot satisfy the no-slip wall boundary conditions but can satisfy this boundary condition when B -C1. The stresses of this analytic slip solution contain eigensolutionsthe coefficients of which are set by the slip coefficient and other wall stress boundary conditions. 5.6 Asymptotic boundary layer analysis Near the separating cylinder, 1, H, and K are all large, and one can drop the constant, p * = -1, and replace I- 1 by I in eqns. (38a)-(38d). The result is K’=
-H--K/&
[F;“-F’/[]I=
-KS+F’K-F(K’+K/E),
H’F = 4K( F’/t - F”) - H& I’F = 2I( F’ - F/e) - I& which has the solution, F=4x
1 eW-loz) - I (
ieZAcl* I = - Ce2xt’, 2t2A
(W
369
0
0.2
0.4
SQUARE
0.6
OF RADIAL
1.0
0.6 COORDINATE,
1
<*
Fig. 7. Semi-logarithmic plot of the similarity functions H, Z(*, K.$, and F’/t against [*, the square of the resealed radial coordinate, for .f,, = 1 on the upper solution branch. The straight-line portions of the lines agree with the asymptotic boundary layer formulas given by eqn. (47).
beyond the separating cylinder. Figure (8) shows F and the stress similarity functions, I, H, and K as functions of 6, for t > &, for &, = 0.8, De = 0.031, on the upper solution branch. The asymptotic formulas, eqn. (47), are almost exact in this region. The super-steepness of the boundary layers is evident: at t = 1.3 the renomalized stresses are all in excess of lo”!
I
0.6
0.9
I
1.0 RENORMALIZED
I I.2
I
I.1 RADIUS,
I.3
c
Fig. 8. Resealed stress and velocity similarity functions beyond the separating cylinder at De = 0.031 on the upper solution branch.
370 Codusions In flows of memory fluids, the upstream boundary of the domain can require special treatment. This is obvious when the memory effects are incorporated in an integral representation, for then a pre-history of the fluid up to the time it crosses the inflow boundary must be provided, but it can also be true for a differential equation, such as the upper-convected Maxwell equation. In the porous tube problem, with Newtonian kinematics specified, the Lodge integral, with pre-history supplied by analytically continuing the flow inside the domain, has simple poles for injection and suction. The differential upper-convected Maxwell equation with Newtonian kinematics and stress boundary conditions only on the centerline agrees with the Lodge integral for suction, but for injection has the Lodge solution plus eigensolutions with undetermined coefficients that are steep at the wall and become singular there as De goes to zero. These eigensolutions and their derivatives are singular at a cylinder outside the tube wall where flow separation occurs and the radial velocity is zero. Thus if one imposes boundary conditions, namely boundedness of stress derivatives, at this upstream separating cylinder, the eigensolutions are eliminated. Thus in either suction or injection, with Newtonian kinematics, eigensolutions are only avoided by applying boundary conditions at the upstream end of the domain. When non-Newtonian kinematics are incorporated by solving the constitutive equation and the momentum balance equation simultaneously, boundary conditions must again be specified at the upstream end of the domain, and, for stability, numerical integration of the equations must proceed in the direction of flow. This integration shows a limit point at De = 0.08. On the upper solution branch, the kinematics are almost Newtonian inside the tube, but are highly non-Newtonian between the tube wail and the separating cylinder. The stresses, which become steep as the tube wall is approached from inside the tube, become even steeper outside the wall, and steeper still out beyond the separating cylinder. The steepness of these stress boundary layers is produced by the biaxial flow outside the tube and the existence of a singularity in the biaxial extensional viscosity for the UCM fluid. An asymptotic solution is derived which is valid for the highly non-Newtonian kinematics and steep stress boundary layers near the separating cylinder and outside of it. On the upper solution branch this boundary layer penetrates into the porous tube. Acknowledgement I gratefully acknowledge stimulating and encouraging discussions with, and insight given by Miral Kim-E, John Brady, Bob Brown, and Eric Shaqfeh.
371 References 1 P.W. Ye&, M.E. Kim-E, RC. Armstrong and R.A. Brown, J. Non-Newtonian Fluid Mech., 16 (1984) 173. 2 M.J. Crochet and K. Walters, Ann. Rev. Fluid Mech., 15 (1983) 241. 3 M.A. Mend&on, P.-W. Yeh, RA. Brown and RC. Armstrong, J. Non-Newtonian Fluid Me&., 10 (1982) 31. 4 J.F. Brady and A. Acrivos, J. Fluid Me&., 112 (1981) 127. 5 M. Kim-E, Ph.D. thesis, Mass. Inst. of Technol., 1984. 6 R.K. Mcnon, M. Kim-E, RC. Armstrong, R.A. Brown and J.F. Brady, J. Non-Newtonian Fluid Mech., 27 (1988) 265. 7 RG. Larson, J. Fluid. Mech., accepted (1988). 8 M.M. Dem, C.J.S. Petrie and P. Avenas, AIChE J., 21 (1975) 791.
ANALYTIC RESULTS IN A POROUS TURE
FOR VISCOELASTIC
349
FLOW
RG. LARSON AT&T Bell Laboratories,
Muray
Hill New Jersey 07974 (U.S.A.)
(Received October 28,1987)
We consider the problem of an upper-convected Maxwell fled that is injected into or sucked from a cylindrical tube with porous walls. Many salient features of the flow are revealed by solving for the stresses assuming Newtonian kinematics. It is thereby found that the stress components have eigensolutions that in suction are eliminated by boundary conditions imposed on the centerline, which is the upstream boundary for suction. For injection, the centerline is the downstream boundary; boundary conditions applied at the centerline then do not eliminate the eigensolutions. As a result, numerical integration of the differential equations starting at the centerline is stable for suction, unstable for injection. This conclusion applies whether or not Newtonian kinematics are assumed. It is shown here that the eigensolution in injection is eliminated, and the integration stabilized, if integration is started at a position of zero radial flow outside the tube wall, in the region of “pre-history”, and carried out in the direction of flow, towards the centerline of the tube. The solution obtained in this way shows steep stress gradients that are driven by a biaxial extensional near-singularity in the region of pre-history.
1. Introduction Calculation of multi-dimensional viscoelastic processing flows, such as mold filling, cable jacketing, or film blowing, requires that one chooses both a constitutive equation that incorporates elastic effects and a numerical technique capable of solving the chosen constitutive equation coupled with the continuity and momentum balance equations. Unfortunately to date the numerical techniques used to solve these sets of nonlinear partial differential 0377-0257/88/$03.50
0 1988 Elsevier Science Publishers B.V.
350 Vr =
-I
k&-__,-___ Vr =$=O Fig. 1. The porous tube geometry with injection; from Kim-E [S].
equations have generally broken down when the elastic effects were made important. Numerical breakdown or loss of convergence is often associated with or accompanied by limit points [l], noisy or osciZZatory stresses or velocities [2,3] and steep boun&v Zuyers [3]. The origin of these effects and even whether they stem from the constitutive equation as opposed to the numerical technique is unclear. There is therefore a great need to develop analytic tools that can distinguishtypes of pathological behavior and trace them to their source in the constitutive equation or in the numerical technique. The behavior of constitutive equations in one-dimensionalflows can be tackled analytically and does yield some useful insight into possible sources of the numerical difficultiesencountered in the process geometries, but these possibilitiesseem now to be nearly exhausted. Needed are analytic results in multiple dimensions, or in axisymmetric geometries so that light can be shed on the difficulties associated with multi-dimensionality in process simulations. An important recent discovery is therefore the axisymmetric similarity structure developedby Brady and Acrivos [4] and applied by Kim-E [5] and by Menon et al. [6] to the porous tube problem with the upper-convected Maxwell constitutiveequation [5,6]. The porous tube geometry is depicted in Fig. 1. Fluid is uniformlyinjected into or sucked from an infinite cylinder of unit radius. A plane of symmetry perpendicular to the centerline at z = 0 separates fluid flowing to the right from that flowing to the left. The similaritystructure also holds in the planar analog of the axisymmetric tube, but this case was not analyzed. The upper-convected Maxwell (UCM) equation has been the most commonly used constitutive equation in numerical calculations, because of its relative simplicity and because it arises from the earliest and simplest molecular models. In dimensionlessform the UCM equation is De ?+T= -20, (1) where V-m T=T-Vv=~T-r~~v
351 is the upper-convected time derivative; i is the substantial time derivative. The Deborah number, De, can be thought of as the material relaxation time. The sign convention chosen for the stress tensor, T, is such that the momentum balance equation without body forces and inertia (the influence of inertia is considered in a forthcoming publication [7]) is vp+
v-7=0,
(2)
with p the pressure. D is the symmetric part of the velocity gradient tensor. The flow is assumed isochoric, i.e., v*u=o.
(3)
The upper convected Maxwell equation is not a particularly realistic equation. The dimensionless shear and first normal stress difference in simple shearing flow as functions of shear rate, 9, are . 712= -Y, (4) 711 -
722 =
-De
q2.
(5)
These dependencies on shear rate are usually only seen in viscoelastic materials at small p. In extensional flows, the equation is even more unrealistic; the normal stresses in uniaxial extension, where z is the direction of extension, are -2; 7ZL =l-2Dei’ where i = au,/az. In biaxial extension the normal stresses are -2i ?&= l-De<;
4i ‘ZZ= l+2Dei’
(7)
where in this case i = Elu,/gr. Therefore each extensional viscosity, (7,, ~,,)/i in uniaxial and (r,, - r,,)/i in biaxial extension, is singular at a critical extension rate. 2. SImIIarity solution Despite its inadequacies, the UCM equation is useful in that it is simple and does show pronounced elastic effects. For the porous tube flow, it is exceptional in that it admits as a solution the similarity form [5,6]: u, =fW/r,
u, = -zf’(
T,== k(r)z,
T,, = g(r) + @(r)z2,
7,, = i(r),
+,
% =_W
P = POW + iPd
(f-0
352
The dependencies of velocity and stress on z exhibited by this similarity solution are consistent with the continuity, momentum balance, and constitutive equations. The forms assumed for u, and u, satisfy the continuity equation for any f(r). 2.1 Ordinary differential equations to be solved The functions of r: f(r), k(r), g(r), h(r), i(r), j(r), p,,(r) and the constant p1 must be determined from a set of 7 equations arising from the momentum balance and constitutive equations of which 4 must be solved as a coupled set for k(r), h(r), i(r), and f(r); namely, r ( z component of momentum balance),
k’=
-PI-h-k/
fll=
fl. + kr+De(fk’--f’k+kf/r) l-Dei
(rz component of constitutive equation), h’=
(9)
i,i,:;
+4kj$)‘)
i’c~{2(~)‘[i--L]--L)
(zz component of constitutive equation),
00) (11)
( rr component of constitutive equation). (12)
Once the functions, f(r), h(r), i(r), and k(r), have been obtained, g(r), and p,(r), can be obtained by solving the remaining three equations, one at a time. For example, g(r) is obtained by solving
j(r),
2.2 Boundary conditions The boundary conditions on the kinematics for uniform injection and suction are f(1) = -1 injection, f (1) = 1 suction.
(13)
The no-slip boundary condition on the tube wall implies that f ‘(1) = 0.
(14) We shah occasionally cite results obtained by Kim-E [5], in which slip at the wall as allowed by replacing eqn. (14) with +f’=(l+)k.
05)
353 The slip coefficient /3, when set to unity, gives the no-slip condition, eqn. (14). Symmetry and continuity of both velocity components at the centerline imply that f(0) =f’(O)
= 0.
06)
If i’ and h’ are bounded at r = 0, the following boundary conditions can be derived from eqns. (9)-(12): k(0) = h(0) = 0, -f “(0) m= l_Def”(0)’
(17)
The linear dependence of u, on z and the z-independence of u, assure that the flow is a combination of uniuxial extension - for injection, biaxial extension - for suction, and shearing that increases linearly with z. From this and the linear dependence of Tagon p and quadratic dependence of 7 - 7,, on q2 in simple shearing flows of upper-convected Maxwell fluids, i;e might surmise that at large z in the porous tube, 7rZought to be linear, and 7,, - T,, quadratic in z. This is exactly the form the similarity solution takes. Other constitutive equations for which a similarity solution exists, namely the Oldroyd-B equation, the lower-convected Maxwell equation, other special cases of the Oldroyd Eight-Constant equation, and the second-order fluid, also have a linear dependence of shear stress and a quadratic dependence of the first normal stress difference on the shear rate. The simple behavior of these equations in shear therefore seems to be responsible for the existence of the similarity solution. More realistic constitutive equations are therefore not expected to have exact similarity solutions, although for other constitutive equations similarity forms that apply at asymptotically large z can be sought. The elongational character of the porous tube flow is exposed if the shear is removed from the problem by replacing the no-slip condition on the tube wall with perfect slip; i.e., zero shear stress. The injection problem is then equivalent to homogeneous uniaxial extension, and the suction problem is equivalent to biaxial extension, both at strain rates, i = 2. According to eqns. (6) and (7), the stresses for these flows exhibit simple poles at De = l/4 and l/2, respectively. 3. Review of the work of Kim-E and Menon et al. 3. I Stress singularity The resealing technique used by Kim-E and Menon et al., and discussed more detailed in Section 5.1, reduces the equations and boundary conditions
354 to an initial value problem. The scaled initial value problem was solved by an integration, using an implicit Gear integration routine, of the set of ordinary differential equations from t = 0, where the boundary conditions given by eqns. (15-(M) were applied, to t = 1, where the boundary conditions were satisfied by appropriately choosing the values of the scaling parameters. For suction, solutions were obtained for De up to 0.38. At De = 0.263, a singularity, namely a simple pole, appears in the solution for pi, and beyond De = 0.263, pi becomes aphysically negative. If perfect slip is introduced into the wall boundary condition, this simple pole shifts to De = 0.5; this is the value of De at which pure biaxial extension becomes singular. Thus the simple pole in the suction problem is closely related to the one in biaxial extension. 3.2 Numerical breakdown and limit points Based on these results, one might expect that in the injection case another. simple pole should be seen, one that shifts to De = 0.25 as slip is introduced. However, the integration from r = 0 fails to give any converged solution for De > 0. The finite difference technique, with boundary conditions applied at both r = 0 and r = 1, gave the same results for the suction problem as were obtained by the integration. For injection, finite-difference solutions were obtained by starting at De = 0 and then increasing De in small increments using arclength continuation, out to De = 0.073, where a limit point appeared. At a limit point, a family of solutions reaches a maximum Deborah number and bends back toward smaller values of De; see Fig. 2. If the constitutive equation, under conditions at which the limit point occurs,
I
!
I
I
-II-
I= -12 5 ; c
-
-13 -
Y g
-14-
# 2
-15-
B g
-16 -
I 0.02
I 0.04 DEBORAH
I 0.06 NUMBER,
I 0.06
D6
Fig. 2. The injection solution family, as computed by Runge-Kutta solid line), and by finite difference (Menon et al. [6] dashed line).
integration (this work,
355 represents real fluid behavior, then the limit point important transition in the flow.
corresponds
to an
3.3 Steep boundary layers and oscillations As the limit point was approached, large stress boundary layers developed near r = 1, especially in h(r). These boundary layers grew as the limit point was passed from below and solutions were obtained for ever lower values of De on the upper branch. On the upper branch, not only were the boundary layers steeper, but there were large oscillations in stresses on coarser meshes [5]; these were eliminated by ultra-refinement of the mesh [6]. Through all these variations in stresses, the velocity function, f(r), remained almost unchanged for the no-slip problem; it remained almost indistinguishable from its zero-Deborah-number limit. These results: limit points, steep boundary layers in stress, oscillations, and numerical breakdown, are banes familiar to practitioners of finite element analysis of similar problems. It is significant that the same pathology as is seen in finite element analysis is also exhibited by the much simpler semi-analytical techniques that were applied to porous tube problem by Kim-E and Menon et al. This suggests that a large part of the numerical difficulty lies in the constitutive equation; here the upper convected Maxwell equation. The simple pole that appears in the suction problem, for example, is clearly related to the constitutive equation. Before abandoning the UCM equation, however, it is important that its pathology be clearly understood so that symptoms of the same afflictions can more readily be spotted in other constitutive equations. In particular, we need to understand the differences between the injection and suction problems. Why can’t the injection problem be solved in resealed variables by the initial value integration? Why does the finite difference analysis of injection show a limit point? What causes the steep boundary layers? 4. Solutions using Newtonian kinematics To address these questions, we simplify the problem still more by taking advantage of the near-invariance of the computed velocity field throughout all calculations described above. Thus we shall compute the stresses from the constitutive equation using the Newtonian velocity field as the assumed kinematics. 4. I Kinematics For injection, the Newtonian solution for
f(r) =r4-2r2,
f,
namely
(18)
356 gives the velocities u, = ;f = -r(2
-2f’
ur = -
- t2),
= 4(1-
r
(1%
T2)Z.
(20)
From these kinematics one obtains explicit formulas for the coordinates, r(t; ro, ro, 0 and z(t; ro, zo, t’) in terms of time and the initial position:
(21)
z =
+z,r,” 1 +
2
- 1 exp(4(t - t’))
[i
I
1
* exp( -4(t
- t’)).
(22)
4.2 Lodge integral equation With the kinematics fixed, the stresses that satisfy the UCM equations can be obtained using the Lodge integral: Q = 1’ exp( - (t - t’)/De)C-‘(t, --cQ
t’) dt’.
(23)
The stress tensor, T, can be obtained from u by T= -(a-s)/De,
(24)
where 8 is the unit tensor. 4.2. I Finger tensor Stresses computed from the Lodge equation must satisfy the differential UCM equation. C- ’ is the Finger strain tensor, whose components are obtained from the kinematics: 2
c-l=
g
i cil=
1
($+
,
(g-t
G’=(g)(gJ
c,-,l=
r ’ . ( r0 1
(25)
357 These components can be calculated with eqns. (21) and (22). The results are Cr;‘=+[r2+(2-r2)exp(-4(t-t’))13exp(8(t-t’)), ~~,‘=16exp(-8(t-t’))[r~+(2-r’)exp(-4(t-t’))]-~ + 2.r2r2 [ r2 + (2 - r2) exp( -4( t - t’))]
(1 - exp(4( t - t’)))2,
Cr.l=+r,[r2+(2-r2)exp(-4(t-t’))12 X (exp(4( t - t’)) - exp(8(t - t’))), C;l=
f[r2 + (2 - r2) exp( -4(t
- t’))].
(26)
Since the kinematics of suction are those of injection with the flow direction reversed, the components of C-l for suction are obtained from eqns. (26) simply by changing each t - t’ to -(t - t’). 4.2.2 Pre-history To compute the stresses, the strains given in eqns. (26) must be multiplied by exp( -( t - t’)/De) and integrated over all times, t’ from - cc to the present time, t; see eqn. (23). For the injection problem, this requires that C-l be evaluated at times t’ when the particle had not yet entered the tube. Thus, for the integral equation, the particle pre-history must be provided before the stresses can be computed. To supply the needed pre-history for the integral equation, we analytically continue the kinematics given by eqns. (19) and (20) into the region outside the tube; that is, we apply these equations to the region r > 1. From eqn. (19), v, goes to zero at r = n. Figure (3) shows the Newtonian streamlines; there is a separating cylinder at r = fi. Thus all needed kinematics come from the region, r I a. With this particular pre-history, the velocity field is analytic inside and outside the tube and across the wall of the tube; the stresses computed from the Lodge integral are therefore analytic functions of De, at least near De = 0.
4.2.3 Injection: pole at De = l/8 For injection, C,S ‘, C;5 ‘, and C,T’ all grow as positive exponentials of 8( t - t’) at large (t - t’). From eqn. (23), the stresses, T,,, T,,, and rrZ therefore all become singular at De = l/8. The singularity in r,, is strongest at the wall. The singularities in T,, and TV=at large z are also strongest at the wall; these are caused by shearing components of the flow, although the elongational singularity in 7,,, which is strongest on the centerline, dominates at small z. In the numerical solution of the injection problem by finite difference with a two-point boundary condition, no simple poles were encountered; De = l/8 was not reached because of the appearance of a limit point at De - l/12.
: 0
I
I
I
1
J
1
2
2
4
5
1rx111l,
COORDINATE,
s
Fig. 3. Streamlines for Newtonian flow; the dashed and dotted lines show the positions of the tube wall and separating cylinder, respectively.
4.2.4 Suction: pole at De = I / 4
In suction, however, the simple pole found at De = 0.263 by both the two-point finite difference calculation and by the scaled initial value integration corresponds closely with the prediction of the Lodge model with Newtonian kinematics, namely a singularity at De = l/4. Exchanging (t - t ‘) for -(t - t ‘) in eqn. (26) gives the components of C-’ for suction. All four of these grow as exp(4( t - t’)) at large (t - t’); thus the corresponding stress components all show poles at De = l/4. For T,, and T,, this singularity is strongest at the wall; for ree and r,,, it is strongest at the centerline. Thus the singularities at De = l/4 and l/2 for biaxial and uniaxial extension occur at exactly twice the values of the Deborah number at which singularities occur in injection and suction, if the Newtonian kinematics are assumed. For suction, assumption of Newtonian kinematics gives almost precisely the behavior exhibited by the numerical calculations of the stresses with coupled non-Newtonian kinematics; for numerical calculation of the injection problem, a limit point intervenes at De - l/12, a limit point that is not predicted by the Lodge integral with Newtonian kinematics. 4.3 Upper-convected
Maxwell equation with Newtonian
kinematics
We again assume the Newtonian kinematics within the tube and solve the upper-convected Maxwell equation for the stress components. Any stress
)t
359
tensor obtained from the Lodge equation must be a solution the UCM equation; the latter equation can, however, have other solutions. The UCM equation can be solved sequentially: with f known,eqn. (12) can be solved for i; with f and i known, eqn. (10) can be solved for k; with f and k known, eqn. (11) can be solved for h. For each stress component, an inhomogeneous linear ordinary differential equation must be solved, for example,
(27)
i’ + P( r)i = Q(r),
(28) Q(r)=
(29)
-k--(f)‘.
4.3. I Eigensolutions The general solution to eqn. (27) is i=exp[-jPdx]/Q(x)
exp[jP(y)
dy] dx=c,
exp[-/Pdx],
(30)
where the integrals are all indefinite. The eigensolution,
=4 _r2)2-1/4Der2+1/2D (31)
comes from the solution to the homogeneous part of eqn. (27): i’ + P( r)i = 0,
(32)
where we are assumingthe Newtonian kinematics for injection given in eqn. (18). The suction case, which is related to injection by a change in the sign of f, can be obtained from eqn. (31) merely by changing De to -De. The unknownconstant, cl, multiplyingthe eigensolution,must be determined by the boundary condition. The boundary condition for i(0) was obtained from differential equation by assuming that i’(0) is bounded. But the eigensolution is always zero at r = 0 in injection, and is infinite in suction unless De is greater than l/4, the value at which the pole appears. Thus in suction, below the simple pole at De = l/4, the eigensolutionand its derivatives are unbounded; its coefficient must therefore be zero, since i’(0) is assumed finite. For injection, however, the coefficient of the eigensolution is undetermined. Furthermore, since for injection the eigensolutionand at least two of its derivatives are zero at r = 0, no reasonable boundary condition
360 applied to i at r = 0 will determine the This behavior is understandable: since in injection the flmd enters the tube at r = 1, it is there that the state of stress needs to be specified, just as with the Lodge intergral the prehistory of fluid entering the tube at r = 1 had to be specified. The eigensolutions for injection can be interpreted as the memory of the fluid as it enters the domain. As De decreases, the memory decays faster and so the eigensolution eigensolution becomes singular at the wall as De approaches zero. Therefore solution techniques that involve analytic or regular away from the Newtonian solution at De = 0 can avoid the indeterminate
(33) The result of the Lodge integration, eqn. (31), is obtained when c1 = 9/8. Added to the Lodge result is an eigensolution, r*/(2 - r’), with undetermined coefficient, which becomes large only at r = 1, and so adds a steep boundary layer there. Using the general solution for i at De = l/12, eqn. (10) can be solved for k: kc
24
(2 - r’)’
(r + cIr9 + c2r7).
(34)
k has a term, c,r9/(2
- r2)2, that comes from the eigensolution for i, as well as its own eigensolution, c2r7/(2 - r2)2. With k known, the linear ordinary
differential for h can be solved: h=
(2~1~~j3 ( r2 - clrlo - 2c2r8 + c&j),
(35)
which has terms coming from the eigensolutions for i and k as well as its own eigensolution, c//(2 - r2)3. This solution of the UCM equation with Newtonian kinematics does not exactly satisfy the momentum balance equation, eqn. (9), for any choice of the constants, cr, c2, and c3. Had we chosen a smaller De, say De = l/16, the eigensolutions would have involved even higher powers of r.
361
4.2.3 Eigensolutions in homogeneous flow Even for homogeneous flows on an unbounded domain there are eigensolutions for the UCM equation. For example, in homogeneous elongation, where the Lodge integral gives the stress components of eqn. (6), the UCM equation has these solutions with eigensolutions added to them: 7rr
-2i
1 *“=l+Dei
i %8=
+ qr
=l-2De1
l+De{
-4+2/Dei
+
c r2+2/Dei 3
+
c
4
+
c
2
=Zr2/Dei
3
,
(36b)
r2+2/Dei
These satisfy momentum balance if c,=Oand
c4=(3+&)c3.
The eigensolutions can be eliminated on a bounded domain by an appropriate choice of boundary conditions, or on an unbounded domain by requiring the stresses to be homogeneous, or by requiring the solution to be an analytic continuation of the zero-Deborah-number solution. Other solutions of the UCM equation that are not obtained from the integral are the aphysical stresses obtained in homogeneous elongation when De i > l/2. (These solutions are not on the solution branch passing through De = 0.) The Lodge integral for these kinematics is unbounded. 4.3.4 Implications for boundary conditions Thus the solution for injection to the UCM equation with Newtonian kinematics contains the results of the integral equation plus eigensolutions with undetermined coefficients. For the suction problem, the coefficients are determined by requiring bounded stress derivatives at r = 0, but for injection, this is not a strong enough condition to determine them. In general, it seems that one needs to impose conditions on the stress at the upstream boundary of the domain. For suction this is at r = 0, where the boundary conditions were applied in the numerical techniques of Kim-E, but for injection the upstream boundary is at r = 1. There is thus a history effect that obviously must be reckoned for in the integral model, but takes on a disguised form of upstream stress boundary conditions in the differential model. This phenomenon has been encountered by others, for example in the problem of drawing a filament from the end of a die [8].
362 5. In~tion
starting upstreamfor
iqjection
If our hypothesis is correct, the scaled integration technique of Kim-E and Menon et al. can be made to work by choosing an upstream starting point, and integrating in the direction of frow. The tube wall is not an acceptable starting point, because the stresses are not known there. However, as noted in the newtonian kinematics, there is a position outside the tube where f = 0. Here, as is true at the centerline, the stress boundary conditions can be obtained from the differential equations themselves, if the stresses are assumed to have bounded first derivatives. We shall assume that such a position exists for the non-Newtonian kinematics, and begin the integration there. 5.1 Resealing Kim-E and Menon et al. found that eqns. (9)-(12) can be resealed if one defines t = De’/*pi/*r, f(r) = De-*p,‘F(
p* =pI/I 0,
k(r) = De-“*p~‘*K(~),
p1 I,
h(r) =p,HW, i(r) = De-‘I(6).
(37)
Note that p* = f 1. Using these definitions of the resealed functions, H(6), K(6), and 40, eons. (9)-(12) become
F( 0,
K’=
-p*-H-K/&
(3W
F,,
F’ : -Kt+ 6
F’K-F(K’+ I-l
K/t)
,
(3W
H’=f(4K(F’,t-F//)-HE), 1’=+{2(1-l)(F’-F,<)-It}. The prime denotes differentiation with respect to 1. The resealing has eliminated De and p1 from the problem. The no-slip boundary condition at the wall transforms to F’(6,)
= 0,
or, if there is slip ,to
(39)
363 where 6 = I, corresponds to r = 1. Once eqns. (38) have been solved subject to appropriate boundary conditions, and a point, & at which I;‘(&) = 0 is found, the Deborah number and p1 can be obtained from the scaling relationships: De=lF(E,)/eZ(;
P15
(41)
r,“/mJ
5.2 Bounaluy conditions We seek a solution to eqn. (38) for which the radial velocity becomes zero for some to > I,; i.e., F(&) = 0. When this condition is inserted into eqns. (38), the boundary conditions at &, on the stresses can be obtained:
%J = to K(I
)
0
=
-
2F’(I‘,) _ 2F’([,)
(424
’
[F”(lO)EO - F'(toW(Io)
- 1)
to(F’t~o)- 60)
H(Eo) = 4~(Io)(~‘(I,)
(42b)
’
- F”(IO)IO)&.
(42c) Determination of these stresses therefore requires specification of F’(Eo) and F”(<,). The two degrees of freedom introduced by the need to specify these values are canceled by the two constraints arising from the boundary conditions at 6 = 0: F(O)=F’(O)
=o.
5.3 Transformation to a set of first or&r ODE’s This problem can be solved by Runge-Kutta formed to a set of first order ODE’s of the form
integration
if it is trans-
(43) where F is a vector function of the set of dependent variables, y, and the independent variable 6. The transformation to a problem of this kind can be carried out by augmenting the set of dependent variables with F’(t) and by using eqn. (38a) to substitute for K’ appearing in eqn. (38b), and eqn. (38b) to substitute for F” in eqn. (38~). 5.4 Perturbation expansions To proceed with the Runge-Kutta (RK) integration, we pick a value of to where we apply the boundary condition, F(6,) = 0. Note that eqns. (38~)
364 and (38d) become singular as I; + 0. To begin the RK integration, therefore expand each function about 6 = & in a Taylor’s series:
we
~=f,+fi(E-Eo)+f2(1-Eo)2+~(~-~~)3+f4(~-Eo)4, H=ho+h,(~-50)+h*(~-Eo)l, I=i,+il(E-~0)+i2(~-50)*, K=
k,
+
k,(5
-
6)
+
k2(6
-
&id*.
w
Using the differential equations, (38), we find the values of the coefficients in terms of fO, fi, and f2. From the boundary condition, F( &,) = 0, we find that f. = 0. There remain two unknown coefficients, fi and f2, corresponding to F’( &J and F”(&,). Th ese must be chosen so that the boundary conditions on the centerline, namely F(0) = 0 and F’(0) = 0, are satisfied. Our procedure, then, is to repeatedly guess values of F’(&) and F”( &,), and integrate from E = [,, towards 6 = 0, until values of F’( &,) and F”( to) are found that allow the centerline boundary conditions to be satisfied. In practice, the differential equations become infinitely stiff as t = 0 is approached, because I; must tend toward zero if a solution is to be found. Therefore, instead of integrating to E = 0, we integrate to a small value of E, and require that I;, H, I, and K match the values obtained from a perturbation expansion about t = 0: j7=
$2
+
p*(t6-
1)
54,
H = h,[* + h4t4, I = -/-y
+ i,.(* + i4i4,
K = k,( + k3t3,
(4%
where the values of the above coefficients are given by Kim-E [S] and by Larson [7].
5.5 Results inside the separating cylinder Starting the integration at the separating cylinder and integrating in the flow direction allows satisfaction of the boundary conditions and gives results that are insensitive to RK step size, A& when the step size is 0.0005 or less. Kim-E was not able to obtain such converged results when he integrated against the flow direction. Figure (4) shows the values of F’(&,) and F”(&) which, when chosen as initial conditions for the RK integration, lead to values of F, H, I, and K
365 ‘.’
0
I
I
I
I
(al
0.4
0.6
1.2
Fig. 4. F’(&,) and F”(&) to, for injection.
1.6
as functions of the resealed position of the separating cylinder,
that satisfy eqns. (45) near the centerline of the tube. A limit point at to = 1.32 is evident. The upper solution branch asymptotically approaches the line, F’(&,) = &,. On this line, I(&,) becomes singular according to eqn. (42a), and therefore K(&) and II(&) also become singular; see eqns. (42b) and (42~). This singularity is associated with the biaxial extensional singularity of the UCM equation, and is approached ever more close as the Deborah number is decreased along the upper solution branch! Note from Fig. (4) that P’(&,) is becoming larger as this singularityis approached. For each solution, there is a point, [, < &,, where we find F’(&) = 0, the no-slip wall boundary condition. From this ES, De and p1 can be obtained from eqn. (41). A plot of p1 against De in Fig. (2) shows a limit point at De - 0.08. The value of &, at this limit point in Deborah number is co - 1.32; thus the limit point in De occurs at or near the limit point shown in Fig. (4) for to. This limit point also occurs at or near the limit point in A = F”(O), which is A - -0.31. Thus, to within numerical accuracy, the limit points in lo, De, and A all occur simultaneously.Kim-E [5,6] using finite difference, found a limit point at a. somewhat different Deborah number, De - 0.073; see Fig. (2). The differencebetween the value obtained by Kim-E and that reported here apparently arises because the discretization in Kim-E’s finite difference technique implicitly set stress boundary conditions at the tube wall that are different from those that are set by the flow outside the tube [6]. Figure (5) shows f(r) and -f’(r)/ r, which is proportional to the axial velocity, both computed by RK integration for &, = 0.8 on the upper
366
I(-
N
-205
0.0
I
0.2
I 0.4
I I
0.6
I
I
0.6 RAolus,
1.0
I
12
I
1.4
6 1.1
r
Fig. 5. (a) Velocity simihuity function, - f(r), and (b) axial velocity at z = 1 for De = 0.031 on the upper solution branch. The dashed lines are the correspondingNewtonian results.
branch. As Kim-E found, the velocity inside the tube is nearly Newtonian. Just outside the tube, however, the velocity becomes highly non-Newtonian. Kim-E observed steep boundary layers in the stress similarity functions, especially h, near the tube wall [5,6]. Figure 6(a) shows that the steep boundary layer in h just inside the tube wall is dwarfed by the growth of h between the tube wall and the separating cylinder. The growth appears steep even on a plot of log h versus r, Fig. 6(b). The existence of two injection solutions for values of -A less than 0.31 raises a question concerning the perturbation expansion about t = 0, the first few terms of which are given in eqn. (45). The series have only one parameter, A. Once it is given, all coefficients of eqn. (45) can be obtained through recursion relations. But for each A, there are two solutions. To which solution, then, does the perturbation expansion converge? Close inspection of the two numerical solutions as I approaches zero shows that they approach each other, but neverthelessare not identical until 6 actually reaches zero. If both solutions are analytic near 6 = 0, there can be only one conclusion: the perturbation expansion about t = 0 must have a zero radius of convergence! In fact, a ratio test performed on successive coefficients of solution
367
I 0.0
I 0.2
I
I
0.4
0.6 RADIUS.
0.6
I
1.0
I!2
I: 4
r
Fig. 6. 7zz normal force hilarity function, h, as a function of radial position, r, on (a) a linear plot, and (b) a semi-log plot, showing the steep boundary layer near the tube wall and between the tube wall and separating cylinder.
the perturbation expansion showed a radius of convergence that tends towards zero as higher order terms are ratioed. Thus, for 6 not equal to zero, the perturbation expansion cannot exceed a limiting accuracy no matter how many terms are included. Likewise, the instability of the numerical integration starting from < = 0 should not be too surprising,since for each value of -A < 0.31 chosen, there are at least two solutions to which the integration could conceivably converge. For each solution to the resealed problem, there is not only a no-slip solution, but also a slip solution for each b < 1. One need only search the no-slip solution for a different stopping condition, the one corresponding to eqn. (40), to find a solution corresponding to slip, that is to j3 < 1. As one decreases I, from the value at which F’(&) equals zero to { = 0; the slip coefficient, j3, varies from unity to zero. Thus each no-slip solution has a set of slip solutionsembeddedin it; the smaller value of the slip coefficient, the smaller, in resealed units, is the diameter of the embedded tube with the slip wall condition. The existence, for each value of -A < 0.31, of two no-slip solutions that remain distinct for t > 0, therefore implies that for any /I < 1,
368 there are again two solutions for each value of -A less than -A,( 8). A,,(p) is the limit point for /3 < 1. In addition to these two slip solutions, Menon et al. [6] found a fully analytic solution for any A. The velocity solution of this solution is a pure homogeneous alongation and so cannot satisfy the no-slip wall boundary conditions but can satisfy this boundary condition when B -C1. The stresses of this analytic slip solution contain eigensolutionsthe coefficients of which are set by the slip coefficient and other wall stress boundary conditions. 5.6 Asymptotic boundary layer analysis Near the separating cylinder, 1, H, and K are all large, and one can drop the constant, p * = -1, and replace I- 1 by I in eqns. (38a)-(38d). The result is K’=
-H--K/&
[F;“-F’/[]I=
-KS+F’K-F(K’+K/E),
H’F = 4K( F’/t - F”) - H& I’F = 2I( F’ - F/e) - I& which has the solution, F=4x
1 eW-loz) - I (
ieZAcl* I = - Ce2xt’, 2t2A
(W
369
0
0.2
0.4
SQUARE
0.6
OF RADIAL
1.0
0.6 COORDINATE,
1
<*
Fig. 7. Semi-logarithmic plot of the similarity functions H, Z(*, K.$, and F’/t against [*, the square of the resealed radial coordinate, for .f,, = 1 on the upper solution branch. The straight-line portions of the lines agree with the asymptotic boundary layer formulas given by eqn. (47).
beyond the separating cylinder. Figure (8) shows F and the stress similarity functions, I, H, and K as functions of 6, for t > &, for &, = 0.8, De = 0.031, on the upper solution branch. The asymptotic formulas, eqn. (47), are almost exact in this region. The super-steepness of the boundary layers is evident: at t = 1.3 the renomalized stresses are all in excess of lo”!
I
0.6
0.9
I
1.0 RENORMALIZED
I I.2
I
I.1 RADIUS,
I.3
c
Fig. 8. Resealed stress and velocity similarity functions beyond the separating cylinder at De = 0.031 on the upper solution branch.
370 Codusions In flows of memory fluids, the upstream boundary of the domain can require special treatment. This is obvious when the memory effects are incorporated in an integral representation, for then a pre-history of the fluid up to the time it crosses the inflow boundary must be provided, but it can also be true for a differential equation, such as the upper-convected Maxwell equation. In the porous tube problem, with Newtonian kinematics specified, the Lodge integral, with pre-history supplied by analytically continuing the flow inside the domain, has simple poles for injection and suction. The differential upper-convected Maxwell equation with Newtonian kinematics and stress boundary conditions only on the centerline agrees with the Lodge integral for suction, but for injection has the Lodge solution plus eigensolutions with undetermined coefficients that are steep at the wall and become singular there as De goes to zero. These eigensolutions and their derivatives are singular at a cylinder outside the tube wall where flow separation occurs and the radial velocity is zero. Thus if one imposes boundary conditions, namely boundedness of stress derivatives, at this upstream separating cylinder, the eigensolutions are eliminated. Thus in either suction or injection, with Newtonian kinematics, eigensolutions are only avoided by applying boundary conditions at the upstream end of the domain. When non-Newtonian kinematics are incorporated by solving the constitutive equation and the momentum balance equation simultaneously, boundary conditions must again be specified at the upstream end of the domain, and, for stability, numerical integration of the equations must proceed in the direction of flow. This integration shows a limit point at De = 0.08. On the upper solution branch, the kinematics are almost Newtonian inside the tube, but are highly non-Newtonian between the tube wail and the separating cylinder. The stresses, which become steep as the tube wall is approached from inside the tube, become even steeper outside the wall, and steeper still out beyond the separating cylinder. The steepness of these stress boundary layers is produced by the biaxial flow outside the tube and the existence of a singularity in the biaxial extensional viscosity for the UCM fluid. An asymptotic solution is derived which is valid for the highly non-Newtonian kinematics and steep stress boundary layers near the separating cylinder and outside of it. On the upper solution branch this boundary layer penetrates into the porous tube. Acknowledgement I gratefully acknowledge stimulating and encouraging discussions with, and insight given by Miral Kim-E, John Brady, Bob Brown, and Eric Shaqfeh.
371 References 1 P.W. Ye&, M.E. Kim-E, RC. Armstrong and R.A. Brown, J. Non-Newtonian Fluid Mech., 16 (1984) 173. 2 M.J. Crochet and K. Walters, Ann. Rev. Fluid Mech., 15 (1983) 241. 3 M.A. Mend&on, P.-W. Yeh, RA. Brown and RC. Armstrong, J. Non-Newtonian Fluid Me&., 10 (1982) 31. 4 J.F. Brady and A. Acrivos, J. Fluid Me&., 112 (1981) 127. 5 M. Kim-E, Ph.D. thesis, Mass. Inst. of Technol., 1984. 6 R.K. Mcnon, M. Kim-E, RC. Armstrong, R.A. Brown and J.F. Brady, J. Non-Newtonian Fluid Mech., 27 (1988) 265. 7 RG. Larson, J. Fluid. Mech., accepted (1988). 8 M.M. Dem, C.J.S. Petrie and P. Avenas, AIChE J., 21 (1975) 791.