This module provides algorithms that allow coupled simulations of LBM for the fluid part and rigid particles from the PE module. This coupling can be achieved in various ways, depending on whether the particles should be fully resolved by the fluid flow around it (as in a direct numerical simulation) or if modelling assumptions are introduced that allow for unresolved particles.
The following coupling techniques are provided by the framework:
This is a coupling method for geometrically fully resolved simulations, i.e. the numerical resolution is significantly finer than the size of the particle. It is based on the work of [19] and [1]. It explicitly maps the particles into the fluid domain by flagging the containing cells as "obstacle". The coupling is then established by a velocity boundary condition for LBM along the particles' surface and by calculating the hydrodynamic force acting on the particles. This is the main coupling technique that we are using in waLBerla.
See [25] for more infos and benchmark tests. See [28] for a large-scale application of this coupling.
This is a second option for fully resolved coupled simulations and is based on the work of [23]. It uses a solid volume fraction field that stores the occupancy of the fluid cells by the solid phase. This scalar value is then used in a special LBM collision operator that models the interaction by a weighting between the fluid collision and the solid collision. This collision step results in the interaction forces that are then set onto the particles.
See [25] for more info and benchmark tests.
This is a completely different class of methods that uses particles that are smaller than an LBM cell. This necessitates the introduction of models for the fluid-particle interaction forces, where the drag and pressure forces are usually the most important ones. In short, fluid quantities like velocity are interpolated to the particle positions and used there to evaluate the empirical formulas for the interaction forces. This forces are then applied to the particles and the corresponding reaction-force is distributed back to the surrounding fluid cells. For denser systems, the equations of fluid motion have to be altered to incorporate displacement effects by the solid phase, yielding the volume-averaged Navier-Stokes equations. For higher Reynolds number flows, also turbulence models might become necessary since the resolution for the fluid motion is usually very coarse. The results are highly dependent on the included models and extensive pre-studies have to be undertaken to incorporate all important effects. The current implementation provides a variety of the mentioned functionalities.
See [26] for more information.
Since it is a coupled simulation, you first of all need both, the LBM parts and the PE parts. So also the restrictions coming from these modules apply (like number of blocks in periodic setups with particles, see Tutorial - Confined Gas: Setting up a simple pe simulation).
We only focus on the additional coupling part here. The functionalities for that are found in src/pe_coupling, and the tests in tests/pe_coupling. A good starting point is tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp as it contains a lot of the mentioned functionality.
Every PE particle can be mapped into the domain as long as its type has the pe::RigidBody::containsPoint() member function implemented. With this mechanism, also a object from the mesh module can be mapped when seen as a mesh::pe::ConvexPolyhedron. Generally, the mapping sets the cells with the cell center contained inside the particle to a given flag that signals that it is not fluid.
Before the simulation starts, the particle have to be mapped into the domain.
Two mechanisms / functionalities are provided: mapBodies() and mapMovingBodies().
The mapBodies() functionality is used to set a constant boundary condition (like lbm::NoSlip), i.e. one does not change in time. With this, e.g. bounding planes can be mapped to set up the outer boundaries of the fluid domain.
The mapMovingBodies() functionality is similar but needs an additional data structure, called the body field, which stores a pointer (pe::BodyID) inside each flagged cell to the containing particle. In this way, e.g. the particle corresponding to a specific boundary cell can be found so that its local surface velocity can be obtained and hydrodynamic forces are applied to the correct particle This is the functionality that is needed for moving particles in your coupled simulation.
Make sure that you map the particle with ONLY ONE of them to avoid overwriting of the already set flags. The mapping functions can be given selector functors that allow you to specify rules on which particles should be mapped. Several can be found in src/pe_coupling/utility/BodySelectorFunctions.h, like pe_coupling::selectAllBodies().
For the moving particles, the pe_coupling::BodyMapping class is available and to be used as a sweep inside the time loop to update the mapping of these particles. It additionally requires a second flag, the formerObstacle flag, whose usage will become clear in the next section. In essence, this class transforms fluid to obstacle cells if required and turns newly uncovered obstacle cells to formerObstacle cells.
As the particles move across the fluid mesh, cells will become covered and uncovered by the particle. Covering simply turns the former fluid cells to obstacle cells (see Update mapping during simulation). Extra work is required for cells that turn from obstacle to fluid before the simulation can continue as the information about the fluid is missing. This is where we need to use the pe_coupling::PDFReconstruction that restores the missing PDF information based on a templated reconstructor.
Available reconstructors are:
The two latter ones require an extrapolation direction. This can be provided in two ways:
After reconstructing the PDFs in a cell, the flag is turned from formerObstacle to fluid and the pointer in the body field is invalidated. Note that this procedure is split into two separate loops (1. reconstruct PDFs in all formerObstacle cells, 2. turn all formerObstacle flags to fluid). The reason is that the reconstructors look for valid fluid cells in the vicinity of the uncovered cell, which is done based on the fluid flag, and we want to avoid that recently reconstructed values are regarded as valid data.
As for the reconstruction part particle information is still required (e.g. the particle's surface velocity), we use the overlap of \(1.5\) cell length to keep the particle available (see arguments in Initial mapping). This originates from \(0.5\) mentioned in Initial mapping, and a totally maximal admissible position change of \(1\) cell per time step, so \(1.5\) in total.
This is the part where the main coupling is happening. The fluid is affected by the moving particle via velocity boundary conditions that use the local surface velocity of the particle. As a result, the hydrodynamic force and torque acting on the particle is obtained. Both parts are computed in the boundary handling routine.
Different variants of different order are available:
The two latter ones require the exact position of the particle surface. See [25] for more details and comparison tests of the boundary conditions. This is calculated as a ray-particle intersection which can be done analytically for a sphere (intersectionRatioSpherePe()), ellipsoid (intersectionRatioEllipsoidPe()) or a plane (intersectionRatioPlanePe()) and is done by a expensive bisection line search in other cases (intersectionRatioBisection()). If you want to add your specialization, do this by extending the cases in the intersectionRatio() variant that takes a PE body.
Once the forces and torques on the particles have been calculated, the physics engine is used to compute the inter-particle collision forces and to update the particle's position and orientation, as well as the linear and angular velocity. This is conveniently achieved by using the pe_coupling::TimeStep class that can be inserted into the time loop. It also features a sub-stepping capability where several PE time steps with smaller time step sizes are carried out which is desirable in closely packed scenarios to avoid large overlaps of the particles, see e.g. [8]. During this sub-stepping, the hydrodynamic forces and torques are kept constant.
The typical coupled simulation looks like the following:
Initially: set up data structures, map particles, set boundary conditions
Time loop:
This is a collection of lessons learned and tips that should be kept in mind when setting up a coupled simulation:
Since the fluid and particle part of the simulation generate different load characteristics, load balancing is non-trivial and different approaches are possible. In waLBerla, we generally want to maintain the same domain partitioning of the fluid and the particle simulation.
See [27] for one possibility. The corresponding application codes can be found in apps/benchmarks/AdaptiveMeshRefinementFluidParticleCoupling/AMRSedimentSettling.cpp to see how the load balancing routines are set up.
Functions | |
| pe::Vec3 | walberla::pe_coupling::LubricationCorrection::compLubricationSphrSphr (real_t gap, const pe::SphereID sphereI, const pe::SphereID sphereJ) const |
| Computes lubrication correction force between spheres. More... | |
| pe::Vec3 | walberla::pe_coupling::LubricationCorrection::compLubricationSphrPlane (real_t gap, const pe::SphereID sphereI, const pe::ConstPlaneID planeJ) const |
| Computes lubrication correction force between sphere and wall. More... | |
|
private |
Computes lubrication correction force between sphere and wall.
Lubrication correction according to Ladd and Verberg, 2001 ("Lattice-Boltzmann Simulations of Particle-Fluid Suspensions" )
Note: Verified quantitatively by computation in spreadsheet and qualitatively by considering direction of force for example setup.
|
private |
Computes lubrication correction force between spheres.
Lubrication correction according to Ladd and Verberg, 2001 ("Lattice-Boltzmann Simulations of Particle-Fluid Suspensions")
Note: Verified quantitatively by computation in spreadsheet and qualitatively by considering direction of force for example setup.