Problems of Resolved Particles with Periodic Boundary

[yueyq]

Dear Community,

I tried to modify the example code dkt2d.cpp to be periodic along x-axis (left and right), but failed. The fluid part is good, but the solid part has two problems:

1. I set bottom boundary to be wall, but it works out that the left and right boundary exert forces to particles as wall when contacting. It shows that the four boundaries share same “wallID”, even though they have different material number. Plus, I have set defaults::periodicity {true, false}.

2. There is no longer momentum exchanged between particle and fluid after the particles left the domain and entered from the other side (i.e. periodic boundary).

Please give me some advice. Thanks a lot.

[jan]

Dear yueyq,

Without having the specific details of your code, I can only make assumptions about the issues you’re facing. Here are a couple of points to consider:

1. It’s possible that you’ve defined the solidBoundaries incorrectly. Ensure that the solidBoundaries are not set for the entire surrounding part, but only for the single wall at the bottom of the container. Updating this might resolve the issue with forces being exerted by the left and right boundaries.

2. Double-check how you’ve set the periodicity. Also, incorrect definition of the walls might cause the force from the wall contact to override the forces from the fluid. This could be the reason why there’s no momentum exchange between the particle and the fluid after the particles leave and re-enter the domain.

I hope these suggestions help you resolve the issues in your code. If you need further assistance or have any other questions, please don’t hesitate to ask.

Best regards,
Jan

[yueyq]

Dear jan,

Thanks for your advice. The code is mostly the same as examples/particles/dkt2d. I will put my code then at another reply post for you to facilitate your reading:)

1. I did get puzzled when I define the “SolidBoundaries”. As the “SolidBoundary” code below, I wonder if “std::unique_ptr<IndicatorF<T, D>> _indPtr” was the whole simulation domain, or just the user-defined boundary wall part? I tried to set the bottom of the container as the “_indPtr”, but particles just went out of the container at the beginning. Plus, why is the cuboid extended larger in the dkt2d example?
struct SolidBoundary {
private:
/// Indicator providing the surface and volume information
std::unique_ptr<IndicatorF<T, D>> _indPtr;
/// Material number on the lattice for identification and handling of boundaries
const unsigned _latticeMaterial;
}

2. I checked again about the periodicity, but the problem still exists. A speculation is that forces of the “ghostPos” are overwritten by “position” according to the line 190-195 of src/functors/lattice/latticeMomentumExchangeForce.hh:
// Calculate force for the ghost particle
particle.template setField<GENERAL,POSITION>(ghostPos);
evaluate(output, particle, iP);
// Calculate force of actual particle
particle.template setField<GENERAL,POSITION>(position);
evaluate(output, particle, iP);
Of course, incorrect definition of the walls could also be the reason since it has not been solved yet.

Really appreciate your suggestions and help. Look forward to your reply.

Sincerely,
yueyq

[yueyq]

#define NEW_FRAMEWORK

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

using namespace olb;
using namespace olb::descriptors;
using namespace olb::graphics;
using namespace olb::particles;
using namespace olb::particles::access;
using namespace olb::particles::contact;
using namespace olb::particles::dynamics;
using namespace olb::util;

using T = FLOATING_POINT_TYPE;

//Define lattice type
typedef PorousParticleWithContactD2Q9Descriptor DESCRIPTOR;

//Define particle type
typedef ResolvedCircleWithContact2D PARTICLETYPE;

// Define particle-particle contact type
typedef ParticleContactArbitraryFromOverlapVolume<T, DESCRIPTOR::d, true>
PARTICLECONTACTTYPE;
// Define particle-wall contact type
typedef WallContactArbitraryFromOverlapVolume<T, DESCRIPTOR::d, true>
WALLCONTACTTYPE;

#define WriteVTK
#define WriteGnuPlot

std::string gnuplotFilename = “gnuplot.dat”;

// Parameters for the simulation setup
int N = 1;
int M = N;

T eps = 0.5; // eps*latticeL: width of transition area

T maxPhysT = 6.; // max. simulation time in s, SI unit
T iTwrite = 0.125; //converter.getLatticeTime(.3);

T lengthX = 0.02;
T lengthY = 0.08;

T centerX1 = 0.01;
T centerY1 = 0.068;
Vector<T,2> center1 = {centerX1,centerY1};
T centerX2 = 0.00999;
T centerY2 = 0.072;
Vector<T,2> center2 = {centerX2,centerY2};

