Numerical Investigation of Freely Falling Objects Using Direct-Forcing Immersed Boundary Method

Abstract

The fluid-structure interaction of solid objects freely falling in a Newtonian fluid was investigated numerically by direct-forcing immersed boundary (DFIB) method. The Navier–Stokes equations are coupled with equations of motion through virtual force to describe the motion of solid objects. Here, we rigorously derived the equations of motion by taking control-volume integration of momentum equation. The method was validated by a popular numerical test example describing the 2D flow caused by the free fall of a circular disk inside a tank of fluid, as well as 3D experimental measurements in the sedimentation of a sphere. Then, we demonstrated the method by a few more 2D sedimentation examples: (1) free fall of two tandem circular disks showing drafting, kissing and tumbling phenomena; (2) sedimentation of multiple circular disks; (3) free fall of a regular triangle, in which the rotation of solid object is significant; (4) free fall of a dropping ellipse to mimic the falling of a leaf. In the last example, we found rich falling patterns exhibiting fluttering, tumbling, and chaotic falling.

1. Introduction

The falling of light objects, such as feathers or leaves in fluids, driven by gravity and hydrodynamic force, poses an intractability in fluid dynamics. The embedded nonlinearity often results in complicated falling patterns and chaotic trajectories. Several experimental and numerical works have been conducted to study the extensive dynamics of freely falling objects. Experiments conducted by Belmonte et al. [1] and Mahadevan et al. [2] focused on thin flat strips and tumbling cards falling through a fluid in a vertical cell. Through experiments and numerical simulations of 2D incompressible Navier–Stokes equations, Pesavento and Wang [3] analyzed the aerodynamics of a freely falling ellipse, a setup akin to a leaf or business card falling in air.

The investigation into the unsteady aerodynamics of fluttering and tumbling plates formed the core work of Andersen et al. [4]. By varying the thickness-to-width ratio and the dimensionless moment of inertia, they were able to capture the transitions among fluttering, tumbling, and steady fall for the free fall of a plate. Jin and Xu [5] studied the unsteady aerodynamics associated with freely falling elliptical and rectangular plates both experimentally and numerically. In their discovery, they showed that the difference in trajectories between elliptical and rectangular plates was a nuance. The distinction, attributable to the geometry of plates, was observed in the angular velocity, in which the rectangular plate was found to rotate at a much slower rate than the elliptical one. Though the falling of solids with a slender shape, like ellipse and rectangular plate exhibits richer falling patterns, here, we managed to study this sedimentation hydrodynamics staring from a circular disk to a triangle and, finally, to an ellipse by decreasing the symmetry in shape using our long-term developing direct-forcing immersed boundary (DFIB) method [6,7,8]. To avoid repetition, the fundamental idea of DFIB method, together with its benchmark testing and error analysis, is referred to Reference [6,7,8]. Besides free falling of a single solid object, sedimentation of multiple solid particles has long been an interesting and fundamental subject bearing many significant industrial applications, such as paper making, specialty chemicals, petroleum, bioengineering, pharmaceuticals, biomass gasification, and combustion. Even today, it is still under extensive and intensive study theoretically, numerically, and experimentally [9,10,11,12,13,14]. This motivated us to study sedimentation of multiple circular disks here in order to understand interaction dynamics among solids, besides fluid-solid interaction, during settling.

The rest of this paper is organized as follows. The basic idea of DFIB approach for solving fluid-structure-interaction (FSI) problems is described in Section 2 with the coupled equations of motion for a freely falling solid object rigorously derived by control-volume integration of momentum equations. In Section 3, we validated DFIB method by a popular numerical test example describing the 2D flow caused by the free fall of a circular disk inside a tank of fluid [15], as well as 3D experimental sedimentation measurements of a sphere [16]. In Section 4, sedimentation of two circular disks in tandem and multiple circular disks in arrays, showing drafting, kissing, and tumbling phenomena, was computed to demonstrate our DFIB method. Section 5 shows the computation of free fall of a regular triangle, where the rotation is as significant as displacement in the motion of solid due to the breaking of rotational symmetry possessed by a circular disk. The falling pattern starts to show minor fluttering and tumbling behaviors. The rotational effect is further studied by the free fall of a dropping ellipse to mimic the falling of a leaf. We found its falling patterns exhibiting fluttering, tumbling, and chaotic falling, which are particularly sensitive to the aspect ratio of ellipse and density ratio between solid and fluid. Finally, a conclusion is given at end.

