Consistency law at Lattice Boltzmann using Palabos

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


[code=“cpp”]

[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 ForcedD3Q19Descriptor
#define DYNAMICS GuoExternalForceBGKdynamics

T channelForce(IncomprFlowParam const& parameters)
{
return (T) 6e-6;
}

void channelSetup(MultiBlockLattice3D<T,DESCRIPTOR>& lattice,
IncomprFlowParam 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::ExternalField::forceBeginsAt,force);

lattice.initialize();

}

int main(int argc, char*argv[]){
plbInit(&argc,&argv);

global::directories().setOutputDir("./tmp/");

IncomprFlowParam 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; iTparameters.getDeltaT()<1.01maxT; ++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;
}