Multiphase flow in 3D using FreeEnergy LBM

[mehrdadvasheghanifarahani]

Hi OpenLB community,

I’m attempting to simulate the process of the multiphase fluid flow in a cubic domain using Free Energy LBM. I’ve successfully simulated the process in 2D; the code works nicely without any warning, error, nan value, etc. However, when attempting to convert the code from 2D to 3D, I had a number of issues.

The first one was a compiler warning related to uninitialization of iC in Coupling objects. To address this issue, I had a look at the other examples implementing Free Energy LB in 3D and noticed that it could be suppressed by using the “shift” function (please see contact angle 3D example).
However, I still have issues when running the code in 3D. That is, after a number of time steps, I get nan values for average rho and average energy and literally the solution diverges. I checked the results and noticed that phi (which is C1 – C2) monotonously increases at the outlet dirichlet boundary and even becomes higher than 1 (which is the constrain for it). I was wondering if anybody could suggest any solution for this issue. I’m afraid that suppressing the compiler warning might cause this issue.

Kind Regards,
Mehrdad

• This topic was modified 3 years, 7 months ago by mehrdadvasheghanifarahani.

[savis]

Hi Mehrdad,

Would it be possible for you to share your code so that I can take a look to see what may be the cause? There are situations where the outlet boundary condition can become unstable, so that may be the case. However, if it is working in 2D then there may just be an error somewhere. Also, the outlet is not a Dirichlet boundary for phi because it must allow different fluids to pass, so some change is to be expected.

Regarding the uninitialized iC warning, I don’t believe this will be the cause of the issue. The warning seems to be because LatticeCouplingGenerator3D has an iC attribute that is uninitialized in the constructor, whereas the 2D version does not. However, I don’t fully understand the lattice coupling interface, so I simply added a shift function to the example. It would be great if someone else with a better understanding could take a look at this.

Best,
Sam

[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

[savis]

I have found the issue, it is a bug in this function. All of the “uPerp = +/- u[0/1/2];” lines have the wrong sign, which causes problems whenever an interface tries to pass the outlet.

Thanks for pointing out this issue. I will fix it for the next release. In the mean time you can also fix it in your own installation.

Best,
Sam

[Author]

Posts