# Problem with Passive Scalar

Hello fellow Palabos users:-)

I am trying to simulate the flow in a microchannel containing obstacles. The latter are supposed to act as vortex promoters in order to enhance mixing of two streams of passive scalars introduced at the inlet and obeying the advection diffusion equation.
Basically, the simulation runs well except for one decisive problem with the passive scalar (PS). When I compute the maximum and minimum PS with

[code=“cpp”]

``````
the minimum value is always zero instead of 1 and worse the maximum value rises to ~2.8 and then stays more or less constant. Analysis of the PS field with ParaView shows that the PS is greater than two at the inlet area and right in front of the obstacles.
If someone could provide any pointer to the cause of the problem I would greatly appreciate it since this is a master thesis and I'm slowly running out of time due to this issue:-( Below you can find the code I'm using.

I use following dynamics and descriptors:
[code="cpp"]
#define NSDYNAMICS GuoExternalForceBGKdynamics
#define NSDESCRIPTOR descriptors::ForcedD3Q19Descriptor

``````

The passive scalar is initialized with following data processor which will be used to set PS to 2 at the top of the inlet and to 1 at the bottom and to 1.5 (average) everywhere else:

[code=“cpp”]
template<typename T, template class nsDescriptor,
struct psFieldInitProcessor3D :
{
: parameters(parameters_)
{ }
{
// the respective Multiblock must know its location
// with respect to the global coordinates

``````for (plint iX=domain.x0; iX<=domain.x1; ++iX) {
for (plint iZ=domain.z0; iZ<=domain.z1; ++iZ) {
for (plint iY=domain.y0; iY<=domain.y1; ++iY) {

plint absoluteY = absoluteOffset.y  +  iY;
plint absoluteX = absoluteOffset.x  +  iX;

// The initial condition is set with respect
// to the global coordinates
T ps;
if(absoluteX == 0){
// in the inlet 50/50
if (absoluteY > (parameters.getNy()-1)/2) ps = parameters.getMaxPS();
else ps = parameters.getMinPS();
}
// else average
else ps = parameters.getMeanPS();
}
}
}
}
// this is needed for any class in Palabos
{
}

virtual void getTypeOfModification(std::vector<modif::ModifT>& modified) const {
modified[0] = modif::staticVariables;
}

virtual BlockDomain::DomainT appliesTo() const {
return BlockDomain::bulkAndEnvelope;
}

private :
``````

};

``````
The channel setup then looks as follows:
[code="cpp"]
void    channelSetup( MultiBlockLattice3D<T,NSDESCRIPTOR>& nsLattice,
OnLatticeBoundaryCondition3D<T,NSDESCRIPTOR>& nsBoundaryCondition,
const plint sections)
{
const T l =  parameters.getResolution();
const plint nx = parameters.getNx();
const plint ny = parameters.getNy();
const plint nz = parameters.getNz();

Box3D inletFluid	(0,	0,	1,   ny-2, 1, nz-2);
Box3D outletFluid	(nx-1,	nx-1,	1,   ny-2, 1, nz-2);
Box3D topFluid	(0,	nx-1,	ny-1,ny-1, 0, nz-1);
Box3D bottomFluid	(0,	nx-1,	0,   0,    0, nz-1);
Box3D leftFluid	(0,	nx-1,	0,   ny-1, 0, 0);
Box3D rightFluid	(0,	nx-1,	0,   ny-1, nz-1,nz-1);

Box3D inletThermal	(0,	0,	0,   ny-1, 0, nz-1);
Box3D outletThermal	(nx-1,	nx-1,	0,   ny-1, 0, nz-1);
Box3D topThermal	(1,	nx-2,	ny-1,ny-1, 1,nz-2);
Box3D bottomThermal	(1,	nx-2,	0,   0, 1,nz-2);
Box3D leftThermal	(1,	nx-2,	0,   ny-1, 0,0);
Box3D rightThermal	(1,	nx-2,	0,   ny-1, nz-1, nz-1);

inletFluid, nsLattice);

outletFluid, nsLattice);
setBoundaryDensity(nsLattice, outletFluid, (T)1.);

bottomFluid,nsLattice);
topFluid, nsLattice);
leftFluid, nsLattice);
rightFluid, nsLattice);

setBoundaryVelocity (
nsLattice, nsLattice.getBoundingBox(),
PoiseuilleVelocity<T>(parameters) );

setBoundaryDensity (
nsLattice, outletFluid,
ConstantDensity<T>(1.) );