T rhoP = 1010.;
T radiusP = 0.001;
Vector<T,2> accExt = {.0, -T(9.81) * (T(1) – T(1000) / rhoP)};

unsigned contactBoxResolutionPerDirection = 8;
unsigned particleContactMaterial = 0;
unsigned wallContactMaterial = 0;
T youngsModulus = 1e6;
T poissonRatio = 0.3;
T coefficientOfRestitution = 0.9;
T coefficientStaticFriction = 0.6;
T coefficientKineticFriction = 0.3;

namespace custom{
template <unsigned D>
const auto customperiodicity = [](){
if constexpr (D==3) {
return Vector<bool,3>(true,true,false);
} else {
return Vector<bool,2>(true,false);
}
};
}

void prepareGeometry(UnitConverter<T,DESCRIPTOR> const& converter,
SuperGeometry<T,2>& superGeometry)
{
OstreamManager clout(std::cout, “prepareGeometry”);
clout << “Prepare Geometry …” << std::endl;

superGeometry.rename(0, 2);
//superGeometry.rename(2, 1, {1, 1});
superGeometry.rename(2, 1, {0, 1});

superGeometry.clean();
superGeometry.innerClean();

superGeometry.checkForErrors();
superGeometry.getStatistics().print();
clout << “Prepare Geometry … OK” << std::endl;
}

