Skip to content

Problems of Resolved Particles with Periodic Boundary

OpenLB – Open Source Lattice Boltzmann Code Forums on OpenLB General Topics Problems of Resolved Particles with Periodic Boundary

Viewing 5 posts - 1 through 5 (of 5 total)
  • Author
    Posts
  • #8895
    yueyq
    Participant

    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.

    #8905
    jan
    Participant

    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

    #8910
    yueyq
    Participant

    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

    #8911
    yueyq
    Participant

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

    #8921
    jan
    Participant

    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 4 days, 2 hours ago by jan.
Viewing 5 posts - 1 through 5 (of 5 total)
  • You must be logged in to reply to this topic.