Welcome! Log In Create A New Profile

Advanced

Microchannel with a pressure gradient at Palabos

Posted by mfduquedmfduqued  
Microchannel with a pressure gradient at Palabos
May 08, 2018 09:04PM
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++
  1. #include"palabos3D.h"
  2. #include"palabos3D.hh"
  3.  
  4. #include <cstdlib>
  5. #include <vector>
  6. #include <cmath>
  7. #include <iostream>
  8. #include <sstream>
  9. #include <fstream>
  10. #include <iomanip>
  11.  
  12. using namespace plb;
  13. using namespace std;
  14. using namespace plb::descriptors;
  15.  
  16. typedef double T;
  17. #define DESCRIPTOR D3Q27Descriptor
  18. #define DYNAMICS BGKdynamics<T, DESCRIPTOR>(parameters.getOmega())
  19. // // Aca definimos el gradiente de presion
  20. class PressureGradient {
  21. public:
  22. PressureGradient(T deltaP_, plint nx_) : deltaP(deltaP_), nx(nx_)
  23. { }
  24. void operator() (plint iX, plint iY, plint iZ, T& density, Array<T,3>& velocity) const
  25. {
  26. velocity.resetToZero();
  27. density = (T)1 - deltaP*DESCRIPTOR<T>::invCs2 / (T)(nx-1) * (T)iX;
  28.  
  29. }
  30. private:
  31. T deltaP;
  32. plint nx;
  33. };
  34.  
  35. void channelSetup(MultiBlockLattice3D<T,DESCRIPTOR>& lattice,
  36. IncomprFlowParam<T> const& parameters,
  37. OnLatticeBoundaryCondition3D<T,DESCRIPTOR>& boundaryCondition, T deltaP)
  38. {
  39. const plint nx = parameters.getNx();
  40. const plint ny = parameters.getNy();
  41. const plint nz = parameters.getNz();
  42.  
  43. // // Microchannel';s Geometry
  44. //==================================================================
  45.  
  46. Box3D bottom(0, nx-1, 0, 0, 0, nz-1);
  47.  
  48. Box3D top (0, nx-1, ny-1, ny-1, 0, nz-1);
  49.  
  50. Box3D right (0, nx-1, 0, ny-1, nz-1, nz-1);
  51.  
  52. Box3D left (0, nx-1, 0, ny-1, 0, 0);
  53.  
  54. Box3D outlet(nx-1, nx-1, 1, ny-2, 1, nz-2);
  55.  
  56. Box3D inlet (0, 0, 1, ny-2, 1, nz-2);
  57. // //boudary====================================================
  58. /
  59.  
  60. boundaryCondition->addPressureBoundary0N(inlet, lattice);
  61. setBoundaryDensity(lattice, inlet, (T) 1.);
  62. // initializeAtEquilibrium(lattice, inlet, PressureGradient(deltaP, nx)
  63. // );
  64.  
  65. boundaryCondition->addPressureBoundary0P(outlet, lattice);
  66. setBoundaryDensity(lattice, outlet, (T) 1. - deltaP*DESCRIPTOR<T>::invCs2);
  67. initializeAtEquilibrium(lattice, lattice.getBoundingBox(), PressureGradient(deltaP, nx)
  68. );
  69.  
  70. lattice.initialize();
  71. delete boundaryCondition;
  72. }
  73.  
  74. template<class BlockLatticeT>
  75. void writeVTK(BlockLatticeT& lattice, IncomprFlowParam<T> const& parameters, plint iter){
  76. T dx = parameters.getDeltaX();
  77. T dt = parameters.getDeltaT();
  78. VtkImageOutput3D<T> vtkOut(createFileName("vtk", iter, 6), dx);
  79. vtkOut.writeData<float>(*computeVelocityNorm(lattice), "velocityNorm", dx/dt);
  80. vtkOut.writeData<float>(*computeDensity(lattice), "density", 1.);
  81. }
  82.  
  83. T computeVmedia (
  84. MultiBlockLattice3D<T,DESCRIPTOR>& lattice,
  85. IncomprFlowParam<T> const& parameters , Box3D domain )
  86. {
  87. T meanU = computeAverage ( *computeVelocityNorm (lattice, domain) );
  88.  
  89. T meanUphys = (parameters.getDeltaX() / parameters.getDeltaT()) * meanU;
  90.  
  91. return meanUphys;
  92. }
  93.  
  94. int main(int argc, char **argv){
  95. plbInit(&argc,&argv);
  96.  
  97. global::directories().setOutputDir("./tmp/");
  98.  
  99. IncomprFlowParam<T> parameters(
  100. (T) 0.63746, //Velocity at physical units
  101. (T) 0.03, //Velocity LB
  102. (T) 0.025, //Knudsen';s number
  103. (T) 1e-6, // caractheristc Lenght
  104. (T) 16., //Resolution
  105. (T) 18, //times lenght x
  106. (T) 1, //times lenght y
  107. (T) 3 //times lenght z
  108. );
  109.  
  110. T gradP = (T)1;
  111. const T maxT = (T)39e-6;
  112. const T imSave = (T)maxT/(T)10;
  113. T deltaP = (T)gradP * parameters.getDeltaX();
  114.  
  115. writeLogFile(parameters, "Parameters_Microchannel");
  116.  
  117. plint nx = parameters.getNx();
  118. plint ny = parameters.getNy();
  119. plint nz = parameters.getNz();
  120.  
  121. MultiBlockLattice3D<T, DESCRIPTOR> lattice(
  122. nx,ny,nz, new DYNAMICS);
  123.  
  124. util::ValueTracer<T> convergen(parameters.getPhysicalU(), parameters.getPhysicalLength(), 1e-5);
  125.  
  126. lattice.periodicity().toggleAll(false);
  127.  
  128. OnLatticeBoundaryCondition3D<T,DESCRIPTOR>*
  129. boundaryCondition = createLocalBoundaryCondition3D<T,DESCRIPTOR>();
  130.  
  131. channelSetup(lattice, parameters, *boundaryCondition, deltaP);
  132.  
  133. plb_ofstream UPvsT("./tmp/VvsT_320.dat");
  134.  
  135.  
  136. Box3D profileSectionFinal(nx/2, nx/2, ny/2, ny/2,0, nz - 1);//A lo largo del eje Z
  137.  
  138. Box3D profileSectionUP1(nx/2, nx/2, ny/2, ny/2, nz/2, nz/2);
  139.  
  140. plint counter = 0;
  141.  
  142. T unit_velo = parameters.getDeltaX() / parameters.getDeltaT() / parameters.getLatticeU();
  143.  
  144.  
  145. for(plint iT=0; iT*parameters.getDeltaT()<1.01*maxT; ++iT){
  146. counter += 1;
  147. double TimE = iT*parameters.getDeltaT();
  148.  
  149. lattice.collideAndStream();
  150. convergencia.takeValue(getStoredAverageEnergy(lattice),true);
  151.  
  152. if(convergen.hasConverged()){
  153. pcout << "Simulation respect to kinetic energy mean" << endl;
  154. pcout << "Simulation';s Time =" << TimE << endl;
  155. break;
  156. }
  157.  
  158. UPvsT << setprecision(16)
  159. << TimE << " " << *multiply(*computeVelocityNorm (lattice, profileSectionUP1),unit_velo) << endl;
  160. }
  161.  
  162. T Vmedia = computeVmedia(lattice, parameters, profileSectionFinal);
  163.  
  164. pcout <<"mean velocity ="<<" "<< Vmedia << endl;
  165.  
  166. return 0;
  167. }



Edited 1 time(s). Last edit at 05/09/2018 02:32AM by mfduqued.
Sorry, you do not have permission to post/reply in this forum.