Similar articles

Abstract: We present the simulation of the dynamics of fluid-cylinder interactions in a narrow three-dimensional channel filled with a Newtonian fluid, using a Lagrange multiplier based fictitious domain methodology combined with a finite element method and an operator splitting technique. As expected, a settling truncated cylinder turns its broadside perpendicular to the main stream direction and the center of mass moves to the central axis of the channel. In the case of two truncated cylinders, they first move around each other for a while and then stay together in a “T” shape. After the “T” shape has been formed for a long enough time, we found no vortex shedding behind the cylinders. When simulating the fluidization of 60 truncated cylinders, we captured the features of interactions among fluidized cylinders as observed in experiments.

.  INTRODUCTION

 The orientation of symmetric long bodies, e.g. ellipsoid and truncated cylinder, settling in liquids of different nature is a fundamental issue in many problems of practical interest (Liu and Joseph, 1993). In the past decade, a fast growing population of researchers developed numerical methods for direct simulation of fluid/particle interaction (Hu et al., 1992; Ladd, 1994a; 1994b; Qi and Luo, 2003; Hu, 1996; Johnson and Tezduyar, 1997; Maury and Glowinski, 1997; Hofler et al., 1998; Glowinski et al., 1999a; 2001; Takagi et al., 2003; Dong et al., 2004; Juárez et al., 2002). In this article we discuss first the generalization of a Lagrange multiplier based fictitious domain method (Glowinski et al., 1999a; 2001) for simulating the motion of particles with edges and corners, such as truncated cylinders, in a Newtonian fluid. Unlike the cases where the particles are spheres, we attach two points to each particle and move them according to the rigid-body motion of the particle. The equations describing the motion of these two points are solved by a distance preserving scheme so that rigidity can be maintained. Then we apply the above methodology to simulate the dynamics of fluid-cylinder interactions in a narrow three-dimensional channel filled with a Newtonian fluid. As expected, a settling truncated cylinder turns its broadside perpendicular to the main stream direction and the center of mass moves toward the central axis of the channel. In the case of two truncated cylinders, they first move around each other for a while after being released from their initial positions and then stay together in a “T” shape, as observed in experiments. We found that there is no vortex shedding behind the cylinders after the “T” shape has been formed for a long enough time. After increasing the density of the cylinders or reducing their length, we still obtained similar interaction between two settling truncated cylinders. In the simulation of fluidization of 60 truncated cylinders, we successfully captured the features of the interactions between fluidized cylinders as observed in experiments.

