
[Sponsors] 
April 12, 2010, 04:40 
Hiting the wall by particle

#1 
Member
Mihails Ščepanskis
Join Date: Jul 2009
Posts: 36
Rep Power: 9 
Hello!
I have a problem with the particles tracking. The problem is to simulate the hit of the particle to the wall (the boundary of the mesh) correctly. There is solidParticle library in the OpenFOAM. The function trackToFace in solidParticle.C should take into account the size of the particle in the case of collision. But if the cell near the wall is smaller then the size of the particle doesn't work! In this case particle can approach the wall closer than the radius of the particle. Can you tell me how can I solve this problem?
__________________
Mihails Ščepanskis Laboratory for Mathematical Modelling of Environmental and Technological Processes University of Latvia 

April 12, 2010, 04:57 

#2 
Super Moderator
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 22 
yup, that is correct.
I think this is a deliberate limitation, because what you want to do is very expensive. Also, particles can only belong to 1 cell at a time, so how do you treat this? This should also be considered, but again, extremely expensive to solve. N 

April 13, 2010, 03:10 

#3 
Senior Member
matej forman
Join Date: Mar 2009
Location: Brno, Czech Republic
Posts: 115
Rep Power: 10 
in more detail, what Niklas is mentioning  the lagrangian particle tracking implemented in FOAM is sitting on the assumption, that the particle is smaller than the cell. Otherwise you would need to integrate e.g. the forces on particle around the particle surface as it expands over more cells. New formulation of coupling to the continuous phase would be needed. The continuous phase would not be continuous any more  force transfer across particle.... You would need to use VOF, or 6DOF or some other treatment...
matej 

April 13, 2010, 03:47 

#4 
Member
Mihails Ščepanskis
Join Date: Jul 2009
Posts: 36
Rep Power: 9 
Niklas and Matej, thank you very much for answer!
Matej, you are right that in the case of big particle we should take into account a lot of forces (for example lift force, that takes into account the circulation of the flow around the particle). I have already implemented this and some other forces in the solidParticle library. The most problem is in the region close to the wall. We should make the size of the mesh elements near the wall much smaller than in the middle zone to simulate the hydrodynamic correctly, especially in the case of EM governed flow (my case). Thus appears the described problem of the collision with the wall. The way how to solve this problem that I want to suggest is so. I propose to make a virtual wall in the distance of the radius of the particle from the real wall. Then make the size of the particle zero for the function wallImpactDistance. Thus we will indicate the collision when the center of the particle reach this virtual wall. Now I can not understand how to implement the check of the collision with this virtual wall in the function trackToFace. If you can give me some advise in this question, I'll be very thankful!
__________________
Mihails Ščepanskis Laboratory for Mathematical Modelling of Environmental and Technological Processes University of Latvia Last edited by shchepan; April 14, 2010 at 09:49. 

April 13, 2010, 09:12 

#5 
Super Moderator
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 22 
This is a really difficult problem.
Let me outline how I would do it and where I see the challenges. Maybe you have a better idea, in which case I would love to hear it. Since we are assuming that wallinteraction can occur even though we are not in a cell close to a wall, we have to move our wallfaces inwards according to the wallImpactDistance and check every wallface before we move the particle to see if we will hit a wall face first. This operation is pretty straightforward: we just add 1 extra check to the original algorithm before we move the particle. Something like this  1. original algorithm. 2. check wall faces. 3. move particle. However...point 2 in the list above is trickier than one would think Maybe a simple example can serve as an illustration to the problem. Assume we have 2 walls as below, that we move inwards and the particle wants to move from A to B. Code:
______________ A > __B________ A > B / 2 / 1 / ====> / / outside / / / Hence, when you check for wallinteraction you have to also check if the intersectionpoint is within the wallface. And, since we have to take into account the size of the particle and move the wallfaces, we are creating small holes in the mesh through which the particle can escape. So in short, wallinteraction has to take into account if the interacting point is within the face (very expensive check) and it's also important when you move the wallfaces that you dont create holes in the mesh. This implies that you cant just move the faces, you have to move the points that constitutes the faces. I hope it is clear now why I find this a bit difficult. N 

April 13, 2010, 12:54 
Query

