[cheng1988sjtu]

Hi,

I had a chance to check what's going on in the alphaEqn.H in twoPhaseEulerFoam, and the basic idea is that they introduce a particle-particle stress to avoid higher volumetric concentration than alphaMax. and this artificial particle-particle stress does appear in the pEqn.H to modify the velocity field in particle phase directly.

However, the author of this code try to modify phia in the alphaEqn.H, well I don't see it necessary because the velocity has already been modified in pEqn.H.

Good thing is that, the rUaAf in alphaEqn.H is zero always ( you can print out the rUaAf value, and because you are calling rUaAf before it's defined in pEqn.H ) and thus ppMagf and phipp in alphaEqn.H is useless in alphaEqn.H, is that correct ?

Any opinions are welcome, Thank you!

Zhen

[alberto]

The code is correct: the idea behind the procedure is to directly include the effect of the particle pressure in the volume fraction equation. You probably noticed that the equation for alpha has an additional laplacian term which depends on ppMag.

If you include this term as done in the code, the flux must be subtracted of the contribution to the flux of the particle pressure computed at the previous time ste
p. If you do not do this, you account for the same term twice (once using the old value, the second time due to the laplacian).

Note that this becomes immediately clear if you perform the derivation: take the volume fraction equation, replace alpha*Ua with alpha_f*phia, as if you were deriving a pressure equation for the phase. Then extract the particle pressure term, and you will get exactly what is coded in OpenFOAM.

[cheng1988sjtu]

Hi, Aberto

Thanks for your reply!

Yes, there is a ppMagf in laplacian term, but if you have a chance to see the value of this term, it's zero, always, even after the pEqn.H has been solved. So my argue is that this particle-particle stress does not appear in alphaEqn.

Which I think is correct, ironically, since if you look at the pEqn, the velocity in particle phase has been modified already if you set g0>0. and I think the basic idea of this particle-particle idea is that through modification of Ua, we hope to constrain the alpha < alphaMax.

Is that correct?

Zhen

Originally Posted by alberto

The code is correct: the idea behind the procedure is to directly include the effect of the particle pressure in the volume fraction equation. You probably noticed that the equation for alpha has an additional laplacian term which depends on ppMag.

If you include this term as done in the code, the flux must be subtracted of the contribution to the flux of the particle pressure computed at the previous time ste
p. If you do not do this, you account for the same term twice (once using the old value, the second time due to the laplacian).

Note that this becomes immediately clear if you perform the derivation: take the volume fraction equation, replace alpha*Ua with alpha_f*phia, as if you were deriving a pressure equation for the phase. Then extract the particle pressure term, and you will get exactly what is coded in OpenFOAM.

[alberto]

In theory, the presence of the particle pressure gradient in the momentum equation should enforce the constraint alpha <= alphaMax. This is however a well known source of numerical instability in multi-fluid solvers, since the particle pressure term has a strongly non-linear dependency on the phase fraction.

In order to address this difficulty, you can include the effect of the particle pressure directly in the equation for the phase fraction. To do this, you proceed as I explained above, and you obtain exactly the form of the equation implemented in OpenFOAM.

There is no reason for ppMagf or for the Laplacian term containing it to be zero everywhere if the phase fraction is high enough. Note that in the case of ppMagf, the term is zero for alpha < alphaMax, so the particle maximum packing will be higher than the specified one. However, the approach is general, and can be extended to the kinetic theory formulation. I have recently done this and other changes to the solver, and published a manuscript on the topic on Powder Technology, if you are interested.

Best,

[cheng1988sjtu]

Hi Alberto,

Thanks for your patient answer!

I agree that there must a numerical trick in alphaEqn. however, let's derive step by step, and focus on where there is g0>0.

first, in pEqn.H, there is a modification in phiDraga, which means that the particle-particle stress is implemented directly in the momentum equation.[CODE]
if (g0.value() > 0.0)
{
phiDraga -= ppMagf*fvc::snGrad(alpha)*mesh.magSf();
}

/CODE]

after solving the pEqn.H, let's take a look at alphaEqn.H, the code modify the phir and phic again, from which we can see that it take into particle-particle stress on Ua twice ( once in pEqn, second here), in the code :
[CODE]
if (g0.value() > 0.0)
{
surfaceScalarField alphaf = fvc::interpolate(alpha);
surfaceScalarField phipp = ppMagf*fvc::snGrad(alpha)*mesh.magSf();
phir += phipp;
phic += fvc::interpolate(alpha)*phipp;
}
/CODE]

then, in the alphaEqn, it eliminate the particle-particle stress once, which is in the laplacian term:

        if (g0.value() > 0.0)
        {
            ppMagf = rUaAf*fvc::interpolate
            (
                (1.0/(rhoa*(alpha + scalar(0.00001))))

               *g0*min(exp(preAlphaExp*(alpha - alphaMax)), expMax)
            );

            alphaEqn -= fvm::laplacian
            (
                (fvc::interpolate(alpha) + scalar(0.00001))*ppMagf,
                alpha,
                "laplacian(alphaPpMag,alpha)"
            );
        }

Until now, the whole solving system has consider particle-particle stress just once ( one in pEqn, one in phia ( in phir and phic), and eliminate one in alphaEqn), and I agree with you that the code is correct, Thanks for pointing that out!

I've added a few lines in the code to print out the value of rUaAf in alphaEqn.H, and found that the value is uniformly 0, so I just want to clarify that, the original idea is to improve numerical stability, however, we found that the terms concerning ppMagf are zero indeed (I'm not to say that they should be zero, but it is zero as the way it was coded).

Can you quickly add the line in the alphaEqn.H and print out the result of rUaAf? if I'm wrong, could you please let me know? Thank you!

