Could not addVelocityBoundary

[wf_guoyl]

Hello developers,

I recently had a problem with the addition of velocity boundary conditions. I simulate the evaporation of droplets on the wall. Rectangles represent solids. All four walls have temperature values, but the temperature on the left wall is higher than the temperature on the right wall. Now I want to add an inlet velocity u on the left wall, but I get an error message “Could not addVelocityBoundary (2,3) …… (2,69)”. I don’t know what causes this problem, who can help me? Here are some codes:

const int materialBase = 0;
const int materialFluid = 1;
const int materialLowerWall = 2;
const int materialUpperWall = 3;
const int materialInflow = 4;
const int materialOutflow = 5;
const int materialObstacle = 6;

void prepareGeometry( SuperGeometry2D<SimType>& superGeometry )
{

  OstreamManager clout( std::cout,"prepareGeometry" );
  clout << "Prepare Geometry ..." << std::endl;

  superGeometry.rename( materialBase, materialLowerWall );

  Vector<SimType,2> extend( lx0, ly0 );
  Vector<SimType,2> origin;

  // create fluid and lower wall
  superGeometry.rename( materialLowerWall, materialFluid, 0, 1 );

  // create upper wall
  origin = { 0, ly0-conversionLength };
  extend = { lx0, conversionLength };
  olb::IndicatorCuboid2D<SimType> upperBC( extend, origin );
  superGeometry.rename( materialLowerWall, materialUpperWall, materialFluid, upperBC );
  
  // create inflow
  origin = { 0, 0 };
  extend = { conversionLength, ly0-conversionLength };
  olb::IndicatorCuboid2D<SimType> inflow( extend, origin );
  superGeometry.rename( materialFluid, materialInflow, inflow ); 

  // create outflow
  origin = { lx0-conversionLength, 0 };
  extend = { conversionLength, ly0-conversionLength };
  olb::IndicatorCuboid2D<SimType> outflow( extend, origin);
  superGeometry.rename( materialFluid, materialOutflow, outflow);

  // creat obstacle
  origin = { lx0/2-50*conversionLength, ly0/2 };
  extend = { 100*conversionLength, 25*conversionLength };
  olb:: IndicatorCuboid2D<SimType> obstacle( extend, origin);
  superGeometry.rename( materialFluid, materialObstacle, obstacle);

  // Removes all not needed boundary voxels outside the surface
  //superGeometry.clean();
  // Removes all not needed boundary voxels inside the surface
  //superGeometry.innerClean();
  superGeometry.checkForErrors();

  superGeometry.print();

  clout << "Prepare Geometry ... OK" << std::endl;
}