#6 
New Member
Sumit Sinha
Join Date: Mar 2009
Location: Champaign, ILLINOIS, U.S.A
Posts: 25
Rep Power: 9 
Hey Guys,
May be you all can help me, I have been trying to play with Lagrangian particle solver, I had two questions in this regard a) Is there some restriction on this solver, when we run it in parallel. b) I have even tried running it on my desktop, every thing keeps on going well for certain time steps but eventually it hangs with message on the screen as follows Time = 0.41 Courant Number mean: 0.0123124 max: 0.234255 Moving Particles 53 Particles moved. 3 walls hit. 0 particles left the model. DILUPBiCG: Solving for Ux, Initial residual = 0.000270292, Final residual = 8.83653e08, No Iterations 1 DILUPBiCG: Solving for Uy, Initial residual = 0.000453004, Final residual = 9.93462e07, No Iterations 1 DILUPBiCG: Solving for Uz, Initial residual = 0.00126948, Final residual = 1.39543e06, No Iterations 1 DICPCG: Solving for p, Initial residual = 0.0139312, Final residual = 0.00119803, No Iterations 3 DICPCG: Solving for p, Initial residual = 0.00276161, Final residual = 0.00026989, No Iterations 56 time step continuity errors : sum local = 6.67705e08, global = 5.95145e09, cumulative = 8.18682e07 DICPCG: Solving for p, Initial residual = 0.00173553, Final residual = 0.000171299, No Iterations 8 DICPCG: Solving for p, Initial residual = 0.00126418, Final residual = 0.000115941, No Iterations 8 time step continuity errors : sum local = 2.8673e08, global = 6.02904e09, cumulative = 8.12653e07 ExecutionTime = 597.83 s ClockTime = 598 s Time = 0.42 Courant Number mean: 0.0123111 max: 0.23442 Moving Particles It just does't proceed beyond this All other residuals seems to be under control and courant number is well below zero. I am also attaching the file cloudProperties, Any help is enormously appreciated. Thanks, Sumit 

April 14, 2010, 02:24 

#7 
Member
Mihails Ščepanskis
Join Date: Jul 2009
Posts: 36
Rep Power: 9 
Niklas,
you are right, there are a lot of problems in your example of 2 walls and obtuse angle between them. Let us initially consider simpler case. Let we have some simple geometry, for example cylinder (there are no obtuse angles). In this case we can build the virtual wall in the distance of particle radius from the real wall. You propose to apply some function to indicate heating this virtual wall before the move function. But now algorithm check heating with the wall inside of move function. There are trackToFace function that indicate the collision and move the particle only to the face of cell or wall. Thus we should implement this new check into the trackToFace function, but not outside of the move function, as you suggest. Am I right? Sumit, where are some problems with parallel calculation when you are working with the big cloud (huge number of particles). I have not solve this problem yet.
__________________
Mihails Ščepanskis Laboratory for Mathematical Modelling of Environmental and Technological Processes University of Latvia 

April 14, 2010, 03:01 

#8  
Super Moderator
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 22 
Quote:
it will all be in the trackToFace routine. It has to be. trackToFace checks all faces in the cell to see if the particle can move to its end position or if it should move the particle to the intersecting face. What I meant was that before we move it (which is actually done inside the trackToFace routine) we have to add an additional check for all, and I do mean all, wall faces before we can go ahead and move it. And I really dont recommend that you implement something that will work only on a special type of geometry. 

April 14, 2010, 04:35 

#9 
Member
Mihails Ščepanskis
Join Date: Jul 2009
Posts: 36
Rep Power: 9 
Niklas,
I am very thankful for your help and ideas! But I want to ask you some more technical questions. There is such code in Particle.C: 00307 // change cell 00308 if (internalFace) // Internal face 00309 { ... 00326 } 00327 else 00328 { 00329 ParticleType& p = static_cast<ParticleType&>(*this); ... 00341 if (!p.hitPatch(patch, td, patchi)) 00342 { ... 00371 else if (isA<wallPolyPatch>(patch)) 00372 { 00373 p.hitWallPatch 00374 ( 00375 static_cast<const wallPolyPatch&>(patch), td 00376 ); 00377 } ... 00382 } 00383 } Here is indicated the collision with the wall and applied the hitWallPatch method. But I cannot understand how to add the new virtual wall to the list of faces and to put flag that it is the wall, not cell face. Can you help me in this question? As a next step I'll think how to patch up the possible holes in the virtual wall.
__________________
Mihails Ščepanskis Laboratory for Mathematical Modelling of Environmental and Technological Processes University of Latvia 

April 14, 2010, 05:00 

