[mehrdadvasheghanifarahani]

Hi Michael,

As far as I know, the relaxation time of the hydrodynamics lattice for the fluids with two different viscosities can be accounted using this equation: 1 / tau_f = c_1 / tau_1 + c_2 / tau_2.

Would you please advise how the dynamics can be implemented using variable relaxation time? Is there any specific example that I can look into to learn how to implement this? Thanks.

Kind regards,
Mehrdad

[mehrdadvasheghanifarahani]

Lovely, thanks very much, Mathias.

Kind Regards,
Mehrdad

[mehrdadvasheghanifarahani]

Hi Mathias,

Thanks very much for your reply.
I checked the OpenLB examples; I think the freeslip boundary implemented in poiseuilleXd examples must give the symmetry boundary, is that correct?

Kind Regards,
Mehrdad

[mehrdadvasheghanifarahani]

Hi Mathias,

Thanks very much for clarifying this.

Kind Regards,
Mehrdad

[mehrdadvasheghanifarahani]

Dear Sam,

Many thanks for your response. Sure, I’d be happy to share the code:

#include “olb3D.h”
#include “olb3D.hh” // use only generic version!

#include <cstdlib>
#include <iostream>
#include <fstream>

using namespace olb;
using namespace olb::descriptors;
using namespace olb::graphics;
using namespace std;

typedef double T;
#define DESCRIPTOR D3Q19<CHEM_POTENTIAL,FORCE>

// Parameters for the simulation setup
const int N = 100;
const T nx = 20.;
const T ny = 100.;
const T nz = 100.;
const T dx = ny / N;

const T inlet1Velocity = 0.00056; // [lattice units]
const T outletDensity = 1.; // [lattice units]
const T alpha = 1.; // Interfacial width [lattice units]

const T kappa1 = 0.001;
const T kappa2 = 0.001;
const T kappa3 = 0.001;

const T gama = 1.; // For mobility of interfaces [lattice units]
const T h1 = 0.0;
const T h2 = 0.0;
const T h3 = 0.0;

const int maxIter = 1000000;
const int vtkIter = 1000;
const int statIter = 2000;