Many thanks!

Zhen

Originally Posted by alberto

In theory, the presence of the particle pressure gradient in the momentum equation should enforce the constraint alpha <= alphaMax. This is however a well known source of numerical instability in multi-fluid solvers, since the particle pressure term has a strongly non-linear dependency on the phase fraction.

In order to address this difficulty, you can include the effect of the particle pressure directly in the equation for the phase fraction. To do this, you proceed as I explained above, and you obtain exactly the form of the equation implemented in OpenFOAM.

There is no reason for ppMagf or for the Laplacian term containing it to be zero everywhere if the phase fraction is high enough. Note that in the case of ppMagf, the term is zero for alpha < alphaMax, so the particle maximum packing will be higher than the specified one. However, the approach is general, and can be extended to the kinetic theory formulation. I have recently done this and other changes to the solver, and published a manuscript on the topic on Powder Technology, if you are interested.

Best,

[alberto]

Hi,

rUaAf is zero only at the first time-step, since rUaA is initialized to zero. If you look at pEqn.H, you have

rUaAf = fvc::interpolate(rUaA);

I added the Info line in

        if (g0.value() > 0.0)
        {
        Info << "min(rUaAf) = " << min(rUaAf).value() << "max(rUaAf) = " <<
          max(rUaAf).value() << endl;
            ppMagf = rUaAf*fvc::interpolate
            (
                (1.0/(rhoa*(alpha + scalar(0.0001))))
               *g0*min(exp(preAlphaExp*(alpha - alphaMax)), expMax)
            );

            alphaEqn -= fvm::laplacian
            (
                (fvc::interpolate(alpha) + scalar(0.0001))*ppMagf,
                alpha,
                "laplacian(alphaPpMag,alpha)"
            );
        }

and the printout of the solver is:

Courant Number mean: 0.0450759 max: 0.0451754
Max Ur Courant Number = 0.0821135
deltaT = 0.00119048
Time = 0.10119

min(rUaAf) = 0max(rUaAf) = 0
DILUPBiCG:  Solving for alpha, Initial residual = 0.00552778, Final residual = 1.6595e-15, No Iterations 2
Dispersed phase volume fraction = 0.33099  Min(alpha) = 0  Max(alpha) = 0.410396
min(rUaAf) = 0max(rUaAf) = 0
DILUPBiCG:  Solving for alpha, Initial residual = 0.000331261, Final residual = 1.05713e-16, No Iterations 2
Dispersed phase volume fraction = 0.33099  Min(alpha) = 0  Max(alpha) = 0.410339
GAMG:  Solving for p, Initial residual = 0.000659113, Final residual = 2.65448e-05, No Iterations 1
time step continuity errors : sum local = 2.53952e-06, global = 3.11192e-07, cumulative = 3.11192e-07
GAMG:  Solving for p, Initial residual = 3.16091e-05, Final residual = 5.42006e-09, No Iterations 24
time step continuity errors : sum local = 5.93454e-07, global = 5.92783e-07, cumulative = 9.03975e-07
DILUPBiCG:  Solving for epsilon, Initial residual = 0.0249787, Final residual = 6.68909e-08, No Iterations 2
DILUPBiCG:  Solving for k, Initial residual = 0.0278856, Final residual = 1.44422e-07, No Iterations 2
ExecutionTime = 0.07 s  ClockTime = 0 s

Calculating averages

Courant Number mean: 0.0536524 max: 0.0537346
Max Ur Courant Number = 0.0999525
deltaT = 0.00141156
Time = 0.102602

min(rUaAf) = 0.00117375max(rUaAf) = 0.00118477
DILUPBiCG:  Solving for alpha, Initial residual = 0.00787966, Final residual = 1.03695e-14, No Iterations 2
Dispersed phase volume fraction = 0.331013  Min(alpha) = 0  Max(alpha) = 0.412684
min(rUaAf) = 0.00117375max(rUaAf) = 0.00118477
DILUPBiCG:  Solving for alpha, Initial residual = 0.000974836, Final residual = 5.82784e-16, No Iterations 2
Dispersed phase volume fraction = 0.331013  Min(alpha) = 0  Max(alpha) = 0.412501
GAMG:  Solving for p, Initial residual = 0.000288249, Final residual = 4.12155e-06, No Iterations 1
time step continuity errors : sum local = 1.20351e-06, global = 5.16087e-07, cumulative = 1.42006e-06
GAMG:  Solving for p, Initial residual = 1.48219e-05, Final residual = 5.379e-09, No Iterations 23
time step continuity errors : sum local = 6.99062e-07, global = 6.9825e-07, cumulative = 2.11831e-06
DILUPBiCG:  Solving for epsilon, Initial residual = 0.0263934, Final residual = 1.0387e-07, No Iterations 2
DILUPBiCG:  Solving for k, Initial residual = 0.0284996, Final residual = 2.31919e-07, No Iterations 2
ExecutionTime = 0.11 s  ClockTime = 0 s

Calculating averages

Courant Number mean: 0.0636171 max: 0.063778
Max Ur Courant Number = 0.118324
deltaT = 0.00167928
Time = 0.104281

