Hi , I am trying to simulate wetting of a 2D porous slice. I need a constant density at the inlet which would simulate a constant head of wetting fluid, that should subsequently fill the vapour space. My code appears to work but with a non constant head of wetting fluid i.e. at iteration 0, the inlet has a density of 200, I would like to preserve this value of density at the inlet for all future iterations, at the moment the density at the inlet gradually decreases as it mixes with the vapour fluid. I tried to change the boundary conditions and use SetPressureCondition for the inlet but this appears to be incompatible with the ShanChen 2D descriptor. I post my code below for anyone to advise how this can be done as I seem to be going around in circles??
Cheers,
Ruth
/* This file is part of the Palabos library.
* Copyright (C) 2009 Jonas Latt
* E-mail contact: jonas@lbmethod.org
* The most recent release of Palabos can be downloaded at
* <http://www.lbmethod.org/palabos/>
*
* The library Palabos is free software: you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3 of the
* License, or (at your option) any later version.
*
* The library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include "palabos2D.h"
#include "palabos2D.hh"
#include <cstdlib>
#include <iostream>
#include <iomanip>
using namespace plb;
using namespace std;
typedef double T;
#define DESCRIPTOR descriptors::ShanChenD2Q9Descriptor
void writeVTK(MultiBlockLattice2D<T,DESCRIPTOR>& lattice,
plint iter)
{
//T dx = parameters.getDeltaX();
// T dt = parameters.getDeltaT();
VtkImageOutput2D<T> vtkOut(createFileName("vtk", iter, 6), 1.0);
vtkOut.writeData<float>(*computeDensity(lattice), "Density", 1.0);
}
int main(int argc, char *argv[])
{
plbInit(&argc, &argv);
global::timer("simTime").start();
global::directories().setOutputDir("./tmp/");
OnLatticeBoundaryCondition2D<T,DESCRIPTOR>* boundaryCondition = createEquilibriumBoundaryCondition2D<T,DESCRIPTOR>();
// For the choice of the parameters G, rho0, and psi0, we refer to the book
// Michael C. Sukop and Daniel T. Thorne (2006),
// Lattice Boltzmann Modeling; an Introduction for Geoscientists and Engineers.
// Springer-Verlag Berlin/Heidelberg.
T tIni = global::timer("simTime").stop();
global::timer("simTime").start();
const T omega = 1.0;
const int nx = 272;
const int ny = 272;
//const int nz = 32;
const T G = -92.0; //setting the force to zero
const int vapour_inlet=10;
const int maxIter = 100001;//100001;
const int saveIter = 10;
const int statIter = 10;
plint evalTime =100;
const T liquid_fluid = 200.0;
const T psi0 = 4.0;
const T vel=0.05;
const T virt_density=liquid_fluid+5.0;
const T vapour_fluid=40.0;
const int inlet=10;
pcout << "Reading the geometry file." << endl;
std::string fNameIn = "SOslice.dat";
MultiScalarField2D<T> geometry(nx,ny);
plb_ifstream geometryFile(fNameIn.c_str());
if (!geometryFile.is_open()) {
pcout << "Error: could not open geometry file " << fNameIn << endl;
return -1;
}
geometryFile >> geometry;
//reset geometry file everything is a pore
//geometry.reset();
MultiBlockLattice2D<T, DESCRIPTOR> lattice ( nx,ny, new ExternalMomentRegularizedBGKdynamics <T, DESCRIPTOR>(omega) );
lattice.periodicity().toggle(0, false);
lattice.periodicity().toggle(1, false);
// Add the data processor which implements the Shan/Chen interaction potential.
plint processorLevel = 1;
integrateProcessingFunctional (
new ShanChenSingleComponentProcessor2D<T,DESCRIPTOR> (
G, new interparticlePotential::PsiShanChen94<T>(psi0,liquid_fluid) ),
lattice.getBoundingBox(),
lattice,
processorLevel );
initializeAtEquilibrium (
lattice, Box2D(0, inlet, 0, ny-1),liquid_fluid,
Array<T,2>(0.0,0.0) );
initializeAtEquilibrium (
lattice, Box2D(inlet+1, nx-1, 0,ny-1),vapour_fluid,
Array<T,2>(0.0, 0.0) );
//OnLatticeBoundaryCondition2D<T,DESCRIPTOR>* boundaryCondition = createEquilibriumBoundaryCondition2D<T,DESCRIPTOR>();
//y boundaries
setBoundaryVelocity(lattice, Box2D(0, inlet, ny-1, ny-1), Array<T,2>(0.0,0.) );
initializeAtEquilibrium (
lattice, Box2D(0, inlet, ny-1, ny-2),liquid_fluid,
Array<T,2>(0.0,0.0) );
setBoundaryVelocity(lattice, Box2D(inlet+1, nx-1, ny-1, ny-1), Array<T,2>(0.0,0.) );
initializeAtEquilibrium (
lattice, Box2D(inlet+1, nx-1, ny-1, ny-2),vapour_fluid,
Array<T,2>(0.0,0.0) );
setBoundaryVelocity(lattice, Box2D(0, inlet, 0, 0), Array<T,2>(0.0,0.) );
initializeAtEquilibrium (
lattice, Box2D(0, inlet, 0, 0),liquid_fluid,
Array<T,2>(0.0, 0.0) );
setBoundaryVelocity(lattice, Box2D(inlet+1, nx-1, 0, 0), Array<T,2>(0.0,0.) );
initializeAtEquilibrium (
lattice, Box2D(inlet+1, nx-1, 0, 0),vapour_fluid,
Array<T,2>(0.0, 0.0) );
//x boundaries
setBoundaryVelocity(lattice, Box2D(nx-1, nx-1, 0, ny-1), Array<T,2>(vel,0.) );
initializeAtEquilibrium (
lattice, Box2D(nx-1, nx-1, 1, ny-2),vapour_fluid,
Array<T,2>(vel,0.0) );//0.00005
setBoundaryVelocity(lattice, Box2D(0, 0, 0,ny-1), Array<T,2>(vel,0.) );
initializeAtEquilibrium (
lattice, Box2D(0, 0, 1,ny-2),liquid_fluid,
Array<T,2>(vel, 0.0) );
// Where "geometry" evaluates to 1, use bounce-back.//0.5 adhesive force
defineDynamics(lattice, geometry, new BounceBack<T, DESCRIPTOR>(virt_density), true);
defineDynamics(lattice, Box2D(Box2D(0, nx-1, 0, 0)), new BounceBack<T, DESCRIPTOR>(virt_density));
defineDynamics(lattice, Box2D(Box2D(0, nx-1, ny-1, ny-1)), new BounceBack<T, DESCRIPTOR>(virt_density));
lattice.initialize();
pcout << "Starting simulation" << endl;
// The value-tracer is used to stop the simulation once is has converged.
// 1st parameter:velocity
// 2nd parameter:size
// 3rd parameters:threshold
// 1st and second parameters ae used for the length of the time average (size/velocity)
util::ValueTracer<T> converge(1.0,1000.0,1.0e-4);
int iT=1;
for (int iT; iT<maxIter; ++iT) {
//initializeAtEquilibrium (
//lattice, Box2D(0, 10, 0, ny-1),200.0,
//Array<T,2>(0.05,0.0) );
//lattice, Box2D(0, 20, 0, ny-1),200.0 );
if (iT%statIter==0) {
auto_ptr<MultiScalarField2D<T> > rho( computeDensity(lattice) );
converge.takeValue(getStoredAverageEnergy(lattice),true);
pcout << "iteration " << iT << endl;
pcout << "Average density: " << computeAverage(*computeDensity(lattice, Box2D(0, nx-1, 0, ny-1 ))) << endl;
pcout << "Minimum density: " << computeMin(*rho) << endl;
pcout << "Maximum density: " << computeMax(*rho) << endl;
}
if (iT %saveIter==0) {
ImageWriter<T>("leeloo.map").writeGif (
createFileName("rho", iT, 6),*computeDensity(lattice, Box2D(0, nx-1, 0, ny-1 )), vapour_fluid, virt_density ); //writeVTK(lattice,iT);
char path[200]="density";//createFileName("density", iT, 6)
char iter[100]="";
sprintf(iter, "%d", iT);
strcat(path, iter);
strcat(path, ".dat");
pcout << path << endl;
plb_ofstream ofile(path);
ofile << setprecision(10) << *computeDensity(lattice) << endl;
}
lattice.collideAndStream();
converge.takeValue(getStoredAverageEnergy(lattice),true);
if (converge.hasConverged())
{
// break;
}
}
//convert(convert -delay 5 rho*.gif animation.gif);
T tEnd = global::timer("simTime").stop();
T totalTime = tEnd-tIni;
pcout << "number of processors: " << global::mpi().getSize() << endl;
pcout << "simulation time: " << totalTime << endl;
}