2D CBC

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.

1
2
3
4
5


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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10


   // 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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36


// 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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11


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):

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38


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:

Poiseuille Flow

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 setVelocityConditionOnBlockBoundaries use boundary::freeslip (replace boundary::outflow)
Convert the initialization and the boundary velocities to a uniform profile. You need to change PoiseuilleVelocity and PoiseuilleVelocityAndDensity to ConstantVelocity and ConstantVelocityAndDensity, respectively.
Finally put a free slip (specular reflection) on the cylinder, instead of plb::BounceBack, changing the dynamics (use plb::SpecularReflection).

The modified function was displayed as follows:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39


// 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:

Free-slip Boundary

2D CBC

CMAKE file

cbc_example_2d.cpp

Poiseuille Flow

Boundary Conditions

相关文章：