Solid boundaries and the fluid‐particle coupling - gnomeCreative/HYBIRD GitHub Wiki

Solid boundaries and the fluid-particle coupling

The particles constituting the granular phase are solved using the DEM (see DEm documentation). However, the computation of the hydrodynamic interaction forces is performed within the LBM framework. The traditional way to handle the fluid-particle coupling is by considering the surface of the particles as a moving no-slip boundary condition. This method is based on the bounce back mechanism. However, this method is unstable if the particles are moving quickly inside the lattice, a fact that motivated the development of alternative methods. Among the many available possibilities, the Immersed Boundary Method with Direct Forcing seems to be the most suitable for debris flow simulations, because of its efficiency and stability.

The discretization of particles in the lattice, regardless of the coupling method, follows the same logic used for the resolution of the free surface. Nodes are divided into two categories: solid and fluid. Fluid nodes can be either of the three categories already described (i.e. liquid, interface, or gas), while solid nodes are all the nodes that are located inside the volume of a grain. The main difference between the "dry" and the "wet" coupling scheme is in the treatment of the solid nodes. The dry scheme works with fluid only in the exterior of the particles, and therefore the solid nodes are inactive and do not contain fluid. In the wet coupling the fluid is also present inside the solid nodes.

The bounce-back condition and the dry coupling scheme

The dry coupling is based on the no-slip boundary condition, which is commonly implemented in LBM through the bounce-back rule. This states that whenever a population is streaming towards a wall, this population is reflected and bounced back in the opposite direction. In LBM notation