min(rUaAf) = 0.00138437max(rUaAf) = 0.00140353
DILUPBiCG:  Solving for alpha, Initial residual = 0.00809664, Final residual = 1.78739e-14, No Iterations 2
Dispersed phase volume fraction = 0.33104  Min(alpha) = 0  Max(alpha) = 0.414937
min(rUaAf) = 0.00138437max(rUaAf) = 0.00140353
DILUPBiCG:  Solving for alpha, Initial residual = 0.00102087, Final residual = 1.24811e-15, No Iterations 2
Dispersed phase volume fraction = 0.33104  Min(alpha) = 0  Max(alpha) = 0.414745
GAMG:  Solving for p, Initial residual = 0.00019334, Final residual = 6.26108e-06, No Iterations 1
time step continuity errors : sum local = 1.57192e-06, global = 5.63326e-07, cumulative = 2.68164e-06
GAMG:  Solving for p, Initial residual = 1.70851e-05, Final residual = 6.02255e-09, No Iterations 23
time step continuity errors : sum local = 8.34896e-07, global = 8.33788e-07, cumulative = 3.51543e-06
DILUPBiCG:  Solving for epsilon, Initial residual = 0.0304074, Final residual = 2.14845e-07, No Iterations 2
DILUPBiCG:  Solving for k, Initial residual = 0.0331744, Final residual = 5.28396e-07, No Iterations 2
ExecutionTime = 0.16 s  ClockTime = 0 s

Calculating averages

Courant Number mean: 0.0756831 max: 0.0758539
Max Ur Courant Number = 0.140544
deltaT = 0.00199414
Time = 0.106275

min(rUaAf) = 0.00164154max(rUaAf) = 0.00166791
DILUPBiCG:  Solving for alpha, Initial residual = 0.00917989, Final residual = 2.60175e-14, No Iterations 2
Dispersed phase volume fraction = 0.331072  Min(alpha) = 0  Max(alpha) = 0.417257
min(rUaAf) = 0.00164154max(rUaAf) = 0.00166791
DILUPBiCG:  Solving for alpha, Initial residual = 0.001015, Final residual = 2.1767e-15, No Iterations 2
Dispersed phase volume fraction = 0.331072  Min(alpha) = 0  Max(alpha) = 0.417063
GAMG:  Solving for p, Initial residual = 0.000238078, Final residual = 9.40541e-06, No Iterations 1
time step continuity errors : sum local = 2.19634e-06, global = 5.81552e-07, cumulative = 4.09698e-06
GAMG:  Solving for p, Initial residual = 2.04994e-05, Final residual = 7.48965e-09, No Iterations 23
time step continuity errors : sum local = 9.87082e-07, global = 9.85411e-07, cumulative = 5.08239e-06
DILUPBiCG:  Solving for epsilon, Initial residual = 0.0347589, Final residual = 4.36944e-07, No Iterations 2
DILUPBiCG:  Solving for k, Initial residual = 0.0383839, Final residual = 1.04614e-06, No Iterations 2
ExecutionTime = 0.2 s  ClockTime = 0 s

which shows rUaAf is not uniformly zero.

Best,

[cheng1988sjtu]

Hi, Alberto,

Great!, thanks for your information, there must be something wrong with my coding.

Thanks again, now my doubt has been removed, cheers!

Best

Zhen

Originally Posted by alberto

Hi,

rUaAf is zero only at the first time-step, since rUaA is initialized to zero. If you look at pEqn.H, you have

rUaAf = fvc::interpolate(rUaA);

I added the Info line in

        if (g0.value() > 0.0)
        {
        Info << "min(rUaAf) = " << min(rUaAf).value() << "max(rUaAf) = " <<
          max(rUaAf).value() << endl;
            ppMagf = rUaAf*fvc::interpolate
            (
                (1.0/(rhoa*(alpha + scalar(0.0001))))
               *g0*min(exp(preAlphaExp*(alpha - alphaMax)), expMax)
            );

            alphaEqn -= fvm::laplacian
            (
                (fvc::interpolate(alpha) + scalar(0.0001))*ppMagf,
                alpha,
                "laplacian(alphaPpMag,alpha)"
            );
        }

and the printout of the solver is:

Courant Number mean: 0.0450759 max: 0.0451754
Max Ur Courant Number = 0.0821135
deltaT = 0.00119048
Time = 0.10119

min(rUaAf) = 0max(rUaAf) = 0
DILUPBiCG:  Solving for alpha, Initial residual = 0.00552778, Final residual = 1.6595e-15, No Iterations 2
Dispersed phase volume fraction = 0.33099  Min(alpha) = 0  Max(alpha) = 0.410396
min(rUaAf) = 0max(rUaAf) = 0
DILUPBiCG:  Solving for alpha, Initial residual = 0.000331261, Final residual = 1.05713e-16, No Iterations 2
Dispersed phase volume fraction = 0.33099  Min(alpha) = 0  Max(alpha) = 0.410339
GAMG:  Solving for p, Initial residual = 0.000659113, Final residual = 2.65448e-05, No Iterations 1
time step continuity errors : sum local = 2.53952e-06, global = 3.11192e-07, cumulative = 3.11192e-07
GAMG:  Solving for p, Initial residual = 3.16091e-05, Final residual = 5.42006e-09, No Iterations 24
time step continuity errors : sum local = 5.93454e-07, global = 5.92783e-07, cumulative = 9.03975e-07
DILUPBiCG:  Solving for epsilon, Initial residual = 0.0249787, Final residual = 6.68909e-08, No Iterations 2
DILUPBiCG:  Solving for k, Initial residual = 0.0278856, Final residual = 1.44422e-07, No Iterations 2
ExecutionTime = 0.07 s  ClockTime = 0 s

Calculating averages

Courant Number mean: 0.0536524 max: 0.0537346
Max Ur Courant Number = 0.0999525
deltaT = 0.00141156
Time = 0.102602

