Direct Forcing Immersed Boundary Method

This section presents the numerical discretization of Fluid Equations. The discretization of the fluid system on a single level is elucidated in Time advancement, followed by an examination of two distinct types of kinematic constraints in Types of kinematic constraints.

Time advancement

Figure 12 Sketch of the two-dimensional semi-staggered grid and variable locations. The blue triangles, red squares, and black circles represent cell-centered variables, node-centered variables, and interface Lagrangian markers.

To solve the partial differential equations of Fluid Equations, the canonical projection [11, 16] is applied to the semi-staggered grid. As shown in Figure 12, the fluid velocity (\(u\) and \(v\)), the Eulerian force \(f\), and particle volume fraction \(\alpha\) are located at the cell center. The pressure \(p\) and level set function \(\phi\) are at the node center. The temporal and spatial discretizations of equations for single-level advancement are considered here. At the time \(t^{n}\), the Eulerian velocity \(\mathbf{u}^{n}\) and pressure \(p^{n-1/2}\) are known. The particle position \(\mathbf{X}^{n}_l\) and velocity \(\mathbf{U}\left(\mathbf{X}^{n}_l\right)\) are also available. The time advancement during the interval \([t^{n}, t^{n+1}]\) proceeds as follows.

Step 1: The intermediate velocity \(\widetilde{\mathbf{u}}^{*,n+1}\) is solved semi-implicitly as

(17)\[\rho_{f}\left(\frac{\widetilde{\mathbf{u}}^{*,n+1}-\mathbf{u}^{n}}{\Delta t}+\nabla \cdot \left(\mathbf{u}\mathbf{u} \right)^{n+\frac{1}{2}}\right)=-\nabla p^{n-\frac{1}{2}}+\frac{1}{2}\left(\nabla \cdot \mu \nabla{\widetilde{\mathbf{u}}^{*,n+1}}+{\nabla} \cdot \mu \nabla {\mathbf{u}}^{n}\right)+\rho_{f}\mathbf{g}\]

where the convective term \({\nabla} \cdot \left(\mathbf{u}\mathbf{u} \right)^{n+\frac{1}{2}}\) is calculated using the second-order Godunov scheme. The Godunov algorithm has been proven to be a very robust method for handling conservation laws, especially for high Reynolds number flows. Many researchers who previously worked on incompressible fluid simulations have applied this method to deal with the convection terms in equations [11, 15, 17]. In our previous work [10, 12], we provided a detailed derivation of the Godunov algorithm. In our recent work [14], we applied this algorithm to simulate two-phase flows at high Reynolds numbers (i.e., \(Re\approx 10^{6}\)) and obtained stable and physically meaningful results. In this step, only the pure fluid system is solved and no particle-related influence is included.

Figure 13 The number of markers change with the \(d/h\) for a single particle scenario, where \(d\) is the particle diameter and \(h\) is the grid spacing.

Step 2: The updated velocity \(\widetilde{\mathbf{u}}^{*, n+1}\) needs to be corrected to satisfy the no-slip boundary condition at the fluid–particle interfaces \(\partial V_{\mathrm{b}}(t)\). This step is divided into four substeps [1, 4] in Algorithm-1. We first interpolate the intermediate Eulerian Velocity obtained from Step 1 to the Lagrangian Velocity of markers. The Lagrangian forces are then calculated based on the desired velocity at the interface and the intermediate velocity. Next, the Eulerian forces are obtained from the spreading of Lagrangian forces by using either the three-point or four-point delta function [1, 4, 10]. As shown in Figure 13, the Eulerian cells, enclosed by red dashed circles, refer to grid areas that are influenced by two blue markers. These two Lagrangian markers also share some intersected areas, which are marked by green arrows. Finally, the Eulerian velocity is corrected by the updated Eulerian Force.

Figure 14 The number of markers change with the \(d/h\) for a single particle scenario, where \(d\) is the particle diameter and \(h\) is the grid spacing.

In Algorithm-1 mentioned above, four points need to be noted. First, the Lagrangian markers only exist on the finest level during the Eulerian-Lagrangian interaction process. This brings the benefits of memory saving since particle-related information does not need to be stored on coarser levels. Second, the Lagrangian markers only distribute on the surface of particles. This is different from the DLM method in others [18, 19], in which the markers also appear inside the particle and there is one marker per Eulerian grid cell. Figure 14 shows how the number of markers changes with the \(d/h\) for a single particle scenario. It is seen that as \(d/h\) increases, the present needs much fewer markers compared with the DLM method. Third, the multi-direct forcing algorithm includes an outer loop with \(m\) ranging from 1 to \(N_{s}\), which controls the degree of coupling between Eulerian and Lagrangian variables. The original method of Uhlmann [20] corresponds to the case of \(N_{s} = 0\). Increasing \(N_{s}\) can enhance their coupling but will also increase the computational load. Based on the experience in the previous work [1, 4] and tests presented in this paper, it is sufficient to set \(N_{s}\) to 2-3 for all cases in Validation Cases. Finally, if a system has multiple particles, each particle goes into 1 to \(N_{s}\) loop sequentially. The corrected Eulerian forces we employ take into account the effects of all particles. This consideration also applies to the calculation of the particle volume fraction (PVF) field in Particle Volume Fraction (PVF).

