
[Sponsors] 
October 9, 2019, 16:57 
Derivation of the PISO algorithm in OF

#1 
New Member
Join Date: Sep 2019
Posts: 18
Rep Power: 2 
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. Code:
{ /**********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 timestep. ); constrainPressure(p, U, phiHbyA, rAUf); //Update the pressure BCs to ensure flux consistency. //Insure mass conservation by adjusting the incoming and outcoming flux if the BC are illdefined while (piso.correctNonOrthogonal()) //Nonorthogonal 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 } } Now, in PISO loop is is splitted up the matrix in two matrices: one with the diagonal elements and the other one with offdiagonal 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)? Last edited by rucky96; October 9, 2019 at 17:58. Reason: I was wrong placing the (2) 

October 19, 2019, 00:33 

#2 
New Member
Weiwen Zhao
Join Date: Dec 2013
Posts: 19
Rep Power: 8 
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. 

October 21, 2019, 03:16 

#3 
New Member
Join Date: Sep 2019
Posts: 18
Rep Power: 2 
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: Code:
fvScalarMatrix pEqn ( fvm::laplacian(rAUf, p) == fvc::div(phiHbyA) ); 

October 23, 2019, 04:45 

#4 
New Member
Weiwen Zhao
Join Date: Dec 2013
Posts: 19
Rep Power: 8 
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). 

Tags 
mhdfoam, piso 
Thread Tools  Search this Thread 
Display Modes  


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