Identification of Geometric Model Errors in Multiphase Flow Simulations

[Bobbie]

Hello Everyone.
When executing the code for a two-dimensional multiphase flow simulation using microfluidics2d, I noticed that for geometric models composed of multiple structures (e.g., the model in bstep2d), the removed sections of the geometry are not properly considered in the multiphase flow simulation. Contrary to the expected behavior where these regions should be treated as solid domains, they are instead interpreted as part of one of the fluid phases. This discrepancy deviates significantly from the intended design. What could be the underlying cause of this issue, and what modifications to the code are necessary to resolve it?

[mathias]

Did you check on the material numbers? Are they set correctly?

[Bobbie]

I have verified that the material number is correctly assigned. Despite attempts to set the material number of the removed geometric model to 0 or treat it as a solid region, the results still indicate that it exhibits fluid-phase characteristics.

[mathias]

Without knowing details, it is not possible to help you. In the spring school we explain step by step how to do it and also in the user guide is a section on it. Individial support we can also offer based on a support contract. If so, please contact the OpenLB team using the contact form.

[Bobbie]

The following complete simulation code combines the examples of cylinder2d and microFluidcs2d. However, despite explicitly setting the cylinder’s material identifier to 0 to designate it as a solid domain, the post-processing results demonstrate unexpected fluid-like behavior within the cylindrical region. Could you please advise on necessary code modifications to properly enforce solid phase characteristics?

#include “olb2D.h”
#include “olb2D.hh”

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

using T = FLOATING_POINT_TYPE;
typedef D2Q9<CHEM_POTENTIAL,FORCE> DESCRIPTOR;

const int N = 100;
const T nx = 800e-6;
const T ny = 100e-6;
const T dx = ny/N;

const T alpha = 1.;
const T kappa1 = 0.005;
const T kappa2 = 0.005;
const T gama = 10.;
const T h1 = 0.0001448;
const T h2 = -0.0001448;
const T inletVelocity = 0.01;
const T outletDensity = 1.;

const int maxIter = 100000;
const int vtkIter = 200;
const int statIter = 200;