.  A MODEL PROBLEM AND A FICTITIOUS DOMAIN FORMULATION FOR 3-D PARTICULATE FLOW

 To perform the direct numerical simulation of the interaction between rigid bodies and fluids, we developed a methodology which combines a distributed Lagrange multiplier based fictitious domain (also called domain embedding) method with operator splitting methods (Glowinski et al., 1998b; 1999a; 1999b; 2001). In (Glowinski et al., 2001; Pan et al., 2001; 2002) we have applied this methodology to simulate the fluidization of a thousand spherical particles in 3-D and the sedimentation of ten thousand circular particles in 2-D. This approach (or closely related ones derived from it) has also been used by other investigators (Dong et al., 2004; Baaijens, 2001; Wagner et al., 2001; Yu et al., 2002; Diaz-Goano et al., 2003). We are going to recall the ideas at the basis of the above methodology by considering the motion of a single particle in a Newtonian viscous incompressible fluid (of density ρ_f and viscosity μ_f) under the effect of gravity. For the situation depicted in Fig.1 (a case where the particle has edges), the flow is modeled by the Navier-Stokes equations, namely, (with obvious notation) (1) (2) (3) (4) where Γ=∂Ω, g denotes gravity and n is the unit normal vector pointing outward to the flow region. We assume a no-slip condition on γ (=∂B). The motion of particle B satisfies the Euler-Newton’s equations, namely (5) (6) (7)




 with the resultant and torque of the hydrodynamical forces given respectively by, with . Eqs.(1)−(7) are completed by the following initial conditions (8)

 and by the no-slip boundary conditions (9)

 M_p, I_p, G, V and ω above are the mass, inertia, center of mass, velocity of the center of mass and angular velocity of particle B, respectively. In Eq.(7) we found it preferable to deal with the kinematic angular momentum I_pω, thus making the formulation more conservative for non-spherical particle cases. To avoid particle-particle and particle-wall penetration which can happen in the numerical simulation, we introduced a artificial short-range repulsion force F ^r in Eq.(6), which becomes active when the shortest distance between two (convex) particles or between (convex) particle and wall is less than a pre-chosen distance (Glowinski et al., 1999a; 2001; Maury, 1997) and then a torque in Eq.(7) acting on the point x_r where F ^r applies on B (see Remark 3 in Section 3.3 for the location of x_r). For non-convex particles, we can apply similar approach to activate the short-range repulsion force F ^r.

 To solve system Eqs.(1)−(9) we can use, for example, Arbitrary Lagrange-Euler (ALE) methods as in ( Hu, 1996; Johnson and Tezduyar, 1997; Maury and Glowinski, 1997), or fictitious domain methods, which allow the flow calculation to be done on a fixed grid, as in (Glowinski et al., 1998b; 1999a; 1999b; 2001). The fictitious domain methods that we advocate have some common features with the immersed boundary method of Peskin (Peskin, 1977; 1981; Peskin and McQueen, 1980), but also some significant differences in the sense that we take advantage of distributed Lagrange multipliers to force the rigid body motion inside the particle. As with the methods in (Peskin, 1977; 1981; Peskin and McQueen, 1980), our approach takes advantage of the fact that the flow can be computed on a grid which does not have to vary in time, a substantial simplification indeed.

 The principle of the fictitious domain method that we employ is rather simple. It consists in

 (1) Filling the particles with a fluid having the same density and viscosity as the surrounding one.

 (2) Compensating the above step by introducing, in some sense, an anti-particle of mass (−1)M_pρ_f /ρ_s and inertia (−1)I_pρ_f /ρ_s, taking into account the fact that any rigid body motion v(x, t) verifies ∇˙v=0 and D(v)=0 (ρ_s: particle density).

 (3) Finally, assuming the fluid contained in B(t) has rigid body velocity, (10)

 via a Lagrange multiplier λ supported by . The vector λ forces rigidity in B(t) in the same way that ∇p forces ∇˙v=0 for incompressible fluids.

 We obtain then an equivalent formulation of Eqs.(1)−(9) defined on the whole domain Ω, namely

 For a.e. t>0, find {u(t), p(t), V(t), G(t), ω(t), λ(t)} such that (11)

 and (12) (13) (14) (15) (16) (17)

 with the following functional spaces Various examples for in Eqs.(12) and (15) are given in (Glowinski et al., 2001; Glowinski, 2003). In Eqs.(11)−(17), only the center of mass, the translation velocity of the center of mass and the angular velocity of the particle are considered, which are sufficient for spherical particle cases. For non-spherical particles, knowing these two velocities and the center of mass of the particle, we can translate and rotate the particle in space by tracking, in each particle, two extra points x₁ and x₂ which follow the rigid body motion, namely (18)

 In practice we track two orthogonal normalized vectors rigidly attached to the body B and originating from the center of mass G.

 Remark 1 The second gravity term on the right-hand-side of Eq.(12) can be combined with the pressure. Hence in the following discussion, this term will not be mentioned anymore.

.  TIME AND SPACE DISCRETIZATION