// Set up the geometry of the simulation
void prepareLattice( olb::ThermalUnitConverter<SimType, NSDESCRIPTOR, TDESCRIPTOR> const& converter,
                     olb::SuperLattice2D<SimType, NSDESCRIPTOR>& NSlattice,
                     olb::SuperLattice2D<SimType, TDESCRIPTOR>& ADlattice,
                     olb::Dynamics<SimType, NSDESCRIPTOR>& bulkDynamics,
                     olb::Dynamics<SimType, TDESCRIPTOR>& advectionDiffusionBulkDynamics,
                     olb::Dynamics<SimType, TDESCRIPTOR>& advectionDiffusionObstacleDynamics,
                     olb::BounceBack<SimType, NSDESCRIPTOR>& bounceBackUpperWall,
                     olb::BounceBack<SimType, NSDESCRIPTOR>& bounceBackLowerWall,
                     olb::sOnLatticeBoundaryCondition2D<SimType,TDESCRIPTOR>& TboundaryCondition,
                     olb::sOnLatticeBoundaryCondition2D<SimType,NSDESCRIPTOR>& NSboundaryCondition,                     
                     SuperGeometry2D<SimType>& superGeometry,
                     SimType TempInitial,
                     SimType rhoVapor, SimType rhoLiquid )
{

  OstreamManager clout( std::cout,"prepareLattice" );
  clout << "Prepare Lattice ..." << endl;

  const SimType omega = converter.getLatticeRelaxationFrequency();
  const SimType Tomega  = converter.getLatticeThermalRelaxationFrequency();

  clout << "NS omega = " << omega << endl;
  clout << "T  omega = " << Tomega << endl;

  auto bulkIndicator = superGeometry.getMaterialIndicator({materialFluid, materialLowerWall, materialUpperWall, materialInflow, materialOutflow, materialObstacle});
  auto fluidIndicator = superGeometry.getMaterialIndicator({materialFluid, materialInflow, materialOutflow});
  auto wallIndicator = superGeometry.getMaterialIndicator({materialUpperWall, materialLowerWall});
  auto obstacleIndicator = superGeometry.getMaterialIndicator({materialObstacle});

  clout << "Define NS Dynamics ..." << endl;
  NSlattice.defineDynamics( superGeometry, materialBase, &olb::instances::getNoDynamics<SimType, NSDESCRIPTOR>() );
  NSlattice.defineDynamics( fluidIndicator, &bulkDynamics );
  NSlattice.defineDynamics( superGeometry, materialUpperWall, &bounceBackUpperWall );
  NSlattice.defineDynamics( superGeometry, materialLowerWall, &bounceBackLowerWall );
  NSlattice.defineDynamics( superGeometry, materialObstacle, &olb::instances::getBounceBack<SimType, NSDESCRIPTOR>() );

  clout << "Define AD Dynamics ..." << endl;
  ADlattice.defineDynamics( superGeometry, materialBase, &olb::instances::getNoDynamics<SimType, TDESCRIPTOR>() );
  ADlattice.defineDynamics( fluidIndicator, &advectionDiffusionBulkDynamics );
  ADlattice.defineDynamics( wallIndicator, &advectionDiffusionBulkDynamics );
  ADlattice.defineDynamics( obstacleIndicator, &advectionDiffusionObstacleDynamics );

  NSboundaryCondition.addVelocityBoundary(superGeometry, materialInflow, omega);
  //NSboundaryCondition.addPressureBoundary(superGeometry, materialOutflow, omega);

  TboundaryCondition.addTemperatureBoundary(superGeometry, materialUpperWall, Tomega);
  TboundaryCondition.addTemperatureBoundary(superGeometry, materialLowerWall, Tomega);
  TboundaryCondition.addTemperatureBoundary(superGeometry, materialInflow, Tomega);
  TboundaryCondition.addTemperatureBoundary(superGeometry, materialOutflow, Tomega);

  clout << "Define Initial Conditions ..." << endl;
  // Fluid Initial conditions
  olb::AnalyticalConst2D<SimType,SimType> zeroVelocity( 0,0 );

  AnalyticalConst2D<SimType,SimType> u( 0.01, 0 );

  AnalyticalConst2D<SimType,SimType> noiseVap( rhoVapor*0.1 );
  AnalyticalConst2D<SimType,SimType> noiseLiq( rhoLiquid*0.1 );
  AnalyticalRandom2D<SimType,SimType> random;

  SimType estimatedAdditionalAverageVapour = 2*(wallRho - rhoVapor)*conversionLength/ly0;
  clout << "estimatedAdditionalAverageVapour: " << estimatedAdditionalAverageVapour << std::endl;
  // this estimate assumes that 2 full fluid nodes equal to the wall density are needed to compensate for the wall density gradient
  // This additional density is then spread over the entire y domain

  olb::AnalyticalConst2D<SimType,SimType> constRhoVapour( rhoVapor+estimatedAdditionalAverageVapour );
  olb::AnalyticalConst2D<SimType,SimType> constRhoLiquid( rhoLiquid * rhoLiquidMultiplier );
  olb::AnalyticalConst2D<SimType,SimType> constRhoUpperBC( rhoVapor );
  olb::AnalyticalConst2D<SimType,SimType> constRhoWall( rhoVapor );

  // Temperature Inital Conditions
  olb::AnalyticalConst2D<SimType,SimType> T_init( converter.getLatticeTemperature(TempInitial));

  // Random density distribution across inlet
  olb::AnalyticalIdentity2D<SimType,SimType> noiseIndicatorVap( random*noiseVap  );
  olb::AnalyticalIdentity2D<SimType,SimType> noiseIndicatorLiq( random*noiseLiq  );
  olb::AnalyticalIdentity2D<SimType,SimType> rhoBase( constRhoVapour );

  olb::SmoothIndicatorCircle2D<SimType,SimType> sphere( sphereOrigin, sphereRadius, sphereInterfaceWidth);

  olb::AnalyticalIdentity2D<SimType,SimType> sphereIndicator( sphere );

  olb::AnalyticalIdentity2D<SimType,SimType> rho3( rhoBase + noiseIndicatorVap + sphereIndicator*( constRhoLiquid - constRhoVapour + noiseIndicatorLiq) );

  //Initialize all values of distribution functions to their local equilibrium

  // Setting initial psi values
  olb::AnalyticalConst2D<SimType,SimType> initialPsi( 1.0 );
  NSlattice.defineField<olb::descriptors::PSI_PSEUDO_RHO>( fluidIndicator, initialPsi ) ; 

  // Fluid init
  NSlattice.defineRhoU( superGeometry, materialFluid, rho3, zeroVelocity );
  NSlattice.iniEquilibrium( superGeometry, materialFluid, rho3, zeroVelocity );
  
  NSlattice.defineRhoU( superGeometry, materialUpperWall, constRhoUpperBC, zeroVelocity );
  NSlattice.iniEquilibrium( superGeometry, materialUpperWall, constRhoUpperBC, zeroVelocity );
  
  NSlattice.defineRhoU( superGeometry, materialLowerWall, constRhoWall, zeroVelocity );
  NSlattice.iniEquilibrium( superGeometry, materialLowerWall, constRhoWall, zeroVelocity ); 

  NSlattice.defineRhoU( superGeometry, materialInflow, rho3, u );
  NSlattice.iniEquilibrium( superGeometry, materialInflow, rho3, u ); 

  NSlattice.defineRhoU( superGeometry, materialOutflow, rho3, zeroVelocity );
  NSlattice.iniEquilibrium( superGeometry, materialOutflow, rho3, zeroVelocity ); 

  // Thermal init
  ADlattice.defineRho( superGeometry, materialFluid, T_init );
  ADlattice.iniEquilibrium( superGeometry, materialFluid, T_init, zeroVelocity );

  ADlattice.defineRho( superGeometry, materialUpperWall, T_init );
  ADlattice.iniEquilibrium( superGeometry, materialUpperWall, T_init, zeroVelocity );

  ADlattice.defineRho( superGeometry, materialLowerWall, T_init );
  ADlattice.iniEquilibrium( superGeometry, materialLowerWall, T_init, zeroVelocity );
  
  ADlattice.defineRho( superGeometry, materialInflow, T_init );
  ADlattice.iniEquilibrium( superGeometry, materialInflow, T_init, zeroVelocity );

  ADlattice.defineRho( superGeometry, materialOutflow, T_init );
  ADlattice.iniEquilibrium( superGeometry, materialOutflow, T_init, zeroVelocity );
  
  ADlattice.defineRho( obstacleIndicator, T_init );
  ADlattice.iniEquilibrium( obstacleIndicator, T_init, zeroVelocity );
  
  // Make the lattice ready for simulation
  NSlattice.initialize();
  ADlattice.initialize();

  clout << "Prepare Lattice ... OK" << std::endl;
}