void prepareGeometry( SuperGeometry<T,2>& superGeometry,
std::shared_ptr<IndicatorF2D<T>> section1
)
{

OstreamManager clout( std::cout,”prepareGeometry” );
clout << “Prepare Geometry …” << std::endl;

superGeometry.rename(0,2,section1);
superGeometry.rename(2,1,{1,1});
IndicatorCuboid2D<T> inlet(dx,ny,{0.0,ny/T(2)});
IndicatorCuboid2D<T> outlet(dx,ny,{nx,ny/T(2)});
superGeometry.rename( 2,3,1,inlet );
superGeometry.rename( 2,4,1,outlet );
IndicatorCircle2D<T> circle({100.0*dx,ny/T(2)}, 12.5*dx);
superGeometry.rename( 1,0,circle );

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

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

void prepareLattice( SuperLattice<T, DESCRIPTOR>& sLattice1,
SuperLattice<T, DESCRIPTOR>& sLattice2,
UnitConverter<T, DESCRIPTOR>& converter,
SuperGeometry<T,2>& superGeometry)
{

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

T omega = converter.getLatticeRelaxationFrequency();

clout << “Prepare Lattice: Define lattice dynamics …” << std::endl;
sLattice1.defineDynamics<ForcedBGKdynamics>(superGeometry, 1);
sLattice2.defineDynamics<FreeEnergyBGKdynamics>(superGeometry, 1);

clout << “Prepare Lattice: Add wall boundary …” << std::endl;
setFreeEnergyWallBoundary(sLattice1, superGeometry, 2, alpha, kappa1, kappa2, h1, h2, 1);
setFreeEnergyWallBoundary(sLattice2, superGeometry, 2, alpha, kappa1, kappa2, h1, h2, 2);

clout << “Prepare Lattice: Add inlet boundaries …” << std::endl;
auto inletIndicator = superGeometry.getMaterialIndicator(3);
setFreeEnergyInletBoundary<T,DESCRIPTOR>(sLattice1, omega, inletIndicator, “velocity”, 1);
setFreeEnergyInletBoundary<T,DESCRIPTOR>(sLattice2, omega, inletIndicator, “velocity”, 2);

clout << “Prepare Lattice: Add outlet boundary …” << std::endl;
auto outletIndicator = superGeometry.getMaterialIndicator(4);
setFreeEnergyOutletBoundary<T,DESCRIPTOR>(sLattice1, omega, outletIndicator, “density”,1);
setFreeEnergyOutletBoundary<T,DESCRIPTOR>(sLattice2, omega, outletIndicator, “density”,2);

clout << “Prepare Lattice: Bulk initial conditions …” << std::endl;
std::vector<T> v( 2,T() );
AnalyticalConst2D<T,T> zeroVelocity( v );

AnalyticalConst2D<T,T> one( 1.0 );
AnalyticalConst2D<T,T> zero( 0.0 );

IndicatorCuboid2D<T> ind1(nx-dx,ny,{dx+(nx-dx)/T(2),ny/T(2)});
SmoothIndicatorCuboid2D<T,T> section1( ind1, 0 );
AnalyticalIdentity2D<T,T> c1( section1 );

AnalyticalIdentity2D<T,T> rho( one ); // c2 + c1
AnalyticalIdentity2D<T,T> phi( one – c1 -c1 );// c2 – c1

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

clout << “Prepare Lattice: Inlet boundary conditions …” << std::endl;
AnalyticalConst2D<T,T> inletU(inletVelocity);
sLattice1.defineU( inletIndicator, inletU );
sLattice2.defineRho( inletIndicator, phi );

clout << “Prepare Lattice: Outlet initial / Boundary conditions …” << std::endl;
AnalyticalConst2D<T,T> rhoOutlet( outletDensity );
AnalyticalIdentity2D<T,T> phiOutlet( zero );
sLattice1.defineRho( outletIndicator, rhoOutlet );
sLattice2.defineRho( outletIndicator, phiOutlet );

sLattice1.setParameter<descriptors::OMEGA>( omega );
sLattice2.setParameter<descriptors::OMEGA>( omega );
sLattice2.setParameter<collision::FreeEnergy::GAMMA>( gama );

clout << “Prepare Lattice: Initialise lattices …” << std::endl;sLattice1.initialize();
sLattice1.initialize();
sLattice2.initialize();

clout << “Prepare Lattice: Communicate …” << std::endl;
sLattice1.communicate();
sLattice2.communicate();

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

void prepareCoupling( SuperLattice<T, DESCRIPTOR>& sLattice1,
SuperLattice<T, DESCRIPTOR>& sLattice2,
SuperGeometry<T,2>& superGeometry )
{

OstreamManager clout( std::cout,”prepareCoupling” );
clout << “Add lattice coupling” << std::endl;

// 主体耦合
FreeEnergyChemicalPotentialGenerator2D<T,DESCRIPTOR> coupling2( alpha, kappa1, kappa2 ); //化学势耦合→coupling2
FreeEnergyForceGenerator2D<T,DESCRIPTOR> coupling3; //外力耦合→coupling3
// 进出口耦合
FreeEnergyDensityOutletGenerator2D<T,DESCRIPTOR> coupling1( outletDensity ); //密度出口耦合→coupling1
FreeEnergyInletOutletGenerator2D<T,DESCRIPTOR> coupling4; //进出口耦合→coupling4

sLattice1.addLatticeCoupling( superGeometry, 4, coupling1, sLattice2 );//密度出口耦合

sLattice1.addLatticeCoupling( superGeometry, 1, coupling2, sLattice2 );//化学势耦合
sLattice2.addLatticeCoupling( superGeometry, 1, coupling3, sLattice1 );//外力耦合

sLattice2.addLatticeCoupling( superGeometry, 3, coupling4, sLattice1 );//进出口耦合
sLattice2.addLatticeCoupling( superGeometry, 4, coupling4, sLattice1 );//进出口耦合

{
auto& communicator = sLattice1.getCommunicator(stage::PostCoupling());
communicator.requestField<CHEM_POTENTIAL>();
communicator.requestOverlap(sLattice1.getOverlap());
communicator.exchangeRequests();
}
{
auto& communicator = sLattice2.getCommunicator(stage::PreCoupling());
communicator.requestField<CHEM_POTENTIAL>();
communicator.requestOverlap(sLattice2.getOverlap());
communicator.exchangeRequests();
}

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

void getResults( SuperLattice<T, DESCRIPTOR>& sLattice1,
SuperLattice<T, DESCRIPTOR>& sLattice2, int iT,
SuperGeometry<T,2>& superGeometry, util::Timer<T>& timer,
UnitConverter<T, DESCRIPTOR> converter )
{

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

if ( iT==0 ) {
// Writes the geometry, cuboid no. and rank no. as vti file for visualization
SuperLatticeGeometry2D<T, DESCRIPTOR> geometry( sLattice1, superGeometry );
SuperLatticeCuboid2D<T, DESCRIPTOR> cuboid( sLattice1 );
SuperLatticeRank2D<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) );
clout << “averageRhoFluidOne=” << sLattice1.getStatistics().getAverageRho();
clout << “; averageRhoFluidTwo=” << sLattice2.getStatistics().getAverageRho() << std::endl;
}

// Writes the VTK files
if ( iT%vtkIter==0 ) {
sLattice1.setProcessingContext(ProcessingContext::Evaluation);
sLattice2.setProcessingContext(ProcessingContext::Evaluation);
SuperLatticeVelocity2D<T, DESCRIPTOR> velocity( sLattice1 );
SuperLatticeDensity2D<T, DESCRIPTOR> rho( sLattice1 );
rho.getName() = “rho”;
SuperLatticeDensity2D<T, DESCRIPTOR> phi( sLattice2 );
phi.getName() = “phi”;

AnalyticalConst2D<T,T> half_( 0.5 );
SuperLatticeFfromAnalyticalF2D<T, DESCRIPTOR> half(half_, sLattice1);
SuperIdentity2D<T,T> c1 (half*(rho+phi));
c1.getName() = “density-fluid-1”;
SuperIdentity2D<T,T> c2 (half*(rho-phi));
c2.getName() = “density-fluid-2”;

vtmWriter.addFunctor( velocity );
vtmWriter.addFunctor( rho );
vtmWriter.addFunctor( phi );
vtmWriter.addFunctor( c1 );
vtmWriter.addFunctor( c2 );
vtmWriter.write( iT );

// BlockReduction2D2D interpolates the data of a SuperF2D functor in a given resolution.
// This is primarily used for exporting GIF images via BlockGifWriter.
// BlockReduction2D2D在给定的分辨率内插一个SuperF2D函子的数据。这主要用于通过BlockGifWriter导出GIF图像。
// olb::BlockReduction2D2D<T>::BlockReduction2D2D(FunctorPtr<SuperF2D<T>>&&f,int resolution = 600,
// BlockDataSyncMode mode=BlockDataSyncMode::ReduceAndBcast)
// BlockReduction2D2D所需参数：FunctorPtr<SuperF2D<T>>&&f, resolution, ReduceAndBcast/ReduceOnly/None
BlockReduction2D2D<T> planeReduction( c1, 600, BlockDataSyncMode::ReduceOnly );
// write output as JPEG
heatmap::write(planeReduction, iT);
}
}

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

// === 第一步：初始化 ===

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) ny, // charPhysLength: reference length of simulation geometry
(T) 0.1, // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__
(T) 1e-5, // physViscosity: physical kinematic viscosity in __m^2 / s__
(T) 1000. // physDensity: physical density in __kg / m^3__
);

