Wall contact in particle simulation (dkt2d example)

[avrachan1]

Dear Community,

I do the following steps starting with the example dkt2d ( sedimentation of 2d particles using porousparticleBGKdynamics).

I switch off gravity.

ParticleManager<T,DESCRIPTOR,PARTICLETYPE> particleManager(
particleSystem, superGeometry, sLattice, converter);

and

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

And then when I run the simulation, I get what I expect – particles don’t move – nothing happens.

Now I change the geometry of the wall. Instead of walls on 4 sides, I would like to just have a wall at the bottom.

Vector<T,2> extend_mL(lengthX,L/2.);
Vector<T,2> origin_mL(0.,-L/2);
IndicatorCuboid2D<T> bottombox(extend_mL,origin_mL);

solidBoundaries.push_back( SolidBoundary<T, DESCRIPTOR::d>(
std::make_unique<IndicInverse<T, DESCRIPTOR::d>>(bottombox), 2, wallContactMaterial));

Ideally nothing should change after this, but I see that the particles are being moved towards the wall.

output after 1 timestep

[ParticleInfo] Particle ID=0
[ParticleInfo] |Circum radius(m)= 0.001025
[ParticleInfo] |Mass(kg)= 0.00317301
[ParticleInfo] |Position(m)= ( 0.01, 0.068)
[ParticleInfo] |Angle(°)= ( 0)
[ParticleInfo] |Velocity(m/s)= ( -5.45708e-05, -145059)
[ParticleInfo] |Ang. Velocity(°/s)= ( -0.0561916)
[ParticleInfo] |Force(N)= ( -0.000692615, -1.8411e+06)
[ParticleInfo] |Acceleration(m/s^2)= ( -0.218283, -5.80237e+08)
[ParticleInfo] |Ang. acc.(°/s^2)= ( -224.767)
[ParticleInfo] |Moment of inertia(kg m^2)= ( 1.5865e-09)
[ParticleInfo] Particle ID=1
[ParticleInfo] |Circum radius(m)= 0.001025
[ParticleInfo] |Mass(kg)= 0.00317301
[ParticleInfo] |Position(m)= ( 0.00999, 0.072)
[ParticleInfo] |Angle(°)= ( 0)
[ParticleInfo] |Velocity(m/s)= ( -5.67728, -18860.8)
[ParticleInfo] |Ang. Velocity(°/s)= ( -5.01562e+06)
[ParticleInfo] |Force(N)= ( -72.0562, -239382)
[ParticleInfo] |Acceleration(m/s^2)= ( -22709.1, -7.54432e+07)
[ParticleInfo] |Ang. acc.(°/s^2)= ( -2.00625e+10)
[ParticleInfo] |Moment of inertia(kg m^2)= ( 1.5865e-09)

What am I doing wrong ?

When inspecting the source code I see that a wall contact is formed in line 421 of src/particles/resolvedblockLatticeInteraction.hh.

