[rucky96]

Hi Foamers,

I am trying to understand the PISO algorithm. I know taht is to work out the momentum equation and obtain the pressure.

I am working with mhdFoam so is the PISO for this solver. The momentum equation in OF is the one for a single component of velocity . For the component one has:

In the solver they sum a divergence with a gradient, so I will close my eyes and think that this is only notation.

        {
        	/**********MOMENTUM*EQUATION**********/
            fvVectorMatrix UEqn //:confused:
            (
                fvm::ddt(U) //Time derivative.
              + fvm::div(phi, U) //Convection.
              - fvm::laplacian(nu, U) //Advection.
              + fvc::div(phiB, 2.0*DBU*B) //Lorentz force term 1.
              - fvc::grad(DBU*magSqr(B)) //Lorentz force term 2.
            );//DBU=1/(2*mu*rho),

            if (piso.momentumPredictor()) //Prediction of the velocity.
            {
                solve(UEqn ==   - fvc::grad(p) //Pressure term.
                                );
            }
            // --- Standard PISO loop starts.
            while (piso.correct())
            {
                volScalarField rAU(1.0/UEqn.A()); //The matrix A^(-1) is saved as rAU. Scalar values at centres.
                surfaceScalarField rAUf //Scalar values at faces.
                (
                    "rAUf",
                    fvc::interpolate(rAU)
                ); 
                volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p)); //Evaluates relation A^(-1)*H and saves to HbyA. Vector values at centres.
                surfaceScalarField phiHbyA //Scalar values at faces.
                (
                    "phiHbyA",
                    fvc::flux(HbyA) //Interpolation to faces and dot product with respective face normal vector
                  + rAUf*fvc::ddtCorr(U, phi) //Correction related to the times scheme used to avoid an unwanted dependency on the time-step.
                );

                constrainPressure(p, U, phiHbyA, rAUf); //Update the pressure BCs to ensure flux consistency.
                //Insure mass conservation by adjusting the in-coming and out-coming flux if the BC are ill-defined
                while (piso.correctNonOrthogonal()) //Non-orthogonal pressure corrector loop.
                {
                    fvScalarMatrix pEqn //Pressure corrector
                    (
                        fvm::laplacian(rAUf, p) == fvc::div(phiHbyA)
                    );

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

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

                #include "continuityErrs.H" //Including header file for continuity error evaluation
                U = HbyA - rAU*fvc::grad(p); //U has been calculated from p.
                U.correctBoundaryConditions();//The BCs do no longer correspond to the ones in 0/U. This function means that the BCs of U must be the ones precised in 0/U.
                //End of current time step
            }
        }

With the discretization of the momentum equation, using the FVM techniques for implicit discretization of the velocity field (the fvm::ddt(U), fvm::div(phi, U) and fvm::laplacian(nu, U)) and explicit discretization schemes for the source terms (fvc::grad(p), fvc::div(phiB, 2.0*DBU*B) and fvc::grad(DBU*magSqr(B))) , one can approximate the momentum equation as three linear system of equations (one for each velocity component) of the kind (with the sub-index I refer to the fact that is the system of equations just for and in is the partial derivative of pressure respect to ). Here one must build a new matrix of coefficients, , that is calculated using an appropriate linear solver. The array contains all the source terms (i.e., the terms due to Lorentz force) except for the pressure gradient. ( is the number of cells) the discretized velocity.

Now, in PISO loop is is splitted up the matrix in two matrices: one with the diagonal elements and the other one with off-diagonal elements, i.e., . Thus:







So far everything is great, my problem is when some people apply the divergence on the previous equation (yes, on a array of scalars) and say that their divergence has to be zero to get
.

The condition that must be met is

i.e.,



not just .

My question is then what equation is OpenFoam solving: (1) or (2)?

[wwzhao]

Hi rucky96,

Eq.1 in your post is wrong. Divergence should be applied to a vector field rather than each component of the vector field.

[rucky96]

Hi wwzhao, thanks for answering.

Of course neither me nor anyone can apply a divergence over a scalar. That is just the justification that some people use. They jump from the previous equation to without explaining anything. What I am asking is if (without saying so) they take the velocity vector, , and apply the divergence on it to get . And, then, if equation (2) is what mhdFoam (and icoFOAM and all solvers with PISO) is really solving here:

                    fvScalarMatrix pEqn
                    (
                        fvm::laplacian(rAUf, p) == fvc::div(phiHbyA)
                    );

[wwzhao]

Hi rucky96,

If you replace \nabla with \nabla \cdot outside the parentheses in (1) and rewrite it in vectorial form, (1) is equivalent to (2).

(2) is easy to calculate in Cartesian grids. However, OpenFOAM adopts unstructured polyhedral mesh. So the procedure for calculating divergence is a bit different from (2).