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 4 years 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

[dirk.neethling]

Dear Savis,

I am trying to simulate a film flow down an inclined (up to vertical) plane.
I think the Free Energy model using 2 components is the best approach, since I can then also model the surface tension and the fluid/air interface.
Using the multiComponent/microFluidics2d example as starting point, I then tried to reduce the domain to an inflow with a given velocity and an outflow with pressure balance. In the domain there would intially be air (component c2) with some of it being filled up the inflowing fluid (component c1).
Once I removed lattice3 from the example, i.e. remove the third component, thus reducing the flow to a binary case, nothing happens anymore, i.e. the inflowing fluid does not advance through the domain. I find this puzzling since the inflow boundary condition would surely force the fluid to move through the domain, right?
Using the approach in the multiComponent/youngLaplace2d example, with

rho = 1 = c1 + c2, => c1 = 1 – c2
and given that phi = c1 – c2 (see §4.4.2.1. Bulk Free Energy Model in the User Guide)
=> phi = (1 – c2) – c2 = 1 – c2 – c2, with psi = 0

However, as I mentioned the fluid is not moving through the domain, compared to the 3 component solution where the confluence of 2 or more inflows forces the column to move to the outlet.
Could you please indicate to me where my thinking is wrong?
Regards, Dirk

• This reply was modified 2 months ago by dirk.neethling.

[stephan]

Dear Dirk,

thank you for posting.
Please also adapt everything else in your modified example to the 2-component case.
Then, the boundary conditions have to made consistent with this reduction.
I rather expect that you have to modify several other sections of the .cpp file to make this work properly. However, you can do this (as you already have started) with looking at the other examples with 2 components.

BR
Stephan

[Author]

Posts