[jan]

Dear Jijo,

I’m not sure what I’m seeing in the pictures. Are the particles passing through the wall or are they bouncing back off it?

If they are passing through it:
– Most likely the contact force is too low, perhaps because the Young’s modulus is too low.
– It could also be that the contact handling in your case is not configured correctly. To debug this, I suggest you check if contacts are detected in your case. If they are detected, then the particleContacts or wallContacts of the contactContainer (see https://www.openlb.net/DoxyGen/html/db/d56/structolb_1_1particles_1_1contact_1_1ContactContainer.html) should contain at least one object. You can add this check after coupleResolvedParticlesToLattice.

If they bounce and you’re wondering about the artifact:
If the “automatic” porosity reset doesn’t work for you, you can try to reset the fields manually with resetSuperParticleField (see https://www.openlb.net/DoxyGen/html/da/d81/namespaceolb_1_1particles.html#a9f93e2392819fc95fc4dc52e282893e0) before writing the VTKs.

Best regards,
Jan

[jan]

Dear Jijo,

I agree, that part is advanced. Let me break it down:
1. solidBoundaries is an std::vector of solid boundaries (walls). In this case we only have one, but we can have as many as we like. This can be useful if the walls are made of different materials, or if we want to make sure that contacts in a corner or on an edge are detected separately.
2. In the call of the push_back method, we create an instance of a SolidBoundary and add it to the above mentioned vector.
3. The SolidBoundary is created from a unique pointer of an IndicatorF, a latticeMaterial, a contactMaterial and an optional enlargement solely for the contact treatment.

With this in mind, it may be a little easier to understand. The last two entries are the material numbers. The 2 is the one used on the lattice, e.g. in prepareLattice to set the fluid boundary conditions. The second, the wallContactMaterial, is used to identify the wall properties during the contact treatment. As indicated above. Walls can be of different materials but you are likely to use the same boundary condition for the fluid interaction. Therefore the latticeMaterial could be the same but the contactMaterial is different.

The first parameter, which is probably the most complex looking one, is a simple unique pointer to an indicator that describes the geometry of the wall. Here we use the inverse of a cuboid, which describes the fluid domain. This means that the wall contains everything except the fluid. This is done by the first argument of IndicInverse. The other two arguments take the new min and max values. If we didn’t set them manually and just used the min and max of the fluid, the wall would have a thickness of 0, which would probably cause problems with the distance calculations during the contact treatment.

Kind regards,
Jan

[jan]

Dear Henderson,

Unfortunately, I don’t know your setup. However, I believe that the Δt and contact forces may be too large, which could cause problems like this.

Best regards,
Jan

[jan]

Dear Jijo,

as stated in [1], it is necessary to enlarge the particle for the consideration of the contact. You can do this using particle.setField<NUMERICPROPERTIES, ENLARGEMENT_FOR_CONTACT>( /* choose a value */ );. Otherwise it is to be expected that the incompressible fluid slows down the particle so that there is no contact with the wall, especially for high resolutions.

Kind regards,
Jan

[1]: 10.1016/j.partic.2022.12.005

[jan]

Dear Henderson,

yes, it should work like that, but if you’re using the ParticleManager, you also have to update the external acceleration stored there. Apparently this is not easy in OpenLB 1.5 because it’s a private variable. In the future release you’ll update it like this: particleManager.setExternalAcceleration(accExt);. This function is missing in 1.5. As a workaround, if waiting for the next release is not an option, you could make the variable _externalAcceleration public in src/olb/particles/particlemanager.h and then modify it directly or write a simple setter function yourself.

Best regards,
Jan

[jan]

Dear Henderson,

You can change the mass in the particle field, that is correct. However, that mass is most likely only used to solve the velocity verlet algorithm in your setup. So in that context the density changes, yes. But it doesn’t affect the external acceleration caused by gravity. In the DKT2D example, this acceleration is quantified by the variable accExt, as you can see, the particle density affects it. If you want to change it, you would have to update it with the new particle density.

Best regards,
Jan

[jan]

Dear Henderson,

Unfortunately, I can’t provide you with a solution to the above problem because I don’t know your simulation setup. Here are a few things you might check:
– Is the periodic boundary applied correctly?
– Maybe also check the “geometry” in ParaView, in the direction of the periodic boundary the material should most likely be 1
– Compare your implementation with examples that use a periodic setup, e.g. Poiseuille2d, are there any differences?

For extensive help with your setup I suggest the Spring School: https://www.openlb.net/spring-school-2023/

Best regards,
Jan

[jan]

Dear Henderson,

The periodic boundaries for particles in OpenLB 1.5 doesn’t work as expected, so it was explicitly removed from the ParticleManager to make it unusable. Therefore there is no simple fix for your problem in OpenLB 1.5. This issue has been fixed and will be controllable via the ParticleManager in the next release, as mentioned by Fedor.

Best regards,
Jan

[jan]

Dear Henderson,

Yes, you can change all the fields provided by the PARTICLETYPE used. The following simple example iterates over all particles and gets and sets the position, velocity and force field. The latter corresponds to the acceleration of the particle according to Newton’s second law of motion.

constexpr unsigned D = 3; // or 2 depending on your setup
using namespace descriptors;

// Iterate over all particles
for(int iP = 0; iP < particleSystem.size(); ++iP) {
    // Get the particle
    auto particle = particleSystem.get(iP);

    // Get and set position
    Vector<T,D> position = particle.template getField<GENERAL,POSITION>();
    particle.template setField<GENERAL,POSITION>( position );

    // Get and set velocity
    Vector<T,D> velocity = particle.template getField<MOBILITY,VELOCITY>();
    particle.template setField<MOBILITY,VELOCITY>( velocity );

    // Get and set force
    Vector<T,D> force = particle.template getField<FORCING,FORCE>();
    particle.template setField<FORCING,FORCE>( force );
}

Best regards,
Jan

[jan]

Dear Henderson,

As of OpenLB 1.5, no explicit particle-wall and particle-particle interactions are supported. All coupling is done through the fluid, but as you have noticed, this is not sufficient in many cases. That’s why we have developed a contact model that works for resolved convex particles (10.1016/j.partic.2022.12.005), which will be included in the next release.

Best regards,
Jan

[jan]

Hello Henderson,

I assume that you mean the two surface resolved particle examples (sedimentation3d and dkt2d).

In OpenLB 1.5 the dynamics have been completely reworked. In the process, most forcing schemes were removed, since studies have shown that the exact difference method of Kupershtokh leads to the most accurate results [1].

If you want to manually change the used forcing scheme, you can modify the apply function of the PorousParticle struct in src/dynamics/porousBGKdynamics.h.

Kind regards,
Jan

[1]: https://doi.org/10.3390/computation9020011

[jan]

Dear Anas,

in the main function of the example a cuboidGeometry is created:

/// Instantiation of an empty cuboidGeometry
CuboidGeometry2D<T> cuboidGeometry(cuboid, converter.getPhysDeltaX(), singleton::mpi().getSize());

Here, you can see that the number of processes is provided as an argument (singleton::mpi().getSize()). Apparently, if this equals 1 or 4, the geometry is falsified. This means for you, that the case would probably work as expected if you run it with a different number of processes, e.g., mpirun -np 3 ./squareCavity2d.

In the meantime, we’re working on a fix.

Hope that answers your question, if not, let me know.

Best regards,
Jan

[jan]

Hello Anas,

thanks for the explanation. I’ve checked our current master, and it seems to work there. However, I’m not aware which commit/change fixes it. Therefore, I’ll ask other contributors and we’ll get back to you within a few days.

Best regards,
Jan

[jan]

Hello Anas,

can you give me some more information about your problem? Is it the default squareCavity2D? Did you change something and what did you change? How do you quantify that the right wall is not isothermal? (And maybe other things that you thing are of interest.)

Best regards,
Jan

[jan]

Hello Anas,

as you’re regarding a thermal flow case, those properties are handled by the ThermalUnitConverter (https://www.openlb.net/DoxyGen/html/d0/de2/classolb_1_1ThermalUnitConverter.html). For example, comparing the constructor with the usage in the squareCavity2D case shows that lx is the characteristic length and is set at the beginning of the main function. Similar is true for other values.
The relaxation time, on other hand, is not a parameter of the used constructor, but is calculated from other values (see the implementation of the constructor: https://www.openlb.net/DoxyGen/html/d0/de2/classolb_1_1ThermalUnitConverter.html#a60ea34c354f56edbee5b720bf2cadac4).

If this does not answer your question, please be more specific.

Best regards,
Jan

[Author]

Posts