.  Lie’s first order operator-splitting scheme

 Many operator-splitting schemes can be applied to problem (11)−(18). One advantage of operator-splitting schemes is that we can decouple difficulties such as (i) the incompressibility condition, (ii) the nonlinear advection term, and (iii) the rigid body motion, so that each one of them can be handled separately, and in principle optimally. Let Δt be a time discretization step and t^n+s=(n+s)Δt. Application of Lie’s first order operator-splitting scheme (Chorin et al., 1978), to problem (11)−(18), yields: (19)

 for n≥0, uⁿ[ u(tⁿ)], Gⁿ, Vⁿ, ωⁿ, and being known, we first compute u^n+1/6, p^n+1/6 via the solution of (20)

 and set u^n+1/6=u(tⁿ⁺¹), p^n+1/6=p(tⁿ⁺¹).

 Next, compute u^n+2/6 via the solution of (21)

 and set u^n+2/6=u(tⁿ⁺¹).

 Then, compute u^n+3/6 via the solution of (22)

 and set u^n+3/6=u(tⁿ⁺¹).

 Now predict the motion of the center of mass and the angular velocity of the particle via the solution of (23) (24) (25) (26) (27)

 and set G^n+4/6=G(tⁿ⁺¹), V^n+4/6=V(tⁿ⁺¹), (I_pω)^n+4/6=

 (I_pω)(tⁿ⁺¹), and Using G^n+4/6, and obtained in the above step, we enforce the rigid body motion on the region B^n+4/6 occupied by the particle (28) (29)

 and set uⁿ⁺¹=u(tⁿ⁺¹), V^n+5/6=V(tⁿ⁺¹), (I_pω)^n+5/6=

 (I_pω)(tⁿ⁺¹).

 Correct the motion of the center of mass and the angular velocity of the particle via the solution of (30) (31) (32) (33) (34)

 Then set Gⁿ⁺¹=G(tⁿ⁺¹), Vⁿ⁺¹=V(tⁿ⁺¹), (I_pω)ⁿ⁺¹=

 (I_pω)(tⁿ⁺¹), and In Eqs.(19)−(34), we have with B^n+4/6 being the region occupied by the particle B according to G^n+4/6, and and 0≤α,β≤1, with α+β=1. In numerical simulations, we always choose α=1 and β=0.

.  Space discretization

 We assume that and is a rectangular parallelepiped. For the finite element approximation of problem (11)−(18), we employ a Bercovier-Pironneau one, namely (Glowinski, 2003) (35) (36) (37) where τ_h is a “tetrahedrization” of Ω, τ_2h is twice coarser than τ_h, and P₁ is the space of the polynomials in three variables of degree≤1. A finite dimensional space approximating Λ(t) is obtained as follows: let be a set of points from which cover (uniformly, for example); we define then (38) where δ(˙) is the Dirac measure at x=0. Then we use <˙,˙>_h defined by (39)

 A typical choice of points for defining Eq.(38) is to take the grid points of the velocity mesh internal to the particle B and whose distance to the boundary of B is greater than, e.g. h/2, and to complete with selected points from the boundary of B(t) (see, e.g., Fig.2 for an example of selected collocation points on the boundary of B(t)).




 Using the above finite dimensional spaces and the backward Euler’s method for most of the subproblems in scheme (19)−(34), we obtain the following scheme after dropping some of the subscripts h [similar schemes are discussed in (Glowinski et al., 1998b; 1999a; 1999b; 2001)]: (40)

 For n≥0, uⁿ( u(tⁿ)), Gⁿ, Vⁿ, ωⁿ, and being known, we first compute u^n+1/6, p^n+1/6 via the solution of (41)

 Next, compute u^n+2/6 via the solution of (42)

 and set u^n+2/6=u(tⁿ⁺¹). Then, compute u^n+3/6 via the solution of (43)

 Now predict the motion of the center of mass and the angular velocity of the particle via the solution of (44) (45) (46) (47)

 G(tⁿ) =G ⁿ, V(tⁿ) =V ⁿ, (I_pω)(tⁿ)=(I_pω)ⁿ, x₁(tⁿ) (48)

 Then set G^n+4/6=G(tⁿ⁺¹), V^n+4/6=V(tⁿ⁺¹), (I_pω)^n+4/6 =(I_pω)(tⁿ⁺¹), and With G^n+4/6, and obtained at the above step, we enforce the rigid body motion on the region B(t^n+4/6) occupied by the particle via the solution of (49) (50)

 Correct the motion of the center of mass and the angular velocity of the particle via the solution of (51) (52) (53) (54) (55)

 Then set Gⁿ⁺¹=G(tⁿ⁺¹), Vⁿ⁺¹=V(tⁿ⁺¹), (I_pω)ⁿ⁺¹=

 (I_pω)(tⁿ⁺¹), and In Eqs.(40)−(55), <0} and = is an approximation of belonging to and verify Remark 2 The rate of convergence (with respect to h) of scheme (40)−(55), when using the Taylor-Hood finite element approximation instead of the Bercovier-Pironneau one, was investigated in (Juárez et al., 2002) for a circular particle in a Stokes flow, the error orders in the L²-norm for velocity and pressure slightly exceed 3 and 2, respectively.

