alphaEqn.H in twoPhaseEulerFoam
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 particleparticle stress to avoid higher volumetric concentration than alphaMax. and this artificial particleparticle 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 
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. 
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 particleparticle 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 particleparticle idea is that through modification of Ua, we hope to constrain the alpha < alphaMax. Is that correct? Zhen Quote:

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 multifluid solvers, since the particle pressure term has a strongly nonlinear 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, 
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 particleparticle 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 particleparticle 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 particleparticle stress once, which is in the laplacian term: Code:
if (g0.value() > 0.0) 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 Quote:

Hi,
rUaAf is zero only at the first timestep, since rUaA is initialized to zero. If you look at pEqn.H, you have Code:
rUaAf = fvc::interpolate(rUaA); Code:
if (g0.value() > 0.0) Code:
Courant Number mean: 0.0450759 max: 0.0451754 Best, 
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 Quote:

You are welcome :)

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.90098e17, Final residual = 3.90098e17, 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.90098e17, Final residual = 3.90098e17, 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) = 1e15 kinTheory: min(nua) = 5.06848e07, max(nua) = 0.04 kinTheory: min(pa) = 0, max(pa) = 2.30439e06 GAMG: Solving for p, Initial residual = 1, Final residual = 9.52665e09, No Iterations 23 time step continuity errors : sum local = 2.15521e12, global = 2.15521e12, cumulative = 2.15521e12 max(rUaAf) = 0 min(rUaAf) = 0 max(ppMagf) = 0 min(ppMagf) = 0 DILUPBiCG: Solving for alpha, Initial residual = 0.00104655, Final residual = 1.46481e11, 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.06684e05, Final residual = 8.19921e12, No Iterations 3 Dispersed phase volume fraction = 0.2015 Min(alpha) = 1.64854e21 Max(alpha) = 0.620052 GAMG: Solving for p, Initial residual = 0.0786421, Final residual = 9.35372e09, No Iterations 17 time step continuity errors : sum local = 2.83579e11, global = 2.83579e11, cumulative = 2.62027e11 DILUPBiCG: Solving for epsilon, Initial residual = 0.135913, Final residual = 2.09582e06, No Iterations 4 DILUPBiCG: Solving for k, Initial residual = 1, Final residual = 1.43273e06, 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.48925e11, No Iterations 3 Dispersed phase volume fraction = 0.2015 Min(alpha) = 8.33627e28 Max(alpha) = 0.620094 max(rUaAf) = 0 min(rUaAf) = 0 max(ppMagf) = 0 min(ppMagf) = 0 DILUPBiCG: Solving for alpha, Initial residual = 4.58883e05, Final residual = 5.27276e12, No Iterations 3 Dispersed phase volume fraction = 0.2015 Min(alpha) = 1.55414e157 Max(alpha) = 0.620094 pressure gradient in UEqn= (0.081 0 0) kinTheory: max(Theta) = 58.918 kinTheory: min(nua) = 1.15627e12, max(nua) = 0.04 kinTheory: min(pa) = 3.88536e169, max(pa) = 2209.46 GAMG: Solving for p, Initial residual = 0.100408, Final residual = 8.11056e09, No Iterations 15 time step continuity errors : sum local = 2.49503e11, global = 2.49503e11, cumulative = 1.25243e12 max(rUaAf) = 0 min(rUaAf) = 0 max(ppMagf) = 0 min(ppMagf) = 0 DILUPBiCG: Solving for alpha, Initial residual = 0.000897803, Final residual = 1.59916e12, No Iterations 4 Dispersed phase volume fraction = 0.2015 Min(alpha) = 1.45549e14 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.37526e11, No Iterations 3 Dispersed phase volume fraction = 0.2015 Min(alpha) = 1.26941e26 Max(alpha) = 0.621929 GAMG: Solving for p, Initial residual = 0.034853, Final residual = 5.40877e09, No Iterations 17 time step continuity errors : sum local = 2.02786e11, global = 2.02786e11, cumulative = 2.15311e11 DILUPBiCG: Solving for epsilon, Initial residual = 0.0187545, Final residual = 8.89096e07, No Iterations 3 DILUPBiCG: Solving for k, Initial residual = 0.800558, Final residual = 3.72508e06, 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.59677e12, No Iterations 4 Dispersed phase volume fraction = 0.2015 Min(alpha) = 6.81511e22 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.52428e11, No Iterations 3 Dispersed phase volume fraction = 0.2015 Min(alpha) = 7.91657e26 Max(alpha) = 0.622328 pressure gradient in UEqn= (0.081 0 0) kinTheory: max(Theta) = 1000 kinTheory: min(nua) = 1.60279e12, max(nua) = 0.04 kinTheory: min(pa) = 3.08412e37, max(pa) = 2045.12 GAMG: Solving for p, Initial residual = 0.122225, Final residual = 7.73269e09, No Iterations 17 time step continuity errors : sum local = 2.63962e11, global = 2.63962e11, cumulative = 4.8651e12 max(rUaAf) = 0 min(rUaAf) = 0 max(ppMagf) = 0 min(ppMagf) = 0 DILUPBiCG: Solving for alpha, Initial residual = 0.000356051, Final residual = 2.09863e11, No Iterations 3 Dispersed phase volume fraction = 0.2015 Min(alpha) = 4.54596e16 Max(alpha) = 0.621036 max(rUaAf) = 0 min(rUaAf) = 0 max(ppMagf) = 0 min(ppMagf) = 0 DILUPBiCG: Solving for alpha, Initial residual = 8.30811e05, Final residual = 7.85124e13, No Iterations 3 Dispersed phase volume fraction = 0.2015 Min(alpha) = 1.08298e25 Max(alpha) = 0.621036 GAMG: Solving for p, Initial residual = 0.016629, Final residual = 6.94633e09, No Iterations 14 time step continuity errors : sum local = 2.18021e11, global = 2.18021e11, cumulative = 1.6937e11 DILUPBiCG: Solving for epsilon, Initial residual = 0.020499, Final residual = 8.62327e07, No Iterations 3 DILUPBiCG: Solving for k, Initial residual = 0.471488, Final residual = 2.08913e06, 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.278e11, No Iterations 3 Dispersed phase volume fraction = 0.2015 Min(alpha) = 6.07852e25 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.213e12, No Iterations 3 Dispersed phase volume fraction = 0.2015 Min(alpha) = 1.04731e25 Max(alpha) = 0.621349 pressure gradient in UEqn= (0.081 0 0) kinTheory: max(Theta) = 1000 kinTheory: min(nua) = 1.60279e12, max(nua) = 0.04 kinTheory: min(pa) = 7.28917e21, max(pa) = 1950.06 GAMG: Solving for p, Initial residual = 0.0905112, Final residual = 9.85605e09, No Iterations 15 time step continuity errors : sum local = 2.99527e11, global = 2.99527e11, cumulative = 4.68897e11 max(rUaAf) = 0 min(rUaAf) = 0 max(ppMagf) = 0 min(ppMagf) = 0 DILUPBiCG: Solving for alpha, Initial residual = 0.000346857, Final residual = 3.04659e12, No Iterations 3 Dispersed phase volume fraction = 0.2015 Min(alpha) = 3.26225e16 Max(alpha) = 0.621475 max(rUaAf) = 0 min(rUaAf) = 0 max(ppMagf) = 0 min(ppMagf) = 0 DILUPBiCG: Solving for alpha, Initial residual = 7.42358e05, Final residual = 1.26536e12, No Iterations 3 Dispersed phase volume fraction = 0.2015 Min(alpha) = 2.41699e25 Max(alpha) = 0.621475 GAMG: Solving for p, Initial residual = 0.0111199, Final residual = 9.2887e09, No Iterations 12 time step continuity errors : sum local = 2.6559e11, global = 2.6559e11, cumulative = 7.34486e11 DILUPBiCG: Solving for epsilon, Initial residual = 0.0216537, Final residual = 8.88015e06, No Iterations 2 DILUPBiCG: Solving for k, Initial residual = 0.272703, Final residual = 7.4626e06, 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.10609e13, No Iterations 4 Dispersed phase volume fraction = 0.2015 Min(alpha) = 4.13607e23 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.92767e12, No Iterations 3 Dispersed phase volume fraction = 0.2015 Min(alpha) = 2.2714e25 Max(alpha) = 0.621452 pressure gradient in UEqn= (0.081 0 0) kinTheory: max(Theta) = 1000 kinTheory: min(nua) = 1.60279e12, max(nua) = 0.04 kinTheory: min(pa) = 1.4887e21, max(pa) = 1860.09 GAMG: Solving for p, Initial residual = 0.0898181, Final residual = 9.60529e09, No Iterations 16 time step continuity errors : sum local = 2.78988e11, global = 2.78988e11, cumulative = 1.01347e10 max(rUaAf) = 0 min(rUaAf) = 0 max(ppMagf) = 0 min(ppMagf) = 0 DILUPBiCG: Solving for alpha, Initial residual = 0.000473163, Final residual = 2.11461e11, No Iterations 3 Dispersed phase volume fraction = 0.2015 Min(alpha) = 3.87245e18 Max(alpha) = 0.621414 max(rUaAf) = 0 min(rUaAf) = 0 max(ppMagf) = 0 min(ppMagf) = 0 DILUPBiCG: Solving for alpha, Initial residual = 6.93413e05, Final residual = 3.96824e12, No Iterations 3 Dispersed phase volume fraction = 0.2015 Min(alpha) = 3.91346e26 Max(alpha) = 0.621414 GAMG: Solving for p, Initial residual = 0.0374737, Final residual = 9.59322e09, No Iterations 13 time step continuity errors : sum local = 2.7584e11, global = 2.75837e11, cumulative = 1.28931e10 DILUPBiCG: Solving for epsilon, Initial residual = 0.0222171, Final residual = 9.84455e06, No Iterations 2 DILUPBiCG: Solving for k, Initial residual = 0.196096, Final residual = 4.79474e06, 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.33643e11, No Iterations 4 Dispersed phase volume fraction = 0.2015 Min(alpha) = 3.90035e25 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.15395e12, No Iterations 3 Dispersed phase volume fraction = 0.2015 Min(alpha) = 1.56767e25 Max(alpha) = 0.621334 pressure gradient in UEqn= (0.081 0 0) kinTheory: max(Theta) = 1000 kinTheory: min(nua) = 1.60279e12, max(nua) = 0.04 kinTheory: min(pa) = 5.56552e23, max(pa) = 1778.39 GAMG: Solving for p, Initial residual = 0.0744853, Final residual = 8.1907e09, No Iterations 17 time step continuity errors : sum local = 2.52555e11, global = 2.52555e11, cumulative = 1.54187e10 max(rUaAf) = 0 min(rUaAf) = 0 max(ppMagf) = 0 min(ppMagf) = 0 DILUPBiCG: Solving for alpha, Initial residual = 0.000528506, Final residual = 3.56838e11, No Iterations 3 Dispersed phase volume fraction = 0.2015 Min(alpha) = 4.19315e17 Max(alpha) = 0.621249 max(rUaAf) = 0 min(rUaAf) = 0 max(ppMagf) = 0 min(ppMagf) = 0 DILUPBiCG: Solving for alpha, Initial residual = 9.13973e05, Final residual = 4.25874e12, No Iterations 3 Dispersed phase volume fraction = 0.2015 Min(alpha) = 4.63227e25 Max(alpha) = 0.621248 GAMG: Solving for p, Initial residual = 0.0460487, Final residual = 7.05426e09, No Iterations 15 time step continuity errors : sum local = 2.20093e11, global = 2.20093e11, cumulative = 1.76196e10 DILUPBiCG: Solving for epsilon, Initial residual = 0.0228629, Final residual = 9.47428e06, No Iterations 2 DILUPBiCG: Solving for k, Initial residual = 0.155566, Final residual = 3.88173e06, 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.31724e12, No Iterations 4 Dispersed phase volume fraction = 0.2015 Min(alpha) = 1.087e24 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.34103e12, No Iterations 3 Dispersed phase volume fraction = 0.2015 Min(alpha) = 2.60664e24 Max(alpha) = 0.621241 pressure gradient in UEqn= (0.081 0 0) kinTheory: max(Theta) = 1000 kinTheory: min(nua) = 1.60279e12, max(nua) = 0.04 kinTheory: min(pa) = 2.89806e19, max(pa) = 1716.94 GAMG: Solving for p, Initial residual = 0.0786497, Final residual = 8.33953e09, No Iterations 16 time step continuity errors : sum local = 2.72705e11, global = 2.72705e11, cumulative = 2.03467e10 max(rUaAf) = 0 min(rUaAf) = 0 max(ppMagf) = 0 min(ppMagf) = 0 DILUPBiCG: Solving for alpha, Initial residual = 0.000359687, Final residual = 2.1303e11, No Iterations 3 Dispersed phase volume fraction = 0.2015 Min(alpha) = 7.1831e18 Max(alpha) = 0.621369 max(rUaAf) = 0 min(rUaAf) = 0 max(ppMagf) = 0 min(ppMagf) = 0 DILUPBiCG: Solving for alpha, Initial residual = 6.24967e05, Final residual = 2.79162e12, No Iterations 3 Dispersed phase volume fraction = 0.2015 Min(alpha) = 6.99314e28 Max(alpha) = 0.621369 GAMG: Solving for p, Initial residual = 0.0277413, Final residual = 6.64361e09, No Iterations 14 time step continuity errors : sum local = 2.22131e11, global = 2.22131e11, cumulative = 2.2568e10 DILUPBiCG: Solving for epsilon, Initial residual = 0.0237308, Final residual = 5.02822e07, No Iterations 3 DILUPBiCG: Solving for k, Initial residual = 0.130643, Final residual = 3.80861e06, 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? :confused: Quote:

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, 
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 Quote:

It turns out they added #include "pimpleControl.H"
Hi,
I've installed the openFOAM2.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 OpenFOAM1.7.1. Best Zhen 
Quote:
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 
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 Quote:

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

All times are GMT 4. The time now is 13:42. 