The code used to work, but I got an error when I added the speed boundary condition.

Many thanks!

[wf_guoyl]

Referring to case “cylinder2d”, the adjustment is as follows:

previous code:

 // create inflow
  origin = { 0, 0 };
  extend = { conversionLength, ly0-conversionLength };
  olb::IndicatorCuboid2D<SimType> inflow( extend, origin );
  superGeometry.rename( materialFluid, materialInflow, inflow ); 

  // create outflow
  origin = { lx0-conversionLength, 0 };
  extend = { conversionLength, ly0-conversionLength };
  olb::IndicatorCuboid2D<SimType> outflow( extend, origin);
  superGeometry.rename( materialFluid, materialOutflow, outflow);

the adjusted code — version 1.0

// create inflow
  origin = { -conversionLength, 0 };
  extend = { conversionLength, ly0-conversionLength };
  olb::IndicatorCuboid2D<SimType> inflow( extend, origin );
  superGeometry.rename( materialFluid, materialInflow, inflow ); 

  // create outflow
  origin = { lx0-conversionLength, 0 };
  extend = { 2*conversionLength, ly0-conversionLength };
  olb::IndicatorCuboid2D<SimType> outflow( extend, origin);
  superGeometry.rename( materialFluid, materialOutflow, outflow);

Error messages disappears.

the adjusted code — version 2.0

// create inflow
  origin = { -conversionLength, 0 };
  extend = { 2*conversionLength, ly0-conversionLength };
  olb::IndicatorCuboid2D<SimType> inflow( extend, origin );
  superGeometry.rename( materialFluid, materialInflow, inflow ); 

  // create outflow
  origin = { lx0-conversionLength, 0 };
  extend = { 2*conversionLength, ly0-conversionLength };
  olb::IndicatorCuboid2D<SimType> outflow( extend, origin);
  superGeometry.rename( materialFluid, materialOutflow, outflow);

Error message appears.

But I compared the cloud images of the three programs, it looks like these three have the same result. And this result feels like a velocity boundary condition has been applied in the left inlet.

So my problems are:
1.This error “Could not addVelocityBoundary (2,3) …… (2,69)” is generally related to which situation, why can the first version solve this problem.

2.When we define geometry, why do we have a little extra region on the left and right?

I’d appreciate it if anyone could answer. If I don’t give a clear description, I can send the program to your email. Thank you very much!

[Adrian]

Warnings of this kind are printed when the boundary setter (the addVelocityBoundary method in this case) fails to instantiate the necessary dynamics / post processors for a particular cell. This means that the boundary for this cell is not well defined. However in practice this is frequently overlooked as single erroneous populations don’t necessarily impact the overall simulation in a noticeable fashion. This seems to also be the case here?

Of course it is good practice to keep an eye on such messages and to fix the underlying geometry problem. In this case I suspect that the discrete wall normal could not be recovered correctly / the wall normal can not be handled by the specific velocity condition.

Did you already verify the material number setup in ParaView?

w.r.t. your second question: Depending on what you mean, this is either an artifact of the specific geometry setup (rounding or otherwise) or the overlap area around each block that is used for communication.

[Fany]

Dear wf_guoyl,

Your code seemed to simulate single component multiphase flow. I am simulating the flow boiling based on Julius Weinmiller’s work. But it always had a problem of “density is NAN”(“uMax=0.523235; avEnergy=-nan; avRho=nan”). I found uMax sharply increased from 0.007 to 0.53 after the endTimeInit. Have you encounterred this problem?
I would like to ask how you set the MaxVelocity. should it be same with the inlet velocity? I appreciate very much if you can share you code to me.

Best regards,
Fany

[Author]

Posts