Problem with Passive Scalar

Palabos General
#1

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”]
computeMax(*computeDensity(adLattice),adLattice.getBoundingBox())
computeMin(*computeDensity(adLattice),adLattice.getBoundingBox())


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 ADYNAMICS AdvectionDiffusionBGKdynamics
#define NSDESCRIPTOR descriptors::ForcedD3Q19Descriptor
#define ADESCRIPTOR descriptors::AdvectionDiffusionD3Q7Descriptor

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,
template class adDescriptor>
struct psFieldInitProcessor3D :
public BoxProcessingFunctional3D_L<T,adDescriptor>
{
psFieldInitProcessor3D(SuperFlowParam<T,nsDescriptor,adDescriptor> parameters_)
: parameters(parameters_)
{ }
virtual void process(Box3D domain, BlockLattice3D<T,adDescriptor>& adLattice)
{
// the respective Multiblock must know its location
// with respect to the global coordinates
Dot3D absoluteOffset = adLattice.getLocation();

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();
	    Array<T,adDescriptor<T>::d> jEq(0.0, 0.0, 0.0);
	    adLattice.get(iX,iY,iZ).defineDensity(ps);
	    iniCellAtEquilibrium(adLattice.get(iX,iY,iZ), ps, jEq);
	}
    }
}
}
// this is needed for any class in Palabos
virtual psFieldInitProcessor3D<T,nsDescriptor,adDescriptor>* clone() const
{
return new psFieldInitProcessor3D<T,nsDescriptor,adDescriptor>(*this);
}


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

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

private :
SuperFlowParam<T,NSDESCRIPTOR,ADESCRIPTOR> parameters;

};


The channel setup then looks as follows:
[code="cpp"]
void    channelSetup( MultiBlockLattice3D<T,NSDESCRIPTOR>& nsLattice,
	MultiBlockLattice3D<T,ADESCRIPTOR>& adLattice,	
	SuperFlowParam<T,NSDESCRIPTOR,ADESCRIPTOR> const& parameters,
	OnLatticeBoundaryCondition3D<T,NSDESCRIPTOR>& nsBoundaryCondition,
	OnLatticeAdvectionDiffusionBoundaryCondition3D<T,ADESCRIPTOR>& adBoundaryCondition,
	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);
    

    nsBoundaryCondition.addVelocityBoundary0N (
	    inletFluid, nsLattice); 

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


    nsBoundaryCondition.addVelocityBoundary1N (
	    bottomFluid,nsLattice);
    nsBoundaryCondition.addVelocityBoundary1P (
	    topFluid, nsLattice);
    nsBoundaryCondition.addVelocityBoundary2N (
	    leftFluid, nsLattice);
    nsBoundaryCondition.addVelocityBoundary2P (
	    rightFluid, nsLattice);  

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

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


    adBoundaryCondition.addTemperatureBoundary0N(inletThermal,  adLattice); 

   
    initializeAtEquilibrium (
	    nsLattice, nsLattice.getBoundingBox(),
	    PoiseuilleVelocityAndDensity<T>(parameters) );
  
     applyProcessingFunctional (
	    new psFieldInitProcessor3D<T,NSDESCRIPTOR,ADESCRIPTOR>(parameters),
	    adLattice.getBoundingBox(),  adLattice );
    	
    /// define adiabatic walls 
    plint processorLevel = 1;
    integrateProcessingFunctional (
	    new FlatAdiabaticBoundaryFunctional3D<T,ADESCRIPTOR,1,1>,
	    topThermal, adLattice, processorLevel );

    integrateProcessingFunctional (
	    new FlatAdiabaticBoundaryFunctional3D<T,ADESCRIPTOR,1,-1>,
	    bottomThermal, adLattice, processorLevel );
   
 
    integrateProcessingFunctional (
	    new FlatAdiabaticBoundaryFunctional3D<T,ADESCRIPTOR,2,-1>,
	    leftThermal, adLattice, processorLevel );
   

    integrateProcessingFunctional (
	    new FlatAdiabaticBoundaryFunctional3D<T,ADESCRIPTOR,2,1>,
	    rightThermal, adLattice, processorLevel );
 
    //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>);

	defineDynamics(adLattice, adLattice.getBoundingBox(),
		new PeriodicGroovedDomain3D<T>(2*l + iM*12*l,parameters),
		new plb::BounceBack<T,ADESCRIPTOR>);
 
   }

    nsLattice.initialize();
    adLattice.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”]
void copyUnknownOnOutlet( MultiBlockLattice3D<T,ADESCRIPTOR>& adLattice,
SuperFlowParam<T,NSDESCRIPTOR,ADESCRIPTOR> const& parameters,
OnLatticeAdvectionDiffusionBoundaryCondition3D<T,ADESCRIPTOR>&
adBoundaryCondition )
{
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(
new CopyUnknownPopulationsFunctional3D<T,ADESCRIPTOR, 0, +1>,
    Box3D(nx-1,nx-1, 0,ny-1,0,nz-1), adLattice);

}



The relevant parts of main() are:

[code="cpp"]
///Generate lattices
    T nsOmega = parameters.getFluidOmega();
    T adOmega = parameters.getPSOmega();
    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); 
  
    MultiBlockLattice3D<T, ADESCRIPTOR> adLattice (
	    nx,ny, nz, new ADYNAMICS<T, ADESCRIPTOR>(adOmega) );
    adLattice.periodicity().toggleAll(false);
 
///Channel Setup
    pcout << "Setting up the channel..." << endl;
    OnLatticeBoundaryCondition3D<T,NSDESCRIPTOR>*
	nsBoundaryCondition = createLocalBoundaryCondition3D<T,NSDESCRIPTOR>();

    OnLatticeAdvectionDiffusionBoundaryCondition3D<T,ADESCRIPTOR>*
	adBoundaryCondition = createLocalAdvectionDiffusionBoundaryCondition3D<T,ADESCRIPTOR>();

    channelSetup(nsLattice,adLattice ,parameters, *nsBoundaryCondition, 
			*adBoundaryCondition, sections);


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

    copyUnknownOnOutlet(adLattice,parameters,*adBoundaryCondition); 
    
    pcout << "Finished setting up channel!" << endl;

where getLatticeGravity() is set to zero.

Regards,

Paul

#2

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

#3

Dear Paul,

It seems most convenient to first answer your last question:

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.

Does this solve your issue?

Cheers,
Jonas

#4

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(
new CopyUnknownPopulationsFunctional3D<T,ADESCRIPTOR, 0, +1>,
Box3D(nx-1,nx-1, 0,ny-1,0,nz-1), adLattice);


in main (and doing nothing else for the adLattice outlet). Do I have to implement a boundary condition with the help of my adBoundaryCondition object for the outlet, too perhaps?


Thanks for your time,

Paul