Body force versus pressure driven flow for porous media

Palabos General
#1

Hi All,

I am using the OpenLB for flow simulation in porous media. I started with body force driven flow and my result was good. However, I decided to use pressure driven flow so that I can since my real porous media is complex and not periodic and also, so that I can try the MRT model but I got weird results. I defined an error function using the average energy (i.e lattice.getStatistic) to determine when to stop my simulation. After about 60 hrs of computation on a HPC cluster, my result was not good and the error never converged within my set tolerance as the average energy in the porous media was still dropping.
I had a look at the image file and I found out that at an early time step, the velocity distribution was large in the system. However, as simulation progressed, it reduced and continued to reduce so that the simulation never converged.

Has anybody experienced this with pressure driven flow?

#2

Have you looked at the tutorial on this topic?

http://www.lbmethod.org/palabos/documentation.tutorial/permeability.html

#3

Hi,

I already looked at that tutorial earlier on but I do not use Palabos. I started with OpenLB and it has worked for me so far. However, I understood what it was trying to do but I have some issues with it.
First it says “use bounding box”. Secondly, I do not see a definition of the inflow and outflow boundary conditions that were used, like say the Zou and He which is what I use.
In any case I have below posted a copy of my code for perusal and help if there is any mistake.

Thanks

T tolerance = 1e-6;

T error = 1;

T rho2 = 1, rho1 = 1.1;

T porosity;

int nb = 200, nres = 1; //nres is used because I might have to increase the simulation size by a factor of nres
char n[200][200][200];

T rhox(int iX, LBunits const& converter){
T Lx = converter.getNx()-1;

return ((rho1 - rho2) / Lx)*iX;

}

void memorycreate(BlockStructure3D<T,DESCRIPTOR>& lattice,
LBunitsconst& converter)
{
const int nx = converter.getNx();
const int ny = converter.getNy();
const int nz = converter.getNz();

char ch;

olb_ifstream in("em1a.vox", ios::in | ios::binary);
if(!in){
	cout << " cannot open input file";
}

for (int iz=0; iz< (nz-2)/nres; ++iz) {
	for (int iy=0; iy< (ny-2)/nres; ++iy) {
		for (int ix=0; ix< (nx-2)/nres; ++ix) {
			in.get(ch);
			
			for(int kk = 0; kk < nres; kk++){
				for(int jj = 0; jj < nres; jj++){
					for(int ii = 0; ii < nres; ii++){
						n[nres*iz + kk][nres*iy + jj][nres*ix + ii] = ch;
			
						
					}
				}
			}	
		}
	}
}

in.close();

}