min(rUaAf) = 0.00117375max(rUaAf) = 0.00118477
DILUPBiCG:  Solving for alpha, Initial residual = 0.00787966, Final residual = 1.03695e-14, No Iterations 2
Dispersed phase volume fraction = 0.331013  Min(alpha) = 0  Max(alpha) = 0.412684
min(rUaAf) = 0.00117375max(rUaAf) = 0.00118477
DILUPBiCG:  Solving for alpha, Initial residual = 0.000974836, Final residual = 5.82784e-16, No Iterations 2
Dispersed phase volume fraction = 0.331013  Min(alpha) = 0  Max(alpha) = 0.412501
GAMG:  Solving for p, Initial residual = 0.000288249, Final residual = 4.12155e-06, No Iterations 1
time step continuity errors : sum local = 1.20351e-06, global = 5.16087e-07, cumulative = 1.42006e-06
GAMG:  Solving for p, Initial residual = 1.48219e-05, Final residual = 5.379e-09, No Iterations 23
time step continuity errors : sum local = 6.99062e-07, global = 6.9825e-07, cumulative = 2.11831e-06
DILUPBiCG:  Solving for epsilon, Initial residual = 0.0263934, Final residual = 1.0387e-07, No Iterations 2
DILUPBiCG:  Solving for k, Initial residual = 0.0284996, Final residual = 2.31919e-07, No Iterations 2
ExecutionTime = 0.11 s  ClockTime = 0 s

Calculating averages

Courant Number mean: 0.0636171 max: 0.063778
Max Ur Courant Number = 0.118324
deltaT = 0.00167928
Time = 0.104281

min(rUaAf) = 0.00138437max(rUaAf) = 0.00140353
DILUPBiCG:  Solving for alpha, Initial residual = 0.00809664, Final residual = 1.78739e-14, No Iterations 2
Dispersed phase volume fraction = 0.33104  Min(alpha) = 0  Max(alpha) = 0.414937
min(rUaAf) = 0.00138437max(rUaAf) = 0.00140353
DILUPBiCG:  Solving for alpha, Initial residual = 0.00102087, Final residual = 1.24811e-15, No Iterations 2
Dispersed phase volume fraction = 0.33104  Min(alpha) = 0  Max(alpha) = 0.414745
GAMG:  Solving for p, Initial residual = 0.00019334, Final residual = 6.26108e-06, No Iterations 1
time step continuity errors : sum local = 1.57192e-06, global = 5.63326e-07, cumulative = 2.68164e-06
GAMG:  Solving for p, Initial residual = 1.70851e-05, Final residual = 6.02255e-09, No Iterations 23
time step continuity errors : sum local = 8.34896e-07, global = 8.33788e-07, cumulative = 3.51543e-06
DILUPBiCG:  Solving for epsilon, Initial residual = 0.0304074, Final residual = 2.14845e-07, No Iterations 2
DILUPBiCG:  Solving for k, Initial residual = 0.0331744, Final residual = 5.28396e-07, No Iterations 2
ExecutionTime = 0.16 s  ClockTime = 0 s

Calculating averages

Courant Number mean: 0.0756831 max: 0.0758539
Max Ur Courant Number = 0.140544
deltaT = 0.00199414
Time = 0.106275

min(rUaAf) = 0.00164154max(rUaAf) = 0.00166791
DILUPBiCG:  Solving for alpha, Initial residual = 0.00917989, Final residual = 2.60175e-14, No Iterations 2
Dispersed phase volume fraction = 0.331072  Min(alpha) = 0  Max(alpha) = 0.417257
min(rUaAf) = 0.00164154max(rUaAf) = 0.00166791
DILUPBiCG:  Solving for alpha, Initial residual = 0.001015, Final residual = 2.1767e-15, No Iterations 2
Dispersed phase volume fraction = 0.331072  Min(alpha) = 0  Max(alpha) = 0.417063
GAMG:  Solving for p, Initial residual = 0.000238078, Final residual = 9.40541e-06, No Iterations 1
time step continuity errors : sum local = 2.19634e-06, global = 5.81552e-07, cumulative = 4.09698e-06
GAMG:  Solving for p, Initial residual = 2.04994e-05, Final residual = 7.48965e-09, No Iterations 23
time step continuity errors : sum local = 9.87082e-07, global = 9.85411e-07, cumulative = 5.08239e-06
DILUPBiCG:  Solving for epsilon, Initial residual = 0.0347589, Final residual = 4.36944e-07, No Iterations 2
DILUPBiCG:  Solving for k, Initial residual = 0.0383839, Final residual = 1.04614e-06, No Iterations 2
ExecutionTime = 0.2 s  ClockTime = 0 s

which shows rUaAf is not uniformly zero.

Best,

[alberto]

You are welcome :-)

[cheng1988sjtu]

Hi Alberto,

I'm sorry to bother again.

actually, I did as you code in the alphaEqn, and compile using Allwmake, and the result is still like this:

[HTML]Starting time loop

Courant Number mean: 0.0110502 max: 0.05
Max Ur Courant Number = 0.05
Reading/calculating field UaMean

Reading/calculating field UbMean

Reading/calculating field alphaMean

Reading/calculating field pMean

fieldAverage: starting averaging at time 0

Time = 0.001