where the following if condition becomes true
if (solidBoundary.getIndicator()->signedDistance(
PhysR<T, D>(physR)) < solidBoundaryDetectionDistance) {
The signed distance,
solidBoundary.getIndicator()->signedDistance(PhysR<T, D>(physR)) is negative.

Shouldn’t we be comparing absolute numbers if wall contact is to become active only if the particles are near wall ?

Thanks.

• This topic was modified 7 months, 4 weeks ago by avrachan1.

[jan]

Dear avrachan,

If you only want to define the bottom wall, using IndicInverse<T, DESCRIPTOR::d> is wrong. This indicator inverts the indicator that is passed to it, i.e. thebottomboxyou define is the hole in a wall that takes up the rest of the indicator’s bounding box (defined byminandmax). Therefore, I suggest that you only use theIndicatorCuboid directly:

Vector<T,2> extend_mL(lengthX,L/2.);
Vector<T,2> origin_mL(0.,-L/2);
solidBoundaries.push_back( SolidBoundary<T, DESCRIPTOR::d>(
std::make_unique<IndicatorCuboid2D<T, DESCRIPTOR::d>>(extend_mL,origin_mL), 2, wallContactMaterial));

The error you saw is caused by the particles being completely immersed in the walls. Therefore, a contact is detected and a force is calculated. Full overlap is unrealistic, and the contact model is not valid in this case either, making the force obviously incorrect and likely causing unexpected behavior.

Regarding the second question. You’re absolutely right, you could wrap that part with another if condition that checks if the corresponding location is near a wall (e.g. using the circumference radius of the geometry). However, in most cases the signed distance functions should be fairly cheap to evaluate, and in my tests the contact detection on the lattice wasn’t the limiting factor, the force calculation was usually more expensive, but that might admittedly change if many complex walls are involved.

Best regards,
Jan

• This reply was modified 7 months, 3 weeks ago by jan.

[avrachan1]

Dear Jan,

Thanks for the explanation. I understand it now.

[avrachan1]

Dear Jan,

Why doesn’t the particles obey the periodic boundary condition ?

I can’t seem to find it anywhere in the documentation.

Thanks.

[avrachan1]

I am using

Vector<T,3> periodicity(true,true,false);
ParticleManager<T,DESCRIPTOR,PARTICLETYPE>particleManager(particleSystem,superGeometry,phaseField,converterPF,externalAcceleration,periodicity);

[avrachan1]

I have had to modify the code at two additional places. ( the last two arguments)

processContacts<T, PARTICLETYPE, PARTICLECONTACTTYPE, WALLCONTACTTYPE, ContactProperties<T, 1>>(
particleSystem, solidBoundaries, contactContainer, contactProperties,
superGeometry, contactBoxResolutionPerDirection, static_cast<T>(4. / (3 * util::sqrt(M_PI))),periodicity_func<DESCRIPTOR::d>);

coupleResolvedParticlesToLattice<T, DESCRIPTOR, PARTICLETYPE, PARTICLECONTACTTYPE, WALLCONTACTTYPE>(
particleSystem, contactContainer, superGeometry, phaseField, converterPF, solidBoundaries, periodicity_func<DESCRIPTOR::d>);

periodicity_func is a function which return (true,true,false) vector.

With this the POROSITY fields behave correctly – it obeys the periodic boundary condition, however, the particle positions when accessed through particle.getField<GENERAL,POSITION> does not.

I have gone through the code and I cannot figure out where the particle positions are updated according to the periodic boundary conditions.

Any help in understanding this issue is appreciated.

• This reply was modified 7 months, 2 weeks ago by avrachan1.

[jan]

Dear avrachan,

you’re on the right track already. It’s important that all functions know about the periodicity. The ParticleManager can only pass it to the functions that have access to it. The functions you mentioned don’t know it, hence you have to manually pass it.

The particle position must be explicitly updated when it moves outside the domain. The example doesn’t include that, so I’m sure it’s difficult to find, sorry about that. Currently, we have the following function that should compute the new position:

const PhysR<T,D> position = particles::communication::applyPeriodicityToPosition(
  cuboidGeometry, periodicity, particle.template getField<GENERAL,POSITION>());

Note that you have to update the position afterwards:
particle.template setField<GENERAL, POSITION>(position);

Also, instead of accessing the fields like that, it’s likely better to use the functions provided in src/particles/functions/dataAccessWrappers.h”, e.g., getPosition and setPosition.

Best regards,
Jan

• This reply was modified 7 months, 2 weeks ago by jan.
• This reply was modified 7 months, 2 weeks ago by jan.

[avrachan1]

Dear Jan,

Thank you for your continued support.

If I understand correctly, I have to loop over all particles and update the positions individually.

[jan]

Dear avrachan,

yes, that should be the easiest way for now:

for (std::size_t iP = 0; iP < particleSystem.size(); ++iP) {
  auto particle = particleSystem.get(iP);
  // Add the above function calls here
}

Best regards,
Jan

[Author]

Posts