// Prints the converter log as console output
converter.print();
converter.write(“FreeEnergy”);

// === 第二步：准备几何模型 ===
#ifdef PARALLEL_MODE_MPI
const int noOfCuboids = singleton::mpi().getSize();
#else
const int noOfCuboids = 1;
#endif

// setup section1
Vector<T,2> extendChannel( nx, ny );
Vector<T,2> originChannel( 0, 0 );
std::shared_ptr<IndicatorF2D<T>> section1 = std::make_shared<IndicatorCuboid2D<T>>( extendChannel, originChannel );

CuboidGeometry2D<T> cGeometry( *(section1), converter.getConversionFactorLength(), noOfCuboids );

// 设置模型周期性边界
//cGeometry.setPeriodicity( false, true );

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

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

prepareGeometry( superGeometry, section1 );

// === 第三步：准备格子 ===
SuperLattice<T, DESCRIPTOR> sLattice1( superGeometry );
SuperLattice<T, DESCRIPTOR> sLattice2( superGeometry );

// 准备格子和边界条件
prepareLattice( sLattice1, sLattice2, converter, superGeometry);

prepareCoupling( sLattice1, sLattice2, superGeometry);

// === 第四步：主循环 ===
int iT = 0;
clout << “starting simulation…” << std::endl;
util::Timer<T> timer( maxIter, superGeometry.getStatistics().getNvoxel() );
timer.start();

for ( iT=0; iT<=maxIter; ++iT ) {

// === 第五步：定义初始化条件及边界条件 ===
//setBoundaryValues( sLatticeOne, sLatticeTwo, force, iT, superGeometry );

// === 第六步：执行碰撞和迁移步骤 ===
sLattice1.collideAndStream();
sLattice2.collideAndStream();

// 执行格子间耦合
sLattice1.executeCoupling();
sLattice2.executeCoupling();

// === 第七步：计算并输出结果 ===
getResults( sLattice1, sLattice2, iT, superGeometry, timer, converter );
}

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

}

[mathias]

You are asking for detailed support which is beyond the scope of the forum. You may consider our next spring school or https://www.openlb.net/consortium/ .

[Author]

Posts