Direct numerical simulations of rigid body dispersions. I. Mobility/friction tensors of assemblies of spheres John J. Molina and Ryoichi Yamamoto Citation: The Journal of Chemical Physics 139, 234105 (2013); doi: 10.1063/1.4844115 View online: http://dx.doi.org/10.1063/1.4844115 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/139/23?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Experimental study on the local turbulence modulation in a horizontal particle laden flow with rigid fence AIP Conf. Proc. 1547, 270 (2013); 10.1063/1.4816876 On the Suitability of Dispersion Models for the Optical Dispersion of Poly(Ethyl Methacrylate) AIP Conf. Proc. 1349, 621 (2011); 10.1063/1.3606011 Letter to the Editor: Wall slip in dispersion rheometry J. Rheol. 54, 1177 (2010); 10.1122/1.3495981 A non-Gaussian stochastic model to describe passive tracer dispersion and its comparison to a direct numerical simulation Phys. Fluids 16, 3006 (2004); 10.1063/1.1760770 Direct numerical simulation of particle dispersion in homogeneous turbulent shear flows Phys. Fluids 13, 3346 (2001); 10.1063/1.1405443

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.24.51.181 On: Mon, 01 Dec 2014 00:01:08

THE JOURNAL OF CHEMICAL PHYSICS 139, 234105 (2013)

Direct numerical simulations of rigid body dispersions. I. Mobility/friction tensors of assemblies of spheres John J. Molinaa) and Ryoichi Yamamotob) Department of Chemical Engineering, Kyoto University, Kyoto 615-8510, Japan

(Received 30 August 2013; accepted 18 November 2013; published online 16 December 2013) An improved formulation of the “Smoothed Profile” method is introduced to perform direct numerical simulations of arbitrary rigid body dispersions in a Newtonian host solvent. Previous implementations of the method were restricted to spherical particles, severely limiting the types of systems that could be studied. The validity of the method is carefully examined by computing the friction/mobility tensors for a wide variety of geometries and comparing them to reference values obtained from accurate solutions to the Stokes-Equation. © 2013 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4844115] I. INTRODUCTION

The dynamics of particles dispersed in a host solvent, how they react to and affect the fluid motion, is a relevant problem for science as well as engineering fields.1 It is necessary to account for the macroscopic properties of suspensions (such as the viscosity, elastic modulus, and thermal conductivity), as well as the mechanics of protein unfolding,2 the kinetics of bio-molecular reactions,3 and the tumbling motion of bacteria.4 As the particles move, they generate long-range disturbances in the fluid, which are transmitted to all other particles. Properly accounting for these hydrodynamic interactions has proven to be a very complicated task due to their nonlinear many-body nature. Several numerical methods have been developed to explicitly include the effect of hydrodynamic interactions in a suspension of particles.5–11 However, their applicability to non-Newtonian host solvents or solvents with internal degrees of freedom is not straightforward, and in some cases not possible. We have proposed an alternative direct numerical simulation (DNS) method, which we refer to as the Smooth Profile (SP) method, that simultaneously solves for the host fluid and the particles. The coupling between the two motions is achieved by introducing a smooth profile at the particle interfaces. This method is similar in spirit to the fluid particle dynamics method,12 in which the particles are modeled as a highly viscous fluid. The main benefit of our model is the ability to use a fixed cartesian grid to solve the fluid equations of motion. The SP method has been successfully used to study the diffusion, sedimentation, and rheology of colloidal dispersions in incompressible fluids.13–17 Recently it has also been extended to include self-propelled swimmers18 and compressible host solvents.19 So far, however, only spherical particles have been considered. In this paper we extend the SP method to be applicable to arbitrary rigid bodies. We show the validity of the method by computing the mobility/friction tensors for a large variety of geometric shapes. The results are compared to a) Electronic mail: [email protected] b) Electronic mail: [email protected]

0021-9606/2013/139(23)/234105/8/$30.00

numerical solutions of the Stokes equation, which are essentially exact, as well as experimental data. The agreement of the DNS results is excellent for all the cases examined here. Future papers in this series will deal with the dynamical properties of rigid body dispersions in detail, with this work we aim to introduce the basics of the model and show its validity. II. SIMULATION MODEL A. Smooth profile method for rigid bodies