void iniGeometry( BlockStructure3D<T,DESCRIPTOR>& lattice,
LBunits const& converter,
Dynamics<T, DESCRIPTOR>& bulkDynamics,
OnLatticeBoundaryCondition3D<T,DESCRIPTOR>& bc )
{
const int nx = converter.getNx();
const int ny = converter.getNy();
const int nz = converter.getNz();

const T omega = converter.getOmega();


lattice.defineDynamics(0, nx-1, 0, ny-1, 0, nz-1, & bulkDynamics);



bc.addPressureBoundary0N(   0,   0,   1,ny-2,   1,nz-2, omega);
bc.addPressureBoundary0P(nx-1,nx-1,   1,ny-2,   1,nz-2, omega);
bc.addVelocityBoundary1N(   1,nx-2,   0,   0,   1,nz-2, omega);
bc.addVelocityBoundary1P(   1,nx-2,ny-1,ny-1,   1,nz-2, omega);
bc.addVelocityBoundary2N(   1,nx-2,   1,ny-2,   0,   0, omega);
bc.addVelocityBoundary2P(   1,nx-2,   1,ny-2,nz-1,nz-1, omega);

bc.addExternalVelocityEdge0NN(   1,nx-2,   0,   0,   0,   0, omega);
bc.addExternalVelocityEdge0NP(   1,nx-2,   0,   0,nz-1,nz-1, omega);
bc.addExternalVelocityEdge0PN(   1,nx-2,ny-1,ny-1,   0,   0, omega);
bc.addExternalVelocityEdge0PP(   1,nx-2,ny-1,ny-1,nz-1,nz-1, omega);

bc.addExternalVelocityEdge1NN(   0,   0,   1,ny-2,   0,   0, omega);
bc.addExternalVelocityEdge1NP(nx-1,nx-1,   1,ny-2,   0,   0, omega);
bc.addExternalVelocityEdge1PN(   0,   0,   1,ny-2,nz-1,nz-1, omega);
bc.addExternalVelocityEdge1PP(nx-1,nx-1,   1,ny-2,nz-1,nz-1, omega);

bc.addExternalVelocityEdge2NN(   0,   0,   0,   0,   1,nz-2, omega);
bc.addExternalVelocityEdge2NP(   0,   0,ny-1,ny-1,   1,nz-2, omega);
bc.addExternalVelocityEdge2PN(nx-1,nx-1,   0,   0,   1,nz-2, omega);
bc.addExternalVelocityEdge2PP(nx-1,nx-1,ny-1,ny-1,   1,nz-2, omega);

bc.addExternalVelocityCornerNNN(   0,   0,   0, omega);
bc.addExternalVelocityCornerNNP(   0,   0,nz-1, omega);
bc.addExternalVelocityCornerNPN(   0,ny-1,   0, omega);
bc.addExternalVelocityCornerNPP(   0,ny-1,nz-1, omega);
bc.addExternalVelocityCornerPNN(nx-1,   0,   0, omega);
bc.addExternalVelocityCornerPNP(nx-1,   0,nz-1, omega);
bc.addExternalVelocityCornerPPN(nx-1,ny-1,   0, omega);
bc.addExternalVelocityCornerPPP(nx-1,ny-1,nz-1, omega);



for (int iZ=1; iZ< nz-1; ++iZ) {
    for (int iY=1; iY< ny-1; ++iY) {
        for (int iX=1; iX<nx-1; ++iX) {
			
			
			if(n[iZ-1][iY-1][iX-1] == '0')
			{
				//T u = converter.getLatticeU();
				T vel[3] = {T(), T(), T() }; 
				T rho = rho1 - rhox(iX, converter);
				lattice.get(iX,iY,iZ).defineRhoU(rho, vel);
			    lattice.get(iX,iY,iZ).iniEquilibrium(rho, vel);
				
			}
			else if(n[iZ-1][iY-1][iX-1] == '1')
			{
				lattice.defineDynamics(iX,iX,iY,iY,iZ,iZ,
									   &instances::getBounceBack<T, DESCRIPTOR> () );
		    }
		}	
	}
}

T vel[3] = {T(), T(), T() }; 



for(int iZ=0; iZ<nz; ++iZ){
	for(int iX=0; iX<nx; ++iX){			
		T rho = rho1 - rhox(iX, converter);
		lattice.get(iX, 0,iZ).defineRhoU(rho, vel);
		lattice.get(iX, 0, iZ).iniEquilibrium(rho, vel);
	}
}


for(int iZ=0; iZ<nz; ++iZ){
	for(int iX=0; iX<nx; ++iX){
		T rho = rho1 - rhox(iX, converter);							  
		lattice.get(iX, ny-1, iZ).defineRhoU(rho, vel); 
		lattice.get(iX, ny-1, iZ).iniEquilibrium(rho, vel);
	}
}

for(int iY=1; iY<ny-1; ++iY){
	for(int iX=0; iX<nx; ++iX){
		T rho = rho1 - rhox(iX, converter);
		lattice.get(iX, iY, 0).defineRhoU(rho, vel); 
		lattice.get(iX, iY, 0).iniEquilibrium(rho, vel);
	}
}

for(int iY=1; iY<ny-1; ++iY){
	for(int iX=0; iX<nx; ++iX){
		T rho = rho1 - rhox(iX, converter);
		lattice.get(iX, iY, nz-1).defineRhoU(rho, vel); 
		lattice.get(iX, iY, nz-1).iniEquilibrium(rho, vel);
	}
}

for(int iZ=1; iZ<nz-1; ++iZ){
	for(int iY=1; iY<ny-1; ++iY){
		lattice.get(0, iY, iZ).defineRhoU(rho1, vel); 
		lattice.get(0, iY, iZ).iniEquilibrium(rho1, vel);
	}
}

for(int iZ=0; iZ<nz; ++iZ){
	for(int iY=1; iY<ny-1; ++iY){
		lattice.get(nx-1, iY, iZ).defineRhoU(rho2, vel); 
		lattice.get(nx-1, iY, iZ).iniEquilibrium(rho2, vel);
	}
}

lattice.initialize();

}

