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

ysun · October 17, 2013, 14:08

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