Hi Everyone,
I finally managed to solve (almost) the advection-diffusion problem. I also managed to couple the concentration field with a velocity field. I thought that maybe can be useful to somebody and I posted the code here. The code solve the concentration field within a poiseuille flow. The coupling is passive and uses a function from the palabos library.
[code=“cpp”]
#include
#include “palabos2D.h”
#include “palabos2D.hh”
using namespace std;
using namespace plb;
typedef double T;
#define NSDES descriptors::D2Q9Descriptor
#define AVDES descriptors::AdvectionDiffusionD2Q5Descriptor
// ************************************************************************
// Part for the NS part of the code. A classic poiseuille flow is simulated.
// Most of the code is taken from Palabos tutorial 1.5
// ************************************************************************
// Poiselle velocity analitical
// Velocity on the parabolic Poiseuille profile
T poiseuilleVelocity(plint iY, IncomprFlowParam const& parameters) {
T y = (T)iY / parameters.getResolution();
return 4.parameters.getLatticeU() * (y-yy);
}
// A functional, used to initialize the velocity for the boundary conditions
template
class PoiseuilleVelocity {
public:
PoiseuilleVelocity(IncomprFlowParam parameters_)
: parameters(parameters_)
{ }
// This version of the operator returns the velocity only,
// to instantiate the boundary condition.
void operator()(plint iX, plint iY, Array<T,2>& u) const {
u[0] = poiseuilleVelocity(iY, parameters);
u[1] = T();
}
// This version of the operator returns also a constant value for
// the density, to create the initial condition.
void operator()(plint iX, plint iY, T& rho, Array<T,2>& u) const {
u[0] = poiseuilleVelocity(iY, parameters);
u[1] = T();
rho = (T)1;
}
private:
IncomprFlowParam parameters;
};
// Function to set up the channel
void channelSetup (
MultiBlockLattice2D<T,NSDES>& lattice,
IncomprFlowParam const& parameters,
OnLatticeBoundaryCondition2D<T,NSDES>& boundaryCondition )
{
// Create Velocity boundary conditions.
boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice);
// Specify the boundary velocity.
setBoundaryVelocity (
lattice, lattice.getBoundingBox(),
PoiseuilleVelocity<T>(parameters) );
// Create the initial condition.
initializeAtEquilibrium (
lattice, lattice.getBoundingBox(), PoiseuilleVelocity<T>(parameters) );
lattice.initialize();
}
// ****************************************************************************
// ****************************************************************************
// ************************************************************************
// Advection-diffusion setup *********************************************
// ************************************************************************
void concentrationSetup(MultiBlockLattice2D<T,AVDES> &lattice,
OnLatticeAdvectionDiffusionBoundaryCondition2D<T, AVDES> &BC,
IncomprFlowParam parameters) {
// Get lattice dimensions
const plint nx = parameters.getNx();
const plint ny = parameters.getNy();
// Get velocity and rho for this lattice
Array<T,2> u0(0.,0.);
// Set-up the geometry
//Box2D topwall(0, nx-1, ny-1, ny-1);
//Box2D bottomwall(0, nx-1, 0, 0);
Box2D leftwall(0,1, 0,ny);
Box2D rightwall(nx,nx, 0,ny);
// Set the boundary-conditions
BC.addTemperatureBoundary1P(rightwall, lattice);
BC.addTemperatureBoundary1N(leftwall, lattice);
// Try to impose a different temperature
setBoundaryDensity(lattice, rightwall, 0.0);
setBoundaryDensity(lattice, leftwall, 1.0);
// Init lattice
initializeAtEquilibrium(lattice, lattice.getBoundingBox(), 0.0,u0);
// Command to init the lattice
lattice.periodicity().toggle(1,true);
lattice.initialize();
}
// ************************************************************************
// ************************************************************************
// ************************************************************************
int main(int argc, char *argv[])
{
// Init simulation
plbInit(&argc, &argv);
// Parameters for the dynamics
IncomprFlowParam<T> parameters(
(T) 1e-2, // Reference velocity (the maximum velocity
// in the Poiseuille profile) in lattice units.
(T) 100., // Reynolds number
100, // Resolution of the reference length (channel height).
2., // Channel length in dimensionless variables
1. // Channel height in dimensionless variables
);
Array<T,2> u2(2.0, 0.0);
// Istantiate the lattice
MultiBlockLattice2D<T, NSDES> nsLattice(parameters.getNx(), parameters.getNy(), new BGKdynamics<T, NSDES>(parameters.getOmega()));
// Istantiate adv-diff lattice
MultiBlockLattice2D<T, AVDES> advLattice(parameters.getNx(), parameters.getNy(), new AdvectionDiffusionBGKdynamics<T, AVDES>(1.0));
// Istantiate the boundary conditions
OnLatticeBoundaryCondition2D<T,NSDES> *nsBC = createLocalBoundaryCondition2D<T, NSDES>();
OnLatticeAdvectionDiffusionBoundaryCondition2D<T, AVDES> *avBC = createLocalAdvectionDiffusionBoundaryCondition2D<T, AVDES>();
// Set-up channel for NS equations
channelSetup(nsLattice, parameters, *nsBC);
//Set-up concentration field
concentrationSetup(advLattice, *avBC, parameters);
// Main LBM cycle
ImageWriter<T> image("leeloo");
for (plint iT=0; iT < 10000; ++iT) {
if(iT%100==0) {
pcout << "Writing GIF... \n";
//image.writeScaledGif(createFileName("conc",iT,2), *computeDensity(advLattice), parameters.getNx(),parameters.getNy());
//image.writeScaledGif(createFileName("vel",iT, 2), *computeVelocityNorm(nsLattice), parameters.getNx(), parameters.getNy());
}
nsLattice.collideAndStream();
}//*/
// First attempt to couple the two physics
latticeToPassiveAdvDiff(nsLattice, advLattice, advLattice.getBoundingBox());
for (plint iT=0; iT < 10000; ++iT) {
if(iT%100==0) {
pcout << "Writing GIF... \n";
//image.writeScaledGif(createFileName("conc",iT,2), *computeDensity(advLattice), parameters.getNx(),parameters.getNy());
//image.writeScaledGif(createFileName("vel",iT, 2), *computeVelocityNorm(nsLattice), parameters.getNx(), parameters.getNy());
}
advLattice.collideAndStream();
}//*/
image.writeScaledGif(createFileName("adv-diff",1,1), *computeDensity(advLattice), parameters.getNx(), parameters.getNy());
pcout << "Hello palabos";
return 0;
}
The last step will be implement the advection diffusion with source. If any of you have some suggestions they will be appreciated.
Cheers,
Salvo