THIS PART IS THE POROUS MEDIA PERMEABILITY CALCULATION
THIS PART IS THE IMAGE FILE OUTPUTS

int main(int argc, char* argv[]) {
olbInit(&argc, &argv);
singleton::directories().setOutputDir("./tmp/");

LBunits<T> converter(
        (T) 1e-2,  // uMax
		(T) 7.6,  // Re
        (nres*nb+1),        // N
        1.,        // lx
        1.,        // ly
        1.         // lz
);
const T logT     = (T)1/(T)10;
const T imSave   = (T)1/(T)5;
const T vtkSave  = (T)1.;
const T maxT     = (T)100.;




writeLogFile(converter, "CT imaged media");

#ifndef PARALLEL_MODE_MPI // sequential program execution
BlockLattice3D<T, DESCRIPTOR> lattice(converter.getNx(), converter.getNy(), converter.getNz() );
#else // parallel program execution
MultiBlockLattice3D<T, DESCRIPTOR> lattice (
createRegularDataDistribution( converter.getNx(), converter.getNy(), converter.getNz() ) );
#endif

const T omega = converter.getOmega();
BGKdynamics<T, DESCRIPTOR> bulkDynamics (
                  omega,
                  instances::getBulkMomenta<T, DESCRIPTOR>()
);

OnLatticeBoundaryCondition3D<T, DESCRIPTOR>*
boundaryCondition = createZouHeBoundaryCondition3D(lattice);



memorycreate(lattice, converter);
iniGeometry(lattice, converter, bulkDynamics, *boundaryCondition);

int iT = 0;

T absvel1 = 0, absvel2 = 1;

for(iT = 0; iT<100000; ++iT){
//do	{
		if (iT%converter.nStep(logT)==0) {
			cout << "step " << iT
				<< "; t=" << iT*converter.getDeltaT()
				<< "; av energy="
				<< lattice.getStatistics().getAverageEnergy()
				<< "; av rho="
				<< lattice.getStatistics().getAverageRho()
				<< " error= " << error << " absvel1= "<< absvel1 << " absvel2= " << absvel2 <<endl << flush;
		}
	
		if (iT%converter.nStep(imSave)==0 && iT>0) {
			cout << "Saving Gif ..." << endl << flush;
			writeGifs(lattice, converter, iT);
		}
	
		if (iT%converter.nStep(vtkSave)==0 && iT>0) {
			cout << "Saving VTK file ..." << endl << flush;
			writeVTK(lattice, converter, iT);
		}
	
		lattice.collideAndStream();
	
		absvel2 = lattice.getStatistics().getAverageEnergy();
		error = fabs(1 - (absvel1/absvel2));
		absvel1 = absvel2;
	//	++iT;
	
	if(error < tolerance)
		break;
}
//while (error > tolerance);

	averagevelocity(lattice, converter, iT);	
	cout << "final time step= " << iT <<" error= "<<error<<" max veloctiy= "<<lattice.getStatistics().getMaxU() << endl;

	delete boundaryCondition;

}

#4

sorry the message were not for you.

ciao
Andrea