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

0000560  OpenFOAM  Bug  public  20120619 11:52  20180710 11:25 
Reporter  tniemi  Assigned To  henry  
Priority  normal  Severity  minor  Reproducibility  always 
Status  resolved  Resolution  fixed  
Summary  0000560: Problems with kinetic theory models  
Description  There seems to be a couple of problems with the kinetic theory models used in eg. compressibleTwoPhaseEulerFoam solver. The first problem is that the solid volume fraction is included too many times in the kinetic theory viscosity terms. In literature there are a couple of ways to write the solid phase stress contribution in the multiphase momentum equations. For example Gidaspow[1,pp. 314] and Syamlal et al.[2] use a formulation, where the stresses are written as div(tau_s), where tau_s=P_s*I+2*mu_s*D+(lambda_s2/3*mu_s)div(v)I. With this notation, for example the bulk viscosity by Lun et al. is lambda_s=4/3*alpha_s^2*rho_s*d_s*g0*(1+e)*sqrt(Theta/pi) and it is implemented in kineticTheoryModel.C in this exact form. However, in the UEgns.H file in compressibleTwoPhaseEulerFoam the stress contribution is written as div(alpha_s*tau_s), which means that in the final equation the alpha_s will have a power of 3 instead of 2 in the bulk viscosity expression. This same thing also affects the shear viscosity terms mu_s and the frictional terms for all the implemented models. I know that compressibleTwoPhaseEulerFoam is just a preliminary release, but I think this problem also affects the earlier solver twoPhaseEulerFoam. There the alpha does not appear in the div terms, but this is because one power of alpha has been divided away. In this case the final power of alpha should be 1, but in the code it is 2. The other problem with the kinetic theory models is related to the radialModels. The radialModels have to implement a function called g0prime, which should calculate the derivative of g0 with respect to the volume fraction. However, the derivatives are calculated incorrectly. For example, the SinclairJacksonRadial gives g0=1.0/(1.0  pow(alpha/alphaMax, 1.0/3.0)) and in this case dg0/dalpha = (1.0/3.0)*pow(alpha/alphaMax, 2.0/3.0)/(alphaMax*sqr(1.0  pow(alpha/alphaMax, 1.0/3.0))). However in code it is currently (1.0/3.0)*pow(alpha/alphaMax, 2.0/3.0)/(alphaMax*sqr(1.0  pow(alpha/alphaMax, 1.0/3.0))) which is not correct (wrong sign). Besides the SinclairJackson, also the other g0primefunctions have errors. I have attached a file, which contains corrected g0Primefunctions for all the radialModels. [1] Gidaspow, D. Multiphase Flow and Fluidization  Continuum and Kinetic Theory Descriptions. Academic Press, 1994. ISBN 0122824709. [2] Syamlal, M., Rogers, W., and O’Brien, T. J. MFIX Documentation: Volume 1, Theory Guide. National Technical Information Service, Springfield, VA, 1993. DOE/METC9411004, NTIS/DE9400087  
Tags  kinetic theory  



Note that twoPhaseEulerFoam is in phase intensive form wheras compressibleTwoPhaseEulerFoam is in conservative form so in the latter the stress contribution should be of the form div(alpha_s*tau_s). Could you start by checing the form of tau_s in twoPhaseEulerFoam and it should be of the same form in compressibleTwoPhaseEulerFoam if the stress term is written div(alpha_s*tau_s). 

Hello, the momentum equation is correct in both the codes. However the stress tensor has to be defined differently if you use the kinetic theory, because what Timo calls tau_s is actually alpha*tau_s (See Enwald, H., Peirano, E., Almstedt, A.E., Eulerian twophase flow theory applied to fluidization, Int. J. Multiphase Flow, Vol. 22, pp 2166, 1996 for a summary). For what concerns the radial distribution functions, I would also suggest to drop the one from Gidaspow. It is a modification of the Sinclair model, however it is not correct from a theoretical point of view, since when alpha = 0, g0 must be 1, and it does not happen with Gidaspow's distribution. 

Just to clarify: the code kinetic theory code should be changed so that the current tau_s is divided by alpha correct? 

Hi Henry, the code kineticTheoryModel.C is correct: tau_s has the correct expression to be used in the equation for Theta. The problem is in UEqns.H, because both nua and lambda have been already multiplied by alpha, so Ra is indeed already alpha*Ra when the kinetic theory is turned on. Probably the easiest way to fix this in twoPhaseEulerFoam is to divide nua by alpha when nuEffa is calculated (line 11), and do the same with lambda when its contribution is added to Rca (line 26). 

It would be good for stability reasons if we could avoid dividing by alpha unnecessarily, would it be possible to reformulate kineticTheory.mua() and kineticTheory.lambda() to provide these already "divided by alpha" as their names would suggest? 

Yes, it should be relatively straightforward to modify the mua and lambda terms. However, this would mean many modifications, because all the different models would have to be modified. (including frictional terms) Also the Theta equation would have to be modified, because currently it is correct. Still, in the end this might be the most correct approach. On the other hand, if the division approach is used, one possible place to put the division is in the end of the kineticTheoryModel.C after the Theta equation is calculated. 

