# Fluid around a flat plate 2D

Posted by Guillermo55
 Fluid around a flat plate 2D September 14, 2016 02:47PM Registered: 3 years ago Posts: 8
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:

```Language: C++#include "palabos2D.h"
#ifndef PLB_PRECOMPILED // Unless precompiled version is used,
#include "palabos2D.hh"   // include full template code
#endif
#include <vector>
#include <cmath>
#include <iostream>
#include <fstream>
#include <iomanip>

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

typedef double T;
#define DESCRIPTOR D2Q9Descriptor

template<typename T>
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<T>* clone() const {
return new LineShapeDomain2D<T>(*this);
}
private:
plint cx;
plint cy;
plint ll;
};

void boundariesAndPlatSetup( MultiBlockLattice2D<T,DESCRIPTOR>& lattice,
IncomprFlowParam<T> 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<class BlockLatticeT>
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;
}```
 Re: Fluid around a flat plate 2D December 05, 2016 10:43AM Registered: 8 years ago Posts: 26
Hello,

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

```Language: C++#include "palabos2D.h"
#include "palabos2D.hh"

using namespace plb;
using namespace std;

typedef double T;
#define DESCRIPTOR descriptors::D2Q9Descriptor

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

template<typename T>
class PoiseuilleSetup {
public:
PoiseuilleSetup(IncomprFlowParam<T> 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<T> parameters;
};

void channelSetup( MultiBlockLattice2D<T,DESCRIPTOR>& lattice,
IncomprFlowParam<T> 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<T> 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
Sorry, you do not have permission to post/reply in this forum.