Algorithm-1 : Multidirect forcing method for fluid-particle interaction

1. \(\mathbf{u}^{(0)} = \widetilde{\mathbf{u}}^{*, n+1}\)

2. For \(m=1\) to \(N_{s}\):

   1. Interpolate Lagrangian Velocity, \(\mathbf{U}^{m-1}\left(\mathbf{X}^{n}_l\right)=\sum_{i=1}^{N_x} \sum_{j=1}^{N_y} \sum_{k=1}^{N_z} \mathbf{u}^{(m-1)}\left(\mathbf{x}_{i, j, k}\right) \delta_h\left(\mathbf{x}_{i, j, k}-\mathbf{X}^{n}_l\right) h^3\)

   2. Calculate Lagrangian Force, \(\mathbf{F}^{m}=\mathbf{F}^{m-1}+\left(\mathbf{U}^d\left(\mathbf{X}^{n}_l\right)-\mathbf{U}^{m-1}\left(\mathbf{X}^{n}_l\right)\right) / \Delta t\)

   3. Spreading Lagrangian Force onto Eulerian Force, \(\mathbf{f}^{m}\left(\mathbf{x}_{i, j, k}\right)=\sum_{l=1}^{N_l} \mathbf{F}^{m}\left(\mathbf{X}_l\right) \delta_h\left(\mathbf{x}_{i, j, k}-\mathbf{X}_l\right) \Delta V_l\)

   4. Correct Eulerian Velocity, \(\mathbf{u}^{(m)}=\mathbf{u}^{(0)}+\Delta t \mathbf{f}^{m}\left(\mathbf{x}_{i, j, k}\right)\)

3. \({\mathbf{u}}^{*, n+1} = \mathbf{u}^{(m)}\)

Step 3: With the updated intermediate velocity \({\mathbf{u}}^{*, n+1}\) in Step 2, a level projection operator is applied to obtain the updated pressure \(p^{n+1/2}\) and velocity \({\mathbf{u}}^{n+1}\) fields. An auxiliary variable \(\boldsymbol{V}\) is first calculated by

(18)\[\boldsymbol{V} = \frac{{\mathbf{u}}^{*,n+1}}{\Delta t} + \frac{1}{\rho_{f}} {\nabla} p^{n-\frac{1}{2}}.\]

Then, \(\boldsymbol{V}\) is projected onto the divergence-free velocity field to obtain the updated pressure \(p^{n+1/2}\) via

(19)\[L^{cc,l}_{\rho_{f}} p^{n+1/2} = {\nabla} \cdot \boldsymbol{V},\]

where \(L^{cc}_{\rho_{f}}p^{n+1/2}\) is the density-weighted Laplacian operator to \({\nabla} \cdot (1/\rho_{f}{\nabla} p^{n+1/2})\) [11, 12]. Finally, the divergence-free velocity \({\mathbf{u}}^{n+1}\) on level \(l\) is obtained as

(20)\[{\mathbf{u}}^{n+1} = \Delta t \left(\boldsymbol{V} - \frac{1}{\rho_{f}} {\nabla} p^{n+1/2}\right).\]

The projection is stable and appears to be well-behaved in various numerical tests [21, 22] and practical applications [13, 15].

Step 4:After completing Step 3, we obtain the divergence-free fluid velocity at \(t^{n+1}\). The particle-related information also needs to be updated from \(t^{n}\) to \(t^{n+1}\). Depending on different kinematic constraints, the particle motion is categorized into prescribed motion and free motion. The specific updates are detailed in Types of kinematic constraints.

Types of kinematic constraints

The approximation to the constrained Lagrangian velocity \((\mathbf{U}_{\mathrm{b}})_{q, m}^{n+1}\) and position \(\mathbf{X}_{q, m}^{n+1}\) depends on the type of kinematic constraint in the FSI. In this work, we consider two types of kinematic constraints: prescribed motion of the particle and free motion of the particle.

Prescribed motion