.  On the solution of subproblems (41), (42), (43), (44)−(48) and (49)−(50)

 The degenerated quasi-Stokes problem (41) is solved by an Uzawa/preconditioned conjugate gradient algorithm as in (Glowinski et al., 1998b), where the discrete elliptic problems used for preconditioning are solved by a matrix-free fast solver from FISHPAK due to Adams et al.(1980). The stopping criterion for the preconditioned conjugate gradient algorithm is ||r^k||/||r⁰||≤ε where r^k is the residue at the kth iteration. It typically takes about 10 iterations in these simulations to achieve convergence when ε =10⁻⁵. The advection problem (42) for the velocity field is solved by a wave-like equation method as in (Dean and Glowinski, 1997; Dean et al., 1998). Problem (43) is a classical discrete elliptic problem which can be solved by the above matrix-free fast solver.

 Systems (44)−(48) and (51)−(55) are systems of ordinary differential equations thanks to operator splitting. For their solutions one can choose a time step smaller than Δt, (i.e. we can divide Δt into smaller steps) to predict the translation velocity of the center of mass, the angular velocity of the particle, the position of the center of mass and the regions occupied by each particle so that the repulsion forces can effectively prevent particle-particle and particle-wall overlapping. At each sub-cycling time step, keeping the distance constant between points x₁ and x₂ in each particle is important since we are dealing with rigid particles. To satisfy the above constraint we applied the following approach:

 1. Translate x₁ and x₂ according to the new position of the mass center at each subcycling time step; 2. Rotate Gx₁ and Gx₂, the relative positions of x₁ and x₂ to the center of mass G, by the following Crank-Nicolson scheme (a Runge-Kutta scheme of order 2, in fact): (56)

 for i=1, 2 with τ as a sub-cycling time step. By Eq.(56), we have , for i=1, 2 and (i.e., scheme (56) is distance and in fact shape preserving).

 Remark 3 In order to activate the short range repulsion force in Eqs.(6) and (7), we have to find the shortest distance between two truncated cylinders. Unlike the cases for spheres, it is not trivial to locate the points on the surface of each truncated cylinder so that their distance is the shortest possible between two truncated cylinders. There is no explicit formula for such distance. The difficulties stem from the fact that a truncated cylinder is not a strictly convex set; it is “flat” in the direction of the central axis and on the disks at the two ends. When two cylinders are parallel to each other, i.e., when the two central axes are parallel to each other, the question which arises is: at which points do we have to apply the repulsion force when it is needed As shown in Fig.3, the shortest distance is between two line segments on the cylinder surfaces on the plane passing through two central axes. We choose the middle point on each line segment as the point at which we apply the repulsion force as shown in Fig.3. When the plane containing one end disk is parallel to the central axis of another cylinder, e.g., two cylinders are in a “T” shape as in Fig.4, we also first obtain the line segment from the circular surface of one cylinder and that from the disk at one end of another cylinder so that the shortest distance is between two line segments. Then we choose the middle point on each line segment as the point at which we apply the repulsion force. When the cylinders are not parallel to each other, we first choose a set of points from the surface of each truncated cylinder. Then we find the point among the chosen points from each surface at which the distance is the shortest. We repeat this (kind of relaxation) process in the neighborhood of the newly located point on each surface of the cylinder until convergence, usually obtained in very few iterations. For the shortest distance between the wall and truncated cylinder, we can take the advantage of the shape of truncated cylinder and only search points on the circles at the two ends and then determine the shortest distance to each wall and whether it is parallel to the wall. When it is parallel to the wall then we use the midpoint on the line segment, which is either parallel to the central axis and on the surface or on the plane containing the disk at the end, as the point for applying the repulsion force. These technical details have to be implemented in order to capture the interactions between cylinders.







 The rigid body motion is enforced in B(t^n+4/6), via Eq.(50). At the same time those hydrodynamical and gravity forces acting on the particles are also taken into account in order to update the translation and angular velocities of the particles. To solve Eqs.(49)−(50), we use a conjugate gradient algorithm as discussed in (Glowinski et al., 1999a). Since we take β=0 in Eq.(49) for the simulation, we actually do not need to solve any non-trivial linear systems for the velocity field; this saves a lot of computing time. To get the angular velocity ωⁿ⁺¹, computed via (57)

 We need to have , the inertia of the particle B(t^n+4/6). We first compute the inertia I₀ in the coordinate system attached to the particle. Then via the center of mass G^n+4/6 and points and , we have the rotation transformation Q (QQ^T=Q^TQ=I_d, detQ=1) which transforms vectors expressed in the particle frame to vectors in the flow domain coordinate system and . Actually in order to update matrix Q we can also use quaternion techniques, as shown, in the review paper (Chou, 1992).

