Skip to content

Advection Diffusion moddelling

OpenLB – Open Source Lattice Boltzmann Code Forums on OpenLB General Topics Advection Diffusion moddelling

Viewing 6 posts - 1 through 6 (of 6 total)
  • Author
    Posts
  • #9472
    HAdi
    Participant

    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.
    #9482
    mathias
    Keymaster

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

    #9485
    HAdi
    Participant

    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

    #9486
    mathias
    Keymaster

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

    #9487
    HAdi
    Participant

    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?

    #9488
    mathias
    Keymaster

    You need to set both!

Viewing 6 posts - 1 through 6 (of 6 total)
  • You must be logged in to reply to this topic.