2. Mathematical Model and Numerical Method

In the present work, we employed the DFIB approach developed in Reference [6,7,8] for all 2D and 3D FSI computations. A virtual force is added to the incompressible Navier–Stokes equations in order to accommodate the interaction between solid and fluid. A solid body, immersed in fluid, is identified by a volume-of-solid (VOS) function, $\eta$, which denotes the volume fraction of solid within a numerical cell. For a cell full of solid, $\eta$ is equal to 1, while it becomes 0 for a fluid-filled cell, as shown in Figure 1. It would be fractional in a cut cell, consisting of both solid and fluid. With this notation, we describe DFIB method as follows.

2.1. DFIB Method

The spirit of DFIB method is treating a rigid solid object immersed in fluid same as ambient fluid but with a prescribed velocity following equations of motion. This idea is achieved by applying an extra virtual force in momentum equations to enforce the fluid within the solid domain to follow prescribed solid-object velocity [7,8]. The mechanics would then be totally equivalent to a rigid solid object interacting with its surrounding fluid. By doing so, we can avoid using traditional body-fitted methods, such as arbitrary Lagrangian-Eulerian (ALE) method [17,18,19], which is generally less easy to implement. More specifically, the incompressible Navier–Stokes equations under this framework are expressed as

$$\begin{array}{cc}\hfill \nabla ·\mathit{u}& =0,\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill \frac{\partial \mathit{u}}{\partial t}+\nabla ·\left(\mathit{u}\mathit{u}\right)& =-\frac{1}{{\rho }_f}\nabla (p+{\rho }_{f}g\mathit{z})+\nu {\nabla }^2\mathit{u}+\eta \mathit{f},\hfill \end{array}$$

(2)

where u, ${\rho }_f$, p, and $\nu$ are velocity, density, pressure, and kinematic viscosity of fluid; g is the gravitational acceleration. Note that, in Equation (2), the virtual force f is exerted only on the solid domain by volume-of-solid function $\eta$.

Traditionally, a time-splitting scheme together with pressure increment projection method, implemented by finite difference method in staggered grids, is used to advance velocity from nth to $(n+1)$th time level with ${\mathit{u}}^{n+1}$ satisfying Equation (1). Here, in DFIB method [6,7,8], we advance velocity from ${\mathit{u}}^n$ to the intermediate time level ${\mathit{u}}^{**}$, instead of ${\mathit{u}}^{n+1}$, with ${\mathit{u}}^{**}$ satisfying Equation (1). This is because velocity in solid domain, ${\Omega }_s$, has not yet satisfied prescribed solid object velocity. To do so, virtual force is recruited for the advancement of ${\mathit{u}}^{**}$ to ${\mathit{u}}^{n+1}$ such that u in ${\Omega }_s$ would comply with prescribed solid object velocity,

$$\frac{{\mathit{u}}^{n+1}-{\mathit{u}}^{**}}{\Delta t}=\eta {\mathit{f}}^{n+1/2}.$$

(3)

Note that, with the virtual force accompanied by VOS, $\eta$, in Equation (3), we know this update of velocity to ${\mathit{u}}^{n+1}$ would only be effective for solid domain. This leaves ${\mathit{u}}^{n+1}={\mathit{u}}^{**}$, satisfying divergence-free condition, for fluid domain. With prescribed velocity for the solid object, denoted by ${\mathit{u}}_s$, known in advance through equations of motion shown later, we can reciprocally calculate the virtual force by letting ${\mathit{u}}^{n+1}={\mathit{u}}_s^{n+1}$ in ${\Omega }_s$ in Equation (3) and obtain

$$\eta {\mathit{f}}^{n+1/2}=\frac{{\mathit{u}}^{n+1}-{\mathit{u}}^{**}}{\Delta t}=\frac{{\mathit{u}}_s^{n+1}-{\mathit{u}}^{**}}{\Delta t}.$$

(4)

Physically, the negative of volume integration of virtual force, $\eta \mathit{f}$, would be the resultant force exerted on the solid object by fluid:

$$\mathit{F}=-{\int \int \int }_{\Omega }\eta {\rho }_f\mathit{f}dV.$$

(5)

This computation of resultant force would request complicated surface integration of pressure and viscous stress if we otherwise use traditional body-fitted methods.

