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

0001563  OpenFOAM  Bug  public  20150310 15:45  20150313 14:41 
Reporter  Assigned To  henry  
Priority  normal  Severity  minor  Reproducibility  N/A 
Status  closed  Resolution  fixed  
Platform  GNU/Linux  OS  Ubuntu  OS Version  14.10 
Summary  0001563: Incorrect momentum predictor for pimple in DPMFoam  
Description  In the DPMFoam solver, UcEqn.H, line 2133 it states: if (pimple.momentumPredictor()) { solve ( UcEqn == fvc::reconstruct ( phicForces/rAUcf  fvc::snGrad(p)*mesh.magSf() ) ); } The right handside should be multiplied by the gas phase fraction, alphacf, i.e. phicForces/rAUcf  fvc::snGrad(p)*mesh.magSf() should be: alphacf * (phicForces/rAUcf  fvc::snGrad(p)*mesh.magSf())  
Additional Information  A valid reference for the momentum equation is found at N.G. Deen, M. Van Sint Annaland, M.A. Van der Hoef, J.A.M. Kuipers, Review of discrete particle modeling of fluidized beds, Chemical Engineering Science, Volume 62, Issues 1–2, January 2007, Pages 2844, ISSN 00092509, http://dx.doi.org/10.1016/j.ces.2006.08.014. (http://www.sciencedirect.com/science/article/pii/S0009250906004830) When looking at equation 8, it can be seen that both the pressure and gravity term need to be multiplied by the gas phase fraction, alphacf.  
Tags  No tags attached.  

In all the references we worked from both for DPM and MPPIC the force and pressure terms are NOT multiplied by the gas phase fraction. If you want to change the code so that the modelling exactly corresponds to the paper you cite you MUST also make the correcsponding changes to the formulation of the pressure equation and the momentum corrector. 

When looking at the MPPIC paper you are correct both pressure and force terms are not multiplied by the gas phase fraction. However, in all MPPIC papers gravity IS multiplied by the gas phase fraction. In OpenFOAM this is not the case, as gravity (g) is directly included in phicForces. 

Agreed, I will change it to multiply gravity by the gas phasefraction. 

Resolved by commit 5c68a98799e67e3f6680f013050d132036ffe232 

Reopened to allow further discussion 

Indeed in original MPPIC papers it seems that only the gravity term is multiplied by gas phase fraction, but in my opinion this is not correct and the equations are inconsistent. For instance if we have a static, packed bed (eg. alphac=0.5, uc=up=0), the gas phase momentum equation would reduce to grad p = alphac*rho_f*g = 0.5*rho_f*g, ie. the hydrostatic pressure would be 50 % too low. In order to have a correct hydrostatic pressure, also the pressure gradient should be multiplied by alphac. Another inconsistency in the MPPIC papers is that the drag term appears to be WenYu drag law, which is divided with alphac. Again to be consistent, also this term should be multiplied with alphac. So I wouldn't put too much trust on those papers. If I read the Deen et al. correctly, it seems to be consistent. However, in OpenFOAM there is currently no pressure gradient term in the solid particle equations, which should be added if the other terms are multiplied with alphac. 

Dear tniemi, According to the original MPPIC paper: Andrews, M.J. and O'Rourke, P.J. (1996). The Multiphase ParticleinCell (MPPIC) Method for Dense Particle Flows. International Journal of Multiphase Flow, 22(2):379–402. It explains in the paragraph under equation 10 why the pressure is not multiplied by the gas phase fraction: it is partially included in the interphase momentum transfer term. Combining the pressure in the interphase momentum transfer term and the expression in the gas momentum equation would yield the correct pressure multiplied by the gas phase fraction. If you are stating that the interphase particle momentum transfer term in OpenFOAM does not include the pressure gradient as in equation 10 in the above mentioned paper, then you are indeed correct. However, I haven't come around checking the particle equations yet. 

Dear kenyi89, Thank you for pointing that out, now the papers make much more sense. I had missed the fact that part of the pressure gradient is "hidden" in the momentum transfer term. 