max(rUaAf) = 0 min(rUaAf) = 0
max(ppMagf) = 0 min(ppMagf) = 0
DILUPBiCG: Solving for alpha, Initial residual = 3.90098e-17, Final residual = 3.90098e-17, No Iterations 0
Dispersed phase volume fraction = 0.2015 Min(alpha) = 0 Max(alpha) = 0.62
max(rUaAf) = 0 min(rUaAf) = 0
max(ppMagf) = 0 min(ppMagf) = 0
DILUPBiCG: Solving for alpha, Initial residual = 3.90098e-17, Final residual = 3.90098e-17, No Iterations 0
Dispersed phase volume fraction = 0.2015 Min(alpha) = 0 Max(alpha) = 0.62
pressure gradient in UEqn= (0.081 0 0)
kinTheory: max(Theta) = 1e-15
kinTheory: min(nua) = 5.06848e-07, max(nua) = 0.04
kinTheory: min(pa) = 0, max(pa) = 2.30439e-06
GAMG: Solving for p, Initial residual = 1, Final residual = 9.52665e-09, No Iterations 23
time step continuity errors : sum local = 2.15521e-12, global = -2.15521e-12, cumulative = -2.15521e-12
max(rUaAf) = 0 min(rUaAf) = 0
max(ppMagf) = 0 min(ppMagf) = 0
DILUPBiCG: Solving for alpha, Initial residual = 0.00104655, Final residual = 1.46481e-11, No Iterations 3
Dispersed phase volume fraction = 0.2015 Min(alpha) = 0 Max(alpha) = 0.620052
max(rUaAf) = 0 min(rUaAf) = 0
max(ppMagf) = 0 min(ppMagf) = 0
DILUPBiCG: Solving for alpha, Initial residual = 4.06684e-05, Final residual = 8.19921e-12, No Iterations 3
Dispersed phase volume fraction = 0.2015 Min(alpha) = -1.64854e-21 Max(alpha) = 0.620052
GAMG: Solving for p, Initial residual = 0.0786421, Final residual = 9.35372e-09, No Iterations 17
time step continuity errors : sum local = 2.83579e-11, global = 2.83579e-11, cumulative = 2.62027e-11
DILUPBiCG: Solving for epsilon, Initial residual = 0.135913, Final residual = 2.09582e-06, No Iterations 4
DILUPBiCG: Solving for k, Initial residual = 1, Final residual = 1.43273e-06, No Iterations 5
pressure gradient = (0.081 0 0)
ExecutionTime = 0.63 s ClockTime = 4 s

Courant Number mean: 0.010322 max: 0.0500001
Max Ur Courant Number = 0.0343635
Calculating averages

Time = 0.002

max(rUaAf) = 0 min(rUaAf) = 0
max(ppMagf) = 0 min(ppMagf) = 0
DILUPBiCG: Solving for alpha, Initial residual = 0.0010924, Final residual = 8.48925e-11, No Iterations 3
Dispersed phase volume fraction = 0.2015 Min(alpha) = -8.33627e-28 Max(alpha) = 0.620094
max(rUaAf) = 0 min(rUaAf) = 0
max(ppMagf) = 0 min(ppMagf) = 0
DILUPBiCG: Solving for alpha, Initial residual = 4.58883e-05, Final residual = 5.27276e-12, No Iterations 3
Dispersed phase volume fraction = 0.2015 Min(alpha) = -1.55414e-157 Max(alpha) = 0.620094
pressure gradient in UEqn= (0.081 0 0)
kinTheory: max(Theta) = 58.918
kinTheory: min(nua) = 1.15627e-12, max(nua) = 0.04
kinTheory: min(pa) = -3.88536e-169, max(pa) = 2209.46
GAMG: Solving for p, Initial residual = 0.100408, Final residual = 8.11056e-09, No Iterations 15
time step continuity errors : sum local = 2.49503e-11, global = -2.49503e-11, cumulative = 1.25243e-12
max(rUaAf) = 0 min(rUaAf) = 0
max(ppMagf) = 0 min(ppMagf) = 0
DILUPBiCG: Solving for alpha, Initial residual = 0.000897803, Final residual = 1.59916e-12, No Iterations 4
Dispersed phase volume fraction = 0.2015 Min(alpha) = -1.45549e-14 Max(alpha) = 0.62193
max(rUaAf) = 0 min(rUaAf) = 0
max(ppMagf) = 0 min(ppMagf) = 0
DILUPBiCG: Solving for alpha, Initial residual = 0.000133328, Final residual = 1.37526e-11, No Iterations 3
Dispersed phase volume fraction = 0.2015 Min(alpha) = -1.26941e-26 Max(alpha) = 0.621929
GAMG: Solving for p, Initial residual = 0.034853, Final residual = 5.40877e-09, No Iterations 17
time step continuity errors : sum local = 2.02786e-11, global = 2.02786e-11, cumulative = 2.15311e-11
DILUPBiCG: Solving for epsilon, Initial residual = 0.0187545, Final residual = 8.89096e-07, No Iterations 3
DILUPBiCG: Solving for k, Initial residual = 0.800558, Final residual = 3.72508e-06, No Iterations 4
pressure gradient = (0.081 0 0)
ExecutionTime = 0.93 s ClockTime = 5 s

Courant Number mean: 0.010183 max: 0.0500034
Max Ur Courant Number = 0.0826851
Calculating averages

Time = 0.003