initializeAtEquilibrium (
nsLattice, nsLattice.getBoundingBox(),
PoiseuilleVelocityAndDensity<T>(parameters) );

applyProcessingFunctional (

plint processorLevel = 1;
integrateProcessingFunctional (

integrateProcessingFunctional (

integrateProcessingFunctional (

integrateProcessingFunctional (

//Obstacles
for(plint iM = 0; iM < sections ; ++iM){
defineDynamics(nsLattice, nsLattice.getBoundingBox(),
new PeriodicGroovedDomain3D<T>(2*l + iM*12*l,parameters),
new plb::BounceBack<T,NSDESCRIPTOR>);

new PeriodicGroovedDomain3D<T>(2*l + iM*12*l,parameters),

}

nsLattice.initialize();
}

``````

The class “PeriodicGroovedDomain3D” is a bool mask which gives true in case of obstacle and false otherwise. The outlet boundary condition for the PS is of Neumann type and implemented by copying the previous populations to the outlet nodes:

[code=“cpp”]
{
const plint nx = parameters.getNx();
const plint ny = parameters.getNy();
const plint nz = parameters.getNz();

``````// On the right boundary, we copy the unknown populations from previous locations
integrateProcessingFunctional(
``````

}

``````

The relevant parts of main() are:

[code="cpp"]
///Generate lattices
T nsOmega = parameters.getFluidOmega();
const plint nx = parameters.getNx();
const plint ny = parameters.getNy();
const plint nz = parameters.getNz();
pcout << "Generating lattices of size\n";
pcout << "nx: " << nx << "\n";
pcout << "ny: " << ny << "\n";
pcout << "nz: " << nz << "\n";
MultiBlockLattice3D<T, NSDESCRIPTOR> nsLattice (
nx,ny,nz,new NSDYNAMICS<T, NSDESCRIPTOR>(nsOmega) );
nsLattice.periodicity().toggleAll(false);

///Channel Setup
pcout << "Setting up the channel..." << endl;
OnLatticeBoundaryCondition3D<T,NSDESCRIPTOR>*
nsBoundaryCondition = createLocalBoundaryCondition3D<T,NSDESCRIPTOR>();

Array<T,NSDESCRIPTOR<T>::d> forceOrientation(T(),(T)1,T());
plint processorLevel = 1;
integrateProcessingFunctional (
parameters.getLatticeGravity(), parameters.getMeanPS(),
parameters.getDeltaPS(),forceOrientation ),
nsLattice.getBoundingBox(),

pcout << "Finished setting up channel!" << endl;

``````

where getLatticeGravity() is set to zero.

Regards,

Paul

Hello,

I suspect that the problem (“unboundedness” of passive scalar, see above post) is due to the boundary conditions for the thermal lattice. I’m using prescribed values at the inlet (data processor, see above post), adiabatic walls (FlatAdiabaticBoundaryFunctional3D) and a zero gradient condition at the outlet (CopyUnknownPopulationsFunctional3D).
Is perhaps something wrong with my definition of the walls/inlet/outlet Box3D objects of the fluid and/or thermal lattice (see above post)?
Does anybody know an other way of implementing adiabatic walls and a Neumann boundary at the outlet for the thermal lattice? In this context I also wonder where the difference between FlatAdiabaticBoundaryFunctional3D and CopyUnknownPopulationsFunctional3D is since they both seem to copy previous values in normal direction.

Appreciate any hints (no matter how small:-)

Regards,

Paul

Dear Paul,

The difference is that, while CopyUnknownPopulationsFunctional3D extrapolates unknown populations (after which the boundary condition is completed), FlatAdiabaticBoundaryFunctional3D extrapolates the value of the scalar, which then should be used to implement a Dirichlet (“fixed Temperature” / “fixed Scalar”) boundary condition.

This directly leads us to what seems to be the issue with your code. It is not sufficient to add the FlatAdiabaticBoundaryFunctional3D data processor to a subdomain of the lattice. You also must implement a Dirichlet boundary condition on the same subdomain with help of your adBoundaryCondition object.

Cheers,
Jonas

Hello Jonas,

thanks for your response, I really appreciate it.
So I issued the “adBoundaryCondition.addTemperatureBoundary” command for the walls right after the one for the inlet (using the same boxes as for the adiabatic functional), the rest of the channel setup was kept as in the first post.
The simulation crashes very fast with the passive scalar increasing exponentially, so I guess I misunderstood you somehow.
What about my outlet boundary condition for the adLattice? Currently I’m calling
[code=“cpp”]
integrateProcessingFunctional(
``````