We solve the dynamics of a rigid body in an incompressible Newtonian host fluid using the SP method. The basic idea behind this method is to replace the sharp boundaries between the particles and the fluid with a smooth interface. This allows us to define all field variables over the entire computational domain, and results in an efficient method to accurately resolve the many-body hydrodynamic interactions. The motion of the host fluid is governed by the incompressible Navier-Stokes equation20 (∂t + uf · ∇)uf = ρf−1 ∇ · σ f ,

(1)

∇ · uf = 0,

(2)

where uf is the fluid velocity, ρ f is the density, and σ f is the Newtonian stress tensor σ f = −pI + σ  ,

(3)

σ f = ηf [∇uf + (∇uf )t ],

(4)

with pf and ηf the pressure and viscosity of the fluid, respectively. The motion of the particles is determined by the Newton-Euler equations,21 ˙ I = VI R MI V˙ I = F I

˙ I = skew(I ) · QI , Q II · I = N I ,

(5) (6)

where RI , V I , I , and QI , are the center of mass (com) position, velocity, angular velocity, and orientation matrix, 139, 234105-1

© 2013 AIP Publishing LLC

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.24.51.181 On: Mon, 01 Dec 2014 00:01:08

234105-2

J. J. Molina and R. Yamamoto

J. Chem. Phys. 139, 234105 (2013)

respectively, of the Ith rigid body (I = 1, . . . , N). The total force and torque experienced by the particles are denoted as F I and N I , respectively, with II the moment of inertia tensor and skew(I ) the skew-symmetric angular velocity matrix ⎞ ⎛ 0 −z y ⎟ ⎜ z 0 −x ⎟ (7) skew() = ⎜ ⎠. ⎝  −y x 0 The force (and torque) on each of the particles F I = F H I + F CI + F Ext is comprised of hydrodynamic contributions I F H , particle-particle interactions (including a core potential to prevent overlap) F C , and a possible external field contribution F Ext . For the present study, we neglect thermal fluctuations, but they are easily included within the SP formalism.22 We use a Eulerian reference frame to describe the fluid variables, while the particles are tracked within a Lagrangian reference frame. The coupling between fluid and particles is obtained by defining a total velocity field u, with respect to the fluid uf and particle up velocity fields, as u(x) = (1 − φ)uf (x) + φup (x), φup (x) =



φI (x; RI , QI )[V I + I × r I ],

(8) (9)

locity field to φu∗p (x) =





. φIn+1 (x) V nI + nI × r n+1 I

Next, we compute the hydrodynamic force and torque exerted by the fluid on the particles. Assuming momentum conservation, the time integrated hydrodynamic force and torque, over a period h, are equal to the corresponding momentum exchange over the particle domain

tn +h = dx ρφIn+1 (u∗ − u∗p ), ds F H (16) I tn



tn +h

ds tn

NH I

=

  dx r n+1 × ρφIn+1 u∗ − u∗p . (17) I

From this, and any other forces on the rigid bodies, we update the velocities of the particles as tn +h

−1 n C Ext V n+1 , (18) = V + M ds F H I I + FI + FI I I tn

= nI + I−1 n+1 I · I



tn +h

tn

(∂t + u · ∇)u = ρ −1 ∇ · σ + φ f p ,

(10)

∇ · uf = 0,

(11)

with ρ = ρ f and η = ηf . The stress tensor is defined as in Eq. (3), but in terms of the total fluid velocity u. The scheme used to solve the equations of motion is the same fractional-step algorithm introduced in Ref. 24, with minor modifications needed to account for the non-spherical geometry of the particles. Let un be the velocity field at time tn = nh (h the time interval). We first solve for the advection and hydrodynamic viscous stress terms, and propagate the particle positions (orientations) using the current particle velocities tn +h ∗ n u =u + ds ∇ · [ρ −1 σ − uu], (12) tn

= RnI + Rn+1 I Qn+1 I

=

QnI

+

tn +h

tn +h

