Wall contact in particle simulation (dkt2d example)
OpenLB – Open Source Lattice Boltzmann Code › Forums › on OpenLB › General Topics › Wall contact in particle simulation (dkt2d example)
- This topic has 8 replies, 2 voices, and was last updated 4 months, 1 week ago by jan.
-
AuthorPosts
-
April 13, 2024 at 11:56 pm #8460avrachan1Participant
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 4 months, 2 weeks ago by avrachan1.
April 16, 2024 at 7:29 pm #8468janParticipantDear 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. thebottombox
you define is the hole in a wall that takes up the rest of the indicator’s bounding box (defined bymin
andmax
). 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 4 months, 2 weeks ago by jan.
April 17, 2024 at 3:08 pm #8476avrachan1ParticipantDear Jan,
Thanks for the explanation. I understand it now.
April 22, 2024 at 7:33 pm #8500avrachan1ParticipantDear Jan,
Why doesn’t the particles obey the periodic boundary condition ?
I can’t seem to find it anywhere in the documentation.
Thanks.
April 22, 2024 at 7:44 pm #8501avrachan1ParticipantI am using
Vector<T,3> periodicity(true,true,false);
ParticleManager<T,DESCRIPTOR,PARTICLETYPE>particleManager(particleSystem,superGeometry,phaseField,converterPF,externalAcceleration,periodicity);April 22, 2024 at 10:28 pm #8502avrachan1ParticipantI 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 4 months, 1 week ago by avrachan1.
April 26, 2024 at 1:05 pm #8528janParticipantDear 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
andsetPosition
.Best regards,
JanApril 26, 2024 at 5:53 pm #8534avrachan1ParticipantDear Jan,
Thank you for your continued support.
If I understand correctly, I have to loop over all particles and update the positions individually.
April 26, 2024 at 7:16 pm #8536janParticipantDear 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 -
AuthorPosts
- You must be logged in to reply to this topic.