Ok, I have now read the original MPPIC paper and the Deen et al. papers more carefully and I have reached the following conclusions: Both papers have equivalent gas phase momentum equations (at least regarding gravity and pressure gradient). The difference is that Deen et al. directly writes alphac*grad p in the gas phase momentum equation, while in the original paper the F term contains (1alphac)*grad p. When summed together F+grad p we get the alphac*grad p term also in the original paper. Before the applied patch, the OF implementation was different than in the references. After the patch we are still missing the (1alphac)*grad p term. Multiplying the grad p term with alphac will work in the gas phase momentum equation, but the (1alphac)*grad p term should also be added to the solid particle forces. kenyi89, do you agree with my conclusions? In addition to above conclusions, I also want to point out, that if we do the changes to the gas phase momentum equation and add the (1alphac)*grad p force to particles, we also need to modify some of the drag models. At least the ErgunWenYuDragForce and WenYuDragForce should be also multiplied with alphac, as the existing versions assume that there is no (1alphac)*grad p term. I have attached modified versions of ErgunWenYuDragForce.C and WenYuDragForce.C. I also attached a short pdf, where I try to explain where the difference comes from. I hope that it does not have huge mistakes, as it just my own notes about the issue. 

WenYuDragForce.C (3,022 bytes)
/**\ =========  \\ / F ield  OpenFOAM: The Open Source CFD Toolbox \\ / O peration  \\ / A nd  Copyright (C) 20132014 OpenFOAM Foundation \\/ M anipulation   License This file is part of OpenFOAM. OpenFOAM is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. OpenFOAM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. \**/ #include "WenYuDragForce.H" #include "volFields.H" // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // template<class CloudType> Foam::scalar Foam::WenYuDragForce<CloudType>::CdRe(const scalar Re) const { if (Re > 1000.0) { return 0.44*Re; } else { return 24.0*(1.0 + 0.15*pow(Re, 0.687)); } } // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template<class CloudType> Foam::WenYuDragForce<CloudType>::WenYuDragForce ( CloudType& owner, const fvMesh& mesh, const dictionary& dict ) : ParticleForce<CloudType>(owner, mesh, dict, typeName, true), alphac_ ( this>mesh().template lookupObject<volScalarField> ( this>coeffs().lookup("alphac") ) ) {} template<class CloudType> Foam::WenYuDragForce<CloudType>::WenYuDragForce ( const WenYuDragForce<CloudType>& df ) : ParticleForce<CloudType>(df), alphac_ ( this>mesh().template lookupObject<volScalarField> ( this>coeffs().lookup("alphac") ) ) {} // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // template<class CloudType> Foam::WenYuDragForce<CloudType>::~WenYuDragForce() {} // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template<class CloudType> Foam::forceSuSp Foam::WenYuDragForce<CloudType>::calcCoupled ( const typename CloudType::parcelType& p, const scalar dt, const scalar mass, const scalar Re, const scalar muc ) const { scalar alphac(alphac_[p.cell()]); return forceSuSp ( vector::zero, (mass/p.rho()) *0.75*CdRe(alphac*Re)*muc*pow(alphac, 2.65)/(sqr(p.d())) ); } // ************************************************************************* // 

ErgunWenYuDragForce.C (3,326 bytes)
/**\ =========  \\ / F ield  OpenFOAM: The Open Source CFD Toolbox \\ / O peration  \\ / A nd  Copyright (C) 20132014 OpenFOAM Foundation \\/ M anipulation   License This file is part of OpenFOAM. OpenFOAM is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. OpenFOAM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. \**/ #include "ErgunWenYuDragForce.H" #include "volFields.H" // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // template<class CloudType> Foam::scalar Foam::ErgunWenYuDragForce<CloudType>::CdRe ( const scalar Re ) const { if (Re > 1000.0) { return 0.44*Re; } else { return 24.0*(1.0 + 0.15*pow(Re, 0.687)); } } // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template<class CloudType> Foam::ErgunWenYuDragForce<CloudType>::ErgunWenYuDragForce ( CloudType& owner, const fvMesh& mesh, const dictionary& dict ) : ParticleForce<CloudType>(owner, mesh, dict, typeName, true), alphac_ ( this>mesh().template lookupObject<volScalarField> ( this>coeffs().lookup("alphac") ) ) {} template<class CloudType> Foam::ErgunWenYuDragForce<CloudType>::ErgunWenYuDragForce ( const ErgunWenYuDragForce<CloudType>& df ) : ParticleForce<CloudType>(df), alphac_ ( this>mesh().template lookupObject<volScalarField> ( this>coeffs().lookup("alphac") ) ) {} // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // template<class CloudType> Foam::ErgunWenYuDragForce<CloudType>::~ErgunWenYuDragForce() {} // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template<class CloudType> Foam::forceSuSp Foam::ErgunWenYuDragForce<CloudType>::calcCoupled ( const typename CloudType::parcelType& p, const scalar dt, const scalar mass, const scalar Re, const scalar muc ) const { scalar alphac(alphac_[p.cell()]); if (alphac < 0.8) { return forceSuSp ( vector::zero, (mass/p.rho()) *(150.0*(1.0  alphac)/alphac + 1.75*Re)*muc/(sqr(p.d())) ); } else { return forceSuSp ( vector::zero, (mass/p.rho()) *0.75*CdRe(alphac*Re)*muc*pow(alphac, 2.65)/(sqr(p.d())) ); } } // ************************************************************************* // 