max(rUaAf) = 0 min(rUaAf) = 0
max(ppMagf) = 0 min(ppMagf) = 0
DILUPBiCG: Solving for alpha, Initial residual = 0.00188765, Final residual = 1.59677e-12, No Iterations 4
Dispersed phase volume fraction = 0.2015 Min(alpha) = -6.81511e-22 Max(alpha) = 0.62233
max(rUaAf) = 0 min(rUaAf) = 0
max(ppMagf) = 0 min(ppMagf) = 0
DILUPBiCG: Solving for alpha, Initial residual = 0.000223479, Final residual = 1.52428e-11, No Iterations 3
Dispersed phase volume fraction = 0.2015 Min(alpha) = -7.91657e-26 Max(alpha) = 0.622328
pressure gradient in UEqn= (0.081 0 0)
kinTheory: max(Theta) = 1000
kinTheory: min(nua) = 1.60279e-12, max(nua) = 0.04
kinTheory: min(pa) = -3.08412e-37, max(pa) = 2045.12
GAMG: Solving for p, Initial residual = 0.122225, Final residual = 7.73269e-09, No Iterations 17
time step continuity errors : sum local = 2.63962e-11, global = -2.63962e-11, cumulative = -4.8651e-12
max(rUaAf) = 0 min(rUaAf) = 0
max(ppMagf) = 0 min(ppMagf) = 0
DILUPBiCG: Solving for alpha, Initial residual = 0.000356051, Final residual = 2.09863e-11, No Iterations 3
Dispersed phase volume fraction = 0.2015 Min(alpha) = -4.54596e-16 Max(alpha) = 0.621036
max(rUaAf) = 0 min(rUaAf) = 0
max(ppMagf) = 0 min(ppMagf) = 0
DILUPBiCG: Solving for alpha, Initial residual = 8.30811e-05, Final residual = 7.85124e-13, No Iterations 3
Dispersed phase volume fraction = 0.2015 Min(alpha) = -1.08298e-25 Max(alpha) = 0.621036
GAMG: Solving for p, Initial residual = 0.016629, Final residual = 6.94633e-09, No Iterations 14
time step continuity errors : sum local = 2.18021e-11, global = 2.18021e-11, cumulative = 1.6937e-11
DILUPBiCG: Solving for epsilon, Initial residual = 0.020499, Final residual = 8.62327e-07, No Iterations 3
DILUPBiCG: Solving for k, Initial residual = 0.471488, Final residual = 2.08913e-06, No Iterations 4
pressure gradient = (0.081 0 0)
ExecutionTime = 1.23 s ClockTime = 6 s

Courant Number mean: 0.0101149 max: 0.0500028
Max Ur Courant Number = 0.126741
Calculating averages

Time = 0.004

max(rUaAf) = 0 min(rUaAf) = 0
max(ppMagf) = 0 min(ppMagf) = 0
DILUPBiCG: Solving for alpha, Initial residual = 0.00167317, Final residual = 6.278e-11, No Iterations 3
Dispersed phase volume fraction = 0.2015 Min(alpha) = -6.07852e-25 Max(alpha) = 0.62135
max(rUaAf) = 0 min(rUaAf) = 0
max(ppMagf) = 0 min(ppMagf) = 0
DILUPBiCG: Solving for alpha, Initial residual = 0.000158778, Final residual = 3.213e-12, No Iterations 3
Dispersed phase volume fraction = 0.2015 Min(alpha) = -1.04731e-25 Max(alpha) = 0.621349
pressure gradient in UEqn= (0.081 0 0)
kinTheory: max(Theta) = 1000
kinTheory: min(nua) = 1.60279e-12, max(nua) = 0.04
kinTheory: min(pa) = -7.28917e-21, max(pa) = 1950.06
GAMG: Solving for p, Initial residual = 0.0905112, Final residual = 9.85605e-09, No Iterations 15
time step continuity errors : sum local = 2.99527e-11, global = 2.99527e-11, cumulative = 4.68897e-11
max(rUaAf) = 0 min(rUaAf) = 0
max(ppMagf) = 0 min(ppMagf) = 0
DILUPBiCG: Solving for alpha, Initial residual = 0.000346857, Final residual = 3.04659e-12, No Iterations 3
Dispersed phase volume fraction = 0.2015 Min(alpha) = -3.26225e-16 Max(alpha) = 0.621475
max(rUaAf) = 0 min(rUaAf) = 0
max(ppMagf) = 0 min(ppMagf) = 0
DILUPBiCG: Solving for alpha, Initial residual = 7.42358e-05, Final residual = 1.26536e-12, No Iterations 3
Dispersed phase volume fraction = 0.2015 Min(alpha) = -2.41699e-25 Max(alpha) = 0.621475
GAMG: Solving for p, Initial residual = 0.0111199, Final residual = 9.2887e-09, No Iterations 12
time step continuity errors : sum local = 2.6559e-11, global = 2.6559e-11, cumulative = 7.34486e-11
DILUPBiCG: Solving for epsilon, Initial residual = 0.0216537, Final residual = 8.88015e-06, No Iterations 2
DILUPBiCG: Solving for k, Initial residual = 0.272703, Final residual = 7.4626e-06, No Iterations 3
pressure gradient = (0.081 0 0)
ExecutionTime = 1.51 s ClockTime = 6 s

Courant Number mean: 0.010065 max: 0.050003
Max Ur Courant Number = 0.145256
Calculating averages

Time = 0.005

max(rUaAf) = 0 min(rUaAf) = 0
max(ppMagf) = 0 min(ppMagf) = 0
DILUPBiCG: Solving for alpha, Initial residual = 0.0014732, Final residual = 3.10609e-13, No Iterations 4
Dispersed phase volume fraction = 0.2015 Min(alpha) = -4.13607e-23 Max(alpha) = 0.621453
max(rUaAf) = 0 min(rUaAf) = 0
max(ppMagf) = 0 min(ppMagf) = 0
DILUPBiCG: Solving for alpha, Initial residual = 0.00013011, Final residual = 3.92767e-12, No Iterations 3
Dispersed phase volume fraction = 0.2015 Min(alpha) = -2.2714e-25 Max(alpha) = 0.621452
pressure gradient in UEqn= (0.081 0 0)
kinTheory: max(Theta) = 1000
kinTheory: min(nua) = 1.60279e-12, max(nua) = 0.04
kinTheory: min(pa) = -1.4887e-21, max(pa) = 1860.09
GAMG: Solving for p, Initial residual = 0.0898181, Final residual = 9.60529e-09, No Iterations 16
time step continuity errors : sum local = 2.78988e-11, global = 2.78988e-11, cumulative = 1.01347e-10
max(rUaAf) = 0 min(rUaAf) = 0
max(ppMagf) = 0 min(ppMagf) = 0
DILUPBiCG: Solving for alpha, Initial residual = 0.000473163, Final residual = 2.11461e-11, No Iterations 3
Dispersed phase volume fraction = 0.2015 Min(alpha) = -3.87245e-18 Max(alpha) = 0.621414
max(rUaAf) = 0 min(rUaAf) = 0
max(ppMagf) = 0 min(ppMagf) = 0
DILUPBiCG: Solving for alpha, Initial residual = 6.93413e-05, Final residual = 3.96824e-12, No Iterations 3
Dispersed phase volume fraction = 0.2015 Min(alpha) = -3.91346e-26 Max(alpha) = 0.621414
GAMG: Solving for p, Initial residual = 0.0374737, Final residual = 9.59322e-09, No Iterations 13
time step continuity errors : sum local = 2.7584e-11, global = 2.75837e-11, cumulative = 1.28931e-10
DILUPBiCG: Solving for epsilon, Initial residual = 0.0222171, Final residual = 9.84455e-06, No Iterations 2
DILUPBiCG: Solving for k, Initial residual = 0.196096, Final residual = 4.79474e-06, No Iterations 3
pressure gradient = (0.081 0 0)
ExecutionTime = 1.79 s ClockTime = 7 s

