CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM (https://www.cfd-online.com/Forums/openfoam/)
-   -   twoPhaseEulerFoam:pEqn.H:"ppDrag" for KTGF confusing: why "alpha1f" omitted (https://www.cfd-online.com/Forums/openfoam/125035-twophaseeulerfoam-peqn-h-ppdrag-ktgf-confusing-why-alpha1f-omitted.html)

ysun October 17, 2013 13:08

twoPhaseEulerFoam:pEqn.H:"ppDrag" for KTGF confusing: why "alpha1f" omitted
 
Dear all,

I've been trying to understand what "twoPhaseEulerFoam" does in detail, but there is one problem that keeps puzzling me. I'm using OpenFOAM-2.2.x, but this problem is the same in OpenFOAM-1.6-ext. The "ppDrag" term accounts for the particle pressure gradient in the momentum equation: -\frac{\nabla p_s}{\rho_s \alpha_s} (van Wachem's thesis, Table 3.1, p. 43). Note that there is the solid volume fraction term in the denominator, which is different from the gas pressure gradient. The powder modulus version of "ppDrag" takes this into account:
Code:

    if (g0.value() > 0.0)
    {
        ppDrag -= ppMagf*fvc::snGrad(alpha1)*mesh.magSf();
    }

where the division of \alpha_s was carried in "ppMagf" defined in "alphaEqn.H"
Code:

            ppMagf =
                rAU1f/(alpha1f + scalar(0.0001))
              *(g0/rho1)*min(exp(preAlphaExp*(alpha1f - alphaMax)), expMax);

However, the KTGF version of "ppDrag" does not include the solid volume fraction in the denominator:
Code:

    if (kineticTheory.on())
    {
        ppDrag -= rAU1f*fvc::snGrad(kineticTheory.pa()/rho1)*mesh.magSf();
    }

I cannot figure out what the justification is for such a different treatment. I tried to add \alpha_s in the denominator and recompile the solver, but it just failed when running on the "bed2" case in the tutorial. The change I made to the above code is:
Code:

    if (kineticTheory.on())
    {
        ppDrag -= rAU1f/(alpha1f + scalar(0.0001))*fvc::snGrad(kineticTheory.pa()/rho1)*mesh.magSf();
    }

When running on "bed2", "max(pa)" suddenly skyrocketed after a few milliseconds and the solver stopped running. I will attach the relevant piece of the result at the end of this post, but before that I think it's worth mentioning that Prof. Passalacqua's "twoPhaseEulerPimpleFoam" also includes \alpha_s in the denominator, although indirectly:
Code:

    if (g0.value() > 0.0 || kineticTheory.on())
    {
        phiDraga -= ppMagf*fvc::snGrad(alpha)*mesh.magSf();
    }

where "ppMagf" is defined for KTGF too:
Code:

        ppMagf = rUaAf*fvc::interpolate
        (
        (1.0/(rhoa*(alpha + scalar(0.0001))))
          *kineticTheory.ppMagf(alpha)
        );

I would greatly appreciate any insight regarding why \alpha_s is omitted in "twoPhaseEulerFoam" and why it doesn't work if added it back. Thanks!

Yujian



Quote:

Time = 0.007

MULES: Solving for alpha1
MULES: Solving for alpha1
MULES: Solving for alpha1
Dispersed phase volume fraction = 0.275 Min(alpha1) = 0 Max(alpha1) = 0.550043
kinTheory: max(Theta) = 0.00609332
kinTheory: min(nu1) = 9.61677e-13, max(nu1) = 2.0919e-05
kinTheory: min(pa) = 0, max(pa) = 363.058
GAMG: Solving for p, Initial residual = 5.87016e-05, Final residual = 5.03626e-06, No Iterations 1
time step continuity errors : sum local = 1.55852e-06, global = 8.16546e-08, cumulative = -8.87394e-06
GAMG: Solving for p, Initial residual = 4.29619e-05, Final residual = 7.37734e-09, No Iterations 11
time step continuity errors : sum local = 2.27498e-09, global = 2.21685e-09, cumulative = -8.87172e-06
DILUPBiCG: Solving for epsilon, Initial residual = 0.0293547, Final residual = 2.26688e-06, No Iterations 2
DILUPBiCG: Solving for k, Initial residual = 0.178011, Final residual = 1.12541e-06, No Iterations 3
ExecutionTime = 1.39 s ClockTime = 1 s

fieldAverage fieldAverage1 output:
Calculating averages

Courant Number mean: 0.0250311 max: 0.0353292
Max Ur Courant Number = 0.0775321
Time = 0.0075

MULES: Solving for alpha1
MULES: Solving for alpha1
MULES: Solving for alpha1
Dispersed phase volume fraction = 0.275 Min(alpha1) = 0 Max(alpha1) = 0.550056
kinTheory: max(Theta) = 31.692
kinTheory: min(nu1) = 9.61677e-13, max(nu1) = 0.000154818
kinTheory: min(pa) = 0, max(pa) = 1.83653e+06
GAMG: Solving for p, Initial residual = 0.00788312, Final residual = 0.000210995, No Iterations 1
time step continuity errors : sum local = 6.48878e-05, global = 1.84827e-10, cumulative = -8.87154e-06
GAMG: Solving for p, Initial residual = 0.0541924, Final residual = 6.9458e-09, No Iterations 18
time step continuity errors : sum local = 2.25176e-09, global = -2.19516e-09, cumulative = -8.87374e-06
DILUPBiCG: Solving for epsilon, Initial residual = 1, Final residual = 7.44946e-06, No Iterations 7
DILUPBiCG: Solving for k, Initial residual = 0.970484, Final residual = 3.50804e-07, No Iterations 1
ExecutionTime = 1.49 s ClockTime = 2 s

fieldAverage fieldAverage1 output:
Calculating averages

Courant Number mean: 0.046579 max: 6.53391
Max Ur Courant Number = 39.2927


All times are GMT -4. The time now is 22:10.