Dear tniemi, I agree with your conclusions. I don't have a full grasp yet on the kinematics of the particle cloud, but it seems the interphase momentum transfer term is NOT including (1alphac)*grad p and it should. @Henry, can you confirm that the interphase momentum transfer term in OpenFOAM is missing the pressure gradient contribution? I.e. the second term in the integral in equation 10 in the original MPPIC paper as stated above. Regarding the drag, I have not looked into that yet in great detail, I'm just on my way of learning OpenFOAM. If I look at your notes, it seems that with the proper (1alphac)*grad p term added, the system as defined in OpenFOAM should be equal to equation (1) in your notes (since grad p in the momentum equation and (1alphac)*grad p in the momentum transfer equation should produce alphac*grad p as equation (1) in your report states). that would mean using K_gs and not \tilde{K_gs}. 

The Kinematic lagrangian library includes the PressureGradientForce which can be applied to DPM or MPPIC. The PressureGradientForce class infers the continuousphase pressure gradient from DDt(Uc) rather than using the gradient of the pressure field directly. However, the direct method could be implemented and Timo has provided a quick hack I am looking into. I can make all the changes discussed here to OpenFOAMdev provided you can validate them at least for the tutorial cases supplied and also ensure that the stability of the solvers are not compromised. 

Okay, so to be clear phicForces in DPMFoam does include the pressure gradient force making the equations consistent? Regarding implementing the direct method I'm willing to have a look once it's implemented, although what would be the benefit of implementing the direct method over the way it's done now? 

> Okay, so to be clear phicForces in DPMFoam does include the pressure gradient > force making the equations consistent? The particlephase pressure gradient force is runtime selected. It is not currently selected in the tutorials. So in the original implementation the pressure gradient and buoyancy were consisent in that they would generate the correct hydrostatic pressure distribution. Now that the phasefraction is in the bouyancy term but not in the pressure gradient term it is inconsistent. So either the change must be reverted or the implementation changes to have the phasefraction in both terms and the pressure gradient forces selected for the particles. > Regarding implementing the direct method I'm willing to have a look once it's > implemented You can of course study the implemention but mor importantly it need to be validated by the cases being run and the results compared to the previous formulation and the results in the papers. > although what would be the benefit of implementing the direct method over the > way it's done now? I do not know, would you like to compare them? 

Dear Henry > So either the change must be reverted or the implementation changes to have the phasefraction in both terms and the pressure gradient forces selected for the particles. Reverting the change means being inconsistent with the equations in the MPPIC papers, so that would not be my preferred choice. I would opt for keeping the change in the buoyancy term, NOT adding the phase fraction in the pressure term in the gas momentum equation and selecting the pressure gradient force for the particles, this is consistent with the equations as these two pressure terms together will reproduce alphac*grad p as explained in the original MPPIC paper. I think this is the least confusing as it follows exactly the MPPIC papers as well as standard approaches for DPM. 

I just found out a good discussion about the different DPM approaches in the PhD thesis by Leboreiro, which is fortunately open access http://gradworks.umi.com/33/03/3303878.html (link to pdf at bottom). A very relevant part to our discussion can be found in section "2.5.1 Decomposition of the Fluid–Solid Interaction Force", where the two different decomposition approaches are presented. The the first one is the same as in the original MPPIC paper and the second one is the same as in OF before the multiplication of the gravity. According to Leboreiro, the two approaches should be mathematically the same, if the DDt(Uc) pressure gradient is selected for parcels. (And gravity has (rho_srho_f) and correct, divided by alphac drag models are used). However, he states that the first approach is more common. > although what would be the benefit of implementing the direct method over the > way it's done now? This is unclear to me and would need testing as Henry said. In theory, as Leboreiro states, the approaches should give the same results, but in practice maybe not? If you are referring whether to compute the "pressure gradient term" as DDt(Uc) or as (1alphac)*grad P then there is a difference. They both can be called "pressure gradient terms", but they are not the same thing. (1alphac)*grad P is used in the first decomposition approach and DDt(Uc) in the second. So if we move to use the first decomposition approach, I think we need to implement (1alphac)*grad P term. > Reverting the change means being inconsistent with the equations in the MPPIC papers, so that would not be my preferred choice. I would opt for keeping the change in the buoyancy term, NOT adding the phase fraction in the pressure term in the gas momentum equation and selecting the pressure gradient force for the particles, this is consistent with the equations as these two pressure terms together will reproduce alphac*grad p as explained in the original MPPIC paper. I think this is the least confusing as it follows exactly the MPPIC papers as well as standard approaches for DPM. On the other hand I would prefer Henry's approach, ie. multiplying the pressure gradient, as it seems a bit counter intuitive to add and subtract the same gradient. Also Deen et al. presents the equations in this form, so it is somewhat a matter of taste I guess. 