2.2. Equations of Motion Governing the Movement of an Immersed Solid

The prescribed velocity at each solid cell ${\mathit{u}}_s$ in Equation (4) has to follow the equations of motion governing the free-fall dynamics of solid objects. The motion of a solid object is tracked in a Lagrangian frame by linear and angular momentum equations, that delivers center-of-mass and angular velocities for solid object. That means we can decompose ${\mathit{u}}_s$ at each solid cell into rigid-body translational and rotational components as follows.

$${\mathit{u}}_s={\mathit{v}}_s+{\omega }_s\times \mathit{r},$$

(6)

where ${\mathit{v}}_s$ is the center-of-mass velocity of solid object, and ${\omega }_s$ is its angular velocity around a rotational axis through center of mass. r is the position vector of a solid cell with respect to center of mass. ${\mathit{v}}_s$ and ${\omega }_s$ further provide position and rotating angle of the solid object by time integration.

To derive an equation for ${\mathit{v}}_s$, a control volume integration over ${\Omega }_s$ is taken on Equation (2). Note that the control volume integration of pressure gradient and viscous term over ${\Omega }_s$ would be equivalent to the surface integration of stress $\tau$ over $\Gamma$ by divergence theorem, and it ends up as

$$\frac{d{\mathit{v}}_s}{dt}{\int \int \int }_{{\Omega }_s}{\rho }_{f}dV={\int \int }_{\Gamma }\tau ·\mathit{n}d\mathit{A}+{\int \int \int }_{{\Omega }_s}{\rho }_f\mathit{g}dV+{\int \int \int }_{{\Omega }_s}{\rho }_f\mathit{f}dV.$$

(7)

With regard to the solid object, its equation of motion for translation takes the form as

$$\frac{d{\mathit{v}}_s}{dt}{\int \int \int }_{{\Omega }_s}{\rho }_{s}dV={\int \int }_{\Gamma }\tau ·\mathit{n}d\mathit{A}+{\int \int \int }_{{\Omega }_s}{\rho }_s\mathit{g}dV.$$

(8)

Since the surface drag term ${\int \int }_{\Gamma }\tau ·\mathit{n}d\mathit{A}$ in Equations (7) and (8) is usually unknown a priori, taking a difference between Equations (7) and (8) would eliminate it and result in

$$\frac{d{\mathit{v}}_s}{dt}{\int \int \int }_{{\Omega }_s}{\rho }_{s}dV={\int \int \int }_{{\Omega }_s}({\rho }_s-{\rho }_f)\mathit{g}dV-{\int \int \int }_{{\Omega }_s}{\rho }_f\mathit{f}dV+\frac{d{\mathit{v}}_s}{dt}{\int \int \int }_{{\Omega }_s}{\rho }_{f}dV,$$

(9)

where gravity and buoyancy are denoted by the first term of the right hand side, and hydrodynamic drag is represented by the negative of virtual force in the second term. The last term denotes a fluid inertia term. We can further express Equation (9) as

$$m_s\frac{d{\mathit{v}}_s}{dt}=(m_s-m_f)\mathit{g}-{\int \int \int }_{\Omega }\eta {\rho }_f\mathit{f}dV+m_f\frac{d{\mathit{v}}_s}{dt},$$

(10)

where $m_s$ and $m_f$ denote the mass of solid object and and its replacement by fluid, respectively, with expressions

$$m_s={\int \int \int }_{{\Omega }_s}{\rho }_{s}dV={\int \int \int }_{\Omega }\eta {\rho }_{s}dV$$

(11)

and

$$m_f={\int \int \int }_{{\Omega }_s}{\rho }_{f}dV={\int \int \int }_{\Omega }\eta {\rho }_{f}dV.$$

(12)

We can then discretize Equation (10) in time with ${\mathit{f}}^{n+1/2}\approx \frac{3}{2}{\mathit{f}}^n-\frac{1}{2}{\mathit{f}}^{n-1}$, which is equivalent to employ 2nd-order-accurate Adams-Bashforth scheme for the virtual force term, and the equation becomes

