View Issue Details

IDProjectCategoryView StatusLast Update
0002666OpenFOAMPatchpublic2017-11-28 11:01
Reportertniemi Assigned Towill  
Status resolvedResolutionfixed 
Product Versiondev 
Fixed in Versiondev 
Summary0002666: Alternative way to calculate momentum transfer from lagrangian parcels
Description(I couldn't add files to, 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, 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.



2017-08-16 15:39


momentum.diff (792 bytes)   
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 =;
momentum.diff (792 bytes)   


2017-08-16 16:12

reporter   ~0008597

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.


2017-08-23 09:51

manager   ~0008624

Last edited: 2017-08-23 09:52

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.


2017-08-23 14:13

reporter   ~0008625

> 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.


2017-08-23 14:48

reporter   ~0008626

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 +*dt


2017-08-25 10:13


splitting.png (8,056 bytes)   
splitting.png (8,056 bytes)   


2017-08-25 10:13

manager   ~0008633

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.


2017-09-01 12:17


derivation.pdf (286,526 bytes)


2017-09-01 12:20

reporter   ~0008651

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

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?


2017-09-25 11:13

manager   ~0008798

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.


2017-11-14 09:46

reporter   ~0009033

Hi, sorry again for very long delay in response, I finally had some time to investigate the changes. I think the code is now quite clean and the splitting approach clearly reduces numerical error in cases where the drag coupling is very strong compared to the integration time step. However, as illustrated in my previous comment, the splitting does not quite work if both coupled and non-coupled terms are important.

To test this, I modified the test case in such a way that I set the intial parcel velocity to zero, carrier velocity to 0.5945 m/s and gravity to -9.81 m/s. This creates a situation, where the gravity and drag are in balance and the parcel does not move. The mass of the parcel is 1 kg and the force is 9.80 N (buoyancy subtracted). By integrating the parcel over 1 s, the momentum transfer should be 9.80 Ns.

Below are the results for different integration time steps (controlled with maxCo):
dt dUTrans(Euler) dUTrans(analytical)
1 0.56 0.59
1e-1 3.70 4.80
1e-2 8.41 9.03
1e-3 9.64 9.72
1e-4 9.78 9.79
1e-5 9.79 9.80

As can be seen, the splitting is not quite correct unless small steps are used. In this case, where u=u0=constant, the explicit method would actually give the correct result and any deviation from it is an error. For the splitting to work accurately, beta*dt << 1. (Here beta is approx 16.5)

This problem also affects heat transfer, as it is quite typical that there can be a reaction source in the parcel and this is balanced by heat transfer out of the parcel. It would be good if the energy was conserved and previously with analytic average the results were exact.

I ran a test with thermoparcels in a similar setup, with the parcels having a fixed energy source of 10 MW/kg. Below are the results for different integration time steps:
dt dHTrans(Euler) dHTrans(analytical)
1 0.303 0.312
1e-1 2.39 2.99
1e-2 7.57 8.56
1e-3 9.69 9.84
1e-4 9.97 9.98
1e-5 9.997 9.998

Again there are sizable errors unless small steps are used.

Based on these experiments I would still propose the original fix, ie. calculate total delta and remove from this the explicit part as alpha_nc*dt. This approach gives exact results for large dt also in these balance cases. The drawback is that it does not allow implicit uncoupled terms, but they are not currently used anyway (especially not for heat transfer). The option 2 in my derivation should also work and enable these terms, but it is more messy and maybe unneccessary.


2017-11-14 11:18


new.diff (2,624 bytes)   
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index e654b56..f97328b 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -181,19 +181,18 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
     // Calculate the integration coefficients
     const vector acp = (Fcp.Sp()*td.Uc() + Fcp.Su())/massEff;
-    const vector ancp = (Fncp.Sp()*td.Uc() + Fncp.Su() + Su)/massEff;
+    const vector ancp = (Fncp.Su() + Su)/massEff;
     const scalar bcp = Fcp.Sp()/massEff;
-    const scalar bncp = Fncp.Sp()/massEff;
+    //const scalar bncp = Fncp.Sp()/massEff;
     // Integrate to find the new parcel velocity
-    const scalar dtEff = cloud.UIntegrator().dtEff(dt, bcp + bncp);
-    const vector deltaUcp = integrationScheme::delta(U_, dtEff, acp, bcp);
-    const vector deltaUncp = integrationScheme::delta(U_, dtEff, ancp, bncp);
+    const scalar dtEff = cloud.UIntegrator().dtEff(dt, bcp);
+    const vector deltaU = integrationScheme::delta(U_, dtEff, acp+ancp, bcp);
     // Calculate the new velocity and the momentum transfer terms
-    vector Unew = U_ + deltaUcp + deltaUncp;
+    vector Unew = U_ + deltaU;
-    dUTrans -= massEff*deltaUcp;
+    dUTrans -= massEff*(deltaU-ancp*dt);
     Spu = dt*Fcp.Sp();
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
index d44526e..3f2d064 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
@@ -287,14 +287,13 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
     // Integrate to find the new parcel temperature
     const scalar dtEff = cloud.TIntegrator().dtEff(dt, bcp);
-    const scalar deltaTcp = integrationScheme::delta(T_, dtEff, acp, bcp);
-    const scalar deltaTncp = integrationScheme::delta(T_, dtEff, ancp, 0);
+    const scalar deltaT = integrationScheme::delta(T_, dtEff, acp+ancp, bcp);
     // Calculate the new temperature and the enthalpy transfer terms
-    scalar Tnew = T_ + deltaTcp + deltaTncp;
+    scalar Tnew = T_ + deltaT;
     Tnew = min(max(Tnew, cloud.constProps().TMin()), cloud.constProps().TMax());
-    dhsTrans -= m*Cp_*deltaTcp;
+    dhsTrans -= m*Cp_*(deltaT-ancp*dt);
     Sph = dt*m*Cp_*bcp;
new.diff (2,624 bytes)   


2017-11-14 11:20

reporter   ~0009035

I uploaded new.diff, which implements option 1 to the current code.


2017-11-28 11:01

manager   ~0009079

I have applied your diff. Commit 9ee3bbc9 in dev. I have tested the velocity case, but not the heat transfer case, as I do not have the files. Please let me know if this does not fix the heat transfer.

I have also corrected (I think) the splitting, but at present it is not being used. It is not the same as your derivation. Your derivation, as I understand it, functions in terms of average velocities. That's fine, but that then makes the issue one of how those averages are computed. In the analytical method, it is relatively straightforward in that an integral average is appropriate. For the Euler method, it is not so clear. What is the average velocity over a step integrated with a Euler scheme?

Another way to think about it is to do the analysis with the combined coefficients to get a profile of velocity over the step. This can then be plugged into the expression for the coupled force only, and integrated to get the coupled velocity difference. That's what I've implemented. The necessary calls are in KinematicParcel::calcVelocity, but they are at present commented in favour of your change. As far as I can tell, the results are the same. I will leave this functionality available but unused in case additional implicit coefficients are dreamed up in future. If you get the chance to test it, that would be great.

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
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
2017-11-14 09:46 tniemi Status resolved => feedback
2017-11-14 09:46 tniemi Resolution fixed => reopened
2017-11-14 09:46 tniemi Note Added: 0009033
2017-11-14 11:18 tniemi File Added: new.diff
2017-11-14 11:20 tniemi Note Added: 0009035
2017-11-14 11:20 tniemi Status feedback => assigned
2017-11-28 11:01 will Status assigned => resolved
2017-11-28 11:01 will Resolution reopened => fixed
2017-11-28 11:01 will Note Added: 0009079