Simulations of settling object using moving domain and immersed-boundary method

Accepted Manuscript

Simulations of settling object using moving domain and immersed-boundary method Sung-Hua Chen, Yen Ku, Chao-An Lin PII: DOI: Reference:

S0045-7930(18)30615-7 https://doi.org/10.1016/j.compfluid.2018.09.007 CAF 4002

To appear in:

Computers and Fluids

Received date: Revised date: Accepted date:

12 March 2018 15 August 2018 6 September 2018

Please cite this article as: Sung-Hua Chen, Yen Ku, Chao-An Lin, Simulations of settling object using moving domain and immersed-boundary method, Computers and Fluids (2018), doi: https://doi.org/10.1016/j.compfluid.2018.09.007

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • An IBM combined with a moving domain is used to simulate settling object.

CR IP T

• The predicted drag coefficients of settling sphere agree reasonably with measurements. • The descending motions of an oblate spheroid are 2-D for the Reynolds numbers investigated • The settling spheroid shows Re dependent steady or fluttering falling.

AC

CE

PT

ED

M

AN US

• The toroidal and hairpin vortices are observed for the settling objects.

1

ACCEPTED MANUSCRIPT

Simulations of settling object using moving domain and immersed-boundary method

CR IP T

Sung-Hua Chen, Yen Ku, Chao-An Lin∗

Department of Power Mechanical Engineering, National Tsing Hua University, Hsinchu 30013, TAIWAN

AN US

Abstract

CE

PT

ED

M

A Cartesian grid-based, immersed-boundary method combined with a moving domain is used to simulate sphere and oblate spheroid settling under gravity. For the moving domain, grids are generated or deleted based on the location of the object centroid, ensuring a constant size of the computational mesh. Sphere settling under gravity is first simulated, and the predicted drag coefficients at different Reynolds numbers, that is, 4.1 to 2000, agree reasonably with measurements. Finally, an oblate spheroid descending under gravity is further simulated with Reynolds numbers at 44.5 to 1617. For the Reynolds numbers examined, the motions are twodimensional, and the fluttering motions exist in all the Reynolds numbers investigated except at Re = 44.5, in which steady falling prevails. This is consistent with previous investigations for falling disks, where the threshold for the steady and fluttering motion is at around Re = 80 for the level of moment of inertia investigated. As the spheroid flutters, different vortex branches are generated in the trailing edge direction, in which hairpin-like vortices occur.

AC

Keywords: immersed boundary method, direct forcing, moving domain, settling object, vortex structure

∗

Corresponding author. Tel:+886-3-5742602 Email address: [email protected] (Chao-An Lin)

Preprint submitted to Computers and Fluids

September 6, 2018

ACCEPTED MANUSCRIPT

1. Introduction

AC

CE

PT

ED

M

AN US

CR IP T

Falling leaves and disks are typical examples of fluid and structure interaction problems. The path instability of a descending object can be due to the interaction of the falling objects wake and geometry. A review of wake-induced oscillatory paths of freely rising or falling bodies in fluids can be found in [1]. Two issues may be considered in numerically investigating the settling object. The first is the embedded objects modelling, and the second is the vast computational domain required to capture the falling bodys transient development after a long-time descent. Although complex moving objects are commonly modelled with bodyfitted or unstructured-grid methods, the transient remeshing further increases computational complexity and hence cost for these approaches. The immersedboundary method (IBM) [2, 3, 4, 5, 6, 7, 8] can be an alternative that not only retains the merit of structured Cartesian grids but also deals with complex geometry, where an external force field is generated along the immersed boundary to mimic the solid borders effect on the fluid. Among different forcing strategies [5], the direct-forcing strategy [6, 7] does not degrade the stability constraint of the time integration scheme. In particular, the direct-forcing approach that explicitly enforces boundary conditions at Eulerian nodes near the immersed boundary [3, 8, 9, 10] provides a better interface representation for many hydrodynamic and thermal problems with time-evolving embedded objects. On the other hand, for a gravity-induced descending object, a significant computational domain is required to capture the transient development after a long-time descent in the object-proceeding direction. Many methods have been adopted to shorten the computational domain. For example, using the stabilized space-time finite-element technique and moving mesh [11], Johnson and Tezduyar [12, 13] were able to calculate the sedimentation of spheres in a Newtonian fluid and studied their arrangements under terminal states. Also, Swaminathan et al. [14] investigated the freely falling ellipsoid in an infinitely long tube using arbitrary Lagrangian and Eulerian (ALE)based finite element method. Another approach is a combination of coarse and fine mesh using a reference frame to compute freely falling disks [15]. The 3-D motion of falling disks was also investigated using the reference-frame approach [16]. In [15, 16], body-fitted grids were used to model the disks. 3

