- Email: [email protected]
The flow of dilute polymer solutions in confined geometries: a consistent numerical approach
Journal of Non-Newtonian Fluid Mechanics, 25 (1987) 347-364 Elsevier Science Publishers B.V., Amsterdam - Printed in The Netherlands
347
THE FLOW OF DILUTE POLYMER SOLUTIONS IN CONFINED GEOMETRIES: A CONSISTENT NUMERICAL APPROACH
P. BILLER
and F. PETRUCCIONE
Fakultiit ftir Physik der Universitiit Freiburg, Hermann-Herder-StraOe Breisgau (West Germany)
3, D-7800 Freiburg im
(Received January 12, 1987)
Starting from rigorous expressions derived from phase space kinetic theory for dumbbell models of polymer solutions, a new numerical approach is presented. It enables one to solve the Langevin equations governing the motion of the dumbbells in a confined geometry consistently with the momentum balance equation. As an example, we discuss the flow of a polymer solution between two parallel shearing planes. For this purpose, we consider linear and nonlinear dumbbell models and investigate typical phenomena such as, for example, the slip effect.
1, Introduction In the flow of dilute polymer solutions in confined geometries, wall effects are known to play an important role [l]. They lead to new phenomena which are expected to be observable when the size of the macromolecules cannot be ignored compared with the dimensions of the flow domain. In the case of, for example, two parallel planes which may move in opposite directions (shear flow) one finds that the concentration of macromolecules decreases towards the walls. The result is a reduction of the viscosity in the wall region. The decrease of the viscosity in this region is also due to the fact that near a wall the macromolecules are more likely in a position parallel to the wall than in the interior of the fluid domain. An experimenter can see this effect (“slip phenomenon”) when the fluid is slipping macroscopically relative to the boundaries. From the balance equation for the momentum 0377-0257/87/$03.50
0 1987 Elsevier Science Publishers B.V.
348 one can then deduce a deviation from the linear velocity profile which would exist for the pure solvent alone. Comparisons of existing studies on the shear flow situation between confining planes reveal that one must be very careful to determine the physical quantities one is interested in and not to take over simple formulae which are known to be correct in unbounded systems. The only studies in which phase space kinetic theory expressions are used rigorously seem to be those of Brunn and Grisafi [2] and Bnmn [3]. In both works, the drawbacks of the wrong definitions on the physical quantities especially for the polymer concentration are explicitly explained. Modelling the macromolecules as dumbbells of two identical beads connected by a Hookean spring, the paper of Brunn and Grisafi [2] treats the equilibrium (no flow) situation which is shown to be analytically solvable and gives the expressions for the local polymer concentration and the stress tensor in dependence on the distance to the wall and for arbitrary channel width. In a subsequent paper [3], the case of an arbitrary flow strength is investigated and analytical results are presented which are built on the assumption that the flow field remains linear but with a changed shear rate due to the wall induced slip effect. In this paper, we seek to determine precisely the flow field which can only be done by a numerical simulation of the process. Numerical studies of the shear flow problem have already been carried out by Goh et al. [4] but they did not only take the wrong definition for the polymer concentration but also they did not take into account the feedback of the polymers on the flow field at the same time moving the dumbbells in the undisturbed shear flow field. We get the exact flow field by fulfilling also the momentum balance equation. This equation determines the flow field if the stress tensor is given, but the stress tensor itself depends on the flow field. The idea is therefore to determine the flow field and the stress tensor consistently. As an example, we will treat the simple geometry of a polymer solution between two parallel shearing planes. In the following section, we will present the mathematical formulation of the problem and the basic formulae. In Section 3, the translation into a computer program is sketched and some general remarks about the simulation technique are given. The subsequent section gives our results for the flow field in the channel and the slip velocity in dependence on the relative size of the macromolecules to the channel width and on the applied shear rate, both for Hookean dumbbells and for a dumbbell with a nonlinear force law where an additional small cubic term is present. The effect of the variation of this cubic term is studied and results for the local polymer concentration and the material functions are given. The final section gives a brief summary.
349 2. Theory
We consider a dilute solution of N macromolecules confined to the volume V of a channel consisting of two.parallel planes at z = + d (Fig. 1). The macromolecules are modelled by dumbbells consisting of two identical beads each with mass m and friction coefficient [ and connected by a spring. Let r, and r, be the position vectors of the two beads and R- r, - r2 their relative vector. The internal force law of the spring is denoted by F(r). Each bead of a dumbbell experiences the force from the other bead, a stochastic force 9 which models the influence of the solvent, and a hydrodynamic drag force which in a flow field o(r) is assumed to be of a Stokes type, i.e. proportional to the relative velocity of the bead and the solvent flow field at the position of the bead. Therefore the equations of motion are mi;, = [[u(q)
- ii] - F(r, - r2) + VI,
(1)
mi;, = {[z)(q)
- i2] - F(r2 - rl) + q,.
(2)
We thus suppose that the solution is sufficiently dilute so that no interaction between different dumbbells has to be taken into account. Then the system is characterized by the single phase space distribution function $~(q, r,, pl, p2, t) which gives the probability density that at time t bead i (i = 1, 2). has a position 5 and a momentum pi. It is convenient to assume that $ can be factorized into a spatial part \c/(rl, r,, t) and a momentum part and that the latter one is identical with the explicitly known equilibrium momentum distribution function (“equilibration in momentum space” [5]). Thus, the system is fully understood when the spatial distribution function 4 is known. Consistent with the above assumption, we neglect the very small acceleration terms in eqns. (1) and (2) and obtain the two Langevin equations i+l =2)(q)
- +F(r,
(3)
- r2) + iv,,
(4)
i2 = z)(r2) - +F(r2 - rl) + +q2.
z=+d
T P
.,..i
________________---
.’
Y x__________________------
z=-d
Fig. 1. Geometry of the problem.
350 The stochastic forces are both modelled by a Gaussian white noise with zero mean and the two time correlations (77,(+?/&‘))
= 2kUS,,G(t
- My,
P =x,
Y, z),
(5)
where k is Boltzmann’s constant and T is the temperature. Because no bead may diffuse through the walls, we have the additional condition ]zi] Id,
(i=l,2).
(6) The Langevin equations (3) and (4) can easily be transformed into a Fokker-Planck equation for the spatial distribution function +( r,, r,, t)
The impenetrability condition for a purely repulsive wall states that the bead flux has no z-component for‘z; = f d, and for a wall which does not move in the z direction it takes on the form [
1
=0
kTg+F,(r;-r,)$ I
forzi=
+d
(i, j=l,
2; i#j).
(8)
The Fokker-Planck equation (7) together with the boundary condition (8) is not expected to be analytically solvable for a nonvanishing flow field u(r). We therefore determine the spatial distribution function \c/ by integrating the stochastic Langevin equations (3) and (4) numerically for a large number of realizations and averaging over these realizations. It is important to note that as soon as I/J is known all functions of interest such as the polymer concentration or the stress tensor can be evaluated as certain moments of this distribution function. But, as indicated in the introduction, one has to take the correct phase space kinetic expressions for the physical quantities. According to our model the masses are located at the positions of the beads (and not at the position of the center of mass) and so the local polymer concentration has to be defined as [5] n(r,
t) = $j-d3R[#(
r, r+R,
t)+#(r+R,
r, t)].
(9)
For identical beads this can be simplified to n(r,
t) = j-d3R$+,
r+R,
t).
00)
The normalization condition for n (and thus for J/) is +,jd3r
n(r,
t) = N/V.
(10
351 The second quantity of interest is the stress tensor 7 which is made up of a contribution 7s of the solvent and a polymer contribution (I, which itself splits up into a kinetic part uk and an intramolecular or spring part (I’ 7 = ,Ts+ ok + UC.
(12)
With ps the solvent pressure, qs the solvent viscosity and qua = 6u,/6x, + SuB/Sxa, i.e. twice the rate of strain tensor, the solvent contribution to the stress tensor is 7s = Psi - II,*-
(13)
The kinetic part of the polymer contribution to the stress tensor stems from the motion of the beads. Since we assumed “equilibration in momentum space” this term is simply [5] uk = 2nkTl.
(W
Having defined the mass density as in eqn. (9) or (10) one consequently has to define the intramolecular contribution to the polymer stress tensor as [5] u’(r,
t)=
-Jd3RJ1dsF(R)R~(r-sR, 0
T+ (1 -s)R,
t).
(15)
So far we did not specify the flow field u(r). In this paper we will treat the situation where the upper plane at z = d is moving with a constant velocity U in the +x direction and the lower plane at z = -d is moving with the same velocity U in the -x direction. We then look for a stationary solution of the form u, = u(z), u,, = u, = 0 of the momentum balance equation which in this context simply reads V l7=0.
(16)
This equation requires that the xz-component as well as the zz-component of the total stress tensor T = T(Z) have to be constant over the whole charmel. If we had only the solvent between the planes then +T= rs and the solution of eqn. (16) would be u(z) = (U/d) z, i.e. simple shear flow with a constant shear rate T0 = U/d. The effect of the macromolecules is to alter the flow field because due to the wall the function u,“, becomes z-dependent. Thus a variation of the local shear rate results as can be seen from the x-component of eqn. (16) d dz 4
-”
d4z) ~ dz
+ u,“,(z)
1
= 0.
(17)
Following Bnmn [6] one can integrate eqn. (17) two times and introduce the slip velocity to obtain y(O)d+u,=
U,
(18)
352
Fig. 2. Definition of the slip velocity. The continuous line is the exact velocity profile. The dotted line represents the imposed shear rate, while the chain-dotted one is the extrapolation of the real velocity profile in the middle of of the channel.
with the slip velocity
(1% where p(O) and u,“,(O) are the corresponding quantities evaluated in the middle of the channel at z = 0. Note that our definition of the slip velocity differs from that of Brurm [6], because he evaluates the above expression with a,“,(O) determined in an infinitely large channel. Consequently, for large channels there is no difference between the two definitions, but our definition has the advantage that it has the intuitive meaning sketched in Fig. 2 for arbitrary channel widths (see also eqn. (18)). For a purely repulsive wall u&(z), which comes out to be negative, increases towards the wall and a positive slip velocity results. Thus the local shear rate p(O) in the middle of the channel is less than the expected value T,, for a pure solvent. It is interesting to see that even if the wall region where u&(z) is mostly influenced is very small, its effects can be felt in the whole fluid because the velocity gradient is changed everywhere. We have to determine the flow field so that it fulfills eqn. (17) but u,“, depends on the flow field through eqns. (15), (7), (3) and (4) in a nontrivial way. The determination of the velocity field has thus to be self-consistent which means that the function u,“l which is generated by the motion of the dumbbells in the unknown velocity field has to fulfil the relation eqn. (17) with this velocity field. 3. The simulation The Brownian dynamic simulation of the flow of a dilute polymer solution between two plates with the consistent determination of the velocity profile requires a great amount of computational effort.
353
The simulations were performed on the vector computer Cray-l/M and the program had to be written in a vectorizable way to reduce computational time as far as possible. Usually [4], one considers the trajectory of a single dumbbell over a long period of time and evaluates the physical quantities of interest by the ergodic assumption as time averages. We do the same but with an ensemble of trajectories (it4 = lOOO),which evolve in time parallel. At the beginning of the simulation the dumbbells are randomly distributed between the two planes, the distance between the two beads being the equilibrium mean spring extension. The initial velocity profile is the linear one that would exist for the pure solvent. Then the dumbbells evolve in time until a stationary state is achieved. During further evolution, all contributions to the interesting physical quantities are determined each time step for each trajectory and at the end one averages over the time and the ensemble. The main part of the program is the integration of the equations of motion. The Langevin equations that determine our system are first-order stochastic differential equations, that have to be solved according to the boundary conditions at the walls. The interaction with the wall is assumed to be a hard-core repulsion and can therefore simply be simulated by reflective boundary conditions at z = +d. The integration algorithm consequently has to be very accurate at the boundaries. Knowing the exact analytical results in the equilibrium (no flow) case [2] we could test different algorithms to choose the most appropriate for our purpose. Let us consider generally the Ito stochastic differential equation ds, = a(~,) dt + b dW,,
(20)
where a is the drift term, b a constant and N: the standard Wiener process [71. The simplest numerical integration algorithm of the above equation is the stochastic analogue of the first order Euler method for ordinary differential equations Xt”+l= xtn + Q(xJ At + bAW&
(21)
where At = tn+l - t, denotes the time increment and AK, = Wt,+, - Wt. the corresponding increment of the Wiener process which is of the order ( At)“2. A second-order Runge-Kutta type integration algorithm is the so-called Heun method [S], which is defined by the following scheme
(22) with the Euler predictor
Xl+,= xtm+ u(xJAt
+ bAe/
(23)
354
x 0.75
k,=o.2
X %4
z?
X
-II
i
eb 0.50
0.25
I
O.ool 0.0
0.1
0.2
0.3
0.4
0.5 4
0.6
0.7
0.8
0.9
1.0
L
Fig. 3. Comparison between two different integration polymer stress tensor versus i in dimensionless notation. exact analytical curve [2], the triangles represent the curve the crosses the curve obtained by the Euler method. Here, size of the marker has the magnitude of the error bars.
algorithms: xx-component of The continous line represents obtained by the Heun method and in all the following figures,
the the and the
The application of both algorithms to our problem revealed that major difficulties arise in the optimal choice of the time step independently of the algorithm used. Significant differences between the two algorithms were detected for large channels. Figure 3 shows the polymer contribution to the xx-component of the stress tensor for a situation where the channel is ten times larger than the root-mean-square equilibrium length of the Hookean dumbbell. The comparison with the exact analytical results [2] shows that the Heun algorithm is more exact than the simple Euler algorithm. We therefore decided to choose the Heun algorithm for all situations. Once the dumbbells have reached a stationary state after a sufficiently long iteration in time-a good check is the constancy of the zz-component of the stress tensor along the whole channel width (see the z-component of eqn. (16))-one can measure the averages of interest in dependence on the distance from the wall. To do so we divide the width of the channel into a large number of cells. For each of the A4 dumbbells, we now have to distribute the contributions to the measured quantities, according to their definitions, into the corresponding cells after each time step. The averaged normalized quantity is then obtained by dividing by the number of trajectories and the number of time steps. Knowing how to compute the stress tensor we can now proceed in the description of the algorithm for the step-by-step improvement of the velocity profile. Having evaluated the xz-component of the stress tensor for the
355 velocity gradient imposed at the beginning (which is not yet constant) we begin to improve the velocity field. In accordance with eqn. (17) the total xz-component of the stress tensor has to be constant over the whole channel. This constant is evaluated in the middle of the channel at z = 0 with the help of eqns. (18) and (19) and the new velocity gradient is now evaluated in each cell by requiring eqn. (17) to be obeyed. The new velocity profile is then obtained by the trivial integration of this gradient. With this new velocity profile we again make an adequate time development of the M trajectories to obtain a corrected xz-component of the stress tensor (which now is more constant than the one before). The above procedure is then repeated until a reasonable convergence condition is satisfied. Finally all material functions of interest can be evaluated. For the numerical simulation, it is useful to introduce dimensionless quantities. To this end, we introduce a length scale d, i.e., half the channel width, a time scale [d*/kT, i.e. the mean time for a dumbbell to diffuse over the channel, and an energy scale kT. All lengths, times, and energies are measured in these units. The dimensionless quantities are denoted with “-“, so for example we have i = z/d. Remembering that N/V is the (global) polymer concentration we measure elements of the stress tensor in units of (N/V)kT. The above scales also fix the units in which all other quantities are measured. So the unit for a viscosity is (N/V){d 2 and for a velocity component it is kT/[d. In our simulation, we prefer to keep the channel width constant and to vary the ratio of the length of the dumbbell to the width of the channel by varying the spring force. For a Hookean dumbbell with the internal force law F(R) = HR with the spring constant H the relevant parameter is the dimensionless quantity fi = Hd2/kT. To determine the interesting values of fi we can use the fact that for Rouse dumbbells there exists an analytical expression which relates the spring constant to the root meanAsquare equilibrium length R, of the dumbbell H = 3kT/Ri. So we get H = 3/& where R, is the dimensionless rootmean-square length of the dumbbell in equilibrium. The A usual situation where the wall effects can be neglected is R, -=K1. For R, x=- 1 the macromolecule is much larger than the channel and we expect that our dilute dumbbell model becomes meaningless. The interesting region for the wall effect studies is thus 10-l ,( &, ,( 1, i.e. 1 s H < 100. The parameter t, which characterizes the solvent is chosen to be 0.1 in all simulations. 4. Results 4. I. Hookean dumbbells The case of Hookean dumbbells in a shear flow between two parallel planes has been thoroughly discussed by Brunn [3]. Nevertheless, his analyti-
356 cal approach is based on the assumption that the velocity gradient is constant in the whole space between the planes, its absolute value being fixed by the exact gradient in the middle of the channel, which may differ from the imposed shear rate. On the other hand, the numerical approach we developed takes into account the influences of the polymers on the velocity field, in which these polymers flow. With the help of the momentum balance equation this velocity field is evaluated self-consistently. Some velocity profiles can be seen in Fig. 4 for difzerent values of the dimensionless root-mean-square equilibrium length R, of the Hookean dumbbell. As expected the flow field is very nearly linear* if the dumbbells are small relative to the channel width. For greater R,, i.e., longer dumbbells or smaller channels, respectively, the deviation from linearity is a considerable effect and great slip velocities may arise. The dependency of the slip velocity on the channel width is shown in Fig. 5 for a constant imposed shear rate. One expects intuitively the slip velocity to increase with I&,. Figure 5 shows that this is indeed the case until the ratio R,/d is of the order of magnitude of about 1. For greater ratios the opposite behaviour is detected. Here, the very large dumbbells are inhibited to take on all orientations and therefore the influence of the walls on the stress tensor is felt not only in a small region near the planes but also in the middle of the channel. Especially the absolute value of the component uiz decreases everywhere in the channel. As this is the quantity which enters in our definition of the slip velocity (eqn. (19)) the effect occurring in Fig. 5 is no longer surprising.
Fig. 4. Velocity profiles versus i for different length ratios I?, of the linear dumbbell for a typical imposed shear rate.
357 0.40 0.35 0.30 0.25 1
0.5
0.0
1.0
j.5
2.0
2.5
3.0
R,
Fig. 5. Slip velocity versus length ratio 8, shear rate in dimensionless notation.
of the Hookean dumbbell for a typical imposed
In Fig. 6, we examined the dependency of the slip velocity on the applied shear rate for different values &, and observed that it is always linear within the limits of numerical errors. Until now we discussed only the velocity field but of course the polymer concentration and stress tensor are determined exactly too. For shear flows the components of the stress tensor that do not vanish are the diagonal elements and the xz- and zx-component, which are equal. We find that a
3.0
*=J= 20 + .
1.0
0.0
0.0
10.0
20.0
30.0
40.0
SC.0
60.0
70.0
80.0
90.0
100.0
U
Fig. 6. Slip velocity versus imposed shear rate for three different values of l?, of the linear dumbbell in dimensionless notation.
358
Fig. 7. Polymer viscosity versus 2 for the linear dumbbell model for two different values of I?, and two different imposed shear rates in dimensionless notation.
nonzero flow field leads to no measurable effects compared to the equilibrium results as far as the function n(z), c$,( z) and u,“,(z) are concerned. Therefore we refer to Ref. 2 where exact analytical expressions for these quantities in the no-flow condition are given. Clearly u,C,(z) and u,“,(z) change compared to the equilibrium situation when a shear rate is applied. To study these variations we define the functions n(z) and Nr( z) through e,C(z> = -q(z)?(z) and
(24)
u,c,(z; U= 0) - u,‘,(z) = N,(z)?*(z).
(2% Thus q(z) is nothing else but the polymer viscosity and A$( z) is the part of the first normal stress function which is induced by the flow. Note that both quantities are z-dependent and are thus defined locally. We find that the so-defined functions are independent of the applied shear rate. This situation is depicted in Figs. 7 and 8 each for two different values of the ratio R ,/d. The linearity of t,he slip velocity with i? and the fact that neither n(z) nor Nr( z) depend on U are direct consequences of the simple linear Hookean dumbbell model we used. Therefore it seems desirable to study the effects of a more realistic model, where the spring between the dumbbells is nonlinear. 4.2. Nonlinear dumbbells The nonlinear internal force law we want to investigate F(R) = HR + gR*R (H, g = const.)
(26)
359 0.125
0.100
0.075
0.050
B 0 B 0 0
0.025
0.000
0.0
0.1
0.2
0.3
0.4
0.5
D
0.6
cl
0.7
0.8
0.9
1.0
1
Fig. 8. Flow induced the same parameters
first normal stress difference as in Fig. 7 in dimensionless
versus i for the linear dumbbell notation.
models for
is the simplest modification of the Hookean force law. As in the FENE-model [5], this force law exhibits a growing stiffness with increasing elongation. For not too high shear or elongational rates this model, that can be discussed in unbounded systems to a great extent analytically, behaves as the FENEmodel [9-111. As the flows we are going to investigate are not strong the chosen model is expected to lead to reasonable results. The additional cubic term enters in the Langevin equations (3) and (4) (and in the corresponding Fokker-Planck equation (7) and the boundary condition (8)) and also in eqn. (15) which determines the intramolecular part of the polymer stress tensor. Before going on to the analysis of the results obtained for the nonlinear dumbbell model a short remark has to be made. In the nonlinear dumbbell model we vary the ratio of the size of the macromolecule to the size of the channel by a variation of A as in the linear model. But now the interpretation in a direct length scale ratio fi, is no longer possible because for the nonlinear dumbbell there exists no explicit analytical relation between the root mean square equilibrium length of the dumbbell and the parameters H and g. As the magnitude of the cubic term is typically of the order of a few percent relative_ to the linear term in weak flows the interesting region for the values of H remains the same as in the Hookean case. Realistic values for the dimensionless parameter g = gd4/kT are obtained if one chooses g^ in the order of magnitude A*/100 [lo]. The slip velocity contains the main information about the velocity field. Let us first consider the dependency of the slip velocity on the strength of
360 4.QA
A
‘I X0-
2.0-
0.0
; 0.0
I
2.5
I
5.0
I
7.5
I
10.0
!
12.5
I
15.0
3ooC& Fig. 9. Slip velocity versus strength of the nonlinearity for two different constant and a typical imposed shear rate in dimensionless notation.
values of the spring
the nonlinearity S for a constant I-? and a constant imposed shear rate. The effect of the nonlinearity is a growing stiffness with increasing elongation and thus the dumbbell is shorter than it would be in the Hookean case especially for high shear rates. Therefore the absolute value of the xz-component of the polymer contribution to the stress tensor u,“, becomes smaller for an increasing g and one expects thf slip velocity to decrease with 2. In Fig. 9 the behaviour for two typical H values is shown. We find that for small values of fi, i.e. a large dumbbell or a small charmel respectively, there has to be an additional effect which leads to an increase of the slip velocity with 8. The interpretation is that for these situations the shortening of the dumbbells leads to a greater absolute value of u,“, because the dumbbell can take on more transversal orientations. Hence, for long dumbbells relative to the channel also an increase of the slip velocity with 2 is possible. As soon as the nonlinearity is strong enough so that the dumbbell gets much shorter than the channel this effect disappears and the expected decrease of the slip velocity is detected. The dependenceAof the slip velocity on the imposed shear rate for constant values of H and 2 is shown in Fig. 10. One immediately notices the difference to the Hookean model results (Fig. 6). The dependence of the slip velocity on the imposed shear rate is no longer linear but the slip velocity increases more slowly than linearly. Turning to the stress tensor it is known [5] that the main effect of a nonlinearity in the internal force law of the dumbbell is to produce a dependence of the material functions on the applied shear rate. The be-
361 7.0
6.0
5.0
11
+ + i=o.75
4.0'3.
+ +
LO-
I3
+ 2.0-
A
+
1.0-
0.0 0.0
+ a3
A
0 I
q
b 0
25.0
A
A q
I
q
I
50.0
75.0
ci=3
A q
0
tb8.303
I
1op.o
I
125.0
I
150.0
I
175.0
200.0
U
Fig. 10. Slip velocity versus imposed shear rate for three different spring constants in dimensionless notation. The 2 values were chosen to be equal to a2/60.
haviour of the polymer viscosity q(z) for different values of the imposed shear rate is shown in Fig. 11. In the middle of the channel the viscosity decreases with 6 while near the wall curves for different 0 converge. Thus the reduction of the viscosity near the wall does not seem to be very sensitive to the imposed shear rate. Figure 12 shows that the function N,(z) behaves in a way similar to the polymer viscosity.
0.10 A
AA
A
1
D
ox 0.04-
0.02-
A
ci=lO
x
L55
0
O=lOO
D :: 0 A
% B b
0.00 0.0
I 0.1
I 0.2
I 0.3
I 0.4
I 0.5
I 0.6
I 0.7
I 0.6
b
m I 0.9
1 1.0
t Fig. 11. Polymer viscosity versus i in the nonlinear model for three different imposed shear rates in dimensionless notation.
362
1
0.0175 0.0150
A
A
A A
h=3
A A
0.0125
kO.15 A
t
x
A
x
x
0 0
x
0
A
x
1 0.0000
A0.0
A
x
ho
x
c-55
q
6=100 I
0.2
o
0
A
I 0.1
x
0
q
0.0025
A
x
x
0
A
x
A
q DX0 x A q X
OX
A
A
"VI
0.3
!
0.4
I
I
0.5
0.6
I
I
0.7
0.6
0.9
1.0
Fig. 12. Flow induced first normal stress difference versus i in the nonlinear same parameters as in Fig. 11 in dimensionless notation.
model for the
Finally, we wish to show two typical polymer c_oncentration profiles over the channel for a small and a large value of H and properly chosen g values. As can be seen in Fig. 13 for small channels the whole region is dominated by wall effects, as the polymer concentration decreases over all the channel width. For larger channels, the polymer concentration decreases only in a smaller wall region, whereas in the middle of the channel it is nearly constant. In this region, the polymer solution behaves as under bulk conditions.
x
0.7-
0
Lo.75 )cj=o.o094
0 0
0.6 0.0
I 0.1
t 0.2
I 0.3
I 0.4
I 0.5
I 0.6
1 0.7
I 0.6
8 0.9
I 1.0
i
Fig. 13. Polymer concentration
versus i in dimensionless
notation.
363 5. Conclusions The study of the flow of dilute polymer solutions in confined geometries requires great care in the definitions of the physical quantities of interest as they might differ significantly from the usual expressions for unbounded systems. The correct expressions are derived in the frame of the phase space kinetic theory [5], as was mentioned by Brunn [2]. In this formalism, the material functions of interest have to be evaluated as well-defined ensemble averages, the ensemble being a sample of trajectories in configuration space governed by the Langevin equations (3) and (4) for our case. Our numerical approach differs from previous ones [4] as we not only solve the Langevin equations but also require the momentum balance equation to be satisfied. This means that the dumbbells flow in the real velocity field that establishes itself in the solution as a consequence of the confining geometry. Our analysis therefore does not require any assumption on the form of the velocity field as it is determined iteratively during the numerical computations. The iterative procedure is based on the interplay of the stress tensor and the velocity field which are related by the momentum balance equation: the stress tensor influences the velocity field and the velocity field enters in the equations of motion from which the stress tensor is derived. One can interpret this concept as a consideration of a kind of hydrodynamic interaction. As a first application of this new approach, we studied the properties of linear and nonlinear dumbbells in a shear flow confined between two parallel planes. First, we presented the real velocity profile between the planes and showed how the slip velocity depends on the shear rate and the ratio of the dumbbell size relative to the channel width. In the Hookean dumbbell model the slip velocity depends linearly on the applied shear rate and the polymer viscosity as well as the flow induced first normal stress difference are independent of the imposed shear rate. This led us to study also a more realistic nonlinear force law. To this end, we examined a force law where in addition to the linear term a small cubic term is present. The effect is that the polymer contribution to the stress tensor and therefore also the slip velocity depend on the shear rate in a nontrivial way. The relevant material functions were studied and typical behaviours were shown. Brownian dynamic simulations appear to be of fundamental importance in the phase space kinetic theory approach in polymer dynamics. As has been shown they are also a powerful tool in the exact handling of nontrivial flow situations, e.g. the flow in confined geometries, and of nonlinear polymer models. Where the complexity of the physical situation would require too crude approximations for an analytical investigation, they offer a possibility to obtain a deep insight in the exact solution.
364 Acknowledgement We would like to thank Professor J. Honerkamp for stimulating our interest in the field of confined geometries ,and for many fruitful discussions. The numerical computations have been performed on the CRAY l/M computer of the Stuttgart University. Both authors are grateful to the Deutsche Forschungsgemeinschaft for financial support. References 1 A.B. Metzner, Y. Cohen and C. Rangel-Nafaile, J. Non-Newtonian Fluid Mech., 5 (1979) 449-462. 2 P.O. Brunn and S. Grisafi, Chem. Eng. Commun., 36 (1985) 367-383. 3 P.O. BEDIM, J. Rheol., 29 (1985) 859-886. 4 C.J. Goh, J.D. Atkinson and N. Phan-Thien, J. Chem Phys., 82 (1985) 988-994. 5 R.B. Bird, 0. Hassager, R.C. Armstrong and C.F. Curtiss, Dynamics of Polymeric Liquids, Vol II: Kinetic Theory, Wiley, New York, 1977. 6 P. Brunn, Rheol. Acta, 15 (1976) 23-29. 7 C.W. Gardiner, Handbook of Stochastic Methods, Springer, Berlin, 1985. 8 W. Riimelin, SIAM, J. Numer. Anal., 19 (1982) 604-613. 9 J. Honerkamp and H.C. Gttinger, J. Non-Newtonian Fluid Mech., 21 (1986) 157-165. 10 J. Honerkamp and H.C. Gttinger, J. Chem. Phys., 84 (1986) 7028-7035. 11 P. Biller, F. Petruccione, J. Honerkamp and H.C. &urger, J. Non-Newtonian Fluid Mech., 22 (1987) 309-324.
347
THE FLOW OF DILUTE POLYMER SOLUTIONS IN CONFINED GEOMETRIES: A CONSISTENT NUMERICAL APPROACH
P. BILLER
and F. PETRUCCIONE
Fakultiit ftir Physik der Universitiit Freiburg, Hermann-Herder-StraOe Breisgau (West Germany)
3, D-7800 Freiburg im
(Received January 12, 1987)
Starting from rigorous expressions derived from phase space kinetic theory for dumbbell models of polymer solutions, a new numerical approach is presented. It enables one to solve the Langevin equations governing the motion of the dumbbells in a confined geometry consistently with the momentum balance equation. As an example, we discuss the flow of a polymer solution between two parallel shearing planes. For this purpose, we consider linear and nonlinear dumbbell models and investigate typical phenomena such as, for example, the slip effect.
1, Introduction In the flow of dilute polymer solutions in confined geometries, wall effects are known to play an important role [l]. They lead to new phenomena which are expected to be observable when the size of the macromolecules cannot be ignored compared with the dimensions of the flow domain. In the case of, for example, two parallel planes which may move in opposite directions (shear flow) one finds that the concentration of macromolecules decreases towards the walls. The result is a reduction of the viscosity in the wall region. The decrease of the viscosity in this region is also due to the fact that near a wall the macromolecules are more likely in a position parallel to the wall than in the interior of the fluid domain. An experimenter can see this effect (“slip phenomenon”) when the fluid is slipping macroscopically relative to the boundaries. From the balance equation for the momentum 0377-0257/87/$03.50
0 1987 Elsevier Science Publishers B.V.
348 one can then deduce a deviation from the linear velocity profile which would exist for the pure solvent alone. Comparisons of existing studies on the shear flow situation between confining planes reveal that one must be very careful to determine the physical quantities one is interested in and not to take over simple formulae which are known to be correct in unbounded systems. The only studies in which phase space kinetic theory expressions are used rigorously seem to be those of Brunn and Grisafi [2] and Bnmn [3]. In both works, the drawbacks of the wrong definitions on the physical quantities especially for the polymer concentration are explicitly explained. Modelling the macromolecules as dumbbells of two identical beads connected by a Hookean spring, the paper of Brunn and Grisafi [2] treats the equilibrium (no flow) situation which is shown to be analytically solvable and gives the expressions for the local polymer concentration and the stress tensor in dependence on the distance to the wall and for arbitrary channel width. In a subsequent paper [3], the case of an arbitrary flow strength is investigated and analytical results are presented which are built on the assumption that the flow field remains linear but with a changed shear rate due to the wall induced slip effect. In this paper, we seek to determine precisely the flow field which can only be done by a numerical simulation of the process. Numerical studies of the shear flow problem have already been carried out by Goh et al. [4] but they did not only take the wrong definition for the polymer concentration but also they did not take into account the feedback of the polymers on the flow field at the same time moving the dumbbells in the undisturbed shear flow field. We get the exact flow field by fulfilling also the momentum balance equation. This equation determines the flow field if the stress tensor is given, but the stress tensor itself depends on the flow field. The idea is therefore to determine the flow field and the stress tensor consistently. As an example, we will treat the simple geometry of a polymer solution between two parallel shearing planes. In the following section, we will present the mathematical formulation of the problem and the basic formulae. In Section 3, the translation into a computer program is sketched and some general remarks about the simulation technique are given. The subsequent section gives our results for the flow field in the channel and the slip velocity in dependence on the relative size of the macromolecules to the channel width and on the applied shear rate, both for Hookean dumbbells and for a dumbbell with a nonlinear force law where an additional small cubic term is present. The effect of the variation of this cubic term is studied and results for the local polymer concentration and the material functions are given. The final section gives a brief summary.
349 2. Theory
We consider a dilute solution of N macromolecules confined to the volume V of a channel consisting of two.parallel planes at z = + d (Fig. 1). The macromolecules are modelled by dumbbells consisting of two identical beads each with mass m and friction coefficient [ and connected by a spring. Let r, and r, be the position vectors of the two beads and R- r, - r2 their relative vector. The internal force law of the spring is denoted by F(r). Each bead of a dumbbell experiences the force from the other bead, a stochastic force 9 which models the influence of the solvent, and a hydrodynamic drag force which in a flow field o(r) is assumed to be of a Stokes type, i.e. proportional to the relative velocity of the bead and the solvent flow field at the position of the bead. Therefore the equations of motion are mi;, = [[u(q)
- ii] - F(r, - r2) + VI,
(1)
mi;, = {[z)(q)
- i2] - F(r2 - rl) + q,.
(2)
We thus suppose that the solution is sufficiently dilute so that no interaction between different dumbbells has to be taken into account. Then the system is characterized by the single phase space distribution function $~(q, r,, pl, p2, t) which gives the probability density that at time t bead i (i = 1, 2). has a position 5 and a momentum pi. It is convenient to assume that $ can be factorized into a spatial part \c/(rl, r,, t) and a momentum part and that the latter one is identical with the explicitly known equilibrium momentum distribution function (“equilibration in momentum space” [5]). Thus, the system is fully understood when the spatial distribution function 4 is known. Consistent with the above assumption, we neglect the very small acceleration terms in eqns. (1) and (2) and obtain the two Langevin equations i+l =2)(q)
- +F(r,
(3)
- r2) + iv,,
(4)
i2 = z)(r2) - +F(r2 - rl) + +q2.
z=+d
T P
.,..i
________________---
.’
Y x__________________------
z=-d
Fig. 1. Geometry of the problem.
350 The stochastic forces are both modelled by a Gaussian white noise with zero mean and the two time correlations (77,(+?/&‘))
= 2kUS,,G(t
- My,
P =x,
Y, z),
(5)
where k is Boltzmann’s constant and T is the temperature. Because no bead may diffuse through the walls, we have the additional condition ]zi] Id,
(i=l,2).
(6) The Langevin equations (3) and (4) can easily be transformed into a Fokker-Planck equation for the spatial distribution function +( r,, r,, t)
The impenetrability condition for a purely repulsive wall states that the bead flux has no z-component for‘z; = f d, and for a wall which does not move in the z direction it takes on the form [
1
=0
kTg+F,(r;-r,)$ I
forzi=
+d
(i, j=l,
2; i#j).
(8)
The Fokker-Planck equation (7) together with the boundary condition (8) is not expected to be analytically solvable for a nonvanishing flow field u(r). We therefore determine the spatial distribution function \c/ by integrating the stochastic Langevin equations (3) and (4) numerically for a large number of realizations and averaging over these realizations. It is important to note that as soon as I/J is known all functions of interest such as the polymer concentration or the stress tensor can be evaluated as certain moments of this distribution function. But, as indicated in the introduction, one has to take the correct phase space kinetic expressions for the physical quantities. According to our model the masses are located at the positions of the beads (and not at the position of the center of mass) and so the local polymer concentration has to be defined as [5] n(r,
t) = $j-d3R[#(
r, r+R,
t)+#(r+R,
r, t)].
(9)
For identical beads this can be simplified to n(r,
t) = j-d3R$+,
r+R,
t).
00)
The normalization condition for n (and thus for J/) is +,jd3r
n(r,
t) = N/V.
(10
351 The second quantity of interest is the stress tensor 7 which is made up of a contribution 7s of the solvent and a polymer contribution (I, which itself splits up into a kinetic part uk and an intramolecular or spring part (I’ 7 = ,Ts+ ok + UC.
(12)
With ps the solvent pressure, qs the solvent viscosity and qua = 6u,/6x, + SuB/Sxa, i.e. twice the rate of strain tensor, the solvent contribution to the stress tensor is 7s = Psi - II,*-
(13)
The kinetic part of the polymer contribution to the stress tensor stems from the motion of the beads. Since we assumed “equilibration in momentum space” this term is simply [5] uk = 2nkTl.
(W
Having defined the mass density as in eqn. (9) or (10) one consequently has to define the intramolecular contribution to the polymer stress tensor as [5] u’(r,
t)=
-Jd3RJ1dsF(R)R~(r-sR, 0
T+ (1 -s)R,
t).
(15)
So far we did not specify the flow field u(r). In this paper we will treat the situation where the upper plane at z = d is moving with a constant velocity U in the +x direction and the lower plane at z = -d is moving with the same velocity U in the -x direction. We then look for a stationary solution of the form u, = u(z), u,, = u, = 0 of the momentum balance equation which in this context simply reads V l7=0.
(16)
This equation requires that the xz-component as well as the zz-component of the total stress tensor T = T(Z) have to be constant over the whole charmel. If we had only the solvent between the planes then +T= rs and the solution of eqn. (16) would be u(z) = (U/d) z, i.e. simple shear flow with a constant shear rate T0 = U/d. The effect of the macromolecules is to alter the flow field because due to the wall the function u,“, becomes z-dependent. Thus a variation of the local shear rate results as can be seen from the x-component of eqn. (16) d dz 4
-”
d4z) ~ dz
+ u,“,(z)
1
= 0.
(17)
Following Bnmn [6] one can integrate eqn. (17) two times and introduce the slip velocity to obtain y(O)d+u,=
U,
(18)
352
Fig. 2. Definition of the slip velocity. The continuous line is the exact velocity profile. The dotted line represents the imposed shear rate, while the chain-dotted one is the extrapolation of the real velocity profile in the middle of of the channel.
with the slip velocity
(1% where p(O) and u,“,(O) are the corresponding quantities evaluated in the middle of the channel at z = 0. Note that our definition of the slip velocity differs from that of Brurm [6], because he evaluates the above expression with a,“,(O) determined in an infinitely large channel. Consequently, for large channels there is no difference between the two definitions, but our definition has the advantage that it has the intuitive meaning sketched in Fig. 2 for arbitrary channel widths (see also eqn. (18)). For a purely repulsive wall u&(z), which comes out to be negative, increases towards the wall and a positive slip velocity results. Thus the local shear rate p(O) in the middle of the channel is less than the expected value T,, for a pure solvent. It is interesting to see that even if the wall region where u&(z) is mostly influenced is very small, its effects can be felt in the whole fluid because the velocity gradient is changed everywhere. We have to determine the flow field so that it fulfills eqn. (17) but u,“, depends on the flow field through eqns. (15), (7), (3) and (4) in a nontrivial way. The determination of the velocity field has thus to be self-consistent which means that the function u,“l which is generated by the motion of the dumbbells in the unknown velocity field has to fulfil the relation eqn. (17) with this velocity field. 3. The simulation The Brownian dynamic simulation of the flow of a dilute polymer solution between two plates with the consistent determination of the velocity profile requires a great amount of computational effort.
353
The simulations were performed on the vector computer Cray-l/M and the program had to be written in a vectorizable way to reduce computational time as far as possible. Usually [4], one considers the trajectory of a single dumbbell over a long period of time and evaluates the physical quantities of interest by the ergodic assumption as time averages. We do the same but with an ensemble of trajectories (it4 = lOOO),which evolve in time parallel. At the beginning of the simulation the dumbbells are randomly distributed between the two planes, the distance between the two beads being the equilibrium mean spring extension. The initial velocity profile is the linear one that would exist for the pure solvent. Then the dumbbells evolve in time until a stationary state is achieved. During further evolution, all contributions to the interesting physical quantities are determined each time step for each trajectory and at the end one averages over the time and the ensemble. The main part of the program is the integration of the equations of motion. The Langevin equations that determine our system are first-order stochastic differential equations, that have to be solved according to the boundary conditions at the walls. The interaction with the wall is assumed to be a hard-core repulsion and can therefore simply be simulated by reflective boundary conditions at z = +d. The integration algorithm consequently has to be very accurate at the boundaries. Knowing the exact analytical results in the equilibrium (no flow) case [2] we could test different algorithms to choose the most appropriate for our purpose. Let us consider generally the Ito stochastic differential equation ds, = a(~,) dt + b dW,,
(20)
where a is the drift term, b a constant and N: the standard Wiener process [71. The simplest numerical integration algorithm of the above equation is the stochastic analogue of the first order Euler method for ordinary differential equations Xt”+l= xtn + Q(xJ At + bAW&
(21)
where At = tn+l - t, denotes the time increment and AK, = Wt,+, - Wt. the corresponding increment of the Wiener process which is of the order ( At)“2. A second-order Runge-Kutta type integration algorithm is the so-called Heun method [S], which is defined by the following scheme
(22) with the Euler predictor
Xl+,= xtm+ u(xJAt
+ bAe/
(23)
354
x 0.75
k,=o.2
X %4
z?
X
-II
i
eb 0.50
0.25
I
O.ool 0.0
0.1
0.2
0.3
0.4
0.5 4
0.6
0.7
0.8
0.9
1.0
L
Fig. 3. Comparison between two different integration polymer stress tensor versus i in dimensionless notation. exact analytical curve [2], the triangles represent the curve the crosses the curve obtained by the Euler method. Here, size of the marker has the magnitude of the error bars.
algorithms: xx-component of The continous line represents obtained by the Heun method and in all the following figures,
the the and the
The application of both algorithms to our problem revealed that major difficulties arise in the optimal choice of the time step independently of the algorithm used. Significant differences between the two algorithms were detected for large channels. Figure 3 shows the polymer contribution to the xx-component of the stress tensor for a situation where the channel is ten times larger than the root-mean-square equilibrium length of the Hookean dumbbell. The comparison with the exact analytical results [2] shows that the Heun algorithm is more exact than the simple Euler algorithm. We therefore decided to choose the Heun algorithm for all situations. Once the dumbbells have reached a stationary state after a sufficiently long iteration in time-a good check is the constancy of the zz-component of the stress tensor along the whole channel width (see the z-component of eqn. (16))-one can measure the averages of interest in dependence on the distance from the wall. To do so we divide the width of the channel into a large number of cells. For each of the A4 dumbbells, we now have to distribute the contributions to the measured quantities, according to their definitions, into the corresponding cells after each time step. The averaged normalized quantity is then obtained by dividing by the number of trajectories and the number of time steps. Knowing how to compute the stress tensor we can now proceed in the description of the algorithm for the step-by-step improvement of the velocity profile. Having evaluated the xz-component of the stress tensor for the
355 velocity gradient imposed at the beginning (which is not yet constant) we begin to improve the velocity field. In accordance with eqn. (17) the total xz-component of the stress tensor has to be constant over the whole channel. This constant is evaluated in the middle of the channel at z = 0 with the help of eqns. (18) and (19) and the new velocity gradient is now evaluated in each cell by requiring eqn. (17) to be obeyed. The new velocity profile is then obtained by the trivial integration of this gradient. With this new velocity profile we again make an adequate time development of the M trajectories to obtain a corrected xz-component of the stress tensor (which now is more constant than the one before). The above procedure is then repeated until a reasonable convergence condition is satisfied. Finally all material functions of interest can be evaluated. For the numerical simulation, it is useful to introduce dimensionless quantities. To this end, we introduce a length scale d, i.e., half the channel width, a time scale [d*/kT, i.e. the mean time for a dumbbell to diffuse over the channel, and an energy scale kT. All lengths, times, and energies are measured in these units. The dimensionless quantities are denoted with “-“, so for example we have i = z/d. Remembering that N/V is the (global) polymer concentration we measure elements of the stress tensor in units of (N/V)kT. The above scales also fix the units in which all other quantities are measured. So the unit for a viscosity is (N/V){d 2 and for a velocity component it is kT/[d. In our simulation, we prefer to keep the channel width constant and to vary the ratio of the length of the dumbbell to the width of the channel by varying the spring force. For a Hookean dumbbell with the internal force law F(R) = HR with the spring constant H the relevant parameter is the dimensionless quantity fi = Hd2/kT. To determine the interesting values of fi we can use the fact that for Rouse dumbbells there exists an analytical expression which relates the spring constant to the root meanAsquare equilibrium length R, of the dumbbell H = 3kT/Ri. So we get H = 3/& where R, is the dimensionless rootmean-square length of the dumbbell in equilibrium. The A usual situation where the wall effects can be neglected is R, -=K1. For R, x=- 1 the macromolecule is much larger than the channel and we expect that our dilute dumbbell model becomes meaningless. The interesting region for the wall effect studies is thus 10-l ,( &, ,( 1, i.e. 1 s H < 100. The parameter t, which characterizes the solvent is chosen to be 0.1 in all simulations. 4. Results 4. I. Hookean dumbbells The case of Hookean dumbbells in a shear flow between two parallel planes has been thoroughly discussed by Brunn [3]. Nevertheless, his analyti-
356 cal approach is based on the assumption that the velocity gradient is constant in the whole space between the planes, its absolute value being fixed by the exact gradient in the middle of the channel, which may differ from the imposed shear rate. On the other hand, the numerical approach we developed takes into account the influences of the polymers on the velocity field, in which these polymers flow. With the help of the momentum balance equation this velocity field is evaluated self-consistently. Some velocity profiles can be seen in Fig. 4 for difzerent values of the dimensionless root-mean-square equilibrium length R, of the Hookean dumbbell. As expected the flow field is very nearly linear* if the dumbbells are small relative to the channel width. For greater R,, i.e., longer dumbbells or smaller channels, respectively, the deviation from linearity is a considerable effect and great slip velocities may arise. The dependency of the slip velocity on the channel width is shown in Fig. 5 for a constant imposed shear rate. One expects intuitively the slip velocity to increase with I&,. Figure 5 shows that this is indeed the case until the ratio R,/d is of the order of magnitude of about 1. For greater ratios the opposite behaviour is detected. Here, the very large dumbbells are inhibited to take on all orientations and therefore the influence of the walls on the stress tensor is felt not only in a small region near the planes but also in the middle of the channel. Especially the absolute value of the component uiz decreases everywhere in the channel. As this is the quantity which enters in our definition of the slip velocity (eqn. (19)) the effect occurring in Fig. 5 is no longer surprising.
Fig. 4. Velocity profiles versus i for different length ratios I?, of the linear dumbbell for a typical imposed shear rate.
357 0.40 0.35 0.30 0.25 1
0.5
0.0
1.0
j.5
2.0
2.5
3.0
R,
Fig. 5. Slip velocity versus length ratio 8, shear rate in dimensionless notation.
of the Hookean dumbbell for a typical imposed
In Fig. 6, we examined the dependency of the slip velocity on the applied shear rate for different values &, and observed that it is always linear within the limits of numerical errors. Until now we discussed only the velocity field but of course the polymer concentration and stress tensor are determined exactly too. For shear flows the components of the stress tensor that do not vanish are the diagonal elements and the xz- and zx-component, which are equal. We find that a
3.0
*=J= 20 + .
1.0
0.0
0.0
10.0
20.0
30.0
40.0
SC.0
60.0
70.0
80.0
90.0
100.0
U
Fig. 6. Slip velocity versus imposed shear rate for three different values of l?, of the linear dumbbell in dimensionless notation.
358
Fig. 7. Polymer viscosity versus 2 for the linear dumbbell model for two different values of I?, and two different imposed shear rates in dimensionless notation.
nonzero flow field leads to no measurable effects compared to the equilibrium results as far as the function n(z), c$,( z) and u,“,(z) are concerned. Therefore we refer to Ref. 2 where exact analytical expressions for these quantities in the no-flow condition are given. Clearly u,C,(z) and u,“,(z) change compared to the equilibrium situation when a shear rate is applied. To study these variations we define the functions n(z) and Nr( z) through e,C(z> = -q(z)?(z) and
(24)
u,c,(z; U= 0) - u,‘,(z) = N,(z)?*(z).
(2% Thus q(z) is nothing else but the polymer viscosity and A$( z) is the part of the first normal stress function which is induced by the flow. Note that both quantities are z-dependent and are thus defined locally. We find that the so-defined functions are independent of the applied shear rate. This situation is depicted in Figs. 7 and 8 each for two different values of the ratio R ,/d. The linearity of t,he slip velocity with i? and the fact that neither n(z) nor Nr( z) depend on U are direct consequences of the simple linear Hookean dumbbell model we used. Therefore it seems desirable to study the effects of a more realistic model, where the spring between the dumbbells is nonlinear. 4.2. Nonlinear dumbbells The nonlinear internal force law we want to investigate F(R) = HR + gR*R (H, g = const.)
(26)
359 0.125
0.100
0.075
0.050
B 0 B 0 0
0.025
0.000
0.0
0.1
0.2
0.3
0.4
0.5
D
0.6
cl
0.7
0.8
0.9
1.0
1
Fig. 8. Flow induced the same parameters
first normal stress difference as in Fig. 7 in dimensionless
versus i for the linear dumbbell notation.
models for
is the simplest modification of the Hookean force law. As in the FENE-model [5], this force law exhibits a growing stiffness with increasing elongation. For not too high shear or elongational rates this model, that can be discussed in unbounded systems to a great extent analytically, behaves as the FENEmodel [9-111. As the flows we are going to investigate are not strong the chosen model is expected to lead to reasonable results. The additional cubic term enters in the Langevin equations (3) and (4) (and in the corresponding Fokker-Planck equation (7) and the boundary condition (8)) and also in eqn. (15) which determines the intramolecular part of the polymer stress tensor. Before going on to the analysis of the results obtained for the nonlinear dumbbell model a short remark has to be made. In the nonlinear dumbbell model we vary the ratio of the size of the macromolecule to the size of the channel by a variation of A as in the linear model. But now the interpretation in a direct length scale ratio fi, is no longer possible because for the nonlinear dumbbell there exists no explicit analytical relation between the root mean square equilibrium length of the dumbbell and the parameters H and g. As the magnitude of the cubic term is typically of the order of a few percent relative_ to the linear term in weak flows the interesting region for the values of H remains the same as in the Hookean case. Realistic values for the dimensionless parameter g = gd4/kT are obtained if one chooses g^ in the order of magnitude A*/100 [lo]. The slip velocity contains the main information about the velocity field. Let us first consider the dependency of the slip velocity on the strength of
360 4.QA
A
‘I X0-
2.0-
0.0
; 0.0
I
2.5
I
5.0
I
7.5
I
10.0
!
12.5
I
15.0
3ooC& Fig. 9. Slip velocity versus strength of the nonlinearity for two different constant and a typical imposed shear rate in dimensionless notation.
values of the spring
the nonlinearity S for a constant I-? and a constant imposed shear rate. The effect of the nonlinearity is a growing stiffness with increasing elongation and thus the dumbbell is shorter than it would be in the Hookean case especially for high shear rates. Therefore the absolute value of the xz-component of the polymer contribution to the stress tensor u,“, becomes smaller for an increasing g and one expects thf slip velocity to decrease with 2. In Fig. 9 the behaviour for two typical H values is shown. We find that for small values of fi, i.e. a large dumbbell or a small charmel respectively, there has to be an additional effect which leads to an increase of the slip velocity with 8. The interpretation is that for these situations the shortening of the dumbbells leads to a greater absolute value of u,“, because the dumbbell can take on more transversal orientations. Hence, for long dumbbells relative to the channel also an increase of the slip velocity with 2 is possible. As soon as the nonlinearity is strong enough so that the dumbbell gets much shorter than the channel this effect disappears and the expected decrease of the slip velocity is detected. The dependenceAof the slip velocity on the imposed shear rate for constant values of H and 2 is shown in Fig. 10. One immediately notices the difference to the Hookean model results (Fig. 6). The dependence of the slip velocity on the imposed shear rate is no longer linear but the slip velocity increases more slowly than linearly. Turning to the stress tensor it is known [5] that the main effect of a nonlinearity in the internal force law of the dumbbell is to produce a dependence of the material functions on the applied shear rate. The be-
361 7.0
6.0
5.0
11
+ + i=o.75
4.0'3.
+ +
LO-
I3
+ 2.0-
A
+
1.0-
0.0 0.0
+ a3
A
0 I
q
b 0
25.0
A
A q
I
q
I
50.0
75.0
ci=3
A q
0
tb8.303
I
1op.o
I
125.0
I
150.0
I
175.0
200.0
U
Fig. 10. Slip velocity versus imposed shear rate for three different spring constants in dimensionless notation. The 2 values were chosen to be equal to a2/60.
haviour of the polymer viscosity q(z) for different values of the imposed shear rate is shown in Fig. 11. In the middle of the channel the viscosity decreases with 6 while near the wall curves for different 0 converge. Thus the reduction of the viscosity near the wall does not seem to be very sensitive to the imposed shear rate. Figure 12 shows that the function N,(z) behaves in a way similar to the polymer viscosity.
0.10 A
AA
A
1
D
ox 0.04-
0.02-
A
ci=lO
x
L55
0
O=lOO
D :: 0 A
% B b
0.00 0.0
I 0.1
I 0.2
I 0.3
I 0.4
I 0.5
I 0.6
I 0.7
I 0.6
b
m I 0.9
1 1.0
t Fig. 11. Polymer viscosity versus i in the nonlinear model for three different imposed shear rates in dimensionless notation.
362
1
0.0175 0.0150
A
A
A A
h=3
A A
0.0125
kO.15 A
t
x
A
x
x
0 0
x
0
A
x
1 0.0000
A0.0
A
x
ho
x
c-55
q
6=100 I
0.2
o
0
A
I 0.1
x
0
q
0.0025
A
x
x
0
A
x
A
q DX0 x A q X
OX
A
A
"VI
0.3
!
0.4
I
I
0.5
0.6
I
I
0.7
0.6
0.9
1.0
Fig. 12. Flow induced first normal stress difference versus i in the nonlinear same parameters as in Fig. 11 in dimensionless notation.
model for the
Finally, we wish to show two typical polymer c_oncentration profiles over the channel for a small and a large value of H and properly chosen g values. As can be seen in Fig. 13 for small channels the whole region is dominated by wall effects, as the polymer concentration decreases over all the channel width. For larger channels, the polymer concentration decreases only in a smaller wall region, whereas in the middle of the channel it is nearly constant. In this region, the polymer solution behaves as under bulk conditions.
x
0.7-
0
Lo.75 )cj=o.o094
0 0
0.6 0.0
I 0.1
t 0.2
I 0.3
I 0.4
I 0.5
I 0.6
1 0.7
I 0.6
8 0.9
I 1.0
i
Fig. 13. Polymer concentration
versus i in dimensionless
notation.
363 5. Conclusions The study of the flow of dilute polymer solutions in confined geometries requires great care in the definitions of the physical quantities of interest as they might differ significantly from the usual expressions for unbounded systems. The correct expressions are derived in the frame of the phase space kinetic theory [5], as was mentioned by Brunn [2]. In this formalism, the material functions of interest have to be evaluated as well-defined ensemble averages, the ensemble being a sample of trajectories in configuration space governed by the Langevin equations (3) and (4) for our case. Our numerical approach differs from previous ones [4] as we not only solve the Langevin equations but also require the momentum balance equation to be satisfied. This means that the dumbbells flow in the real velocity field that establishes itself in the solution as a consequence of the confining geometry. Our analysis therefore does not require any assumption on the form of the velocity field as it is determined iteratively during the numerical computations. The iterative procedure is based on the interplay of the stress tensor and the velocity field which are related by the momentum balance equation: the stress tensor influences the velocity field and the velocity field enters in the equations of motion from which the stress tensor is derived. One can interpret this concept as a consideration of a kind of hydrodynamic interaction. As a first application of this new approach, we studied the properties of linear and nonlinear dumbbells in a shear flow confined between two parallel planes. First, we presented the real velocity profile between the planes and showed how the slip velocity depends on the shear rate and the ratio of the dumbbell size relative to the channel width. In the Hookean dumbbell model the slip velocity depends linearly on the applied shear rate and the polymer viscosity as well as the flow induced first normal stress difference are independent of the imposed shear rate. This led us to study also a more realistic nonlinear force law. To this end, we examined a force law where in addition to the linear term a small cubic term is present. The effect is that the polymer contribution to the stress tensor and therefore also the slip velocity depend on the shear rate in a nontrivial way. The relevant material functions were studied and typical behaviours were shown. Brownian dynamic simulations appear to be of fundamental importance in the phase space kinetic theory approach in polymer dynamics. As has been shown they are also a powerful tool in the exact handling of nontrivial flow situations, e.g. the flow in confined geometries, and of nonlinear polymer models. Where the complexity of the physical situation would require too crude approximations for an analytical investigation, they offer a possibility to obtain a deep insight in the exact solution.
364 Acknowledgement We would like to thank Professor J. Honerkamp for stimulating our interest in the field of confined geometries ,and for many fruitful discussions. The numerical computations have been performed on the CRAY l/M computer of the Stuttgart University. Both authors are grateful to the Deutsche Forschungsgemeinschaft for financial support. References 1 A.B. Metzner, Y. Cohen and C. Rangel-Nafaile, J. Non-Newtonian Fluid Mech., 5 (1979) 449-462. 2 P.O. Brunn and S. Grisafi, Chem. Eng. Commun., 36 (1985) 367-383. 3 P.O. BEDIM, J. Rheol., 29 (1985) 859-886. 4 C.J. Goh, J.D. Atkinson and N. Phan-Thien, J. Chem Phys., 82 (1985) 988-994. 5 R.B. Bird, 0. Hassager, R.C. Armstrong and C.F. Curtiss, Dynamics of Polymeric Liquids, Vol II: Kinetic Theory, Wiley, New York, 1977. 6 P. Brunn, Rheol. Acta, 15 (1976) 23-29. 7 C.W. Gardiner, Handbook of Stochastic Methods, Springer, Berlin, 1985. 8 W. Riimelin, SIAM, J. Numer. Anal., 19 (1982) 604-613. 9 J. Honerkamp and H.C. Gttinger, J. Non-Newtonian Fluid Mech., 21 (1986) 157-165. 10 J. Honerkamp and H.C. Gttinger, J. Chem. Phys., 84 (1986) 7028-7035. 11 P. Biller, F. Petruccione, J. Honerkamp and H.C. &urger, J. Non-Newtonian Fluid Mech., 22 (1987) 309-324.