2017-10-22 00:11 BST

View Issue Details Jump to Notes ]
IDProjectCategoryView StatusLast Update
0002666OpenFOAMPatchpublic2017-09-25 11:13
Reportertniemi 
Assigned Towill 
PrioritylowSeveritytweakReproducibilitysometimes
StatusresolvedResolutionfixed 
Product Versiondev 
Target VersionFixed in Versiondev 
Summary0002666: Alternative way to calculate momentum transfer from lagrangian parcels
Description(I couldn't add files to https://bugs.openfoam.org/view.php?id=2282, so I opened this)

The problem:
Currently, the momentum transfer in KinematicParcel.C is calculated as dUTrans=dt*(Feff.Sp()*(Ures.average() - Uc_)), where Ures.average() is the mean of the particle velocity during the time step. With Euler integration Ures.average() is just the linear average of the velocities before and after the integration. This is fine if the velocity change has been roughly linear, ie. the integration time step is smaller than the response time of a parcel. On the other hand, if the velocity change is not linear, the error can be very large.
In the testcase in https://bugs.openfoam.org/view.php?id=2282, the error in momentum transfer is ~13% with 150 micron particle and ~0.02s step size. If the particle size is changed to 1 micron, the momentum source is 2385.25.

Proposed fix:
Calculate the momentum of the parcel before and after the integration. The momentum sink to continuous phase is then the difference between these two values when the contribution of possible non-coupled forces is subtracted.
TagsNo tags attached.
Attached Files
  • diff file icon momentum.diff (792 bytes) 2017-08-16 15:39 -
    diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
    index 572cd97..df10890 100644
    --- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
    +++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
    @@ -185,8 +185,7 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
     
         vector Unew = Ures.value();
     
    -    // note: Feff.Sp() and Fc.Sp() must be the same
    -    dUTrans += dt*(Feff.Sp()*(Ures.average() - Uc_) - Fcp.Su());
    +    dUTrans -= (Unew-U_)*massEff-dt*(Fncp.Su()+Su);
     
         // Apply correction to velocity and dUTrans for reduced-D cases
         const polyMesh& mesh = td.cloud().pMesh();
    
    diff file icon momentum.diff (792 bytes) 2017-08-16 15:39 +
  • png file icon splitting.png (8,056 bytes) 2017-08-25 10:13 -
    png file icon splitting.png (8,056 bytes) 2017-08-25 10:13 +
  • pdf file icon derivation.pdf (286,526 bytes) 2017-09-01 12:17

-Relationships
+Relationships

-Notes

~0008597

tniemi (reporter)

Note: using analytical integration will help a lot, it works well up to 1e-12 parcel size or so in the test case. However, I guess that Euler is quite often used for U and the behavior may be surprising.

~0008624

will (manager)

Last edited: 2017-08-23 09:52

View 2 revisions

Sorry for the delay here. I think what you are saying makes sense. I don't think the average of the velocity over the integration step multiplied by the force coefficient is a particularly good way of calculating the momentum change.

Beyond that, the average velocity over the Euler integration step is not (u_0 + u_1)/2. Euler has the velocity profile with time as asymptotic, so the average (as determined by integration) is a logarithmic function. This clearly removes the advantage of Euler integration as the "cheap" option. Taking the difference is a much better idea.

The complication is that the uncoupled forces should not be included, so your solution subtracts them from the final momentum difference. This is probably OK, so long as there aren't any implicit (Sp term) uncoupled forces. It is currently assumed that this is the case. If there are implicit uncoupled forces, then I don't know how to isolate coupled momentum change from uncoupled.

Also, it is worth noting that this mechanism is in use for energy, too. I suppose all the same points apply.

~0008625

tniemi (reporter)

> If there are implicit uncoupled forces, then I don't know how to isolate coupled momentum change from uncoupled.

Yes, I don't know if there is a robust solution for this without resorting to using average velocity. But as you said, currently there are no models using non-coupled Sp and this is also assumed in the current implementation.

> Also, it is worth noting that this mechanism is in use for energy, too. I suppose all the same points apply.

Yes, probably so. I haven't been that concerned with heat transfer as I have always used analytical for that, which works ok because it gives the correct average. I would guess that analytical is a more common and reasonable choice here, because with a constant htc it gives the exact solution. This makes the possible performance penalty more justified.

Also, the problems might not be that bad with heat transfer. With drag, if there is eg. turbulent dispersion, the slip velocity can be temporarily quite large. Even with not so small parcels this can lead to a huge coupling factor/small timescales. The heat transfer coeff tends to react more mildly to slip velocity.

~0008626

tniemi (reporter)

Actually, now that I thought about it, it might be possible to calculate the momentum transfer without knowing the velocity. If one also removes the contribution of coupled Su, the remaining momentum change has to be explained by the Sp terms. By using the ratio of uncoupled/coupled Sp one could calculate the coupled portion:

coupled mom = Fcp.Sp/(Fcp.Sp+Fncp.Sp)*(mom without Fcp.su) + Fcp.su*dt

~0008633

will (manager)

I don't understand where that expression has come from. Could you elaborate?

I have an alternative proposal. All the integration routines can be rewritten as a fully explicit contribution multiplied by a factor, F, that depends only on the implicit coefficient. The explicit contribution can be trivially split into it's coupled and non-coupled components. See the uploaded image for the equations.

I can't justify the split physically (i.e., why the explicit contribution splits, but F does not), but it's straightforward, conservative, and it tends to the correct form in all the limiting cases I can think of.

~0008651

tniemi (reporter)

Sorry for the long delay.

I have attached a derivation for my proposal. Case 1 corresponds to original suggestion, case 2 for nonzero non-coupled Sp.

You proposal also looks reasonable for splitting the change in velocity to coupled and non-coupled parts, but I'm not sure how to use it to get the momentum transfer. Maybe I'm just thinking this wrong, but let's consider the following:

Assume that there is a particle in gas with a drag force and gravity and that the slip velocity is such that drag and gravity are in balance. Then physically delta up should be 0 and the momentum source should be equal to Sp_c*(ug-up)*dt (where ug=gas velocity, up=particle velocity).

Now in this case ac=Sp_c*ug and bc=Sp_c

Using the split formulation and Euler,
du_c=(Sp_c*(ug-u0))F=(Sp_c(*ug-u0))*[dt/(1+Sp_c* dt)]

If dt or Sp_c is large, then dt/(1+Sp_c* dt) = 1/Sp_c
and
du_c=(Sp_c(ug-u0))/Sp_c=(ug-u0)

This tells correctly that given infinite time, the coupled force tries to undo the velocity difference. However, how to get the momentum source from this?

~0008798

will (manager)

I've pushed some changes to dev which I think resolves this. Commits a5806207 and fa3ed3b8.

The transfer terms are now calculated from the change in particle property, not from evaluating the rate expressions with an average value over the integration step. The integration methods have been rewritten to allow for the split between coupled and non-coupled, and to take advantage of the simplification of not calculating the step average.

The Euler scheme is, as far as I can tell, now fully conservative. The difference between UTrans values between Euler and analytical schemes in the #2282 test case is also not as big as before. This has made me confident enough to push. I'm aware that we didn't fully conclude how to handle the coupled/non-coupled split, though. Please reopen the report if you think the calculation still isn't right.
+Notes

-Issue History
Date Modified Username Field Change
2017-08-16 15:39 tniemi New Issue
2017-08-16 15:39 tniemi File Added: momentum.diff
2017-08-16 16:12 tniemi Note Added: 0008597
2017-08-23 09:51 will Note Added: 0008624
2017-08-23 09:52 will Note Edited: 0008624 View Revisions
2017-08-23 14:13 tniemi Note Added: 0008625
2017-08-23 14:48 tniemi Note Added: 0008626
2017-08-25 10:13 will File Added: splitting.png
2017-08-25 10:13 will Note Added: 0008633
2017-09-01 12:17 tniemi File Added: derivation.pdf
2017-09-01 12:20 tniemi Note Added: 0008651
2017-09-25 11:13 will Assigned To => will
2017-09-25 11:13 will Status new => resolved
2017-09-25 11:13 will Resolution open => fixed
2017-09-25 11:13 will Fixed in Version => dev
2017-09-25 11:13 will Note Added: 0008798
+Issue History