[shadabdyn]

Hello everyone,

I need to add the gravity settling effect delrho*g*phi to the existing self-made code. and I am using the Piso algorithm. so I added delrho*g*phi at the solve momentum predictor step. here I choose g as a volume vector defined in the create field properties simultaneously additional initial property in 0 folder of the case directory and density difference as constant transport properties ( dimensional scalar in create property.h file). when I compile the Source code it is running fine without any error..when I run twice it is showing up to date.
So I am asking that what I have done is correct or not. I am adding my source code also below. I am attaching a picture of my governing equation also.

// Pressure-velocity PISO corrector
for (int corr=0; corr<nCorr; corr++)
{
// Momentum predictor

tmp<fvVectorMatrix> UEqn
(

- fvm::laplacian(nu1, U) -fvc::div(nu1*dev(fvc::grad(U)().T())) -fvc::div(SigmaL)
);

UEqn().relax();

solve(UEqn() == -fvc::grad(p)+delrhog*c); //+delrhog*c

p.boundaryField().updateCoeffs();
volScalarField rUA = 1.0/UEqn().A();
U = rUA*UEqn().H();
UEqn.clear();
phi = fvc::interpolate(U) & mesh.Sf();
adjustPhi(phi, U, p);

// Store pressure for under-relaxation
p.storePrevIter();

// Non-orthogonal pressure corrector loop
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(rUA, p) == fvc::div(phi)
);

pEqn.setReference(pRefCell, pRefValue);
pEqn.solve();

if (nonOrth == nNonOrthCorr)
{
phi -= pEqn.flux();
}
}

# include "continuityErrs.H"

// Explicitly relax pressure for momentum corrector
p.relax();

// Momentum corrector
U -= rUA*fvc::grad(p);
U.correctBoundaryConditions();
}


up to best of my knowledge we should'nt add I source term at momentum corrector step?
Thanks in advance.

[pbachant]

Looks like your image isn't showing up. Also, please put your code in the [CODE] blocks so it's in monospace font with indentation.

[shadabdyn]

Thank u pbachant for your reply. I hope now this image and code could help u to understand my problem.

        // Pressure-velocity SIMPLE corrector
        for (int corr=0; corr<nCorr; corr++)
        {



            // Momentum predictor

            tmp<fvVectorMatrix> UEqn
            (
    
               - fvm::laplacian(nu1, U) -fvc::div(nu1*dev(fvc::grad(U)().T())) -fvc::div(SigmaL)
            );

            UEqn().relax();

            solve(UEqn() == -fvc::grad(p)+delrhog*c);  //+delrhog*c

            p.boundaryField().updateCoeffs();
            volScalarField rUA = 1.0/UEqn().A();
            U = rUA*UEqn().H();
            UEqn.clear();
            phi = fvc::interpolate(U) & mesh.Sf();
            adjustPhi(phi, U, p);

            // Store pressure for under-relaxation
            p.storePrevIter();

            // Non-orthogonal pressure corrector loop
            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
            {
                fvScalarMatrix pEqn
                (
                    fvm::laplacian(rUA, p) == fvc::div(phi)
                );

                pEqn.setReference(pRefCell, pRefValue);
                pEqn.solve();

                if (nonOrth == nNonOrthCorr)
                {
                    phi -= pEqn.flux();
                }
            }

#           include "continuityErrs.H"

            // Explicitly relax pressure for momentum corrector
            p.relax();

            // Momentum corrector
            U -= rUA*fvc::grad(p);
            U.correctBoundaryConditions();
        }