Viscoelastic flow simulations through an array of cylinders
J. J. J. GillissenDepartment of Chemical Engineering, Delft University of Technology, Julianalaan 136, 2628 BL Delft, The Netherlands
(Received 27 June 2012; revised manuscript received 14 January 2013; published 7 February 2013) Polymer solution flow is studied numerically in a periodic, hexagonal array of cylinders as a model for a porous medium. We use a lattice Boltzmann method supplemented by a polymer stress, where the polymers are modeled as finitely extensible, nonlinear, elastic dumbbells. The simulated, nonmonotonic behavior of the effective viscosity μeffas a function of the Weissenberg number We is in qualitative agreement with experiments in the literature. An analytical model, which replaces the flexible polymers by rods and that replaces the flow field in the porous medium by a superposition of shear and elongation, correctly reproduces the simulated μeffas a function of the polymer extensibility parameter b in the limit of large We.
DOI:10.1103/PhysRevE.87.023003 PACS number(s): 47.56.+r, 47.57.Ng, 47.50.Cd, 46.35.+z
I. INTRODUCTION
Darcy’s law relates the pressure gradient −dp/dx in a porous medium to the fluid velocity u:
−dp
dx = μu
k . (1)
Here k is the permeability of the porous medium, μ is the fluid dynamic viscosity, and· · · denotes a volume average over the void space. Equation(1)is valid when the flow is Newtonian and the Reynolds number Re= u√k/ν 1, where√kis the typical length scale of the pores, ν= μ/ρ is the fluid kinematic viscosity, and ρ is the fluid mass density.
Viscoelastic flows of polymer solutions in porous media involve a complex interplay of polymer chains and irregular flow patterns. For these liquids pressure versus flow rate experiments have been conducted, not only in porous media and packed beds, but also in artificially made geometries that mimic the conditions in a porous medium [1–4]. From these studies, the following picture has emerged. Below a critical flow speed Darcy’s law [Eq.(1)] holds, with μ replaced by an effective viscosity μeff, which corresponds to the shear rate-dependent viscosity of the polymer solution, as can be measured using a shear rheometer. Above the critical flow speed, the pressure rises above the effective Darcy pressure. This effect is believed to be due to the coil-stretch transition, which occurs when the Weissenberg number We= λγ 1, which corresponds to the departure of the polymers from their equilibrium, coiled configuration [2]. Here λ is the polymer relaxation time, and γ is the typical fluid shear rate.
The goal of the present work is to understand the relation between the porous geometry, the flow topology, and the pressure loss, induced by the polymers. For this purpose we use numerical simulations of a viscoelastic flow in a repeat unit of an array of cylinders, which represents a small portion of a porous medium.
In the simulations the polymer stress is computed by evolving the time- and space-dependencies of the polymer configuration. Simulations as such have previously been performed in undulated channels [5,6] and in repeat units of arrays of cylinders [7,8]. These previous simulations have been limited, however, to modest values of the Weissenberg number We 5. Consequently the coil-stretch transition has not yet been fully resolved numerically in a porous medium.
In the present work simulations have been conducted up to We∼ 103, which is sufficient to fully stretch the polymer chains.
II. GOVERNING EQUATIONS
Polymer solution flow at Re 1 is governed by the Stokes’ equation, supplemented with the divergence of the polymer stress tensorτP:
0= ∇ · (−pδ + 2μS + τP). (2) Here p is the fluid pressure,δ is the unit tensor, S = 1
2(∇u +
∇uT
) is the rate of strain tensor, and u is the fluid velocity vector.
We model the flexible polymers as finitely extensible, nonlinear elastic (FENE) dumbbells [9]. As drawn in Fig.1, the FENE dumbbell consists of two spheres of diameter d, which are separated by a vector l and connected by a FENE spring, which accounts for the tendency of Brownian fluctuations to drive the polymer to its equilibrium, coiled configuration. The time scale of the spring λ= 3π4 μdleq2/kBT = 3π4μd/H is referred to as the polymer relaxation time, with kB the Boltzmann constant, T the temperature, and H the spring constant at zero extension. The maximum and equilibrium lengths of the spring are lmax and leq. The dumbbells are assumed massless, noninteracting, and much smaller than the typical pore size√k. The dumbbell stress component τijPis the
ith component of the spring force summed over the dumbbells that cross the unit plane with a normal in the j direction:
τP = μα λ c 1−trcb − δ . (3a)
The stress depends on the configuration tensor c= ll/leq2, where· · · is an average over polymers contained in a (small compared to√k) volume centered at the point, where the stress is to be determined. The dumbbell stress is determined by three dimensionless parameters: the concentration parameter
α= 94φl2eq/d2, the extensibility parameter b= lmax2 / leq2, and the Weissenberg number We= γ λ. Here φ is the dumbbell volume fraction and γ = u/√k is the typical fluid strain rate.
The FENE-dumbbell has extensively been used to model viscoelastic flows [9]. The dumbbell is based on a single
FIG. 1. The FENE model, where l is the polymer length, which is leqin zero flow and extends to lmax=√bleqin strong flow, where
bis the extensibility parameter.
relaxation time λ. Since a real polymer possesses a broad spectrum of relaxation times, the model provides explanations at a qualitative level only. Quantitatively, the model is limited. In transient, extensional flow, for instance, the FENE dumbbell underpredicts the polymer stress by several orders of magnitude [10]. In this context the FENE chain with N beads provides improvement over the FENE dumbbell (N= 2). In the present work we qualitatively explain the pressure loss in viscoelastic, porous media flows. For this purpose we use the FENE dumbbell, since it captures the qualitative trends at reasonable computational expenses.
The time and space dependencies of model configurations, such as the FENE dumbbell or the FENE chain, have been computed in the literature using various numerical methods. We roughly classify these into field methods and particle methods, where the latter are usually more exact and computa-tionally more extensive. Particle methods compute individual polymer configurations, while field methods compute the polymer configuration statistics. The statistics are described by field variables, that evolve according to partial differential equations. Similar to the polymers, the solvent molecules can also be treated as particles or as fields. In particle-solvent methods the solvent is replaced by individual particles, which collide and interact with each other and with the polymer particles; see, e.g., Ref. [11]. In field-solvent methods the solvent is treated as a continuum that evolves according to the (Navier-)Stokes or Boltzman equation, and the polymer interacts with the solvent using empirical drag laws and stochastic forcing; see, e.g., Ref. [12]. To minimize the computational expenses in the present work, we use field representations for both the solvent as well as for the polymers. The dumbbell configuration tensor evolves according to [9]:
Dc Dt − ∇u T · c − c · ∇u = 1 λ δ − c 1−trcb + κ∇2c. (3b) Here κ is the polymer mass diffusivity due to Brownian motion. Equation(3b)describes how the polymers rotate and stretch due to the fluid velocity gradient and contract due to the spring force. For an extensive derivation of Eqs.(3a)and (3b)the reader is referred to Ref. [9].
III. POLYMER VISCOSITY
To account for the polymers we replace the viscosity μ in Darcy’s law [Eq.(1)] by an effective viscosity μeff, which contains the contributions of both the solvent and the polymers:
μeff = μ + μP:
−dp
dx =
(μ+ μP)u
k . (4)
We base the polymer viscosity μP on the integral energy equation:
−dp
dxu= , (5)
where the energy dissipation rate is due to Newtonian stress τNand due to polymer stressτP:
= N+ P =S : (τN+ τP)= S : (2μS + τP). (6) By defining the polymer viscosity as
μP = μ
P
N, (7)
the energy balance [Eq.(5)] becomes −dp dxu= N 1+μ P μ . (8)
To proceed we use that in Newtonian flow μP = 0, and Eq.(8)reduces to the Newtonian Darcy law [Eq.(1)], which is satisfied when N = μu2/k. By inserting this expression for N into Eq. (8) we recover Eq. (4). It is noted that our definitions of μP and k are based on the Newtonian value for N. Therefore, this μP also reflects changes in the Newtonian stress due to changes in the flow patterns induced by viscoelasticity.
To demonstrate some basic properties of μP, we consider the following, alternative, the so-called Giesekus expression for the polymer stress, which is obtained by combining Eqs.(3a)and(3b):
τP = μα ∇uT· c + c · ∇u − Dc Dt − κ∇ 2c . (9)
For We 1 Eq.(3b)is dominated by the first term on the right-hand side. Equating this term to zero and using that b 1 gives an isotropic configuration tensor, c= δ, which corresponds to randomly coiled polymers. Inserting c= δ into Eq.(9)gives a viscous polymer stress:τP = α2μS and P = α2μS : S. This exercise demonstrates that for We 1 the polymer viscosity [Eq. (7)] equals μP = αμ. In the rest of this paper we will scale the polymer viscosity with this zero shear value. The scaled viscosity is referred to as the intrinsic viscosity f :
f = μ P μα = P α N = S :τP 2α S : S. (10)
For illustration purposes, we present in Figs.2(a)and2(b) the intrinsic viscosity for b= 100 in steady elongational flow,
∇u = γ (δxδx− δyδy), and in steady shear flow,∇u = γ δyδx,
whereδx andδy are the unit vectors in the x and y direction. In elongational flow, the polymers extend along the axis of positive strain rateδx, resulting in an increasing f as a function
10−2 10−1 100 101 102 100 101 102 We f b 2 (a) 10−2 10−1 100 101 102 10−1 100 We f (b)
FIG. 2. The intrinsic polymer viscosity f for b= 100 as a function of the Weissenberg number We in (a) steady elongation and (b) steady shear.
of We, referred to as shear thickening. In the limit We→ ∞, the polymers are fully extended, c= bδxδx, corresponding to
the maximum value of f = b2. In shear flow the polymers rotate away from the axes of deformation, giving a decreasing
f with increasing We, referred to as shear thinning.
In this work we study the behavior of f in porous media flow. To this end we examine the numerical solution to Eqs.(2) and(3), where the flow is driven by a constant pressure gradient in a complex geometry, which mimics a small portion of a porous medium.
IV. NUMERICAL METHODS
We numerically approximate the solution to Eqs. (2) and(3). We discretize the polymer evolution equation [Eq. (3b)] using the finite volume method, with the first order upwind scheme for the convection term and the second order, central difference scheme for the diffusion term. Time integration is achieved with the explicit, second order Adams Bashforth scheme for the advection, diffusion, and polymer stretching terms and the implicit, second order Crank Nicolson scheme for the nonlinear spring term. On the solid walls we set n· ∇c = 0, with n the wall-normal unit vector, which means that there is no diffusion flux of polymers through the walls.
The solution to the Stokes equation with the polymer stress [Eq. (2)] is approximated using a lattice Boltzmann (LB) scheme [13]. The main advantage of using LB is that no-slip conditions on complexly shaped boundaries are easily and efficiently implemented using the bounce-back scheme, with the wall located halfway between a fluid node and a solid node [14]. The LB method is based on discretizing the Boltzmann equation, which governs the probability distribution function
gof the velocityv and the position x of the solvent molecules. The key of the LB method is to discretize the velocity space into a minimum set of velocities vα, which is still large enough to represent the essential features of g that play a role in the (Navier-)Stokes limit. In this work we use nine velocity directions vα on a two-dimensional (2D) lattice, referred to as D2Q9 [13]. The velocity directions vα are such that tvα equal the distances between neigh-boring lattice points, where t is the computational time step.
We use the following lattice Boltzmann scheme, which is extended with the polymer stress [15]:
gα(t+ t,x + tvα)− gα(t,x) t = g eq α (t,x)− gα(t,x) τ − ωα c2 vα− u −vαvα· u c2 · dp dxδx +ωα 2 vαvα c2 − 3δ :τP. (11a)
Here gα equals g evaluated at vα multiplied with the corresponding weight ωα of the Gauss-Hermite quadrature [16], and τ = ν/c2+ t/2 is the relaxation time, with c the sound speed, which in the present LB scheme is related to
t and the grid spacing x by c=√1 3
x
t. The left-hand side of Eq. (11a) represents the streaming of the solvent molecules, which is numerically integrated over one t by simply shifting gα between neighboring lattice points. The first term on the right-hand side (r.h.s.) of Eq. (11a) is the change of gαdue to collisions between the solvent molecules, which is modeled as a relaxation process toward the Maxwell Boltzmann distribution geq[13]: gαeq= ωαρ 1+vα· u c2 + 1 2 (vα· u)2 c4 − 1 2 |u|2 c2 , (11b)
which in LB is truncated to second order in the Mach number Ma= |u|c. The second and the third terms on the r.h.s. of Eq.(11a)are the change of gαdue to the mean pressure gradient and due to the polymer stress, respectively.
To find the steady state solution to Eqs.(2)and(3), we nu-merically integrate Eqs.(3)and(11)in time until a steady state is reached. The initial conditions are zero flow and equilibrium polymer configuration, c= δ. The numerical procedure for one time step starts with computing the velocity field from the distribution function using u=αgαvα/
αgα [13]. Then the velocity gradient∇u is computed using the second order, central differencing scheme, and this∇u is fed into the polymer configuration equation(3b), which is then integrated over one time step. Finally the polymer stressτPis computed from Eq.(3a)and fed into the LB equation(11), which is then integrated over one time step.
H R W
FIG. 3. The computational geometry corresponds to a repeat unit in a periodic, hexagonal array of cylinders, with W/H≈√3 and
R/H≈ 0.42.
V. SIMULATION PARAMETERS
We compute the pressure-driven, viscoelastic flow in a 2D geometry, which represents a small portion of a porous medium. Simulating in two instead of three dimensions allows studying wide ranges of parameters using reasonable amounts of computational resources. The flow geometry, which is shown in Fig. 3, corresponds to a repeat unit of a periodic, hexagonal array of cylinders of radius R, which resembles a fibrous medium. In fibrous media viscoelastic liquids show similar resistance versus flow rate curves as in porous media or packed beds; see, e.g., Ref. [17]. The pressure gradient is applied in the horizontal direction. The domain width W to domain height H ratio equals W/H ≈√3. The cylinder radius to domain height ratio equals R/H ≈ 0.42, corresponding to a solid volume fraction of 0.63. The number of grid points over the domain width and the domain height are 111 and 64, respectively, which corresponds to 27 grid points
over the cylinder radius. The solvent kinematic viscosity is
ν= x2
6 t.
We conduct one simulation of Newtonian flow (α= 0) and 72 simulations of viscoelastic flow, using three different polymer solutions, which are defined by the concentration parameter and the extensibility parameter: (α,b)= (100,10), (10,100), and (1,1000). In each of these solution the maximum ratio of the polymer stress to the Newtonian stress is αb= 1000. For each polymer solution, we simulate the flow through the porous medium at 24 values for the Weissenberg number:
We= √λu
k, (12)
varied between 2× 10−2 and 3× 103, by adjusting the value of λ. The permeability k in Eq.(12)is obtained by applying Eq.(1)to the Newtonian flow simulation: k≈ 9.5 × 10−3R2. Figure 4(a)shows Newtonian flow simulations [Eq. (11) withτP = 0] at different pressure gradients. Stokes’ condition is satisfied for Re= u√k/ν 10−1, meaning that the perme-ability is independent on Re [Eq.(1)]. In the remainder of the paper the pressure gradient is chosen such that the Reynolds number is Re≈ 10−2. For the 111× 64 grid this corresponds to a dimensionless time step of u t/ x≈ 6 × 10−4.
Figure 4(b) shows the permeability in Newtonian flow simulations [Eq.(1)] as a function of the grid resolution at fixed
10−2 10−1 100 7 8 9 10x 10 −3
Re
k R2 (a) 101 102 0.009 0.0095 0.01 0.0105 0.011 R Δ x k R2 (b) 10−1 100 101 10 20 30 40 50Pe
f
(c)FIG. 4. (a) Permeability [Eq.(1)] of Newtonian flow simulations as a function of the Reynolds number Re= u√k/ν. (b) Permeability [Eq.(1)] of Newtonian flow simulations as a function of the grid resolution. (c) Intrinsic polymer viscosity [Eq.(10)] in one way coupled simulations as a function of the P´eclet number Pe= u√k/κ, for R/ x= 13.4 (circles), 26.9 (squares), and 53.8 (triangles).
Re≈ 1 × 10−2. The 111× 64 grid (R/ x = 26.9) provides a permeability within 5% of the grid converged solution.
To study the grid convergence of the viscoelastic flow we numerically solve Eqs.(3) and(11)on three grids: 55× 32, 111× 64, and 222 × 128, corresponding to R/ x = 13.4, 26.9, and 53.8 grid points per cylinder radius, respectively. On each grid, we use a Weissenberg number of We= 102, an extensibility parameter of b= 103, a concentration parameter of α= 0 and various values for the polymer mass diffusivity
κ. These simulations are one-way coupled, meaning that the polymer dynamics is affected by the fluid flow, while the fluid flow is unaffected by the polymer stress, i.e., τP = 0 in Eq.(11). Figure4(c)shows the intrinsic polymer viscosity [Eq.(10)], as a function of the P´eclet number Pe= u√k/κ. The results from the 111× 64 grid are within a few percent of those of the 222× 128 grid.
In a realistic, experimental situation: √k= 10−6 m,
u= 10−2m s−1, κ= 10−11m2s−1, and Pe= 103. To resolve the corresponding c gradients requires the P´eclet number based on the grid spacing Pe x= u x/κ 1. Using Pe = 103,
R/H= 0.42, and √k/R≈ 0.1, the Pe x requirement is equivalent to H / x 104, which is several orders of mag-nitude larger than what we presently consider as numerically feasible. Instead of using such large grids, we resort to the 111× 64 grid. On this grid the condition of resolvedness reads Pe≈ Pe x 1. Numerical tests showed that for the 111 × 64 grid and We∼ 103, numerical stability requires Pe≈ 0.1, which is the value used in the viscoelastic flow simulations reported in the remainder of this paper.
VI. NUMERICAL RESULTS
The simulated intrinsic polymer viscosity f as a function of the Weissenberg number We is shown in Fig.5for different values of the extensibility parameter b. Quantity f is computed by combining Eqs.(4)and(10):
f = −dp dx k αμu− 1 α.
With increasing We, f first decreases and then it increases. Similar, nonmonotonic behavior has been observed experi-mentally; see, e.g., Ref. [2]. To shed light on this behavior, we decompose the energy dissipation rate [Eq.(6)] over the flow
10−2 10−1 100 101 102 103 10−2 10−1 100 101 We f
FIG. 5. The intrinsic polymer viscosity f versus the Weissenberg number We for b= 10 (circles), b = 100 (squares), and b = 1000 (triangles). δx δy δx δy δx δy (a) (b) (c)
FIG. 6. Three limiting cases of a 2D velocity gradient, cor-responding to (a) rotation (Q= −1), (b) shear (Q = 0), and (c) elongation (Q= 1).
topology parameter Q. This parameter is defined as the second invariant of the normalized velocity gradient tensor:
Q= S
2− 2
S2+ 2,
where S2 and 2 are the second invariants of the rate of strain tensor S= 1
2(∇u
T + ∇u) and of the vorticity tensor = 1
2(∇u
T− ∇u), respectively:
S2= 12(S : S), 2= −12( : ).
As illustrated in Fig.6, Q= −1, Q = 0, and Q = 1 corre-spond to pure rotational flow, pure shearing flow and pure elongational flow, respectively.
To determine the relative contributions of the different flow topologies to the energy dissipation, we write the energy dissipation as the following integral over Q space:
=
1 −1p(Q)
QN+ QPdQ,
where p(Q) is the probability density function of Q and QNand
P
Qare the Newtonian dissipation and the polymer dissipation, conditionally averaged on Q. Figures7(a)and7(b)show p(Q) in the Newtonian simulation and in the simulation with the strongest viscoelastic effects, b= 1000 and We ≈ 3 × 103. The p’s are centered around Q= 0, which corresponds to shear flow. The secondary peak at Q= 1, which corresponds to elongational flow, is observed to increase with increasing We.
We divide the Q domain into three subdomains: (−1, −1 3), (− 1 3, 1 3), and ( 1
3,1), corresponding to flows dominated by rotation, shear, and elongation, respectively. Accordingly we decompose the energy dissipation into rotational dissipation
N
rot+ rotP, shearing dissipation shearN + P
shear, and elongational
−10 −0.5 0 0.5 1 0.1 0.2 0.3 Q p (a) −1 −0.5 0 0.5 1 0 0.1 0.2 0.3 Q p (b)
FIG. 7. The probability density function p of the flow topology parameter Q for (a) Newtonian flow and (b) viscoelastic flow with
10−2 10−1 100 101 102 103 10−4 10−3 10−2 10−1 100 101 We
FIG. 8. The intrinsic polymer viscosity f versus the Weisenberg number We for b= 1000, decomposed over the different flow topologies, corresponding to rotational flow (frot, circles), shear flow (fshear, squares) and elongational flow (felong, triangles).
dissipation elongN + P elong:
= Nrot+ rotP + shearN + Pshear+ elongN + elongP .
These quantities are defined as integrals over the respective Q subdomains: = −1 3 −1 p(Q) QN+ QPdQ N rot+ Prot + 1 3 −1 3 p(Q) QN+ PQdQ N shear+ shearP + 1 1 3 p(Q) QN+ QPdQ N elong+ elongP .
Similar as outlined in Sec.III[Eqs.(4)–(8)], we define intrinsic polymer viscosities for the different flow topologies frot, fshear, and felong, such that
−dp
dx =
μ[1+ α(frot+ fshear+ felong)]u
k . (13) Assuming that N = N rot+ N shear+ N elong= μu 2/k, we can define fα= μP α μα = P α Nα = P αk μu2α,
where α∈ {rot, shear, elong}. As given in Eq.(13), the intrinsic viscosities fα correspond to the contribution in the pressure drop−dp/dx by the polymer stress in flow regions, which are dominated by rotation, shear, and elongation, respectively.
Figure 8 shows the intrinsic viscosities fα as functions of We for b= 1000. As expected, polymers in rotational flow regions do not contribute significantly to−dp/dx. The observed negative values for frot at We∼ 10 (absent data points in Fig. 8) correspond to an energy transfer from the polymers to the solvent, which can occur when the polymers relax in rotation dominated regions. In the Newtonian limit (We 1), the polymer contribution to −dp/dx is concen-trated in the shear regions, which reflects the dominance of this flow topology, as shown in Fig.7(a). Similar as in steady shear flow [Fig.2(b)], the decrease in f with We for small We is related to a decrease of the polymer viscosity in the
shear regions fshear. Furthermore the rise of felong with We for larger We is in line with the behavior in steady elongation [Fig.2(a)]. However, the increase of felong is much smaller than b2, that corresponds to steady elongational flow. Another unexpected result is that for large We, the fshear increases significantly, which is opposite to the situation in steady shear flow. The anomalous increase of fshearwith We reflects that the polymers in the shear-dominated regions in a porous medium have significant alignment with the strain rate. This is opposite to steady shear flow, where the polymers are predominantly orientated along an axis of zero strain rate.
VII. ANALYTICAL MODEL
Here we propose an analytical model for f at large We. Figure9shows that at large We the flexible polymers are fully extended, with the trace of the configuration tensor reaching the maximum value of b. Fully extended polymers behave as rigid rods and the corresponding stress equals [18]
τP = 2μαbS : δrδrδrδr,
(14) whereδr= cos θδx + sin θδy is the polymer orientation unit
vector, which is defined by the polar angle θ betweenδr and δx . Hereδx andδy span the local Cartesian coordinate system
that is orientated such that the local velocity gradient can be written as Eq.(17). In two dimensions the rotation of the rod equals [18]
˙ δr = ∇uT
:δrδθδθ, (15)
where the dot represents time differentiation and δθ= − sin θδx + cos θδy is the unit vector in the θ direction.
Com-bining Eqs.(10)and(14)gives the following expression for f :
f = bS :δrδrδrδr : S
S : S . (16)
In the model we compute the polymer dynamics in a reference frame, which moves with the center of mass of the polymer and which rotates such that in this frame the velocity gradient tensor can be written as the superposition of shear and elongation: ∇u = Gδyδx + E[δxδx− δyδy]. (17) 10−2 10−1 100 101 102 103 100 101 102 103 We tr c
FIG. 9. The squared, scaled polymer length trc versus the Weissenberg number We for b= 10 (circles), b = 100 (squares), and b= 1000 (triangles).
101 102 103 10−1 100 101 b f∞
FIG. 10. The intrinsic polymer viscosity f∞in the limit of large Weissenberg numbers We versus the polymer extensibility parameter
b. The circles are simulated values at We≈ 3 × 103. The dashed line is a fit of our model [Eq.(19)] to the simulation data: f∞= 1.3× 10−2b.
In a porous medium the shear rate G and the elongation rate
E are rather complicated functions of time. By combining Eqs. (15) and (17) we obtain the polymer rotation in the comoving and corotating reference frame:
˙
θ= −G2/2+ G2/4+ E2cos(2θ+ arctan[2E/G]). (18) Equation (18) has a stable point defined by ˙θ = 0 and ∂ ˙θ /∂θ <0. Inserting this stable point into Eq. (16)
gives
f =b
2
E2
E2+ G2/4, (19)
where the operator· · · averages a quantity in time as “seen” by a fluid element that moves through the porous medium.
Figure10shows the simulated intrinsic viscosities in the limit of large Weissenberg numbers f∞ as a function of b, which are taken from the simulations at We≈ 3 × 103. The data follow a linear relationship, in agreement with the model [Eq.(19)]. The dashed line is a fit of the model to the simulation data, which gives a shear to elongation ratio of G2/E2≈ 1.5× 102.
VIII. CONCLUSIONS
In agreement with experiments in the literature the present viscoelastic flow simulations in porous media predict shear thinning for small Weissenberg numbers We and shear thick-ening for large We. Due to the combination of shear and elongation, the shear thickening in the porous medium is a few orders of magnitude smaller than in steady elongational flow.
In the limit of large We the simulation results are rational-ized by modeling the flexible chains as rods and modeling the flow field in the porous medium by a superposition of shear and elongation. The relative magnitude of these flow types, which is a geometrical property, is the only free parameter of the model. Future work involves the study of this parameter for different porous geometries.
[1] R. Haas and F. Durst,Rheol. Acta 21, 566 (1982).
[2] A. Magueur, G. M. Moan, and G. Chaveteau, Chem. Eng. Commun. 36, 351 (1985).
[3] R. P. Chhabra, J. Comiti, and I. Machac,Chem. Eng. Sci. 56, 1 (2001).
[4] F. Galindo-Rosales, L. Campo-Dea˜no, F. Pinho, E. van Bokhorst, P. Hamersma, M. Oliveira, and M. Alves,Microfluid. Nanofluid. 12, 485 (2012).
[5] T. Koshiba, N. Mori, S. Sugiyama, and K. Nakamura,Rheol. Acta 38, 375 (1999).
[6] S. H. Momeni-Masuleh and T. N. Phillips,Comput. Fluids 33, 1075 (2004).
[7] K. K. Talwar and B. Khomami,J. Non-Newtonian Fluid Mech. 57, 177 (1995).
[8] A. W. Liu, D. E. Bornside, R. C. Armstrong, and R. A. Brown, J. Non-Newtonian Fluid Mech. 77, 153 (1998).
[9] R. B. Bird, R. C. Armstrong, O. Hassager, C. F. Curtiss, and S. Middleman, Dynamics of Polymeric Liquids, Vols. 1 and 2 (John Wiley and Sons, New York, 1987).
[10] I. Ghosh, G. H. McKinley, R. A. Brown, and R. C. Armstrong, J. Rheol. 45, 721 (2001).
[11] M. Ellero, P. Espa˜nol, and E. G. Flekkøy, Phys. Rev. E 68, 041504 (2003).
[12] O. B. Usta, J. E. Butler, and A. J. C. Ladd,Phys. Rev. Lett. 98, 098301 (2007).
[13] D. A. Wolf-Gladrow, Lattice-Gas Cellular Automata and Lattice
Boltzmann Models: An Introduction (Springer-Verlag, Berlin,
2000).
[14] X. He, Q. Zou, L.-S. Luo, and M. Dembo,J. Stat. Phys. 87, 115 (1997).
[15] J. Onishi, Y. Chen, and H. Ohashi,Prog. Comput. Fluid Dyn. 5, 75 (2004).
[16] X. He and L.-S. Luo, Phys. Rev. E 55, 6333 (1997).
[17] L. Skartsis, B. Khomami, and J. L. Kardos,J. Rheol. 36, 589 (1992).
[18] M. Doi and S. F. Edwards, The Theory of Polymer Dynamics (Clarendon Press, Oxford, 1986).
