[Houthuys]

Hi All,

I've read several discussions in this forum about including the gravity body force in pisoFoam or similar solvers. The most successful approach seems to be the one implemented in interFoam, i.e. using the following term in UEqn:

 ghf*fvc::grad(rho)

However, this approach obviously not suitable for the case of single phase incompressible flow with constant density.

So, I tried to get through this problem and I've modified icoFoam using two approaches and I want to discuss them here and see if there's any better approach or if we can modify them to get better.

In the first approach I've just add the gravity body force to the UEqn and I didn't change any other parts in the code. The new UEqn will be like:

        fvVectorMatrix UEqn 
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)==fvc::grad(gh)
        );

I didn't change the pressure equation since the H(U) now includes the contribution from the gravity body force. Then I test using cavity tutorial by setting the moving wall boundary conditions as p: fixedValue = 0, U: fixedValue = 0 0 0, since I want to reproduce hydro-static pressure.

The results for the pressure was almost like the analytical solution (hydro-static pressure distribution) while unfortunately I get some velocities near the upper and lower walls while velocities should be zero everywhere. You can see the results in the attachment.

In the second approach I just combined the gravity force with pressure and called it p_gh where:

p_gh = p-gh
-grad(p)+grad(gh) = -grad(p-gh)=-fvc::grad(p_gh)

So, basically I didn't change anything in the code. I only I dded one line at the end of the time loop to calculate p from P_gh. The code looks like:

 while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        #include "CourantNo.H"  

        // Momentum predictor
        fvVectorMatrix UEqn 
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
        );

        if (piso.momentumPredictor()) 
        {
            solve(UEqn == -fvc::grad(p_gh)); 
        } 

        // --- PISO loop starts here
        while (piso.correct())
        {

	   volScalarField rAU(1.0/UEqn.A()); 
           
            volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_gh));
            surfaceScalarField phiHbyA
            (
                "phiHbyA",
                fvc::flux(HbyA)
              + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
            );

            adjustPhi(phiHbyA, U, p_gh);

            // Update the pressure BCs to ensure flux consistency
            constrainPressure(p_gh, U, phiHbyA, rAU);

            // Non-orthogonal pressure corrector loop
            while (piso.correctNonOrthogonal())
            {
                // Pressure corrector

                fvScalarMatrix pEqn
                (
                    fvm::laplacian(rAU, p_gh) == fvc::div(phiHbyA)
                );

                pEqn.setReference(pRefCell, pRefValue);

                pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));

                if (piso.finalNonOrthogonalIter())
                {
                    phi = phiHbyA - pEqn.flux();
                }
            }

            #include "continuityErrs.H"

            U = HbyA - rAU*fvc::grad(p_gh);
            U.correctBoundaryConditions();
	    p = p_gh + gh;
        }

The result in this case for p and U are perfect (now U is zero everywhere), you can see this results also in the attachments.
The only problem with this approach is that you have to set the pressure b.c carefully for p_gh.

So, in your opinion, is the second approach is suitable in general? does it need improvement? or should we use totally different approach?

[jherb]

Did you look at the implementation of buoyantBoussinesqSimpleFoam (buoyantBoussinesqPimpleFoam)? This should do what you want. Or perhaps twoLiquidMixingFoam.

[Houthuys]

Originally Posted by jherb

Did you look at the implementation of buoyantBoussinesqSimpleFoam (buoyantBoussinesqPimpleFoam)? This should do what you want. Or perhaps twoLiquidMixingFoam.

Thanks Joachim for the reply.

Actually both solvers use density gradient to introduce the gravity body force, e.g.:

- ghf*fvc::snGrad(rho)

In my case I have constant density. So, how can we add that effect in this case?

[jherb]

Sorry for the late reply. But if you have constant density, what should be the effect of gravity?

[Santiago]

I don't see the point of putting gravity if you are not solving a variable density, or buoyancy-driven, flow

Sent from my GT-I8190L using CFD Online Forum mobile app

[fedez91]

Hi rick,

I've been trying to do the same. I don't really understand what do you mean by "gh". Is it only g? is g*z (if z is your vertical direction)?

I have implemented g in icoFoam adding g to transport properties as a vector and in the U equation as follows,

 
        fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
          - g
        );

        if (piso.momentumPredictor())
        {
            solve(UEqn == -fvc::grad(p));
        }

but I got the same error as you do (some velocities near the wall in the cavity test).

Where you able to solve this problem? Could you please explain me what did you mean by gh and how did you defined it in the case?

(Sorry for the silly questions, first time trying to programming in OF)

Best,

[Santiago]

Originally Posted by fedez91

Hi rick,

I've been trying to do the same. I don't really understand what do you mean by "gh". Is it only g? is g*z (if z is your vertical direction)?

I have implemented g in icoFoam adding g to transport properties as a vector and in the U equation as follows,

 
        fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
          - g
        );

        if (piso.momentumPredictor())
        {
            solve(UEqn == -fvc::grad(p));
        }

but I got the same error as you do (some velocities near the wall in the cavity test).

Where you able to solve this problem? Could you please explain me what did you mean by gh and how did you defined it in the case?

(Sorry for the silly questions, first time trying to programming in OF)

Best,

Youre explicitly adding a source to a trabsport eqn that, in principle, is time-discretized in an implicit form. You should add the term to the diagonal of the resulting linear system, via the fvm::Sp method

[Mahmoud Abbaszadeh]

Originally Posted by jherb

Sorry for the late reply. But if you have constant density, what should be the effect of gravity?

to see the effects of the hydrostatic pressure