CFD Online Logo CFD Online URL
Home > Forums > Software User Forums > OpenFOAM

Derivation of the PISO algorithm in OF

Register Blogs Members List Search Today's Posts Mark Forums Read

LinkBack Thread Tools Search this Thread Display Modes
Old   October 9, 2019, 15:57
Default Derivation of the PISO algorithm in OF
New Member
Join Date: Sep 2019
Posts: 4
Rep Power: 2
rucky96 is on a distinguished road
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 \mathbf{u}. For the x component one has:
\dfrac{\partial U_x}{\partial t}+\boldsymbol\nabla\left(\mathbf u U_x\right)-\boldsymbol\nabla^2\left(\nu U_x\right)+\boldsymbol\nabla\left(\dfrac{\mathbf B B}{\rho\mu}\right)-\dfrac{\partial }{\partial x}\left(\dfrac{\mathbf B^2}{2\rho\mu}\right)=-\dfrac{\partial p}{\partial x}.
In the solver they sum a divergence with a gradient, so I will close my eyes and think that this is only notation.
            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.

            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.
                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.
                    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);//

                    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 M_x\cdot U_x-S_x=\boldsymbol\nabla_xp (with the sub-index I refer to the fact that is the system of equations just for x and in \boldsymbol\nabla_xp is the partial derivative of pressure respect to x). Here one must build a new matrix of coefficients, M_x, that is calculated using an appropriate linear solver. The array S contains all the source terms (i.e., the terms due to Lorentz force) except for the pressure gradient. U_x=\left(U_{x,1}, ..., U_{x,i},...,U_{x,n}\right) (n is the number of cells) the discretized velocity.

Now, in PISO loop is is splitted up the matrix M_x in two matrices: one with the diagonal elements and the other one with off-diagonal elements, i.e., M_x=A_x+H_x^\prime. Thus:
\left(A_x+H_x^\prime\right)\cdot U_x-S_x=\boldsymbol\nabla_xp,

A_x\cdot U_x+\underbrace{H_x^\prime\cdot U_x-S_x}_{H_x}=\boldsymbol\nabla_xp,

A_x\cdot U_x=\boldsymbol\nabla_xp-H_x,

U_x=A_x^{-1}\cdot\boldsymbol\nabla_xp-A_x^{-1}\cdot H_x,
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
\boldsymbol\nabla\left(A_x^{-1}\cdot\boldsymbol\nabla_xp\right)=\boldsymbol\nabla\left(A_x^{-1}\cdot H_x\right)\ \ \ (1).

The condition that must be met is
\boldsymbol\nabla\cdot\mathbf u=0=\dfrac{\partial U_x}{\partial x}+\dfrac{\partial U_y}{\partial y}+\dfrac{\partial U_z}{\partial z},
\dfrac{\partial}{\partial x}\left(A_x^{-1}\cdot\boldsymbol\nabla_xp-A_x^{-1}\cdot H_x\right)
+\dfrac{\partial}{\partial y}\left(A_y^{-1}\cdot\boldsymbol\nabla_yp-A_y^{-1}\cdot H_y\right)
+\dfrac{\partial}{\partial z}\left(A_z^{-1}\cdot\boldsymbol\nabla_zp-A_z^{-1}\cdot H_z\right)=0,\ \ \ (2)
not just \dfrac{\partial U_x}{\partial x}=0.

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

Last edited by rucky96; October 9, 2019 at 16:58. Reason: I was wrong placing the (2)
rucky96 is offline   Reply With Quote

Old   October 18, 2019, 23:33
New Member
Weiwen Zhao
Join Date: Dec 2013
Posts: 18
Rep Power: 7
wwzhao is on a distinguished road
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.
wwzhao is offline   Reply With Quote

Old   Today, 02:16
New Member
Join Date: Sep 2019
Posts: 4
Rep Power: 2
rucky96 is on a distinguished road
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 \left(1\right) without explaining anything. What I am asking is if (without saying so) they take the velocity vector, \mathbf u= U_x\boldsymbol\imath+U_y\boldsymbol\jmath+U_z\boldsymbol k, and apply the divergence on it to get \left(2\right). 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)
rucky96 is offline   Reply With Quote


mhdfoam, piso

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On

Similar Threads
Thread Thread Starter Forum Replies Last Post
Development of PISO Algorithm Roger9318 Main CFD Forum 1 September 12, 2018 07:51
Eulerian multiphase model with PISO algorithm Shehryar Fluent Multiphase 0 July 6, 2017 02:13
Doubt Regarding the Predictor step in PISO algorithm shadabdyn OpenFOAM Programming & Development 0 February 12, 2017 01:41
Time step size sensitivity of piso algorithm TypeR OpenFOAM 1 February 1, 2017 16:00
Non-linearity Pressure Equation -- PISO algorithm gdeneyer OpenFOAM Programming & Development 1 August 23, 2012 05:19

All times are GMT -4. The time now is 08:09.