View Issue Details
| ID | Project | Category | View Status | Date Submitted | Last Update |
|---|---|---|---|---|---|
| 0002282 | OpenFOAM | Bug | public | 2016-10-04 13:43 | 2017-08-16 10:44 |
| Reporter | tniemi | Assigned To | will | ||
| Priority | low | Severity | minor | Reproducibility | always |
| Status | resolved | Resolution | fixed | ||
| Summary | 0002282: Lagrangian parcels, small issue in setCellValues | ||||
| Description | In KinematicParcel.C, line 291, the current cell of a parcel is cached in case the parcel moves to another cell during the current time step. The cached cell index is then passed on to setCellValues and calc functions. However, in setCellValues the carrier phase properties, such as density and velocity, are interpolated to the new position of the parcel and the cached index is ignored. If the parcel has moved to another cell and eg. cell interpolation is used, the carrier phase values are taken incorrectly from the new cell instead of the old cell. I would guess that this issue is probably not very significant in practice and I only noticed it by chance. However, it is slightly inconsistent that the sources to the carrier phase are applied to the old cell, but the carrier phase properties are taken from the new cell. If the properties differ much between the two cells, the source will be inaccurate. The issue could be fixed by caching also the particle position and tetIndices before movement and passing these to setCellValues. This way the carrier phase properties would be interpolated to the same cell that the sources are applied to. The same issue affects also the other parcels that implement setCellValues, such as ThermoParcel and ReactingParcel. | ||||
| Tags | No tags attached. | ||||
|
|
Can you provide a patch to resolve this issue? |
|
|
|
|
|
While playing around I noticed another, perhaps a more serious bug that affects the parcels when they hit patches, such as walls and outlets. The patchInteraction-models are called during the trackToFace call and the interaction models immediately alter the parcel velocity. This velocity is then used to calculate the momentum sources. What this means is that when eg. a parcel flies towards a wall, the velocity that is used to calculate the momentum sources is the velocity after the bounce, ie. velocity away from the wall. Also if the boundary has an escape condition, the parcel velocity is set to zero in the patchInteraction model and a large momentum source can get created if the parcel velocity was high. I have attached an icoUncoupledKinematicParcelFoam test case that demonstrates both issues. The case has 3 cells: the velocity of the first cell is 5 m/s and the velocity of the other cells is 6 m/s. A 1 kg parcel is injected into the first cell with a velocity of 5 m/s. What should happen is that the momentum source in the first cell should be zero, because there is no slip velocity. Then in the second cell the parcel should accelerate to ~5.9 m/s due to drag and the source should be about -0.9 kgm/s. In the final cell there should be about 0.1 kgm/s momentum transfer left. However, when run with maxCo=0.1, the momentum sources are -0.273738,-0.77907 and -1.59833. In the first cell momentum transfer happens, because the velocity is interpolated wrongly in between the cells. In the last cell momentum transfer occurs due to escape boundary condition which sets the velocity to zero. I can take a look about making a patch to the interpolation issue when I have a bit more time, because it should be simple to fix. However, the patch interaction issue can be more difficult to fix. |
|
|
To clarify the patch interaction issue, what should happen when a parcel is going to hit a patch: 1. Move particle from it's initial position to the patch with its initial velocity U0 2. Calculate the velocity change during the travel from the initial position to the patch. Calculate momentum transfer. 3. Apply patch interaction, ie. in rebound boundary flip the new velocity computed in step 2. What currently happens is: 1. Move particle from it's initial position to the patch with its initial velocity U0. 2. Apply patch interaction, eg. flip U0 or set it to zero. 3. Calculate the velocity change to modified U0 using deltaT from step 1. Update momentum transfer based on this computation. |
|
|
I agree with your assessment. |
|
|
Hmm, properly fixing the interpolation issue is not as straightforward as I thought. Besides in setCellValues, intepolation is also used in several other places, such as in some force models and in thermoParcel calcSurfaceValues. As such, many function signatures would need to be changed in order to pass around the cached tetIs and position. I guess a more elegant fix would be to have a different type of trackToFace-function, which would not move the parcel, but just would compute the dt. The actual movement (and patch interaction) could then be done after all the other computations. |
|
|
What would the "different type of trackToFace-function" do if not track to the face? What is the "dt" needed for? |
|
|
What I meant is that it would be handy if the trackToFace-function would not alter the data of the parcel so that the position() and currentTetIndices() etc. would still refer to the original values in calc and setCellValues functions. This way it would not be necessary to worry about using cached values in the various places. Maybe the trackToFace could create a separate data structure that holds the necessary tracking info, such as the final position of the parcel, dt it took to travel there and a possible patch that was hit at the end? Then we could use the dt in calc and setCellValues and after that move the parcel to its final position and optionally call hit patch. |
|
|
Now that the the lagrangian tracking functions are very clean, I was wondering if you have had time to consider the issues presented here? As the hitFace and trackToFace are now separate functions, it looks to me that the both issues could be solved by simply doing the calc phase in between trackToFace and hitFace? I have added a simple patch that does this and it seems to work in simple test cases. I haven't yet tested it much and there may be side-effects that I'm missing. I could come up with two corner-cases: 1. Reacting particle gets used up during calc -> check keepParticle before hitFace 2. Velocity of the particle has turned away from the wall -> at least most rebound interactions already have a check for this. What do you think about this fix, could it work or are there some other corner cases/fundamental problems with it that I'm not taking into account? With cellPoint interpolation the calc still uses the end position, but I would guess it is as good an approximation as using the start position. |
|
|
patch.diff (1,027 bytes)
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index b0c5418..bf56198 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -291,8 +291,8 @@ bool Foam::KinematicParcel<ParcelType>::move
f = min(f, maxCo*l/max(SMALL*l, mag(s)));
if (p.active())
{
- // Track to and hit the next face
- p.trackToAndHitFace(f*s, f, td);
+ // Track to the next face
+ p.trackToFace(f*s, f);
}
else
{
@@ -332,6 +332,10 @@ bool Foam::KinematicParcel<ParcelType>::move
if (p.onFace())
{
+ if (td.keepParticle)
+ {
+ p.hitFace(f*s, td);
+ }
td.cloud().functions().postFace(p, p.face(), td.keepParticle);
}
|
|
|
I think the change makes a lot of sense. I'm not sure corner case #1 is even an issue. I can't envisage a case in which a face/patch interaction with a to-be-deleted parcel is problematic. A little wasteful maybe, but I can't see it failing. Rebound interactions would need to be careful in the case that the velocity has changed, I agree, and I think most of them are already. There's also an ambiguity for stick and escape; if the velocity has changed so that the particle is now leaving the patch, should it still stick or escape? I would tend towards the opinion that if rebounds do nothing when the particle is leaving the patch, then the other interactions should follow suit. If we're doing this, then we should probably get rid of the celli arguments to calc and similar in order to remove any ambiguity. These methods should get the cell index from the particle. |
|
|
> I think the change makes a lot of sense. I'm not sure corner case #1 is even an issue. I can't envisage a case in which a face/patch interaction with a to-be-deleted parcel is problematic. A little wasteful maybe, but I can't see it failing. Yes, it is probably a hypothetical concern for most cases. I was mainly thinking more complex models like ThermoSurfaceFilm-model, where maybe something bad could happen. However, I haven't ever used the model, so I'm not sure if there is any reason to be concerned. However, in ThermoSurfaceFilm, in contrast to other rebound interactions, there is no checks for the direction of the velocity. > I would tend towards the opinion that if rebounds do nothing when the particle is leaving the patch, then the other interactions should follow suit. Yes, I think it would be logical. If the velocity is away from the wall, then it is likely that if smaller steps would have been taken, the particle wouldn't touch the patch at all. One solution could be to have a common pre-check that patch interaction is only applied if the (relative) velocity is still towards the wall. |
|
|
I've run a smattering of big cases with the change and haven't had any errors, so I'm happy to push this to dev. Commit e7fc53c0. I've kept the change to the bare minimum for now. Your example case now records y-momentum sources of 0, -1.05063, and -0.0882667, which I think is more in keeping with what you were expecting. I haven't modified any of the patch interactions, as I don't think there's really any chance of failure. Even on a rebound (of whatever sort, thermoSurfaceFilm or otherwise), if the direction was wrong, it would simply result in two hits and the parcel carrying along on it's original trajectory. I think this is fine for the moment. If we discover significant inconsistencies in the interactions stemming from the fact that the parcel now has the newly calculated velocity, then we can fix them as they arise. Thank you for the report. This was a significant inconsistency in the way in which parcels were evolved, and it's a good thing to have resolved. |
| Date Modified | Username | Field | Change |
|---|---|---|---|
| 2016-10-04 13:43 | tniemi | New Issue | |
| 2016-10-04 14:48 | henry | Note Added: 0006959 | |
| 2016-10-05 12:42 | tniemi | File Added: testcase.tar.gz | |
| 2016-10-05 12:56 | tniemi | Note Added: 0006961 | |
| 2016-10-05 14:11 | tniemi | Note Added: 0006962 | |
| 2016-10-05 14:23 | henry | Note Added: 0006963 | |
| 2016-10-05 17:58 | tniemi | Note Added: 0006965 | |
| 2016-10-05 18:09 | henry | Note Added: 0006966 | |
| 2016-10-06 07:22 | tniemi | Note Added: 0006972 | |
| 2017-08-14 07:37 | tniemi | Note Added: 0008568 | |
| 2017-08-14 07:37 | tniemi | File Added: patch.diff | |
| 2017-08-14 11:01 | will | Note Added: 0008572 | |
| 2017-08-14 11:31 | tniemi | Note Added: 0008573 | |
| 2017-08-16 10:44 | will | Assigned To | => will |
| 2017-08-16 10:44 | will | Status | new => resolved |
| 2017-08-16 10:44 | will | Resolution | open => fixed |
| 2017-08-16 10:44 | will | Note Added: 0008593 |