If the motion of the particle is prescribed, then its velocity and position are known a priori and not influenced by the surrounding fluid. Thus, the centroid position \(\mathbf{X}_{r}^{n}\), centroid velocity \(\mathbf{U}_{r}^{n}\) at \(t^{n}\), centroid velocity \(\mathbf{U}_{r}^{n+1}\) at \(t^{n+1}\), and angular velocity \(\mathbf{W}_{r}^{n}\) of the body are given. The desired velocity \(\mathbf{U}^d\left(\mathbf{X}^{n}_l\right)\) of the markers in Algorithm-1 is calculated as

(21)\[\mathbf{U}^d\left(\mathbf{X}^{n}_l\right)=\mathbf{U}_{r}^{n}+\mathbf{W}_{r}^{n} \times \mathbf{R}_l^n,\]

where \(\mathbf{R}_l^n = \left(\mathbf{X}^{n}_{l} - \mathbf{X}_{r}^{n}\right)\). The new position of the centroid of the particle \(\mathbf{X}_{r}^{n+1}\) is updated using the midpoint scheme as

(22)\[\mathbf{X}_{r}^{n+1}=\mathbf{X}_{r}^{n}+\frac{\Delta t}{2}(\mathbf{U}_{r}^{n+1}+\mathbf{U}_{r}^{n}).\]

Free motion

In contrast to the prescribed kinematics case, the motion of a freely moving particle is influenced by the surrounding fluid. To account for this two-way interaction, the following governing equations of the particle systems are solved [1, 4].

(23)\[\begin{split}\begin{equation} \label{particle1} \rho_p V_p \frac{d \mathbf{U}_r}{d t} \approx-\rho_f \sum_{l=1}^{N_L} \mathbf{F}_l^{n+1/2} \Delta V_l+\rho_f \frac{d}{d t}\left(\int_{V_p} \mathbf{u} d V\right) \\ +\left(\rho_p-\rho_f\right) V_p \mathbf{g}+\mathbf{F}_c^{n+1 / 2} \end{equation}\end{split}\]

(24)\[\begin{split}\begin{equation} \label{particle2} I_p \frac{d \mathbf{W}_r}{d t} \approx-\rho_f \sum_{l=1}^{N_L} \mathbf{R}_l^n \times \mathbf{F}_l^{n+1 / 2} \Delta V_l \\ +\rho_f \frac{d}{d t}\left(\int_{V_p} \mathbf{r} \times \mathbf{u} d V\right)+\mathbf{T}_c^{n+1 / 2} \end{equation}\end{split}\]

On the right-hand side of (23)-(24), the term \(\mathbf{F}_l^{n+1/2}\) refers to the Lagrangian Force, coming from the final value of \(\mathbf{F}^{m}\) in Algorithm-1. The time derivatives of momentum integration \(\rho_f \frac{d}{d t}\left(\int_{V_p} \mathbf{u} d V\right)\) and angular momentum integration \(\rho_f \frac{d}{d t}\left(\int_{V_p} \mathbf{r} \times \mathbf{u} d V\right)\) within the particle are also included. These two integrated terms account for flow unsteadiness by using the PVF field (Particle Volume Fraction (PVF)). The term \(\left(\rho_p-\rho_f\right) V_p \mathbf{g}\) considers the buoyancy effects. The terms \(\mathbf{F}_c^{n+1 / 2}\) and \(\mathbf{T}_c^{n+1 / 2}\) refer to the induced force and torque generated by the particle collision, respectively. If there is only one single particle in the system, both \(\mathbf{F}_c^{n+1 / 2}\) and \(\mathbf{T}_c^{n+1 / 2}\) are set to be zero. On the left-hand side of (23)-(24), we use the second-order mid-point scheme to integrate particle motions [3]. After updating the particle centroid velocity \(\mathbf{U}_{r}^{n+1}\) and angular velocity \(\mathbf{W}_{r}^{n+1}\) at \(t^{n+1}\), we go back to (22) to update new position of the particle centroid \(\mathbf{X}_{r}^{n+1}\).

The time advancement scheme in this work is not fully implicit [3], yet it can deal with the free motion applies to particles either with a large density ratio (i.e., \(\frac{\rho_p}{\rho_f} \geq 10\)) or a small density ratio (i.e., \(\frac{\rho_p}{\rho_f} \approx 1.5 - 2\)) [1, 4]. We found it is robust and fast enough to handle all the testing cases in Validation Cases. Before ending this Section, we also emphasize that our method is similar to the “weak coupling” method used in the sharp-interfaced immersed boundary method, which requires only one solution for fluid and solid solver during each time step and no iterations are needed between these two solvers [23, 24]. It is easier to extend the current portable solver to the “strong coupling” method, which then re-projects the flow part, re-updates the solid particle, and performs a convergence checking between the fluid solver and the solid solver during each sub-iteration [10, 25].