Blasius Flat Plate

Guillermo55 · September 14, 2016, 2:44pm

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;

}




Guillermo Izquierdo

kk1 · September 15, 2016, 6:54am

Dear Guillermo Izquierdo,

for others to understand what the problem is, you should give some information about what exactly went wrong (error messages or details about the output of the simulation).

Concerning the periodic boundary conditions, the syntax should be

[code=“cpp”]
lattice.periodicity().toggle(0, true);



instead of

[code="cpp"]
lattice.periodicity.toggle(0, true);

Regards,

kk

Guillermo55 · September 20, 2016, 11:01am

Thanks for your response.
I give for details of the problem and let see if you can help me because I do not know what else to do.

The blasius flat plate solution (analytical one) gives:
n=ysqrt(vx/U)
f’(n) = u/U
Being u=velocity in the boundary layer; U=velocity outside boundary layer (inlet velocity); v=kinematic viscosity; x=characteristic length
I introduce those parameters in the Incompressible Flow Parameters of Palabos:
Re=25
uMax=0,003
N=1600
lx=1
ly=1
which gives me a v (kinematic viscosity) of 0,192
And then I take the velocities in the middle of the flat plate, which is situated in the bottom wall. I have then, u and y for LBM.
I compare both results by changing Blasius flat Plate units to LBM by:
y=nsqrt(U/(vx))=nsqrt(0,003/(0,192(1600/2)))
u=f’(n)*U=f’(n)*0,003
And I have u and y for Blasius.
I obtain the graph y vs u in both cases and they do not look at all similar.

I would be very pleased if you could help me, because I have tried everything by nothing changes.

I post the code of LBM for Blasius Flat Plate.