$$m_s\frac{{\mathit{v}}_s^{n+1}-{\mathit{v}}_s^n}{\Delta t}=(m_s-m_f)\mathit{g}-\left({\displaystyle \frac{3}{2}}{\int \int \int }_{{\Omega }_s}{\rho }_f{\mathit{f}}^{n}dV-{\displaystyle \frac{1}{2}}{\int \int \int }_{{\Omega }_s}{\rho }_f{\mathit{f}}^{n-1}dV\right)+m_f\frac{{\mathit{v}}_s^n-{\mathit{v}}_s^{n-1}}{\Delta t}.$$

(13)

Note that the fluid inertia term is discretized by lagging with a time step on purpose. Without this lagging, the initial acceleration felt by solid at a quiescent environment would be g instead of the correct $\frac{m_s-m_f}{m_s}g$ on account of buoyancy. This lagging is comprehensible, since flow motion is driven by movement of solid first, as described by Equation (8), and then the response of fluid follows by Equation (7). Equation (10) can be rewritten in short as

$$\frac{d\left(m_s{\mathit{v}}_s\right)}{dt}={\mathit{f}}_{gb}+{\mathit{f}}_d+\frac{d\left(m_f{\mathit{v}}_s\right)}{dt},$$

(14)

where ${\mathit{f}}_d$ and ${\mathbf{f}}_{gb}$ represent the hydrodynamic drag and gravity/buoyancy, respectively.

Similarly, to derive an equation for ${\mathbf{\omega }}_s$, we can first take a moment around the axis of rotation through center of mass on Equation (2). Following similar procedures as Equations (7)–(10), we can obtain

$$\frac{d({\mathit{I}}_s·{\mathbf{\omega }}_s)}{dt}={\mathit{T}}_s+\frac{d({\mathit{I}}_f·{\mathbf{\omega }}_s)}{dt},$$

(15)

where ${\mathit{I}}_s$ and ${\mathit{I}}_f$ are the moment of inertia for the solid object and its fluid replacement, respectively. The torque ${\mathit{T}}_s$ acting upon the solid is determined by

$${\mathit{T}}_s=-{\int \int \int }_{\Omega }\eta {\rho }_f\mathit{r}\times \mathit{f}dV.$$

(16)

Equation (16) can be further discretized similar to Equation (13). Once ${\mathit{v}}_s$ and ${\omega }_s$ are found, ${\mathit{u}}_s$ can then be determined by Equation (6). The trajectory of every solid cell can be further determined by

$$\frac{d{\mathit{X}}_s}{dt}={\mathit{v}}_s,$$

(17)

and then ${\Omega }_s\left(t\right)$ can be determined, as well.

2.3. Collision Model

When we compute the free fall of multiple solid objects immersed in fluid, the collision among solid objects or between solid objects and wall are often inevitable. In computation, the interference (overlapping) among solid objects or between solid objects and wall during falling could happen, though in reality this is often prevented by a thin in-between lubrication fluid layer. Unless our mesh can resolve this thin lubrication layer, which would request a rather exhaustive resolution, collision model is necessary in practice to avoid possible interferences. Glowinski et al. [15] proposed a repulsive force model, in which an artificial short-range repulsive force was introduced to make sure the separation among solid objects and between solid object and wall. Here, we introduce a simple collision model by designing a repulsion force, which would be activated when body-body or body-wall collision is in progress.

For simplicity, we only consider solid objects as circular disks here. Assuming there are $N_p$ disks immersed in the fluid and labeled as $1,2,\cdots ,N_p$, we denote the total repulsive force exerted on disk i by

$${\mathit{F}}_i^{co}:={\mathit{F}}_i^p+{\mathit{F}}_i^w:=\sum _{j=1,j\ne i}^{N_p}{\mathit{F}}_{ij}^p+\sum _{j=1}^{N_w}{\mathit{F}}_{ij}^w,$$

(18)

where ${\mathit{F}}_i^p$ is the repulsive force between disk i and the others, and ${\mathit{F}}_i^w$ is the repulsive force between disk i and walls (labelled as $1,2,\cdots ,N_w$). ${\mathit{F}}_{ij}^p$ is the disk-to-disk repulsive force that disk i bears from disk j. Similarly, ${\mathit{F}}_{ij}^w$ is the disk-to-wall repulsive force that disk i bears from wall j.

In the present model, the repulsive force is activated whenever a disk is close enough to another disk or walls by the introduction of small positive parameters ${\epsilon }_p,{\epsilon }_w$ to control the size of the repulsive force and a small distance tolerance $\delta >0$ for the activation and adjustment of the force. More specifically, we define the disk-to-disk repulsive force ${\mathit{F}}_{ij}^p$, acting along center-to-center direction on disk i due to the collision with disk j, by [20]

