3D CBC

The goal of this exercise is to familiarize with off-lattice boundary conditions (3D Curved Boundary Conditions) in Palabos

Drag coefficient

Check the available methods of OffLatticeBoundaryCondition3D<Real, Descriptor, Array<Real, 3>>

The definition of drag coefficient is:

$$c_d = \frac{2F_d}{\rho\mu^2A}$$

Fd is the drag force, which is by definition the force component in the direction of the flow velocity

\rho is the mass density of the fluid

\mu is the flow speed of the object relative to the fluid,

A is the reference area.

In the Palabos program, the codes:

1
2
3
4
5
6
7
8
9


        if (iT % parameters.nStep(logT) == 0) {
            cd = boundaryoff->getForceOnObject()[0] /
                      (0.5 * 1.0 * lu.getUlb() * lu.getUlb() * lu.getLref() *
                       lu.getLref() * M_PI);
            std::string fullName =
                global::directories().getLogOutDir() + "Cd.dat";
            plb_ofstream ofile(fullName.c_str(), std::ostream::app);
            ofile << cd << std::endl;
        }

Define the Dynamics

Use the TRT dynamics for low Reynolds and Smagorinsky for the turbulent case

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12


    // 4. Define the dynamics
    if (parameters.getRe() > 2000.0) {
        defineDynamics(*lattice_ptr, lattice_ptr->getBoundingBox(),
                       new SmagorinskyBGKdynamics<Real, Descriptor>(
                           parameters.getOmega(), 0.1));
        pcout << "Using Smagorinsky BGK dynamics." << std::endl;
    } else {
        defineDynamics(
            *lattice_ptr, lattice_ptr->getBoundingBox(),
            new BGKdynamics<Real, Descriptor>(parameters.getOmega()));
        pcout << "Using BGK dynamics." << std::endl;
    }

Change Boundary Conditions

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12


    // 6. Define outer boundary conditions
    OnLatticeBoundaryCondition3D<Real, Descriptor> *onlatt_boundary_condition =
        inject_on_lattice_bc(lattice_ptr, parameters);

    // 7. Define inner-offlattice boundary conditions
    OffLatticeBoundaryCondition3D<Real, Descriptor, Array<Real, 3>>
        *offlatt_boundary_condition =
            inject_off_lattice_bc(lattice_ptr, voxelized_domain);
    lattice_ptr->initialize();
    // return lattice and boundary conditions
    return std::tuple{lattice_ptr, onlatt_boundary_condition,
                      offlatt_boundary_condition};

Check the available methods in palabos/src/offLattice/, available Boundary Conditions for the task:

Guo boundary condition (GZS) , function guoOffLatticeModel3D
Bouzidi boundary condition (BFL), function bouzidiOffLatticeModel3D
Filippova-Hanel boundary condition (FH), function filippovaHaenelOffLatticeModel3D
Mei-Luo-Shyy (MLS), function meiLuoShyyOffLatticeModel3D

Use Sponge Zones

1
2
3


    // 5. Set sponge zones to reduce pressure waves reflections
        createSpongeZones(parameters, Array<plint,6>(5,20,0,0,0,0),
        lattice_ptr);

Add a second sphere

In shape.h file, add generateTwoEllipsoid function:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18


template <typename Real>
TriangleSet<Real> *generateTwoEllipsoid(Array<Real, 3> &center, Real a, Real b,
                                     Real c, Real dx = 1.) {
    Real p = 2.0;
    Real approxSurface =
        4. * M_PI *
        pow(1. / 3. * (pow(a * b, p) + pow(a * c, p) + pow(b * c, p)), 1. / p);
    approxSurface /= (dx * dx);
    TriangleSet<Real> sphere =
        constructSphere<Real>(Array<Real, 3>(0, 0, 0), 1., approxSurface);
    auto& sphere2 = *sphere.clone();
    sphere.scale(a, b, c);
    sphere2.scale(a, b, c);
    sphere.translate(center);
    sphere2.translate(center+Array<Real,3>(5.0*a,.1,.0));
    sphere.append(sphere2);
    return sphere.clone();
}