Skip to content

Liquid-Solid Conjugate Heat Transfer

OpenLB – Open Source Lattice Boltzmann Code Forums on OpenLB General Topics Liquid-Solid Conjugate Heat Transfer

Viewing 8 posts - 1 through 8 (of 8 total)
  • Author
    Posts
  • #5950
    wf_guoyl
    Participant

    Hello, everyone! I have some trouble with the liquid-solid conjugate heat transfer. I have defined two descriptors: NSDESCRIPTOR and TDESCRIPTOR. During thermal simulations, I define two different “omega” to simulate transfer of heat in fluids and solids. The solid part is called “obstacle”. Here is part of my code:

    /***************************************************/
    prepare Lattice:

    Vector<SimType,2> center( lx0/2, ly0/2 );
    IndicatorCircle2D<SimType> obstacle( center, 25*conversionLength );

    auto obstacleIndicator = superGeometry.getMaterialIndicator({materialObstacle});

    // define NS Dynamics
    NSlattice.defineDynamics( superGeometry, materialObstacle, &olb::instances::getNoDynamics<SimType, NSDESCRIPTOR>() );
    offBC.addZeroVelocityBoundary( superGeometry, materialObstacle, obstacle );

    // define AD Dynamics
    ADlattice.defineDynamics( fluidIndicator, &advectionDiffusionBulkDynamics );
    ADlattice.defineDynamics( obstacleIndicator, &advectionDiffusionObstacleDynamics );

    // fluid init
    There’s no need to change it.

    // thermal init
    ADlattice.defineRho( obstacleIndicator, T_init );
    ADlattice.iniEquilibrium( obstacleIndicator, T_init, zeroVelocity );

    /*************************************************/

    But the simulation didn’t work. The error message is “density is NAN”. Specifically, four white blocks appear on the top, bottom, left and right endpoints of the circle obstacle. I don’t know what caused it. Does anyone want to help me? If someone can provide some information on fluid-solid conjugate heat transfer, I would appreciate it very much. My email is wfguoyl@163.com .

    #5953
    Adrian
    Keymaster

    The reason for this is likely an issue in the boundary setup, i.e. not all links between undefined no-dynamics nodes and fluid nodes are reconstructed correctly. This is a relatively common problem that can occur when using Bouzidi boundaries. One trick that helps surprisingly often is to change the discretization resolution by e.g. a single cell. Do you get any warnings printed during simulation setup?

    If you can post your full code I’d be happy to take a closer look.

    #5988
    wf_guoyl
    Participant

    Hello, Adrian:

    Thank you very much for your reply!
    Because there are a lot of documents involved, I have sent the documents to your email.
    My model is modified based on Julius Weinmiller’s work, which can be found on https://www.openlb.net/forum/users/julius-weinmiller/.

    In the code below, I adjust the “Bouzidi” boundary condition to “BounceBack” boundary condition. But the solid part “obstacle” still does not conduct heat.

    The following is the main code:

    /*****************************************************/

    /* OpenLB Extension:
    * Phase-change multiphase flow
    * Written by Julius Weinmiller for his thesis
    * https://www.openlb.net/forum/users/julius-weinmiller/
    *
    * The most recent release of OpenLB can be downloaded at
    * <http://www.openlb.net/&gt;
    */

    /* D2Law.cpp:
    * This simulation will combine both AD and NS lattice for a multiphase problem
    * A droplet is spended without body forces
    * periodic boundaries are applied in all three axis
    * The termpature is increased and evaporation rate recorded
    * The results of the simulations show that the simulations follow the
    * well known D2 evaporative law.
    */

    #include “olb2D.h”
    #include “olb2D.hh” // include full template code
    #include “phaseChangeExtension.h”
    #include “phaseChangeExtension.hh”

    using namespace olb;
    #define ThermalMRT

    typedef double SimType;
    #define NSDESCRIPTOR olb::descriptors::D2Q9<olb::descriptors::VELOCITY, \
    olb::descriptors::FORCE, \
    olb::descriptors::EXTERNAL_FORCE, \
    olb::descriptors::PSI_PSEUDO_RHO>

    #ifdef ThermalMRT
    #define TDESCRIPTOR olb::descriptors::D2Q5<olb::descriptors::tag::MRT, \
    olb::descriptors::VELOCITY, \
    olb::descriptors::PHI_THERMAL, \
    olb::descriptors::PREV_T_V, \
    olb::descriptors::PREV_PHI>
    #else
    #define TDESCRIPTOR olb::descriptors::D2Q5<olb::descriptors::VELOCITY, \
    olb::descriptors::PHI_THERMAL, \
    olb::descriptors::PREV_T_V, \
    olb::descriptors::PREV_PHI>
    #endif

    // Parameters for the simulation setup

    // Resolutions
    const int Nx = 100; // spacial resolution of the model
    const SimType Nt = 7; // temporal resolution of the model

    const SimType physLength = 100e-6; // m
    const SimType lx0 = 2*physLength; // length of channel
    const SimType ly0 = 2*physLength; // height of channel
    const SimType conversionLength = physLength/Nx ;

    const SimType sphereRadius = physLength/5.;
    const SimType sphereInterfaceWidth = sphereRadius/4;
    const SimType sphereOrigin2[2] = {lx0/2., ly0/2.};
    const SimType sphereOrigin[2] = {lx0/2., sphereRadius+sphereInterfaceWidth};
    const SimType sphereRadius2 = 25*conversionLength;
    const SimType sphereInterfaceWidth2 = sphereRadius2/4;

    const SimType rhoVapourMultiplier = 1.0;
    const SimType rhoLiquidMultiplier = 0.99;

    const SimType physViscosity = 1e-5; // m2 / s

    const SimType maxVelocity = 20;
    const SimType charTime = physLength/maxVelocity;

    const SimType endTimeInit = 10*charTime; // Time to allow droplet to come to rest
    const SimType endTimeEvap = 2010.0*charTime; // max. simulation time in s, SI unit

    const SimType dx = physLength/ Nx;
    const SimType dt = charTime/Nx/Nt;

    ////
    // Thermal constants
    const SimType TempCrit = 647; // [K] critical temperature of water in Kelvin

    const SimType specificGasConstantVapour = 461.52;
    const SimType propertyTemperatureRatio = 0.6;
    const SimType physWaterTemp = TempCrit * propertyTemperatureRatio – 273.15;
    // http://www.thermopedia.com/content/1150/
    // Table 8 for thermal conductivity
    // Table 4 for specific heat at constant pressure
    // Table 1 for density
    const SimType physDensityLiquid = 947.13;
    const SimType specificHeatCapacityVolumeLiquid = 3.8e3; // J / kg K
    const SimType thermalConductivityLiquid = 682e-3; // W / m K
    const SimType physDensityVapour = 0.9643;
    const SimType specificHeatCapacityVolumeVapour = 2.3e3 – specificGasConstantVapour; // Cv = Cp – Rsp REAL
    const SimType thermalConductivityVapour = 31e-3; // W / m K

    const SimType liquidDensityRatio = physDensityLiquid/8.72518;

    const SimType specificHeatCapacityVolume = 4e3; // J / kg K
    const SimType thermalConductivity = 4; // W / m K

    const SimType wallRho = 5.0;
    const SimType initLiquidHeightFraction = 0.4;

    const int iTsmoothUpdate = 10;
    const int rhoRampIter = 5000;

    const int materialBase = 0;
    const int materialFluid = 1;
    const int materialLowerWall = 2;
    const int materialUpperWall = 3;
    const int materialInflow = 4;
    const int materialOutflow = 5;
    const int materialObstacle = 6;

    void prepareGeometry( SuperGeometry2D<SimType>& superGeometry )
    {

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

    superGeometry.rename( materialBase, materialLowerWall );

    Vector<SimType,2> extend( lx0, ly0 );
    Vector<SimType,2> origin;

    // create fluid and lower wall
    superGeometry.rename( materialLowerWall, materialFluid, 0, 1 );

    // create upper wall
    origin = { 0, ly0-conversionLength };
    extend = { lx0, conversionLength };
    olb::IndicatorCuboid2D<SimType> upperBC( extend, origin );
    superGeometry.rename( materialLowerWall, materialUpperWall, materialFluid, upperBC );

    // create inflow
    origin = { 0, 0 };
    extend = { conversionLength, ly0-conversionLength };
    olb::IndicatorCuboid2D<SimType> inflow( extend, origin );
    superGeometry.rename( materialFluid, materialInflow, inflow );

    // create outflow
    origin = { lx0-conversionLength, 0 };
    extend = { conversionLength, ly0-conversionLength };
    olb::IndicatorCuboid2D<SimType> outflow( extend, origin);
    superGeometry.rename( materialFluid, materialOutflow, outflow);

    // creat obstacle
    origin = { lx0/2-50*conversionLength, ly0/2 };
    extend = { 100*conversionLength, 25*conversionLength };
    olb:: IndicatorCuboid2D<SimType> obstacle( extend, origin);
    superGeometry.rename( materialFluid, materialObstacle, obstacle);

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

    superGeometry.print();

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

    // Set up the geometry of the simulation
    void prepareLattice( olb::ThermalUnitConverter<SimType, NSDESCRIPTOR, TDESCRIPTOR> const& converter,
    olb::SuperLattice2D<SimType, NSDESCRIPTOR>& NSlattice,
    olb::SuperLattice2D<SimType, TDESCRIPTOR>& ADlattice,
    olb::Dynamics<SimType, NSDESCRIPTOR>& bulkDynamics,
    olb::Dynamics<SimType, TDESCRIPTOR>& advectionDiffusionBulkDynamics,
    olb::Dynamics<SimType, TDESCRIPTOR>& advectionDiffusionObstacleDynamics,
    olb::BounceBack<SimType, NSDESCRIPTOR>& bounceBackUpperWall,
    olb::BounceBack<SimType, NSDESCRIPTOR>& bounceBackLowerWall,
    olb::sOnLatticeBoundaryCondition2D<SimType,TDESCRIPTOR>& TboundaryCondition,
    olb::sOnLatticeBoundaryCondition2D<SimType,NSDESCRIPTOR>& NSboundaryCondition,
    SuperGeometry2D<SimType>& superGeometry,
    SimType TempInitial,
    SimType rhoVapor, SimType rhoLiquid )
    {

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

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

    clout << “NS omega = ” << omega << endl;
    clout << “T omega = ” << Tomega << endl;

    auto bulkIndicator = superGeometry.getMaterialIndicator({materialFluid, materialLowerWall, materialUpperWall, materialInflow, materialOutflow, materialObstacle});
    auto fluidIndicator = superGeometry.getMaterialIndicator({materialFluid, materialInflow, materialOutflow});
    auto wallIndicator = superGeometry.getMaterialIndicator({materialUpperWall, materialLowerWall});
    auto obstacleIndicator = superGeometry.getMaterialIndicator({materialObstacle});

    clout << “Define NS Dynamics …” << endl;
    NSlattice.defineDynamics( superGeometry, materialBase, &olb::instances::getNoDynamics<SimType, NSDESCRIPTOR>() );
    NSlattice.defineDynamics( fluidIndicator, &bulkDynamics );
    NSlattice.defineDynamics( superGeometry, materialUpperWall, &bounceBackUpperWall );
    NSlattice.defineDynamics( superGeometry, materialLowerWall, &bounceBackLowerWall );
    NSlattice.defineDynamics( superGeometry, materialObstacle, &olb::instances::getBounceBack<SimType, NSDESCRIPTOR>() );

    clout << “Define AD Dynamics …” << endl;
    ADlattice.defineDynamics( superGeometry, materialBase, &olb::instances::getNoDynamics<SimType, TDESCRIPTOR>() );
    ADlattice.defineDynamics( fluidIndicator, &advectionDiffusionBulkDynamics );
    ADlattice.defineDynamics( wallIndicator, &advectionDiffusionBulkDynamics );
    ADlattice.defineDynamics( obstacleIndicator, &advectionDiffusionObstacleDynamics );

    TboundaryCondition.addTemperatureBoundary(superGeometry, materialUpperWall, Tomega);
    TboundaryCondition.addTemperatureBoundary(superGeometry, materialLowerWall, Tomega);
    TboundaryCondition.addTemperatureBoundary(superGeometry, materialInflow, Tomega);
    TboundaryCondition.addTemperatureBoundary(superGeometry, materialOutflow, Tomega);

    clout << “Define Initial Conditions …” << endl;
    // Fluid Initial conditions
    olb::AnalyticalConst2D<SimType,SimType> zeroVelocity( 0,0 );

    AnalyticalConst2D<SimType,SimType> noiseVap( rhoVapor*0.1 );
    AnalyticalConst2D<SimType,SimType> noiseLiq( rhoLiquid*0.1 );
    AnalyticalRandom2D<SimType,SimType> random;

    SimType estimatedAdditionalAverageVapour = 2*(wallRho – rhoVapor)*conversionLength/ly0;
    clout << “estimatedAdditionalAverageVapour: ” << estimatedAdditionalAverageVapour << std::endl;
    // this estimate assumes that 2 full fluid nodes equal to the wall density are needed to compensate for the wall density gradient
    // This additional density is then spread over the entire y domain

    olb::AnalyticalConst2D<SimType,SimType> constRhoVapour( rhoVapor+estimatedAdditionalAverageVapour );
    olb::AnalyticalConst2D<SimType,SimType> constRhoLiquid( rhoLiquid * rhoLiquidMultiplier );
    olb::AnalyticalConst2D<SimType,SimType> constRhoUpperBC( rhoVapor );
    olb::AnalyticalConst2D<SimType,SimType> constRhoWall( rhoVapor );

    // Temperature Inital Conditions
    olb::AnalyticalConst2D<SimType,SimType> T_init( converter.getLatticeTemperature(TempInitial));

    // Random density distribution across inlet
    olb::AnalyticalIdentity2D<SimType,SimType> noiseIndicatorVap( random*noiseVap );
    olb::AnalyticalIdentity2D<SimType,SimType> noiseIndicatorLiq( random*noiseLiq );
    olb::AnalyticalIdentity2D<SimType,SimType> rhoBase( constRhoVapour );

    olb::SmoothIndicatorCircle2D<SimType,SimType> sphere( sphereOrigin, sphereRadius, sphereInterfaceWidth);

    olb::AnalyticalIdentity2D<SimType,SimType> sphereIndicator( sphere );

    olb::AnalyticalIdentity2D<SimType,SimType> rho3( rhoBase + noiseIndicatorVap + sphereIndicator*( constRhoLiquid – constRhoVapour + noiseIndicatorLiq) );

    //Initialize all values of distribution functions to their local equilibrium

    // Setting initial psi values
    olb::AnalyticalConst2D<SimType,SimType> initialPsi( 1.0 );
    NSlattice.defineField<olb::descriptors::PSI_PSEUDO_RHO>( fluidIndicator, initialPsi ) ;

    // Fluid init
    NSlattice.defineRhoU( superGeometry, materialFluid, rho3, zeroVelocity );
    NSlattice.iniEquilibrium( superGeometry, materialFluid, rho3, zeroVelocity );

    NSlattice.defineRhoU( superGeometry, materialUpperWall, constRhoUpperBC, zeroVelocity );
    NSlattice.iniEquilibrium( superGeometry, materialUpperWall, constRhoUpperBC, zeroVelocity );

    NSlattice.defineRhoU( superGeometry, materialLowerWall, constRhoWall, zeroVelocity );
    NSlattice.iniEquilibrium( superGeometry, materialLowerWall, constRhoWall, zeroVelocity );

    NSlattice.defineRhoU( superGeometry, materialInflow, rho3, zeroVelocity );
    NSlattice.iniEquilibrium( superGeometry, materialInflow, rho3, zeroVelocity );

    NSlattice.defineRhoU( superGeometry, materialOutflow, rho3, zeroVelocity );
    NSlattice.iniEquilibrium( superGeometry, materialOutflow, rho3, zeroVelocity );

    // Thermal init
    ADlattice.defineRho( superGeometry, materialFluid, T_init );
    ADlattice.iniEquilibrium( superGeometry, materialFluid, T_init, zeroVelocity );

    ADlattice.defineRho( superGeometry, materialUpperWall, T_init );
    ADlattice.iniEquilibrium( superGeometry, materialUpperWall, T_init, zeroVelocity );

    ADlattice.defineRho( superGeometry, materialLowerWall, T_init );
    ADlattice.iniEquilibrium( superGeometry, materialLowerWall, T_init, zeroVelocity );

    ADlattice.defineRho( superGeometry, materialInflow, T_init );
    ADlattice.iniEquilibrium( superGeometry, materialInflow, T_init, zeroVelocity );

    ADlattice.defineRho( superGeometry, materialOutflow, T_init );
    ADlattice.iniEquilibrium( superGeometry, materialOutflow, T_init, zeroVelocity );

    ADlattice.defineRho( obstacleIndicator, T_init );
    ADlattice.iniEquilibrium( obstacleIndicator, T_init, zeroVelocity );

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

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

    // Slowly guides wall density to wanted density
    void setBoundaryValues( olb::ThermalUnitConverter<SimType, NSDESCRIPTOR, TDESCRIPTOR> &converter,
    olb::SuperLattice2D<SimType, NSDESCRIPTOR>& NSlattice,
    olb::SuperLattice2D<SimType, TDESCRIPTOR>& ADlattice,
    int iT,
    SuperGeometry2D<SimType>& superGeometry,
    olb::BounceBackVariableRho<SimType, NSDESCRIPTOR> &bounceBackWall,
    SimType wallRhoMax, SimType rhoVapour,
    SimType newTempRatio )
    {
    olb::OstreamManager clout( std::cout,”setBoundaryValues” );

    if ( iT%iTsmoothUpdate==0 ) {

    SimType frac[1]= {};

    // if below iTmaxStart, get percentage ramp up, else 1
    if (iT<= rhoRampIter){
    // Smooth start curve, polynomial
    olb::PolynomialStartScale<SimType,int> startScale( rhoRampIter, SimType( 1 ) );
    int iTvec[1]= {iT};
    startScale( frac,iTvec );
    }
    else {
    frac[0] = 1;
    }
    bounceBackWall.setDensity( frac[0]*(wallRhoMax-rhoVapour) + rhoVapour );
    }

    }

    // Output to console and files
    void getResults( olb::SuperLattice2D<SimType, NSDESCRIPTOR>& NSlattice,
    olb::SuperLattice2D<SimType, TDESCRIPTOR>& ADlattice,
    olb::SuperGeometry2D<SimType>& superGeometry,
    olb::ThermalUnitConverter<SimType, NSDESCRIPTOR, TDESCRIPTOR> &converter,
    int iT,
    Timer<SimType>& timer,
    bool forcedSave=false)
    {
    olb::OstreamManager clout( std::cout,”getResults” );

    olb::SuperVTMwriter2D<SimType> vtmWriterNS( “fluidLatticeResults” );
    olb::SuperLatticePhysVelocity2D<SimType, NSDESCRIPTOR> velocity( NSlattice, converter );
    olb::SuperLatticeDensity2D<SimType, NSDESCRIPTOR> density( NSlattice );
    olb::SuperLatticeVelocity2D<SimType, NSDESCRIPTOR> latticeVelocity( NSlattice );
    olb::SuperLatticePhysPressure2D<SimType, NSDESCRIPTOR> pressure( NSlattice, converter );

    olb::SuperVTMwriter2D<SimType> vtmWriterThermal( “thermalLatticeResults” );
    olb::SuperLatticeDensity2D<SimType, TDESCRIPTOR> latticeTemperature( ADlattice );
    olb::SuperLatticePhysTemperature2D<SimType, NSDESCRIPTOR, TDESCRIPTOR> temperature(ADlattice, converter);

    vtmWriterNS.addFunctor( velocity );
    vtmWriterNS.addFunctor( latticeVelocity );
    vtmWriterNS.addFunctor( pressure );
    vtmWriterNS.addFunctor( density );
    vtmWriterThermal.addFunctor( latticeTemperature );
    vtmWriterThermal.addFunctor( temperature );

    const int vtkIter = converter.getLatticeTime( 2500*charTime/Nx/Nt );
    const int statIter = converter.getLatticeTime( 5000*charTime/Nx/Nt );
    const int ppmIter = converter.getLatticeTime( 2500*charTime/Nx/Nt );

    if ( iT==0 ) {
    // Writes the geometry, cuboid no. and rank no. as vti file for visualization
    olb::SuperLatticeGeometry2D<SimType, NSDESCRIPTOR> geometry( NSlattice, superGeometry );
    olb::SuperLatticeCuboid2D<SimType, NSDESCRIPTOR> cuboid( NSlattice );
    olb::SuperLatticeRank2D<SimType, NSDESCRIPTOR> rank( NSlattice );
    vtmWriterNS.write( geometry );
    vtmWriterNS.write( cuboid );
    vtmWriterNS.write( rank );
    vtmWriterNS.createMasterFile();
    vtmWriterThermal.createMasterFile();
    }

    // Writes the vtk files
    if ( iT%vtkIter==0 || forcedSave) {
    vtmWriterNS.write( iT );
    vtmWriterThermal.write( iT );
    }

    // Writes the ppm files
    if ( iT%ppmIter==0 || (iT<=50 && (iT%10==0) ) || forcedSave) {
    olb::BlockReduction2D2D<SimType> planeReductionRho( density );
    olb::BlockReduction2D2D<SimType> planeReductionT( temperature );
    olb::heatmap::plotParam<SimType> thermalPlotParam;
    thermalPlotParam.name = “physTemperature”;
    thermalPlotParam.minValue = converter.getCharPhysLowTemperature()-10;
    thermalPlotParam.maxValue = converter.getCharPhysHighTemperature()+10;
    // write output as JPEG
    olb::heatmap::write(planeReductionRho, iT);
    olb::heatmap::write(planeReductionT, iT, thermalPlotParam);
    }

    // Writes output on the console
    if ( iT%statIter==0 || forcedSave ) {
    // Timer console output
    timer.update( iT );
    timer.printStep();

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

    void writeGeometry( olb::SuperLattice2D<SimType, NSDESCRIPTOR>& NSlattice,
    olb::SuperGeometry2D<SimType>& superGeometry)
    {
    olb::SuperVTMwriter2D<SimType> vtmWriter( “d2law2D” );
    // Writes the geometry, cuboid no. and rank no. as vti file for visualization
    olb::SuperLatticeGeometry2D<SimType, NSDESCRIPTOR> geometry( NSlattice, superGeometry );
    olb::SuperLatticeCuboid2D<SimType, NSDESCRIPTOR> cuboid( NSlattice );
    olb::SuperLatticeRank2D<SimType, NSDESCRIPTOR> rank( NSlattice );
    vtmWriter.write( geometry );
    vtmWriter.write( cuboid );
    vtmWriter.write( rank );
    //vtmWriter.createMasterFile();
    }

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

    // === 1st Step: Initialization ===
    bool verbosity = false;
    olb::olbInit( &argc, &argv, verbosity );
    olb::OstreamManager clout( std::cout,”main” );

    // == Define simulations constants
    const SimType sigma = 0.35;
    const SimType eosA = 0.5/49.0;
    std::vector<SimType> rho0 ={1.0};
    const SimType G = -1;

    SimType molarMass = 18.0152;

    const SimType temperatureRatio = 0.7;
    const SimType newTempRatio = 0.8;

    SimType physThermalConductivity;
    for( physThermalConductivity = 100.0e-3; physThermalConductivity<=500.01e-3; physThermalConductivity+=200.0e-3 ){

    std::ostringstream outputDirStream;
    outputDirStream << std::setprecision(3);
    outputDirStream << “./results2D/”;
    mkdir( outputDirStream.str().c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH );
    outputDirStream << “lambda_” << physThermalConductivity << “/”;
    olb::singleton::directories().setOutputDir( outputDirStream.str() );

    // === 2nd Step: Prepare Geometry ===
    olb::Vector<SimType,2> extend( lx0, ly0 );
    olb::Vector<SimType,2> origin;
    olb::IndicatorCuboid2D<SimType> cuboid( extend, origin );

    // Instantiation of a cuboidGeometry with weights
    #ifdef PARALLEL_MODE_MPI
    const int noOfCuboids = olb::singleton::mpi().getSize();
    #else
    const int noOfCuboids = 1;
    #endif
    olb::CuboidGeometry2D<SimType> cuboidGeometry( cuboid, conversionLength, noOfCuboids );
    cuboidGeometry.setPeriodicity( false, false );

    // Instantiation of a loadBalancer
    olb::HeuristicLoadBalancer<SimType> loadBalancer( cuboidGeometry );

    // Instantiation of a superGeometry
    olb::SuperGeometry2D<SimType> superGeometry( cuboidGeometry, loadBalancer, 2 );

    prepareGeometry( superGeometry );

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

    // Pseudopotential multiphase single component dynamics
    // A dimensionless R is needed to recover correct latent heat
    SimType Rdimless = 8314.0/molarMass /(dx*dx/dt/dt) * TempCrit;

    olb::PengRobinsonExt<SimType,SimType> interactionPotential( G, 0.3443, eosA, 2./21., temperatureRatio, Rdimless );
    interactionPotential.write();

    SimType rho[2];
    olb::MaxwellConstruction<SimType, decltype(interactionPotential)> maxwellConstruction( interactionPotential );
    maxwellConstruction(rho, temperatureRatio); //temperatureRatio
    SimType densityLatticeVapour = rho[1];
    SimType densityLatticeLiquid = rho[0];
    clout << “Density liquid: ” << densityLatticeLiquid << std::endl;
    clout << “Density vapour: ” << densityLatticeVapour << std::endl;

    // == Define converter
    olb::ThermalUnitConverter<SimType, NSDESCRIPTOR, TDESCRIPTOR> converter(
    (SimType) dx, // physDeltaX: spacing between two lattice cells in __m__
    (SimType) dt, // physDeltaT = charLatticeVelocity / charPhysVelocity * dx
    (SimType) physLength, // charPhysLength: reference length of simulation geometry
    (SimType) maxVelocity, // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__
    (SimType) physViscosity, // physViscosity: physical kinematic viscosity in __m^2 / s__
    (SimType) liquidDensityRatio, // physDensity: physical density in __kg / m^3__
    (SimType) thermalConductivity, // physThermalConductivity
    (SimType) specificHeatCapacityVolume, // specificHeatCapacity
    (SimType) 0.0, // physThermalExpansionCoefficient
    (SimType) temperatureRatio*TempCrit, // charPhysLowTemperature
    (SimType) newTempRatio*TempCrit // charPhysHighTemperature
    );
    // Prints the converter log as console output
    converter.print();
    converter.write();

    olb::FluidProperties<SimType> fluidProperties(
    densityLatticeVapour, densityLatticeLiquid,
    converter.getLatticeSpecificHeatCapacity( specificHeatCapacityVolume ),
    converter.getLatticeSpecificHeatCapacity( specificHeatCapacityVolume ),
    converter.getLatticeThermalConductivity( physThermalConductivity ),
    converter.getLatticeThermalConductivity( physThermalConductivity ),
    converter.getCharPhysLowTemperature(), converter.getCharPhysHighTemperature(), TempCrit );

    fluidProperties.write();

    // Modified Guo forcing dynamics
    olb::ForcedModifiedGuoBGKdynamics<SimType, NSDESCRIPTOR> fluidSimulationDynamics (
    converter.getLatticeRelaxationFrequency(),
    olb::instances::getBulkMomenta<SimType,NSDESCRIPTOR>(), sigma
    );

    // Modified Advection Diffusion dynamics
    olb::Dynamics<SimType, TDESCRIPTOR>* thermalSimulationDynamics;
    olb::Dynamics<SimType, TDESCRIPTOR>* obstacleThermalSimulationDynamics;
    #ifdef ThermalMRT
    thermalSimulationDynamics = new olb::PseudopotentialAdvectionDiffusionMRTdynamics<SimType, TDESCRIPTOR>(
    converter.getLatticeThermalRelaxationFrequency(),
    olb::instances::getBulkMomenta<SimType,TDESCRIPTOR>()
    );

    obstacleThermalSimulationDynamics = new olb::PseudopotentialAdvectionDiffusionMRTdynamics<SimType, TDESCRIPTOR>(
    converter.getLatticeThermalRelaxationFrequency()-0.2,
    olb::instances::getBulkMomenta<SimType,TDESCRIPTOR>()
    );

    #else
    thermalSimulationDynamics = new olb::PseudopotentialAdvectionDiffusionBGKdynamics<SimType, TDESCRIPTOR>(
    converter.getLatticeThermalRelaxationFrequency(),
    olb::instances::getBulkMomenta<SimType,TDESCRIPTOR>()
    );
    #endif

    // boundary condition
    olb::sOnLatticeBoundaryCondition2D<SimType,TDESCRIPTOR> TboundaryCondition(ADlattice);
    olb::createAdvectionDiffusionBoundaryCondition2D<SimType,TDESCRIPTOR>(TboundaryCondition);

    // A bounce-back node with fictitious density
    // which is experienced by the partner fluid
    olb::BounceBackVariableRho<SimType, NSDESCRIPTOR> bounceBackUpperWall(densityLatticeVapour);

    // A bounce-back node with fictitious density
    // which is experienced by the partner fluid
    olb::BounceBackVariableRho<SimType, NSDESCRIPTOR> bounceBackLowerWall(densityLatticeVapour);

    olb::sOnLatticeBoundaryCondition2D<SimType,NSDESCRIPTOR> NSboundaryCondition(NSlattice);

    prepareLattice( converter,
    NSlattice, ADlattice,
    fluidSimulationDynamics, *thermalSimulationDynamics,
    *obstacleThermalSimulationDynamics,
    bounceBackUpperWall, bounceBackLowerWall,
    TboundaryCondition, NSboundaryCondition,
    superGeometry,
    TempCrit * temperatureRatio,
    densityLatticeVapour, densityLatticeLiquid );

    // === 4th Step: Prepare Lattice Coupling
    olb::LatticeCouplingGenerator2D<SimType, NSDESCRIPTOR> *thermalMultiphaseCoupling;
    #ifdef ThermalMRT
    thermalMultiphaseCoupling = new olb::CombinedShanChenThermalMRTCouplingGenerator2D<SimType, NSDESCRIPTOR>(
    0, converter.getLatticeLength(lx0),
    0, converter.getLatticeLength(ly0),
    G, converter.getLatticeThermalRelaxationFrequency(),
    0, 0, {0,1},
    fluidProperties, interactionPotential);
    #else
    thermalMultiphaseCoupling = new olb::CombinedShanChenThermalCouplingGenerator2D<SimType, NSDESCRIPTOR>(
    0, converter.getLatticeLength(lx0),
    0, converter.getLatticeLength(ly0),
    G, converter.getLatticeThermalRelaxationFrequency(),
    0, 0, {0,1},
    fluidProperties, interactionPotential);
    #endif

    NSlattice.addLatticeCoupling(*thermalMultiphaseCoupling, ADlattice);

    // === 5th Step: Main Loop with Timer ===
    clout << “starting simulation…” << endl;
    olb::util::Timer<SimType> timer( converter.getLatticeTime( endTimeEvap ), superGeometry.getStatistics().getNvoxel() );
    timer.start();
    int latticeInitTime = converter.getLatticeTime( endTimeInit );

    //setBoundaryValues( converter, NSlattice, ADlattice, 0, superGeometry, temperatureRatio);
    getResults(NSlattice, ADlattice, superGeometry, converter, 0, timer );

    olb::AnalyticalConst2D<SimType,SimType> zeroVelocity( 0,0 );
    olb::AnalyticalConst2D<SimType,SimType> T_const( converter.getLatticeTemperature(TempCrit * temperatureRatio));

    // ******************the first loop******************
    for ( int iT = 0; iT < latticeInitTime; ++iT ){

    // === 6th Step: Definition of Initial and Boundary Conditions ===
    setBoundaryValues( converter, NSlattice, ADlattice, iT, superGeometry, bounceBackLowerWall, wallRho, densityLatticeVapour, temperatureRatio);
    // === 7th Step: Collide and Stream Execution ===

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

    // === Reset temperature field
    ADlattice.defineRho(superGeometry, materialFluid, T_const);
    ADlattice.iniEquilibrium(superGeometry, materialFluid, T_const, zeroVelocity);

    ADlattice.defineRho(superGeometry, materialInflow, T_const);
    ADlattice.iniEquilibrium(superGeometry, materialInflow, T_const, zeroVelocity);

    ADlattice.defineRho(superGeometry, materialOutflow, T_const);
    ADlattice.iniEquilibrium(superGeometry, materialOutflow, T_const, zeroVelocity);

    // === 8th Step: Computation and Output of the Results ===
    getResults(NSlattice, ADlattice, superGeometry, converter, iT+1, timer );
    if ( std::isnan( NSlattice.getStatistics().getAverageRho() ) ) {
    clout << “Terminating this run at step ” << iT+1 << ” => density is NAN” << std::endl;
    getResults(NSlattice, ADlattice, superGeometry, converter, iT+1, timer, true );
    break;
    }
    }

    olb::AnalyticalConst2D<SimType,SimType> T_new( converter.getLatticeTemperature(TempCrit * newTempRatio));
    olb::AnalyticalConst2D<SimType,SimType> T_newHigher( converter.getLatticeTemperature(TempCrit * 0.9));

    ADlattice.defineRho(superGeometry, materialUpperWall, T_new);
    ADlattice.defineRho(superGeometry, materialLowerWall, T_new);
    ADlattice.defineRho(superGeometry, materialInflow, T_newHigher);
    ADlattice.defineRho(superGeometry, materialOutflow, T_new);

    // ******************the second loop******************
    for ( int iT = latticeInitTime; iT < converter.getLatticeTime( endTimeEvap ); ++iT ) {

    // === 6th Step: Definition of Initial and Boundary Conditions ===

    // === 7th Step: Collide and Stream Execution ===

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

    // === 8th Step: Computation and Output of the Results ===
    getResults(NSlattice, ADlattice, superGeometry, converter, iT+1, timer );
    if ( std::isnan( NSlattice.getStatistics().getAverageRho() ) ) {
    clout << “Terminating this run at step ” << iT+1 << ” => density is NAN” << std::endl;
    getResults(NSlattice, ADlattice, superGeometry, converter, iT+1, timer, true );
    break;
    }
    }
    delete thermalSimulationDynamics; // free previous alloc to avoid memory leak
    delete thermalMultiphaseCoupling;

    timer.stop();
    timer.printSummary();
    } // end for loop
    }

    #5989
    wf_guoyl
    Participant

    Now let me recap the question, I define two different “omega” to simulate transfer of heat in fluids and solids(obstacle). That is “obstacleThermalSimulationDynamics” and “thermalSimulationDynamics”. The program can run. But heat transfer can not occur in the solid part (obstacle). If anyone could help me, I would be very grateful.

    #5990
    Adrian
    Keymaster

    I received your mail and was able to compile and run the cuboid case.

    Contrary to what you want to achieve (as I understand it) the advectionDiffusionObstacleDynamics are not actually applied to the entire obstacle cuboid on the AD lattice. The reason for this is that a material indicator for number 6 is used but this is only set on the outer boundary of the obstacle due to the default geometry cleaning. You can see this in the geometry VTK output.

    I removed the clean calls and the obstacle now takes part in the AD dynamics. Alternatively one could use the original analytical obstacle indicator instead of a material indicator.

    Temperature field

    • This reply was modified 11 months, 1 week ago by Adrian.
    #5994
    wf_guoyl
    Participant

    Hello, Adrian

    Thank you very much for your answer! I delete “superGeometry.clean()” and “superGeometry.innerClean()”, and the program can run very well! But I don’t understand where the “original analytical obstacle indicator” is. Do you mean “SmoothIndicatorCuboid2D”? Are there any examples of this in OpenLB?

    Thanks again for your kind help!

    #5996
    Adrian
    Keymaster

    Good to hear!

    I mean the indicator you use to set the material number of the obstacle:

    
    olb:: IndicatorCuboid2D<SimType> obstacle( extend, origin);
    superGeometry.rename( materialFluid, materialObstacle, obstacle);
    

    You could pass this indicator directly to defineDynamics. But it makes sense to maintain a 1:1 relation between material and dynamics here, so no need to change it if it works.

    #5997
    wf_guoyl
    Participant

    Thank you for your quick reply!

    OpenLB is an excellent piece of software, thanks again to you and all the developers!

    Best wishes!

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