#10 
Super Moderator
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 22 
That part of the code is just responsible for the walltreatment and will remain unchanged.
What you want to do is to edit the Particle.Cfile and the trackToFaceroutine. just stick an includefile there and add the necessary checks and the rest should work, since boundaryfaces are moved according wallImpactDistance. Just remove the check in the lambdaroutine that wallImpactDistance is shorter than cellsize also of course. Hitting a face is signalled with facei_ not being equal to 1. So when you have moved the particle to the virtual wall, setting facei_ accordingly, the tracking will think it is on a wall and apply the boundary treatment. Code:
template<class ParticleType> template<class TrackData> Foam::scalar Foam::Particle<ParticleType>::trackToFace ( const vector& endPosition, TrackData& td ) { const polyMesh& mesh = cloud_.polyMesh_; DynamicList<label>& faces = cloud_.labels_; findFaces(endPosition, faces); // check for all boundary faces here and if a face has been hit // add it to the faces variable, the rest of the code should remain unmodified #include "myWallCheck.H" facei_ = 1; scalar trackFraction = 0.0; if (faces.empty()) // inside cell { trackFraction = 1.0; position_ = endPosition; } else // hit face { scalar lambdaMin = GREAT; 

April 14, 2010, 08:58 

#11 
Member
Mihails Ščepanskis
Join Date: Jul 2009
Posts: 36
Rep Power: 9 
Niklas,
I understand, that in myWallCheck.H I should check the wallImpactDistance in direction of normal vector of face for all faces. If this wallImpactDistance for the faces[j] is less than d_/2.0 then we should redefine this faces[j] as a boundary face. And all other algorithm will work. But 1) how can I get the normal vector of faces[j]? 2) how can I mark this faces[j] as a boundary face? Thank you a lot, Niklas! P.S. I think I have found answer on the 1st question: mesh.Sf()[facei_]/mesh.magSf()[facei_] Am I right?
__________________
Mihails Ščepanskis Laboratory for Mathematical Modelling of Environmental and Technological Processes University of Latvia Last edited by shchepan; April 14, 2010 at 09:20. 

April 14, 2010, 09:55 

#12  
Super Moderator
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 22 
Quote:
You dont have to mark a face as a boundary face. FOAM knows if the face is a internal face or not. You just have to loop over every patch and check if its a wall or not, then you can just grab the faces of an entire patch and loop over those. not sure of the exact syntax, but to grab them should be something like this const surfaceVectorField& Sfbi = mesh.Sf().boundaryField()[patchi]; or const faceList& bcfaces = mesh.cells().boundaryMesh()[patchi] alternatively mesh.boundaryMesh()[patchi].cells() (not sure how it is) Also, remember that cells() returns the indeces to the faces, while faces() returns the indeces to the points. if it doesnt work try looking in the wallFunctions directory for the turbulenceModels should be something useful there If you cant figure it out I will try to find out the exact expression. 

April 14, 2010, 10:28 

#13  
Member
Mihails Ščepanskis
Join Date: Jul 2009
Posts: 36
Rep Power: 9 
Sorry, Niklas, I don't understand your algorithm:
Quote:
What I proposed: Quote:
Niklas, would you be so kind to explain me the algorithm you suggest and/or estimate my one?
__________________
Mihails Ščepanskis Laboratory for Mathematical Modelling of Environmental and Technological Processes University of Latvia 

April 15, 2010, 04:55 

#14 
Super Moderator
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 22 
I am away at a conference at the moment, but when I get back
I will write you a piece of code that I think will make it clear. N 

April 15, 2010, 05:25 

#15 
Member
Mihails Ščepanskis
Join Date: Jul 2009
Posts: 36
Rep Power: 9 
Thank you very much, Niklas! I'll be waiting for code.
Have a good time on the conference!
__________________
Mihails Ščepanskis Laboratory for Mathematical Modelling of Environmental and Technological Processes University of Latvia 

April 15, 2010, 15:12 

