Reply To: Identification of Geometric Model Errors in Multiphase Flow Simulations

[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();

}