- Email: [email protected]
Direct numerical simulations of hard particle suspensions in planar elongational flow
J. Non-Newtonian Fluid Mech. 136 (2006) 167–178
Direct numerical simulations of hard particle suspensions in planar elongational flow Wook Ryol Hwang a,∗ , Martien A. Hulsen b a
School of Mechanical and Aerospace Engineering, Research Center for Aircraft Parts Technology (ReCAPT), Gyeongsang National University, Jinju 660-701, South Korea b Department of Mechanical Engineering, Eindhoven University of Technology, P.O. Box 513, 5600MB Eindhoven, The Netherlands Received 28 February 2006; received in revised form 13 April 2006; accepted 13 April 2006
Abstract We present a new direct simulation technique of inertialess particle suspensions in planar elongational flow of a Newtonian fluid. The extensional bi-periodic domain concept is introduced such that a single cell problem with a small number of particles may represent a large number of repeated structures of such a cell in planar elongational flow. For implicit treatment of the hydrodynamic interaction between particles and fluid, we employ a finite-element/fictitious-domain method similar to the distributed Lagrangian multipliers (DLM) method together with a rigid-ring description of the particle. The extensional bi-periodic frame is incorporated by constraint equations with Lagrangian multipliers and is implemented by the mortar element method. In our formulation, the bulk stress is evaluated by simple boundary integrals. Concentrating on 2D circular disk particles, we present numerical examples of single-particle, two-particle and 100-particle problems in the extensional bi-periodic frame. We discuss effects of solid fraction and particle configuration on the elongational viscosity of the suspension, in comparison with simple shear flow. We found that, at zero strain, the relative elongational viscosity is almost the same as the relative shear viscosity in simple shear for moderately concentrated suspensions. There is a small increase in elongational viscosity for large strains, which is related to an anisotropic distribution of the particles. © 2006 Elsevier B.V. All rights reserved. Keywords: Direct numerical simulation; Suspension rheology; Extensional rheology; Elongational viscosity; Extensional bi-periodic frame
1. Introduction It is quite common to add particles to fluids and polymer fluids in particular. The reasons can be quite diverse, e.g. for the manufacturing of ceramics [1] where the polymer is removed after processing, changing the product properties (such as adding flame retardants) or modifying the properties during processing. Therefore, studying the rheological properties of particle-filled polymer is very important. In most studies so far only shear deformation has been considered, whereas in processing both shear and elongational deformation components are present. Studying particle-filled polymers in elongational flow is therefore very relevant for processing. Since our final objective is to study flow of particles in viscoelastic fluids in general flows and in elongational flows in particular, we need to use direct simulation methods that model
∗
Corresponding author. Tel.: +82 55 751 6295; fax: +82 55 757 5622. E-mail addresses: [email protected] (W.R. Hwang), [email protected] (M.A. Hulsen). 0377-0257/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.jnnfm.2006.04.004
the interaction between particles by a direct simulation of the fluid in between the particles. Studying systems with many particles using the finite-element method with boundary-fitted meshes would be too complicated, especially in 3D. Therefore, we employ the fictitious domain method (see, for example Ref. [2]) in which internal boundaries are created using Lagrangian multiplier techniques on simple domains, such as a rectangular region, that can be meshed rather easily. In this paper, we present a new direct numerical simulation technique for 2D hard particle suspensions in a Newtonian fluid under planar elongational flow. The technique is based on the fictitious domain method for solid particles in fluids in bi-periodic domains under shear as developed in [3] and [4], where the domain has been fixed. Here we extend this approach to biperiodic regions that deform in time according to the average deformation in planar elongational flow. In this way we are able to obtain bulk properties for (the start-up of) elongational flows. The methods employed can relatively easy be extended to viscoelastic flows and 3D domains. The paper is structured as follows. First in Section 2 we define the problem and explain how it is modeled. Then in Section 3 we
168
W.R. Hwang, M.A. Hulsen / J. Non-Newtonian Fluid Mech. 136 (2006) 167–178
describe the numerical methods employed and the implementation of these methods. In Section 4 we show results of three examples: a flow with a single particle, two particles and 100 particles, respectively. Finally we end with some conclusions. 2. Modeling 2.1. Problem definition We consider the suspension with a Newtonian fluid consisting of a large number of circular disk particles subject to planar elongational flow. The hydrodynamic interaction between particles and fluid leads to complicated particle motion and microstructural development for both fluid and particles. In order to deal with such a problem at reasonable computational cost, a welldefined bi-periodic domain concept needs to be introduced that may transform the suspension problem with a large number of particles into a particulate flow problem in a unit cell. In simple shear flow, the Lees–Edwards boundary condition (LEbc), proposed by Lees and Edwards [5] for Molecular Dynamics, satisfies the above requirement exactly with the help of the sliding boundary concept. In the previous studies, the authors extended the LEbc concept to the continuous field problems by combining it with the finite-element/mortar-element method, for particle suspension problems [3,4]. In non-equilibrium molecular dynamics simulation of elongational flow Heyes [6] introduced the concept of the deforming simulation cell, which is aligned in flow direction and stretches in that direction. The cross-flow dimensions of the cell decrease such that the volume of the cell remains constant. A disadvantage of this approach is that only relatively small strains can be reached: the contracting dimensions eventually become of the same order of magnitude as the interaction range of a single molecule. For complex molecules an equilibrium condition (steady state) has usually not been reached yet at this point. For planar elongational flows it is possible to use so-called Kraynik–Reivelt boundary conditions [7] and simulations can be performed indefinitely [8,9]. However this approach is difficult to use and cannot be extended to uniaxial elongational or bi-axial stretching flows. In this paper we adopt the simple deforming cell approach of Heyes [6] and propose the extensional bi-periodic frame as a
computational domain for the representative unit cell problem in planar elongational flow of particle filled fluids. Fig. 1 shows the extensional bi-periodic frames with the elongation rate ˙ . (The coordinate system is indicated therein as well.) The extensional frame is a material frame in that it translates and deforms along the flow kinematics. The hatched region in Fig. 1 indicates a single frame for the illustration. At an arbitrary instance, say t = 0, the unbounded domain is regularly divided into an infinite number of cells (or frames) of the width L0 and the height H0 . As time goes on, frames elongate in the horizontal direction and contract in the vertical direction, while translating along the flow. Given the elongation rate ˙ as in Fig. 1, the time-dependent changes of the width L and the height H of the frame can be expressed as L(t) = L0 exp(˙t),
H(t) = H0 exp(−˙t).
(1)
Although each frame translates at different velocity, the frame dimensions L(t) and H(t) is the same for all the extensional frames and for all time t. A position vector x˜ measured from an arbitrary global reference frame OXY can be decomposed as the sum of the position of the frame center XF and the relative position vector x based on the reference frame OF xy (Fig. 2): x˜ = XF + x Note that the relative position vector x is periodic, in contrast to the global position vector x˜ . Similarly, the velocity field can be expressed as ˜ x) = U F (XF ) + u(x), u(˜
(2)
˜ UF and u are the velocity of the material point with where u, respect to an arbitrary global reference, the frame velocity and the relative velocity inside the frame, respectively. The frame velocity UF , the velocity at the center of the extensional biperiodic frame, with respect to the global reference is given by: U F (XF ) = (˙XF , −˙YF ).
(3)
We note that, in the absence of particles, the material points inside a frame will deform in the same way as the frame and thus the relative velocity inside the frame u(x) is the same for all the extensional frames, regardless of the global location of
Fig. 1. Extensional bi-periodic frames in planar elongational flow. The movement and the deformation of a single frame is indicated for the illustration (hatched).
W.R. Hwang, M.A. Hulsen / J. Non-Newtonian Fluid Mech. 136 (2006) 167–178
169
Fig. 4. The extensional bi-periodic frame with a four-particle configuration at a certain instance, which is the computational domain in this study.
Fig. 2. The material position x˜ -based on an arbitrary global reference can be expressed as the vector sum of the frame center coordinate XF and the relative coordinate x.
the frame, i.e. it is periodic and it is given by: u(x) = (˙x, −˙y)
(4)
Indeed, the position of the frame with respect to the global reference affects only the frame velocity UF . It is essential to realize that u(x) will be periodic also for non-homogeneous fluid systems subject to planar elongational flow, e.g. particle suspensions or droplet emulsions, if the initial configuration of the particles or the droplets is periodic from the beginning. Fig. 3 shows a set of extensional bi-periodic frames at
a certain instance. There are an infinite number of such frames in the unbounded domain. Each frame contains three circular disk-like particles and two of them cross the frame boundary. Note that the particle configuration in all frames is the same. As shown in Fig. 3, the velocity field in each frame can be decomposed into the frame velocity UF , which has been given in Eq. (3), and the relative velocity u(x). The relative velocity in the frame 1 in Fig. 3 can be obtained by specifying the global zero velocity at OF1 . In the same way, by assigning zero velocity at OF2 (or OF3 ), one gets the relative velocity in the frame 2 (or 3). Then the relative velocities in the frames 1–3 can be considered the same, since there are an infinite number of repeated structures. This means that the velocity field in each frame can be decomposed into the frame velocity UF , as in Eq. (3), and the relative velocity u(x) will be the same for all frames and be periodic. In other words, all relative events in one frame happen in the other. However, care should be exercised in determining the size of the frame: the length scale of the frame should be large enough in comparison with the scale of the microstructure. In this study, we solve a single extensional bi-periodic frame with a number of rigid particles. The frame can be periodically extended to fill the full 2D space by adding the proper frame velocity Eq. (3). The word ‘frame’ seems to be adequate for the unit cell domain, since it is a small window through which one observes what happens inside the suspension. An extensional bi-periodic frame, denoted by Ω, is the computational domain of the work and the four boundaries of the frame are denoted by Γ i (i = 1, 2, 3 and 4) (Fig. 4). We use a symbol Γ for ∪4i=1 Γi . The Cartesian x and y coordinates are selected for the direction of the stretching and that of the contraction, respectively. Particles are denoted by Pi (t) (i = 1,. . ., N) and N is the number of particles in a single frame. We use a symbol P(t) for ∪N i=1 Pi (t), a collective region occupied by particles at a certain time t. For a particle Pi , Xi = (Xi ,Yi ), Ui = (Ui ,Vi ), ωi = ωi k and i = Θi k are used for the coordinates of the particle center, the translational velocity, the angular velocity and the angular rotation, respectively; and k is the unit vector in the direction normal to the plane. 2.2. Governing sets of equations The set of equations for the fluid domain is given by:
Fig. 3. A set of the extensional bi-periodic frames containing three particles at a certain instance. Some particles cross the frame boundary.
∇ · = 0,
in Ω/P(t),
(5)
170
W.R. Hwang, M.A. Hulsen / J. Non-Newtonian Fluid Mech. 136 (2006) 167–178
∇ · u = 0,
in Ω/P(t),
= −pI + 2ηD,
(6)
in Ω/P(t),
u = U i + ωi k(x − Xi )
on ∂Pi (t),
(7) (i = 1, . . . , N).
(8)
Eqs. (5)–(8) are equations for the momentum balance, the continuity, the constitutive relation, and rigid-body conditions on particle boundaries, respectively. u, σ, p, I, D and η are the velocity, the stress, the pressure, the identity tensor, the rate of deformation tensor and the viscosity, respectively. The unknown rigid-body motions (Ui ,ωi ) in Eq. (8) will be determined by the hydrodynamic interaction that will be discussed later. In the absence of inertia, initial conditions are not necessary for the velocity of the fluid as well as the particles. In addition, there is no explicit boundary condition on the domain boundary Γ . Instead we will introduce constraint equations on Γ in Section 2.3 to assign a planar elongational flow condition. Following the work by Hwang et al. [3], we consider the circular particle as a rigid ring, which is filled with the same fluid as in the fluid domain and the rigid body condition is imposed on the particle boundary only. Thus, one needs the discretization of the particle along the boundary only. From the rigid-ring description, the set of governing equations for a region occupied by a particle Pi at a certain time t is: ∇ · = 0,
in Pi (t),
(9)
∇ · u = 0,
in Pi (t),
(10)
= −pI + 2ηD,
in Pi (t),
u = U i + ωi k × (x − Xi ),
on ∂Pi (t).
(11) (12)
Eqs. (9)–(12) are the equations for the momentum balance, the continuity, the constitutive relation, and the boundary condition, respectively, which are exactly the same as the fluid domain equations in Eqs. (5)–(8). In addition, the movement of the particle is given by the kinematic equations: dXi = U i, dt
Xi|t=0 = Xi,0 ,
(13)
dΘi = i , dt
Θi|t=0 = Θi,0 .
(14)
Eq. (14) is completely decoupled from the other equations, since the particle is circular. To determine the unknown rigid body motions (U i , i )’s of the particles, one needs balance equations for drag forces and torques on particle boundaries. In the absence of inertia and external forces or torques, particles are force-free and torquefree: Fi = · n ds = 0, (15) Ti =
∂Pi (t)
∂Pi (t)
(x − Xi ) × ( · n)ds = 0,
(16)
where Ti = Ti k and n is a normal vector on ∂Pi pointing out of the particle (i = 1,. . ., N). We did not use an artificial particleparticle collision scheme [2], because the particle overlap and
particle/wall collision could be avoided for the multiple-particle problems we studied in this paper by taking a relatively small time-step and a sufficiently refined particle boundary discretization [3]. 2.3. Extensional bi-periodic frame constraints Here we discuss the mathematical description of the biperiodicity in the extensional frame. The continuity of the velocity field and the force balance need to be satisfied across the domain boundary Γ . The conditions in the horizontal direction (between Γ 2 and Γ 4 ) at a certain time t can be written as follows: u(0.5L(t), y) = u(−0.5L(t), y) + f h , y ∈ [−0.5H(t), 0.5H(t)]
(17)
t(0.5L(t), y) = −t(−0.5L(t), y), y ∈ [−0.5H(t), 0.5H(t)],
(18)
where the vector t is the traction force on the boundary and f h = (˙L(t), 0). With the elongation rate e as shown in Fig. 2, the velocity in the x direction on the right boundary (Γ 2 ) is faster than that on the left boundary (Γ 4 ) by the amount of the velocity difference ˙ L(t). Similarly, the conditions for the periodicity in the vertical direction between Γ 1 and Γ 3 are u(x, 0.5H(t)) = u(x, −0.5H(t)) + f v , x ∈ [−0.5L(t), 0.5L(t)] t(x, 0.5H(t)) = −t(x, −0.5H(t)),
(19) x ∈ [−0.5L(t), 0.5L(t)]. (20)
The velocity difference in the vertical direction is specified by f v = (0, −˙H(t)) in Eq. (19). In the weak formulation of the finite element method, the kinematic constraint is usually combined with the Lagrangian multipliers and then the force balance can be satisfied implicitly through the multipliers. In this regard, we use the kinematic equations Eqs. (17) and (19) only in the derivation of the weak formulation and we call the two equations the extensional biperiodic frame constraints. 2.4. Boundary-crossing particles When a particle crosses the boundary of the computational domain, the particle or parts of the particle need to be relocated into the domain. The relocation of the particle always involves the change of the translational velocity of the particle in the extensional bi-periodic frame to satisfy the continuity of the velocity field. An example situation is described in Fig. 5: the part of the particle which is located outside the computational domain is relocated into the domain according to the bi-periodicity. As indicated, the relocation involves the change of the particle velocity. In order to describe the relocation, we consider two sets of the coordinates: the unprimed set for the original coordinate (before relocation) and the primed set for the relocated coordinate. The relocated coordinate x = (x ,y ) for the
W.R. Hwang, M.A. Hulsen / J. Non-Newtonian Fluid Mech. 136 (2006) 167–178
171
in the vertical direction (Eq. (19)) and the rigid-ring constraint equation along the i-th particle boundary (Eqs. (8) and (12)): h = (hx , hy ) ∈ L2 (Γ4 )2 ,
v = (vx , vy ) ∈ L2 (Γ3 )2 ,
p,i 2 2 p,i = (p,i x , y ) ∈ L (∂Pi (t)) .
Fig. 5. An illustration of a particle crossing the frame boundary: the part of the particle which is present outside the domain is relocated into the domain according to the bi-periodicity.
original coordinate x on the particle boundary is determined as follows:
In Ref. [3], we showed that the multipliers h and v are the traction force on the domain boundary Γ and that the multiplier p,i is related to the traction force acting on the particle boundary ∂Pi (t). The weak form for the problem can be stated as follows:At a certain time t with L(t), H(t) and for a given particle configuration Xi (i = 1,. . ., N), find u ∈ H1 (Ω)2 , U i ∈ R2 , ωi ∈ R, p,i ∈ L2 (∂Pi (t)), p ∈ L2 (Ω), h ∈ L2 (Γ 4 ), and v ∈ L2 (Γ 3 ) (i = 1,. . ., N) such that − p∇ · v dΩ + 2ηD[u] : D[v] dΩ Ω
x = mod(x + 0.5L, L) − 0.5L, y = mod(x + 0.5H, H) − 0.5H
+
(21)
where ‘mod’ is the modular function: e.g. mod (1.7X,X) = 0.7X and mod (−0.7X,X) = 0.3X for a positive real number X. The accompanying modification for the translational velocity of the particle associated with the relocated point x , denoted by U = (U ,V ), can be given as follows: ⎧ ⎪ ⎨ U − ˙ L(t), x > 0.5L(t), U = U + ˙ L(t), x < −0.5L(t), ⎪ ⎩ U, otherwise. ⎧ ⎪ ⎨ V + ˙ H(t), y > 0.5H(t), V = V − ˙ H(t), y < −0.5H(t), (22) ⎪ ⎩ V, otherwise. where U = (U,V) denotes the translational velocity of the particle associated with the original coordinate x. Note that the angular velocity does not change during the relocation. 3. Numerical methods 3.1. Combined weak formulation
Ω
N i=1
((p,i , v − [V i − χi k × (x − Xi )])∂Pi
+ (h , v(0.5L, y) − v(0.5L, y))Γ4 + (v , v(x, 0.5H) − v(x, −0.5H))Γ3 = 0
(23)
q∇ · u dx = 0
(24)
Ω p,i
(µ , u − [U i − ωi k × (x − Xi )])∂Pi = 0,
(i = 1, . . . , N) (25)
h
(µ , u(−0.5L, y) − u(0.5L, y))Γ4 = (µ , f h )Γ4 , v
f h = (˙L(t), 0)
(26)
(µv , u(x, 0.5H) − u(x, −0.5H)Γ3 = (µv , f v )Γ3 , f v = (0, −˙H(t))
(27)
for all v ∈ H 1 (Ω)2 , V i ∈ R2 , χi ∈ R, µp,i ∈ L2 (∂Pi (t)), q ∈ L2 (Ω), µh ∈ L2 (Γ 4 ), and µv ∈ L2 (Γ3 ). The inner product (·, ·)Γj is the standard inner product in L2 (Γ j ): (µ, v)Γj = µ · v ds. Γj
We have several remarks on the weak form Eqs. (23)–(27): Following the combined weak formulation of Glowinski et al. [2] in which the hydrodynamic force and torque acting on the particle canceled exactly, we derived the weak form for the Newtonian system together with the rigid-ring description of the particle and with the sliding bi-periodic frame in the simple shear flow in our previous study [3]. The only change in this work from the previous work is the extensional bi-periodic frame constraint (Eqs. (17) and (19)). Therefore, we present the final weak form without the detailed derivation that is available in Ref. [3]. We introduce three different kinds of Lagrangian multipliers h , v and p,i which are associated with the kinematic constraint equation for the periodicity in the horizontal direction (Eq. (17)), the constraint equation for the sliding periodicity
• The rigid-body motion on the particle boundary is treated via the rigid-ring constraint, which is satisfied in a weak sense with the multiplier λp (Eqs. (23) and (25)); the extensional bi-periodic frame constraint is satisfied again weakly with the multipliers λh and λv (Eqs. (23), (26) and (27)). • Since the extensional bi-periodic constraint (Eqs. (17) and (19)) impose only the relative relation in the velocity, one needs to specify the reference velocity at a single point in the domain. The choice is arbitrary. We specify the zero velocity at the center of the domain (x,y) = (0,0) throughout the study. • The pressure level of the fluid domain is determined by specifying one of the normal component of the Lagrangian mul-
172
W.R. Hwang, M.A. Hulsen / J. Non-Newtonian Fluid Mech. 136 (2006) 167–178
tipliers on Γ (λhx or λvy ) since the multiplier can be identified as the traction force. • For a given particle configuration, L(t) and H(t), one can solve Eqs. (23)–(27) to get all the solution and then determine the next particle configuration and the domain configuration (dimension) using Eq. (13) and Eq. (1), respectively. • When a particle crosses the domain boundary, a part of the particle needs to be relocated as described in Eqs. (21) and (22). In this case, Eqs. (23) and (25) need to be modified as follows: (µp,i (x ), u(x ) − {U i + i × (x − Xi )})∂Pi = (µp,i (x ), f )∂Pi ,
(28)
where f = (fx , fy ) and it depends on the original coordinate (before relocation): ⎧ ⎪ ⎨ −˙L(t), x > 0.5L(t), x < −0.5L(t), fx = ˙ L(t), ⎪ ⎩ 0, otherwise, ⎧ y > 0.5H(t), ⎪ ⎨ ˙ H(t), fy = −˙H(t), y < −0.5H(t), (29) ⎪ ⎩ 0, otherwise.
Fig. 6. The regular rectangular discretization is used for the entire computational domain and the particle is described by the collocation points along the particle boundary.
where A = L(t)H(t) is the area of the extensional frame that is constant. The first terms in the right-hand side of Eq. (30) are the contribution from the boundary traction (i.e. λh or λv ) and the second ones are the contribution from the boundary-crossing particles. The latter vanishes when all particles are completely immersed in the extensional frame: i.e., x = x. 3.3. Implementation
3.2. Bulk stress The bulk stress, the average stress over the domain, can be expressed, for a volume V, as the sum of the fluid contribution and the particle contribution [10]: 1 1 1 σ = σ dV = dV + σ dV, V V V Vf V Vp where · denotes the averaged quantity in V, Vf and Vp are the volume occupied by the fluid and the particle, respectively. In the weak formulation of the sliding bi-periodic frame, we have the boundary integral expression for the bulk stress using the identity between the traction force and the Lagrangian multipliers λh , λv and λp which was discussed in detail in Ref. [3]. Such a relationship is still valid in the extensional bi-periodic frame. In this regard, we present the final expression of the bulk stress using the boundary integral without detailed proof: σ11 =
1 H(t)
Γ4
A
∂Pi
i=1
N
1 λvy ds + A Γ3
1 σ22 = − L(t) 1 σ12 = L(t)
λhx ds +
N 1
i=1
Γ3
N
λhy
1 ds + A i=1
Mi k=1
µk · {u(xk ) − (U i + ωi × (xk − Xi ))}, p,i
(y − y )λp,i y (x ) ds,
(x − x )λp,i y (x ) ds,
(30)
(31)
where Mi , xk , xk and µk are the number of the collocation points on ∂Pi , the original (before relocation) of the k-th collocation point, the relocated coordinate of the k-th collocation point and the collocated Lagrangian multiplier at xk , respectively. We define equally distributed collocation points on the particle boundary and the number of the collocation points M is taken to be proportional to the radius of the particle. In our p,i
∂Pk
∂Pi
≈
(x − x )λp,i x (x ) ds,
The spatial discretization for the weak form in Eqs. (23)–(27) is basically the same as our previous work on the sliding biperiodic frame in simple shear flow [3], except for the extensional bi-periodic frame constraint equations in Eqs. (26) and (27). However, for the completeness of the paper, we briefly summarize the method. For the discretization of the momentum equation in the weak form, we use regular quadrilateral elements with continuous bi-quadratic interpolation (Q2 ) for the velocity u, discontinuous linear interpolation (P1 ) for the pressure p for the entire computational domain, including the interior of the particle. An example mesh for the two-particle problem is presented in Fig. 6. To discretize the weak form of the rigid-ring constraint (Eqs. (23), (25) or (28)), we use the point collocation method. For example, the integral in Eq. (28) has been approximated as follows: (µp,i (x ), u(x ) − {U i + ωi × (x − Xi )})∂Pi
W.R. Hwang, M.A. Hulsen / J. Non-Newtonian Fluid Mech. 136 (2006) 167–178
previous work, we analyzed effects of the number of collocation points and found that approximately one collocation point in an element appears to give the most accurate solution [3]. A small change in M from the optimal number does not affect solutions much, but it gives additional control to avoid particle collision, although it hardly occurs in our simulation with a reasonable time step size. To discretize the boundary integrals for the horizontal periodicity in Eq. (26), we use the mortar element method for the integral in Eq. (27) using the continuous linear interpolation of the λh . For the vertical periodicity, we employed the nodal collocation method, point collocation at all nodes. Since the facing elements between Γ 2 and Γ 4 or between Γ 1 and Γ 3 are always conforming due to the regular discretization of the computational domain, one could use the nodal collocation method to implement the periodicity constraint for both directions. The latter method found to generate identical results with the method used in this work. For each time step, we have the change in the particle configuration Xi (i = 1,. . ., N) and the change of the computational domain, i.e. L(t) and H(t). We remesh the computational domain to keep the aspect ratio of the element almost one. For given Xi (tn ) (i = l,. . ., N), L(tn ) and H(tn ), we performed the following procedures: • Step 1: solve the weak form at tn in Eqs. (23)–(27) to get (u, p, Ui , ωi , λp,i , λh , λv ) for the n-th step. • Step 2: update the particle configuration Xi (tn+1 ) by solving the kinematic equation of Eq. (13) with the explicit Euler method. • Step 3: update the computational domain configuration L(tn+1 ) and H(tn+1 ) using Eq. (1). In the first step, an equation with a sparse symmetric matrix with many zeroes on the diagonal appears as a result of the above discretizations, which has been solved by a direct method based on the sparse multi frontal variant of the Gaussian elimination (HSL2002/MA41) [11]. In the step 3, the remeshing of the computational domain is done by the simple scheme below: L H Nx = int , Ny = int , (32) le le
173
where Nx , Ny , le and int are the number of discretization in the x and y direction, the predetermined element size and the integer truncation operator, respectively. 4. Example problems In this section, we present three example problems to demonstrate the feasibility of our formulation and the numerical schemes: single particle, two particle, and 100 particle problems. Because of the bi-periodicity of the computational domain, the result from a single frame computation with a number of particles can be considered as the problem with the periodically extended configuration of particles in an unbounded domain. 4.1. The single particle problems The first problem is stated as follows: a single particle of radius r is suspended freely at the center of an extensional biperiodic frame with L0 = H0 = 1 under the given elongation rate ˙ = 1 in a Newtonian fluid having a viscosity η = 1. We used a 50 × 50 fluid mesh for the single particle problem, except for the r = 0.05 and 0.075 cases in which a 100 × 100 mesh is employed. The time step t is 0.02. As shown in Fig. 7, this problem represents a regular configuration of an infinite number of particles. It is an extreme case such that the particles in the neighboring frames form vertical strings as time goes on. Therefore one can expect excessively high flow resistance (elongational viscosity) when the particles approach each other (when H = 2r, the elongational viscosity is infinite). Although the problem is an artificial one, it will show that the elongation viscosity increases for an anisotropic particle distribution. This result can be used to explain the increase in viscosity for the large particle systems that will be analyzed later. Fig. 8 shows an example result for a r = 0.15 particle, plotted with the distribution of the total shear rate (II2D ). A high shear rate region is found both on the top and the bottom of the particle due to the existence of other particles in the upper and the lower frames, when the extensional strain is large (˙t = 0.8). The bulk elongational viscosity η¯ 1 can be defined as η¯ 1 =
σ11 − σ22 ˙
Fig. 7. The single particle problem in the extensional bi-periodic frame represents a regular configuration of an infinite number of particles.
(33)
174
W.R. Hwang, M.A. Hulsen / J. Non-Newtonian Fluid Mech. 136 (2006) 167–178
Fig. 8. The single particle problem result in case of r = 0.15. The plotted quantity is the distribution of the shear rate II2D , the second invariant of the rate-ofdeformation tensor 2D: ˙ t = 0, 0.4, 0.8 (from the top to the bottom).
in the planar elongational flow. For the pure Newtonian fluid in 2D, the elongational viscosity becomes four times the shear viscosity η : η¯ 1 = 4η (The factor is 3, for uniaxial elongational flow in 3D, which is known as the Trouton ratio.) In this regard, we could define the relative bulk elongational viscosity as follows: η¯ 1 η¯ 1r = . (34) 4η The elongational viscosity increases with the presence of the particle. For the dilute system where the interaction between particles is neglected, it increases linearly with the solid area fraction φ with the factor of 2 [12]. The same linear relationship is valid not only for simple shear flow but also for planar elongational flow. In Fig. 9, we present the relative bulk elongational viscosity as a function of the extensional strain = ˙ t from the single particle problem from a very dilute system with φ = 0.079 (r = 0.05) to a concentrated system with φ = 0.38 (r = 0.35). The viscosity increases with , since the vertical distances between particles decrease exponentially that leads to large flow resistance, as mentioned above. Also the viscosity increases with φ. Since the governing equations and the boundary condition are linear, one can expect a symmetric viscosity result for the negative strain < 0, i.e. contraction in the horizontal direction but stretching in the vertical direction. Thus, in the single particle problem, the elongational viscosity for zero strain = 0 reaches a minimum value, where both horizontal and vertical distances between the particle are the same. The elongational viscosity at zero strain is of interest, since the comparison between elongational and shear viscosities can be made for the same particle configuration in the same domain size L = H = 1. Fig. 10 shows the comparison of the relative bulk shear viscosity ηr and the relative elongational viscosity η¯ 1r at zero strain. The simple shear flow computation has been performed by using our previous code on particle suspensions in sliding bi-periodic frame [3].
Fig. 9. The relative elongational viscosity η¯ 1r as a function of the extensional strain = ˙ t for various solid area fraction φ.
We remark that care should be exercised to interpret the results from the single particle problem with an artificial configuration of the particle. The excessive increase of the elongational viscosity for the high φ value cases is largely due to the regular configuration of the particle.
Fig. 10. The comparison of the relative bulk elongational viscosity η¯ 1r and the relative shear viscosity ηr at zero strain as a function of the solid area fraction φ from the single particle problem. The relative shear viscosity has been adopted from Ref. [3].
W.R. Hwang, M.A. Hulsen / J. Non-Newtonian Fluid Mech. 136 (2006) 167–178
175
Fig. 11. The initial configuration of the two-particle problem in the extensional bi-periodic frame.
4.2. The two-particle problem The two-particle problem is constructed to investigate the effect of the hydrodynamic interaction between particles on the bulk suspension behavior and is stated as follows: two identical particles of radius r = 0.1 are suspended initially in an extensional bi-periodic frame of L0 = H0 = 1 under the elongation rate ˙ = 1 and the fluid viscosity η = 1. The fluid domain is discretized by a 50 × 50 mesh and the time step is 0.02. As shown in Fig. 11, the two particles are located symmetrically about the center of the domain (0, 0) in the beginning: one particle center coordinate is at (−d, 0.25) and the other is at (d, 0.25). The parameter d is the offset from the center line and it is roughly a measure of (head-to-head) interaction between two particles. When d = 0, the two particle align vertically and they will end up with head-to-head collision as the single particle problem, whereas such a interaction will be minimized for the case d = 0.25. Fig. 12 shows the particle trajectories of the two-particle problem. The trajectories for the zero offset case appears to be just straight vertical lines and the two particles approach gradually as time goes on. For the other value of d, the trajectories roughly follow the given flow field of the planar elongational flow. Fig. 13 is the bulk relative elongational viscosity η¯ 1r of the two-particle problem as a function of the extensional strain for different offset d’s. The elongational viscosity for d = 0 increases the fastest, while that for d = 0.25 increases the slowest with respect to the extensional strain. These behaviors can be understood easily, since the height of the frame H(t) that leads to the infinite viscosity is twice of the particle diameter (4r) for the zero offset case and about the particle diameter (2r) for 0 < d ≤ 0.25, leading to asymptotes near = ln 2.5 ≈ 0.92 and ln 5 ≈ 1.61. For the intermediate values of d (0 < d < 0.25, so excluding d = 0 and d = 0.25), the elongational viscosity increases for small strain but later decreases with the strain under further strain. When the two particles mostly align in the vertical direction, the viscosity increases with the strain as is the case of d = 0; but, as the particles separate from each other
Fig. 12. The particle trajectories for various d values in the two-particle problem of r = 0.1 and ˙ = 1 (the solid fraction φ = 0.0628).
along the flow, the interaction between the particle becomes less and the viscosity behavior becomes similar to that of the case d = 0.25. It is interesting that, if there is a non-zero offset, even though it is small, the bulk viscosity behaves independent of the offset value after a certain value of the strain and that the effect of the offset parameter d is found to be negligible, since the non-zero offset guarantees the exponential separation of particles. This implies that there would not be a sharp increase of the bulk viscosity due to the (vertical) alignment of particles, for the suspension with the random particle
Fig. 13. The bulk relative elongational viscosity η¯ 1r as a function of the extensional strain = ˙ t in the two-particle problem of r = 0.1 for various d values φ = 0.0628).
176
W.R. Hwang, M.A. Hulsen / J. Non-Newtonian Fluid Mech. 136 (2006) 167–178
Fig. 14. An example 100 particle problem result for r = 0.025. The plotted quantity is the distribution of the shear rate II2D : ˙ t = 0, 0.5 and 1 (from the top to the bottom).
configuration, which will be illustrated through the next example. 4.3. Many particle problems The two previous problems are more or less restricted in investigating the suspension behavior because of their regular particle configuration, in the last example problem, we attempt randomly distributed many particles in the extensional bi-periodic frame in order to observe evolution of the particle structures and bulk suspension properties. The problem is stated as follows: randomly distributed 100 particles of the same radius r are suspended initially in an extensional bi-periodic frame of L0 = H0 = 1 under the given elongation rate ˙ = 1 and the fluid viscosity η = 1. The fluid domain is discretized by a 100 × 100 mesh and the time step is 0.02. Fig. 14 shows an example result with a particle radius r = 0.025 (φ = 0.196) for the strain ˙ t = 0, 0.5 and 1 with the distribution of the total shear rate (II2D ).1 In this case, the vertical distances between particles decrease exponentially, while the horizontal distance increases, due to the planar elongational flow. However, there is no vertical alignment of particles, since randomly distributed particles and their interactions always keep the particles from being aligned: particles continuously fill in the space between the exponentially separating other particles. Therefore, one can expect a much smaller increase of the bulk elongational viscosity with respect to the applied strain in many particle problems than in the single- or two-particle problems. Plotted in Fig. 15 is the relative bulk elongational viscosity computed in the 100 particle problems of several different solid area fraction φ with respect to the extensional strain . Some results are readily comparable with the results in the single particle problem presented in Fig. 9, since they have exactly the same solid fraction. The comparison demonstrates the effect of the particle configuration on 1 The movie has been supplemented with this article in the online version, at: doi:10.1016/j.jnnfm.2006.04.004.
Fig. 15. The bulk relative elongational viscosity η¯ 1r of the 100 particle problems as a function of the extensional strain = ˙ t for several values of the solid area fraction φ.
the bulk behavior. The relative elongational viscosity η¯ 1r of the 100 particle problem increases much slower with respect to e than the single particle problem: e.g., for φ = 0.196, the viscosity goes up to 5 at ≈ 0.6 in the single particle problem in Fig. 9, while it remains around 1.6 in the 100 particle problem in Fig. 15. We see consistent increase in the elongational viscosity around the strain ≈ 1 for all cases in Fig. 15. The increase might be related to the squeezing particles in the vertical direction, like in the one and two particle problems. However here, the alignment can be avoided, since particles fill in the space generated by the exponentially separating other particles. Due to this settling of particles, the elongation viscosity does not increase excessively and, sometimes, a decrease in the viscosity even occurs afterward (φ = 0.159 and 0.126 in Fig. 15). Plotted in Fig. 16 are the average horizontal and vertical distances between two nearby particles for φ = 0.196 (r = 0.025). The horizontal distance for each particle to other particles is defined by taking the horizontal distance to the closest particle within a angular window of within ±45◦ about the x-axis. The vertical distance is defined similarly but now using ±45◦ about the y-axis. It can clearly be seen in Fig. 16 that the average horizontal distance slightly increases and the average vertical distance slightly decreases as a function of strain. This means that the particle distribution becomes slightly anisotropic, which can also be noticed in Fig. 14. We have shown for the single particle case that the viscosity increases when the particle distribution becomes anisotropic. It seems likely that this is true for the many particle case also, so we believe that the anisotropic distribution here explains the slight increase in elongational viscosity as a function of strain.
W.R. Hwang, M.A. Hulsen / J. Non-Newtonian Fluid Mech. 136 (2006) 167–178
Fig. 16. The average vertical and horizontal distances between nearby particles for φ = 0.196 in the extensional bi-periodic frame.
As we did in the single particle problem, it would be worthwhile to compare the relative elongational viscosity of planar elongational flow and the relative shear viscosity of simple shear flow for the 100-particle problem. To get the results for the simple shear, we have performed the 100 particle problem set in a sliding bi-periodic frame of the size 1 × 1 at t = 0, using the same fluid and particle meshes as the elongational flow [3]. The imposed shear rate γ˙ = 1 and the fluid viscosity η = 1. Again we emphasize that the comparison here has been made only at the initial particle configuration at zero strain at which the computational domain size and the particle configuration for both problems are the same. Therefore, this comparison does not reflect the evolutionary effects of particle configurational structures. Plotted in Fig. 17 is the comparison result of the relative bulk elongational viscosity and the relative bulk shear viscosity of the 100 particle problem, along with the single particle results for both flows that have been presented in Fig. 10. The result show that both the relative viscosities for the 100 particle problems are almost the same for 0.07 < φ < 0.25 for the given particle configuration at zero strain. Similar results has been reported experimentally in Greener and Evans [1] such that the Trouton ratio is close to the theoretical value of 3 for moderate solid volume fractions in uniaxial elongational flow. In addition, Fig. 17 shows that the prediction from the single particle problem of planar elongational flow is excessive. Small differences in the relative viscosity for elongational and shear flows are believed to be caused by the difference in hydrodynamic interaction. In Fig. 18, we present an example 100 particle problem result for r = 0.025 in simple shear flow for comparison of the hydrodynamic interaction with planar elongational flow. Again, the plotted quantity is the distribution of the shear rate II2D . Comparing Fig. 18 and the top figure in Fig. 14, one
177
Fig. 17. The comparison of the bulk relative elongational viscosity η¯ 1r as a function of the solid fraction φ for the 100 particle problems in the planar elongational flow (denoted ‘pef’) at the zero strain and in the simple shear flow (denoted ‘ssf’) along with the single particle problem results in Fig. 10.
finds that the direction of the high sheared region in between particles is different for simple shear and planar elongational flows. In simple shear flow, the high shear region is formed in the diagonal directions, whereas it is oriented in the vertical and horizontal direction in planar elongational flow. This might be responsible for the difference in the relative viscosity in Fig. 17. Further investigation on the particle microstructural effects on the bulk suspension properties is necessary but this will be part of our future work.
Fig. 18. An example 100 particle problem result for r = 0.025 in a sliding biperiodic frame for comparison with that of planar elongational flow (the top figure in Fig. 14). The plotted quantity is the distribution of the shear rate II2D .
178
W.R. Hwang, M.A. Hulsen / J. Non-Newtonian Fluid Mech. 136 (2006) 167–178
5. Conclusions
Acknowledgements
In this work, we developed a new bi-periodic domain concept, the extensional bi-periodic frame, for the planar elongational flow, which is quite well suited for studying rheology and microstructural development in non-homogeneous multiphase materials such as droplet emulsions, particle suspensions, etc. We combined the extensional bi-periodic frame with particle suspensions in 2D and developed numerical techniques for implementation using the finite-element/fictitiousdomain/mortar-element methods. Our method is able to deal with particles crossing the computational boundary and the bulk stress can be expressed by simple boundary integrals. Concentrating on disk-like particle suspensions in 2D, we studied numerical examples of single-particle, two-particle and 100-particle problems in the ex-tensional bi-periodic frame. Through the examples, we have discussed effects of the solid fraction and the particle configuration on the elongational viscosity of particle suspensions in comparison with that in simple shear flow. In the single particle problem with a regular particle configuration, we obtained the bulk elongational viscosity for a wide range of the solid fraction and reported that the elongational viscosity increases much faster than the shear viscosity with respect to the solid fraction, on account of the head-tohead collision. Also the elongational viscosity has been found to have a minimum value at the instance when both horizontal and vertical distances are the same. In the two-particle problem, we discussed the effect of the particle interaction on the bulk elongational viscosity and found that the bulk elongational viscosity behaves almost the same irrespective of the initial particle configuration after a certain amount of the strain, unless the particles are aligned exactly on the contraction axis. Finally, from the 100 particle problems, we found that the dependence of the bulk elongational viscosity on the amount of the strain is minor in comparison to the single particle problem. We also reported that the relative bulk shear and elongational viscosities are found almost the same for moderately concentrated suspensions at zero strain, by comparison with the suspension simulation with the sliding bi-periodic frame in simple shear flow. We remarked that the anisotropy in the particle microstructure is responsible for a small increase in the elongational viscosity with increasing strain.
The authors thank Mr. Michael Baltussen for his early work on this topic in his Bachelor final project. This work was supported by the Korea Research Foundation Grant (KRF-2005005-J09902) funded by the Korean government (MOEHRD) and partially by the Dutch Polymer Institute (DPI). Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.jnnfm.2006.04.004. References [1] J. Greener, J.R.G. Evans, Uniaxial elongational flow of particle-filled polymer melts, J. Rheol. 42 (1998) 697–709. [2] R. Glowinski, T.-W. Pan, T.I. Hesla, D.D. Joseph, A distributed Lagrangian multipliers/fictitious domain method for participate flows, Intern. J. Multiphase Flow 25 (1999) 755–794. [3] W.R. Hwang, M.A. Hulsen, H.E.H. Meijer, Direct simulation of particle suspensions in sliding bi-periodic frames, J. Comput. Phys. 194 (2004) 742–772. [4] W.R. Hwang, M.A. Hulsen, H.E.H. Meijer, Direct simulation of particle suspensions in a viscoelastic fluid in sliding bi-periodic frames, J. NonNewtonian Fluid Mech. 121 (2004) 15–33. [5] A.W. Lees, S.F. Edwards, The computer study of transport processes under extreme conditions, J. Phys. C: Solid State Phys. 5 (1972) 1921–1929. [6] D.M. Heyes, Molecular dynamics simulations of extensional, sheet and unidirectional flow, Chem. Phys. 98 (1985) 15–27. [7] A.M. Kraynik, D.A. Reinelt, Extensional motions of spatially periodic lattices, Intern. J. Multiphase Flow 18 (1992) 1045– 1059. [8] B.D. Todd, P.J. Daivis, Nonequilibrium molecular dynamics simulation of planar elongational flow with spatially and temporally periodic boundary conditions, Phys. Rev. Lett. 81 (1998) 1118–1121. [9] A. Baranyai, P.T. Cummings, Steady state simulation of planar elongation flow by nonequilibirum molecular dynamics, J. Chem. Phys. 110 (1999) 42–45. [10] G.K. Batchelor, The stress system in a suspension of force-free particles, J. Fluid Mech. 41 (1970) 545–570. [11] HSL2002, A collection of Fortran codes for large scale scientific computation, 2002, http://www.numerical.rl.ac.uk/hsl. [12] J.F. Brady, The Einstein viscosity correction in n dimensions, Intern. J. Multiphase Flow 10 (1984) 113–114.
Direct numerical simulations of hard particle suspensions in planar elongational flow Wook Ryol Hwang a,∗ , Martien A. Hulsen b a
School of Mechanical and Aerospace Engineering, Research Center for Aircraft Parts Technology (ReCAPT), Gyeongsang National University, Jinju 660-701, South Korea b Department of Mechanical Engineering, Eindhoven University of Technology, P.O. Box 513, 5600MB Eindhoven, The Netherlands Received 28 February 2006; received in revised form 13 April 2006; accepted 13 April 2006
Abstract We present a new direct simulation technique of inertialess particle suspensions in planar elongational flow of a Newtonian fluid. The extensional bi-periodic domain concept is introduced such that a single cell problem with a small number of particles may represent a large number of repeated structures of such a cell in planar elongational flow. For implicit treatment of the hydrodynamic interaction between particles and fluid, we employ a finite-element/fictitious-domain method similar to the distributed Lagrangian multipliers (DLM) method together with a rigid-ring description of the particle. The extensional bi-periodic frame is incorporated by constraint equations with Lagrangian multipliers and is implemented by the mortar element method. In our formulation, the bulk stress is evaluated by simple boundary integrals. Concentrating on 2D circular disk particles, we present numerical examples of single-particle, two-particle and 100-particle problems in the extensional bi-periodic frame. We discuss effects of solid fraction and particle configuration on the elongational viscosity of the suspension, in comparison with simple shear flow. We found that, at zero strain, the relative elongational viscosity is almost the same as the relative shear viscosity in simple shear for moderately concentrated suspensions. There is a small increase in elongational viscosity for large strains, which is related to an anisotropic distribution of the particles. © 2006 Elsevier B.V. All rights reserved. Keywords: Direct numerical simulation; Suspension rheology; Extensional rheology; Elongational viscosity; Extensional bi-periodic frame
1. Introduction It is quite common to add particles to fluids and polymer fluids in particular. The reasons can be quite diverse, e.g. for the manufacturing of ceramics [1] where the polymer is removed after processing, changing the product properties (such as adding flame retardants) or modifying the properties during processing. Therefore, studying the rheological properties of particle-filled polymer is very important. In most studies so far only shear deformation has been considered, whereas in processing both shear and elongational deformation components are present. Studying particle-filled polymers in elongational flow is therefore very relevant for processing. Since our final objective is to study flow of particles in viscoelastic fluids in general flows and in elongational flows in particular, we need to use direct simulation methods that model
∗
Corresponding author. Tel.: +82 55 751 6295; fax: +82 55 757 5622. E-mail addresses: [email protected] (W.R. Hwang), [email protected] (M.A. Hulsen). 0377-0257/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.jnnfm.2006.04.004
the interaction between particles by a direct simulation of the fluid in between the particles. Studying systems with many particles using the finite-element method with boundary-fitted meshes would be too complicated, especially in 3D. Therefore, we employ the fictitious domain method (see, for example Ref. [2]) in which internal boundaries are created using Lagrangian multiplier techniques on simple domains, such as a rectangular region, that can be meshed rather easily. In this paper, we present a new direct numerical simulation technique for 2D hard particle suspensions in a Newtonian fluid under planar elongational flow. The technique is based on the fictitious domain method for solid particles in fluids in bi-periodic domains under shear as developed in [3] and [4], where the domain has been fixed. Here we extend this approach to biperiodic regions that deform in time according to the average deformation in planar elongational flow. In this way we are able to obtain bulk properties for (the start-up of) elongational flows. The methods employed can relatively easy be extended to viscoelastic flows and 3D domains. The paper is structured as follows. First in Section 2 we define the problem and explain how it is modeled. Then in Section 3 we
168
W.R. Hwang, M.A. Hulsen / J. Non-Newtonian Fluid Mech. 136 (2006) 167–178
describe the numerical methods employed and the implementation of these methods. In Section 4 we show results of three examples: a flow with a single particle, two particles and 100 particles, respectively. Finally we end with some conclusions. 2. Modeling 2.1. Problem definition We consider the suspension with a Newtonian fluid consisting of a large number of circular disk particles subject to planar elongational flow. The hydrodynamic interaction between particles and fluid leads to complicated particle motion and microstructural development for both fluid and particles. In order to deal with such a problem at reasonable computational cost, a welldefined bi-periodic domain concept needs to be introduced that may transform the suspension problem with a large number of particles into a particulate flow problem in a unit cell. In simple shear flow, the Lees–Edwards boundary condition (LEbc), proposed by Lees and Edwards [5] for Molecular Dynamics, satisfies the above requirement exactly with the help of the sliding boundary concept. In the previous studies, the authors extended the LEbc concept to the continuous field problems by combining it with the finite-element/mortar-element method, for particle suspension problems [3,4]. In non-equilibrium molecular dynamics simulation of elongational flow Heyes [6] introduced the concept of the deforming simulation cell, which is aligned in flow direction and stretches in that direction. The cross-flow dimensions of the cell decrease such that the volume of the cell remains constant. A disadvantage of this approach is that only relatively small strains can be reached: the contracting dimensions eventually become of the same order of magnitude as the interaction range of a single molecule. For complex molecules an equilibrium condition (steady state) has usually not been reached yet at this point. For planar elongational flows it is possible to use so-called Kraynik–Reivelt boundary conditions [7] and simulations can be performed indefinitely [8,9]. However this approach is difficult to use and cannot be extended to uniaxial elongational or bi-axial stretching flows. In this paper we adopt the simple deforming cell approach of Heyes [6] and propose the extensional bi-periodic frame as a
computational domain for the representative unit cell problem in planar elongational flow of particle filled fluids. Fig. 1 shows the extensional bi-periodic frames with the elongation rate ˙ . (The coordinate system is indicated therein as well.) The extensional frame is a material frame in that it translates and deforms along the flow kinematics. The hatched region in Fig. 1 indicates a single frame for the illustration. At an arbitrary instance, say t = 0, the unbounded domain is regularly divided into an infinite number of cells (or frames) of the width L0 and the height H0 . As time goes on, frames elongate in the horizontal direction and contract in the vertical direction, while translating along the flow. Given the elongation rate ˙ as in Fig. 1, the time-dependent changes of the width L and the height H of the frame can be expressed as L(t) = L0 exp(˙t),
H(t) = H0 exp(−˙t).
(1)
Although each frame translates at different velocity, the frame dimensions L(t) and H(t) is the same for all the extensional frames and for all time t. A position vector x˜ measured from an arbitrary global reference frame OXY can be decomposed as the sum of the position of the frame center XF and the relative position vector x based on the reference frame OF xy (Fig. 2): x˜ = XF + x Note that the relative position vector x is periodic, in contrast to the global position vector x˜ . Similarly, the velocity field can be expressed as ˜ x) = U F (XF ) + u(x), u(˜
(2)
˜ UF and u are the velocity of the material point with where u, respect to an arbitrary global reference, the frame velocity and the relative velocity inside the frame, respectively. The frame velocity UF , the velocity at the center of the extensional biperiodic frame, with respect to the global reference is given by: U F (XF ) = (˙XF , −˙YF ).
(3)
We note that, in the absence of particles, the material points inside a frame will deform in the same way as the frame and thus the relative velocity inside the frame u(x) is the same for all the extensional frames, regardless of the global location of
Fig. 1. Extensional bi-periodic frames in planar elongational flow. The movement and the deformation of a single frame is indicated for the illustration (hatched).
W.R. Hwang, M.A. Hulsen / J. Non-Newtonian Fluid Mech. 136 (2006) 167–178
169
Fig. 4. The extensional bi-periodic frame with a four-particle configuration at a certain instance, which is the computational domain in this study.
Fig. 2. The material position x˜ -based on an arbitrary global reference can be expressed as the vector sum of the frame center coordinate XF and the relative coordinate x.
the frame, i.e. it is periodic and it is given by: u(x) = (˙x, −˙y)
(4)
Indeed, the position of the frame with respect to the global reference affects only the frame velocity UF . It is essential to realize that u(x) will be periodic also for non-homogeneous fluid systems subject to planar elongational flow, e.g. particle suspensions or droplet emulsions, if the initial configuration of the particles or the droplets is periodic from the beginning. Fig. 3 shows a set of extensional bi-periodic frames at
a certain instance. There are an infinite number of such frames in the unbounded domain. Each frame contains three circular disk-like particles and two of them cross the frame boundary. Note that the particle configuration in all frames is the same. As shown in Fig. 3, the velocity field in each frame can be decomposed into the frame velocity UF , which has been given in Eq. (3), and the relative velocity u(x). The relative velocity in the frame 1 in Fig. 3 can be obtained by specifying the global zero velocity at OF1 . In the same way, by assigning zero velocity at OF2 (or OF3 ), one gets the relative velocity in the frame 2 (or 3). Then the relative velocities in the frames 1–3 can be considered the same, since there are an infinite number of repeated structures. This means that the velocity field in each frame can be decomposed into the frame velocity UF , as in Eq. (3), and the relative velocity u(x) will be the same for all frames and be periodic. In other words, all relative events in one frame happen in the other. However, care should be exercised in determining the size of the frame: the length scale of the frame should be large enough in comparison with the scale of the microstructure. In this study, we solve a single extensional bi-periodic frame with a number of rigid particles. The frame can be periodically extended to fill the full 2D space by adding the proper frame velocity Eq. (3). The word ‘frame’ seems to be adequate for the unit cell domain, since it is a small window through which one observes what happens inside the suspension. An extensional bi-periodic frame, denoted by Ω, is the computational domain of the work and the four boundaries of the frame are denoted by Γ i (i = 1, 2, 3 and 4) (Fig. 4). We use a symbol Γ for ∪4i=1 Γi . The Cartesian x and y coordinates are selected for the direction of the stretching and that of the contraction, respectively. Particles are denoted by Pi (t) (i = 1,. . ., N) and N is the number of particles in a single frame. We use a symbol P(t) for ∪N i=1 Pi (t), a collective region occupied by particles at a certain time t. For a particle Pi , Xi = (Xi ,Yi ), Ui = (Ui ,Vi ), ωi = ωi k and i = Θi k are used for the coordinates of the particle center, the translational velocity, the angular velocity and the angular rotation, respectively; and k is the unit vector in the direction normal to the plane. 2.2. Governing sets of equations The set of equations for the fluid domain is given by:
Fig. 3. A set of the extensional bi-periodic frames containing three particles at a certain instance. Some particles cross the frame boundary.
∇ · = 0,
in Ω/P(t),
(5)
170
W.R. Hwang, M.A. Hulsen / J. Non-Newtonian Fluid Mech. 136 (2006) 167–178
∇ · u = 0,
in Ω/P(t),
= −pI + 2ηD,
(6)
in Ω/P(t),
u = U i + ωi k(x − Xi )
on ∂Pi (t),
(7) (i = 1, . . . , N).
(8)
Eqs. (5)–(8) are equations for the momentum balance, the continuity, the constitutive relation, and rigid-body conditions on particle boundaries, respectively. u, σ, p, I, D and η are the velocity, the stress, the pressure, the identity tensor, the rate of deformation tensor and the viscosity, respectively. The unknown rigid-body motions (Ui ,ωi ) in Eq. (8) will be determined by the hydrodynamic interaction that will be discussed later. In the absence of inertia, initial conditions are not necessary for the velocity of the fluid as well as the particles. In addition, there is no explicit boundary condition on the domain boundary Γ . Instead we will introduce constraint equations on Γ in Section 2.3 to assign a planar elongational flow condition. Following the work by Hwang et al. [3], we consider the circular particle as a rigid ring, which is filled with the same fluid as in the fluid domain and the rigid body condition is imposed on the particle boundary only. Thus, one needs the discretization of the particle along the boundary only. From the rigid-ring description, the set of governing equations for a region occupied by a particle Pi at a certain time t is: ∇ · = 0,
in Pi (t),
(9)
∇ · u = 0,
in Pi (t),
(10)
= −pI + 2ηD,
in Pi (t),
u = U i + ωi k × (x − Xi ),
on ∂Pi (t).
(11) (12)
Eqs. (9)–(12) are the equations for the momentum balance, the continuity, the constitutive relation, and the boundary condition, respectively, which are exactly the same as the fluid domain equations in Eqs. (5)–(8). In addition, the movement of the particle is given by the kinematic equations: dXi = U i, dt
Xi|t=0 = Xi,0 ,
(13)
dΘi = i , dt
Θi|t=0 = Θi,0 .
(14)
Eq. (14) is completely decoupled from the other equations, since the particle is circular. To determine the unknown rigid body motions (U i , i )’s of the particles, one needs balance equations for drag forces and torques on particle boundaries. In the absence of inertia and external forces or torques, particles are force-free and torquefree: Fi = · n ds = 0, (15) Ti =
∂Pi (t)
∂Pi (t)
(x − Xi ) × ( · n)ds = 0,
(16)
where Ti = Ti k and n is a normal vector on ∂Pi pointing out of the particle (i = 1,. . ., N). We did not use an artificial particleparticle collision scheme [2], because the particle overlap and
particle/wall collision could be avoided for the multiple-particle problems we studied in this paper by taking a relatively small time-step and a sufficiently refined particle boundary discretization [3]. 2.3. Extensional bi-periodic frame constraints Here we discuss the mathematical description of the biperiodicity in the extensional frame. The continuity of the velocity field and the force balance need to be satisfied across the domain boundary Γ . The conditions in the horizontal direction (between Γ 2 and Γ 4 ) at a certain time t can be written as follows: u(0.5L(t), y) = u(−0.5L(t), y) + f h , y ∈ [−0.5H(t), 0.5H(t)]
(17)
t(0.5L(t), y) = −t(−0.5L(t), y), y ∈ [−0.5H(t), 0.5H(t)],
(18)
where the vector t is the traction force on the boundary and f h = (˙L(t), 0). With the elongation rate e as shown in Fig. 2, the velocity in the x direction on the right boundary (Γ 2 ) is faster than that on the left boundary (Γ 4 ) by the amount of the velocity difference ˙ L(t). Similarly, the conditions for the periodicity in the vertical direction between Γ 1 and Γ 3 are u(x, 0.5H(t)) = u(x, −0.5H(t)) + f v , x ∈ [−0.5L(t), 0.5L(t)] t(x, 0.5H(t)) = −t(x, −0.5H(t)),
(19) x ∈ [−0.5L(t), 0.5L(t)]. (20)
The velocity difference in the vertical direction is specified by f v = (0, −˙H(t)) in Eq. (19). In the weak formulation of the finite element method, the kinematic constraint is usually combined with the Lagrangian multipliers and then the force balance can be satisfied implicitly through the multipliers. In this regard, we use the kinematic equations Eqs. (17) and (19) only in the derivation of the weak formulation and we call the two equations the extensional biperiodic frame constraints. 2.4. Boundary-crossing particles When a particle crosses the boundary of the computational domain, the particle or parts of the particle need to be relocated into the domain. The relocation of the particle always involves the change of the translational velocity of the particle in the extensional bi-periodic frame to satisfy the continuity of the velocity field. An example situation is described in Fig. 5: the part of the particle which is located outside the computational domain is relocated into the domain according to the bi-periodicity. As indicated, the relocation involves the change of the particle velocity. In order to describe the relocation, we consider two sets of the coordinates: the unprimed set for the original coordinate (before relocation) and the primed set for the relocated coordinate. The relocated coordinate x = (x ,y ) for the
W.R. Hwang, M.A. Hulsen / J. Non-Newtonian Fluid Mech. 136 (2006) 167–178
171
in the vertical direction (Eq. (19)) and the rigid-ring constraint equation along the i-th particle boundary (Eqs. (8) and (12)): h = (hx , hy ) ∈ L2 (Γ4 )2 ,
v = (vx , vy ) ∈ L2 (Γ3 )2 ,
p,i 2 2 p,i = (p,i x , y ) ∈ L (∂Pi (t)) .
Fig. 5. An illustration of a particle crossing the frame boundary: the part of the particle which is present outside the domain is relocated into the domain according to the bi-periodicity.
original coordinate x on the particle boundary is determined as follows:
In Ref. [3], we showed that the multipliers h and v are the traction force on the domain boundary Γ and that the multiplier p,i is related to the traction force acting on the particle boundary ∂Pi (t). The weak form for the problem can be stated as follows:At a certain time t with L(t), H(t) and for a given particle configuration Xi (i = 1,. . ., N), find u ∈ H1 (Ω)2 , U i ∈ R2 , ωi ∈ R, p,i ∈ L2 (∂Pi (t)), p ∈ L2 (Ω), h ∈ L2 (Γ 4 ), and v ∈ L2 (Γ 3 ) (i = 1,. . ., N) such that − p∇ · v dΩ + 2ηD[u] : D[v] dΩ Ω
x = mod(x + 0.5L, L) − 0.5L, y = mod(x + 0.5H, H) − 0.5H
+
(21)
where ‘mod’ is the modular function: e.g. mod (1.7X,X) = 0.7X and mod (−0.7X,X) = 0.3X for a positive real number X. The accompanying modification for the translational velocity of the particle associated with the relocated point x , denoted by U = (U ,V ), can be given as follows: ⎧ ⎪ ⎨ U − ˙ L(t), x > 0.5L(t), U = U + ˙ L(t), x < −0.5L(t), ⎪ ⎩ U, otherwise. ⎧ ⎪ ⎨ V + ˙ H(t), y > 0.5H(t), V = V − ˙ H(t), y < −0.5H(t), (22) ⎪ ⎩ V, otherwise. where U = (U,V) denotes the translational velocity of the particle associated with the original coordinate x. Note that the angular velocity does not change during the relocation. 3. Numerical methods 3.1. Combined weak formulation
Ω
N i=1
((p,i , v − [V i − χi k × (x − Xi )])∂Pi
+ (h , v(0.5L, y) − v(0.5L, y))Γ4 + (v , v(x, 0.5H) − v(x, −0.5H))Γ3 = 0
(23)
q∇ · u dx = 0
(24)
Ω p,i
(µ , u − [U i − ωi k × (x − Xi )])∂Pi = 0,
(i = 1, . . . , N) (25)
h
(µ , u(−0.5L, y) − u(0.5L, y))Γ4 = (µ , f h )Γ4 , v
f h = (˙L(t), 0)
(26)
(µv , u(x, 0.5H) − u(x, −0.5H)Γ3 = (µv , f v )Γ3 , f v = (0, −˙H(t))
(27)
for all v ∈ H 1 (Ω)2 , V i ∈ R2 , χi ∈ R, µp,i ∈ L2 (∂Pi (t)), q ∈ L2 (Ω), µh ∈ L2 (Γ 4 ), and µv ∈ L2 (Γ3 ). The inner product (·, ·)Γj is the standard inner product in L2 (Γ j ): (µ, v)Γj = µ · v ds. Γj
We have several remarks on the weak form Eqs. (23)–(27): Following the combined weak formulation of Glowinski et al. [2] in which the hydrodynamic force and torque acting on the particle canceled exactly, we derived the weak form for the Newtonian system together with the rigid-ring description of the particle and with the sliding bi-periodic frame in the simple shear flow in our previous study [3]. The only change in this work from the previous work is the extensional bi-periodic frame constraint (Eqs. (17) and (19)). Therefore, we present the final weak form without the detailed derivation that is available in Ref. [3]. We introduce three different kinds of Lagrangian multipliers h , v and p,i which are associated with the kinematic constraint equation for the periodicity in the horizontal direction (Eq. (17)), the constraint equation for the sliding periodicity
• The rigid-body motion on the particle boundary is treated via the rigid-ring constraint, which is satisfied in a weak sense with the multiplier λp (Eqs. (23) and (25)); the extensional bi-periodic frame constraint is satisfied again weakly with the multipliers λh and λv (Eqs. (23), (26) and (27)). • Since the extensional bi-periodic constraint (Eqs. (17) and (19)) impose only the relative relation in the velocity, one needs to specify the reference velocity at a single point in the domain. The choice is arbitrary. We specify the zero velocity at the center of the domain (x,y) = (0,0) throughout the study. • The pressure level of the fluid domain is determined by specifying one of the normal component of the Lagrangian mul-
172
W.R. Hwang, M.A. Hulsen / J. Non-Newtonian Fluid Mech. 136 (2006) 167–178
tipliers on Γ (λhx or λvy ) since the multiplier can be identified as the traction force. • For a given particle configuration, L(t) and H(t), one can solve Eqs. (23)–(27) to get all the solution and then determine the next particle configuration and the domain configuration (dimension) using Eq. (13) and Eq. (1), respectively. • When a particle crosses the domain boundary, a part of the particle needs to be relocated as described in Eqs. (21) and (22). In this case, Eqs. (23) and (25) need to be modified as follows: (µp,i (x ), u(x ) − {U i + i × (x − Xi )})∂Pi = (µp,i (x ), f )∂Pi ,
(28)
where f = (fx , fy ) and it depends on the original coordinate (before relocation): ⎧ ⎪ ⎨ −˙L(t), x > 0.5L(t), x < −0.5L(t), fx = ˙ L(t), ⎪ ⎩ 0, otherwise, ⎧ y > 0.5H(t), ⎪ ⎨ ˙ H(t), fy = −˙H(t), y < −0.5H(t), (29) ⎪ ⎩ 0, otherwise.
Fig. 6. The regular rectangular discretization is used for the entire computational domain and the particle is described by the collocation points along the particle boundary.
where A = L(t)H(t) is the area of the extensional frame that is constant. The first terms in the right-hand side of Eq. (30) are the contribution from the boundary traction (i.e. λh or λv ) and the second ones are the contribution from the boundary-crossing particles. The latter vanishes when all particles are completely immersed in the extensional frame: i.e., x = x. 3.3. Implementation
3.2. Bulk stress The bulk stress, the average stress over the domain, can be expressed, for a volume V, as the sum of the fluid contribution and the particle contribution [10]: 1 1 1 σ = σ dV = dV + σ dV, V V V Vf V Vp where · denotes the averaged quantity in V, Vf and Vp are the volume occupied by the fluid and the particle, respectively. In the weak formulation of the sliding bi-periodic frame, we have the boundary integral expression for the bulk stress using the identity between the traction force and the Lagrangian multipliers λh , λv and λp which was discussed in detail in Ref. [3]. Such a relationship is still valid in the extensional bi-periodic frame. In this regard, we present the final expression of the bulk stress using the boundary integral without detailed proof: σ11 =
1 H(t)
Γ4
A
∂Pi
i=1
N
1 λvy ds + A Γ3
1 σ22 = − L(t) 1 σ12 = L(t)
λhx ds +
N 1
i=1
Γ3
N
λhy
1 ds + A i=1
Mi k=1
µk · {u(xk ) − (U i + ωi × (xk − Xi ))}, p,i
(y − y )λp,i y (x ) ds,
(x − x )λp,i y (x ) ds,
(30)
(31)
where Mi , xk , xk and µk are the number of the collocation points on ∂Pi , the original (before relocation) of the k-th collocation point, the relocated coordinate of the k-th collocation point and the collocated Lagrangian multiplier at xk , respectively. We define equally distributed collocation points on the particle boundary and the number of the collocation points M is taken to be proportional to the radius of the particle. In our p,i
∂Pk
∂Pi
≈
(x − x )λp,i x (x ) ds,
The spatial discretization for the weak form in Eqs. (23)–(27) is basically the same as our previous work on the sliding biperiodic frame in simple shear flow [3], except for the extensional bi-periodic frame constraint equations in Eqs. (26) and (27). However, for the completeness of the paper, we briefly summarize the method. For the discretization of the momentum equation in the weak form, we use regular quadrilateral elements with continuous bi-quadratic interpolation (Q2 ) for the velocity u, discontinuous linear interpolation (P1 ) for the pressure p for the entire computational domain, including the interior of the particle. An example mesh for the two-particle problem is presented in Fig. 6. To discretize the weak form of the rigid-ring constraint (Eqs. (23), (25) or (28)), we use the point collocation method. For example, the integral in Eq. (28) has been approximated as follows: (µp,i (x ), u(x ) − {U i + ωi × (x − Xi )})∂Pi
W.R. Hwang, M.A. Hulsen / J. Non-Newtonian Fluid Mech. 136 (2006) 167–178
previous work, we analyzed effects of the number of collocation points and found that approximately one collocation point in an element appears to give the most accurate solution [3]. A small change in M from the optimal number does not affect solutions much, but it gives additional control to avoid particle collision, although it hardly occurs in our simulation with a reasonable time step size. To discretize the boundary integrals for the horizontal periodicity in Eq. (26), we use the mortar element method for the integral in Eq. (27) using the continuous linear interpolation of the λh . For the vertical periodicity, we employed the nodal collocation method, point collocation at all nodes. Since the facing elements between Γ 2 and Γ 4 or between Γ 1 and Γ 3 are always conforming due to the regular discretization of the computational domain, one could use the nodal collocation method to implement the periodicity constraint for both directions. The latter method found to generate identical results with the method used in this work. For each time step, we have the change in the particle configuration Xi (i = 1,. . ., N) and the change of the computational domain, i.e. L(t) and H(t). We remesh the computational domain to keep the aspect ratio of the element almost one. For given Xi (tn ) (i = l,. . ., N), L(tn ) and H(tn ), we performed the following procedures: • Step 1: solve the weak form at tn in Eqs. (23)–(27) to get (u, p, Ui , ωi , λp,i , λh , λv ) for the n-th step. • Step 2: update the particle configuration Xi (tn+1 ) by solving the kinematic equation of Eq. (13) with the explicit Euler method. • Step 3: update the computational domain configuration L(tn+1 ) and H(tn+1 ) using Eq. (1). In the first step, an equation with a sparse symmetric matrix with many zeroes on the diagonal appears as a result of the above discretizations, which has been solved by a direct method based on the sparse multi frontal variant of the Gaussian elimination (HSL2002/MA41) [11]. In the step 3, the remeshing of the computational domain is done by the simple scheme below: L H Nx = int , Ny = int , (32) le le
173
where Nx , Ny , le and int are the number of discretization in the x and y direction, the predetermined element size and the integer truncation operator, respectively. 4. Example problems In this section, we present three example problems to demonstrate the feasibility of our formulation and the numerical schemes: single particle, two particle, and 100 particle problems. Because of the bi-periodicity of the computational domain, the result from a single frame computation with a number of particles can be considered as the problem with the periodically extended configuration of particles in an unbounded domain. 4.1. The single particle problems The first problem is stated as follows: a single particle of radius r is suspended freely at the center of an extensional biperiodic frame with L0 = H0 = 1 under the given elongation rate ˙ = 1 in a Newtonian fluid having a viscosity η = 1. We used a 50 × 50 fluid mesh for the single particle problem, except for the r = 0.05 and 0.075 cases in which a 100 × 100 mesh is employed. The time step t is 0.02. As shown in Fig. 7, this problem represents a regular configuration of an infinite number of particles. It is an extreme case such that the particles in the neighboring frames form vertical strings as time goes on. Therefore one can expect excessively high flow resistance (elongational viscosity) when the particles approach each other (when H = 2r, the elongational viscosity is infinite). Although the problem is an artificial one, it will show that the elongation viscosity increases for an anisotropic particle distribution. This result can be used to explain the increase in viscosity for the large particle systems that will be analyzed later. Fig. 8 shows an example result for a r = 0.15 particle, plotted with the distribution of the total shear rate (II2D ). A high shear rate region is found both on the top and the bottom of the particle due to the existence of other particles in the upper and the lower frames, when the extensional strain is large (˙t = 0.8). The bulk elongational viscosity η¯ 1 can be defined as η¯ 1 =
σ11 − σ22 ˙
Fig. 7. The single particle problem in the extensional bi-periodic frame represents a regular configuration of an infinite number of particles.
(33)
174
W.R. Hwang, M.A. Hulsen / J. Non-Newtonian Fluid Mech. 136 (2006) 167–178
Fig. 8. The single particle problem result in case of r = 0.15. The plotted quantity is the distribution of the shear rate II2D , the second invariant of the rate-ofdeformation tensor 2D: ˙ t = 0, 0.4, 0.8 (from the top to the bottom).
in the planar elongational flow. For the pure Newtonian fluid in 2D, the elongational viscosity becomes four times the shear viscosity η : η¯ 1 = 4η (The factor is 3, for uniaxial elongational flow in 3D, which is known as the Trouton ratio.) In this regard, we could define the relative bulk elongational viscosity as follows: η¯ 1 η¯ 1r = . (34) 4η The elongational viscosity increases with the presence of the particle. For the dilute system where the interaction between particles is neglected, it increases linearly with the solid area fraction φ with the factor of 2 [12]. The same linear relationship is valid not only for simple shear flow but also for planar elongational flow. In Fig. 9, we present the relative bulk elongational viscosity as a function of the extensional strain = ˙ t from the single particle problem from a very dilute system with φ = 0.079 (r = 0.05) to a concentrated system with φ = 0.38 (r = 0.35). The viscosity increases with , since the vertical distances between particles decrease exponentially that leads to large flow resistance, as mentioned above. Also the viscosity increases with φ. Since the governing equations and the boundary condition are linear, one can expect a symmetric viscosity result for the negative strain < 0, i.e. contraction in the horizontal direction but stretching in the vertical direction. Thus, in the single particle problem, the elongational viscosity for zero strain = 0 reaches a minimum value, where both horizontal and vertical distances between the particle are the same. The elongational viscosity at zero strain is of interest, since the comparison between elongational and shear viscosities can be made for the same particle configuration in the same domain size L = H = 1. Fig. 10 shows the comparison of the relative bulk shear viscosity ηr and the relative elongational viscosity η¯ 1r at zero strain. The simple shear flow computation has been performed by using our previous code on particle suspensions in sliding bi-periodic frame [3].
Fig. 9. The relative elongational viscosity η¯ 1r as a function of the extensional strain = ˙ t for various solid area fraction φ.
We remark that care should be exercised to interpret the results from the single particle problem with an artificial configuration of the particle. The excessive increase of the elongational viscosity for the high φ value cases is largely due to the regular configuration of the particle.
Fig. 10. The comparison of the relative bulk elongational viscosity η¯ 1r and the relative shear viscosity ηr at zero strain as a function of the solid area fraction φ from the single particle problem. The relative shear viscosity has been adopted from Ref. [3].
W.R. Hwang, M.A. Hulsen / J. Non-Newtonian Fluid Mech. 136 (2006) 167–178
175
Fig. 11. The initial configuration of the two-particle problem in the extensional bi-periodic frame.
4.2. The two-particle problem The two-particle problem is constructed to investigate the effect of the hydrodynamic interaction between particles on the bulk suspension behavior and is stated as follows: two identical particles of radius r = 0.1 are suspended initially in an extensional bi-periodic frame of L0 = H0 = 1 under the elongation rate ˙ = 1 and the fluid viscosity η = 1. The fluid domain is discretized by a 50 × 50 mesh and the time step is 0.02. As shown in Fig. 11, the two particles are located symmetrically about the center of the domain (0, 0) in the beginning: one particle center coordinate is at (−d, 0.25) and the other is at (d, 0.25). The parameter d is the offset from the center line and it is roughly a measure of (head-to-head) interaction between two particles. When d = 0, the two particle align vertically and they will end up with head-to-head collision as the single particle problem, whereas such a interaction will be minimized for the case d = 0.25. Fig. 12 shows the particle trajectories of the two-particle problem. The trajectories for the zero offset case appears to be just straight vertical lines and the two particles approach gradually as time goes on. For the other value of d, the trajectories roughly follow the given flow field of the planar elongational flow. Fig. 13 is the bulk relative elongational viscosity η¯ 1r of the two-particle problem as a function of the extensional strain for different offset d’s. The elongational viscosity for d = 0 increases the fastest, while that for d = 0.25 increases the slowest with respect to the extensional strain. These behaviors can be understood easily, since the height of the frame H(t) that leads to the infinite viscosity is twice of the particle diameter (4r) for the zero offset case and about the particle diameter (2r) for 0 < d ≤ 0.25, leading to asymptotes near = ln 2.5 ≈ 0.92 and ln 5 ≈ 1.61. For the intermediate values of d (0 < d < 0.25, so excluding d = 0 and d = 0.25), the elongational viscosity increases for small strain but later decreases with the strain under further strain. When the two particles mostly align in the vertical direction, the viscosity increases with the strain as is the case of d = 0; but, as the particles separate from each other
Fig. 12. The particle trajectories for various d values in the two-particle problem of r = 0.1 and ˙ = 1 (the solid fraction φ = 0.0628).
along the flow, the interaction between the particle becomes less and the viscosity behavior becomes similar to that of the case d = 0.25. It is interesting that, if there is a non-zero offset, even though it is small, the bulk viscosity behaves independent of the offset value after a certain value of the strain and that the effect of the offset parameter d is found to be negligible, since the non-zero offset guarantees the exponential separation of particles. This implies that there would not be a sharp increase of the bulk viscosity due to the (vertical) alignment of particles, for the suspension with the random particle
Fig. 13. The bulk relative elongational viscosity η¯ 1r as a function of the extensional strain = ˙ t in the two-particle problem of r = 0.1 for various d values φ = 0.0628).
176
W.R. Hwang, M.A. Hulsen / J. Non-Newtonian Fluid Mech. 136 (2006) 167–178
Fig. 14. An example 100 particle problem result for r = 0.025. The plotted quantity is the distribution of the shear rate II2D : ˙ t = 0, 0.5 and 1 (from the top to the bottom).
configuration, which will be illustrated through the next example. 4.3. Many particle problems The two previous problems are more or less restricted in investigating the suspension behavior because of their regular particle configuration, in the last example problem, we attempt randomly distributed many particles in the extensional bi-periodic frame in order to observe evolution of the particle structures and bulk suspension properties. The problem is stated as follows: randomly distributed 100 particles of the same radius r are suspended initially in an extensional bi-periodic frame of L0 = H0 = 1 under the given elongation rate ˙ = 1 and the fluid viscosity η = 1. The fluid domain is discretized by a 100 × 100 mesh and the time step is 0.02. Fig. 14 shows an example result with a particle radius r = 0.025 (φ = 0.196) for the strain ˙ t = 0, 0.5 and 1 with the distribution of the total shear rate (II2D ).1 In this case, the vertical distances between particles decrease exponentially, while the horizontal distance increases, due to the planar elongational flow. However, there is no vertical alignment of particles, since randomly distributed particles and their interactions always keep the particles from being aligned: particles continuously fill in the space between the exponentially separating other particles. Therefore, one can expect a much smaller increase of the bulk elongational viscosity with respect to the applied strain in many particle problems than in the single- or two-particle problems. Plotted in Fig. 15 is the relative bulk elongational viscosity computed in the 100 particle problems of several different solid area fraction φ with respect to the extensional strain . Some results are readily comparable with the results in the single particle problem presented in Fig. 9, since they have exactly the same solid fraction. The comparison demonstrates the effect of the particle configuration on 1 The movie has been supplemented with this article in the online version, at: doi:10.1016/j.jnnfm.2006.04.004.
Fig. 15. The bulk relative elongational viscosity η¯ 1r of the 100 particle problems as a function of the extensional strain = ˙ t for several values of the solid area fraction φ.
the bulk behavior. The relative elongational viscosity η¯ 1r of the 100 particle problem increases much slower with respect to e than the single particle problem: e.g., for φ = 0.196, the viscosity goes up to 5 at ≈ 0.6 in the single particle problem in Fig. 9, while it remains around 1.6 in the 100 particle problem in Fig. 15. We see consistent increase in the elongational viscosity around the strain ≈ 1 for all cases in Fig. 15. The increase might be related to the squeezing particles in the vertical direction, like in the one and two particle problems. However here, the alignment can be avoided, since particles fill in the space generated by the exponentially separating other particles. Due to this settling of particles, the elongation viscosity does not increase excessively and, sometimes, a decrease in the viscosity even occurs afterward (φ = 0.159 and 0.126 in Fig. 15). Plotted in Fig. 16 are the average horizontal and vertical distances between two nearby particles for φ = 0.196 (r = 0.025). The horizontal distance for each particle to other particles is defined by taking the horizontal distance to the closest particle within a angular window of within ±45◦ about the x-axis. The vertical distance is defined similarly but now using ±45◦ about the y-axis. It can clearly be seen in Fig. 16 that the average horizontal distance slightly increases and the average vertical distance slightly decreases as a function of strain. This means that the particle distribution becomes slightly anisotropic, which can also be noticed in Fig. 14. We have shown for the single particle case that the viscosity increases when the particle distribution becomes anisotropic. It seems likely that this is true for the many particle case also, so we believe that the anisotropic distribution here explains the slight increase in elongational viscosity as a function of strain.
W.R. Hwang, M.A. Hulsen / J. Non-Newtonian Fluid Mech. 136 (2006) 167–178
Fig. 16. The average vertical and horizontal distances between nearby particles for φ = 0.196 in the extensional bi-periodic frame.
As we did in the single particle problem, it would be worthwhile to compare the relative elongational viscosity of planar elongational flow and the relative shear viscosity of simple shear flow for the 100-particle problem. To get the results for the simple shear, we have performed the 100 particle problem set in a sliding bi-periodic frame of the size 1 × 1 at t = 0, using the same fluid and particle meshes as the elongational flow [3]. The imposed shear rate γ˙ = 1 and the fluid viscosity η = 1. Again we emphasize that the comparison here has been made only at the initial particle configuration at zero strain at which the computational domain size and the particle configuration for both problems are the same. Therefore, this comparison does not reflect the evolutionary effects of particle configurational structures. Plotted in Fig. 17 is the comparison result of the relative bulk elongational viscosity and the relative bulk shear viscosity of the 100 particle problem, along with the single particle results for both flows that have been presented in Fig. 10. The result show that both the relative viscosities for the 100 particle problems are almost the same for 0.07 < φ < 0.25 for the given particle configuration at zero strain. Similar results has been reported experimentally in Greener and Evans [1] such that the Trouton ratio is close to the theoretical value of 3 for moderate solid volume fractions in uniaxial elongational flow. In addition, Fig. 17 shows that the prediction from the single particle problem of planar elongational flow is excessive. Small differences in the relative viscosity for elongational and shear flows are believed to be caused by the difference in hydrodynamic interaction. In Fig. 18, we present an example 100 particle problem result for r = 0.025 in simple shear flow for comparison of the hydrodynamic interaction with planar elongational flow. Again, the plotted quantity is the distribution of the shear rate II2D . Comparing Fig. 18 and the top figure in Fig. 14, one
177
Fig. 17. The comparison of the bulk relative elongational viscosity η¯ 1r as a function of the solid fraction φ for the 100 particle problems in the planar elongational flow (denoted ‘pef’) at the zero strain and in the simple shear flow (denoted ‘ssf’) along with the single particle problem results in Fig. 10.
finds that the direction of the high sheared region in between particles is different for simple shear and planar elongational flows. In simple shear flow, the high shear region is formed in the diagonal directions, whereas it is oriented in the vertical and horizontal direction in planar elongational flow. This might be responsible for the difference in the relative viscosity in Fig. 17. Further investigation on the particle microstructural effects on the bulk suspension properties is necessary but this will be part of our future work.
Fig. 18. An example 100 particle problem result for r = 0.025 in a sliding biperiodic frame for comparison with that of planar elongational flow (the top figure in Fig. 14). The plotted quantity is the distribution of the shear rate II2D .
178
W.R. Hwang, M.A. Hulsen / J. Non-Newtonian Fluid Mech. 136 (2006) 167–178
5. Conclusions
Acknowledgements
In this work, we developed a new bi-periodic domain concept, the extensional bi-periodic frame, for the planar elongational flow, which is quite well suited for studying rheology and microstructural development in non-homogeneous multiphase materials such as droplet emulsions, particle suspensions, etc. We combined the extensional bi-periodic frame with particle suspensions in 2D and developed numerical techniques for implementation using the finite-element/fictitiousdomain/mortar-element methods. Our method is able to deal with particles crossing the computational boundary and the bulk stress can be expressed by simple boundary integrals. Concentrating on disk-like particle suspensions in 2D, we studied numerical examples of single-particle, two-particle and 100-particle problems in the ex-tensional bi-periodic frame. Through the examples, we have discussed effects of the solid fraction and the particle configuration on the elongational viscosity of particle suspensions in comparison with that in simple shear flow. In the single particle problem with a regular particle configuration, we obtained the bulk elongational viscosity for a wide range of the solid fraction and reported that the elongational viscosity increases much faster than the shear viscosity with respect to the solid fraction, on account of the head-tohead collision. Also the elongational viscosity has been found to have a minimum value at the instance when both horizontal and vertical distances are the same. In the two-particle problem, we discussed the effect of the particle interaction on the bulk elongational viscosity and found that the bulk elongational viscosity behaves almost the same irrespective of the initial particle configuration after a certain amount of the strain, unless the particles are aligned exactly on the contraction axis. Finally, from the 100 particle problems, we found that the dependence of the bulk elongational viscosity on the amount of the strain is minor in comparison to the single particle problem. We also reported that the relative bulk shear and elongational viscosities are found almost the same for moderately concentrated suspensions at zero strain, by comparison with the suspension simulation with the sliding bi-periodic frame in simple shear flow. We remarked that the anisotropy in the particle microstructure is responsible for a small increase in the elongational viscosity with increasing strain.
The authors thank Mr. Michael Baltussen for his early work on this topic in his Bachelor final project. This work was supported by the Korea Research Foundation Grant (KRF-2005005-J09902) funded by the Korean government (MOEHRD) and partially by the Dutch Polymer Institute (DPI). Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.jnnfm.2006.04.004. References [1] J. Greener, J.R.G. Evans, Uniaxial elongational flow of particle-filled polymer melts, J. Rheol. 42 (1998) 697–709. [2] R. Glowinski, T.-W. Pan, T.I. Hesla, D.D. Joseph, A distributed Lagrangian multipliers/fictitious domain method for participate flows, Intern. J. Multiphase Flow 25 (1999) 755–794. [3] W.R. Hwang, M.A. Hulsen, H.E.H. Meijer, Direct simulation of particle suspensions in sliding bi-periodic frames, J. Comput. Phys. 194 (2004) 742–772. [4] W.R. Hwang, M.A. Hulsen, H.E.H. Meijer, Direct simulation of particle suspensions in a viscoelastic fluid in sliding bi-periodic frames, J. NonNewtonian Fluid Mech. 121 (2004) 15–33. [5] A.W. Lees, S.F. Edwards, The computer study of transport processes under extreme conditions, J. Phys. C: Solid State Phys. 5 (1972) 1921–1929. [6] D.M. Heyes, Molecular dynamics simulations of extensional, sheet and unidirectional flow, Chem. Phys. 98 (1985) 15–27. [7] A.M. Kraynik, D.A. Reinelt, Extensional motions of spatially periodic lattices, Intern. J. Multiphase Flow 18 (1992) 1045– 1059. [8] B.D. Todd, P.J. Daivis, Nonequilibrium molecular dynamics simulation of planar elongational flow with spatially and temporally periodic boundary conditions, Phys. Rev. Lett. 81 (1998) 1118–1121. [9] A. Baranyai, P.T. Cummings, Steady state simulation of planar elongation flow by nonequilibirum molecular dynamics, J. Chem. Phys. 110 (1999) 42–45. [10] G.K. Batchelor, The stress system in a suspension of force-free particles, J. Fluid Mech. 41 (1970) 545–570. [11] HSL2002, A collection of Fortran codes for large scale scientific computation, 2002, http://www.numerical.rl.ac.uk/hsl. [12] J.F. Brady, The Einstein viscosity correction in n dimensions, Intern. J. Multiphase Flow 10 (1984) 113–114.