Shan-Chen

[xinhui]

After a period of self-study, I still find it somewhat difficult to use the ShanChen model. Could you provide any relevant simulation cases, preferably using ShanChen93, for single-component two-phase flow or immiscible two-component flow for learning purposes?

[mathias]

There will be an update in the next release. Yet, we need your and the other users and developeres support, see https://www.openlb.net/consortium/.

[xinhui]

Thank you for your reply! I am very much looking forward to the update of the example.

Currently, I want to simulate gas-liquid two-phase flow (without phase change), and I need to set the contact angle at the boundary. In the existing examples, the contact angle settings are all based on the free energy model. However, I am not sure where to find information about the relevant settings for the ShanChen model (such as setting the desity of soilds). Do you have any recommendations for learning resources?

[luizeducze]

Dear xinhui,

Are you using single-component or multi-component ShanChen method?

You can create a layer with a different material number than the fluid to be your solid.
For example, your fluid can be material number 1 and the solid 5.

In your prepare lattice function you can set use the sLattice.defineRho() to set the density in different material numbers. You can find this in many examples in OpenLB.

Also, in your COUPLING, you can choose the material numbers where the operations will be applied, as an example: coupling.restrictTo(superGeometry.getMaterialIndicator({1}));

Since you want to compute the interaction in the solid for the contact angle, you can do that:
coupling.restrictTo(superGeometry.getMaterialIndicator({1}));
coupling.restrictTo(superGeometry.getMaterialIndicator({5}));

Check the binary shear flow example to see how it is applied.

In this way, you can prescribe a density to your solid in prepare lattice, then tell your code to compute the interaction from the solid. Remember to use the bounce-back or other boundary condition for your solid (material number 5) in your prepare lattice.

If you have more doubts just ask.

Kind regards,
Luiz

[Bobbie]

Dear Luiz,

In our current investigation of multiphase flow systems, a critical issue persists regarding solid material numbers. The following code integrates the example of cylinder2d and microFluidcs2d. While the cylindrical region has been specifically assigned material number to zero through phase tagging to represent solid boundaries, subsequent visualization through ParaView reveals liquid density(rho) inside cylindrical region.The results is wrong obviously.So, What fundamental adjustments would ensure the correct running of this code.

The codes are as follows:

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

}

Thank you for your sincere help.

Best regards,
Bobbie

[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