ds V I ,

(13)

C Ext ds N H . I + NI + NI

(19)

Finally, the particle rigidity is imposed on the total fluid velocity through the body force φ f p in the Navier-Stokes equation

I

where φI (x; RI , QI ) is a suitably defined SP function (0 ≤ φ I ≤ 1) that interpolates between fluid and particle domains (as described below), and r I = x − RI . The modified Navier-Stokes equation which governs the evolution of the total fluid velocity field (host fluid + particle) is given by23

(15)

I

un+1 = u∗ +



tn +h

ds φ f p ,

(20)

tn



tn +h tn

  h ds φ f p = φ n+1 un+1 − u∗ − ∇pp , p ρ

(21)

with the pressure due to the rigidity constraint obtained from the incompressibility condition ∇ · un+1 = 0. We note that H neither the hydrodynamic forces (F H I and N I ) nor the constraint forces (φ f p ) need to be computed explicitly. Only their time integrals are required, as expressed in Eqs. (16), (17), and (21), and these are easily solved by evaluating the change in momentum over the particle domains. B. Rigid body representation

For computational simplicity, we consider each particle (I) as being composed of a rigid collection of nI spherical beads (see Fig. 1), with position, velocity, and angular velocity given by Ri , V i , and i . We use upper and lowercase variables to differentiate between rigid body particles and the spherical beads used to construct them, as well as the shorthand i ∈ I to refer to the nI beads belonging to the rigid body I. The rigidity constraint on the bead velocities is given by V i = V I + I × G i ,

(22)

i = I ,

(23)

tn

ds skew(I ) · QI .

(14)

tn

Given the dependence of the profile function on the particle position and orientation, we must also update the particle ve-

where G i = Ri − RI is the distance vector from the com of the rigid body (RI ) to the bead’s com (Ri ). The rigidity constraint on the position of the beads requires that the relative distances between any two of them remain constant. Thus,

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.24.51.181 On: Mon, 01 Dec 2014 00:01:08

234105-3

J. J. Molina and R. Yamamoto

J. Chem. Phys. 139, 234105 (2013)

particle geometries, because the constituent beads are free to overlap with each other. III. STOKES DRAG

If the dispersion under consideration is such that the inertial forces are negligible compared to the viscous forces, i.e., when the particle Reynolds number is vanishingly small Re = ρU L/η  1 (U and L being the characteristic velocity and length scales), the incompressible Navier-Stokes equation (Eq. (1)) reduces to the linear Stokes equation. In the overdamped stationary limit, relevant to the systems studied here, this is given by20 η∇ 2 u = ∇p + f Ext , FIG. 1. Rigid body representation as arbitrary collection of spherical beads.

(30)

with f ext any external forces on the fluid. Due to the linear nature of this equation, the force (torque) exerted by the fluid on the particles is also linear in the velocities25, 26

the G i vectors, expressed within the reference frame of the particle  G i , are constants of motion,

F = Z · U,

(31)

 G i = QtI · G i = constant,

U = M · F,

(32)

(24)

t

where A is the matrix transpose of A. The individual positions of the beads can be directly obtained from the position and orientation of the rigid body to which they belong through Ri = RI + G i ,

(25)

Gi . G i = QI · 

(26)

We note that alternative representations of the rigid body, which maintain the smooth boundary, are possible. The advantage of the spherical-bead representation we have chosen is the ease with which the smooth profile function φ I of an arbitrary rigid body can be defined. This is relevant when evaluating the body integrals over the particle domains, as we need an efficient procedure to determine whether or not a given grid point x will contribute to the integral, i.e., φI (x) > 0. We start with the profile function for spherical particles introduced in Ref. 23, φi (x) =

h[(a + ξ/2) − ri ] , h[(a + ξ/2) − ri ] + h[ri − (a − ξ/2)]  exp (−2 /x 2 ) x ≥ 0 h(x) = , 0 x

friction tensors of assemblies of spheres.

An improved formulation of the "Smoothed Profile" method is introduced to perform direct numerical simulations of arbitrary rigid body dispersions in ...
1MB Sizes 0 Downloads 0 Views