Hi All,
I was trying out some palabos feature and I tried to simulate poiseuille flow in a channell without following the palabos tutorial. Basically I imposed a velocity in the inlet-outlet walls and a bouce back dynamic on the top and bottom wall. It should be straight forward but gives me strange results.
Here the code:
[code=“cpp”]
#include
#include “palabos2D.h”
#include “palabos2D.hh”
using namespace std;
using namespace plb;
typedef double T;
#define NSDES descriptors::D2Q9Descriptor
int main(int argc, char *argv[])
{
// Init simulation
plbInit(&argc, &argv);
// Set omegas and constant velocity
const T nsOmega = 1.0;
Array<T,2> u0(1.0,0.0);
// Allocate lattice
MultiBlockLattice2D<T, NSDES> nsLattice(100, 10, new BGKdynamics<T, NSDES>(nsOmega));
// Create box2D for the boundaries
const plint nx = nsLattice.getNx();
const plint ny = nsLattice.getNy();
Box2D topWall(0, nx, ny-1,ny);
Box2D bottomWall(0,nx, 0,1);
Box2D inlet(0,1, 1,ny -1);
Box2D outlet(nx-1,nx, 1,ny-1);
// Boundary conditions for NS dynamics
OnLatticeBoundaryCondition2D<T, NSDES> *nsBC = createLocalBoundaryCondition2D<T, NSDES>();
// Set the boundary conditions type
nsBC->setVelocityConditionOnBlockBoundaries(nsLattice, inlet);
nsBC->setVelocityConditionOnBlockBoundaries(nsLattice, outlet);
defineDynamics(nsLattice, topWall, new BounceBack<T, NSDES>);
defineDynamics(nsLattice, bottomWall, new BounceBack<T, NSDES>);
// Set the velocity value
setBoundaryVelocity(nsLattice, inlet, u0);
setBoundaryVelocity(nsLattice, outlet, u0);
// Init lattice
initializeAtEquilibrium(nsLattice, nsLattice.getBoundingBox(), 1.0, u0);
nsLattice.initialize();
//Write first image
ImageWriter<T> image("leeloo");
image.writeScaledGif(createFileName("vel",1,1), *computeVelocityNorm(nsLattice), nx,ny);
// Main LBM cycle
for (plint iT=0; iT < 10000; ++iT) {
if(iT%100==0) {
pcout << "Writing GIF...";
image.writeScaledGif(createFileName("vel",iT,2), *computeVelocityNorm(nsLattice), nx,ny);
}
nsLattice.collideAndStream();
}
pcout << "Hello palabos";
return 0;
}
Any suggestions? Is there smth wrong with the code?