ACCEPTED MANUSCRIPT

AN US

CR IP T

To the best knowledge of the authors, the combination of the immersed-boundary method and moving mesh has not been commonly used to model the transient development of a moving object. In the present study, the immersed-boundary method proposed by Liao et al. [8, 17] is used to mimic the embedded objects movement within the flow. A finite mesh moves with the embedded object to examine the transient development of the solid after a long-time descent. For the moving domain, grids are generated or deleted based on the location of the object centroid, thus ensuring the constant size of the computational mesh. Also, a PETSc-based [18] parallel-iterative pressure Poisson solver is adopted to speed up the 3-D computations. The present methodology is first validated with flows induced by a sphere settling under gravity in a small container. Further, the vertical structures of settling spheres at different Reynolds numbers are presented. Finally, the descents of oblate spheroids under gravity at different Reynolds numbers are also simulated to investigate path instability.

M

2. Methodology of the immersed-boundary technique

AC

CE

PT

ED

2.1. Mathematical formulation Consider the problem of a viscous incompressible fluid in a three-dimensional rectangular domain containing an immersed massless boundary, where a noslip boundary condition is imposed on the boundary. The influence of the immersed boundary on the fluid is represented by forces exerted on the fluid and causing it to move by the prescribed velocity (via the no-slip condition) at the immersed boundary. Thus, the governing equations of this fluid-structure interaction system are the usual incompressible N.–S. equations, 1 ∂u + ∇·(uu) = −∇p + ∆u + f , ∂t Re ∇·u = 0 .

(1) (2)

Here, x = (x, y, z), u(x, t) is the fluid velocity with components u, v and w, p(x, t) is the fluid pressure, and Re is the Reynolds number; ∇ and ∆ are the usual Cartesian-coordinate gradient and Laplacian. Note that the discrete momentum forcing f (x, t) is applied at the fluid nodes adjacent to 4

ACCEPTED MANUSCRIPT

mrb

CR IP T

the immersed boundary and the solid nodes in an Eulerian grid system to satisfy the no-slip condition on the immersed boundaries. The dynamics of a rigid body is described by Newton’s law, where the translational velocity vrb and the angular velocity ωrb of the object can be determined as follows: dvrb = Fg + FB + FD = Fnet , dt dω rb Irb =T dt

(3) (4)

AN US

Here, mrb and Irb are the mass and the moment of inertia of the rigid body, and the latter is a second rank tensor. The net imposed force Fnet on the object is the summation of gravitational force Fg , buoyancy FB , and drag R force FD which equals to −ρ fdV ol. T is the torque on the object which f R equals to −ρf r × fdV ol[19, 20], and r is the position vector from the object geometric center to the forcing point where the velocity is to be determined.

AC

CE

PT

ED

M

