Advection Diffusion moddelling

[HAdi]

Dear OpenLB community,

I want to simulate solute transport with OpenLB. I start with the most simple problem in which I want to see the transport of a spherical solute in a pipe with Poiseuille velocity profile. I use SmoothIndicatorSphere3D to create a sphere in the middle of domain and I use sLatticeAD.defineField to add the Poiseuille velocity profile to it. But When I simulate it I can see from the results that the solute is only diffusing, it does not move with the velocity. Below is my code. Thank you in advance for your help.

#include “olb3D.h”
#include “olb3D.hh”

using namespace olb;
using namespace olb::descriptors;
using namespace olb::graphics;
using namespace olb::util;

using T = FLOATING_POINT_TYPE;
typedef D3Q7<VELOCITY> ADDESCRIPTOR;

// Parameters for simulation setup

const T length = 2.4e-3; // length of the pipe
const T diameter = 0.4e-3; // diameter of the pipe
const int N = 40; // resolution of the model
const T physU = 2.95e-4; // physical velocity
const T Re = 7.08e-2; // Reynolds number
const T physRho = 1.; // physical density
const T tau = 0.8; // lattice relaxation time
const T maxPhysT = 20.; // max. simulation time in s, SI unit

// Scaled Parameters
const T radius = diameter/2.; // radius of the pipe
const T physInterval = 0.0125*maxPhysT; // interval for the convergence check in s
const T diffusion = 1.e-8; // diffusion coefficient for advection-diffusion equation in m2/s

// These below constants are for developing
// SS flow field before particle injection

const T residuum = 3.e-1; // residdum for the convergence check
int iTinjection = 0; // injection time in s

const int iTperiod = 50; // amount of timesteps when new boundary conditions are reset and results are visualized

void prepareGeometry( UnitConverter<T,ADDESCRIPTOR> const& converter,
SuperGeometry<T,3>& superGeometry )
{

OstreamManager clout( std::cout, “prepareGeometry” );

clout << “Prepare Geometry …” << std::endl;

// Material 0 –> Isolated
// Material 1 –> Fluid Domain
// Material 2 –> Bounce Back

Vector<T, 3> center0(-converter.getPhysDeltaX() * 0.2, radius, radius);
Vector<T, 3> center1(length, radius, radius);

center0[0] -= 3.*converter.getPhysDeltaX();
center1[0] += 3.*converter.getPhysDeltaX();

IndicatorCylinder3D<T> pipe(center0, center1, radius);

superGeometry.rename(0, 2);
superGeometry.rename(2, 1, pipe);

superGeometry.clean();
superGeometry.innerClean();
superGeometry.checkForErrors();

superGeometry.print();

clout << “Prepare Geometry … OK” << std::endl;
return;
}

void prepareLattice( SuperLattice<T, ADDESCRIPTOR>& sLatticeAD,
UnitConverter<T,ADDESCRIPTOR> const& converter,
SuperGeometry<T,3>& superGeometry,
T omegaAD )
{

OstreamManager clout( std::cout, “prepareLattice” );
clout << “Prepare Lattice …” << std::endl;

// Material=1 –> bulk dynamic (AD-BGK)
sLatticeAD.defineDynamics<AdvectionDiffusionBGKdynamics>(superGeometry, 1);

// Material=2 –> zero distribution dynamic
sLatticeAD.defineDynamics<ZeroDistributionDynamics>(superGeometry, 2);

// AD boundary conditions
setZeroDistributionBoundary<T,ADDESCRIPTOR>(sLatticeAD, superGeometry, 2);

/// Initial conditions
// Setting the Poiseuille flow field
std::vector<T> origin = { length, radius, radius};
std::vector<T> axis = { 1, 0, 0 };
CirclePoiseuille3D<T> poiseuilleU(origin, axis, converter.getCharLatticeVelocity(), radius);
sLatticeAD.defineField<descriptors::VELOCITY>(superGeometry, 1, poiseuilleU);

AnalyticalConst3D<T,T> rho0( 0 );
AnalyticalConst3D<T,T> rho1( 1. );

Vector<T, 3> centerp(length/4, diameter/2, diameter/2);
IndicatorSphere3D<T> particle(centerp, radius/3);
SmoothIndicatorSphere3D<T,T> section1( particle, 0. );

AnalyticalIdentity3D<T,T> rho(section1 – rho0);

// Initialize all values of distribution functions to their local equilibrium
sLatticeAD.defineRho( superGeometry.getMaterialIndicator(1), rho );
sLatticeAD.defineU(superGeometry, 1, poiseuilleU);
sLatticeAD.iniEquilibrium( superGeometry.getMaterialIndicator({1}), rho, poiseuilleU );

sLatticeAD.setParameter<descriptors::OMEGA>(omegaAD);

// Lattice initialize
sLatticeAD.initialize();

{
auto& communicator = sLatticeAD.getCommunicator(stage::PostCoupling());
communicator.requestField<descriptors::VELOCITY>();
communicator.requestOverlap(sLatticeAD.getOverlap());
communicator.exchangeRequests();
}

clout << “Prepare Lattice … OK” << std::endl;
return;
}