Courant Number mean: 0.010023 max: 0.0500031
Max Ur Courant Number = 0.150859
Calculating averages

Time = 0.006

max(rUaAf) = 0 min(rUaAf) = 0
max(ppMagf) = 0 min(ppMagf) = 0
DILUPBiCG: Solving for alpha, Initial residual = 0.00166162, Final residual = 2.33643e-11, No Iterations 4
Dispersed phase volume fraction = 0.2015 Min(alpha) = -3.90035e-25 Max(alpha) = 0.621333
max(rUaAf) = 0 min(rUaAf) = 0
max(ppMagf) = 0 min(ppMagf) = 0
DILUPBiCG: Solving for alpha, Initial residual = 0.000128721, Final residual = 5.15395e-12, No Iterations 3
Dispersed phase volume fraction = 0.2015 Min(alpha) = -1.56767e-25 Max(alpha) = 0.621334
pressure gradient in UEqn= (0.081 0 0)
kinTheory: max(Theta) = 1000
kinTheory: min(nua) = 1.60279e-12, max(nua) = 0.04
kinTheory: min(pa) = -5.56552e-23, max(pa) = 1778.39
GAMG: Solving for p, Initial residual = 0.0744853, Final residual = 8.1907e-09, No Iterations 17
time step continuity errors : sum local = 2.52555e-11, global = 2.52555e-11, cumulative = 1.54187e-10
max(rUaAf) = 0 min(rUaAf) = 0
max(ppMagf) = 0 min(ppMagf) = 0
DILUPBiCG: Solving for alpha, Initial residual = 0.000528506, Final residual = 3.56838e-11, No Iterations 3
Dispersed phase volume fraction = 0.2015 Min(alpha) = -4.19315e-17 Max(alpha) = 0.621249
max(rUaAf) = 0 min(rUaAf) = 0
max(ppMagf) = 0 min(ppMagf) = 0
DILUPBiCG: Solving for alpha, Initial residual = 9.13973e-05, Final residual = 4.25874e-12, No Iterations 3
Dispersed phase volume fraction = 0.2015 Min(alpha) = -4.63227e-25 Max(alpha) = 0.621248
GAMG: Solving for p, Initial residual = 0.0460487, Final residual = 7.05426e-09, No Iterations 15
time step continuity errors : sum local = 2.20093e-11, global = 2.20093e-11, cumulative = 1.76196e-10
DILUPBiCG: Solving for epsilon, Initial residual = 0.0228629, Final residual = 9.47428e-06, No Iterations 2
DILUPBiCG: Solving for k, Initial residual = 0.155566, Final residual = 3.88173e-06, No Iterations 3
pressure gradient = (0.081 0 0)
ExecutionTime = 2.1 s ClockTime = 8 s

Courant Number mean: 0.00998571 max: 0.0500031
Max Ur Courant Number = 0.150683
Calculating averages

Time = 0.007

max(rUaAf) = 0 min(rUaAf) = 0
max(ppMagf) = 0 min(ppMagf) = 0
DILUPBiCG: Solving for alpha, Initial residual = 0.00188168, Final residual = 7.31724e-12, No Iterations 4
Dispersed phase volume fraction = 0.2015 Min(alpha) = -1.087e-24 Max(alpha) = 0.62124
max(rUaAf) = 0 min(rUaAf) = 0
max(ppMagf) = 0 min(ppMagf) = 0
DILUPBiCG: Solving for alpha, Initial residual = 0.000164972, Final residual = 7.34103e-12, No Iterations 3
Dispersed phase volume fraction = 0.2015 Min(alpha) = -2.60664e-24 Max(alpha) = 0.621241
pressure gradient in UEqn= (0.081 0 0)
kinTheory: max(Theta) = 1000
kinTheory: min(nua) = 1.60279e-12, max(nua) = 0.04
kinTheory: min(pa) = -2.89806e-19, max(pa) = 1716.94
GAMG: Solving for p, Initial residual = 0.0786497, Final residual = 8.33953e-09, No Iterations 16
time step continuity errors : sum local = 2.72705e-11, global = 2.72705e-11, cumulative = 2.03467e-10
max(rUaAf) = 0 min(rUaAf) = 0
max(ppMagf) = 0 min(ppMagf) = 0
DILUPBiCG: Solving for alpha, Initial residual = 0.000359687, Final residual = 2.1303e-11, No Iterations 3
Dispersed phase volume fraction = 0.2015 Min(alpha) = -7.1831e-18 Max(alpha) = 0.621369
max(rUaAf) = 0 min(rUaAf) = 0
max(ppMagf) = 0 min(ppMagf) = 0
DILUPBiCG: Solving for alpha, Initial residual = 6.24967e-05, Final residual = 2.79162e-12, No Iterations 3
Dispersed phase volume fraction = 0.2015 Min(alpha) = -6.99314e-28 Max(alpha) = 0.621369
GAMG: Solving for p, Initial residual = 0.0277413, Final residual = 6.64361e-09, No Iterations 14
time step continuity errors : sum local = 2.22131e-11, global = 2.22131e-11, cumulative = 2.2568e-10
DILUPBiCG: Solving for epsilon, Initial residual = 0.0237308, Final residual = 5.02822e-07, No Iterations 3
DILUPBiCG: Solving for k, Initial residual = 0.130643, Final residual = 3.80861e-06, No Iterations 3
pressure gradient = (0.081 0 0)
ExecutionTime = 2.4 s ClockTime = 9 s