2.2. Numerical scheme The numerical procedure used herein consists of a finite-volume method discretized in Cartesian coordinates with a staggered-grid arrangement of dependent variables. This procedure is based on the integration of the transport equations over arbitrary control volumes, leading to the conservation of mass and the balance of momentum and any scalar flow property over each such volume. Spatial derivatives are approximated using second-order centered differencing, and a fractional-step (projection) method implemented with a combination of the Adams–Bashforth and Crank–Nicolson methods for advective and diffusive terms, respectively, is used for temporal discretization. Here, the four-step time-advancement formulation[8, 18, 21, 22] is employed for the momentum equations due to its storage economy and straightforward implementation of boundary conditions. The discrete form of this algorithm is embodied in the following sequence of equations: e − un 3 1 u = − ∇h · (uu)n + ∇h · (uu)n−1 ∆t 2 2 1 e ) + f n+1 , − ∇h pn + ∆h (un + u 2Re 5

(5)

ACCEPTED MANUSCRIPT

(6) (7)

CR IP T

e u∗ − u = ∇h pn , ∆t ∇h · u∗ ∆h pn+1 = , ∆t un+1 − u∗ = −∇h pn+1 , ∆t

(8)

AN US

e and u∗ might be viewed as intermediate velocities occurring between where u the time steps n and n + 1 although they are, in fact, non-divergence-free approximations to un+1 . The discrete operators ∇h and ∆h represent the usual centered-difference approximations for the gradient and Laplace operators, respectively. To evaluate Eq. (5), the forcing f n+1 must be determined in advance to satisfy the no-slip boundary condition on any immersed boundary at the advanced time level. This will be treated in the next subsection. The immersed-boundary location at step n+1 is determined by summing all forces on the sphere at previous time step. The translation and rotation (Eqs. (3) and (4)) of the rigid body at each time step can be obtained by the first-order Euler method: Fnnet ∆t mrb n = ω nrb + I−1,n rb T ∆t

M

vn+1 = vnrb + rb

(10)

ED

ω n+1 rb

(9)

CE

PT

It should be noted that other discretization schemes can also be applied. Thus, the velocity, position and angle at an arbitrary point on the rigid body surface or inside the solid is computed as n+1 Umove = vn+1 rb + ωrb × r n+1 xrb = xnrb + Umove ∆t θ n+1 = θ nrb + ω n+1 rb rb ∆t

(11) (12) (13)

AC

The second moment of mass can calculate the moment of inertia about an axis, i.e., Z (14) Ii = ri2 dm

where ri is the perpendicular distance to the specified axis xi . For regular shapes such as a sphere or a spheroid, there is moment of inertia formula available along their geometric principal axis. However, for an ellipsoid, if the object moves in three-dimensional space, then the moment of inertia has to 6

ACCEPTED MANUSCRIPT

be modified, as shown in Fig. 1. Based on θ and x, the non-inertia reference frame (X0 ) can be constructed, and the transformation to the inertia frame is through the transformation matrix M cosθy0 cosθz0

  0 0 M =  −cosθy sinθz  sinθy0

sinθx0 sinθy0 cosθz0 −cosθx0 sinθy0 cosθz0 +cosθx0 sinθz0 +sinθx0 sinθz0 −sinθx0 sinθy0 sinθz0 cosθx0 sinθy0 sinθz0 0 0 +cosθx cosθz +sinθx0 cosθz0 −sinθx0 cosθy0 cosθx0 cosθy0



CR IP T



(15)

PT

ED

M

AN US

where θ 0 = −θ in Eq. 13.

    

CE

Figure 1: Relationship between the inertial reference frame X0 Y0 Z0 and the non-inertial reference frame X 0 Y 0 Z 0 .

AC

If the moment inertia of an object in the non-inertia frame can be expressed in its principal axis, i.e.,   Ix0 x0 0 0 I rb,ni =  0 Iy0 y0 0  (16) 0 0 0 0 Iz z Then, the moment of inertia in the inertia frame can be expressed as, I rb = M I rb,ni M T 7

(17)

ACCEPTED MANUSCRIPT

f n+1 =

b uF − u , ∆t

CR IP T

2.3. Forcing strategies The force applied on the computational node due to the presence of a solid body must be determined at the advanced time level before the solution procedure presented in Eqs. (5)–(8) can be started. This force (per unit mass) is merely a Newtonian acceleration of any node near the solid body. Thus, (18)

M

AN US

b is an estimate of the time-level n + 1 velocity field in a neighborhood where u of the solid boundary in the absence of forcing, and uF is the velocity (at b is estimated the same location) when forcing is taken into account. Here, u with a simple, explicit Adams–Bashforth scheme where Eq. (1) is provisionally discretized explicitly in time to derive the momentum forcing value, as described in Refs.[8], i.e.,  n  1 n 3 b = u − ∆t ∇h · (uu) − ∆h u u 2 Re  n−1  1 n 1 − 2 ∇h · (uu) − ∆h u + ∇h p , (19) Re

AC

CE

PT

ED

The issue of computing the forcing function f of Eq. (18), and hence uF , when the interface does not coincide with grid nodes, is discussed here. For the present applications, sphere or oblate spheroid are considered. Here, we define the inertia coordinate ri and the body coordinate rb relative to the object, and the coordinate transformation matrix is rb = Mi→b ri . The general governing equation of the sphere or oblate spheroid at the body coordinate can be expressed as, (xb )2 (yb )2 (zb )2 + 2 + 2 =1 (20) a2 b a Then, a cubic box extended 0.52 primary diameter from the object center is used to search for the forcing points and solid points. Any point in the cubic box ri has its corresponding body coordinate counterpart, i.e. rb = (xb , yb , zb ). The following function is used to determine whether the node is solid or fluid, i.e., f (rb ) =

(xb )2 (yb )2 (zb )2 + 2 + 2 −1 a2 b c

Thus, 8

(21)

ACCEPTED MANUSCRIPT

f (rb ) =



≤ 0, > 0,

∈ solid region , ∈ fluid region .

M

AN US

CR IP T

Forcing point is defined as the node located in the fluid region and with at least one neighboring solid node. A schematic of forcing points in the fluid region used to calculate uF is shown in Fig. 2, which is carried out via linear interpolation in the present study, but more accurate interpolation methods could be employed if needed. From the figure, forcing nodes are Cartesian grid points outside the solid body, but immediately adjacent to it. The no-slip or Dirichlet condition is imposed at points of the boundary, Γ, of the solid body at points xB lying on the same grid line as the forcing point xF . Then uF at this point is obtained by (linear) interpolation between points xB and xA , with the latter being the next point outward from xF on the grid line containing points xB and xF . Thus, for the u velocity component at the point xF shown in Fig. 2,  
uˆA − uˆB yA − yF , (22) uF = uˆA − yA − yB

PT

ED

This leads to direct calculation of f n+1 via substitution of (19) and (22) into Eq. (18):  n 3 1 uF − un n+1 + ∇h · (uu) − ∆h u f = ∆t 2 Re n−1  1 1 − ∇h · (uu) − ∆h u + ∇pn . (23) 2 Re

AC

CE

The preceding approach works well if the forcing point is associated with only one solid node, as is true for forcing nodes shown in Fig. 2. On the other hand, if the forcing point is associated with two (or more in 3-D) solid nodes, then ambiguity occurs as indicated in the figure. To remove this ambiguity, Liao et al.[8] proposed utilizing simple linear interpolation between the diagonal node xA and the node xB on the boundary, as shown in Fig. 2. 2.4. Moving domain Consider a two-dimensional fixed computational domain with uniform mesh, as shown in Fig. 3, where the centroid of the object rcen is located at 9

ACCEPTED MANUSCRIPT

(a)

(b Solid Domain Fluid Domain

A

forcing nodes

ambiguous node

F3

B ambiguous node

F2 B

B

CR IP T

F1

A

Fluid Domain

Solid Domain

AN US

A

Figure 2: Schematic diagram for the interpolation procedure for (a) 2-D and (d) 3-D geometry.

CE

PT

ED

M

control volume ∆Vij with grid density being (xn ; 1 : N x × 1 : N y). If, at the next time step, the centroid moves out of the control volume to ∆Vi+1,j−1 , then computational domain is readjusted to (xn+1 ; 2 : N x + 1 × 0 : N y − 1). Thus, grids will be generated or deleted based on the location of the object centroid, which could ensure that the size of the computational mesh remains constant. In the upwind direction of the newly created mesh, i.e. AB and BC, zero inlet velocities uuw n−1,n,∗,n+1 = 0 are prescribed. In the downwind udw = direction, i.e. CD and DA, zero Neumann boundary condition, i.e. ∇˜ 0 is imposed, where u ˜ is computed from Eq. (5). To ensure the global continuity, a uniform correction undc to the wall normal velocity is calculated as, Z Z udw + undwc ) · dADA = 0 (24) (˜ udw + undwc ) · dACD + (˜

AC

Here, the corrections at boundaries CD and DA are assumed to be the same. For simplicity, the following can also be adopted, i.e., Z (25) (˜ udw + undwc ) · dACD = 0 Z (˜ udw + undwc ) · dADA = 0 (26) Therefore, the boundary velocities at the downwind boundary are udw ∗,n+1 = 10

ACCEPTED MANUSCRIPT

xn

8

xn+1 C

D 6

rcen 4

AN US

rcen

2

0

CR IP T

u ˜dw +undc and ∇˜ udw = 0 The boundary condition implementation for threedimensional domain is similar.

A

B

4

6

8

10

M

2

ED

Figure 3: Demonstration of the moving domain method.

PT

2.5. Complete solution procedure Suppose n time steps have completed. The solution at the time level n+1 can be obtained by executing the following steps.

AC

CE

1. Determine the immersed-boundary location at step n + 1 with Eqs. (12) and (13). 2. Adjust the computational mesh if needed based on the location of the centroid of the object. 3. Determine the fluid and solid points on Eulerian grid. 4. Establish locations of forcing points. b. 5. Use Eq. (19) to determine u 6. Calculate the forcing term f n+1 from Eq. (23) with uF obtained via linear interpolation using Eq. (22). e implicitly with force f n+1 7. Compute non-divergence-free velocity field u M by solving Eq. (5). 11

ACCEPTED MANUSCRIPT

CR IP T

8. Obtain estimation of the intermediate velocity u∗ from Eq. (6). 9. Compute pressure at time level n + 1, pn+1 , by solving pressure Poisson equation, Eq. (7). 10. Compute fluid velocity un+1 with divergence-free on time level n + 1 via projection, Eq. (8). 3. Numerical results

AC

CE

PT

ED

M

AN US

3.1. Simulation of a sphere settling under gravity The present formulation is applied to predict three-dimensional sphere settling under gravity. Simulations are first used to predict sphere settling in a rectangular container with fluid, where the experimental study of Cate et al. [23] is used to validate the accuracy of the moving domain. The dimensions of the rectangular container are 6.67D × 10.67D × 6.67D and those of the moving computational domain are 6.67D ×6.67D ×6.67D, which is smaller than the experimental setup. The sphere is released at 2.67D from the top of the box, and the Reynolds numbers investigated are 4.1, 11.6, and 32.2, respectively. Vortex shedding would not occur for Re ≤ 200 [24]. At these small Reynolds numbers, previous results [8] show that ∆x/D ≤ 20 is adequate; thus, a uniform grid of 122 × 122 × 122 cells is used. Simulations with a stationary computational domain of 6.67D × 10.67D × 6.67D are also considered to assess the effectiveness of the results using the moving domain. The measurements contrasted with predicted sphere velocity and position are shown in Fig. 4, where good agreement is achieved. Also, both the computational domain and the moving domain produce the same results, indicating the effectiveness of the latter. As noted, terminal velocity is reached for the lower Reynolds number, that is, Re = 4.1, whereas at higher Reynolds number (Re = 32.2), because of the short distance of the container, terminal velocity is not achieved before the sphere impinges on the bottom wall. The movement of the computational domain and vorticity contours at Re = 32.2 can be observed from Fig. 5, where the dashed line is the size of the original container and the solid line the dimension of the moving domain. When the computational domain reaches the bottom container wall, the sphere is allowed to descend within the moving domain. The development of the vortical structure can be clearly observed in Fig. 5, where, as the sphere travels farther downward, the wake region expands since it has not reached the terminal state. 12

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Figure 4: The comparisons of the simulation results using the moving domain and the experimental results from Cate et al. [23]. (a) the falling velocity of the sphere Vcen (b) the trajectory along the y=direction of the centre of the sphere.

Ref erence Wu and Shu (2010) Johnson and Patel (1999) Wang et al. (2017) Present (D/30) Present (D/40) Present (D/50)

Cd 0.8 0.79 0.8 0.793 0.789 0.788

Table 1: Comparison of drag coefficients Cd for flow over a sphere at Re=200

13

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

Figure 5: The isosurface and slices of instantaneous vortices of the moving domain using vorticity magnitude (|ωv |) at Re=32.2. (a) t=0.51 s (b) t=1.024 s.

14

ACCEPTED MANUSCRIPT

M

AN US

CR IP T

A grid sensitivity test is performed at Reynolds number 200 to examine the influence of grid density on the predicted drag coefficient. Three different grids are adopted, namely, ∆x = D/30, ∆x = D/40 and ∆x = D/50. Table 1 shows the predicted drag coefficients. Despite the variations in grid density, the difference is within 1%. The present predictions are also compatible with those in the literature [25, 26]. Also, in the experimental study of Cate et al. [23], domain size is constant. However, for an unbounded descending object, the size of the moving domain is an important parameter for evaluation. Fig. 6(a) shows the computational domain with an embedded object, where the object centroid is at the center of the horizontal plane (X − Z), and the distances to the vertical (Y ) inlet and exit planes are 1 and 3 Ly , respectively. Three different domain sizes are tested Lx : Ly : Lz = 6.67D : 6.67D : 6.67, 6.67D : 13.3D : 6.67D and 13.3D : 13.3D : 13.3D, respectively, at ∆x = D/40. A spheres descent under gravity is simulated, and the Reynolds number is Re ∼ 2000. Fig. 6(b) shows the descending velocity of the sphere using different domain sizes, and Lx : Ly : Lz = 6.67D : 13.3D : 6.67D and 13.3D : 13.3D : 13.3D produce the same results. Thus, ∆x = D/40 and 6.67D × 13.3D × 6.67D computational domain are used to simulate flows with the descending sphere. Y

X

AC

CE

PT

ED

Z

Figure 6: (a) The domain setting. (b) The time history of the y velocity

Fig. 7 shows the comparisons of the predicted drag coefficients with measurements at different Reynolds numbers, that is, 4.1 to 2000. Symbols are predictions while the solid line shows the measurements. Although some discrepancies are found, reasonably good results are obtained over the investi15

AN US

CR IP T

ACCEPTED MANUSCRIPT

M

Figure 7: The comparisons of the simulation results with the experimental results.

AC

CE

PT

ED

gated range of Reynolds numbers. Detailed terminal-state vortical structures are demonstrated here at two Reynolds numbers (75 and 2000), where the former and the latter are in the steady and unsteady states, respectively. Figs. 8 and 9 show the iso-surface and vorticity contours at two different Reynolds numbers, where unsteady motion at Re = 2000 can be observed. For Re = 75, a regular vortical structure resides behind the sphere, showing a steady vortex tube. On the other hand, at Re = 2000, because of vortex shedding and a higher Reynolds number, different scales of vortical structures exist at Re = 2000. At this Reynolds number, the vortex structure at a high Reynolds number is more complicated than at a low one, where vortices remain behind the sphere, delaying the separation behind it. To further visualize the vortical structure, the flows are presented using the λ2 -criterion [27, 28], where λ2 < 0 represents the vortex core. Figs. 10 and 11 show the iso-surface λ2 at two different Reynolds numbers. At Re = 75, as shown in Fig. 10, the hollow space right behind the sphere is the lowpressure region attributed to the flow separation. The iso-surface of the λ2 coloring with the vertical velocity shows a clear vortex tube behind the sphere. The petal- and ring-shaped vortex cores are typical in the region of 16

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

Figure 8: The iso-surface (a) and contours kωv k = 20 (b) of vorticity of a falling sphere at Re = 75.

17

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

Figure 9: The iso-surface (a) and contours kωv k = 100 (b) of vorticity of a falling sphere at Re = 2000.

18

AN US

CR IP T

ACCEPTED MANUSCRIPT

M

Figure 10: The iso-surface λ2 = −100 of a falling sphere at Re = 75 (color: vertical velocity).

AC

CE

PT

ED

the low Reynolds number flow . Fig. 11 shows the instantaneous iso-surface of λ2 = −2000 at Re = 2000 coloring with the vertical velocity. Toroidal vortices VT resulting from the KelvinHelmholtz instability here form around the sphere edges. From the horizontal plane, the vortices are concentrated in the four corners around the sphere and gradually break up to the ear-shaped vortices VE and then the four hairpin-like vortices VH , which still have small rings around them. The legs of the hairpin vortices are immersed in the inner part of the wake. Other similar hairpin vortices are also observed by Shenoy and Kleinstruer [29] at Re = 180 and Tian et al. [30] at Re=1.5×104 . Because of viscous dissipation, the vortices gradually diminish as the starfishshaped VS and finally the ring-shaped VR . In the last stage of the vortices, the energy on the vortex decays, and hence, these separated vortices reconnect and form a small ring.

3.2. Falling of an oblate spheroid Here, an oblate spheroid descending under gravity (negative ydirection) is simulated, and the aspect ratio is 6. The Reynolds ρ V a numbers defined as f µavg are 44.5 to 1617. Vavg is the average descending velocity at the last two fluttering cycles, and a is the 19

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 11: The iso-surface λ2 = −2000 of a falling sphere at Re = 2000 (color: vertical velocity).

AC

CE

PT

ED

M

radius of the oblate spheroids major axis. The size of the computational domain is 14a(x)×10a(y)×14a(z) with a grid size of ∆x = a/30. Walls along the X and Z directions are stationary. As the computational domain descends in the Y direction, zero momentum is prescribed at the inlet and a Neumann boundary with global mass correction imposed at the exit. The oblate spheroid is allowed to move in the X − Z plane. However, the computational domain moves in the Y direction to ensure that the bottom boundary of the spheroid centroids residing control volume is 2a above the bottom of the computational domain. The simulation starts from a stationary spheroid at the center of the X − Z plane with a π/3 inclined angle. Fig. 12[31, 32, 33] shows the phase diagram indicating the likely pattern of the falling ellipsoid. Note that the phase diagram is for a falling disk. Here, I ∗ is defined as Iellipsoid /ρf d5 . Depending on I ∗ and Re, the path can be steadily falling, fluttering, chaotic, and tumbling. The triangles represent the cases examined. Figs. 13(a) and (b) show the predicted paths of the falling ellipsoid at different Reynolds numbers. For the Reynolds number examined, the motion is two-dimensional as demonstrated in Fig. 13(a), which is consistent with the previous observation of the falling disk. Also, the fluttering motion prevails in all Reynolds

20

Figure 12: Phase diagram.

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

21

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 13: Trajectories of the center of the ellipsoid at different Reynolds numbers.

22

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 14: Transient development of the falling ellipsoid at three Reynolds numbers.

23

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

Figure 15: Transient development of vortical structure at three Reynolds numbers-isosurfaces of λ2 = −3000..

24

ACCEPTED MANUSCRIPT

M

AN US

CR IP T

numbers investigated except Re = 44.5. This phenomenon is also consistent with the results from the phase diagram shown in Fig. 12, where the critical Reynolds number is about 80 for the I ∗ investigated. Fig. 13(b) also shows that the fluttering motion strengthens as the Reynolds number increases. Fig. 14 shows a 3-D view of the falling oblate spheroid at three Reynolds numbers, where the steady falling (Fig. 14(a)) and fluttering motions (Figs. 14(b) and (c)) exist. The horizontal movements of the spheroid have been demonstrated in Figs. 13(a) and (b). On the other hand, snapshots of the iso-surface of λ2 are shown in Figs. 15(a)-(c). At Reynolds number Re = 44.5 shown in Fig. 15(a), steady falling produces a stable toroidal vortex around the spheroid edges. Figs. 15(b) and (c) indicate that the spheroid flutters and that the hairpin-like vortices appear when the Reynolds number increases. During fluttering, different vortex branches develop at the downwind region, and a complex hairpin-like vortical structure exists as shown in Fig. 15(c) at a higher Reynolds number. 4. Conclusion

AC

CE

PT

ED

In the present study, three-dimensional rigid objects descending under gravity are simulated using the Cartesian grid-based, immersed-boundary method and a moving computational domain. As the centroid of the object moves across different control volumes, the computational domain is readjusted, with grids being generated or deleted based on the location of the object centroid, ensuring a constant size of the computational mesh. Threedimensional sphere settling under gravity is first simulated. The predicted drag coefficients at the terminal state at different Reynolds numbers, that is, 4.1 to 2000, are in reasonably good agreement with measurements. The iso-surface of λ2 is used to visualize the vortical structure at Re = 75 and Re = 2000. Petal- and ring-shaped vortex cores are observed in the flow with a Reynolds number of 75. On the other hand, a more complex vortical structure is observed for Re = 2000, where the KelvinHelmholtz instabilityinduced toroidal vortices and hairpin-like vortices are present. Finally, an oblate spheroid descending under gravity is further simulated at Reynolds numbers 44.5 to 1617. For the Reynolds numbers examined, the motion is two-dimensional, and a fluttering motion 25

ACCEPTED MANUSCRIPT

CR IP T

exists in all the investigated Reynolds numbers except at Re = 44.5. This observation is consistent with those in previous investigations of falling disks, where the threshold for the steady and fluttering motion is at around Re = 80 for the level of moment of inertia investigated. As the spheroid flutters, different vortex branches develop in the trailing direction, where the hairpin-like vortices develop. Also, the previously generated vortices decay because of the lack of energy supply. 5. ACKNOWLEDGMENTS

AN US

This present research was supported by Ministry of Science and Technology of Taiwan under project MOST 105-2221-E-007-016-MY3. We are thankful for the computational facilities and resources provided by the National Centre for High-Performance Computing of Taiwan.

References

M

References

ED

[1] P. Ern, F. Risso, D. Fabre and J. Magnaudet, Wake-induced oscillatory paths of bodies freely rising or falling in fluids, Annu. Rev. Fluid Mech. 44 (2012) 97–121.

PT

[2] C.S. Peskin, Flow patterns around heart valves: a numerical method, J. Comput. Phys. 10 (1972) 252–271.

CE

[3] J. Mohd-Yusof, Combined immersed boundary/B-Spline method for simulations of flows in complex geometries CTR annual research briefs, NASA Ames/Stanford University, 1997.

AC

[4] M.C. Lai, C.S. Peskin, An immersed boundary method with formal second-order accuracy and reduced numerical viscocity, J. Comput. Phys. 160 (2000) 705–719. [5] R. Mittal, G. Iaccarino, Immersed boundary methods, Annu. Rev. Fluid Mech. 37 (2005) 239–261. [6] M. Uhlmann, An immersed boundary method with direct forcing for the simulation of particulate flows, J. Fluid Mech. 209 (2005) 448–476. 26

ACCEPTED MANUSCRIPT

[7] S.W. Su, M.C. Lai, C.A. Lin, A simple immersed boundary technique for simulating complex flows with rigid boundary, Comput. Fluids 36 (2007) 313–324.

CR IP T

[8] C.C. Liao, Y.W. Chang, C.A. Lin, J.M. McDonough, Simulating flows with moving rigid boundary using immersed-boundary method, Comput. Fluids 39 (2010) 152–167.

[9] C.C. Liao, C.A. Lin, Simulations of natural and forced convection flows with moving embedded object using immersed boundary method, Comput. Methods Appl. Mech. Eng. 58 (2012) 213–216.

AN US

[10] C.C. Liao, C.A. Lin, Influences of a confined elliptic cylinder at different aspect ratios and inclinations on the laminar natural and mixed convection flows, Int. J. Heat Mass Transf. 55 (2012) 6638–6650

M

[11] T.E. Tezduyar, M. Behr, J. Liou, A new strategy for finite element computations involvingmoving boundaries and interfacesthe deformingspatial-domain/space-time procedure: I. The concept and the preliminary numerical tests. Comput Methods Appl Mech Eng 94 (1992) 339– 351.

ED

[12] A. Johnson, T.E. Tezduyar, Simulation of multiple spheres falling in a liquid-filled tube. Comput Methods Appl Mech Eng 134 (1996) 351–373.

PT

[13] A. Johnson, T.E. Tezduyar, 3D simulation of fluid-particle interactions with the number of particles reaching 100. Comput Methods Appl Mech Eng 145 (1997) 301–321.

CE

[14] T. N. Swaminathan, K. Mukundakrishnan and H. H. Hu, Sedimentation of an ellipsoid inside an infinitely long tube at low and intermediate Reynolds numbers, J. Fluid Mech. 551 (2006) 357–385.

AC

[15] M. Chrust, G. Bouchet and J. Dusek, Numerical simulations of the dynamics of freely falling discs, Phys. Fluids, 25 (2013) 044102. [16] F. Auguste, J. Magnaudet and D. Fabre, Falling styles of disks, J. Fluid Mech, 719 (2013) 388–405.

27

ACCEPTED MANUSCRIPT

CR IP T

[17] C.C. Liao, C.A. Lin, Influence of Prandtl number on the instability of natural convection flows within a square enclosure containing an embedded heated cylinder at moderate Rayleigh number. Phys of Fluids 27 (2015) 013603. [18] S.W. Hsu, F.N. Hwang, Z.H. Wei, S.H. Lai, C.A. Lin, A parallel multilevel preconditioned iterative pressure Poisson solver for the large-eddy simulation of turbulent flow inside a duct. Comput Fluids 45 (2012) 138–146.

AN US

[19] J. Yang, F. Stern, A simple and efficient direct forcing immersed boundary framework for fluid-strture interactions. J. Comput. Phys. 231 (2012) 5029–5061. [20] C. C. Liao, W. W. Hsiao, Y. T. Lin, C. A. Lin, Simulations of two sedimenting-interacting spheres with different sizes and initial configurations using immersed boundary method, Computational Mechanics, 55 (2015) 1191–1200.

M

[21] S.W. Hsu, J.B. Hsu, W. Lo, C.A. Lin, Large eddy simulations of turbulent Couette-Poiseuille and Couette flows inside a square duct. J Fluid Mech 702 (2012) 89–101.

ED

[22] B. E. Owolabi, C. A. Lin, Marginally turbulent Couette flow in a spanwise confined passage of square cross section, Phys of Fluids 30 (2018) 075102.

CE

PT

[23] A. ten Cate, C.H. Nieuwstad, J.J. Derksen, H.E. A Van den Akker, Particle imaging velocimetry experiments and lattice-Boltzmann simulations on a single sphere settling under gravity. Physics of Fluids 14 (2002) 4012–4025.

AC

[24] T.A. Johnson, V.C. Patel, Flow past a sphere up to a Reynolds number of 300 , J. Fluid Mech. 378 (1999) 19–70. [25] T. A. Johnson, V. C. Patel, Flow past a sphere up to a Reynolds number of 300, J. Fluid. Mech. (1999), vol. 378, pp. 19-70. [26] Y. Wang, C. Shu, L. M. Yang, C. J. Teo, An immersed boundary-lattice Boltzmann flux solver in a moving frame to study three-dimensional 28

ACCEPTED MANUSCRIPT

freely falling rigid bodies, Journal of Fluids and Structures (2017), vol. 68, pp. 444-465.

CR IP T

[27] J.C.R. Hunt, A.A. Wray, P. Moin, Eddies, Stream, and Convergence Zones in Turbulent Flows. Proceedings of the 1988 CTR Summer Program, NASA Ames/Stanford University, Stanford, CA, (1988) 193–208. [28] J. Jeong, F. Hussain, On the Identification of a Vortex, J. Fluid Mech. 285 (1995) 69–94.

AN US

[29] A. R. Shenoy, C. Kleinstruer, Flow over a thin circular disk at low to moderate Reynolds numbers, J. Fluid Mech (2008), vol 605, pp. 253-262. [30] X. Tian, M. C. Ong, J. Yang and D. Myrhaug, Large-eddy simulations of flow normal to a circular disk at Re=1.5 × 105 , Computers and Fluids (2016), vol 140, pp. 422-434.

M

[31] W. W. Willmarth, N. E. Hawk, R. L. Harvey, Steady and Unsteady Motions and Wakes of freely falling disks, Phys. Fluids, 7 (1964) 197– 208.

ED

[32] G. E. Stringham, D. B. Simons, H. P. Guy, The Behavior of Large Particles Falling in Quiescent Liquids, Prof. Pap. US Geol. Surv. 562-C (1969).

AC

CE

PT

[33] S. Field, M. Klaus, M. G. Moore, Franco Nori, Chaotic dynamics of falling disks, Nature., 388 (1997) 17.

29