View Issue Details
ID  Project  Category  View Status  Date Submitted  Last Update 

0002666  OpenFOAM  Patch  public  20170816 15:39  20171128 11:01 
Reporter  tniemi  Assigned To  will  
Priority  low  Severity  tweak  Reproducibility  sometimes 
Status  resolved  Resolution  fixed  
Product Version  dev  
Fixed in Version  dev  
Summary  0002666: 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 noncoupled forces is subtracted.  
Tags  No tags attached.  

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 = (UnewU_)*massEffdt*(Fncp.Su()+Su); // Apply correction to velocity and dUTrans for reducedD cases const polyMesh& mesh = td.cloud().pMesh(); 

Note: using analytical integration will help a lot, it works well up to 1e12 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. 

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. 

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

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 



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



Sorry for the long delay. I have attached a derivation for my proposal. Case 1 corresponds to original suggestion, case 2 for nonzero noncoupled Sp. You proposal also looks reasonable for splitting the change in velocity to coupled and noncoupled 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*(ugup)*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*(ugu0))F=(Sp_c(*ugu0))*[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(ugu0))/Sp_c=(ugu0) 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? 

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 noncoupled, 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/noncoupled split, though. Please reopen the report if you think the calculation still isn't right. 

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 noncoupled 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 1e1 3.70 4.80 1e2 8.41 9.03 1e3 9.64 9.72 1e4 9.78 9.79 1e5 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 1e1 2.39 2.99 1e2 7.57 8.56 1e3 9.69 9.84 1e4 9.97 9.98 1e5 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. 

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*(deltaUancp*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_*(deltaTancp*dt); Sph = dt*m*Cp_*bcp; 

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

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. 
Date Modified  Username  Field  Change 

20170816 15:39  tniemi  New Issue  
20170816 15:39  tniemi  File Added: momentum.diff  
20170816 16:12  tniemi  Note Added: 0008597  
20170823 09:51  will  Note Added: 0008624  
20170823 09:52  will  Note Edited: 0008624  
20170823 14:13  tniemi  Note Added: 0008625  
20170823 14:48  tniemi  Note Added: 0008626  
20170825 10:13  will  File Added: splitting.png  
20170825 10:13  will  Note Added: 0008633  
20170901 12:17  tniemi  File Added: derivation.pdf  
20170901 12:20  tniemi  Note Added: 0008651  
20170925 11:13  will  Assigned To  => will 
20170925 11:13  will  Status  new => resolved 
20170925 11:13  will  Resolution  open => fixed 
20170925 11:13  will  Fixed in Version  => dev 
20170925 11:13  will  Note Added: 0008798  
20171114 09:46  tniemi  Status  resolved => feedback 
20171114 09:46  tniemi  Resolution  fixed => reopened 
20171114 09:46  tniemi  Note Added: 0009033  
20171114 11:18  tniemi  File Added: new.diff  
20171114 11:20  tniemi  Note Added: 0009035  
20171114 11:20  tniemi  Status  feedback => assigned 
20171128 11:01  will  Status  assigned => resolved 
20171128 11:01  will  Resolution  reopened => fixed 
20171128 11:01  will  Note Added: 0009079 