void prepareGeometry( UnitConverter<T,DESCRIPTOR> const& converter,
SuperGeometry3D<T>& superGeometry ) {
OstreamManager clout( std::cout,”prepareGeometry” );
clout << “Prepare Geometry …” << std::endl;

superGeometry.rename(0, 2);
superGeometry.rename(2, 1, 1, 1, 1);

// Inlets and outlet

Vector<T,3> origin = superGeometry.getStatistics().getMinPhysR( 2 );
origin[1] += converter.getConversionFactorLength()/2.;
origin[2] += converter.getConversionFactorLength()/2.;

Vector<T,3> extend = superGeometry.getStatistics().getMaxPhysR( 2 );
extend[1] = extend[1]-origin[1]-converter.getConversionFactorLength()/2.;
extend[2] = extend[2]-origin[2]-converter.getConversionFactorLength()/2.;

// Set material number for inflow
origin[0] = superGeometry.getStatistics().getMinPhysR( 2 )[0]-converter.getConversionFactorLength();
extend[0] = 2*converter.getConversionFactorLength();
IndicatorCuboid3D<T> inflow( extend,origin );
superGeometry.rename( 2,3,1,inflow );

// Set material number for outflow
origin[0] = superGeometry.getStatistics().getMaxPhysR( 2 )[0]-converter.getConversionFactorLength();
extend[0] = 2*converter.getConversionFactorLength();
IndicatorCuboid3D<T> outflow( extend,origin );
superGeometry.rename( 2,8,1,outflow );

superGeometry.innerClean();
superGeometry.checkForErrors();
superGeometry.print();

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

void prepareLattice( SuperLattice3D<T, DESCRIPTOR>& sLattice1,
SuperLattice3D<T, DESCRIPTOR>& sLattice2,
SuperLattice3D<T, DESCRIPTOR>& sLattice3,
Dynamics<T, DESCRIPTOR>& bulkDynamics1,
Dynamics<T, DESCRIPTOR>& bulkDynamics2,
Dynamics<T, DESCRIPTOR>& bulkDynamics3,
sOnLatticeBoundaryCondition3D<T,DESCRIPTOR>& sOnBC1,
sOnLatticeBoundaryCondition3D<T,DESCRIPTOR>& sOnBC2,
sOnLatticeBoundaryCondition3D<T,DESCRIPTOR>& sOnBC3,
UnitConverter<T, DESCRIPTOR>& converter,
SuperGeometry3D<T>& superGeometry ) {

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

// define lattice dynamics
sLattice1.defineDynamics( superGeometry, 0, &instances::getNoDynamics<T, DESCRIPTOR>() );
sLattice2.defineDynamics( superGeometry, 0, &instances::getNoDynamics<T, DESCRIPTOR>() );
sLattice3.defineDynamics( superGeometry, 0, &instances::getNoDynamics<T, DESCRIPTOR>() );

sLattice1.defineDynamics( superGeometry, 1, &bulkDynamics1 );
sLattice2.defineDynamics( superGeometry, 1, &bulkDynamics2 );
sLattice3.defineDynamics( superGeometry, 1, &bulkDynamics3 );

sLattice1.defineDynamics( superGeometry, 2, &instances::getNoDynamics<T, DESCRIPTOR>() );
sLattice2.defineDynamics( superGeometry, 2, &instances::getNoDynamics<T, DESCRIPTOR>() );
sLattice3.defineDynamics( superGeometry, 2, &instances::getNoDynamics<T, DESCRIPTOR>() );

sLattice1.defineDynamics( superGeometry, 3, &instances::getNoDynamics<T, DESCRIPTOR>() );
sLattice2.defineDynamics( superGeometry, 3, &instances::getNoDynamics<T, DESCRIPTOR>() );
sLattice3.defineDynamics( superGeometry, 3, &instances::getNoDynamics<T, DESCRIPTOR>() );

sLattice1.defineDynamics( superGeometry, 8, &instances::getNoDynamics<T, DESCRIPTOR>() );
sLattice2.defineDynamics( superGeometry, 8, &instances::getNoDynamics<T, DESCRIPTOR>() );
sLattice3.defineDynamics( superGeometry, 8, &instances::getNoDynamics<T, DESCRIPTOR>() );

// add wall boundary
sOnBC1.addFreeEnergyWallBoundary( superGeometry, 2, alpha, kappa1, kappa2, kappa3, h1, h2, h3, 1 );
sOnBC2.addFreeEnergyWallBoundary( superGeometry, 2, alpha, kappa1, kappa2, kappa3, h1, h2, h3, 2 );
sOnBC3.addFreeEnergyWallBoundary( superGeometry, 2, alpha, kappa1, kappa2, kappa3, h1, h2, h3, 3 );

// add inlet boundaries
T omega = converter.getLatticeRelaxationFrequency();
auto inlet1Indicator = superGeometry.getMaterialIndicator(3);
sOnBC1.addFreeEnergyInletBoundary( inlet1Indicator, omega, “velocity”, 1 );
sOnBC2.addFreeEnergyInletBoundary( inlet1Indicator, omega, “velocity”, 2 );
sOnBC3.addFreeEnergyInletBoundary( inlet1Indicator, omega, “velocity”, 3 );

// add outlet boundary
auto outletIndicator = superGeometry.getMaterialIndicator(8);
sOnBC1.addFreeEnergyOutletBoundary( outletIndicator, omega, “density”, 1 );
sOnBC2.addFreeEnergyOutletBoundary( outletIndicator, omega, “density”, 2 );
sOnBC3.addFreeEnergyOutletBoundary( outletIndicator, omega, “density”, 3 );

// bulk initial conditions
std::vector<T> v( 3,T() );
AnalyticalConst3D<T,T> zeroVelocity( v );

AnalyticalConst3D<T,T> zero ( 0. );
AnalyticalConst3D<T,T> one ( 1. );
SmoothIndicatorCuboid3D<T,T> section1( {nx * 0.25, ny/2., nz/2.}, nx * 0.5, ny, nz, 0.);
SmoothIndicatorCuboid3D<T,T> section2( {nx * 0.65, ny/2., nz/2.}, nx * 0.3, ny, nz, 0.);

AnalyticalIdentity3D<T,T> c1( section1 );
AnalyticalIdentity3D<T,T> c2( section2 );
AnalyticalIdentity3D<T,T> rho( one );
AnalyticalIdentity3D<T,T> phi( c1 – c2 );
AnalyticalIdentity3D<T,T> psi( rho – c1 – c2 );

auto allIndicator = superGeometry.getMaterialIndicator({1, 2, 3});
sLattice1.iniEquilibrium( allIndicator, rho, zeroVelocity );
sLattice2.iniEquilibrium( allIndicator, phi, zeroVelocity );
sLattice3.iniEquilibrium( allIndicator, psi, zeroVelocity );

// inlet boundary conditions

std::vector<T> maxVelocity( 3,0 );
maxVelocity[0] = 1.5*inlet1Velocity;

T distance2Wall = converter.getConversionFactorLength()/2.;
RectanglePoiseuille3D<T> poiseuilleU( superGeometry, 3, maxVelocity, distance2Wall, distance2Wall, distance2Wall );

sLattice1.defineU( inlet1Indicator, poiseuilleU);
sLattice2.defineRho( inlet1Indicator, phi);
sLattice3.defineRho( inlet1Indicator, psi);

// outlet initial / boundary conditions
AnalyticalConst3D<T,T> rhoOutlet( outletDensity );
AnalyticalIdentity3D<T,T> phiOutlet( zero );
AnalyticalIdentity3D<T,T> psiOutlet( rhoOutlet );

sLattice1.defineRho( outletIndicator, rhoOutlet );
sLattice2.defineRho( outletIndicator, phiOutlet );
sLattice3.defineRho( outletIndicator, psiOutlet );

// initialise lattices
sLattice1.initialize();
sLattice2.initialize();
sLattice3.initialize();

sLattice1.communicate();
sLattice2.communicate();
sLattice3.communicate();

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

void prepareCoupling(SuperLattice3D<T, DESCRIPTOR>& sLattice1,
SuperLattice3D<T, DESCRIPTOR>& sLattice2,
SuperLattice3D<T, DESCRIPTOR>& sLattice3,
SuperGeometry3D<T>& superGeometry) {
OstreamManager clout( std::cout,”prepareCoupling” );
clout << “Add lattice coupling” << endl;

// Bulk couplings
FreeEnergyChemicalPotentialGenerator3D<T,DESCRIPTOR> coupling2( alpha, kappa1, kappa2, kappa3 );
FreeEnergyForceGenerator3D<T,DESCRIPTOR> coupling3;

// Inlet / outlet couplings
FreeEnergyDensityOutletGenerator3D<T,DESCRIPTOR> coupling1( outletDensity );
FreeEnergyInletOutletGenerator3D<T,DESCRIPTOR> coupling4;

// Suppress compiler warnings
coupling2.shift(0, 0, 0);
coupling3.shift(0, 0, 0);
coupling1.shift(0, 0, 0);
coupling4.shift(0, 0, 0);

// The DensityOutlet coupling must be added to the first lattice and come before the ChemicalPotential coupling
// The InletOutlet couplings must be added to the second lattice and come after the Force coupling.
sLattice1.addLatticeCoupling<DESCRIPTOR>( superGeometry, 8, coupling1, {&sLattice2, &sLattice3} );

sLattice1.addLatticeCoupling<DESCRIPTOR>( superGeometry, 1, coupling2, {&sLattice2, &sLattice3} );
sLattice2.addLatticeCoupling<DESCRIPTOR>( superGeometry, 1, coupling3, {&sLattice1, &sLattice3} );

sLattice2.addLatticeCoupling<DESCRIPTOR>( superGeometry, 3, coupling4, {&sLattice1, &sLattice3} );

sLattice2.addLatticeCoupling<DESCRIPTOR>( superGeometry, 8, coupling4, {&sLattice1, &sLattice3} );

clout << “Add lattice coupling … OK!” << endl;
}

void getResults( SuperLattice3D<T, DESCRIPTOR>& sLattice1,
SuperLattice3D<T, DESCRIPTOR>& sLattice2,
SuperLattice3D<T, DESCRIPTOR>& sLattice3, int iT,
SuperGeometry3D<T>& superGeometry, Timer<T>& timer,
UnitConverter<T, DESCRIPTOR> converter) {

OstreamManager clout( std::cout,”getResults” );
SuperVTMwriter3D<T> vtmWriter( “microFluidics2d” );

if ( iT==0 ) {
// Writes the geometry, cuboid no. and rank no. as vti file for visualization
SuperLatticeGeometry3D<T, DESCRIPTOR> geometry( sLattice1, superGeometry );
SuperLatticeCuboid3D<T, DESCRIPTOR> cuboid( sLattice1 );
SuperLatticeRank3D<T, DESCRIPTOR> rank( sLattice1 );
vtmWriter.write( geometry );
vtmWriter.write( cuboid );
vtmWriter.write( rank );
vtmWriter.createMasterFile();
}

// Get statistics
if ( iT%statIter==0 ) {
// Timer console output
timer.update( iT );
timer.printStep();
sLattice1.getStatistics().print( iT, converter.getPhysTime(iT) );
sLattice2.getStatistics().print( iT, converter.getPhysTime(iT) );
sLattice3.getStatistics().print( iT, converter.getPhysTime(iT) );
}

// Writes the VTK files
if ( iT%vtkIter==0 ) {
SuperLatticeVelocity3D<T, DESCRIPTOR> velocity( sLattice1 );
SuperLatticeDensity3D<T, DESCRIPTOR> density1( sLattice1 );
density1.getName() = “rho”;
SuperLatticeDensity3D<T, DESCRIPTOR> density2( sLattice2 );
density2.getName() = “phi”;
SuperLatticeDensity3D<T, DESCRIPTOR> density3( sLattice3 );
density3.getName() = “density-fluid-3”;

AnalyticalConst3D<T,T> half_( 0.5 );
SuperLatticeFfromAnalyticalF3D<T, DESCRIPTOR> half(half_, sLattice1);

SuperIdentity3D<T,T> c1 (half*(density1+density2-density3));
c1.getName() = “density-fluid-1”;
SuperIdentity3D<T,T> c2 (half*(density1-density2-density3));
c2.getName() = “density-fluid-2”;

vtmWriter.addFunctor( velocity );
vtmWriter.addFunctor( density1 );
vtmWriter.addFunctor( density2 );
vtmWriter.addFunctor( density3 );
vtmWriter.addFunctor( c1 );
vtmWriter.addFunctor( c2 );
vtmWriter.write( iT );
}
}

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

// === 1st Step: Initialization ===

olbInit( &argc, &argv );
singleton::directories().setOutputDir( “./tmp/” );
OstreamManager clout( std::cout,”main” );

UnitConverterFromResolutionAndRelaxationTime<T,DESCRIPTOR> converter(
(T) N, // resolution
(T) 1., // lattice relaxation time (tau)
(T) nz, // charPhysLength: reference length of simulation geometry
(T) 1.e-6, // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__
(T) 0.1, // physViscosity: physical kinematic viscosity in __m^2 / s__
(T) 1. // physDensity: physical density in __kg / m^3__
);

// Prints the converter log as console output
converter.print();

// === 2nd Step: Prepare Geometry ===
std::vector<T> extend = {nx, ny, nz};
std::vector<T> origin = {0, 0, 0};
IndicatorCuboid3D<T> cuboid(extend,origin);
#ifdef PARALLEL_MODE_MPI
CuboidGeometry3D<T> cGeometry( cuboid, converter.getPhysDeltaX(), singleton::mpi().getSize() );
#else
CuboidGeometry3D<T> cGeometry( cuboid, converter.getPhysDeltaX() );
#endif

// Instantiation of loadbalancer
HeuristicLoadBalancer<T> loadBalancer( cGeometry );
loadBalancer.print();

// Instantiation of superGeometry
SuperGeometry3D<T> superGeometry( cGeometry,loadBalancer );

prepareGeometry( converter, superGeometry );

// === 3rd Step: Prepare Lattice ===
SuperLattice3D<T, DESCRIPTOR> sLattice1( superGeometry );
SuperLattice3D<T, DESCRIPTOR> sLattice2( superGeometry );
SuperLattice3D<T, DESCRIPTOR> sLattice3( superGeometry );

ForcedBGKdynamics<T, DESCRIPTOR> bulkDynamics1 (
converter.getLatticeRelaxationFrequency(),
instances::getBulkMomenta<T,DESCRIPTOR>() );

FreeEnergyBGKdynamics<T, DESCRIPTOR> bulkDynamics23 (
converter.getLatticeRelaxationFrequency(), gama,
instances::getBulkMomenta<T,DESCRIPTOR>() );

sOnLatticeBoundaryCondition3D<T, DESCRIPTOR> sOnBC1( sLattice1 );
sOnLatticeBoundaryCondition3D<T, DESCRIPTOR> sOnBC2( sLattice2 );
sOnLatticeBoundaryCondition3D<T, DESCRIPTOR> sOnBC3( sLattice3 );
createLocalBoundaryCondition3D<T, DESCRIPTOR> (sOnBC1);
createLocalBoundaryCondition3D<T, DESCRIPTOR> (sOnBC2);
createLocalBoundaryCondition3D<T, DESCRIPTOR> (sOnBC3);

prepareLattice( sLattice1, sLattice2, sLattice3, bulkDynamics1, bulkDynamics23,
bulkDynamics23, sOnBC1, sOnBC2, sOnBC3, converter, superGeometry );

prepareCoupling( sLattice1, sLattice2, sLattice3, superGeometry );

SuperExternal3D<T,DESCRIPTOR,CHEM_POTENTIAL> sExternal1 (superGeometry, sLattice1, sLattice1.getOverlap() );
SuperExternal3D<T,DESCRIPTOR,CHEM_POTENTIAL> sExternal2 (superGeometry, sLattice2, sLattice2.getOverlap() );
SuperExternal3D<T,DESCRIPTOR,CHEM_POTENTIAL> sExternal3 (superGeometry, sLattice3, sLattice3.getOverlap() );

// === 4th Step: Main Loop with Timer ===
int iT = 0;
clout << “starting simulation…” << endl;
Timer<T> timer( maxIter, superGeometry.getStatistics().getNvoxel() );
timer.start();

for ( iT=0; iT<maxIter; ++iT ) {
// Computation and output of the results
getResults( sLattice1, sLattice2, sLattice3, iT, superGeometry, timer, converter );

// Collide and stream execution
sLattice1.collideAndStream();
sLattice2.collideAndStream();
sLattice3.collideAndStream();

// MPI communication for lattice data
sLattice1.communicate();
sLattice2.communicate();
sLattice3.communicate();

// Execute coupling between the two lattices
sLattice1.executeCoupling();
sExternal1.communicate();
sExternal2.communicate();
sExternal3.communicate();
sLattice2.executeCoupling();
}

timer.stop();
timer.printSummary();

}

As I mentioned, my results show that “phi” is monotonously increasing in the outlet boundary and I believe this might cause the issue.

Kind Regards,
Mehrdad

[mehrdadvasheghanifarahani]

Hi Mathias,

Many thanks.
Yes, that’s a good idea. I’ll go for it.

Kind Regards,
Mehrdad

[mehrdadvasheghanifarahani]

Hi Mathias,

Many thanks for your kind response.

We have addFreeEnergyWallBoundary function for Free-energy-based LBM. For SC LBM, the interaction between the fluids is captured by Cohesive and Repulsive forces. It would be great if we could consider the interactions between the solid boundary and the fluids as well (I mean Adhesive force). Then, we would be able to see the phenomenon.

Kind Regards,
Mehrdad

[Author]

Posts