PHYSICAL REVIEW E 88, 063006 (2013)
Dynamics of finite-symmetry and general-shaped objects under shear and shear alignment of uniaxial objects at finite temperatures Peilong Chen Department of Physics, Center for Complex Systems, National Central University, Chungli, Taiwan (Received 23 May 2013; published 4 December 2013) We prove that, for an object with a finitefold rotational symmetry (except for a twofold one) around an axis and mirror symmetries (such as a square rod or pentagonal slab, etc.), dynamics of the symmetry axis in low Reynolds number shear flow exactly follows the same form as that of a uniaxial object (e.g., a circular rod or symmetric ellipsoid) as the so-called Jeffery orbits. We use the formulation in which the dynamics of the rigid body follows first-order ordinary differential equations in time [Phys. Rev. E 84, 056309 (2011)]. Interaction between the object and the shear flow enters through a set of scalar coefficients, and the flow field does not need to be solved dynamically. Results of numerical simulations for general-shaped objects also are discussed. In the second part, Brownian dynamics of a uniaxial object is studied numerically. With D as the rotational diffusion constant, α as a parameter characterizing the aspect ratio, and γ as the shear rate, the object starts to align with the flow when the value of D/(γ α) decreases near 1. At large α (the long object limit), the results suggest much lower flow alignment when D/(γ α) > 1. DOI: 10.1103/PhysRevE.88.063006
PACS number(s): 47.15.G−, 47.57.Qk, 82.70.Dd
I. INTRODUCTION
Dynamics of rigid or deformable bodies in flow fields have been important subjects in fluid dynamics, soft matter, and biophysics. In solutions, bulk properties often depend crucially on how the solutes (microscopic particles or macromolecules) respond to the flow, owing to their distinct particle shapes. Vast varieties of particles have recently been developed in many materials [1]. Microfluidic systems [2] also, very often, incorporate solid objects or macromolecules in the working fluids. In living organisms, many of these flow-response behaviors are essential for the functions of life processes. Some examples are how the various blood cells move with the blood flows in different blood vessels [3] and how bacteria, with their huge variety of body shapes, move under shear flow [4,5]. Many objects in soft matter and biophysics are in small scales such that the relevant Reynolds numbers are much less than 1. The fluid is in the creeping flow regime where the Stokes equation applies [6]. Although these objects can encounter vastly different and complicated flows, as they move with the local mean flow, what they first experience is the local velocity gradient and most notably the shear flows. Thus, one of the fundamental properties is how objects respond to shear flows, which naturally arises in complicated flow patterns or in confined geometries. There have been many experimental and theoretical studies on how macromolecules or colloids respond to a shear flow [7] or even a viscoelastic fluid [8]. However, the first and possibly the only complete analytic result for dynamics of a rigid body inside a sheared Stokes flow is the Jeffery orbits [9] for uniaxial objects (e.g., circular rods or symmetric ellipsoids). For other, “less-symmetric” shapes, numerical simulations are mostly the only way to obtain particle dynamics. Usually, the Oseen tensor and surface integral formulation [10] were used. Such simulations involve, at each time, a first step using the surface integral (on the object’s surface) to numerically calculate the flow field by simultaneously solving for the stress distribution on the object surface as well as linear and rotational velocities of the object, typically with an iterative procedure. 1539-3755/2013/88(6)/063006(7)
In the second step, the computed object velocities are used to advance the position and orientation of the object a small time step forward. By repeating the steps, the time evolution is obtained. Because the second step numerically is easy, especially for nondeformable objects, most of the computing time is spent on the first step of solving the flow field. In this paper, using a formulation previously detailed [11] (although it mostly dealt with analytic results for chiral objects), we will demonstrate that the time-consuming first step of solving the flow field can be avoided. (So simulations can be performed very efficiently.) We then analytically prove a surprising result that, for a rod with a cross section being a regular polygon, the dynamics of the central axis is independent of the orientation of the other two axes. The symmetry axis moves in exactly the same way a uniaxial object as the Jeffery orbits. For general-shaped objects, we will discuss their dynamical patterns using results obtained from simulations. By adding random effects due to finite temperatures, we can quickly obtain statistical properties of the Brownian dynamics. Specifically, we will give detailed results of shear alignment for a uniaxial object. Transition to flow alignment is found to be determined by a particular combination of the shear rate, the rotational diffusion constant, and an object geometry parameter.
II. STOKES FLOW FORMULATION FOR AN OBJECT UNDER SHEAR
In this section, we summarize the formulation, which has been derived previously [11], to study dynamics of a rigid body under shear, and as a special case, the resulting so-called Jeffery-orbit solutions for uniaxial objects. Consider an object moving through a stationary fluid with the Reynolds number Re 1 such that the momentum is negligible and the Stokes equation [6] applies. Denote the instantaneous linear and angular velocities of the object as v and ω, respectively. The induced flow field u1 (r) is zero at
063006-1
©2013 American Physical Society
PEILONG CHEN
PHYSICAL REVIEW E 88, 063006 (2013)
infinity and matches no-slip conditions on the (moving) surface of the object. The total force f 1 and torque τ 1 experienced by the object from the fluid can be written as, due to the linearity of the Stokes equation,
on inversion, the resistance tensor b is zero. Orientation of the object is specified completely by a unit vector nˆ along the symmetry axis. By symmetry, the tensors a and c are in the form of (with I as the second-rank unit tensor)
f 1 = a · v + b · ω, τ 1 = b · v + c · ω.
ˆ c = −c1 nˆ nˆ − c2 (I − nˆ n).
Here (a,b,c) are the second-rank resistance tensors depending on the shape and orientation of the object. The two cross term b tensors are the same due to the Lorentz reciprocal theorem [6]. Next we consider the same object fixed, both in space and in orientation, under an external simple shear. The shear induces a flow u2 (r), which is zero on the object surface (which is fixed) and matches the simple shear flow at infinity. Again the total force f 2 and torque τ 2 experienced by the object from the fluid can be written as f 2 = g:e, τ 2 = h:e. Here g and h are the third-rank shear-response tensors also depending on the shape and orientation. e ≡ ∇us is the second-rank velocity gradient tensor of the simple shear flow us (for example, γ y xˆ and γ is the shear rate). When a free object is subjected to an external shear flow and moves with v and ω due to the shear, u1 + u2 will be the resulting flow field. The combination is the solution of the Stokes equation due to the linearity and has the correct no-slip boundary condition on the object’s surface and simple shear at infinity. Of course, v and ω of the object are still unknown. Now because the momentum is negligible in the zero Reynolds number limit, total force and torque on the object are zero, a · v + b · ω + g:e + f r = 0,
(1)
b · v + c · ω + h:e + τ r = 0,
(2)
which determine v and ω. Here we also add the random force and torque terms to model the effects of finite temperatures. Since the resistance and shear-response tensors are orientation dependent, so are the resulting v and ω. The tensors (a,b,c,g,h) depend on the shape and orientation of the object. For a rigid body, if we define a body coordinate on the object (nˆ 1 ,nˆ 2 ,nˆ 3 ), then the second-rank tensors (a,b,c) can be written in the form of a=
3
aαβ nˆ α nˆ β .
αβ=1
Meanwhile, the third-rank tensors g and h are g=
3
gαβγ nˆ α nˆ β nˆ γ .
(3)
αβγ =1
Given a specific object shape, all these scalar coefficients can be determined by solving the Stokes equation. There are well established numerical methods, which usually are based on the Oseen tensors and surface integrals [10]. Once these coefficients are obtained, the time evolution of the position and orientation of the objects can be obtained by calculating v and ω of the object using Eqs. (1) and (2) without the need to solve the flow field at each time step. For a uniaxial object, such as a circular rod, which is rotationally symmetric with respect to an axis and symmetric
(4)
For shear response tensors, g = 0 and hij k
3 = (h1 ilk nj nl + h2 lj k ni nl + h3 ij l nk nl ).
(5)
l=1
Physically, h1 , h2 , and h3 are the torque experienced by the object under shear when the symmetry axis nˆ is fixed and pointing in the velocity, velocity-gradient, and vorticity directions, respectively. (In this paper, in components, the inner products between tensors are defined, such as (a · h)ij k = l ali hj kl and (h:e)k = l hij k ej i as in Ref. [11].) Now the rotation ω can be obtained from Eq. (2), and the dynamics of nˆ is given by d nˆ = ω × nˆ = −(c−1 · h:e) × nˆ dt ˆ × n. ˆ = γ (h1 /c2 )nx (nˆ × ˆy) × nˆ + γ (h2 /c2 )ny ( xˆ × n) (6) For the last equation, we have assumed, without the loss of generality, the external simple shear flow us = −γ y xˆ . Straightforward algebra would lead to the second line, Eq. (6). So the dynamics is determined by the ratio between h1 and h2 , characterizing the geometry of the object. It is well known that Eq. (6) has closed curve solutions when nˆ is plotted on the surface of a unit sphere as the socalled Jeffery orbits [9] (e.g., the black curves in Fig. 3) with α 2 ≡ h2 / h1 , tan2 θ (α 2 sin2 φ + cos2 φ) = const.
(7)
Here θ and φ are the polar and azimuthal angles of the spherical coordinate with the z axis along the vorticity direction of the shear flow and x in the flow direction. III. DYNAMICS OF GENERAL-SHAPED OBJECTS
In this section, we first present a surprising result for objects with a finitefold rotational symmetry (except for the twofold one) around an axis and mirror symmetries. Typical examples are rods with cross sections of regular polygons. It will be shown that, for such objects besides b always being zero, c and h are again in the forms of Eqs. (4) and (5), respectively. The rotational dynamics of the symmetry axis is, hence, also like Eq. (6) and follows the Jeffery orbits as the uniaxial objects. Two examples of such objects are illustrated in Figs. 1(a) and 1(b) with (nˆ 1 ,nˆ 2 ,nˆ 3 ) as defined. The symmetry axis nˆ 3 will also be denoted as n. Three symmetry operations, which leave these objects unchanged, will be used for the proof. The first is a π rotation around the bisection line nˆ 1 + nˆ 2 such that nˆ 1 ↔ nˆ 2 ,
nˆ 3 → −nˆ 3 .
(8)
The second is the mirror reflection of the plane spanned by the bisection line nˆ 1 + nˆ 2 and nˆ 3 ,
063006-2
nˆ 1 ↔ nˆ 2 ,
nˆ 3 → nˆ 3 .
(9)
DYNAMICS OF FINITE-SYMMETRY AND GENERAL- . . .
PHYSICAL REVIEW E 88, 063006 (2013)
n1
So we only have h123 nˆ 1 nˆ 2 nˆ 3 + permutations, and Eq. (8) requires that
(a)
(b )
n2 n1
n
h = h123 (nˆ 1 nˆ 2 − nˆ 2 nˆ 1 )nˆ 3 + h312 nˆ 3 (nˆ 1 nˆ 2 − nˆ 2 nˆ 1 ) + h132 (nˆ 1 nˆ 3 nˆ 2 − nˆ 2 nˆ 3 nˆ 1 ). (12)
n
(c)
(d ) 1
1
0.5
0.5
z
n2
z
0
0 -0.5
-0.5
1
1 -1
0.5 -1
-0.5
x
0 0
0.5
y
-1
0.5 -1
-0.5 1-1
-0.5
x
0 0
0.5
y
-0.5 1-1
FIG. 1. (Color online) (a) A rod with a triangular cross section. (b) A slab with a pentagonal cross section. (c) One example of the dynamics of nˆ for (b) as a Jeffery orbit. For clarity, other Jeffery orbits obtained with different initial conditions are not shown. (d) The corresponding dynamics of nˆ 1 to (c).
The third is a π rotation about nˆ 1 , nˆ 1 → nˆ 1 ,
nˆ 2 → −nˆ 2 + a nˆ 1 ,
nˆ 3 → −nˆ 3 ,
(10)
with a = 2nˆ 1 · nˆ 2 = 2 cos θ12 and θ12 as the angle between nˆ 1 and nˆ 2 . All the tensors (a,b,c,g,h) should be invariant under Eqs. (8) and (10). However, for reflection Eq. (9), we need to recognize that, although (a,c,g) are tensors, (b,h) are pseudotensors. These can be seen from Eqs. (1) and (2) that (v, f r ,e) are vectors and tensors (also called polar vectors for v and f r ) and (ω,τ r ) are pseudovectors (also called axial vectors). A pseudovector or tensor will acquire an extra negative sign under reflection [12]. When pseudotensors are expanded in terms of (nˆ 1 ,nˆ 2 ,nˆ 3 ), combining the requirements of Eqs. (8) and (9) plus a negative sign, we get the following requirement: b and h have no terms without nˆ 3 . (11) Thus, for b = − 3αβ=1 bαβ nˆ α nˆ β being a pseudotensor, the term nˆ 3 nˆ 3 is excluded by the requirement of Eq. (9) plus a negative sign, (nˆ 1 + nˆ 2 )nˆ 3 by Eq. (8), and 2αβ=1 nˆ α nˆ β by Eq. (11). We are left with (nˆ 1 − nˆ 2 )nˆ 3 , which becomes [(a − 1)nˆ 1 − nˆ 2 ]nˆ 3 under Eq. (10). However, a cannot be 2 since nˆ 1 and nˆ 2 are not parallel to each other. So b = 0. The rotation ω under deterministic dynamics (no random motion) is, thus, determined by c and h in Eq. (2). For h = 3αβγ =1 hαβγ nˆ α nˆ β nˆ γ being a pseudotensor, nˆ 3 nˆ 3 nˆ 3 is excluded by Eq. (10), (nˆ 1 + nˆ 2 )nˆ 3 nˆ 3 is excluded by Eq. (9) plus a negative sign, (nˆ 1 − nˆ 2 )nˆ 3 nˆ 3 is excluded 2 ˆ α nˆ β nˆ γ is excluded by Eq. (11). by Eq. (8), and αβγ =1 n The remaining terms have only one nˆ 3 each. Of these, (nˆ 1 nˆ 1 + nˆ 2 nˆ 2 )nˆ 3 is excluded by Eq. (8). For (nˆ 1 nˆ 1 − nˆ 2 nˆ 2 )nˆ 3 , after Eq. (10), we get a positive sign for nˆ 2 nˆ 2 nˆ 3 contradicting the original, and this term cannot appear from the remaining terms h123 nˆ 1 nˆ 2 nˆ 3 + permutations under Eq. (10). Thus, (nˆ 1 nˆ 1 − nˆ 2 nˆ 2 )nˆ 3 also is excluded.
Now because nˆ 1 × nˆ 2 = sin θ12 nˆ 3 , the h tensor is in exactly the same form as that of a uniaxial object Eq. (5), which only consists of nˆ 3 . (Actually, using a similar analysis, we can get g = 0. Nevertheless, this is not needed for our results as b = 0.) For c = − 3αβ=1 cαβ nˆ α nˆ β being a (regular) tensor, the term (nˆ 1 + nˆ 2 )nˆ 3 is excluded by Eq. (8) and (nˆ 1 − nˆ 2 )nˆ 3 by Eq. (9). So, c = −c11 (nˆ 1 nˆ 1 + nˆ 2 nˆ 2 ) − c12 (nˆ 1 nˆ 2 + nˆ 2 nˆ 1 ) − c33 nˆ 3 nˆ 3 . (13) However, when the third symmetry operation of rotating around nˆ 1 Eq. (10) is applied, we get c = −[(a 2 + 1)c11 + 2ac12 ]nˆ 1 nˆ 1 − c11 nˆ 2 nˆ 2 − (−ac11 − c12 )(nˆ 1 nˆ 2 + nˆ 2 nˆ 1 ) − c33 nˆ 3 nˆ 3 .
(14)
For Eq. (14) to be equal to Eq. (13), the coefficient c12 needs to be proportional to c11 as c12 = −(a/2)c11 . (15) √ Now √ by defining nˆ 1 = (nˆ 1 + nˆ 2 )/ 2 + a and nˆ 2 = (nˆ 1 − nˆ 2 )/ 2 − a [with (nˆ 1 ,nˆ 2 ,nˆ 3 ) as an orthogonal coordinate], we can diagonalize (nˆ 1 ,nˆ 2 ) terms in c. Furthermore, due to Eq. (15), both nˆ 1 nˆ 1 and nˆ 2 nˆ 2 have the same coefficient ca = c11 (2 + a)(2 − a)/4, and c is again in the form of Eq. (4), c = −ca (nˆ 1 nˆ 1 + nˆ 2 nˆ 2 ) − c33 nˆ 3 nˆ 3 = −ca (I − nˆ 3 nˆ 3 ) − c33 nˆ 3 nˆ 3 . With h and c in the same forms as in the case of uniaxial objects, dynamics of the symmetry axis nˆ = nˆ 3 would again follow Eq. (6), d nˆ = ω × nˆ = −(c−1 · h:e) × nˆ dt ˆ × n. ˆ = γ (h132 /ca )nx (nˆ × ˆy) × nˆ + γ (h312 /ca )ny ( xˆ × n) (16) So we finally conclude that the dynamics of nˆ will follow Jeffery orbits, regardless of the orientation of nˆ 1 –nˆ 2 relative ˆ to n. We also perform numerical simulations for such objects from randomly chosen initial conditions to check the dynamics of nˆ 1 and nˆ 2 by solving Eqs. (1) and (2) without the random force and torque terms. We construct the object by placing a small sphere on each of the vortices of Figs. 1(a) and 1(b) [six in total for (a) and ten for (b)]. The hydrodynamic interaction from the Oseen tensor between two well-separated spheres then is used to calculate the tensors (a,b,c,g,h). This approximation, although affecting the exact numerical values of the components of these tensors, does not change the symmetry properties of the tensors. The tensors are checked to conform to the forms obtained in the symmetry analysis above.
063006-3
PEILONG CHEN
PHYSICAL REVIEW E 88, 063006 (2013)
(a)
(b )
1
z 0
x
0
1 -1
x (f)
1
0
0
1 -1
x (g)
1
x
0
1 -1
x
0
x (h)
1
0
1 0 y 1 -1
1
z 0
1 -1 -1 0 y 1 -1
0
1 -1 -1 0 y
z 0
-1 -1
0
1 -1 -1 0 y
z
1
z 0
1 -1 -1 0 y
z
(d )
1
z 0
-1 -1
(e)
(c)
1
z
0
1 -1 -1 0 y 1 -1
x
0
1 -1 -1 0 y 1 -1
x
0
1 0 y 1 -1
FIG. 2. (Color online) (a)–(f) Dynamics of the central axis of a long rod with a rectangular cross section. (g)–(h) For a long rod with a triangular cross section. See the text for detailed geometries.
In Fig. 1(c), we plot one example of the dynamics of nˆ for the pentagonal object, which nicely follows one of the Jeffery orbits, as predicted from previous analysis. Different initial conditions, indeed, lead to different Jeffery orbits. In the meantime, the trajectories of nˆ 1 and nˆ 2 follow complex patterns with the corresponding dynamics of nˆ 1 to Fig. 1(c) shown in Fig. 1(d). Next we discuss the dynamics of general-shaped objects. Because of the infinite possibilities of shapes, we first consider two limiting categories: One dimension of the object is as follows: (1) much larger than the other two, i.e., a long rod and (2) much smaller than the others, i.e., a plate. For objects of category (1), Figs. 2(a)–2(f) show the dynamics of a rectangular cuboid with its three edge lengths denoted as L, W , and H . We use the same method to construct the objects by placing spheres on the vertices of the geometry shape and calculate the tensors (approximately) as above. With L = 200, W = 12, and H = 10, the dynamics nˆ L (t) are shown with six different initial conditions on (nˆ L ,nˆ W ,nˆ H ). We find that the trajectories mostly are similar to the Jeffery orbits when nˆ L is near the vorticity direction as illustrated in Figs. 2(a) and 2(b). When nˆ L is closer to the velocity (v)– velocity gradient (∇v) plane, its trajectory is no longer a closed Jeffery-like orbit and gradually covers a larger extended region as in Figs. 2(c) and 2(d). Near the v-∇v plane, the covered region usually is the largest as shown in Fig. 2(e). However, among the complex extended trajectories, we often see closed trajectories as in Fig. 2(f) with various degrees of complexity. We find that these general behaviors are rather representative of objects resembling a long rod. For rectangular cuboids, when the ratio W/H is closer to 1, we see a larger range of Jeffery-orbit-like trajectories. At W/H = 1, we are back to the case in Fig. 1 with all Jeffery orbits for nˆ L . Meanwhile, if the cross section of the long rod is less symmetric, for example, a triangle with edge lengths of 13, 9, and 7, not only a smaller region with simple Jeffery-orbit-like trajectories as in Fig. 2(g) with a rod length of 100 is there, but also most of the complex trajectories cover a large portion of the whole sphere. One example is shown in Fig. 2(h). For objects of category (2) as thin plates, the dynamics of nˆ H along the shortest axis are very similar to those of a long
rod as discussed. Again when L/W (for rectangular plates) approaches unity, there is a larger range of Jeffery-like orbits. Dynamics for less symmetric cross sections also follow the same pattern for the long rod. Besides the two limiting shapes, for objects with all three dimensions being comparable, we find that there are much less Jeffery-orbit-like trajectories. Furthermore, most of the trajectories could cover a large portion of the sphere surface similar to Fig. 2(h). IV. BROWNIAN DYNAMICS OF A UNIAXIAL OBJECT
In physical situations, especially in microscopic length scales, the orientation of an object will be perturbed by random forces due to finite temperatures and, hence, will not follow the Jeffery orbits exactly, even for a uniaxial object. The random effect typically is quantified by the P´eclet number Pe ≡ γ /D as the ratio between the shear rate γ and the rotational diffusion constant D, which is proportional to the temperature. For a uniform circular rod as an example, in the limit of Pe 1 (high temperatures), shear effects will be small, and the rod will orient more or less uniformly in all directions. At lower temperatures, although the rod would follow Jeffery orbits closely, the random force provides a transition mechanism between different Jeffery orbits. Thus, with random forces, there could be a steady-state distribution on the unit sphere at long times. The time needed to reach the steady-state distribution will increase with decreasing temperature. To find these steady-state distributions at finite temperatures, we perform numerical simulations of Eq. (6) with an ˆ This random rotation additional random rotation term ωr × n. always is chosen on the plane perpendicular to the direction of the symmetry axis at each time since the component in the direction of the symmetry axis has no effect on the orientation. In this way, the magnitude of ωr is assumed to follow a Gaussian distribution with the variance |ωr |2 equal to the diffusion constant D. The temperature T and D follow the Einstein relation D = BkB T /ζ with kB as the Boltzmann constant and ζ as the fluid viscosity. B is a factor depending on the exact object geometry, in general, with a highly nontrivial dependence. For our discussions, D is used as the control
063006-4
DYNAMICS OF FINITE-SYMMETRY AND GENERAL- . . .
PHYSICAL REVIEW E 88, 063006 (2013)
1
1
α2 = 3 10 30 100 300 1000
0
λ
Vorticity
0.5
0.5
-0.5
-1 0 -1
0 -0.5 0 Flow
0.5
-0.5
1 1
0 0.5 Velocity gradient
-1
Average density (arb. units)
parameter. We will relate the results to T later for a long circular rod. For α 2 = 10, D = 0.1, and γ = 1, the dynamics of nˆ is shown in Fig. 3 by plotting blue points (the upper part) on the ˆ i ), i = 1 · · · N unit sphere representing the orientation of n(t with ti = ts + i t, t as a chosen suitable time interval, and ts as a transient time interval such that the results do not depend on the initial condition. The Jeffery orbits of Eq. (7) at this value of α also are drawn as the black curves, although for clarity, only those on the upper semisphere are drawn. The concentration of blue dots on the two polar regions (of the flow directions) shows the preference of the rod to align near the flow direction under shear. The alignment can further be seen in the projection of the blue points onto the v-∇v plane (the red points on the lower part). To perform a quantitative comparison between different values of D, in Fig. 4, we plot the average densities as functions of the polar angle θf with respect to the flow direction. The two concentrated regions in Fig. 3 correspond to θf = 0 and π at which the density has peak magnitudes as in Fig. 4. Higher peak values and narrower widths indicate stronger alignment in the flow direction. In Fig. 4, we see the degree of alignment increasing with lower diffusion constant D. Furthermore, the distributions have a low-D limit for D below 0.2. This low-D limit of f (r) can actually be calculated using a lowtemperature expansion of the corresponding Fokker-Planck
10
2
FIG. 5. (Color online) The alignment index λ as a function of D/(γ α).
FIG. 3. (Color online) Brownian dynamics of nˆ showing shear alignment for α 2 = 10, D = 0.1, and γ = 1. See the text for details.
15
1 D/(γ α)
D = 0.02 0.05 0.1 0.2 0.5 1 2 4
5
0
equation for f (r) [13] (or see the Appendix). The contour lines of the calculated f (r) are plotted in Fig. 3 as the red (gray) contour curves (the mostly vertical ones). They are seen to match the distribution of the blue points from simulations well. A more quantitative comparison is seen by drawing the average density as a function of the polar angle in the flow direction in Fig. 4 as the red solid line, which perfectly agrees with the simulation results. For reference, our formulations of the low-T expansion of the Fokker-Planck equation are detailed in the Appendix, which are the same, in principle, with those in Ref. [13]. The peaks in Fig. 4 signify the degree of flow alignment of the rod under the competition between the random motion and the shear flow. To quantify the alignment, we calculate an alignment index λ, which is the portion of distribution above the “background” value at θf = π/2, λ = [f (θf ,φf ) − f (θf = π/2,φf ) φf ]d f , with φf denoting the average over φf and f denoting the solid angle for (θf ,φf ). A uniform distribution without flow alignment has λ = 0, and λ = 1 is the maximum. The results are shown in Fig. 5. It shows the object only starts to have significant alignment with the flow when D/(γ α) is near 1. Furthermore, the data suggest that, at the limit of long rods (α → ∞), the flow alignment remains quite insignificant when D/(γ α) > 1. More accurate numerical calculations could be needed to say whether λ → 0 at α → ∞ and D/(γ α) > 1. On the other hand, the maximum flow alignment λm at D → 0 significantly depends on α. We find that they closely follow a scaling relation 1 − λm ≈ α −0.43 . As an example for long circular rods, it is known that the rotational diffusion constant D goes mostly as ∼ T /L3 plus another correction term involving ln(L/a) [14] with L and a being the length and radius of the rod, respectively. On the other hand, for a rod in Stokes flow, we expect h1 ∼ L2 and h2 ∼ L, thus, α 2 ∼ L. So Fig. 5 tells us that the transition to significant flow alignment, at least for long rods, is determined approximately by T L−3.5 .
0 0.5 1 Polar angle in flow direction (units of π)
V. CONCLUSIONS
FIG. 4. (Color online) Average density as a function of the polar angle in the flow direction for α 2 = 10 and γ = 1.
To summarize, we prove that the central axis of a rod with a cross section of a regular polygon moves exactly in the same
063006-5
PEILONG CHEN
PHYSICAL REVIEW E 88, 063006 (2013)
patterns with those of circular rods or symmetric ellipsoids under external shear flow. Using our formulation, temporal dynamics can be simulated very efficiently for general-shaped objects, which possess different degrees of Jeffery-orbit-like trajectories near the vorticity direction. For uniaxial objects, we obtain shear alignment distributions at finite temperatures depending on the object shapes. Transition to flow alignment is found to mostly be controlled by the combination D/(γ α), and the maximum flow alignment at low temperature scales with α. ACKNOWLEDGMENT
The support of the National Science Council of Taiwan is acknowledged.
with k(u) to be determined. On the order of , we have ∂ hu ∂f0 ∂ hv ∂f0 ∂ hu Kf1 = + . ∂v ∂u hu ∂u ∂v hv ∂v Now if we integrate the equation by requirement of periodicity on v, we have ∂ ∂u
In this Appendix, we outline the low-temperature expansion of the Fokker-Planck equation [13]. It is well established that the Brownian dynamics of nˆ as a trajectory on the unit sphere has a corresponding Fokker-Planck equation [15] for the distribution function f (r) with r being the coordinate on the surface of the sphere, ∂f = −∇ · (K f ) + ∇ 2 f. ∂t The drift velocity K is the deterministic velocity induced by the shear and corresponds to the right hand side of Eq. (6). The second term is the diffusion of f (r), and the diffusion constant is proportional to the temperature. We want to solve for steady-state solutions: ∂f/∂t = 0 at 1. The distribution is expanded in the orders of , f = f0 + f1 + f2 + · · · . 2
Collecting terms with the same order of , we obtain ∇ · (K f0 ) = 0,
∇ · (K fi ) = ∇ 2 fi−1 , i 1.
It is noted that, with the small term having the highest order of derivative, an expansion perturbation often becomes singular [16], leading to, for example, boundary layers. However, in the current case, the space on the unit sphere is finite and has no boundaries. There also are no sink or source terms. So the perturbation could remain regular. Assume (u,v) is a curvilinear coordinate on the surface of the sphere. Note that the drift velocity K is tangent to the Jeffery orbits everywhere. Thus, if we choose u as the Jeffery orbits, u = tan2 θ (α 2 sin2 φ + cos2 φ), then K = K(u,v)ˆv . The coordinate v traces around the closed Jeffery orbits and can always be defined as 0 v 2π with v = 0 and v = 2π being the same. All functions are periodic in v. We have, on the order of 0 with hu and hv , the scale factors of the (u,v) coordinate,
dv 0
dv, from the
hv ∂f0 = 0. hu ∂u
(A2)
dk(u) + B(u)k(u) = C1 , du
with
2π
A(u) ≡ 0
hv dv G, B(u) ≡ hu
2π
dv 0
hv ∂G . hu ∂u
(A3)
C1 equals the integration in Eq. (A2) and, physically, is exactly the total flux across a contour line of a fixed u (Jeffery orbits), the same for all values of u for a steady-state distribution. Because the space is finite and there are no sink or source contributions, the constant flux is necessarily zero. So we have the lowest-order solution f0 (u,v) = [hu (u,v)K(u,v)]−1 k(u) with B(u) du . (A4) k(u) = C2 exp − A(u) C2 is determined by the normalization condition for f0 . Because we have an analytic form for the coordinate v as detailed below, we have an analytic expression for the integrands in Eq. (A3). We, thus, can obtain k(u) by numerically integrating Eq. (A3) to get A(u) and B(u) and then using them in Eq. (A4). So we choose the u coordinate of (u,v) as the Jeffery orbits, u(θ,φ) = tan2 θ (α 2 sin2 φ + cos2 φ) ≡ a(θ )b(φ). Here θ and φ are the angular parts of the spherical coordinate. To find v(θ,φ), we require 0 = ∇u · ∇v =
1 ∂u ∂v ∂u ∂v + . ∂θ ∂θ sin2 θ ∂φ ∂φ
Assuming v(θ,φ) = f (θ )g(φ) with some straightforward algebra, we obtain with some particular choices of integration constants (the choices will not affect the final calculated distribution function), u(θ,φ) = tan2 θ (α 2 sin2 φ + cos2 φ),
1 ∂ hu Kf0 = 0. hu hv ∂v
v(θ,φ) = sinα
Thus, hu Kf0 is a function of u only, f0 (u,v) = (hu K)−1 k(u) ≡ G(u,v)k(u),
2π
0
The integration in the above equation will equal a constant C1 . Thus, in the order, we obtain an equation for f0 (r) to determine k(u). Using Eq. (A1), we get A(u)
APPENDIX: LOW-TEMPERATURE SOLUTIONS OF THE FOKKER-PLANCK EQUATION FOR A UNIAXIAL OBJECT
2π
(A1)
2
−1
2
θ cosα φ/ sin φ.
The scale factors hu and hv can also be found as functions of θ and φ.
063006-6
DYNAMICS OF FINITE-SYMMETRY AND GENERAL- . . .
PHYSICAL REVIEW E 88, 063006 (2013)
[1] See, e.g., B. S. John, A. Stroock, and F. A. Escobedo, J. Chem. Phys. 120, 9383 (2004); C. J. Murphy et al., J. Phys. Chem. B 109, 13857 (2005); L. Rossi et al., Soft Matter 7, 4139 (2011); Y. Wang et al., Nature (London) 491, 51 (2012). [2] Microfluids, edited by S. Colin (Wiley, London, 2010); B. J. Kirby, Micro- and Nanoscale Fluid Mechanics: Transport in Microfluidic Devices (Cambridge University Press, Cambridge, UK, 2010). [3] R. D. Kamm, Annu. Rev. Fluid Mech. 34, 211 (2002). [4] J. T. Locsei and T. J. Pedley, Bull. Math. Biol. 71, 1089 (2009). [5] T. Ahmed, T. S. Shimizu, and R. Stocker, Integr. Biol. 2, 604 (2010). [6] J. Happel and H. Brenner, Low Reynolds Number Hydrodynamics with Special Applications to Particulate Media (PrenticeHall, Upper Saddle River, 1965).
[7] See, e.g., Marcos, H. C. Fu, T. R. Powers, and R. Stocker, Phys. Rev. Lett. 102, 158103 (2009); W. Tang and S. G. Advani, J. Chem. Phys. 126, 144711 (2007), and references therein. [8] Z. Fan and S. G. Advani, Polymer 46, 5232 (2005). [9] G. B. Jeffery, Proc. R. Soc. London, Ser. A 102, 161 (1922). [10] See, e.g., M. Makino and M. Doi, Phys. Fluids 17, 103605 (2005). [11] P. Chen and Q. Zhang, Phys. Rev. E 84, 056309 (2011). [12] G. B. Arfken, Mathematical Methods for Physicists, 3rd ed. (Academic Press, Orlando, FL, 1985), pp. 128–137. [13] L. G. Leal and E. J. Hinch, J. Fluid Mech. 46, 685 (1971). [14] S. J. Broersma, Chem. Phys. 32, 1626 (1960); 32, 1632 (1960). [15] R. Zwanzig, Nonequilibrium Statistical Mechanics (Oxford University Press, New York, 2001). [16] E. J. Hinch, Perturbation Methods (Cambridge University Press, Cambridge, UK, 1991).
063006-7
Dynamics of finite-symmetry and general-shaped objects under shear and shear alignment of uniaxial objects at finite temperatures Peilong Chen Department of Physics, Center for Complex Systems, National Central University, Chungli, Taiwan (Received 23 May 2013; published 4 December 2013) We prove that, for an object with a finitefold rotational symmetry (except for a twofold one) around an axis and mirror symmetries (such as a square rod or pentagonal slab, etc.), dynamics of the symmetry axis in low Reynolds number shear flow exactly follows the same form as that of a uniaxial object (e.g., a circular rod or symmetric ellipsoid) as the so-called Jeffery orbits. We use the formulation in which the dynamics of the rigid body follows first-order ordinary differential equations in time [Phys. Rev. E 84, 056309 (2011)]. Interaction between the object and the shear flow enters through a set of scalar coefficients, and the flow field does not need to be solved dynamically. Results of numerical simulations for general-shaped objects also are discussed. In the second part, Brownian dynamics of a uniaxial object is studied numerically. With D as the rotational diffusion constant, α as a parameter characterizing the aspect ratio, and γ as the shear rate, the object starts to align with the flow when the value of D/(γ α) decreases near 1. At large α (the long object limit), the results suggest much lower flow alignment when D/(γ α) > 1. DOI: 10.1103/PhysRevE.88.063006
PACS number(s): 47.15.G−, 47.57.Qk, 82.70.Dd
I. INTRODUCTION
Dynamics of rigid or deformable bodies in flow fields have been important subjects in fluid dynamics, soft matter, and biophysics. In solutions, bulk properties often depend crucially on how the solutes (microscopic particles or macromolecules) respond to the flow, owing to their distinct particle shapes. Vast varieties of particles have recently been developed in many materials [1]. Microfluidic systems [2] also, very often, incorporate solid objects or macromolecules in the working fluids. In living organisms, many of these flow-response behaviors are essential for the functions of life processes. Some examples are how the various blood cells move with the blood flows in different blood vessels [3] and how bacteria, with their huge variety of body shapes, move under shear flow [4,5]. Many objects in soft matter and biophysics are in small scales such that the relevant Reynolds numbers are much less than 1. The fluid is in the creeping flow regime where the Stokes equation applies [6]. Although these objects can encounter vastly different and complicated flows, as they move with the local mean flow, what they first experience is the local velocity gradient and most notably the shear flows. Thus, one of the fundamental properties is how objects respond to shear flows, which naturally arises in complicated flow patterns or in confined geometries. There have been many experimental and theoretical studies on how macromolecules or colloids respond to a shear flow [7] or even a viscoelastic fluid [8]. However, the first and possibly the only complete analytic result for dynamics of a rigid body inside a sheared Stokes flow is the Jeffery orbits [9] for uniaxial objects (e.g., circular rods or symmetric ellipsoids). For other, “less-symmetric” shapes, numerical simulations are mostly the only way to obtain particle dynamics. Usually, the Oseen tensor and surface integral formulation [10] were used. Such simulations involve, at each time, a first step using the surface integral (on the object’s surface) to numerically calculate the flow field by simultaneously solving for the stress distribution on the object surface as well as linear and rotational velocities of the object, typically with an iterative procedure. 1539-3755/2013/88(6)/063006(7)
In the second step, the computed object velocities are used to advance the position and orientation of the object a small time step forward. By repeating the steps, the time evolution is obtained. Because the second step numerically is easy, especially for nondeformable objects, most of the computing time is spent on the first step of solving the flow field. In this paper, using a formulation previously detailed [11] (although it mostly dealt with analytic results for chiral objects), we will demonstrate that the time-consuming first step of solving the flow field can be avoided. (So simulations can be performed very efficiently.) We then analytically prove a surprising result that, for a rod with a cross section being a regular polygon, the dynamics of the central axis is independent of the orientation of the other two axes. The symmetry axis moves in exactly the same way a uniaxial object as the Jeffery orbits. For general-shaped objects, we will discuss their dynamical patterns using results obtained from simulations. By adding random effects due to finite temperatures, we can quickly obtain statistical properties of the Brownian dynamics. Specifically, we will give detailed results of shear alignment for a uniaxial object. Transition to flow alignment is found to be determined by a particular combination of the shear rate, the rotational diffusion constant, and an object geometry parameter.
II. STOKES FLOW FORMULATION FOR AN OBJECT UNDER SHEAR
In this section, we summarize the formulation, which has been derived previously [11], to study dynamics of a rigid body under shear, and as a special case, the resulting so-called Jeffery-orbit solutions for uniaxial objects. Consider an object moving through a stationary fluid with the Reynolds number Re 1 such that the momentum is negligible and the Stokes equation [6] applies. Denote the instantaneous linear and angular velocities of the object as v and ω, respectively. The induced flow field u1 (r) is zero at
063006-1
©2013 American Physical Society
PEILONG CHEN
PHYSICAL REVIEW E 88, 063006 (2013)
infinity and matches no-slip conditions on the (moving) surface of the object. The total force f 1 and torque τ 1 experienced by the object from the fluid can be written as, due to the linearity of the Stokes equation,
on inversion, the resistance tensor b is zero. Orientation of the object is specified completely by a unit vector nˆ along the symmetry axis. By symmetry, the tensors a and c are in the form of (with I as the second-rank unit tensor)
f 1 = a · v + b · ω, τ 1 = b · v + c · ω.
ˆ c = −c1 nˆ nˆ − c2 (I − nˆ n).
Here (a,b,c) are the second-rank resistance tensors depending on the shape and orientation of the object. The two cross term b tensors are the same due to the Lorentz reciprocal theorem [6]. Next we consider the same object fixed, both in space and in orientation, under an external simple shear. The shear induces a flow u2 (r), which is zero on the object surface (which is fixed) and matches the simple shear flow at infinity. Again the total force f 2 and torque τ 2 experienced by the object from the fluid can be written as f 2 = g:e, τ 2 = h:e. Here g and h are the third-rank shear-response tensors also depending on the shape and orientation. e ≡ ∇us is the second-rank velocity gradient tensor of the simple shear flow us (for example, γ y xˆ and γ is the shear rate). When a free object is subjected to an external shear flow and moves with v and ω due to the shear, u1 + u2 will be the resulting flow field. The combination is the solution of the Stokes equation due to the linearity and has the correct no-slip boundary condition on the object’s surface and simple shear at infinity. Of course, v and ω of the object are still unknown. Now because the momentum is negligible in the zero Reynolds number limit, total force and torque on the object are zero, a · v + b · ω + g:e + f r = 0,
(1)
b · v + c · ω + h:e + τ r = 0,
(2)
which determine v and ω. Here we also add the random force and torque terms to model the effects of finite temperatures. Since the resistance and shear-response tensors are orientation dependent, so are the resulting v and ω. The tensors (a,b,c,g,h) depend on the shape and orientation of the object. For a rigid body, if we define a body coordinate on the object (nˆ 1 ,nˆ 2 ,nˆ 3 ), then the second-rank tensors (a,b,c) can be written in the form of a=
3
aαβ nˆ α nˆ β .
αβ=1
Meanwhile, the third-rank tensors g and h are g=
3
gαβγ nˆ α nˆ β nˆ γ .
(3)
αβγ =1
Given a specific object shape, all these scalar coefficients can be determined by solving the Stokes equation. There are well established numerical methods, which usually are based on the Oseen tensors and surface integrals [10]. Once these coefficients are obtained, the time evolution of the position and orientation of the objects can be obtained by calculating v and ω of the object using Eqs. (1) and (2) without the need to solve the flow field at each time step. For a uniaxial object, such as a circular rod, which is rotationally symmetric with respect to an axis and symmetric
(4)
For shear response tensors, g = 0 and hij k
3 = (h1 ilk nj nl + h2 lj k ni nl + h3 ij l nk nl ).
(5)
l=1
Physically, h1 , h2 , and h3 are the torque experienced by the object under shear when the symmetry axis nˆ is fixed and pointing in the velocity, velocity-gradient, and vorticity directions, respectively. (In this paper, in components, the inner products between tensors are defined, such as (a · h)ij k = l ali hj kl and (h:e)k = l hij k ej i as in Ref. [11].) Now the rotation ω can be obtained from Eq. (2), and the dynamics of nˆ is given by d nˆ = ω × nˆ = −(c−1 · h:e) × nˆ dt ˆ × n. ˆ = γ (h1 /c2 )nx (nˆ × ˆy) × nˆ + γ (h2 /c2 )ny ( xˆ × n) (6) For the last equation, we have assumed, without the loss of generality, the external simple shear flow us = −γ y xˆ . Straightforward algebra would lead to the second line, Eq. (6). So the dynamics is determined by the ratio between h1 and h2 , characterizing the geometry of the object. It is well known that Eq. (6) has closed curve solutions when nˆ is plotted on the surface of a unit sphere as the socalled Jeffery orbits [9] (e.g., the black curves in Fig. 3) with α 2 ≡ h2 / h1 , tan2 θ (α 2 sin2 φ + cos2 φ) = const.
(7)
Here θ and φ are the polar and azimuthal angles of the spherical coordinate with the z axis along the vorticity direction of the shear flow and x in the flow direction. III. DYNAMICS OF GENERAL-SHAPED OBJECTS
In this section, we first present a surprising result for objects with a finitefold rotational symmetry (except for the twofold one) around an axis and mirror symmetries. Typical examples are rods with cross sections of regular polygons. It will be shown that, for such objects besides b always being zero, c and h are again in the forms of Eqs. (4) and (5), respectively. The rotational dynamics of the symmetry axis is, hence, also like Eq. (6) and follows the Jeffery orbits as the uniaxial objects. Two examples of such objects are illustrated in Figs. 1(a) and 1(b) with (nˆ 1 ,nˆ 2 ,nˆ 3 ) as defined. The symmetry axis nˆ 3 will also be denoted as n. Three symmetry operations, which leave these objects unchanged, will be used for the proof. The first is a π rotation around the bisection line nˆ 1 + nˆ 2 such that nˆ 1 ↔ nˆ 2 ,
nˆ 3 → −nˆ 3 .
(8)
The second is the mirror reflection of the plane spanned by the bisection line nˆ 1 + nˆ 2 and nˆ 3 ,
063006-2
nˆ 1 ↔ nˆ 2 ,
nˆ 3 → nˆ 3 .
(9)
DYNAMICS OF FINITE-SYMMETRY AND GENERAL- . . .
PHYSICAL REVIEW E 88, 063006 (2013)
n1
So we only have h123 nˆ 1 nˆ 2 nˆ 3 + permutations, and Eq. (8) requires that
(a)
(b )
n2 n1
n
h = h123 (nˆ 1 nˆ 2 − nˆ 2 nˆ 1 )nˆ 3 + h312 nˆ 3 (nˆ 1 nˆ 2 − nˆ 2 nˆ 1 ) + h132 (nˆ 1 nˆ 3 nˆ 2 − nˆ 2 nˆ 3 nˆ 1 ). (12)
n
(c)
(d ) 1
1
0.5
0.5
z
n2
z
0
0 -0.5
-0.5
1
1 -1
0.5 -1
-0.5
x
0 0
0.5
y
-1
0.5 -1
-0.5 1-1
-0.5
x
0 0
0.5
y
-0.5 1-1
FIG. 1. (Color online) (a) A rod with a triangular cross section. (b) A slab with a pentagonal cross section. (c) One example of the dynamics of nˆ for (b) as a Jeffery orbit. For clarity, other Jeffery orbits obtained with different initial conditions are not shown. (d) The corresponding dynamics of nˆ 1 to (c).
The third is a π rotation about nˆ 1 , nˆ 1 → nˆ 1 ,
nˆ 2 → −nˆ 2 + a nˆ 1 ,
nˆ 3 → −nˆ 3 ,
(10)
with a = 2nˆ 1 · nˆ 2 = 2 cos θ12 and θ12 as the angle between nˆ 1 and nˆ 2 . All the tensors (a,b,c,g,h) should be invariant under Eqs. (8) and (10). However, for reflection Eq. (9), we need to recognize that, although (a,c,g) are tensors, (b,h) are pseudotensors. These can be seen from Eqs. (1) and (2) that (v, f r ,e) are vectors and tensors (also called polar vectors for v and f r ) and (ω,τ r ) are pseudovectors (also called axial vectors). A pseudovector or tensor will acquire an extra negative sign under reflection [12]. When pseudotensors are expanded in terms of (nˆ 1 ,nˆ 2 ,nˆ 3 ), combining the requirements of Eqs. (8) and (9) plus a negative sign, we get the following requirement: b and h have no terms without nˆ 3 . (11) Thus, for b = − 3αβ=1 bαβ nˆ α nˆ β being a pseudotensor, the term nˆ 3 nˆ 3 is excluded by the requirement of Eq. (9) plus a negative sign, (nˆ 1 + nˆ 2 )nˆ 3 by Eq. (8), and 2αβ=1 nˆ α nˆ β by Eq. (11). We are left with (nˆ 1 − nˆ 2 )nˆ 3 , which becomes [(a − 1)nˆ 1 − nˆ 2 ]nˆ 3 under Eq. (10). However, a cannot be 2 since nˆ 1 and nˆ 2 are not parallel to each other. So b = 0. The rotation ω under deterministic dynamics (no random motion) is, thus, determined by c and h in Eq. (2). For h = 3αβγ =1 hαβγ nˆ α nˆ β nˆ γ being a pseudotensor, nˆ 3 nˆ 3 nˆ 3 is excluded by Eq. (10), (nˆ 1 + nˆ 2 )nˆ 3 nˆ 3 is excluded by Eq. (9) plus a negative sign, (nˆ 1 − nˆ 2 )nˆ 3 nˆ 3 is excluded 2 ˆ α nˆ β nˆ γ is excluded by Eq. (11). by Eq. (8), and αβγ =1 n The remaining terms have only one nˆ 3 each. Of these, (nˆ 1 nˆ 1 + nˆ 2 nˆ 2 )nˆ 3 is excluded by Eq. (8). For (nˆ 1 nˆ 1 − nˆ 2 nˆ 2 )nˆ 3 , after Eq. (10), we get a positive sign for nˆ 2 nˆ 2 nˆ 3 contradicting the original, and this term cannot appear from the remaining terms h123 nˆ 1 nˆ 2 nˆ 3 + permutations under Eq. (10). Thus, (nˆ 1 nˆ 1 − nˆ 2 nˆ 2 )nˆ 3 also is excluded.
Now because nˆ 1 × nˆ 2 = sin θ12 nˆ 3 , the h tensor is in exactly the same form as that of a uniaxial object Eq. (5), which only consists of nˆ 3 . (Actually, using a similar analysis, we can get g = 0. Nevertheless, this is not needed for our results as b = 0.) For c = − 3αβ=1 cαβ nˆ α nˆ β being a (regular) tensor, the term (nˆ 1 + nˆ 2 )nˆ 3 is excluded by Eq. (8) and (nˆ 1 − nˆ 2 )nˆ 3 by Eq. (9). So, c = −c11 (nˆ 1 nˆ 1 + nˆ 2 nˆ 2 ) − c12 (nˆ 1 nˆ 2 + nˆ 2 nˆ 1 ) − c33 nˆ 3 nˆ 3 . (13) However, when the third symmetry operation of rotating around nˆ 1 Eq. (10) is applied, we get c = −[(a 2 + 1)c11 + 2ac12 ]nˆ 1 nˆ 1 − c11 nˆ 2 nˆ 2 − (−ac11 − c12 )(nˆ 1 nˆ 2 + nˆ 2 nˆ 1 ) − c33 nˆ 3 nˆ 3 .
(14)
For Eq. (14) to be equal to Eq. (13), the coefficient c12 needs to be proportional to c11 as c12 = −(a/2)c11 . (15) √ Now √ by defining nˆ 1 = (nˆ 1 + nˆ 2 )/ 2 + a and nˆ 2 = (nˆ 1 − nˆ 2 )/ 2 − a [with (nˆ 1 ,nˆ 2 ,nˆ 3 ) as an orthogonal coordinate], we can diagonalize (nˆ 1 ,nˆ 2 ) terms in c. Furthermore, due to Eq. (15), both nˆ 1 nˆ 1 and nˆ 2 nˆ 2 have the same coefficient ca = c11 (2 + a)(2 − a)/4, and c is again in the form of Eq. (4), c = −ca (nˆ 1 nˆ 1 + nˆ 2 nˆ 2 ) − c33 nˆ 3 nˆ 3 = −ca (I − nˆ 3 nˆ 3 ) − c33 nˆ 3 nˆ 3 . With h and c in the same forms as in the case of uniaxial objects, dynamics of the symmetry axis nˆ = nˆ 3 would again follow Eq. (6), d nˆ = ω × nˆ = −(c−1 · h:e) × nˆ dt ˆ × n. ˆ = γ (h132 /ca )nx (nˆ × ˆy) × nˆ + γ (h312 /ca )ny ( xˆ × n) (16) So we finally conclude that the dynamics of nˆ will follow Jeffery orbits, regardless of the orientation of nˆ 1 –nˆ 2 relative ˆ to n. We also perform numerical simulations for such objects from randomly chosen initial conditions to check the dynamics of nˆ 1 and nˆ 2 by solving Eqs. (1) and (2) without the random force and torque terms. We construct the object by placing a small sphere on each of the vortices of Figs. 1(a) and 1(b) [six in total for (a) and ten for (b)]. The hydrodynamic interaction from the Oseen tensor between two well-separated spheres then is used to calculate the tensors (a,b,c,g,h). This approximation, although affecting the exact numerical values of the components of these tensors, does not change the symmetry properties of the tensors. The tensors are checked to conform to the forms obtained in the symmetry analysis above.
063006-3
PEILONG CHEN
PHYSICAL REVIEW E 88, 063006 (2013)
(a)
(b )
1
z 0
x
0
1 -1
x (f)
1
0
0
1 -1
x (g)
1
x
0
1 -1
x
0
x (h)
1
0
1 0 y 1 -1
1
z 0
1 -1 -1 0 y 1 -1
0
1 -1 -1 0 y
z 0
-1 -1
0
1 -1 -1 0 y
z
1
z 0
1 -1 -1 0 y
z
(d )
1
z 0
-1 -1
(e)
(c)
1
z
0
1 -1 -1 0 y 1 -1
x
0
1 -1 -1 0 y 1 -1
x
0
1 0 y 1 -1
FIG. 2. (Color online) (a)–(f) Dynamics of the central axis of a long rod with a rectangular cross section. (g)–(h) For a long rod with a triangular cross section. See the text for detailed geometries.
In Fig. 1(c), we plot one example of the dynamics of nˆ for the pentagonal object, which nicely follows one of the Jeffery orbits, as predicted from previous analysis. Different initial conditions, indeed, lead to different Jeffery orbits. In the meantime, the trajectories of nˆ 1 and nˆ 2 follow complex patterns with the corresponding dynamics of nˆ 1 to Fig. 1(c) shown in Fig. 1(d). Next we discuss the dynamics of general-shaped objects. Because of the infinite possibilities of shapes, we first consider two limiting categories: One dimension of the object is as follows: (1) much larger than the other two, i.e., a long rod and (2) much smaller than the others, i.e., a plate. For objects of category (1), Figs. 2(a)–2(f) show the dynamics of a rectangular cuboid with its three edge lengths denoted as L, W , and H . We use the same method to construct the objects by placing spheres on the vertices of the geometry shape and calculate the tensors (approximately) as above. With L = 200, W = 12, and H = 10, the dynamics nˆ L (t) are shown with six different initial conditions on (nˆ L ,nˆ W ,nˆ H ). We find that the trajectories mostly are similar to the Jeffery orbits when nˆ L is near the vorticity direction as illustrated in Figs. 2(a) and 2(b). When nˆ L is closer to the velocity (v)– velocity gradient (∇v) plane, its trajectory is no longer a closed Jeffery-like orbit and gradually covers a larger extended region as in Figs. 2(c) and 2(d). Near the v-∇v plane, the covered region usually is the largest as shown in Fig. 2(e). However, among the complex extended trajectories, we often see closed trajectories as in Fig. 2(f) with various degrees of complexity. We find that these general behaviors are rather representative of objects resembling a long rod. For rectangular cuboids, when the ratio W/H is closer to 1, we see a larger range of Jeffery-orbit-like trajectories. At W/H = 1, we are back to the case in Fig. 1 with all Jeffery orbits for nˆ L . Meanwhile, if the cross section of the long rod is less symmetric, for example, a triangle with edge lengths of 13, 9, and 7, not only a smaller region with simple Jeffery-orbit-like trajectories as in Fig. 2(g) with a rod length of 100 is there, but also most of the complex trajectories cover a large portion of the whole sphere. One example is shown in Fig. 2(h). For objects of category (2) as thin plates, the dynamics of nˆ H along the shortest axis are very similar to those of a long
rod as discussed. Again when L/W (for rectangular plates) approaches unity, there is a larger range of Jeffery-like orbits. Dynamics for less symmetric cross sections also follow the same pattern for the long rod. Besides the two limiting shapes, for objects with all three dimensions being comparable, we find that there are much less Jeffery-orbit-like trajectories. Furthermore, most of the trajectories could cover a large portion of the sphere surface similar to Fig. 2(h). IV. BROWNIAN DYNAMICS OF A UNIAXIAL OBJECT
In physical situations, especially in microscopic length scales, the orientation of an object will be perturbed by random forces due to finite temperatures and, hence, will not follow the Jeffery orbits exactly, even for a uniaxial object. The random effect typically is quantified by the P´eclet number Pe ≡ γ /D as the ratio between the shear rate γ and the rotational diffusion constant D, which is proportional to the temperature. For a uniform circular rod as an example, in the limit of Pe 1 (high temperatures), shear effects will be small, and the rod will orient more or less uniformly in all directions. At lower temperatures, although the rod would follow Jeffery orbits closely, the random force provides a transition mechanism between different Jeffery orbits. Thus, with random forces, there could be a steady-state distribution on the unit sphere at long times. The time needed to reach the steady-state distribution will increase with decreasing temperature. To find these steady-state distributions at finite temperatures, we perform numerical simulations of Eq. (6) with an ˆ This random rotation additional random rotation term ωr × n. always is chosen on the plane perpendicular to the direction of the symmetry axis at each time since the component in the direction of the symmetry axis has no effect on the orientation. In this way, the magnitude of ωr is assumed to follow a Gaussian distribution with the variance |ωr |2 equal to the diffusion constant D. The temperature T and D follow the Einstein relation D = BkB T /ζ with kB as the Boltzmann constant and ζ as the fluid viscosity. B is a factor depending on the exact object geometry, in general, with a highly nontrivial dependence. For our discussions, D is used as the control
063006-4
DYNAMICS OF FINITE-SYMMETRY AND GENERAL- . . .
PHYSICAL REVIEW E 88, 063006 (2013)
1
1
α2 = 3 10 30 100 300 1000
0
λ
Vorticity
0.5
0.5
-0.5
-1 0 -1
0 -0.5 0 Flow
0.5
-0.5
1 1
0 0.5 Velocity gradient
-1
Average density (arb. units)
parameter. We will relate the results to T later for a long circular rod. For α 2 = 10, D = 0.1, and γ = 1, the dynamics of nˆ is shown in Fig. 3 by plotting blue points (the upper part) on the ˆ i ), i = 1 · · · N unit sphere representing the orientation of n(t with ti = ts + i t, t as a chosen suitable time interval, and ts as a transient time interval such that the results do not depend on the initial condition. The Jeffery orbits of Eq. (7) at this value of α also are drawn as the black curves, although for clarity, only those on the upper semisphere are drawn. The concentration of blue dots on the two polar regions (of the flow directions) shows the preference of the rod to align near the flow direction under shear. The alignment can further be seen in the projection of the blue points onto the v-∇v plane (the red points on the lower part). To perform a quantitative comparison between different values of D, in Fig. 4, we plot the average densities as functions of the polar angle θf with respect to the flow direction. The two concentrated regions in Fig. 3 correspond to θf = 0 and π at which the density has peak magnitudes as in Fig. 4. Higher peak values and narrower widths indicate stronger alignment in the flow direction. In Fig. 4, we see the degree of alignment increasing with lower diffusion constant D. Furthermore, the distributions have a low-D limit for D below 0.2. This low-D limit of f (r) can actually be calculated using a lowtemperature expansion of the corresponding Fokker-Planck
10
2
FIG. 5. (Color online) The alignment index λ as a function of D/(γ α).
FIG. 3. (Color online) Brownian dynamics of nˆ showing shear alignment for α 2 = 10, D = 0.1, and γ = 1. See the text for details.
15
1 D/(γ α)
D = 0.02 0.05 0.1 0.2 0.5 1 2 4
5
0
equation for f (r) [13] (or see the Appendix). The contour lines of the calculated f (r) are plotted in Fig. 3 as the red (gray) contour curves (the mostly vertical ones). They are seen to match the distribution of the blue points from simulations well. A more quantitative comparison is seen by drawing the average density as a function of the polar angle in the flow direction in Fig. 4 as the red solid line, which perfectly agrees with the simulation results. For reference, our formulations of the low-T expansion of the Fokker-Planck equation are detailed in the Appendix, which are the same, in principle, with those in Ref. [13]. The peaks in Fig. 4 signify the degree of flow alignment of the rod under the competition between the random motion and the shear flow. To quantify the alignment, we calculate an alignment index λ, which is the portion of distribution above the “background” value at θf = π/2, λ = [f (θf ,φf ) − f (θf = π/2,φf ) φf ]d f , with φf denoting the average over φf and f denoting the solid angle for (θf ,φf ). A uniform distribution without flow alignment has λ = 0, and λ = 1 is the maximum. The results are shown in Fig. 5. It shows the object only starts to have significant alignment with the flow when D/(γ α) is near 1. Furthermore, the data suggest that, at the limit of long rods (α → ∞), the flow alignment remains quite insignificant when D/(γ α) > 1. More accurate numerical calculations could be needed to say whether λ → 0 at α → ∞ and D/(γ α) > 1. On the other hand, the maximum flow alignment λm at D → 0 significantly depends on α. We find that they closely follow a scaling relation 1 − λm ≈ α −0.43 . As an example for long circular rods, it is known that the rotational diffusion constant D goes mostly as ∼ T /L3 plus another correction term involving ln(L/a) [14] with L and a being the length and radius of the rod, respectively. On the other hand, for a rod in Stokes flow, we expect h1 ∼ L2 and h2 ∼ L, thus, α 2 ∼ L. So Fig. 5 tells us that the transition to significant flow alignment, at least for long rods, is determined approximately by T L−3.5 .
0 0.5 1 Polar angle in flow direction (units of π)
V. CONCLUSIONS
FIG. 4. (Color online) Average density as a function of the polar angle in the flow direction for α 2 = 10 and γ = 1.
To summarize, we prove that the central axis of a rod with a cross section of a regular polygon moves exactly in the same
063006-5
PEILONG CHEN
PHYSICAL REVIEW E 88, 063006 (2013)
patterns with those of circular rods or symmetric ellipsoids under external shear flow. Using our formulation, temporal dynamics can be simulated very efficiently for general-shaped objects, which possess different degrees of Jeffery-orbit-like trajectories near the vorticity direction. For uniaxial objects, we obtain shear alignment distributions at finite temperatures depending on the object shapes. Transition to flow alignment is found to mostly be controlled by the combination D/(γ α), and the maximum flow alignment at low temperature scales with α. ACKNOWLEDGMENT
The support of the National Science Council of Taiwan is acknowledged.
with k(u) to be determined. On the order of , we have ∂ hu ∂f0 ∂ hv ∂f0 ∂ hu Kf1 = + . ∂v ∂u hu ∂u ∂v hv ∂v Now if we integrate the equation by requirement of periodicity on v, we have ∂ ∂u
In this Appendix, we outline the low-temperature expansion of the Fokker-Planck equation [13]. It is well established that the Brownian dynamics of nˆ as a trajectory on the unit sphere has a corresponding Fokker-Planck equation [15] for the distribution function f (r) with r being the coordinate on the surface of the sphere, ∂f = −∇ · (K f ) + ∇ 2 f. ∂t The drift velocity K is the deterministic velocity induced by the shear and corresponds to the right hand side of Eq. (6). The second term is the diffusion of f (r), and the diffusion constant is proportional to the temperature. We want to solve for steady-state solutions: ∂f/∂t = 0 at 1. The distribution is expanded in the orders of , f = f0 + f1 + f2 + · · · . 2
Collecting terms with the same order of , we obtain ∇ · (K f0 ) = 0,
∇ · (K fi ) = ∇ 2 fi−1 , i 1.
It is noted that, with the small term having the highest order of derivative, an expansion perturbation often becomes singular [16], leading to, for example, boundary layers. However, in the current case, the space on the unit sphere is finite and has no boundaries. There also are no sink or source terms. So the perturbation could remain regular. Assume (u,v) is a curvilinear coordinate on the surface of the sphere. Note that the drift velocity K is tangent to the Jeffery orbits everywhere. Thus, if we choose u as the Jeffery orbits, u = tan2 θ (α 2 sin2 φ + cos2 φ), then K = K(u,v)ˆv . The coordinate v traces around the closed Jeffery orbits and can always be defined as 0 v 2π with v = 0 and v = 2π being the same. All functions are periodic in v. We have, on the order of 0 with hu and hv , the scale factors of the (u,v) coordinate,
dv 0
dv, from the
hv ∂f0 = 0. hu ∂u
(A2)
dk(u) + B(u)k(u) = C1 , du
with
2π
A(u) ≡ 0
hv dv G, B(u) ≡ hu
2π
dv 0
hv ∂G . hu ∂u
(A3)
C1 equals the integration in Eq. (A2) and, physically, is exactly the total flux across a contour line of a fixed u (Jeffery orbits), the same for all values of u for a steady-state distribution. Because the space is finite and there are no sink or source contributions, the constant flux is necessarily zero. So we have the lowest-order solution f0 (u,v) = [hu (u,v)K(u,v)]−1 k(u) with B(u) du . (A4) k(u) = C2 exp − A(u) C2 is determined by the normalization condition for f0 . Because we have an analytic form for the coordinate v as detailed below, we have an analytic expression for the integrands in Eq. (A3). We, thus, can obtain k(u) by numerically integrating Eq. (A3) to get A(u) and B(u) and then using them in Eq. (A4). So we choose the u coordinate of (u,v) as the Jeffery orbits, u(θ,φ) = tan2 θ (α 2 sin2 φ + cos2 φ) ≡ a(θ )b(φ). Here θ and φ are the angular parts of the spherical coordinate. To find v(θ,φ), we require 0 = ∇u · ∇v =
1 ∂u ∂v ∂u ∂v + . ∂θ ∂θ sin2 θ ∂φ ∂φ
Assuming v(θ,φ) = f (θ )g(φ) with some straightforward algebra, we obtain with some particular choices of integration constants (the choices will not affect the final calculated distribution function), u(θ,φ) = tan2 θ (α 2 sin2 φ + cos2 φ),
1 ∂ hu Kf0 = 0. hu hv ∂v
v(θ,φ) = sinα
Thus, hu Kf0 is a function of u only, f0 (u,v) = (hu K)−1 k(u) ≡ G(u,v)k(u),
2π
0
The integration in the above equation will equal a constant C1 . Thus, in the order, we obtain an equation for f0 (r) to determine k(u). Using Eq. (A1), we get A(u)
APPENDIX: LOW-TEMPERATURE SOLUTIONS OF THE FOKKER-PLANCK EQUATION FOR A UNIAXIAL OBJECT
2π
(A1)
2
−1
2
θ cosα φ/ sin φ.
The scale factors hu and hv can also be found as functions of θ and φ.
063006-6
DYNAMICS OF FINITE-SYMMETRY AND GENERAL- . . .
PHYSICAL REVIEW E 88, 063006 (2013)
[1] See, e.g., B. S. John, A. Stroock, and F. A. Escobedo, J. Chem. Phys. 120, 9383 (2004); C. J. Murphy et al., J. Phys. Chem. B 109, 13857 (2005); L. Rossi et al., Soft Matter 7, 4139 (2011); Y. Wang et al., Nature (London) 491, 51 (2012). [2] Microfluids, edited by S. Colin (Wiley, London, 2010); B. J. Kirby, Micro- and Nanoscale Fluid Mechanics: Transport in Microfluidic Devices (Cambridge University Press, Cambridge, UK, 2010). [3] R. D. Kamm, Annu. Rev. Fluid Mech. 34, 211 (2002). [4] J. T. Locsei and T. J. Pedley, Bull. Math. Biol. 71, 1089 (2009). [5] T. Ahmed, T. S. Shimizu, and R. Stocker, Integr. Biol. 2, 604 (2010). [6] J. Happel and H. Brenner, Low Reynolds Number Hydrodynamics with Special Applications to Particulate Media (PrenticeHall, Upper Saddle River, 1965).
[7] See, e.g., Marcos, H. C. Fu, T. R. Powers, and R. Stocker, Phys. Rev. Lett. 102, 158103 (2009); W. Tang and S. G. Advani, J. Chem. Phys. 126, 144711 (2007), and references therein. [8] Z. Fan and S. G. Advani, Polymer 46, 5232 (2005). [9] G. B. Jeffery, Proc. R. Soc. London, Ser. A 102, 161 (1922). [10] See, e.g., M. Makino and M. Doi, Phys. Fluids 17, 103605 (2005). [11] P. Chen and Q. Zhang, Phys. Rev. E 84, 056309 (2011). [12] G. B. Arfken, Mathematical Methods for Physicists, 3rd ed. (Academic Press, Orlando, FL, 1985), pp. 128–137. [13] L. G. Leal and E. J. Hinch, J. Fluid Mech. 46, 685 (1971). [14] S. J. Broersma, Chem. Phys. 32, 1626 (1960); 32, 1632 (1960). [15] R. Zwanzig, Nonequilibrium Statistical Mechanics (Oxford University Press, New York, 2001). [16] E. J. Hinch, Perturbation Methods (Cambridge University Press, Cambridge, UK, 1991).
063006-7