$f_{i'} (x,t+1) = f_{i} (x,t)$,

where $x$ is the position of the fluid node close to the wall, $i$ the streaming direction that points at the wall, and $i'$ the opposite direction. If the wall is moving, the reflection has to take into account the momentum transfer between population and wall. The correction is implemented by imposing

$f_{i'} (x,t+1) = f_{i} (x,t) - 6 w_i \rho_\textrm{f} u_\textrm{w} \cdot c_i$,

where $u_\textrm{w}$ is the velocity of the wall at the bounce-back location. To achieve second-order accuracy, the bounce-back location should be located exactly halfway between the solid node and fluid node. In a debris-flow simulation, this rule is applied whenever the fluid is at contact with a wall, e.g. the channel bed. The same rule can also be applied to grains (see Ladd, 1994), with the additional complexity of calculating the velocity of the grain surface at every time step. The grains never undergo large deformations, and can therefore be safely assumed to be rigid bodies for what concerns the determination of $u_\textrm{w}$. If $u_\textrm{p}$ and $\omega_\textrm{p}$ are the translational and rotational velocities of the particle, respectively, the velocity of the particle at the bounce back-location can be determined as

$u_\textrm{w}=u_\textrm{p} + r_\textrm{BB} \times \omega_\textrm{p}$,

where $r_\textrm{BB}$ is the vector connecting the centre of the grain to the bounce-back location.

The momentum exchange described by the bounce-back equation results in a force applied to the wall, in direction $i$. The resultant of all forces given by momentum exchange on a wall can be written as

$F_\textrm{hydro}=\sum \left( 2 f_{i} (x,t) - 6 w_i \rho_f u_w \cdot c_i \right) c_i$,

and the resultant moment

$M_\textrm{hydro}=\sum \left( 2 f_{i} (x,t) - 6 w_i \rho_f u_w \cdot c_i \right) c_i \times r_\textrm{BB}$.

This implementation was initially suggested by Aidun & Lu (1995), who published a scheme that improved the previous method by Ladd (1994) removing the requirement that particles are filled with fluid. This method converges at the limit of $\Delta x \rightarrow 0$, since the zig-zag description of the particle surface becomes closer to the exact shape for finer lattice resolutions. This is not an important problem in a debris-flow simulation because the force on the grains does not need to be computed exactly. In fact, the spherical shape of the particles is anyway only a crude approximation of the actual shape of the grains, which is irregular.

The motion of the particles in the lattice requires a continuous update of the state of the nodes. A moving particle creates new solid nodes at its front and new fluid nodes at its rear. This is problematic for two reasons. The newly created fluid nodes need to be initialized, and their velocity, density, and liquid fraction are in general unknown. The easiest way to deal with this is by computing the new variables using their averaged value over the surrounding liquid and interface nodes. Another problem is the appearance and disappearance of mass whenever a fluid node is created or destroyed. While the total mass is likely to only slightly oscillate around its initial value, the transfer of mass between wake and front can strongly affects the dynamics of the flow. These problems motivated the search for more precise coupling method, and many alternatives have been developed in the last two decades. Particularly successful are the application that use the Immersed Boundary Method (e.g. Leonardi et al., 2012) or the fictitious domain technique (e.g., Glowinski et al., 2001). The drawback of these methods is that, the higher their precision, the more complex and computationally expensive their implementation becomes. For this reason, the next section presents a technique, the simplified direct forcing, which is a compromise between precision and speed. Note that even though the dry scheme described here is unsuitable for the coupling with the particle, it is still used for any other stationary object, like the bottom of the flow and in general any other stationary wall.

The wet coupling scheme: the direct forcing

The Immersed Boundary Method (IBM) and similar schemes have enjoyed lately increasing popularity because they allow a more precise resolution of the particle shapes without a small lattice spacing. The method proposed by Feng & Michaelides (2005) uses an additional discretization of the particles into a set of segments, each characterized by its center of mass $x_\textrm{seg}$, velocity $u_\textrm{seg}$ and volume $V_\textrm{seg}$. For every segment, the velocity of the fluid calculated at the position $x_\textrm{seg}$ is computed with an interpolating function $\mathcal{F}$ running over every lattice node $n$ in the neighborhood

$u_\textrm{f}(x_\textrm{seg},t)=\sum_n \mathcal{F}(u_ {\textrm{f},n}(x,t))$.

This velocity is used to compute a force $F_\textrm{hydro,seg}$ for every segment as

$F_\textrm{hydro,seg}(x_\textrm{seg},t)=\rho_\textrm{f}(x,t)\left[ u_\textrm{f}(x_\textrm{seg},t) - u_\textrm{seg}(x_\textrm{seg},t)\right] V_\textrm{seg}$.

This force is transmitted to the DEM for solution of the particle dynamics. Its counterpart for the LBM is a volumetric force $e(x,t)$, whose determination is carried on using the same interpolation function $\mathcal{F}$, this time running over every $m$ segment $x_{\textrm{seg},m}$ in the proximity of a lattice node locations $x$:

$e(x,t)=-\sum_m \mathcal{F}(F_\textrm{hydro,seg}(x_{\textrm{seg},m},t))$.

This force $e(x,t)$ is then applied to the fluid together with the other body forces. Note that this approach requires the interior of the particle to be filled with fluid, whose velocity quickly relaxes to the velocity of the particle. To fasten the relaxation, the interior of the particle can be assigned a higher viscosity (e.g. $\mu_\textrm{max}$).

The direct forcing described so far has the great advantage of allowing the resolution of particles of arbitrary shapes, and has been used with success for the simulation of fiber-reinforced fresh concrete by Svec et al, (2012). For the simulations of debris flows, where particles are mostly rounded, the additional complexity given by the generation of the segments might be unjustified. For this reason, a simplified version is used, where the segments correspond to the cubic cells of the lattice, and their center of mass to the lattice nodes $x_\textrm{cell}=x$. This results in a simpler and faster algorithm, at the price of a less precise resolution of the particle shape, and a less smooth evolution in time of the force. No interpolation is required, the volume of a cubic cell is unitary, and the force can be computed directly for the particle as

$F_\textrm{hydro,cell}(x,t)=\rho_\textrm{f}(x,t)\left[ u_\textrm{f}(x,t) - u_\textrm{p}(x,t)\right]$,

and for the fluid as

$e(x,t)=-F_\textrm{hydro,cell}(x,t)$.

This approach works well if the particle is large enough to cover a significant number of lattice nodes. The resultant force on the particle is then computed by simple summation over each of the $l$ covered nodes

$F_\textrm{hydro}=\sum_l F_ {\textrm{hydro},l}$,

together with the resultant moment

$M_\textrm{hydro}=\sum_l F_{\textrm{hydro},l}\times r_{\textrm{DF},l}$,

where $r_{\textrm{DF},l}$ is the vector connecting the center of mass of the particle and the node location $x_l$.