[ramirofreile]

Hi Fedor,

Thank you for your reply on this.
I am doing exactly what you recommend, keeping the lattice u at 0.1 and adding an extra tau_advection as a function of tau_fluid and turbulent Pr number. I have also tried refining the mesh of course.

I feel like the double distribution formulation hits a limit at a sufficiently high Reynolds number. For example, at a Re 10,000 it is very hard to simulate Prandtl numbers around 10. The thermal field gets too unstable.

Thank you,

Ramiro Freile

[ramirofreile]

Dear Stephan,

I have gone over the implementation of the wall functions in OpenLB and the succesful application of this postprocessor for cylindrical geometries is not very straightforward to me.

Even though there is a function to calculate wall distances with respect to an indicator, there are some sections of the code that seem more suitable to rectangular geometries, like the van driest calculation of Tau_eff (which uses a wall distance of 0.5 if none is provided). Also, some velocity and gradient reconstructions seem to be suitable for rectangular geometries, where the normal of the wall is aligned with the lattice nodes. That is why I asked this particular question.

If I am completely wrong in my analysis, please feel free to tell me so.

Thank you,

Ramiro

[ramirofreile]

Hi Stephan,

If you take a look at SuperLatticeTimeAveraged3D.hh, in line 64 (current Doxygen), you will see that the expression matches the expression I wrote in my first bug report message.

Mathematically, the cross correlation of two fields ends up being equal to

CrossCorr = (uv)_av – u_av v_av ,

where _av denotes a time average. Each time average is the sum of N observations divided by _ensembles.
Thus, the first term requires a single division by _ensembles and the second term requires a double division by _ensembles.

In my simulations I notices that the Reynolds stresses never converged and kept increasing in value.

Thank you,

Ramiro Freile

[ramirofreile]

Hi Adrian,

Thanks for your reply.
I am working in the 1.4 version, mainly because I have several custom functions for my problem already in the 1.4 format, which would require some adaptation to the new version.

I have a thermal turbulent application, and I was looking for a periodic BC in the first portion of my domain, so that the flow develops. Basically a periodic BC between the pipe inlet and a pipe slice of the domain at a certain distance (both planes would have the same normal).

So far I implemented PostProcessor which attempts a periodic BC between these two slices, using two material indicators. The problem is that when I identify the corresponding lattice coordinates of the node (upstream or downstream), I am unable to access the cell data.

I am currently looking at blockreduction functors as you suggested to see if I can solve the issue.

Thank you,

Ramiro Freile

[ramirofreile]

Hello Anas,

There is a coupling for the total enthalpy in 3D. I suggest you always look for classes in the Doxygen search engine:

https://www.openlb.net/DoxyGen/html/index.html

The conversion should be similar than how you would convert some of the examples in the laminar folder, e.g. cavity2D->3D, poiseuille 2D->3D.

Ramiro Freile

[ramirofreile]

Hello Antonio,

I have been trying to implement a Neumann boundary condition for the AD equation as well. Is there any chance you could share your way of implementing it to have it as a benchmark?
I would greatly appreciate it,

Thank you,

Ramiro

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

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

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

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]

Adrian,

I greatly appreciate your help.
I could solve the issue following your suggestions. I would just like to mention that in getResults I added NSlattice.communicate() and also ADlattice.communicate() for it to output the results properly.
Thank you very much,

Ramiro

[ramirofreile]

Hello,

Just to add more evidence about the issue.
Here I attach an image of “rayleighbernard2d” sample case.
I also have the issue while viewing the results in Paraview. It feels as if two contiguous blocks overlap a distance of deltax.
https://ibb.co/wCGrRPW

Thank you,

Ramiro Freile

[ramirofreile]

Hello Nicolas,

Thank you for your answer.
The image of the previous post is the temperature in a Poiseuille flow, and the way I am outputting the temperature is the following :

SuperLatticeField3D<T, TDESCRIPTOR, TEMPERATURE> temperature(ADlattice);
 temperature.getName() = "temperature";
 vtmWriter.addFunctor( temperature );

To add clarity to the issue, I ran a case that is only divided into two subblocks. https://ibb.co/sqC8Jpc
https://ibb.co/RbdgF1z
https://ibb.co/Jy0wHLh

In one figure, a slice of the pipe is shown. You can see that at the midlength of the pipe, the temperature has a strange output, so does the enthalpy. The velocity field and pressure do not have any issue. In the remaining two figures I show each of the subblocks separately. You can still see the issue in each block at the merging point.

I tried doing Resample to Image but I still have issues when I do line plots at the centerline of the pipe and on the slices.
Thank you in advance,

Ramiro Freile

[ramirofreile]

Julius,

Thank you very much for the discussion. To finish understanding your point I made a simple sketch which I attach. https://ibb.co/gmv52zY
In the figure, the orange thick line represents a geometry indicator for the boundary. We can say that the red nodes which coincide with the orange line are at the wall boundary. The blue nodes are fluid nodes and the black nodes are virtual nodes.

That said, for the Temperature lattice, after collision and streaming, the perpendicular population entering the fluid domain is set accordingly to impose the Dirichlet temperature there. However for the velocity, the halfway or the fullway bounce back would set the velocity to zero, deltax/2 in the upwards direction of the figure.

What I can say though, is that looking at the results of some examples in OpenLB, at those boundary nodes, the velocity is zero and the temperature is set as it is supposed to (rayleighBernard2d example). I am just having a struggle figuring out how it is done.
Thanks,

Ramiro

[ramirofreile]

Julius,

The fact that in AdvectionDiffusionBoundariesDynamics for D3Q7 we set the missing populations by T difference = dirichletTemperature - (T) 1 - sum; makes me think we are at a wet node. Also, for D3Q19 a bounce back of the non-equilibrium populations is done, like the Zou-He wet node BC approach.

While it is true that the supergeometry is created before the boundary conditions and the nodes locations are not influenced by this step, if there is a deltax/2 of difference between lattice node locations the coupling between the lattices could be affected.

I may be missing some details of the code..

Thank you,
Ramiro

• This reply was modified 2 years, 7 months ago by ramirofreile.

[Author]

Posts