|
|
momentum.diff (792 bytes) |
|
|
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. |
|
|
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 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. |
|
|
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 |
|
|
splitting.png (8,056 bytes) |
|
|
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. |
|
|
|
|
|
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? |
|
|
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. |
|
|
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. |
|
|
new.diff (2,624 bytes) |
|
|
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. |
- Anonymous