.  NUMERICAL EXPERIMENTS

.  Settling of a truncated cylinder in a 2D-like channel

 In the first test case, we consider the simulation of the motion of a truncated cylinder settling in an infinite length narrow channel filled with a Newtonian fluid. The computational domain is Ω=(0, 1)×(0, 0.25)×(0, 4) initially, then it moves down with the mass center of the cylinder (see, e.g., (Hu et al., 1992) for adjusting the computational domain according to the position of the particle). The fluid density is ρ_f=1 and the fluid viscosity is μ_f=0.01. The flow field initial condition is u=0. The density of the cylinder is ρ_s= 1.1. The radius of the cylinder is 0.1 and its length is 0.4. Initially its central axis is in the vertical direction. The initial position of the mass center is at (0.5, 0.125, 1). The initial translation and angular velocities of the cylinder are 0. To check the convergence, we chose the following two pairs of the mesh size of the velocity field and the time step: {h_v, Δt}={1/80, 0.001}, {1/112, 0.001}. The mesh size of the pressure is always h_p=2h_v.

 The snapshots of the cylinder position and its orientation for {h_v, Δt}={1/112, 0.001} at different times in the channel are shown in Fig.5. The computations were performed in a moving frame of reference, so the cylinder does not appear moving downward. In Fig.6 we show the histories versus time of the x₁-coordinate of the mass center, the x₁- and x₃-components of the translation velocity and the x₂-component of the angular velocity. We obtained a good agreement between these results. The movement of the cylinder is very fast at the beginning when it moves toward the side wall after being released from its initial position. Later on the oscillations were damped out (Fig.6). As expected, the truncated cylinder turns its broadside perpendicular to the main stream direction and the mass center moves back to the central axis of the channel. The velocity fields projected on the x₁x₃-plane passing through the mass center are shown in Fig.7 for h_v=1/112 and Δt=0.001.










 For h_v=1/112 and Δt=0.001, the averaged particle velocity when the cylinder reaches a stable position is about 2.255 (the maximal velocity is about 4.681) so the averaged particle Reynolds number with the longest axis as characteristic length is 90.2. The number of nodes for the velocity field is 546021 (resp., 1471373) for h_v=1/80 (resp., 1/112). The memory used in the simulation is about 69 MB (resp., 184 MB) for h_v=1/80 (resp., 1/112). The simulation takes about 14 s (resp., 38 s) per time step for {h_v, Δt}={1/80, 0.001} (resp., {1/112, 0.001}) on a Linux based PC with 1.6 GHz AMD Athlon CPU.