#16 
Super Moderator
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 22 
On a train back home, lets hope the shaky connection works
Here's how you access all the stuff you need. I'm leaving it up to you to fill out the details, i.e, how you check if the facehit is within the faceedges myWallCheck.H Code:
const polyBoundaryMesh& patches = mesh.boundaryMesh(); const pointField& points = mesh.points(); forAll(patches, i) { const polyPatch& patchi = patches[i]; if (isA<wallPolyPatch>(patchi)) { const vectorField& Sfbc = patches[i].faceAreas(); const vectorField& Cfbc = patches[i].faceCentres(); const label nStartFaceIndex = patches[i].start(); forAll(patchi, j) { const labelList& facePoints = patchi[j]; //const vector& Cf = Cfbc[j]; //vector Sf = Sfbc[j]; //Sf /= mag(Sf); label facei = nStartFaceIndex + j; // lambda already accounts for wallImpactDistance scalar lam = lambda(position_, endPosition, facei, stepFraction_); // check if trajectory is intersected by face if ( (lam <= 1.0) && (lam >= 0.0)) { // face intersects trajectory, now we need to check if the face hit // is withing the face vector pHit = position_ + lam*(endPosition  position_); forAll(facePoints, fi) { label fi1 = (fi + 1) % facePoints.size(); vector p1 = points[facePoints[fi1]]; vector p0 = points[facePoints[fi]]; vector edge = p1  p0; // check if position is within all edges... } } } } } 

April 16, 2010, 03:00 

#17 
Member
Mihails Ščepanskis
Join Date: Jul 2009
Posts: 36
Rep Power: 9 
Niklas, thank you. But I want to ask you to explain some places in your code to understand it:
myWallCheck.H Code:
const polyBoundaryMesh& patches = mesh.boundaryMesh(); const pointField& points = mesh.points(); forAll(patches, i) { const polyPatch& patchi = patches[i]; //1) What does this check mean? if (isA<wallPolyPatch>(patchi)) { //2) You don't use this variables further, do you? const vectorField& Sfbc = patches[i].faceAreas(); const vectorField& Cfbc = patches[i].faceCentres(); const label nStartFaceIndex = patches[i].start(); forAll(patchi, j) { const labelList& facePoints = patchi[j]; //const vector& Cf = Cfbc[j]; //vector Sf = Sfbc[j]; //Sf /= mag(Sf); label facei = nStartFaceIndex + j; // lambda already accounts for wallImpactDistance scalar lam = lambda(position_, endPosition, facei, stepFraction_); // check if trajectory is intersected by face if ( (lam <= 1.0) && (lam >= 0.0)) { // face intersects trajectory, now we need to check if the face hit // is withing the face //3) This is the vector of the coordinate of the point where the particle intersects the face, isn't it? vector pHit = position_ + lam*(endPosition  position_); forAll(facePoints, fi) { label fi1 = (fi + 1) % facePoints.size(); vector p1 = points[facePoints[fi1]]; vector p0 = points[facePoints[fi]]; vector edge = p1  p0; // check if position is within all edges... //4) Here I should check if pHit is in the face. But it will be //within anyway, because we have checked that the //particle intersects this face! } } } } } 1) to remove the check that the distance between the cell center and face center is larger than wallImpactDistance in lambda function in ParticleI.H 2) to add this code in Particle.C After this modification will been made, FOAM will check if the particle hit the wall anyway, but not only when the particle is in the boundary cell. Have I understood your idea correct?
__________________
Mihails Ščepanskis Laboratory for Mathematical Modelling of Environmental and Technological Processes University of Latvia 

April 16, 2010, 03:14 

#18 
Super Moderator
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 22 
1. I thought that was obvious. It checks if the patch is a wall or not.
2. You will (probably ) need them later when you move the face centre point, to check if the face hit is withing the faceedges 3. yes. 4. nope. the lambda calculation assumes that the face is infinite, since it only has a face centre point and a face normal. I would suggest that you decompose the face into triangles using the facecenter and the points on the edges and then check if the face hit is within one of these triangles. general. 1. yes. 2. yes. and if the face is hit you need to append the face index to the facesvariable in the trackToFaceroutine 

April 16, 2010, 07:33 

#19 
Member
Mihails Ščepanskis
Join Date: Jul 2009
Posts: 36
Rep Power: 9 
I have written the code to check if pHit is within the edges of the face.
Code:
vector pHit = position_ + lam*(endPosition  position_); scalar S = 0.0; forAll(facePoints, fi) { label fi1 = (fi + 1) % facePoints.size(); vector p1 = points[facePoints[fi1]]; vector p0 = points[facePoints[fi]]; vector lHit1 = p1  pHit; vector lHit2 = p0  pHit; scalar S += mag(lHit1 ^ lHit2) / 2.0; } if (S == mag(Sfbc)) { //the face should be added to facesvariable }
__________________
Mihails Ščepanskis Laboratory for Mathematical Modelling of Environmental and Technological Processes University of Latvia 

April 16, 2010, 07:50 

#20 
Super Moderator
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 22 
faces.append(facei)
...I seriously doubt that check though. I think you should come up with something else 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Fluent3DMeshToFoam  simvun  OpenFOAM Other Meshers: ICEM, Star, Ansys, Pointwise, GridPro, Ansa, ...  48  May 14, 2012 05:20 
DPM UDF particle position using the macro P_POS(p)[i]  dm2747  FLUENT  0  April 17, 2009 01:29 
Errors running allwmake in OpenFOAM141dev with WM_COMPILE_OPTION%3ddebug  unoder  OpenFOAM Installation  11  January 30, 2008 21:30 
Multicomponent fluid  Andrea  CFX  2  October 11, 2004 05:12 
Iclude particle diameter effect in wall?  yueroo  FLUENT  0  April 14, 2001 03:56 