# How to impose flux through boundaries.

I’m trying to simulate an imbibition-drainage cycle of a tank in which some obstacles are included inside the tank. In order to perform this sumlation I have modified the example code damBreak3d. I have managed to “create” and “remove” mass by changing the flag of some nodes (fluid flag or empty flag). So , I’m able to reproduce the imbibition and drainage cycle, in fact, you can take a look at the code to better understand what I’m talking about. However, I would like to solve the problem in a more neat way. Rather than introducing a soure of mass, I would prefer imposing a flux through one of the boundaries. Is that possible? I have seen examples in which flow is injected (i.e. rectangularChannel3D.cpp), nevertheless, no flags are involved. I don’t know how can I melt both problems, I mean, use flag decriptions and impose the boundary conditions I want. Thank you!

[code=“cpp”]

#include “palabos3D.h”
#include “palabos3D.hh”

using namespace plb;
using namespace std;

#define DESCRIPTOR descriptors::ForcedD3Q19Descriptor
#define NMAX 40

typedef double T;

// Smagorinsky constant for LES model.
const T cSmago = 0.14;

// Physical dimensions of the system (in meters).
const T lx = 1.0;
const T ly = 1.0;
const T lz = 1.0;
const T rhoEmpty = T(1);
Array<T,3> forceOrientation(T(),T(),(T)1);

plint writeImagesIter = 100;
plint getStatisticsIter = 20;

plint nmin;
plint maxIter;
plint N;
T delta_t, delta_x;
Array<T,3> externalForce;
T nuPhys, nuLB, tau, omega, Bo, surfaceTensionLB, contactAngle;

std::string outDir;
// plint beginWaterReservoir, waterReservoirHeight;
// plint waterLevelOne, waterLevelTwo, waterLevelThree, waterLevelFour;

void setupParameters() {
delta_x = lz / N;
nx = util::roundToInt(lx / delta_x);
ny = util::roundToInt(ly / delta_x);
nz = util::roundToInt(lz / delta_x);
// Gravity in lattice units.
T gLB = 9.8 * delta_t * delta_t/delta_x;
externalForce = Array<T,3>(0., 0., -gLB);
tau = (nuPhysDESCRIPTOR::invCs2delta_t)/(delta_x*delta_x) + 0.5;
omega = 1./tau;
nuLB = (tau-0.5)*DESCRIPTOR::cs2; // Viscosity in lattice units.
nmin = 100;
surfaceTensionLB = rhoEmpty * gLB * N * N / Bo;

}

int initialFluidFlags(plint iX, plint iY, plint iZ) {
// Place an obstacle on the left end, which is hit by the fluid.
Array<T,3> pos(iX, iY, iZ);

``````Array<T,3> center_part1(nx/2, ny/2.0, radius_particleLB);
T r_part1 = norm(pos-center_part1);

T r_part3 = norm(pos-center_part3);

T r_part2 = norm(pos-center_part2);

return twoPhaseFlag::wall;
}

else if (r_part3 <= radius_particleLB) {
return twoPhaseFlag::wall;
}

else if (r_part2 <= radius_particleLB) {
return twoPhaseFlag::wall;
}

else if ( iZ <= 0.) {
return twoPhaseFlag::fluid;
}

else {
return twoPhaseFlag::empty;
}
``````

}

int secondFluidFlags(plint iX, plint iY, plint iZ) {

``````if ((iX <= nx/30) && (iX >= 1) && (iY <= ny-1) && (iY >= 1) && (iZ <= nz/4) && (iZ >= 1) ) {
return twoPhaseFlag::fluid;
}
``````

}

int thirdFluidFlags(plint iX, plint iY, plint iZ) {

``````if ((iX <= nx-1) && (iX >= nx - nx/30) && (iY <= ny-1) && (iY >= 1) && (iZ <= nz/4) && (iZ >= 1) ) {
return twoPhaseFlag::fluid;
}
``````

}

int forthFluidFlags(plint iX, plint iY, plint iZ) {

``````if ((iX <= nx/30) && (iX >= 1) && (iY <= ny-1) && (iY >= 1) && (iZ <= nz/4) && (iZ >= 1) ) {
return twoPhaseFlag::empty;
}
``````

}

int fifthFluidFlags(plint iX, plint iY, plint iZ) {

``````if ((iX <= nx-1) && (iX >= nx - nx/30) && (iY <= ny-1) && (iY >= 1) && (iZ <= nz/4) && (iZ >= 1) ) {
return twoPhaseFlag::empty;
}
``````

}

