Welcome! Log In Create A New Profile

Advanced

Consistency law at Lattice Boltzmann

Posted by mfduquedmfduqued  
Consistency law at Lattice Boltzmann
October 30, 2017 01:31PM
I am a student PhD at this moment I work with palabos for modeling a channel rectangular tridimensional, whose initial condition is V=0cm/s, and affected for a force. I used a model D3Q19, and did several cases.

My question is, when I modified resolution number's without changed my problem physical (physical velocity, Reynolds number and characteristic length) my velocity profile changed.

According to consistency of problem, this doesn't must change the solution because I would to change the resolution number, but that the effect of this change is increasing accuracy.

Please community someone can tell me why do it happen this, or where can I found literature about of this topic.

Thank you
Re: Consistency law at Lattice Boltzmann using palabos
October 30, 2017 01:34PM
I am a student PhD at this moment I work with palabos for modeling a channel rectangular tridimensional, whose initial condition is V=0cm/s, and affected for a force. I used a model D3Q19, and did several cases.

My question is, when I modified resolution number's without changed my problem physical (physical velocity, Reynolds number and characteristic length) my velocity profile changed.

According to consistency of problem, this doesn't must change the solution because I would to change the resolution number, but that the effect of this change is increasing accuracy.

Please community someone can tell me why do it happen this, or where can I found literature about of this topic.

Thank you