Dear tniemi, Excellent reference, thank you! To give my opinion on this, I would then prefer multiplying the pressure gradient in the momentum equation by alphac, that is the first option in the thesis, following equation 2.59 and 2.61. As an extra reason I would like to add that the viscosity part of the stress is already multiplied by alphac, which is the divDevRhoReff(Uc) term in UcEqn. If I understand things correctly moving to the second option would also mean removing the alphac factor in this term. We should make sure though that equation 2.70 is satisfied in this case, which describes the particle velocity updates for the first option. So to conclude, I would vote for the first approach, is this is most intuitive (and avoids adjusting gravity and adding DDt(Uc) terms to the particle forces). The only outstanding confustion for me is that the stress (S_f) is differently defined in equations 2.61 and 2.70, since in the latter the viscous part seems to be multiplied by alpha/epsilon already. 

Sorry, my claim about divDevRhoReff(Uc) is not true, it is not multiplied by alphac. Then the first option would mean also multiplying divDevRhoReff(Uc) by alphac to be consistent. 

> The only outstanding confustion for me is that the stress (S_f) is differently defined in equations 2.61 and 2.70, since in the latter the viscous part seems to be multiplied by alpha/epsilon already. Yes, that is odd. However, as the viscous component is anyhow neglected, it shouldn't matter much. I also can't recall any other reference, where there was a gas viscosity tensor in the particle/solid phase equations. 

I think it is inconsistent to include the viscuous stess in the UcEqn, but then neglect it when treating the particle velocity update. I would also like to add that in the thesis it seems that epsilon is constant, in which case his derivations for option 2 makes sense. The derivations are not that straightforward if epsilon has to be included in the material derivative. Tniemi, do you agree option 1 is the most logical one given what is already in OpenFOAM? So this would include:  Multiplying pressure, viscuous stress and gravity with alpha in the momentum equation.  Multiplying pressure with alpha in the momentum predictor and pEqn.  Multiplying pEqn.flux() with alpha when calculating the new fluxes and velocity.  Making sure particle velocity updates are according to equation 2.70, including the viscuous stress term (and excluding buoyancy and DDt(Uc) terms) 

My understanding is that the original implementation was correct if the particle pressure gradient force is included as currently defined. What you are proposing is a significant change and will require careful testing and validation which is very timeconsuming. If you would like to implement and test the alternative strategy outlined above and provide a set of patches I would be happy to apply them if the testing is successful. Otherwise I would need to secure funding for this development. In the meantime I propose to revert the change to the buoyancyterm so that the implementation corresponds to the second option in Leboreiro as described by Timo above. 

>I would also like to add that in the thesis it seems that epsilon is constant, in which case his derivations for option 2 makes sense. The derivations are not that straightforward if epsilon has to be included in the material derivative. It may be so. This is also the reason, why I would lean towards option 1, as I'm not 100 % sure that including the implemented "DDt(Uc)"term makes the two options exactly equal. However, I also don't know what are the practical differences, ie. does "DDt(Uc)"term really matter and how easy it is to implement the (1alphac)*grad P term for instance. In my experince gradient terms tend to produce problems... > Tniemi, do you agree option 1 is the most logical one given what is already in OpenFOAM? In addition to previous, option 1 is more logical in the sense, that other MPPIC references seem to favour it. But as Henry said, the modifications require work and non trivial amount of testing. >  Multiplying pressure, viscuous stress and gravity with alpha in the momentum equation. >  Multiplying pressure with alpha in the momentum predictor and pEqn. > Multiplying pEqn.flux() with alpha when calculating the new fluxes and velocity. I guess so. However, I'm always nervous when dealing with the pressure correction equations, so I don't want to say anything too definitive :) > I think it is inconsistent to include the viscuous stess in the UcEqn, but then neglect it when treating the particle velocity update. >  Making sure particle velocity updates are according to equation 2.70, including the viscuous stress term (and excluding buoyancy and DDt(Uc) terms) Velocity update should have the (1alphac)*grad P term, but I'm not convinced about the viscous stess. As I stated in my previous comment, I don't recall other papers, which do it like that. Also if it was included, then we would need to further modify the drag expressions. In EulerianEulerian approach only the pressure gradient is separated, not the viscous stress, and the traditional, empirical drag expressions rely on that assumption. > In the meantime I propose to revert the change to the buoyancyterm so that the implementation corresponds to the second option in Leboreiro as described by Timo above. Yes, if we don't have enough muscle to do all the changes now, we should return to the previous, stable state. 

