# Fluid around a flat plate 2D

Dear community of Palabos,

I am in the process of understanding your code for academic purpose. I want to compare the result that gives Palabos with the analytical solution (that of Blasius) in the case of an external fluid around a Flat Plate. I am doing it in 2D. Theoretically it should be quite simple, but there is something I am not doing well and I don’t know why. Perhaps boundaries conditions and how to make them periodic, don’t know. I would be grateful if you could help me. Thank you in advance, here is the code:

[code=“cpp”]
#include “palabos2D.h”
#ifndef PLB_PRECOMPILED // Unless precompiled version is used,
#include “palabos2D.hh” // include full template code
#endif
#include
#include
#include
#include
#include

using namespace plb;
using namespace plb::descriptors;
using namespace std;

typedef double T;
#define DESCRIPTOR D2Q9Descriptor

template
class LineShapeDomain2D : public DomainFunctional2D {
public:
LineShapeDomain2D(plint cx_, plint cy_, plint ll_)
: cx(cx_),
cy(cy_),
ll(ll_)
{ }
virtual bool operator() (plint iX, plint iY) const {
return (iY==cy) and (iX>=cx and iX<=(cx+ll));
}
virtual LineShapeDomain2D* clone() const {
return new LineShapeDomain2D(*this);
}
private:
plint cx;
plint cy;
plint ll;
};

void boundariesAndPlatSetup( MultiBlockLattice2D<T,DESCRIPTOR>& lattice,
IncomprFlowParam const& parameters,
OnLatticeBoundaryCondition2D<T,DESCRIPTOR>& boundaryCondition )
{
const plint nx = parameters.getNx();
const plint ny = parameters.getNy();
T u = parameters.getLatticeU();

``````lattice.periodicity.toggle(0, true);

// inlet
Box2D inlet(0,0,0,ny-1);
boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice,inlet);
setBoundaryVelocity(lattice,inlet,Array<T,2>(u,0.));

// outlet
Box2D outlet(nx-1,nx-1,0,ny-1);
boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice,outlet,boundary::outflow);

// top and bottom wall
Box2D bottomWall(1,nx-2,0,0);
Box2D topWall(1,nx-2,ny-1,ny-1);
boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice,topWall,boundary::freeslip);
boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice,bottomWall,boundary::freeslip);

initializeAtEquilibrium(lattice, lattice.getBoundingBox(), 1., Array<T,2>(u,0.) );

defineDynamics(lattice,lattice.getBoundingBox(), new LineShapeDomain2D<T>(nx/4,ny/2,nx/2), new BounceBack<T,DESCRIPTOR>);

lattice.initialize();
``````

}

template
void writeGif(BlockLatticeT& lattice, plint iter)
{
const plint imSize = 600;

``````ImageWriter<T> 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/");

IncomprFlowParam<T> parameters(
(T) 0.05,  // uMax
(T) 10000.,  // Re
1.,       // N
500.,        // lx
500.         // ly
);
const T logT     = (T)0.1;
const T imSave   = (T)1.;
const T vtkSave  = (T)1.;
const T maxT     = (T)20.2;

writeLogFile(parameters, "solo 1 pared");

MultiBlockLattice2D<T, DESCRIPTOR> lattice (
parameters.getNx(), parameters.getNy(),
new BGKdynamics<T,DESCRIPTOR>(parameters.getOmega()) );

lattice.periodicity().toggle(0,false);

OnLatticeBoundaryCondition2D<T,DESCRIPTOR>*
boundaryCondition = createZouHeBoundaryCondition2D<T,DESCRIPTOR>();

boundariesAndPlatSetup(lattice, parameters, *boundaryCondition);

T previousIterationTime = T();

// Main loop over time iterations.
for (plint iT=0; iT<5100; ++iT) {
global::timer("mainLoop").restart();

if (iT%parameters.nStep(imSave)==0 && iT>0) {
pcout << "Saving Gif ..." << endl;
writeGif(lattice, iT);
pcout << endl;
}

lattice.collideAndStream();

}
writeGif(lattice, 1);

// velocities at 0.50 plat
const plint nx = parameters.getNx();
const plint ny = parameters.getNy();
plb_ofstream ofile("velocities.dat");
for (plint i=0;i<ny/2-1;i++){
Array<T,2> velocity;
lattice.get(nx/2,i+ny/2).computeVelocity(velocity);
ofile << setprecision(3) << velocity[0] << endl;
}
pcout<<"nx: "<<nx<<"   ny: "<<ny<<endl;

delete boundaryCondition;
``````

}

