View Issue Details

Jump to NotesJump to History
IDProjectCategoryView StatusDate SubmittedLast Update
0000560OpenFOAMBugpublic2012-06-19 11:522018-07-10 11:25
Reportertniemi Assigned Tohenry  
PrioritynormalSeverityminorReproducibilityalways
Status resolvedResolutionfixed 
Summary0000560: Problems with kinetic theory models
DescriptionThere 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_s-2/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 g0prime-functions have errors.

I have attached a file, which contains corrected g0Prime-functions for all the radialModels.

[1] Gidaspow, D. Multiphase Flow and Fluidization - Continuum and Kinetic Theory
Descriptions. Academic Press, 1994. ISBN 0-12-282470-9.

[2] Syamlal, M., Rogers, W., and O’Brien, T. J. MFIX Documentation: Volume 1,
Theory Guide. National Technical Information Service, Springfield, VA, 1993.
DOE/METC-9411004, NTIS/DE9400087
Tagskinetic theory

Activities

tniemi

2012-06-19 11:52

reporter  

radialModel_fixed.zip (29,414 bytes)

henry

2012-06-19 14:08

manager   ~0001398

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

albertop

2012-06-19 20:05

reporter   ~0001400

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 two-phase flow theory applied to fluidization, Int. J. Multiphase Flow, Vol. 22, pp 21-66, 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.

henry

2012-06-19 21:33

manager   ~0001401

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

albertop

2012-06-19 23:04

reporter   ~0001402

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

henry

2012-06-20 09:33

manager   ~0001403

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?

tniemi

2012-06-20 10:03

reporter   ~0001404

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.

albertop

2012-06-20 19:38

reporter   ~0001406

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 sub-models (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.

henry

2012-06-20 23:07

manager   ~0001407

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

albertop

2012-06-21 02:15

reporter   ~0001408

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 non-standard form.

henry

2012-06-21 07:43

manager   ~0001409

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.

albertop

2012-06-21 08:40

reporter   ~0001410

OK.

The inconsistency exists because of the two-fluid 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 Boltzmann-Enskog 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.

henry

2012-06-21 10:03

manager   ~0001411

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.

henry

2012-06-21 16:47

manager   ~0001413

The radial models are updated in commit dbf4cf0b0efd27c96a8e9cddf44ebe663ea163bd

albertop

2012-06-21 19:34

reporter   ~0001414

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

henry

2013-09-02 13:16

manager   ~0002448

All these issues are now resolved in our development line and will be included in the next release.

Issue History

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