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
[code=“cpp” number]
#include"palabos3D.h"
#include"palabos3D.hh"
#include
#include
#include
#include
#include
#include
#include
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 presion
class 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::invCs2 / (T)(nx-1) * (T)iX;
}
private:
T deltaP;
plint nx;
};
void channelSetup(MultiBlockLattice3D<T,DESCRIPTOR>& lattice,
IncomprFlowParam 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
void writeVTK(BlockLatticeT& lattice, IncomprFlowParam const& parameters, plint iter){
T dx = parameters.getDeltaX();
T dt = parameters.getDeltaT();
VtkImageOutput3D vtkOut(createFileName(“vtk”, iter, 6), dx);
vtkOut.writeData(*computeVelocityNorm(lattice), “velocityNorm”, dx/dt);
vtkOut.writeData(*computeDensity(lattice), “density”, 1.);
}
T computeVmedia (
MultiBlockLattice3D<T,DESCRIPTOR>& lattice,
IncomprFlowParam 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 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; iTparameters.getDeltaT()<1.01maxT; ++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;
}