A possible workaround to avoid divisions by alpha in the case of compressibleTwoPhaseEulerFoam would be to make a small change in UEqns.H. We can define the U1Eqn without the viscous stresses, and then add a check: if (kineticTheory.on()) { U1Eqn +=  fvm::laplacian(nuEff1, U1) + fvc::div(Rc1); } else { U1Eqn +=  fvm::laplacian(alpha1*nuEff1, U1) + fvc::div(alpha1*Rc1); } This would not require any change in the submodels (except a correction to the frictional viscosity, which must be multiplied by alpha before adding it to mua in kineticTheoryModel.C). Unfortunately this does not resolve the problem in twoPhaseEulerFoam. However, dividing by max(alpha, 0.0001) should not cause troubles, since it's done already elsewhere in more critical terms like ppMagf. 

> Yes, it should be relatively straightforward to modify the mua and lambda > terms. However, this would mean many modifications, because all the different > models would have to be modified. (including frictional terms) Also the Theta > equation would have to be modified, because currently it is correct. Still, in > the end this might be the most correct approach. I agree, I think it would be better if these terms are evaluated without alpha as this is more physical and consistent with the other stress terms. > On the other hand, if the division approach is used, one possible place to > put the division is in the end of the kineticTheoryModel.C after the Theta > equation is calculated. But this is still a hack and I would not wish to see OpenFOAM degenerate into layers of interdependent hacks. 

Unfortunately you cannot avoid divisions completely: there is no explicit power of alpha in the expressions of the frictional terms. Also, even though it's a minor issue, modifying the expressions from lambda and mua would make the implementation inconsistent with respect to what is found in the literature. The code would have the granular conductivity expressed as in the literature, but the other quantities in a nonstandard form. 

Inconsistency with the literature is not as bad as inconsintency with the rest of multiphase continuum mechancics. I will look into this further when I get some time to do so. 

OK. The inconsistency exists because of the twofluid model and the kinetic theory model are derived from different starting points. The formulation used by authors working with the kinetic theory is fully consistent with the derivation of the transport equations from the BoltzmannEnskog equation, which leads to a momentum equation with div(tau_s), then modeled in analogy to what done for fluids. However, most of the difficulties in making the code consistent affect only one specific solver, twoPhaseEulerFoam, due to the form of the momentum equation it implements. In compressibleTwoPhaseEulerFoam the fix is trivial, and does not not require divisions or modifications to the kinetic theory library (see above). Also, I noticed that the term depending on Ct is present in the stress tensor volTensorField Rc1 ( "Rc1", ((2.0/3.0)*I)*(sqr(Ct)*k + nuEff1*tr(gradU1T))  nuEff1*gradU1T ); but it is not included in nuEffa, when the kinetic theory is used. Also this should be made consistent. 

The Ct modelling was created origininally for bubbly flows; are you sure it is appropriate for particulate flows? I am not even convinced it is appropriate for bubbly flows and we plan to replace it with a second set of turbulence property equations which will be more appropriate for real systems which almost universally "suffer" from phase inversion. 

The radial models are updated in commit dbf4cf0b0efd27c96a8e9cddf44ebe663ea163bd 

I agree with you about the Ct modelling. I don't think it is appropriate for particulate flows, where the effect of fluid turbulence typically is accounted for into the equation for Theta through a model for the fluidparticle velocity correlation. I have implemented some algebraic model to do this, which I can contribute. For more complex models, based on transport equations, most of the work was done in Prof. Olivier Simonin's group in Tolouse. Based on this, I would drop the term containing Ct for particulate flows, and I agree with you about implementing an alternative turbulence model for bubbly flows too. 

All these issues are now resolved in our development line and will be included in the next release. 
Date Modified  Username  Field  Change 

20120619 11:52  tniemi  New Issue  
20120619 11:52  tniemi  File Added: radialModel_fixed.zip  
20120619 14:08  henry  Note Added: 0001398  
20120619 20:05  albertop  Note Added: 0001400  
20120619 21:33  henry  Note Added: 0001401  
20120619 23:04  albertop  Note Added: 0001402  
20120620 09:33  henry  Note Added: 0001403  
20120620 10:03  tniemi  Note Added: 0001404  
20120620 19:38  albertop  Note Added: 0001406  
20120620 23:07  henry  Note Added: 0001407  
20120621 02:15  albertop  Note Added: 0001408  
20120621 07:43  henry  Note Added: 0001409  
20120621 08:40  albertop  Note Added: 0001410  
20120621 10:03  henry  Note Added: 0001411  
20120621 16:47  henry  Note Added: 0001413  
20120621 19:34  albertop  Note Added: 0001414  
20130902 13:16  henry  Note Added: 0002448  
20130902 13:16  henry  Status  new => resolved 
20130902 13:16  henry  Resolution  open => fixed 
20130902 13:16  henry  Assigned To  => henry 
20170803 13:47  Juho  Tag Attached: kinetic theory  
20170803 13:48  Juho  Tag Attached: kineticTheory  
20180710 11:25  administrator  Tag Detached: kineticTheory 