Reply To: Thermal Poiseuille 2D.
OpenLB – Open Source Lattice Boltzmann Code › Forums › on OpenLB › General Topics › Thermal Poiseuille 2D. › Reply To: Thermal Poiseuille 2D.
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***************************
***************************************************************