void writeResults(MultiBlockLattice3D<T,DESCRIPTOR>& lattice, MultiScalarField3D& volumeFraction, plint iT)
{

``````Box3D slice(0, nx-1, 0, ny-1, 0, nz-1);
ImageWriter<T> imageWriter("leeloo");
imageWriter.writeScaledPpm(createFileName("u", iT, 6),
*computeVelocityNorm(lattice, slice));

imageWriter.writeScaledPpm(createFileName("rho", iT, 6),
*computeDensity(lattice, slice));

imageWriter.writeScaledPpm(createFileName("volumeFraction", iT, 6), *extractSubDomain(volumeFraction, slice));

// Use a marching-cube algorithm to reconstruct the free surface and write an STL file.
std::vector<T> isoLevels;
isoLevels.push_back((T) 0.5);

typedef TriangleSet<T>::Triangle Triangle;

std::vector<Triangle> triangles;

isoSurfaceMarchingCube(triangles, volumeFraction, isoLevels, volumeFraction.getBoundingBox());

TriangleSet<T>(triangles).writeBinarySTL(createFileName(outDir+"/interface", iT, 6)+".stl");
``````

// pcout << "aqui 5: " << std::endl;
VtkImageOutput3D vtkOut(createFileName(“volumeFraction”, iT, 6), 1.);
vtkOut.writeData(volumeFraction, “vf”, 1.);
vtkOut.writeData(*computeVelocityNorm(lattice, slice),
“velocityNorm”, delta_x/delta_t );
vtkOut.writeData<3,float>(*computeVelocity(lattice, slice), “velocity”, delta_x/delta_t);

``````Array<T,3> center_part1(nx/2, ny/2.0, radius_particleLB);

TriangleSet<T> triangleSet;

// The next few lines of code are typical. They transform the surface geometry of the
//   tube to more efficient data structures that are internally used by palabos.
//   The TriangleBoundary3D structure will be later used to assign proper boundary conditions.
DEFscaledMesh<T> defMesh(triangleSet);
defMesh.getMesh().inflate();
TriangleBoundary3D<T> boundary(defMesh);
boundary.getMesh().writeAsciiSTL(outDir+"/sphere1.stl");

TriangleSet<T> triangleSet2;

// The next few lines of code are typical. They transform the surface geometry of the
//   tube to more efficient data structures that are internally used by palabos.
//   The TriangleBoundary3D structure will be later used to assign proper boundary conditions.
DEFscaledMesh<T> defMesh2(triangleSet2);
defMesh2.getMesh().inflate();
TriangleBoundary3D<T> boundary2(defMesh2);
boundary2.getMesh().writeAsciiSTL(outDir+"/sphere2.stl");
``````

//

``````Array<T,3> center_part3(nx/2, ny/2.0+2*radius_particleLB, radius_particleLB);
TriangleSet<T> triangleSet3;

// The next few lines of code are typical. They transform the surface geometry of the
//   tube to more efficient data structures that are internally used by palabos.
//   The TriangleBoundary3D structure will be later used to assign proper boundary conditions.
DEFscaledMesh<T> defMesh3(triangleSet3);
defMesh3.getMesh().inflate();
TriangleBoundary3D<T> boundary3(defMesh3);
boundary3.getMesh().writeAsciiSTL(outDir+"/sphere3.stl");
``````

}

void writeStatistics(TwoPhaseFields3D<T,DESCRIPTOR>& fields) {
pcout << " ------------------------------- " << endl;
T averageMass = freeSurfaceAverageMass<T,DESCRIPTOR>(fields.twoPhaseArgs, fields.lattice.getBoundingBox());
pcout << "Average Mass: " << averageMass << endl;
T averageDensity = freeSurfaceAverageDensity<T,DESCRIPTOR>(fields.twoPhaseArgs, fields.lattice.getBoundingBox());
pcout << "Average Density: " << setprecision(12) << averageDensity << endl;

``````T averageVolumeFraction = freeSurfaceAverageVolumeFraction<T,DESCRIPTOR>(fields.twoPhaseArgs, fields.lattice.getBoundingBox());
pcout << "Average Volume-Fraction: " << setprecision(12) << averageVolumeFraction  << endl;

pcout << " -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- " << endl;
``````

}

