Thermal Poiseuille 2D.

[ramirofreile]

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

[FBukreev]

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

[ramirofreile]

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

[ramirofreile]

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 2 years, 2 months ago by ramirofreile.

[ramirofreile]

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.

[stephan]

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

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

[stephan]

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

[Author]

Posts