Ok I did some calculations and I think what Henry said checks out. If we revert the change and modify the interaction term to include: F = <drag_term> + rho_f*V_p*DDt(Uc) (per particle) the expressions should be consistent. Implementing the grad P term instead will lead to all the changes I listed previously (option 1 in Leboreiro), and hence is not feasible right now indeed, even though this is the implementation method described in the MPPIC papers... My only current (big) concern is that particle velocities should be updated according to 2.71, i.e. it should include the (rho_f  rho_p)*g term and DDt(Uc) terms, and not fluid stress and explicit pressure gradients. 

> My only current (big) concern is that particle velocities should be updated according to 2.71, i.e. it should include the (rho_f  rho_p)*g term and DDt(Uc) terms, and not fluid stress and explicit pressure gradients. It already works like this. If the user turns on "PressureGradientForce", it will add the "DDt(Uc)" term. PressureGradientForce is a "calcCoupled"force, so it affects both phases. The Gravityforce is also written as mass*g_*(1.0  p.rhoc()/p.rho()) = mass/p.rho()*g_*(p.rho()  p.rhoc()) There are no fluid stress nor explicit pressure gradients in particle equations in the present implementation. 

Ok great, then I suppose reverting is the solution. The PressureGradientForce should be added though in the tutorials related to DPM/MPPIC, as otherwise we have inconsistencies. 

Change reverted. 
Date Modified  Username  Field  Change 

20150310 15:45 

New Issue  
20150310 15:59  henry  Note Added: 0004069  
20150310 15:59  henry  Status  new => closed 
20150310 15:59  henry  Assigned To  => henry 
20150310 15:59  henry  Resolution  open => no change required 
20150310 16:35 

Note Added: 0004071  
20150310 16:35 

Status  closed => feedback 
20150310 16:35 

Resolution  no change required => reopened 
20150310 16:48  henry  Note Added: 0004072  
20150311 10:53  henry  Note Added: 0004079  
20150311 10:53  henry  Status  feedback => resolved 
20150311 10:53  henry  Resolution  reopened => fixed 
20150311 13:09  henry  Note Added: 0004083  
20150311 13:09  henry  Status  resolved => feedback 
20150311 13:09  henry  Resolution  fixed => reopened 
20150311 15:27  tniemi  Note Added: 0004085  
20150311 16:14 

Note Added: 0004086  
20150311 16:14 

Status  feedback => assigned 
20150311 17:07  tniemi  Note Added: 0004088  
20150312 13:12  tniemi  Note Added: 0004091  
20150312 13:14  tniemi  File Added: WenYuDragForce.C  
20150312 13:14  tniemi  File Added: ErgunWenYuDragForce.C  
20150312 13:19  tniemi  File Added: drag_explanation.pdf  
20150312 14:47 

Note Added: 0004092  
20150312 16:21  henry  Note Added: 0004093  
20150312 17:07 

Note Added: 0004097  
20150312 17:14  henry  Note Added: 0004098  
20150313 07:52 

Note Added: 0004100  
20150313 09:05  tniemi  Note Added: 0004101  
20150313 10:21 

Note Added: 0004102  
20150313 10:33 

Note Added: 0004103  
20150313 11:10  tniemi  Note Added: 0004104  
20150313 11:25 

Note Added: 0004105  
20150313 11:53  henry  Note Added: 0004106  
20150313 12:08  tniemi  Note Added: 0004109  
20150313 12:11  tniemi  Note Edited: 0004109  
20150313 12:15  tniemi  Note Edited: 0004109  
20150313 12:41 

Note Added: 0004111  
20150313 13:11  tniemi  Note Added: 0004113  
20150313 13:19 

Note Added: 0004114  
20150313 14:40  henry  Note Added: 0004115  
20150313 14:40  henry  Status  assigned => closed 
20150313 14:41  henry  Resolution  reopened => fixed 
20150324 00:17  liuhuafei  Issue cloned: 0001614 