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;
}