void getResults( SuperLattice<T, ADDESCRIPTOR>& sLatticeAD,
UnitConverter<T,ADDESCRIPTOR> const& converter, int iT,
SuperGeometry<T,3>& superGeometry,
Timer<double>& timer )
{
OstreamManager clout( std::cout, “getResults” );
SuperVTMwriter3D<T> vtmWriterAD( “NSAD_particle” );

SuperLatticeDensity3D<T, ADDESCRIPTOR> particles( sLatticeAD );
SuperLatticePhysField3D<T, ADDESCRIPTOR, descriptors::VELOCITY> extField(
sLatticeAD, converter.getConversionFactorVelocity());
vtmWriterAD.addFunctor( particles );
vtmWriterAD.addFunctor( extField );

if ( iT == 0 ) {
SuperLatticeGeometry3D<T, ADDESCRIPTOR> geometry( sLatticeAD, superGeometry );
SuperLatticeCuboid3D<T, ADDESCRIPTOR> cuboid( sLatticeAD );
SuperLatticeRank3D<T, ADDESCRIPTOR> rank( sLatticeAD );
vtmWriterAD.write( geometry );
vtmWriterAD.write( cuboid );
vtmWriterAD.write( rank );
vtmWriterAD.createMasterFile();

}

if ( iT % iTperiod == 0 ) {
// Writes the vtk files
vtmWriterAD.write( iT );

// Writes output on the console
timer.update( iT );
timer.printStep();
sLatticeAD.getStatistics().print( iT, iT * converter.getCharLatticeVelocity() / T(converter.getResolution()) );

}
}

int main( int argc, char* argv[] )
{

// === 1st Step: Initialization ===

olbInit( &argc, &argv );
singleton::directories().setOutputDir( “./tmp/” );
OstreamManager clout( std::cout, “main” );

UnitConverterFromResolutionAndRelaxationTime<T,ADDESCRIPTOR> const converter(
int {N}, // resolution: number of voxels per charPhysL
(T) tau, // latticeRelaxationTime: relaxation time, have to be greater than 0.5!
(T) diameter, // charPhysLength: reference length of simulation geometry
(T) physU, // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__
(T) diameter*physU/Re, // physViscosity: physical kinematic viscosity in __m^2 / s__
(T) physRho // physDensity: physical density in __kg / m^3__
);

// Prints the converter log as console output
converter.print();
// Writes the converter log in a file
converter.write(“bifurcation3d”);

// compute relaxation parameter to solve the advection-diffusion equation in the lattice Boltzmann context
T omegaAD = converter.getLatticeRelaxationFrequencyFromDiffusivity<ADDESCRIPTOR>( diffusion );
clout << “tau AD: ” << (1.0 / omegaAD) << std::endl;

// === 2nd Step: Prepare Geometry ===

Vector<T, 3> center0(0, radius, radius);
Vector<T, 3> center1(length + 0.5 * converter.getPhysDeltaX(), radius, radius);
IndicatorCylinder3D<T> pipe(center0, center1, radius);
IndicatorLayer3D<T> extendedDomain(pipe, converter.getPhysDeltaX());

// Instantiation of an empty cuboidGeometry
int noOfCuboids = util::max( 20, singleton::mpi().getSize() );
CuboidGeometry3D<T> cuboidGeometry( extendedDomain, converter.getConversionFactorLength(),
noOfCuboids, “weight” );
cuboidGeometry.printExtended();
clout << “min / max ratio (volume) = ” << (T) cuboidGeometry.getMinLatticeVolume()
/ cuboidGeometry.getMaxLatticeVolume() << std::endl;
clout << “min / max ratio (weight) = ” << (T) cuboidGeometry.getMinLatticeWeight()
/ cuboidGeometry.getMaxLatticeWeight() << std::endl;

cuboidGeometry.setPeriodicity( true, false, false );

// Instantiation of an empty loadBalancer
HeuristicLoadBalancer<T> loadBalancer( cuboidGeometry );

// Default instantiation of superGeometry
SuperGeometry<T,3> superGeometry( cuboidGeometry, loadBalancer, 2);

prepareGeometry( converter, superGeometry );

// === 3rd Step: Prepare Lattice ===
// SuperLattice<T, NSDESCRIPTOR> sLatticeNS( superGeometry );

SuperLattice<T, ADDESCRIPTOR> sLatticeAD( superGeometry );

//prepareLattice and setBoundaryConditions
prepareLattice(sLatticeAD, converter, superGeometry, omegaAD );

// === 4th Step: Main Loop with Timer ===
clout << “starting simulation…” << std::endl;
Timer<double> timer( converter.getLatticeTime( maxPhysT ),
superGeometry.getStatistics().getNvoxel() );
util::ValueTracer<T> converge (converter.getLatticeTime (physInterval), residuum);
timer.start();

for (std::size_t iT = 0; iT <= converter.getLatticeTime(maxPhysT); ++iT) {
sLatticeAD.setParameter<descriptors::LATTICE_TIME>(iT);
getResults(sLatticeAD, converter, iT, superGeometry, timer );

sLatticeAD.collideAndStream();
}
timer.stop();
timer.printSummary();

}

• This topic was modified 1 month ago by HAdi.

[mathias]

A good starting point is the bifukation example in OpenLB..

[HAdi]

I worked with it, and I could run my problem with this example successfully, but when I do not want to solve the fluid flow (Navier Stokes) and defied the velocity field, I face problem. Could you please help me in this regard.

Thanks,
Hadi

[mathias]

It might be that you need a second internal velocity also be initialized..

[HAdi]

You are right if I use ParticleAdvectionDiffusionBGK(VELOCITY, VELOCITY2) (D3Q7 <VELOCITY, VELOCITY2>) with 2 velocity profile. The second velocity profile as long as I understand is used to determine drag force on the particles. However, in my case I do not want to consider this, so I used AdvectionDiffusionBGK(VELOCITY) (D3Q7 <VELOCITY>) with only one velocity field. Do I miss something?

[mathias]

You need to set both!

[Author]

Posts