int main(int argc, char **argv)
{
plbInit(&argc, &argv);
global::directories().setInputDir("./");

``````if (global::argc() != 8) {
pcout << "Error missing some input parameter\n";
}

try {
global::directories().setOutputDir(outDir+"/");

}
catch(PlbIOException& except) {
pcout << except.what() << std::endl;
pcout << "The parameters for this program are :\n";
pcout << "1. Output directory name.\n";
pcout << "2. kinematic viscosity in physical Units (m^2/s) .\n";
pcout << "3. Bond number (Bo = rho * g * L^2 / gamma).\n";
pcout << "4. Contact angle (in degrees).\n";
pcout << "5. number of lattice nodes for lz .\n";
pcout << "6. delta_t .\n";
pcout << "7. maxIter .\n";
pcout << "Reasonable parameters on a desktop computer are: " << (std::string)global::argv(0) << " tmp 1.e-5 200 80.0 40 1.e-3 80000\n";
pcout << "Reasonable parameters on a parallel machine are: " << (std::string)global::argv(0) << " tmp 1.e-6 200 80.0 100 1.e-4 80000\n";
exit (EXIT_FAILURE);
}

setupParameters();

pcout << "delta_t= " << delta_t << endl;
pcout << "delta_x= " << delta_x << endl;
pcout << "delta_t*delta_t/delta_x= " << delta_t*delta_t/delta_x << endl;
pcout << "externalForce= " << externalForce[2] << endl;
pcout << "relaxation time= " << tau << endl;
pcout << "omega= " << omega << endl;
pcout << "kinematic viscosity physical units = " << nuPhys << endl;
pcout << "kinematic viscosity lattice units= " << nuLB << endl;

global::timer("initialization").start();

Dynamics<T,DESCRIPTOR>* dynamics
= new SmagorinskyBGKdynamics<T,DESCRIPTOR>(omega, cSmago);

SparseBlockStructure3D blockStructure(createRegularDistribution3D(nx, ny, nz));

// If surfaceTensionLB is 0, then the surface tension algorithm is deactivated.
// If contactAngle is less than 0, then the contact angle algorithm is deactivated.
TwoPhaseFields3D<T,DESCRIPTOR> fields( blockStructure, dynamics->clone(), rhoEmpty,
surfaceTensionLB, contactAngle, externalForce );

// Set all outer-wall cells to "wall" (here, bulk-cells are also set to "wall", but it
// doesn't matter, because they are overwritten on the next line).

setToConstant(fields.flag, fields.flag.getBoundingBox(), (int)twoPhaseFlag::wall);

// In the bulk (all except outer wall layer), initialize the flags as specified by
// the function "initialFluidFlags".

setToFunction(fields.flag, fields.flag.getBoundingBox().enlarge(-1), initialFluidFlags);

fields.defaultInitialize();
``````

/*
OnLatticeBoundaryCondition3D<T,DESCRIPTOR>* boundaryCondition
= createLocalBoundaryCondition3D<T,DESCRIPTOR>();

``````SetUpBC(lattice, parameters, *boundaryCondition);*/

pcout << "Time spent for setting up lattices: "
<< global::timer("initialization").stop() << endl;
T lastIterationTime = T();

Box3D slice2(1, nx/30, 1, ny, 1, nz/4);
Box3D slice3(nx-nx/30, nx-1, 1, ny, 1, nz/4);

for (plint iT = 0; iT <= maxIter; ++iT) {
global::timer("iteration").restart();

T sum_of_mass_matrix = T();
T lost_mass = T();
if (iT % getStatisticsIter==0) {
pcout << endl;
pcout << "ITERATION = " << iT << endl;
pcout << "Time of last iteration is " << lastIterationTime << " seconds" << endl;
writeStatistics(fields);
sum_of_mass_matrix = fields.lattice.getInternalStatistics().getSum(0);
pcout << "Sum of mass matrix: " << sum_of_mass_matrix << std::endl;
lost_mass = fields.lattice.getInternalStatistics().getSum(1);
pcout << "Lost mass: " << lost_mass << std::endl;
pcout << "Total mass: " << sum_of_mass_matrix + lost_mass << std::endl;
pcout << "Interface cells: " << fields.lattice.getInternalStatistics().getIntSum(0) << std::endl;
}

if (iT % writeImagesIter == 0) {
global::timer("images").start();
writeResults(fields.lattice, fields.volumeFraction, iT);
pcout << "Total time spent for writing images: "
<< global::timer("images").stop() << endl;

if (iT <= 50000) {
setToFunction(fields.flag, slice2, secondFluidFlags);
setToFunction(fields.flag, slice3, thirdFluidFlags);

fields.partiallyDefaultInitialize();
}

if (iT > 50000) {
setToFunction(fields.flag, slice2, forthFluidFlags);
setToFunction(fields.flag, slice3, fifthFluidFlags);

fields.partiallyDefaultInitialize();
}

}

// This includes the collision-streaming cycle, plus all free-surface operations.

fields.lattice.executeInternalProcessors();
fields.lattice.evaluateStatistics();
fields.lattice.incrementTime();

lastIterationTime = global::timer("iteration").stop();
}
``````

}

I think flow injection is possible by adding the following processing functionals. I also noticed that the PouringLiquid3D seems to inject fluid only if the coressponding cells are initialized as fluid cells and you need the DESCRIPTOR descriptors::ForcedShanChenD3Q19Descriptor that the functionals work.

[code=“cpp”]
// For Inlets
integrateProcessingFunctional(new PouringLiquid3D<T,DESCRIPTOR>(dynamics->clone(),Array<T,3>(parameters.getLatticeU(),0.,0.)), inlet, fields.twoPhaseArgs, 0);
// For Outlets
integrateProcessingFunctional(new RemoveMass3D<T,DESCRIPTOR>(), outlet, fields.twoPhaseArgs, 0);