void prepareLattice(
SuperLattice<T, DESCRIPTOR>& sLattice, UnitConverter<T,DESCRIPTOR> const& converter,
SuperGeometry<T,2>& superGeometry)
{
OstreamManager clout(std::cout, “prepareLattice”);
clout << “Prepare Lattice …” << std::endl;

/// Material=0 –>do nothing
sLattice.defineDynamics<PorousParticleBGKdynamics>(superGeometry, 1);
setBounceBackBoundary(sLattice, superGeometry, 2);

sLattice.setParameter<descriptors::OMEGA>(converter.getLatticeRelaxationFrequency());

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

void setBoundaryValues(SuperLattice<T, DESCRIPTOR>& sLattice,
UnitConverter<T,DESCRIPTOR> const& converter,
SuperGeometry<T,2>& superGeometry)
{
OstreamManager clout(std::cout, “setBoundaryValues”);

AnalyticalConst2D<T, T> one(1.);
sLattice.defineField<POROSITY>(superGeometry.getMaterialIndicator({1,2}), one);

// Set initial condition
AnalyticalConst2D<T, T> ux(0.);
AnalyticalConst2D<T, T> uy(0.);
AnalyticalConst2D<T, T> rho(1.);
AnalyticalComposed2D<T, T> u(ux, uy);

//Initialize all values of distribution functions to their local equilibrium
sLattice.defineRhoU(superGeometry, 1, rho, u);
sLattice.iniEquilibrium(superGeometry, 1, rho, u);

// Make the lattice ready for simulation
sLattice.initialize();
}

void getResults(SuperLattice<T, DESCRIPTOR>& sLattice,
UnitConverter<T,DESCRIPTOR> const& converter, int iT,
SuperGeometry<T,2>& superGeometry, Timer<double>& timer,
ParticleSystem<T,PARTICLETYPE>& particleSystem )
{
OstreamManager clout(std::cout, “getResults”);

#ifdef WriteVTK
SuperVTMwriter2D<T> vtkWriter(“sedimentation”);
SuperLatticePhysVelocity2D<T, DESCRIPTOR> velocity(sLattice, converter);
SuperLatticePhysPressure2D<T, DESCRIPTOR> pressure(sLattice, converter);
SuperLatticePhysExternalPorosity2D<T, DESCRIPTOR> externalPor(sLattice, converter);
SuperLatticeMomentumExchangeForceLocal<T, DESCRIPTOR, PARTICLETYPE> momentumExchange(
sLattice, converter, superGeometry, particleSystem);

vtkWriter.addFunctor(velocity);
vtkWriter.addFunctor(pressure);
vtkWriter.addFunctor(externalPor);
vtkWriter.addFunctor(momentumExchange);

if (iT == 0) {
converter.write(“dkt”);
SuperLatticeGeometry2D<T, DESCRIPTOR> geometry(sLattice, superGeometry);
SuperLatticeCuboid2D<T, DESCRIPTOR> cuboid(sLattice);
SuperLatticeRank2D<T, DESCRIPTOR> rank(sLattice);
vtkWriter.write(geometry);
vtkWriter.write(cuboid);
vtkWriter.write(rank);
vtkWriter.createMasterFile();
}

if (iT % converter.getLatticeTime(iTwrite) == 0) {
vtkWriter.write(iT);
}
#endif

auto particleA = particleSystem.get( 0 );
auto particleB = particleSystem.get( 1 );

#ifdef WriteGnuPlot
if (iT % converter.getLatticeTime(iTwrite) == 0) {
if (singleton::mpi().getRank() == 0) {

std::ofstream myfile;
myfile.open (gnuplotFilename.c_str(), std::ios::app);
T p2PosY = particleB.getField<GENERAL,POSITION>()[1];
T p1PosY = particleA.getField<GENERAL,POSITION>()[1];
T p2PosX = particleB.getField<GENERAL,POSITION>()[0];
T p1PosX = particleA.getField<GENERAL,POSITION>()[0];
myfile
<< converter.getPhysTime(iT) << ” ”
<< std::setprecision(9)
<< p2PosY << ” ”
<< p1PosY << ” ”
<< p2PosX << ” ”
<< p1PosX << std::endl;
myfile.close();
}
}
#endif

/// Writes output on the console
if (iT % converter.getLatticeTime(iTwrite) == 0) {
timer.update(iT);
timer.printStep();
sLattice.getStatistics().print(iT, converter.getPhysTime(iT));
for (std::size_t iP=0; iP<particleSystem.size(); ++iP) {
auto particle = particleSystem.get(iP);
io::printResolvedParticleInfo(particle);
}
}

return;
}

int main(int argc, char* argv[])
{
/// === 1st Step: Initialization ===
olbInit(&argc, &argv);
singleton::directories().setOutputDir(“./tmp/”);
OstreamManager clout(std::cout, “main”);

UnitConverter<T,DESCRIPTOR> converter(
( T ) 0.0001/ N, //physDeltaX
( T ) 5.e-4/(N*M), //physDeltaT,
( T ) .002, //charPhysLength
( T ) 0.2, //charPhysVelocity
( T ) 1E-6, //physViscosity
( T ) 1000. //physDensity
);
converter.print();

/// === 2nd Step: Prepare Geometry ===
std::vector<T> extend(2, T());
extend[0] = lengthX;
extend[1] = lengthY;
std::vector<T> origin(2, T());
IndicatorCuboid2D<T> cuboid(extend, origin);

#ifdef PARALLEL_MODE_MPI
CuboidGeometry2D<T> cuboidGeometry(cuboid, converter.getConversionFactorLength(), singleton::mpi().getSize());
#else
CuboidGeometry2D<T> cuboidGeometry(cuboid, converter.getConversionFactorLength(), 1);
#endif

cuboidGeometry.setPeriodicity(true, false);

HeuristicLoadBalancer<T> loadBalancer(cuboidGeometry);
SuperGeometry<T,2> superGeometry(cuboidGeometry, loadBalancer, 2);
prepareGeometry(converter, superGeometry);

/// === 3rd Step: Prepare Lattice ===
SuperLattice<T, DESCRIPTOR> sLattice(superGeometry);

prepareLattice(sLattice, converter, superGeometry);

/// === 4th Step: Main Loop with Timer ===
Timer<double> timer(converter.getLatticeTime(maxPhysT), superGeometry.getStatistics().getNvoxel());
timer.start();

// Create ParticleSystem
ParticleSystem<T,PARTICLETYPE> particleSystem;

Vector<bool,DESCRIPTOR::d> period (true, false);

//Create particle manager handling coupling, gravity and particle dynamics
ParticleManager<T,DESCRIPTOR,PARTICLETYPE> particleManager(
particleSystem, superGeometry, sLattice, converter, accExt, period);

// Create and assign resolved particle dynamics
particleSystem.defineDynamics<
VerletParticleDynamics<T,PARTICLETYPE>>();

// Create solid boundaries for particle interaction
std::vector<SolidBoundary<T, DESCRIPTOR::d>> solidBoundaries;
solidBoundaries.push_back( SolidBoundary<T, DESCRIPTOR::d>(
std::make_unique<IndicInverse<T, DESCRIPTOR::d>>(
cuboid, cuboid.getMin() – 5 * converter.getPhysDeltaX(),
cuboid.getMax() + 5 * converter.getPhysDeltaX()), 2, wallContactMaterial));

// Create objects for contact treatment
ContactContainer<T, PARTICLECONTACTTYPE, WALLCONTACTTYPE> contactContainer;
// Generate lookup table for contact properties
ContactProperties<T, 1> contactProperties;
contactProperties.set(particleContactMaterial, wallContactMaterial,
evalEffectiveYoungModulus(youngsModulus, youngsModulus,
poissonRatio, poissonRatio),
coefficientOfRestitution, coefficientKineticFriction, coefficientStaticFriction);

T epsilon = eps * converter.getConversionFactorLength();
T radius = radiusP;

// Create Particle 1
creators::addResolvedCircle2D( particleSystem, center1,
radius, epsilon, rhoP );

// Create Particle 2
creators::addResolvedCircle2D( particleSystem, center2,
radius, epsilon, rhoP );

//Check ParticleSystem
particleSystem.checkForErrors();

// Set contact material
for (std::size_t iP = 0; iP < particleSystem.size(); ++iP) {
auto particle = particleSystem.get(iP);
setContactMaterial(particle, particleContactMaterial);
}

/// === 5th Step: Definition of Initial and Boundary Conditions ===
setBoundaryValues(sLattice, converter, superGeometry);

{
auto& communicator = sLattice.getCommunicator(stage::PostPostProcess());
communicator.requestOverlap(sLattice.getOverlap());
communicator.requestFields<POROSITY,VELOCITY_NUMERATOR,VELOCITY_DENOMINATOR>();
communicator.exchangeRequests();
}

clout << “MaxIT: ” << converter.getLatticeTime(maxPhysT) << std::endl;
for (std::size_t iT = 0; iT < converter.getLatticeTime(maxPhysT)+10; ++iT) {

// Execute particle manager
particleManager.execute<
couple_lattice_to_particles<T,DESCRIPTOR,PARTICLETYPE>,
apply_gravity<T,PARTICLETYPE>
>();

const T k = static_cast<T>(4. / (3 * util::sqrt(M_PI)));

decltype(custom::customperiodicity<2>) customperiod = custom::customperiodicity<2>;

// Calculate and apply contact forces
processContacts<T, PARTICLETYPE, PARTICLECONTACTTYPE, WALLCONTACTTYPE, ContactProperties<T, 1>>(
particleSystem, solidBoundaries, contactContainer, contactProperties,
superGeometry, contactBoxResolutionPerDirection, k, customperiod);

// Solve equations of motion
particleManager.execute<process_dynamics<T,PARTICLETYPE>>();

// Couple particles to lattice (with contact detection)
coupleResolvedParticlesToLattice<T, DESCRIPTOR, PARTICLETYPE, PARTICLECONTACTTYPE, WALLCONTACTTYPE>(
particleSystem, contactContainer, superGeometry, sLattice, converter, solidBoundaries, customperiod);

// Get Results
getResults(sLattice, converter, iT, superGeometry, timer, particleSystem);

// Collide and stream
sLattice.collideAndStream();
}

// Run Gnuplot
if (singleton::mpi().getRank() == 0) {
if (!system(NULL)) {
exit (EXIT_FAILURE);
}
int ret = system(“gnuplot dkt.p”);
if (ret == -1) {
clout << “Writing Gnuplot failed!” << std::endl;
}
}

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

[jan]

Dear yueyq,

solidBoundaries.push_back( SolidBoundary<T, DESCRIPTOR::d>(
std::make_unique<IndicInverse<T, DESCRIPTOR::d>>(
cuboid, cuboid.getMin() – 5 * converter.getPhysDeltaX(),
cuboid.getMax() + 5 * converter.getPhysDeltaX()), 2, wallContactMaterial));

the above part takes the inverse of the whole computation domain, which means each side is regarded as a wall. If you solely want the bottom to act as a wall, then I suggest that you introduce a IndicatorCuboid instead and place it at the bottom wall. Be sure to place it exactly where the bottom wall is and give it a sufficient thickness, because the particles will overlap it.

That’s also part of the reason why the SolidBoundary in the example is extended. If you’d have a wall with a infinitesimal thickness, it would lead to problems during the contact treatment, as we want to calculate the overlap volume.

Best regards,
Jan

• This reply was modified 1 month, 3 weeks ago by jan.

[Author]

Posts