# Microchannel with a pressure gradient at Palabos

Posted by mfduqued
 Microchannel with a pressure gradient at Palabos May 09, 2018 02:34AM Registered: 3 years ago Posts: 20
Dear all,
I am a PhD student in mechanical engineering, working on a problem that requires a microchannel flow field.

In these problems, compressibility effects are present even at low speeds. In such a way that the density is not constant, therefore the pressure gradient can not be represented as a force body.

On the other hand, the relaxation time is not related to the viscosity (by means of Reynolds' number), but to the Knudsen's number, additionally there is a slip condition of the fluid with respect to the wall.

Therefore, I modified a file in Palabos (unit.h), so that the relaxation time is obtained by means of Knudsen's number. Already verify this change, with benchmark problems. Observing how to work with the concept of pressure gradient in the example permeability.cpp, I implemented it to work with the Poiseuille problem and a pressure gradient, but it does not work.

Thanks to everyone for reading, I appreciate any contribution. It can be very short, that is, "read this document" or "looking at the source code that I append".

Thank you very much,

Mauricio

`Language: C++#include"palabos3D.h"#include"palabos3D.hh" #include <cstdlib>#include <vector>#include <cmath>#include <iostream>#include <sstream>#include <fstream>#include <iomanip> using namespace plb;using namespace std;using namespace plb::descriptors; typedef double T;#define DESCRIPTOR D3Q27Descriptor#define DYNAMICS BGKdynamics<T, DESCRIPTOR>(parameters.getOmega())// // Aca definimos el gradiente de presionclass PressureGradient {public:    PressureGradient(T deltaP_, plint nx_) : deltaP(deltaP_), nx(nx_)    { }    void operator() (plint iX, plint iY, plint iZ, T& density, Array<T,3>& velocity) const    {        velocity.resetToZero();        density = (T)1 - deltaP*DESCRIPTOR<T>::invCs2 / (T)(nx-1) * (T)iX;     }private:    T deltaP;    plint nx;}; void channelSetup(MultiBlockLattice3D<T,DESCRIPTOR>& lattice,                  IncomprFlowParam<T> const& parameters,                  OnLatticeBoundaryCondition3D<T,DESCRIPTOR>& boundaryCondition, T deltaP){      	  const plint nx = parameters.getNx();  	  const plint ny = parameters.getNy();  	  const plint nz = parameters.getNz(); 	  // // Microchannel';s Geometry	//================================================================== 	        Box3D bottom(0, nx-1, 0,    0,      0,   nz-1);         	Box3D top   (0, nx-1, ny-1, ny-1,   0,   nz-1);  	        Box3D right (0,  nx-1, 0,   ny-1,  nz-1, nz-1);          	Box3D left  (0,  nx-1, 0,   ny-1,  0,   0);  	        Box3D outlet(nx-1, nx-1, 1,   ny-2,  1,   nz-2);          	Box3D inlet (0,   0,    1,  ny-2,  1,   nz-2);// 	//boudary====================================================/     boundaryCondition->addPressureBoundary0N(inlet, lattice);    setBoundaryDensity(lattice, inlet, (T) 1.);//    initializeAtEquilibrium(lattice, inlet, PressureGradient(deltaP, nx)//    );     boundaryCondition->addPressureBoundary0P(outlet, lattice);    setBoundaryDensity(lattice, outlet, (T) 1. - deltaP*DESCRIPTOR<T>::invCs2);    initializeAtEquilibrium(lattice, lattice.getBoundingBox(), PressureGradient(deltaP, nx)    );     lattice.initialize();    delete boundaryCondition;} template<class BlockLatticeT>void writeVTK(BlockLatticeT& lattice, IncomprFlowParam<T> const& parameters, plint iter){    T dx = parameters.getDeltaX();    T dt = parameters.getDeltaT();    VtkImageOutput3D<T> vtkOut(createFileName("vtk", iter, 6), dx);    vtkOut.writeData<float>(*computeVelocityNorm(lattice), "velocityNorm", dx/dt);    vtkOut.writeData<float>(*computeDensity(lattice), "density", 1.);} T computeVmedia (    MultiBlockLattice3D<T,DESCRIPTOR>& lattice,                  IncomprFlowParam<T> const& parameters , Box3D domain ){    T meanU =  computeAverage ( *computeVelocityNorm (lattice, domain) );       T meanUphys = (parameters.getDeltaX() / parameters.getDeltaT()) *  meanU;     return meanUphys;} int main(int argc, char **argv){  plbInit(&argc,&argv);   global::directories().setOutputDir("./tmp/");   IncomprFlowParam<T> parameters(    (T)  0.63746,               //Velocity at physical units    (T)  0.03,                      //Velocity LB    (T)  0.025,                   //Knudsen';s number      (T)  1e-6,                     // caractheristc Lenght    (T)  16.,                       //Resolution    (T)  18,                        //times lenght x     (T)  1,			        //times lenght y    (T)  3			        //times lenght z    );   T gradP		= (T)1;  const T maxT          = (T)39e-6;  const T imSave        = (T)maxT/(T)10;T deltaP		= (T)gradP * parameters.getDeltaX();   writeLogFile(parameters, "Parameters_Microchannel");         plint nx = parameters.getNx();        plint ny = parameters.getNy();	plint nz = parameters.getNz(); 	MultiBlockLattice3D<T, DESCRIPTOR> lattice(				      nx,ny,nz, new DYNAMICS); 	util::ValueTracer<T> convergen(parameters.getPhysicalU(), parameters.getPhysicalLength(), 1e-5);       lattice.periodicity().toggleAll(false);   OnLatticeBoundaryCondition3D<T,DESCRIPTOR>*        boundaryCondition = createLocalBoundaryCondition3D<T,DESCRIPTOR>();     channelSetup(lattice, parameters, *boundaryCondition, deltaP);        plb_ofstream UPvsT("./tmp/VvsT_320.dat");         Box3D profileSectionFinal(nx/2, nx/2, ny/2, ny/2,0, nz - 1);//A lo largo del eje Z 	Box3D profileSectionUP1(nx/2, nx/2, ny/2, ny/2, nz/2, nz/2); 	plint counter = 0; 	T unit_velo = parameters.getDeltaX() / parameters.getDeltaT() / parameters.getLatticeU();    for(plint iT=0; iT*parameters.getDeltaT()<1.01*maxT; ++iT){	counter += 1;        double TimE = iT*parameters.getDeltaT();             lattice.collideAndStream(); 	convergencia.takeValue(getStoredAverageEnergy(lattice),true); 	if(convergen.hasConverged()){	  pcout << "Simulation respect to kinetic energy mean" << endl;	  pcout << "Simulation';s Time =" << TimE << endl;	  break; 	}         UPvsT << setprecision(16)	  	      << TimE << "    " << *multiply(*computeVelocityNorm (lattice, profileSectionUP1),unit_velo) << endl;  }	      T Vmedia = computeVmedia(lattice, parameters, profileSectionFinal);   pcout <<"mean velocity ="<<"   "<< Vmedia << endl;   return 0;}`
Sorry, you do not have permission to post/reply in this forum.