[qi.yang@polimi.it]

It is so weird. When I used the mesh with larger boundary layer thickness, the result is reasonable and accurate. After I decreased the boundary layer thickness, it was crashed cause the bounding k and epsilon were infinite. Who can give me some suggestions? Thanks!

[Cyp]

Originally Posted by qi.yang@polimi.it

I think I solved the problem with the following code:
- fvc::div(fvc::grad(alpha1)*rho1*phase2.turbulence( ).nut()/sigma*U1)
it was compiled successfully and I ran the test case. The results seem resonable.
However, I doubt that whether I need write the code like:
- fvc::Sp(div(fvc::grad(alpha1)*rho1*phase2.turbulen ce().nut()/sigma, U1))

Because I found the following two codes will lead different results.
solve( fvm::laplacian(k, T) + fvc::Sp(sp, T) );
solve( fvm::laplacian(k, T) + sp*T );

However, when I wrote - fvm::Sp(div(fvc::grad(alpha1)*rho1*phase2.turbulen ce().nut()/sigma, U1)), it failed to be compiled.

it failed to compile because the operator fvc::div() needs fluxes at the face center (so, a surfaceScalarField and not a volVectorField). When you do fvc::grad(alpha1), you create a volVectorField. when you write phase2.turbulence().nut(), you generate a volScalarField, where your values are computed at the cell centers. This is why you have to interpolate.

[qi.yang@polimi.it]

Originally Posted by Cyp

it failed to compile because the operator fvc::div() needs fluxes at the face center (so, a surfaceScalarField and not a volVectorField). When you do fvc::grad(alpha1), you create a volVectorField. when you write phase2.turbulence().nut(), you generate a volScalarField, where your values are computed at the cell centers. This is why you have to interpolate.

Thanks. In your opinion,
- fvc::div(fvc::grad(alpha1)*rho1*phase2.turbulence( ).nut()/sigma*U1) this code is reasonable or not? After I added this code in the equation, it can be compiled and run well for the corse mesh. However, when I use a lower thickness of boundary layer, the simulation crashed.

[qi.yang@polimi.it]

-fvm::div(fvc::snGrad(alpha1)*mesh.magSf()*fvc::int erpolate(rho1*phase2.turbulence().nut())/sigma, U1)

I changed it to be like this, it compiled. But the simulation is also crashed for the fine boundary layer mesh.

[qi.yang@polimi.it]

Originally Posted by qi.yang@polimi.it

-fvm::div(fvc::snGrad(alpha1)*mesh.magSf()*fvc::int erpolate(rho1*phase2.turbulence().nut())/sigma, U1)

I changed it to be like this, it compiled. But the simulation is also crashed for the fine boundary layer mesh.

My obejective is to add the phase dispersion term in both U and p equations
Please help me to verify the codes I wrote are right or not?
In Uequation I wrote like this:
U1Eqn =
(
fvm::div(alphaRhoPhi1, U1) - fvm::Sp(fvc::div(alphaRhoPhi1), U1)
+ MRF.DDt(alpha1*rho1, U1)
+ - fvm::laplacian(alpha1*rho1*phase2.turbulence().nuE ff(), U1)
- fvc::div(fvc::grad(alpha1)*rho1*phase2.turbulence( ).nut()/sigma*U1)

While in p equation I added one code in p equacomp1 and equacomp2:
pEqnComp1 =
(
contErr1
- fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
)/rho1
+ (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh)) - fvc::laplacian(phase2.turbulence().nut()/sigma,alpha1);

[sharonyue]

The equation looks not correct. Perhaps the first term should be a \nabla\cdot instead of \nabla. Meanwhile,
1) to add a diffusion term, you better to add that after MULES to ensure boundedness. See how it does in driftFluxFoam.
2) There is already a turbulent dispersion force in E-E model and it was implemented already, see "turbulent dispersion force".
3) The momentum interfacial exchange term needs to be addressed in UEqn, not pEqn.
4) div(grad()) employs an extended stencil instead of a compact stencil, which means it may introduce possible oscillations.

[qi.yang@polimi.it]

Thanks a lot Dongyue. The equations that I wrote here are not accurate and you are right.
We will not use the dispersion force defined in the openfoam because our group provided a new two fluid model that add two dispersion terms in both continuity and momentum equations which makes numerical stability and fast calculating.
In terms of "div(grad()) employs an extended stencil instead of a compact stencil", how to solve the problem you mentioned. However we have to add the div(grad(alpha))