Courant Number mean: 0.00995146 max: 0.050003
Max Ur Courant Number = 0.147262
Calculating averages/HTML]

the code is like this:
[CODE] if (g0.value() > 0.0)
{
Info<<"max(rUaAf) = "<<max(rUaAf).value()<<" min(rUaAf) = "<<min(rUaAf).value()<<endl;
ppMagf = rUaAf*fvc::interpolate
(
(1.0/(rhoa*(alpha + scalar(0.0001))))
*g0*min(exp(preAlphaExp*(alpha - alphaMax)), expMax)
);
// Info<<"max(rUaAf) = "<<max(rUaAf).value()<<"min(rUaAf) = "<<min(rUaAf).value()<<endl;
Info<<"max(ppMagf) = "<<max(ppMagf).value()<<" min(ppMagf) = "<<min(ppMagf).value()<<endl;
alphaEqn -= fvm::laplacian
(
(fvc::interpolate(alpha) + scalar(0.0001))*ppMagf,
alpha,
"laplacian(alphaPpMag,alpha)"
);
}
/CODE]

so what I did wrong? do you have any idea?

Originally Posted by alberto

You are welcome :-)

[alberto]

I use OpenFOAM 2.0.x, and tested on the tutorial for twoPhaseEulerFoam called "bed". Maybe you want to test on the same tutorial.

Best,

[cheng1988sjtu]

Hi Alberto,

Yeah, I'm using 1.7.1, maybe that's the reason. I'll try to implement 2.0, and see what's the difference.

Thanks!

Zhen

Originally Posted by alberto

I use OpenFOAM 2.0.x, and tested on the tutorial for twoPhaseEulerFoam called "bed". Maybe you want to test on the same tutorial.

Best,

[cheng1988sjtu]

Hi,

I've installed the openFOAM-2.0.1, and had a look at the twoPhaseEulerFoam, it turns out that they are using p.storePrevIter() to store the information obtained in pEqn, which has solved my concern. And solve the problem for OpenFOAM-1.7.1.

Best

Zhen

[alberto]

Originally Posted by cheng1988sjtu

Hi,

I've installed the openFOAM-2.0.1, and had a look at the twoPhaseEulerFoam, it turns out that they are using p.storePrevIter() to store the information obtained in pEqn, which has solved my concern. And solve the problem for OpenFOAM-1.7.1.

Best

Zhen

The presence of p.storePrevIter() is only there because in OpenFOAM 2.0.x the twoPhaseEulerFoam code uses the PIMPLE algorithm, which allows relaxation of the pressure.

This has nothing to do with your former observation about rUaAf being zero. The twoPhaseEulerFoam in both 1.7.x and 2.0.x does not have rUaAf uniformly zero, since this would mean the reciprocal of the central coefficient rUaA (= 1/A) would be zero everywhere too, which is clearly not the case.

Did you modify the code you were using in 1.7.x?

Best,
Alberto

[cheng1988sjtu]

Hi,

Yes, I've modified the code in 1.7.1 a little bit, but it's just in UEqns.H and twoPhaseEulerFoam.C, and I've tried to print out the rUaAf in pEqn.H. it's not zero, but in alphaEqn.H, it's zero, I've no idea what's going on.

Did you try to see the value of rUaAf in alphaEqn.H in 1.7.1?

Best

Zhen

Originally Posted by alberto

The presence of p.storePrevIter() is only there because in OpenFOAM 2.0.x the twoPhaseEulerFoam code uses the PIMPLE algorithm, which allows relaxation of the pressure.

This has nothing to do with your former observation about rUaAf being zero. The twoPhaseEulerFoam in both 1.7.x and 2.0.x does not have rUaAf uniformly zero, since this would mean the reciprocal of the central coefficient rUaA (= 1/A) would be zero everywhere too, which is clearly not the case.

Did you modify the code you were using in 1.7.x?

Best,
Alberto

[alberto]

I switched some month ago to 2.0.x, however I did all my implementation in 1.7.x, and I do not meet the problem. Anyways, I would suggest to use 2.0.x, which has many improvements :-)

[openfoammaofnepo]

Dear Charlie,

For twoPhaseEulerFoam, did you try any LES for particle laden flows, i.e. gas is the continuous phase and particle is the dispersed phase? Do you have any comments about running LES with twoPhaseEulerFOAM? Thanks.

Originally Posted by cheng1988sjtu

Hi,

I've installed the openFOAM-2.0.1, and had a look at the twoPhaseEulerFoam, it turns out that they are using p.storePrevIter() to store the information obtained in pEqn, which has solved my concern. And solve the problem for OpenFOAM-1.7.1.

Best

Zhen