I am a new user of palabos.
I try to solve sudden Expansion flow problem.
In order to solve it I took the tutorial code of Poiseuille flow and changed the geometry.
However when running the simulation it seems that the walls are not solid and I still recieve a tube tunnel in the result.
The code is attached and I will be greatfull for any asistance.
Thank you.
[code=“cpp”]
#include “palabos2D.h”
#include “palabos2D.hh”
#include
#include
#include
/* Code 1.5 in the Palabos tutorial
*/
using namespace plb;
using namespace std;
typedef double T;
#define DESCRIPTOR plb::descriptors::D2Q9Descriptor
const int nx = 1200;
const int ny = 400;
/// Velocity on the parabolic Poiseuille profile
T poiseuilleVelocity(plint iY, IncomprFlowParam const& parameters) {
T y = (T)iY / parameters.getResolution();
return 4.parameters.getLatticeU() * (y-yy);
}
/// A functional, used to initialize the velocity for the boundary conditions
template
class PoiseuilleVelocity {
public:
PoiseuilleVelocity(IncomprFlowParam parameters_)
: parameters(parameters_)
{ }
/// This version of the operator returns the velocity only,
/// to instantiate the boundary condition.
void operator()(plint iX, plint iY, Array<T,2>& u) const {
u[0] = poiseuilleVelocity(iY, parameters);
u[1] = T();
}
/// This version of the operator returns also a constant value for
/// the density, to create the initial condition.
void operator()(plint iX, plint iY, T& rho, Array<T,2>& u) const {
u[0] = poiseuilleVelocity(iY, parameters);
u[1] = T();
rho = (T)1;
}
private:
IncomprFlowParam parameters;
};
void channelSetup (
MultiBlockLattice2D<T,DESCRIPTOR>& lattice,
IncomprFlowParam const& parameters,
OnLatticeBoundaryCondition2D<T,DESCRIPTOR>& boundaryCondition )
{
Box2D topLow(0, nx/2-1, ny*5/8-1, ny*5/8-1);
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, topLow, boundary::freeslip );
Box2D topSide(nx/2-1, nx/2-1, ny*5/8-1 , ny-1);
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, topSide, boundary::freeslip );
Box2D topWall(nx/2-1, nx-1, ny-1, ny-1);
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, topWall, boundary::freeslip );
Box2D bottomLow(0, nx/2-1, ny*3/8-1, ny*3/8-1);
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, bottomLow, boundary::freeslip );
Box2D bottomSide(nx/2-1, nx/2-1, ny*3/8-1, 0);
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, bottomSide, boundary::freeslip );
Box2D bottomWall(nx/2-1, nx-1, 0, 0);
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, bottomWall, boundary::freeslip );
Box2D inlet(0, 0, ny*3/8+2, ny*5/8+2);
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, inlet );
Box2D outlet(nx-1, nx-1, 2, ny-2);
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, outlet); //, boundary::outflow );
boundaryCondition.addInternalVelocityCornerPN( nx/2, ny*5/8, lattice, boundary::freeslip );
boundaryCondition.addInternalVelocityCornerPP( nx/2, ny*3/8, lattice, boundary::freeslip );
boundaryCondition.addExternalVelocityCornerPN( nx/2, ny, lattice, boundary::freeslip );
boundaryCondition.addExternalVelocityCornerPP( nx/2, 0, lattice, boundary::freeslip );
// Specify the boundary velocity.
setBoundaryVelocity (
lattice, lattice.getBoundingBox(),
PoiseuilleVelocity<T>(parameters) );
// Create the initial condition.
initializeAtEquilibrium (
lattice, lattice.getBoundingBox(), PoiseuilleVelocity<T>(parameters) );
lattice.initialize();
}
void writeGifs(MultiBlockLattice2D<T,DESCRIPTOR>& lattice, plint iter)
{
const plint imSize = 600;
ImageWriter imageWriter(“leeloo”);
imageWriter.writeScaledGif(createFileName(“u”, iter, 6),
*computeVelocityNorm(lattice),
imSize, imSize );
}
int main(int argc, char* argv[]) {
plbInit(&argc, &argv);
global::directories().setOutputDir("./tmp/");
// Use the class IncomprFlowParam to convert from
// dimensionless variables to lattice units, in the
// context of incompressible flows.
IncomprFlowParam<T> parameters(
(T) 1e-1, // Reference velocity (the maximum velocity
// in the Poiseuille profile) in lattice units.
(T) 100., // Reynolds number
300, // Resolution of the reference length (channel height).
3., // Channel length in dimensionless variables
1. // Channel height in dimensionless variables
);
const T imSave = (T)0.1; // Time intervals at which to save GIF
// images, in dimensionless time units.
const T maxT = (T)3.1; // Total simulation time, in dimensionless
// time units.
writeLogFile(parameters, "Poiseuille flow");
MultiBlockLattice2D<T, DESCRIPTOR> lattice (
parameters.getNx(), parameters.getNy(),
new BGKdynamics<T,DESCRIPTOR>(parameters.getOmega()) );
OnLatticeBoundaryCondition2D<T,DESCRIPTOR>*
boundaryCondition = createLocalBoundaryCondition2D<T,DESCRIPTOR>();
channelSetup(lattice, parameters, *boundaryCondition);
// Main loop over time iterations.
for (plint iT=0; iT*parameters.getDeltaT()<maxT; ++iT) {
if (iT%parameters.nStep(imSave)==0 && iT>0) {
pcout << "Saving Gif at time step " << iT << endl;
writeGifs(lattice, iT);
}
// Execute lattice Boltzmann iteration.
lattice.collideAndStream();
}
delete boundaryCondition;
}