.  Two truncated cylinder settling in a 2D-like channel

 Fig.10.4.9 in (Joseph, 1992) shows experimental results; we can observe many types of interactions among truncated cylinders when fluidized in a 2D-like bed (there is only one layer of particles in a narrow channel). Even though they are fluidized, the features of interaction between two truncated cylinders are similar to those of two truncated cylinders settling in a 2D-like channel. To reproduce some of those features computationally, we consider the following test case. The initial mass center positions are (0.45, 0.125, 1) and (0.45, 0.125, 1.3), respectively. The central axes are in the horizontal direction initially so the frames rigidly attached to the cylinders are {(1, 0, 0), (0, 1, 0), (0, 0, 1)} initially. The radii are 0.1 and their lengths are 0.46; we have taken {1/80, 0.001} and {1/112, 0.001} for {h_v, Δt}. All other parameters are as in the previous case.




 Snapshots of the positions and orientations of two cylinders for {h_v, Δt}={1/112, 0.001} are shown in Fig.8. The cylinders first move around each other and then stay together in a “T” shape, as had been observed in experiments. Snapshots of the velocity field for {h_v, Δt}={1/112, 0.001} projected on the x₁x₃-plane passing through the mass center of the upper cylinder are shown in Fig.9. As observed in the experiments shown in (Joseph, 1992), we observed that in Fig.9 the wake behind the lower cylinder makes it and the other cylinder stay together to form a “T” shape configuration. We also observed that there is no vortex shedding behind the cylinders after the “T” shape configuration has been formed for a sufficiently long time (see snapshots at t=5 and 10). For {h_v, Δt}={1/80, 0.001} and {1/112, 0.001}, the histories of the x₁-coordinates, the x₁- and x₃-components of the translation velocities of the mass center, and the x₂-component of the angular velocities are shown in Fig.10 for 0≤t≤2. We obtained a good agreement between results, considering the violence of the motion during the first two time units. For the last two time units (8≤t≤10), the averaged particle speed is 2.794 so that the averaged particle Reynolds number with the cylinder length as characteristic length is 128.524. The maximal averaged speed of the two cylinders is 2.971 for 0≤t≤10.







 To see the effect of the length and the density of the cylinder, we considered the two following tests. We only either increase the density of cylinders to 1.25 or reduce the lengths to 0.4 and keep all other parameters the same as above with {h_v, Δt}={1/112, 0.001}. The numerical results of both cases are similar to those of the above case. Histories of the translation velocity of the mass center are shown in Fig.11.




 The memory used in the simulation was about 71 MB (resp., 190 MB) for h_v=1/80 (resp., 1/112) for the first test case. The simulation takes about 5.76 s (resp., 14.34 s) per time step for {h_v, Δt}={1/80, 0.001} (resp., {1/112, 0.001}) on a Linux based PC with 2.2 GHz Opteron 248 CPU.

.  Fluidization of 60 truncated cylinders in a 2D-like channel

 To simulate the fluidization of truncated cylinders similar to the one shown in Fig.10.4.9 in (Joseph, 1992), we considered the following test case. The computational domain is Ω=(0, 3)×(0, 0.25)×(0, 8). The initial mass center positions and orientations are shown in Fig.13 with randomly chosen angle. The radii are 0.1, their lengths are 0.4 and {h_v, Δt}={1/80, 0.001}. All other parameters are as in the case in Section 4.1. At the bottom and the top of the channel, the velocity field is (0, 0, U(t)) with U(t)=1.5(1−e^−20t). In Figs.12, 13 and 14 the snapshots of the position and orientation of the cylinders and the velocity fields projected to the x₁x₃ central plane are shown. We successfully captured the features shown in Fig.10.4.9 in (Joseph, 1992). For 0≤t≤3 the averaged speed of particles is 1.109 and the maximal averaged speed of 60 cylinders is 1.476. But for a single cylinder, the maximal speed is 4.931 for 0≤t≤3. Hence the particle Reynolds numbers based on the length of the cylinders are 44.36, 59.04 and 197.24, respectively. The memory used in the simulation was about 507 MB and it takes about 103.13 s per time step on a Linux-based PC with 2.2 GHz Opteron 248 CPU.

Acknowledgements

We acknowledge the helpful comments and suggestions of R. Bai, E.J. Dean, J. He, H.H. Hu, P.Y. Huang, G.P. Galdi, Y. Kuznetsov and J. Periaux.

---

* Project supported by NSF (Nos. ECS-9527123, CTS-9873236, DMS-9973318, CCR-9902035, DMS-0209066), and DOE/LASCI (No. R71700K-292-000-99), USA