Skip to content

Thermal Poiseuille 2D.

Viewing 8 posts - 1 through 8 (of 8 total)
  • Author
    Posts
  • #6331
    ramirofreile
    Participant

    Hello everyone,

    I am currently working on a Thermal Poiseuille 2D and I am observing a strange behavior in the temperature field, which initiates at the corners of the channel inlet.

    For simplicity, I set the temperature everywhere to a value of 1.5. However, higher temperatures appear at the corners and lower temperatures appear as in the nearby lattices. In the next 2 figures you can see both the geometry and the temperature at the inlet corners and bulk.

    https://ibb.co/YZvfMzB
    https://ibb.co/d6pDbVs

    For Dynamics I am using ForcedBGK for fluid and AdvectiondiffusionBGK for temperature.
    I am using the latest version of the code.

    To try and solve the issue I have tried:
    -Modifying dimensionless numbers (currently Re 100, and Pr 13). Lowering either makes the temperature deviation from 1.5 lower, but the issue is still there.
    -Playing with relaxation times according to Krueger book.
    -Modifying collision schemes, (BGK,TRT,MRT)
    -Modifying lattices Descriptors (D2Q5,D2Q9)

    Boundary conditions for the fluid flow eqs. are velocity inlet, pressure outlet and no-slip. In OpenLB:

    auto bulkIndicator = superGeometry.getMaterialIndicator({1, 3, 4});
    NSlattice.defineDynamics( bulkIndicator, &bulkDynamics );
    NSlattice.defineDynamics( superGeometry, 2, &instances::getBounceBack<T, NSDESCRIPTOR>() );
    setInterpolatedVelocityBoundary<T,NSDESCRIPTOR>(NSlattice, omega, superGeometry, 3);
    setInterpolatedPressureBoundary<T,NSDESCRIPTOR>(NSlattice, omega, superGeometry, 4);
    NSlattice.iniEquilibrium( bulkIndicator, rhoF, u );
    NSlattice.defineRhoU( bulkIndicator, rhoF, u );
    NSlattice.defineRhoU( superGeometry, 2, rhoF, uzero );
    NSlattice.iniEquilibrium( superGeometry, 2, rhoF, uzero );

    Boundary conditions for the energy Eq. are constant temperature everywhere. In Open LB:

    ADlattice.defineDynamics(superGeometry, 0, &instances::getNoDynamics<T, TDESCRIPTOR>());
    ADlattice.defineDynamics(superGeometry.getMaterialIndicator({1,2,3,4}), &advectionDiffusionBulkDynamics);
    setAdvectionDiffusionTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry, 2);
    setAdvectionDiffusionTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry, 3);
    setAdvectionDiffusionTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry, 4);

    ADlattice.defineRhoU(superGeometry, 2 , T_hot,uzero);
    ADlattice.defineRhoU(superGeometry, 3 , T_hot,u);
    ADlattice.defineRho(superGeometry, 4 , T_hot);
    ADlattice.iniEquilibrium(superGeometry, 1 , T_hot, u);
    ADlattice.iniEquilibrium(superGeometry, 3 , T_hot, u);
    ADlattice.iniEquilibrium(superGeometry, 4 , T_hot, u);
    ADlattice.iniEquilibrium(superGeometry, 2 , T_hot, uzero);
    NSlattice.initialize();
    ADlattice.initialize();

    I would appreciate any input on this problem.

    Thank you very much in advance,

    Ramiro

    #6335
    FBukreev
    Keymaster

    Hi Ramiro,

    by ADE Lattice it would be better to use defineRho() instead of defineRhoU(), because the velocity values are taken automatically from the NSE Lattice by the code itself.

    By definition of the initial velocities you could try to set 0 at the start. At the inlet the necessary velocity should be reached gradually by using the SinusStartScale functor. Such initialization should be better for convergency. If not, then write us again and we will look for an other solution.

    Greetings
    Fedor

    #6343
    ramirofreile
    Participant

    Dear Fedor,

    Thank you very much for your time.

    I spent these days trying both of your suggestions, defineRho() for the ADE lattice and using a SinusStartScale based off of the cylinder2d.cpp.

    Unfortunately I am still seeing the strange behavior in the temperature field. These are two images at different timesteps.

    https://ibb.co/3SFjSWL
    https://ibb.co/z6GpJwX

    For these images, the Reynolds number is 50; the Prandtl number is 13; the maximum lattice velocity is 0.0375; the relaxation time 0.55; Courant number 0.04.

    The problem seems to be starting at the cell right next to the inlet. The inlet is constantly kept at 1.5 (as it should be). However, in the cell directly next to it downstream, the temperature begins increasing to non-sense values. The temperature should be bounded between 1.5 and 0.5, as those are the imposed boundary conditions at the inlet and a portion of the wall, respectively.
    Thank you again,

    Ramiro Freile

    #6344
    ramirofreile
    Participant

    Dear Fedor,

    I want to inform you that changing the ADE lattice dynamics from AdvectionDiffusionTRTDynamics to AdvectionDiffusionBGKDynamics solved the problem. Do you have a lead on why this would happen?

    Ramiro Freile

    • This reply was modified 5 months, 4 weeks ago by ramirofreile.
    #6348
    ramirofreile
    Participant

    I have also noticed that in AdvectionDiffusionTRTDynamics the collision function is:

    cell[iPop] -= _omega2 * (fPlus[iPop] – fEqPlus[iPop]) + this->_omega * (fMinus[iPop] – fEqMinus[iPop]);

    whereas in TRTDynamics:

    cell[iPop] -= _omega * (fPlus[iPop] – fEqPlus[iPop]) + _omega2 * (fMinus[iPop] – fEqMinus[iPop]);

    The omega related to the viscosity and the free parameter omega are switched in the collision.

    #6387
    stephan
    Moderator

    Hi Ramiro,

    thanks for your message.
    That is interesting.
    We will have a look.

    I would be very happy if you could also give it a try, switch the relaxation times as in the NSE case, simulate again with AdvectionDifusionTRTDynamics, and post here again about the effect.

    However, the switch should be correct, since the NSE collision conserves moments up to order one, whereas the ADE does this only up to order zero, which in turn gives reason to use diffusion relaxation for the symmetric populations in the former and anti-symmetric ones in the latter application case.

    Please also note that the TRT scheme has variable numerical properties. Dependent on the second relaxation parameter you used, accuracy or stability could be increased; but both can also be decreased. Hence using the plain BGK model is just the right choice for your setting.

    BR
    Stephan

    #6396
    ramirofreile
    Participant

    Hello Stephan,

    I greatly appreciate your help and input on this.

    I tried switching the relaxation times for the ADE-TRTDynamics following the NSE-TRTDynamics, and the results were better. Here are some figures of the old and new results:

    Old Results without the switching:
    https://ibb.co/QQdmXP4

    New Results with the switching:
    https://ibb.co/ypP7hnF

    New Results decreasing the magic parameter
    https://ibb.co/JQVQ2Gt

    As a reminder, the boundary conditions for the temperature are set to 1.5 everywhere.
    In addition, you can find the code to recreate the results (It runs in ~20 seconds in 1 core).

    Thank you,

    Ramiro Freile

    ***************************************************************
    *********************CODE STARTS HERE**************************
    ***************************************************************

    #include “olb2D.h”
    #ifndef OLB_PRECOMPILED // Unless precompiled version is used,
    #include “olb2D.hh” // include full template code
    #endif
    #include <vector>
    #include <cmath>
    #include <iostream>
    #include <fstream>

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

    typedef double T;
    typedef D2Q9<FORCE> NSDESCRIPTOR;
    typedef D2Q5<VELOCITY> TDESCRIPTOR;

    //typedef D2Q9<tag::MRT,FORCE> NSDESCRIPTOR;
    //typedef D2Q5<tag::MRT,VELOCITY> TDESCRIPTOR;

    // Parameters for the simulation setup
    const T lengthX = 10.;
    const T lengthY = 1.;
    const T hydraulicD = 2.*lengthY;

    int N = 50.; // resolution of the model
    const T physDeltaX = lengthY/N; // latticeL

    T tau = 0.6; // relaxation time

    const T Re = 100;
    const T rho = 2050.;
    const T mu = 5.7e-3;
    const T nu = mu/rho;
    const T u_inlet = Re*nu/hydraulicD;

    const T maxPhysT = 15000.; // simulation time
    const T output_frequency = 300.;

    //Thermal Sector
    const T Tcold = 0.5; //in K
    const T Thot = 1.5; //in K
    const T lambda_l = 0.83; // W / m K
    const T cp_real = 1900.; // J / kg K

    // Stores geometry information in form of material numbers
    void prepareGeometry( ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> const& converter,
    SuperGeometry2D<T>& superGeometry )
    {

    OstreamManager clout( std::cout,”prepareGeometry” );
    clout << “Prepare Geometry …” << std::endl;

    Vector<T,2> extend( lengthX,lengthY );
    Vector<T,2> origin;

    superGeometry.rename( 0,2 );

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

    extend[0] = physDeltaX;
    extend[1] = lengthY;
    origin[0] = 0.;
    origin[1] = -physDeltaX;
    IndicatorCuboid2D<T> inflow( extend, origin );
    superGeometry.rename( 2,3,1,inflow );
    // Set material number for
    extend[0] = physDeltaX;
    extend[1] = lengthY+physDeltaX/2.;
    origin[0] = lengthX-physDeltaX/2.;
    origin[1] = -physDeltaX/2.;

    IndicatorCuboid2D<T> outflow( extend, origin );
    superGeometry.rename( 2,4,1,outflow );

    // Removes all not needed boundary voxels outside the surface
    superGeometry.clean();
    superGeometry.checkForErrors();

    superGeometry.print();

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

    // Set up the geometry of the simulation
    void prepareLattice( SuperLattice2D<T, NSDESCRIPTOR>& NSlattice,
    SuperLattice2D<T, TDESCRIPTOR>& ADlattice,
    ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> const& converter,
    Dynamics<T, NSDESCRIPTOR>& bulkDynamics,
    Dynamics<T, TDESCRIPTOR>& advectionDiffusionBulkDynamics,
    SuperGeometry2D<T>& superGeometry )
    {

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

    const T omega = converter.getLatticeRelaxationFrequency();
    T Tomega = converter.getLatticeThermalRelaxationFrequency();

    // Material=0 –>do nothing
    NSlattice.defineDynamics( superGeometry, 0, &instances::getNoDynamics<T, NSDESCRIPTOR>() );

    auto bulkIndicator = superGeometry.getMaterialIndicator({1, 3, 4});
    NSlattice.defineDynamics( bulkIndicator, &bulkDynamics );

    NSlattice.defineDynamics( superGeometry, 2, &instances::getBounceBack<T, NSDESCRIPTOR>() );

    //if boundary conditions are chosen to be local
    setInterpolatedVelocityBoundary<T,NSDESCRIPTOR>(NSlattice, omega, superGeometry, 3);
    setInterpolatedPressureBoundary<T,NSDESCRIPTOR>(NSlattice, omega, superGeometry, 4);

    // Initial conditions
    AnalyticalConst2D<T,T> rhoF( 1 );
    std::vector<T> velocity( 2,T( 0 ) );
    AnalyticalConst2D<T,T> uzero( velocity );

    velocity[0] = converter.getCharLatticeVelocity();
    AnalyticalConst2D<T,T> u( velocity );

    NSlattice.iniEquilibrium( bulkIndicator, rhoF, u );
    NSlattice.defineRhoU( bulkIndicator, rhoF, u );
    NSlattice.defineRhoU( superGeometry, 2, rhoF, uzero );
    NSlattice.iniEquilibrium( superGeometry, 2, rhoF, uzero );

    ////Thermal Begins here
    ADlattice.defineDynamics(superGeometry, 0, &instances::getNoDynamics<T, TDESCRIPTOR>());

    ADlattice.defineDynamics(superGeometry.getMaterialIndicator({1,2,3,4}), &advectionDiffusionBulkDynamics);

    setAdvectionDiffusionTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry, 2);
    setAdvectionDiffusionTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry, 3);
    setAdvectionDiffusionTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry, 4);

    AnalyticalConst2D<T,T> T_hot(Thot);

    ADlattice.defineRho(superGeometry, 2 , T_hot);
    ADlattice.defineRho(superGeometry, 3 , T_hot);
    ADlattice.defineRho(superGeometry, 4 , T_hot);

    ADlattice.iniEquilibrium(superGeometry, 1 , T_hot, u);
    ADlattice.iniEquilibrium(superGeometry, 3 , T_hot, u);
    ADlattice.iniEquilibrium(superGeometry, 4 , T_hot, u);
    ADlattice.iniEquilibrium(superGeometry, 2 , T_hot, uzero);

    // Make the lattice ready for simulation
    NSlattice.initialize();
    ADlattice.initialize();

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

    // Computes the pressure drop between the voxels before and after the cylinder
    void getResults( SuperLattice2D<T, NSDESCRIPTOR>& NSlattice,
    SuperLattice2D<T, TDESCRIPTOR>& ADlattice,
    ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> const& converter, int iT,
    SuperGeometry2D<T>& superGeometry, Timer<T>& timer )
    {

    SuperVTMwriter2D<T> vtmWriter( “channel2d” );
    SuperLatticeGeometry2D<T, NSDESCRIPTOR> geometry(NSlattice, superGeometry);
    SuperLatticePhysVelocity2D<T, NSDESCRIPTOR> velocity( NSlattice, converter );
    SuperLatticePhysPressure2D<T, NSDESCRIPTOR> pressure( NSlattice, converter );
    SuperLatticePhysTemperature2D<T, NSDESCRIPTOR, TDESCRIPTOR> temperature(ADlattice, converter);
    vtmWriter.addFunctor( velocity );
    vtmWriter.addFunctor( pressure );
    vtmWriter.addFunctor( temperature );
    vtmWriter.addFunctor( geometry );

    const int vtkIter = int (output_frequency/converter.getPhysDeltaT() );

    if ( iT == 0 ) {
    // Writes the geometry, cuboid no. and rank no. as vti file for visualization
    SuperLatticeGeometry2D<T, NSDESCRIPTOR> geometry( NSlattice, superGeometry );
    SuperLatticeCuboid2D<T, NSDESCRIPTOR> cuboid( NSlattice );
    SuperLatticeRank2D<T, NSDESCRIPTOR> rank( NSlattice );
    vtmWriter.write( geometry );
    vtmWriter.write( cuboid );
    vtmWriter.write( rank );

    vtmWriter.createMasterFile();
    }

    /// Writes the VTK files
    if (iT % vtkIter == 0) {

    timer.update(iT);
    timer.printStep();

    /// NSLattice statistics console output
    NSlattice.getStatistics().print(iT,converter.getPhysTime(iT));

    /// ADLattice statistics console output
    ADlattice.getStatistics().print(iT,converter.getPhysTime(iT));

    if ( NSlattice.getStatistics().getAverageRho() != NSlattice.getStatistics().getAverageRho() or ADlattice.getStatistics().getAverageRho() != ADlattice.getStatistics().getAverageRho() ) {
    cout << “simulation diverged! stopping now.” << endl;
    exit(1);
    }
    vtmWriter.write(iT);
    }

    }

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

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

    const T physDeltaT = 1 / nu / descriptors::invCs2<T,NSDESCRIPTOR>() * (tau – 0.5) * physDeltaX * physDeltaX ;

    ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> const converter(
    (T) physDeltaX, // physDeltaX
    (T) physDeltaT, // physDeltaT
    (T) hydraulicD, // charPhysLength
    (T) u_inlet, // charPhysVelocity
    (T) nu, // physViscosity
    (T) rho, // physDensity
    (T) lambda_l, // physThermalConductivity
    (T) cp_real, // physSpecificHeatCapacity
    (T) 1.0e-4, // physThermalExpansionCoefficient
    (T) Tcold, // charPhysLowTemperature
    (T) Thot // charPhysHighTemperature
    );
    converter.print();

    // === 2rd Step: Prepare Geometry ===
    Vector<T,2> extend( lengthX,lengthY );
    Vector<T,2> origin;
    IndicatorCuboid2D<T> cuboid( extend, origin );

    // Instantiation of a cuboidGeometry with weights
    #ifdef PARALLEL_MODE_MPI
    const int noOfCuboids = singleton::mpi().getSize();
    #else
    const int noOfCuboids = 7;
    #endif
    CuboidGeometry2D<T> cuboidGeometry( cuboid, physDeltaX, noOfCuboids );

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

    // Instantiation of a superGeometry
    SuperGeometry2D<T> superGeometry( cuboidGeometry, loadBalancer, 2 );

    prepareGeometry( converter, superGeometry );

    // === 3rd Step: Prepare Lattice ===
    SuperLattice2D<T, NSDESCRIPTOR> NSlattice( superGeometry );
    SuperLattice2D<T, TDESCRIPTOR> ADlattice(superGeometry);

    ForcedBGKdynamics<T, NSDESCRIPTOR> NSbulkDynamics(
    converter.getLatticeRelaxationFrequency(),
    instances::getBulkMomenta<T,NSDESCRIPTOR>()
    );

    AdvectionDiffusionTRTdynamics<T, TDESCRIPTOR> TbulkDynamics (
    converter.getLatticeThermalRelaxationFrequency(),
    instances::getAdvectionDiffusionBulkMomenta<T,TDESCRIPTOR>(),
    1./4..
    );

    std::vector<T> dir{0.0, 1.0};

    NavierStokesAdvectionDiffusionCouplingGenerator2D<T,NSDESCRIPTOR>
    coupling(0, converter.getLatticeLength(lengthX),
    0, converter.getLatticeLength(lengthY),
    0.0, Tcold, 1., dir);

    NSlattice.addLatticeCoupling(coupling, ADlattice);

    prepareLattice( NSlattice, ADlattice, converter, NSbulkDynamics, TbulkDynamics, superGeometry );

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

    for ( std::size_t iT = 0; iT < converter.getLatticeTime( maxPhysT ); ++iT ) {

    NSlattice.executeCoupling();
    NSlattice.collideAndStream();
    ADlattice.collideAndStream();

    // === 7th Step: Computation and Output of the Results ===
    getResults( NSlattice, ADlattice, converter, iT, superGeometry, timer );
    }

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

    ***************************************************************
    **********************CODE ENDS HERE***************************
    ***************************************************************

    #6541
    stephan
    Moderator

    Dear Ramiro,

    thanks for reporting back.
    I tested your example on OpenLB r1.4 and can reproduce the results.
    The appearing effect could also be a result of boundary/inlet conditions or dynamics coupling.

    Since we have realized a new dynamics concept in r1.5, I would suggest to update your code according to the new release.
    If this does not change anything, please let me know again.

    Thank you again for contributing!

    BR
    Stephan

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