The goal of this exercise is to familiarize with off-lattice boundary conditions (2D Curved Boundary Conditions) in Palabos
CMAKE file
At first, you need compile the CMAKE file, and produce a EXE file. For this file, you need to know where to change the project name and define the executable_name, when necessary.
project(cavity2d)
enable_language(CXX)
set(EXECUTABLE_NAME "cavity2d_incomplete_with_bug")
set (CMAKE_RUNTIME_OUTPUT_DIRECTORY "../")
Most importantly, you need set the PALABOS_ROOT. Otherwise you will face a compile error after the command cmake ..
# # # If the palabos root folder is in your path, you don't need to modify anything.
# # # In the other case you have 2 options: 1. add PALABOS_ROOT=path/to/palabos to the path of your system
# # # 2. you comment the next line and uncomment the one after, explicitly giving the path/to/palabos
# # # (NB you can use ${CMAKE_CURRENT_SOURCE_DIR} to give paths relative to the CMakeLists.txt location)
#file(TO_CMAKE_PATH $ENV{PALABOS_ROOT} PALABOS_ROOT)
set(PALABOS_ROOT "../../.")
cbc_example_2d.cpp
// Define simulation parameters: we use two ad-hoc units helper check
// simulation_parameters.h for the documentation.
const Real re = 250; // NB: Obstacle reynolds!
sp::Numerics<Real, Int> lu;
sp::NonDimensional<Real> dimless;
dimless.initReLxLyLz(re, 24, 4, 4).printParameters();
IncomprFlowParam<Real> parameters =
lu.initLrefluNodim(50, &dimless, 0.01, false) //25 0.02
.printParameters()
.getIncomprFlowParam(); // be careful with this casting, check the documentation
Definition for the initLrefluNodim function
// LOCATION: Line #209 of File (simulationParameters.cpp)
/// Generates the parameters in lattice units from an object with the non
/// dimensional parameters and the lattice resolution, i.e. l_ref_lu_. l_ref_lu_
/// is Real, but the generated box dimension are rounded to int! \param
/// l_ref_lu_ resolution (reference length in lattice units) \param dimless_
/// NonDimensional<Real> object that has already been initialized
/// @param u_lb_ if given > 0, it overrides the Ma number in the dimless_
/// @param tau_ if given > 0, it overrides the Ma number in the dimless_
/// @return *this
template <typename Real, typename Int>
Numerics<Real, Int>& Numerics<Real, Int>::initLrefluNodim(
Int l_ref_lu_, NonDimensional<Real>* dimless_, Real u_lb_, Real tau_) {
assert(dimless_->initialized());
l_ref_lu = l_ref_lu_;
if (u_lb_ > getEpsilon(1)) {
setUlb(u_lb_);
setTau(3. * (u_lb_ * l_ref_lu_ / dimless_->getRe()) + 0.5);
if (tau_ > 0.499)
pcout << "initLrefluNodim() Waring: you cannot set both u_lb_ "
"and tau_. Giving priority to u_lb_...."
<< std::endl;
} else if (tau_ > 0.499) {
setUlb(dimless_->getRe() * getNulu(tau_) / l_ref_lu_);
setTau(tau_);
} else {
setUlb(dimless_->getMa() * getCs());
setTau(getCs2Inv() * (u_lb * l_ref_lu_ / dimless_->getRe()) + 0.5);
}
lx_domain = util::roundToInt(dimless_->getLx() * l_ref_lu_);
ly_domain = util::roundToInt(dimless_->getLy() * l_ref_lu_);
lz_domain = util::roundToInt(dimless_->getLz() * l_ref_lu_);
is_initialized = true;
delete dimless;
dimless = dimless_;
return *this;
}
Output of plbLog.dat
Poiseuille flow
Velocity in lattice units: u=0.02
Reynolds number: Re=250
Lattice resolution: N=25
Relaxation frequency: omega=1.97628
Extent of the system: lx=24
Extent of the system: ly=4
Extent of the system: lz=4
Grid spacing deltaX: dx=0.04
Time step deltaT: dt=0.0008
Poiseuille Flow
Poiseuille flow is pressure-induced flow (Channel Flow) in a long duct, usually a pipe.
The current Poiseuille setup, in the definition of setup_boundary_conditions_2d() (file: setup.h):
template <typename Real, template <typename U> class Descriptor>
void setup_boundary_conditions_2d(
MultiBlockLattice2D<Real, Descriptor>& lattice,
IncomprFlowParam<Real> const& parameters,
OnLatticeBoundaryCondition2D<Real, Descriptor>& boundaryCondition) {
const plint nx = parameters.getNx();
const plint ny = parameters.getNy();
Box2D outlet(nx - 1, nx - 1, 1, ny - 2);
// Create Velocity boundary conditions everywhere
boundaryCondition.setVelocityConditionOnBlockBoundaries(
lattice, Box2D(0, 0, 1, ny - 2));
boundaryCondition.setVelocityConditionOnBlockBoundaries(
lattice, Box2D(0, nx - 1, 0, 0));
boundaryCondition.setVelocityConditionOnBlockBoundaries(
lattice, Box2D(0, nx - 1, ny - 1, ny - 1));
// .. except on right boundary, where we prefer an outflow condition
// (zero velocity-gradient).
boundaryCondition.setVelocityConditionOnBlockBoundaries(
lattice, Box2D(nx - 1, nx - 1, 1, ny - 2), boundary::outflow);
setBoundaryVelocity(lattice, lattice.getBoundingBox(),
PoiseuilleVelocity<Real>(parameters));
setBoundaryDensity(lattice, outlet, ConstantDensity<Real>(1.));
initializeAtEquilibrium(
lattice, lattice.getBoundingBox(),
PoiseuilleVelocityAndDensity<Real, Descriptor>(parameters));
plint cx = nx / 4;
plint cy = ny / 2 + 2; // cy is slightly offset to avoid full symmetry,
// and to get a Von Karman Vortex street.
plint radius = parameters.getResolution() / 2.0;
defineDynamics(lattice, lattice.getBoundingBox(),
new CylinderShapeDomain2D<Real, Descriptor>(cx, cy, radius),
new plb::BounceBack<Real, Descriptor>);
lattice.initialize();
}
The GIF animation for Poiseuille flow looks like:
Boundary Conditions
No-slip boundary condition: at the interface between a moving fluid and a stationary wall, both the normal and tangential components of the fluid velocity field are equal to zero.
Free-slip boundary condition: at the interface between a moving fluid and a stationary wall, the normal component of the fluid velocity field is equal to zero, but the tangential component is unrestricted.
-
Slip: the normal component of the velocity is zero, i.e., no flux across the boundary
-
Free: the tangential force is zero.
To change to Free-slip boundary condition, three aspects to modify:
- In function
setVelocityConditionOnBlockBoundariesuseboundary::freeslip(replaceboundary::outflow) - Convert the initialization and the boundary velocities to a uniform profile. You need to change
PoiseuilleVelocityandPoiseuilleVelocityAndDensitytoConstantVelocityandConstantVelocityAndDensity, respectively. - Finally put a free slip (specular reflection) on the cylinder, instead of
plb::BounceBack, changing the dynamics (useplb::SpecularReflection).
The modified function was displayed as follows:
// Free-slip boundary condition:
template <typename Real, template <typename U> class Descriptor>
void setup_boundary_conditions_2d(
MultiBlockLattice2D<Real, Descriptor>& lattice,
IncomprFlowParam<Real> const& parameters,
OnLatticeBoundaryCondition2D<Real, Descriptor>& boundaryCondition) {
const plint nx = parameters.getNx();
const plint ny = parameters.getNy();
Box2D outlet(nx - 1, nx - 1, 1, ny - 2);
// Create Velocity boundary conditions everywhere
boundaryCondition.setVelocityConditionOnBlockBoundaries(
lattice, Box2D(0, 0, 1, ny - 2));
boundaryCondition.setVelocityConditionOnBlockBoundaries(
lattice, Box2D(0, nx - 1, 0, 0));
boundaryCondition.setVelocityConditionOnBlockBoundaries(
lattice, Box2D(0, nx - 1, ny - 1, ny - 1));
// .. except on right boundary, where we prefer an freeslip condition
// (zero velocity-gradient).
boundaryCondition.setVelocityConditionOnBlockBoundaries(
lattice, Box2D(nx - 1, nx - 1, 1, ny - 2), boundary::freeslip); //Change#1
setBoundaryVelocity(lattice, lattice.getBoundingBox(),
ConstantVelocity<Real>(parameters)); //Change#2.1
setBoundaryDensity(lattice, outlet, ConstantDensity<Real>(1.));
initializeAtEquilibrium(
lattice, lattice.getBoundingBox(),
ConstantVelocityAndDensity<Real, Descriptor>(parameters));//Change#2.2
plint cx = nx / 4;
plint cy = ny / 2 + 2; // cy is slightly offset to avoid full symmetry,
// and to get a Von Karman Vortex street.
plint radius = parameters.getResolution() / 2.0;
defineDynamics(lattice, lattice.getBoundingBox(),
new CylinderShapeDomain2D<Real, Descriptor>(cx, cy, radius),
new plb::SpecularReflection<Real, Descriptor>);//Change#3
lattice.initialize();
}
The GIF animation for Channel flow with free-slip boundary looks like:
