Reply To: Thermal Poiseuille 2D.

[ramirofreile]

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***************************
***************************************************************