Hello,

I just saw your post. I have already written a simple code to simulate flat plate. I hope it helps.

[code=“cpp”]
#include “palabos2D.h”
#include “palabos2D.hh”

using namespace plb;
using namespace std;

typedef double T;
#define DESCRIPTOR descriptors::D2Q9Descriptor

T poiseuilleVelocity(plint iY, IncomprFlowParam const& parameters) {
T y = (T)iY / parameters.getResolution();
return 4.parameters.getLatticeU() * (y-yy);
}

template
class PoiseuilleSetup {
public:
PoiseuilleSetup(IncomprFlowParam parameters_)
: parameters(parameters_) { }
void operator()(plint iX, plint iY, T& rho, Array<T,2>& u) const {
rho = 1.0;
u[0] = poiseuilleVelocity(iY, parameters);
u[1] = 0.0;
}
void operator()(plint iX, plint iY, Array<T,2>& u) const {
u[0] = poiseuilleVelocity(iY, parameters);
u[1] = 0.0;
}

private:
IncomprFlowParam parameters;
};

void channelSetup( MultiBlockLattice2D<T,DESCRIPTOR>& lattice,
IncomprFlowParam const& parameters,
OnLatticeBoundaryCondition2D<T,DESCRIPTOR>& boundaryCondition)
{
const plint nx = parameters.getNx();
const plint ny = parameters.getNy();

``````boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice);
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, Box2D(nx-1, nx, 1, ny-2), boundary::outflow );

setBoundaryVelocity( lattice, lattice.getBoundingBox(),
PoiseuilleSetup<T>(parameters) );

initializeAtEquilibrium ( lattice, lattice.getBoundingBox(),
PoiseuilleSetup<T>(parameters) );

lattice.initialize();
``````

}

/// Produce a GIF snapshot of the velocity-norm.
void writeGif(MultiBlockLattice2D<T,DESCRIPTOR>& lattice, plint iter)
{
ImageWriter imageWriter(“leeloo”);
imageWriter.writeScaledGif(createFileName(“u”, iter, 6),
*computeVelocityNorm(lattice));
}

int main(int argc, char* argv[]) {
plbInit(&argc, &argv);

``````global::directories().setOutputDir("./tmp/");

IncomprFlowParam<T> parameters(
(T) 1e-2,  // uMax
(T) 600.,  // Re
60,       // N
5.,        // lx
1.         // ly
);
const T logT     = (T)0.02;
const T imSave   = (T)0.06;
const T maxT     = (T)20.1;

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>();

// Flat Plate Geometry //////////////
plint X0 = parameters.getNx()/10, W = 20;
plint Y0 = parameters.getNy()/2+2 , H = 2;
Box2D FlatPlate( X0, X0+W, Y0-H/2, Y0+H/2);
defineDynamics(lattice, FlatPlate, new BounceBack<T,DESCRIPTOR>);

channelSetup(lattice, parameters, *boundaryCondition);

// Main loop over time iterations.
for (plint iT=0; iT*parameters.getDeltaT()<maxT; ++iT) {
if (iT%parameters.nStep(logT)==0) {
pcout << "step " << iT
<< "; t=" << iT*parameters.getDeltaT()<<endl;
}
if (iT%parameters.nStep(imSave)==0) {
pcout << "Saving Gif ...\n" << endl;
writeGif(lattice, iT);
}
// Lattice Boltzmann iteration step.
lattice.collideAndStream();
}

delete boundaryCondition;
``````

}

``````

BTW, please create a topic  in an appropriate place, all the posts here should be related to Palabos software bugs.

The Best,
Arman``````