Welcome! Log In Create A New Profile

Advanced

Fluid around a flat plate 2D

Posted by Guillermo55Guillermo55  
Fluid around a flat plate 2D
September 14, 2016 02:47PM
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
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.