I attache my code
******************************************************************
Language: C++
[code="cpp" number] #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 ForcedD3Q19Descriptor #define DYNAMICS GuoExternalForceBGKdynamics   T channelForce(IncomprFlowParam<T> const& parameters) { return (T) 6e-6; }   void channelSetup(MultiBlockLattice3D<T,DESCRIPTOR>& lattice, IncomprFlowParam<T> const& parameters, OnLatticeBoundaryCondition3D<T,DESCRIPTOR>& boundaryCondition) { const plint nx = parameters.getNx(); const plint ny = parameters.getNy(); const plint nz = parameters.getNz();   ============================ 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, 1, ny-2, nz-1, nz-1); Box3D left (0, nx-1, 1, ny-2, 0, 0); Box3D outlet(nx-1, nx-1, 1, ny-2, 0, nz-1); Box3D inlet (0, 0, 1, ny-2, 0, nz-1); ==================================== boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice, inlet, boundary::outflow); initializeAtEquilibrium(lattice, inlet, (T)1., Array<T,3>((T)0.,(T)0.,(T)0.) );   boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice, outlet, boundary::outflow); initializeAtEquilibrium(lattice, outlet, (T)1., Array<T,3>((T)0.,(T)0.,(T)0.) );   boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice, bottom ); setBoundaryVelocity(lattice, bottom, Array<T,3>((T)0.,(T)0.,(T)0.)); initializeAtEquilibrium(lattice, bottom, (T)1., Array<T,3>((T)0.,(T)0.,(T)0.) );   boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice, left, boundary::outflow); initializeAtEquilibrium(lattice, left, (T)1., Array<T,3>((T)0.,(T)0.,(T)0.) );   boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice, right, boundary::outflow); initializeAtEquilibrium(lattice, right, (T)1., Array<T,3>((T)0.,(T)0.,(T)0.) );   boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice, top ); setBoundaryVelocity(lattice, top, Array<T,3>((T)0.,(T)0., (T)0.) ); initializeAtEquilibrium(lattice, top, (T)1., Array<T,3>((T)0,(T)0., (T) 0.) );   Array<T,3> force(channelForce(parameters), (T)0.,(T)0); setExternalVector(lattice,lattice.getBoundingBox(),DESCRIPTOR<T>::ExternalField::forceBeginsAt,force);   lattice.initialize(); }   int main(int argc, char*argv[]){ plbInit(&argc,&argv);   global::directories().setOutputDir("./tmp/");   IncomprFlowParam<T> parameters( (T) 3.350e-2, //Characteristic Velocity (S.I. units) (T) 0.095, //Lattice Boltzmann Velocity (T) 1000., //Reynolds (T) 0.3, //Characteristic Lenght (T) 48., //resolution number';s (T) 3., //Lx (T) 1., //Ly (T) 1. //Lz );   const T maxT = (T)4500; const T imSave = (T)maxT/(T)100; const T promedio = (T)1110;   writeLogFile(parameters, "Parameters_Channel MRT");   plint nx = parameters.getNx(); plint ny = parameters.getNy(); plint nz = parameters.getNz();   MultiBlockLattice3D<T, DESCRIPTOR> lattice ( nx, ny, nz, new DYNAMICS<T, DESCRIPTOR>(parameters.getOmega()));   lattice.periodicity().toggle(0,true); lattice.periodicity().toggle(2,true);   OnLatticeBoundaryCondition3D<T,DESCRIPTOR>* boundaryCondition = createLocalBoundaryCondition3D<T,DESCRIPTOR>();   channelSetup(lattice, parameters, *boundaryCondition);   Box3D Cavity(0, nx-1, 0, ny-1, 0, nz-1);   plb_ofstream UPvsT("./tmp/VvsT_320.dat");     Box3D profileSectionHalf(0, nx-1, ny/2, ny/2,nz/2, nz/2); //axis X Box3D profileSectionMit(nx/2, nx/2, 0, ny - 1,nz/2, nz/2); //axis Y Box3D profileSectionFinal(nx/2, nx/2, ny/2, ny/2,0, nz - 1);//axis Z   Box3D profileSectionUP1(nx/2, nx/2, ny/2, ny/2, nz/2, nz/2);     Box3D profileSectionIniz(0, nx-1, ny - 1, ny - 1, 0, nz-1); Box3D profileSectionMitz(0, nx-1, ny - 1, ny - 1, 0, nz-1); Box3D profileSectionFinalz(0, nx-1, ny - 1, ny - 1, 0, nz-1); Box3D profileSectionDownz(0, nx-1, ny/4, ny/4, 0, nz-1); Box3D profileSectionHalfz(0, nx-1, ny/2, ny/2, 0, nz-1); Box3D profileSectionUpz(0, nx-1, 3*ny/4, 3*ny/4, 0, nz-1);   Box3D StrengthTapUp(0,nx-1,ny-1,ny-1,0,nz-1); Box3D StrengthTapDown(0,nx-1,0,0,0,nz-1);   plint counter = 0;   MultiScalarField3D<T> Velocity_iter(nx,ny,nz); MultiScalarField3D<T> Velocidad_iter(nx,ny,nz); MultiScalarField3D<T> ShearUP_iter(nx,ny,nz); MultiScalarField3D<T> EsfuerzoUP_iter(nx,ny,nz); MultiScalarField3D<T> ShearDown_iter(nx,ny,nz); MultiScalarField3D<T> EsfuerzoDown_iter(nx,ny,nz);   setToConstant(Velocidad_iter, Cavity, (T)0.); setToConstant(EsfuerzoUP_iter, Cavity, (T)0.); setToConstant(EsfuerzoDown_iter, Cavity, (T)0.);   T unit_Nu = parameters.getLatticeNu() * parameters.getDeltaX() * parameters.getDeltaX() / parameters.getDeltaT();   for(plint iT=0; iT*parameters.getDeltaT()<1.01*maxT; ++iT){ counter += 1; double TimE = iT*parameters.getDeltaT();   Velocity_iter = *multiply (unit_velo, *computeVelocity(lattice)); Velocidad_iter = *add(Velocidad_iter, Velocity_iter);   ShearUP_iter = *multiply (unit_Nu, *extractComponent( *computeStrainRateFromStress( lattice, StrengthTapUp), StrengthTapUp, 1)); EsfuerzoUP_iter = *add(ShearUP_iter,EsfuerzoUP_iter);   ShearDown_iter = *multiply (unit_Nu, *extractComponent( *computeStrainRateFromStress ( lattice, StrengthTapDown), StrengthTapDown, 1)); EsfuerzoDown_iter = *add(ShearDown_iter,EsfuerzoDown_iter);   if(counter%parameters.nStep(promedio)==0 && TimE > 2750){ T prom = counter; pcout << prom << endl; EsfuerzoUP_iter = *divide(EsfuerzoUP_iter, prom); std::string str=createFileName("./tmp/TapUp_Re_320_", iT, 6); const char * fileName = str.c_str(); plb_ofstream perfilUvsY1(fileName); perfilUvsY1 << setprecision(16) << EsfuerzoUP_iter << endl; //setToConstant(EsfuerzoUP_iter, Cavity, (T)0.); }   if(counter%parameters.nStep(promedio)==0 && TimE > 2750){ T prom = counter; pcout << prom << endl; EsfuerzoDown_iter = *divide(EsfuerzoDown_iter,prom); std::string str=createFileName("./tmp/TapDown_Re_320_", iT, 6); const char * fileName = str.c_str(); plb_ofstream perfilUvsY2(fileName); perfilUvsY2 << setprecision(16) << EsfuerzoDown_iter << endl; //setToConstant(EsfuerzoDown_iter, Cavity, (T)0.); }   if(counter%parameters.nStep(promedio)==0 && TimE > 2750){ T prom = counter; pcout << prom << endl; Velocidad_iter = *divide(*divide(Velocidad_iter,prom),parameters.getPhysicalU()); std::string str=createFileName("./tmp/ProfilesU_Re_320_", iT, 6); const char * fileName = str.c_str(); plb_ofstream perfilUvsY4(fileName); perfilUvsY4 << setprecision(16)<< Velocidad_iter << endl; //setToConstant(EsfuerzoRight_iter, Cavity, (T)0.); }   UPvsT << setprecision(16) << TimE << " " << *computeVelocityNorm (lattice, profileSectionUP1) << endl;   lattice.collideAndStream(); }   successiveProfilesHalf << setprecision(16) << *computeVelocityNorm(lattice, profileSectionIni) << endl;   successiveStrainUp << setprecision(16) << *multiply (unit_Nu, *extractComponent(*computeStrainRateFromStress(lattice, StrengthTapUp), StrengthTapUp, 1)) << endl;   successiveStrainDown << setprecision(16) << *multiply (unit_Nu, *extractComponent(*computeStrainRateFromStress(lattice, StrengthTapDown), StrengthTapDown, 1)) << endl;   delete boundaryCondition; }

[/code]
Sorry, you do not have permission to post/reply in this forum.