$${\mathit{F}}_{ij}^p=\left\{\begin{array}{cc}\mathbf{0},\hfill & d_{ij}>R_i+R_j+\delta ,\hfill \\ {\displaystyle \frac{c_{ij}}{{\epsilon }_p}}\left({\displaystyle \frac{{\mathit{X}}_c^{\left(i\right)}-{\mathit{X}}_c^{\left(j\right)}}{d_{ij}}}\right){\left({\displaystyle \frac{R_i+R_j+\delta -d_{ij}}{\delta }}\right)}^2,\hfill & d_{ij}\le R_i+R_j+\delta ,\hfill \end{array}\right.$$

(19)

where ${\mathit{X}}_c^{\left(i\right)}$ and $R_i$ are center coordinates and radius of disk i; $d_{ij}$ is the center-to-center distance between disks i and j; $c_{ij}$ is a scaling factor with a dimension of force. In this simulation, $c_{ij}$ is chosen to be the buoyant force acting on disk i. Similarly, the repulsive force with walls ${\mathit{F}}_{ij}^w$, acting normal to wall j and through center of disk i is given by [20]

$${\mathit{F}}_{ij}^w=\left\{\begin{array}{cc}\mathbf{0},\hfill & d_{ij}^{\prime }>2R_i+\delta ,\hfill \\ {\displaystyle \frac{c_{ij}}{{\epsilon }_w}}\left({\displaystyle \frac{{\mathit{X}}_c^{\left(i\right)}-{\mathit{X}}_{ij}^{\prime }}{d_{ij}}}\right){\left({\displaystyle \frac{2R_i+\delta -d_{ij}^{\prime }}{\delta }}\right)}^2,\hfill & d_{ij}^{\prime }\le 2R_i+\delta ,\hfill \end{array}\right.$$

(20)

where ${\mathit{X}}_{ij}^{\prime }$ is the center coordinates of disk i’s mirror image with respect to the wall j; $d_{ij}^{\prime }$ is the center-to-center distance between disk i and its mirror image against wall j.

2.4. Numerical Procedure

The numerical procedure used to compute fluid flow and trace the position and orientation of falling object at each time step can be summarized in the following steps:

• Calculate ${\mathit{v}}_s$ and ${\omega }_s$ via Equations (14) and (15) and then determine ${\mathit{u}}_s$ for each solid cell via Equation (6).
• Integrate ${\mathit{v}}_s$ and ${\omega }_s$ to obtain center-of-mass position ${\mathit{X}}_s$ and angle of rotation for solid object and then determine solid domain ${\Omega }_s$.
• Determine volume of solid function $\eta$ through ${\Omega }_s$.
• Calculate ${\mathit{u}}^{**}$ via projection method.
• Advance to ${\mathit{u}}^{n+1}$ and compute the virtual force by ${\mathit{u}}_s$ via Equations (3) and (4).

3. Numerical Validations

This section presents benchmark testing of current method by sedimentation of a circular disk for 2D flows [15,20,21] and a sphere for 3D flows [16].

3.1. Sedimentation of a Circular Disk

The numerical simulation describing sedimentation of a circular disk released from rest inside a rectangular enclosure was conducted following the physical setting and geometric configuration of Glowinski et al. [15]. More specifically, the computational domain is $\Omega =(0\mathrm{cm},2\mathrm{cm})\times (0\mathrm{cm},6\mathrm{cm})$ and the disk is located at $(1\mathrm{cm},4\mathrm{cm})$ initially. The densities of fluid and solid are ${\rho }_f=1$ g/cm${}^3$ and ${\rho }_s=1.25$ g/cm${}^3$, respectively. The diameter of the disk is $D=0.25\mathrm{cm}$ and the kinematic viscosity of the fluid is $\nu =0.1 {\mathrm{cm}}^2$/s. The gravity is $g=-981$ cm/s${}^2$. In current simulations, the computational domain is discretized by a uniform mesh with the mesh size $\Delta x=\Delta y=0.025$ cm. To capture the trajectory and orientation of a free-fall solid object accurately, we deliberately choose a rather small time step in current study, which is much smaller than the time step requested by CFL constraint. For the current case, we chose $\Delta t={10}^{-4}$ s.

Simulated flow field featured by vorticity contours at various times is shown in Figure 2. Time traces of disk’s vertical position and velocity were particularly calculated and shown in Figure 3. Figure 3a,b shows the grid-independence of current computations. Basically, results under different mesh sizes coincide with each other except the velocity near the hitting of floor in Figure 3b. Here, we further adopt the result with the finest mesh size in Figure 3a,b to compare with [15,20,21] and show it in Figure 3c,d. Our result generally bears a good agreement with references mentioned above, especially with Glowinski et al. [15]. This can be told from the immediate velocity near the hitting of floor in Figure 3d, where actually none of [15,20,21] is close to one another. From Figure 3d, the disk is first accelerated by gravity, and then the fluid drag increases as disk’s velocity increases. When gravity is finally balanced by drag, the velocity of disk reaches its maximum terminal velocity, which is steady and not altered until hitting the floor. Note that the peak followed by oscillations in Figure 3d, happening near the hitting of floor in our computation, is because of the usage of collision model, described in Section 2.3. The peak and following oscillations imply the damping rebounds from the floor. The collision model is a necessary evil. Without it, falling solid object would generally penetrate the floor numerically.

3.2. Sedimentation of a 3D Sphere

Here, the numerical simulation of sedimentation of a 3D sphere inside a fluid tank was conducted with physical parameters and geometry settings same as the experiment reported by Cate et al. [16]. The computational domain is set to $\Omega =(-5\mathrm{cm},5\mathrm{cm})\times (-5\mathrm{cm},5\mathrm{cm})\times (-12\mathrm{cm},4\mathrm{cm})$ based on the rectangular fluid tank of size $10\mathrm{cm}\times 10\mathrm{cm}\times 16\mathrm{cm}$ mentioned in experiment. A sphere of diameter $D=1.5$ cm is initially centered at $(0\mathrm{cm},0\mathrm{cm},0.75\mathrm{cm})$ and released from rest with its density being ${\rho }_s=1120$ kg/m${}^3$, slightly larger than the density of ambient fluid. Following the experiment, various fluid density, ${\rho }_f$, and dynamic viscosity, ${\mu }_f$, were employed for simulations and are listed in

Table 1. Physical parameters used to simulate the sedimentation of a 3D sphere based on [16].

	${\mathit{\rho }}_{\mathit{f}}$ [kg/m${}^3$]	${\mathit{\mu }}_{\mathit{f}}$ [Ns/ m${}^2$]	${\mathit{U}}_{\mathit{\infty }}$ [m/s]	$\mathit{Re}$
Case 1	970	3.73	3.8	1.5
Case 2	965	2.12	6.0	4.1
Case 3	962	1.13	11.5	11.6
Case 4	960	0.58	12.8	32.2

with Reynolds number defined as $Re={\rho }_f{\mathit{U}}_{\infty }D/{\mu }_f$, in which ${\mathit{U}}_{\infty }$ is the terminal velocity. We conducted the numerical simulations in a uniform Cartesian grid with $\Delta x=\Delta y=\Delta z=0.05$ cm, time step $\Delta t={10}^{-5}$ s, and the results are depicted in Figure 4, where the time traces of sphere’s vertical position and velocity for each case are shown in Figure 4a,b, and a snap shot of vorticity magnitude contours for the falling sphere is shown in Figure 4c. We find that our numerical results are in excellent agreements with the experimental data in Reference [16].

4. More Sedimentation Simulations by DFIB Method

Here, we further computed the sedimentation of (1) two tandem circular disks and (2) multiple circular disks to demonstrate current method.

4.1. Sedimentation of Two Tandem Circular Disks

The sedimentation of two tandem circular disks inside a rectangular enclosure is studied here with the densities of disk and fluid being ${\rho }_s=1.5$ g/cm${}^3$ and ${\rho }_f=1$ g/cm${}^3$, respectively. The kinematic viscosity of fluid is set as $\nu =0.01$ cm${}^2$/s. The rectangular enclosure is 2 cm wide and 6 cm long, and the diameter of the disks is $D=0.25$ cm. These two tandem disks are located initially at the top of the rectangular enclosure with the distance between centers of disks being 0.5 cm, and so the gap between two disks is 0.25 cm. A uniform mesh of size $\Delta x=\Delta y=1/128$ cm, and time step $\Delta t={10}^{-5}$ s are used in current simulation. The collision parameter are set as ${\epsilon }_w={\epsilon }_p={10}^{-5}$, and the distance tolerance $\delta =2\Delta x$.

Figure 5 displays movements of two tandem disks with vorticity contours used for flow visualization. From the outset, these two disks fall in tandem. Eventually, the upper disk falls faster than the lower one due to the wake of the lower disk, since the low pressure in lower disk’s wake accelerates the upper one. The hydrodynamic force affects the lower disk as the upper one approaches, and causes the lower disk to drift a little sideways. Ultimately, the upper disk catches up with the lower one, and the two disks collide with each other in what is referred to as the kissing stage. Afterwards, the two disks separate, and then tumble and fall in a side-by-side way. The formerly lower disk is then pushed sideways and upwards by the tumbling motion, and swaps place with the formerly upper disk, which now becomes the lower one. If the rectangular enclosure is long enough, these events would be repeated. The modes of motion delineated above are classified as drafting, kissing and tumbling, which was first observed experimentally by Fortes et al. [22] and simulated in two dimensions by Hu et al. [23], Glowinski et al. [15], and Choi [24]. These well-known phenomena are also predicted by the current DFIB method as shown in Figure 5.

4.2. Sedimentation of Multiple Circular Disks

In order to demonstrate the ability of handling complicated collisions among multiple solid objects by current DFIB method, we also computed the sedimentation of multiple circular disks here. A set of $10\times 10$ identical circular disks are initially placed at the top of a rectangular cavity filled with incompressible Newtonian fluid, released from rest, fall by gravity, and settle at bottom of cavity at last. The parameters and geometry configuration are all the same as the previous example in Section 4.1, except the diameter of disk reduced to be $D=0.15$ cm. The complicated sedimentation process of multiple disks is successfully simulated here, and the result is shown in Figure 6, from which we can observe those drafting, kissing, and tumbling phenomena mentioned above keep happening among disks during sedimentation.

5. Sedimentation of Solid Objects without Rotational Symmetry

5.1. Sedimentation of a Regular Triangle

In this example, we consider the sedimentation of a equilateral triangle with edge length $\sqrt{3}D/2$, where D is the diameter of the circumscribed circle of triangle, in an incompressible viscous fluid. D is 1 cm here. The regular triangle is initially placed at rest with the center-of-mass position at $(15\mathrm{cm},28\mathrm{cm})$ in a $30\mathrm{cm}\times 30\mathrm{cm}$ rectangular domain, with an initial angle of attack (AOA), ${\alpha }_{att}$ = 5${}^{\circ }$. Note that this small AOA would not affect the asymptotic falling pattern of triangle, but is meant to break the symmetry in flow field in order to accelerate the development of its falling pattern. Otherwise, a regular triangle, placed with symmetry initially, would just steadily fall all the way down, since the vertical span is not long enough for disturbances to kick in and break flow symmetry. All the physical parameters are set to be the same as the previous example in Section 4.1. A uniform mesh of size $\Delta x=\Delta y=1/64$ cm, and time step $\Delta t={10}^{-4}$ s are employed in current simulation.

The numerical results are displayed in Figure 7. Figure 7a shows triangle’s trajectory and orientation, and Figure 7b–d shows the snap shots of flow field at certain times during the falling. From that, we can tell the falling pattern is fluttering at the beginning, followed by tumbling at some tipping point. A MATLAB code for current sedimentation simulation of a regular triangle is attached at Appendix A for readers further interested in this problem.

5.2. A Freely Falling Ellipse

It is a common fact that not all objects fall straight down by gravity. For example, a leaf, feather and paper card all fall in a seemingly unpredictable manner. From time to time, a leaf or feather may reverse its falling direction and momentarily rise against the gravity as it flutters or tumbles through air. This rich dynamical behavior has been investigated in many experimental and modeling works, e.g., Reference [1,3,5,25,26]. Numerical simulations regarding a freely falling rectangular or elliptical plate, to mimic a falling leaf or feather, have been studied before in Reference [3,5,25,26]. Taking 2D ellipse as an example, its dynamics can be characterized by the following three dimensionless parameters [3,25]:

• Reynolds number, $R{e}^{*}={\mathit{u}}_t^{*}a/\nu$, with ${\mathit{u}}_t^{*}=\sqrt{2b\mathit{g}({\rho }_s/{\rho }_f-1)}$,
• the dimensionless moment of inertia, $I^{*}=\left({\rho }_{s}b(a^2+b^2)\right)/\left(2{\rho }_f{a}^3\right)$,
• the aspect ratio of the ellipse, $e=a/b$, where a and b are the width and thickness of the ellipse, respectively.

$R{e}^{*}$, $I^{*}$, and e can be recast to I, $Re$, ${\rho }_s/{\rho }_f$, and e, where $I=\left(b(a^2+b^2)\right)/\left(2a^3\right)$, $Re={\mathit{u}}_{t}b/\nu$, with ${\mathit{u}}_t=\sqrt{b\mathit{g}}$. The latter set of parameters is more physically intuitive. Note that the dimensionless moment of inertia without density ratio, I, is actually equivalent to $(e^2+1)/\left(2e^3\right)$, which means I solely depends on e.

Here, we use the current DFIB method to study the falling patterns of a 2D ellipse. The ellipse is initially placed at rest with its center of mass positioned at $(5\mathrm{cm},28\mathrm{cm})$ in a $30\mathrm{cm}\times 30\mathrm{cm}$ rectangular domain, with an initial AOA, ${\alpha }_{att}$ = 45${}^{\circ }$. A uniform mesh of size $\Delta x=\Delta y=1/100$ cm, and time step $\Delta t={10}^{-5}$ s were employed in current simulation. We fixed the density of fluid to be ${\rho }_f=1.0\mathrm{g}/{\mathrm{cm}}^3$ and changed ${\rho }_s$ to vary the density ratio between solid and fluid. We also fixed b = 0.162 cm and varied a to generate various aspect ratios e. Through exhaustive surveys, three combinations of dimensionless parameters as shown in

Table 2. Physical parameters employed for three falling patterns of ellipse.

	${\mathit{Re}}^{*}$	${\mathit{I}}^{*}$	e	${\mathit{\rho }}_{\mathit{s}}/{\mathit{\rho }}_{\mathit{f}}$	$\mathit{\nu }$ [cm${}^2$/s]	$\mathit{Re}$
fluttering	1000	0.08	8	1.2603	0.0118	173.2413
chaotic falling	1000	0.16	8	2.5206	0.0285	71.6779
tumbling	1000	0.25	8	3.9385	0.0396	51.5626

were employed in simulations to demonstrate three typical falling patterns: fluttering, chaotic falling, and tumbling [3,25]. The trajectory of ellipse together with its orientation for these three falling patterns was calculated and is shown in Figure 8. The associated vortex shedding and wake structure of these three falling patterns are shown respectively in Figure 9, Figure 10 and Figure 11.

In Table 2, e is equally set to 8. This is because we found a falling ellipse is more probable to exhibit all three typical falling patterns at a larger e. In addition, the equal setting of e implies an equal setting of I, as well. We further discovered the falling pattern is more sensitive to density ratio, ${\rho }_s/{\rho }_f$, and less sensitive to viscosity $\nu$. By increasing the density ratio gradually, we explored fluttering, chaotic falling, and tumbling, respectively. Generally larger density ratio implies larger driving force and torque, which tends to rotate the ellipse more severely during falling. This can be seen from maximum AOA hardly reaching 90${}^{\circ }$ in fluttering as shown in Figure 8a, but its counterparts in chaotic falling and tumbling can pass 90${}^{\circ }$ during falling as shown in Figure 8b,c.

6. Summary and Conclusions

The current DFIB approach successfully predicted the interaction between falling solid objects and fluid through simulation cases demonstrated above, with some of them validated with the literature. The method is efficient to calculate the resultant force and torque exerted on the solid object and determine its displacement and orientation at each time step. Among all sedimentation cases studied here, rotation becomes important for the sedimentation of a regular triangle and an ellipse due to the breaking of rotational symmetry. The sedimentation of a regular triangle, which is rarely discussed in the literature, shows a falling pattern of fluttering at the beginning, followed by tumbling starting at some tipping point. The sedimentation of an ellipse has the richest dynamics here. With various combinations of physical parameters, our simulation exhibits three typical falling patterns: fluttering, tumbling, and chaotic falling. Animations of sedimentation simulations in current paper can be found at https://www.youtube.com/channel/UC9Qz6kOB-WVFnwijBHZRfVA/.