CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

phi -= pEqn.flux() vs. linearInterpolate(U) & mesh.Sf()

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

Like Tree17Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   January 5, 2011, 22:57
Default phi -= pEqn.flux() vs. linearInterpolate(U) & mesh.Sf()
  #1
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 14
santiagomarquezd will become famous soon enough
Hi all FOAMers, I'm wondering why final phi at the end of each timestep in icoFoam is calculated as:

a)
Code:
phi -= pEqn.flux();
instead of doing:

b)
Code:
phi=linearInterpolate(U) & mesh.Sf();
particularly having in mind that at this point U is already corrected by PISO.

I've checked the values of both methods and they are slightly different, maybe with a difference of 10-30%. It's interesting that, if you are running continued, phi for the next timestep is calculated like a), but if you stop running and then restart phi for the next timestep then phi is calculated following b).


It should be equal? If not, What is correct and why? Why one is used for continued runnings and the other one for restarting?

Thx in advance.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   January 6, 2011, 06:08
Default
  #2
Senior Member
 
Steven van Haren
Join Date: Aug 2010
Location: The Netherlands
Posts: 149
Rep Power: 6
stevenvanharen is on a distinguished road
Hi Santiago,

please keep in mind that -= means subtraction.

So, in step a phi is not completely recalculated, but subtracted with pEqn.flux();

For a complete understanding of the procedure for solving the Navier-Stokes equation I suggest you read pages 143 till and with 146 of Jasak's thesis:

http://powerlab.fsb.hr/ped/kturbo/Op...jeJasakPhD.pdf

Cheers,

Steven
pyt likes this.
stevenvanharen is offline   Reply With Quote

Old   January 6, 2011, 06:33
Default
  #3
Cyp
Senior Member
 
Cyprien
Join Date: Feb 2010
Location: Stanford University
Posts: 227
Rep Power: 9
Cyp is on a distinguished road
I was wondering the same question. In my own experiment the method with linearInterpolate(U)&mesh.Sf() can yields in very erroneous results...

Regards,
Cyp
Cyp is offline   Reply With Quote

Old   January 6, 2011, 11:58
Default
  #4
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 14
santiagomarquezd will become famous soon enough
Quote:
Originally Posted by stevenvanharen View Post
Hi Santiago,

please keep in mind that -= means subtraction.

So, in step a phi is not completely recalculated, but subtracted with pEqn.flux();

For a complete understanding of the procedure for solving the Navier-Stokes equation I suggest you read pages 143 till and with 146 of Jasak's thesis:

http://powerlab.fsb.hr/ped/kturbo/Op...jeJasakPhD.pdf

Cheers,

Steven
Tnx for the answers. Yes I know that part, the question is in CyP's way, because in both cases you're trying to assemble a corrected flux, moreover, as I pointed out if you do a continued running a) method is used, but in restart b) method is used instead. The question is if the two methods are different and why and why each method have to used for restarting/continuing. As CyP mentioned using method b in continued runnings leads to wrong results.

My advisor is implementing icoFoam in MatLab and is assembling phi at the beginning of each timestep by method b), things aren't going well that way.

Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   January 6, 2011, 12:18
Default
  #5
Senior Member
 
Steven van Haren
Join Date: Aug 2010
Location: The Netherlands
Posts: 149
Rep Power: 6
stevenvanharen is on a distinguished road
Hi Santiago,

Ok, but why do you think method a) is used in continued running and method b) in restart?

If I look in icoFoam I don't see any differences for restarting or not.

at line 72:

Code:
phi = (fvc::interpolate(U) & mesh.Sf())
                + fvc::ddtPhiCorr(rUA, U, phi);
and at line 89:

Code:
phi -= pEqn.flux();
Steven
stevenvanharen is offline   Reply With Quote

Old   January 6, 2011, 12:41
Default
  #6
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 14
santiagomarquezd will become famous soon enough
Steven, thanks for the interest in the topic. When you restart or start from /0 UEqn is assembled using phi calculated by createFields.H line 43 of icoFoam.C. This is calculated using method b). This phi surfaceScalarField affects the values of H operator, U predictor, etc. used afterwards. In case of continued running UEqn is assembled using phi calculated at the end of previous timestep using method a).

That's the question.

Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar

Last edited by santiagomarquezd; January 6, 2011 at 17:41. Reason: Spelling
santiagomarquezd is offline   Reply With Quote

Old   January 6, 2011, 16:08
Default
  #7
Senior Member
 
Steven van Haren
Join Date: Aug 2010
Location: The Netherlands
Posts: 149
Rep Power: 6
stevenvanharen is on a distinguished road
Hi Santiago,

Ok, at start method b is used in createFields.H indeed. This is just to get the simulation started. The non-linearity of the momentum equation is lagged using phi from the previous timestep. But since phi is usually not available at start this is created in createFields.H.

But if phi is available (usual in a restart unless you delete phi) phi is read from the timefolder. See line 46 of createFields.H:

Code:
IOobject::READ_IF_PRESENT
Now during a simulation both methods are used, first:

Code:
phi = (fvc::interpolate(U) & mesh.Sf())
               + fvc::ddtPhiCorr(rUA, U, phi);
and then:

Code:
 phi -= pEqn.flux();
So this is what your advisor should do as well.

This set of equations is not easy to understand. Especially if you compare to the original paper of Issa describing the PISO method. Reading the four pages in Jasak's thesis helped me a lot so I really urge you to do the same.

Steven
stevenvanharen is offline   Reply With Quote

Old   January 6, 2011, 17:41
Default
  #8
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 14
santiagomarquezd will become famous soon enough
Aha, it's true, phi is read if it is present, I had forgotten this point, and as you said if you delete this file and approximated phi is assembled to restarting.

Respect of using both methods while running it's true too, but the phi stored for the next time is what was calculated using method a). This method resembles Eqn. 144 from Jasak thesis, which comes from Eqn.140 & Sf. Eqn. 140 is the same of Eqn. 139 but interpolated at faces.

Checking the notation is worthy to note that Eqn. 140 is written as:

a) U_f=[H(U)/a_P]_f-(1/a_P)_f (\nabla p)_f

(the last has each factor previously interpolated)

and not

b) U_f=[H(U)/a_P]_f-[(1/a_P) (\nabla p)]_f

which one would obtain applying the interpolation the each term of the RHS in Eqn. 129. Checking the code in laplacianScheme.C we have:

Code:
00108     return fvmLaplacian(tinterpGammaScheme_().interpolate(gamma)(), vf);
so that code is implemented as in Jasak's thesis (option a). Assembling phi as

phi=U_f & S.f

leads to

phi=(U_f=[H(U)/a_P]_f-[(1/a_P) (\nabla p)]_f)&Sf

which uses b) and not a), this explains the differences we've found, but it's not clear for me and my advisor why a) is correct and not b).

Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   January 7, 2011, 08:28
Default
  #9
Senior Member
 
Steven van Haren
Join Date: Aug 2010
Location: The Netherlands
Posts: 149
Rep Power: 6
stevenvanharen is on a distinguished road
Quote:
Originally Posted by santiagomarquezd View Post
phi stored for the next time is what was calculated using method a).
I don't see this, what is actually saved is:

Code:
 phi = (fvc::interpolate(U) & mesh.Sf())
                + fvc::ddtPhiCorr(rUA, U, phi)-pEqn.flux();
So I think still a combination of method a and b.

Concerning the difference between:

[(1/a_P) (\nabla p)]_f and (1/a_P)_f (\nabla p)_f


This is the way OpenFoam discretizes a laplacian operator (see equation 3.24 in Jasak):

(1/a_P)_f (\nabla p)_f:

What about:

[(1/a_P) (\nabla p)]_f

How would you discretize this? Since now you need [(1/a_P) (\nabla p)]_f. The program would need to use 3.26 of Jasak instead of 3.25 of Jasak. This method is thus not incorrect, but undesirable. Read section 3.3.1.3 of Jasak.
stevenvanharen is offline   Reply With Quote

Old   January 7, 2011, 12:23
Default
  #10
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26
alberto will become famous soon enoughalberto will become famous soon enough
Quote:
Originally Posted by santiagomarquezd View Post
Hi all FOAMers, I'm wondering why final phi at the end of each timestep in icoFoam is calculated as:

a)
Code:
phi -= pEqn.flux();
instead of doing:

b)
Code:
phi=linearInterpolate(U) & mesh.Sf();
particularly having in mind that at this point U is already corrected by PISO.
Because a) enforces the mass conservation principle, while the second does not. Remember that the pressure equation is assembled with the face fluxes, and the flux is the conservative variable, not the velocity.

You use pEqn.flux() to obtain the exact term you have to subtract from the matrix of the pEqn, in order not to introduce conservation errors.

Best,
su_junwei, fumiya and abtinansari like this.
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image.
OpenQBMM - An open-source implementation of quadrature-based moment methods
alberto is offline   Reply With Quote

Old   January 7, 2011, 12:35
Default
  #11
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26
alberto will become famous soon enoughalberto will become famous soon enough
Quote:
Originally Posted by stevenvanharen View Post
Concerning the difference between:

[(1/a_P) (\nabla p)]_f and (1/a_P)_f (\nabla p)_f


This is the way OpenFoam discretizes a laplacian operator (see equation 3.24 in Jasak):

(1/a_P)_f (\nabla p)_f:

What about:

[(1/a_P) (\nabla p)]_f

How would you discretize this? Since now you need [(1/a_P) (\nabla p)]_f. The program would need to use 3.26 of Jasak instead of 3.25 of Jasak. This method is thus not incorrect, but undesirable. Read section 3.3.1.3 of Jasak.
I agree. Additionally, notice that OpenFOAM and many papers in the literature make the assumption

(a*b*c)_f = a_f * b_f * c_f

This is not true in general, and it holds only for a smooth functional. However it often simplifies the development of the solution procedure.

For example, if you take a look at compressible solvers, you find

phi = rho_f*U_f

while from the theory it should be

phi = (rho*U)_f

In practice, in many applications, the approximation is acceptable.

Best,
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image.
OpenQBMM - An open-source implementation of quadrature-based moment methods
alberto is offline   Reply With Quote

Old   January 7, 2011, 13:25
Default
  #12
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 14
santiagomarquezd will become famous soon enough
Hi,

Quote:
Originally Posted by stevenvanharen View Post
I don't see this, what is actually saved is:

Code:
 phi = (fvc::interpolate(U) & mesh.Sf())
                + fvc::ddtPhiCorr(rUA, U, phi)-pEqn.flux();
So I think still a combination of method a and b.
yes, talking in general both methods are used, but I'm referring to final lines of each timestep in icoFoam.C, particularly:

Code:
00088                 if (nonOrth == nNonOrthCorr)
00089                 {
00090                     phi -= pEqn.flux();
00091                 }
00092             }
00093 
00094             #include "continuityErrs.H"
00095 
00096             U -= rUA*fvc::grad(p);
00097             U.correctBoundaryConditions();
00098         }
once U in cells centers is updated you have U and phi ready for the next timestep, at this point you can generate phi using method b). Using gdb to debug an example you have, for phi calculated by a) (line 90):

phi = {-3.4645991399177079e-06, 3.4646001127746356e-06, -6.3903343759551748e-06, 2.9257355721740963e-06, ...}

and for phi by method b)

$62 = {-3.531581793505539e-06, 4.1013330984044824e-06, -5.9540249950781068e-06, 3.2691627424804398e-06, ...}

values are slightly different. I think the answer is was given by Alberto in #10. My advisor pointed me in the necessity of have a look at fluxes conservation.

Respect of:

Quote:
Originally Posted by stevenvanharen View Post
Concerning the difference between:

[(1/a_P) (\nabla p)]_f and (1/a_P)_f (\nabla p)_f


This is the way OpenFoam discretizes a laplacian operator (see equation 3.24 in Jasak):

(1/a_P)_f (\nabla p)_f:

What about:

[(1/a_P) (\nabla p)]_f

How would you discretize this? Since now you need [(1/a_P) (\nabla p)]_f. The program would need to use 3.26 of Jasak instead of 3.25 of Jasak. This method is thus not incorrect, but undesirable. Read section 3.3.1.3 of Jasak.
very good point, nevertheless this is the "algorithmic" reason, basically reduce the computational molecule. The theoretical problem is how to pass from r.h.s in the first line of 3.24 of Jasak to the second one. My concern was related to Alberto's explanation:

Quote:
Originally Posted by alberto View Post
Additionally, notice that OpenFOAM and many papers in the literature make the assumption

(a*b*c)_f = a_f * b_f * c_f

This is not true in general, and it holds only for a smooth functional. However it often simplifies the development of the solution procedure.

For example, if you take a look at compressible solvers, you find

phi = rho_f*U_f

while from the theory it should be

phi = (rho*U)_f

In practice, in many applications, the approximation is acceptable.
Alberto, do you have any references of those papers in order to go deeper? and another question, assuming this linearization What's the background difference between a) and b)?, or in other words Why b) don't "enforces the mass conservation principle" ? I think these are two different methods of assembling 3.144 of Jasak.

Thanks to you guys, discussion was really useful.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   January 7, 2011, 13:44
Default
  #13
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26
alberto will become famous soon enoughalberto will become famous soon enough
I have seen the approximation (a*b)_f = a_f*b_f used quite often in Issa papers (I have in mind a paper on multiphase flows (Oliveira & Issa, 2003, Int. J. Num. Meth. Fluids), but also elsewhere), and in extension to the Rhie-Chow formula (S. Zhang, X. Zhao, General formulation for Rhie-Chow Interpolation, ASME Heat Transfer/FLuids Engineering Summer Conference, HT-FED04-56274, Charlotte, North Carolina, July 11-15 2004).

Keep in mind that this assumption is absolutely arbitrary. If you want to be picky, you should not make that assumption, since it introduces an approximation.

I have to say it is a very convenient approximation however, and I have used it myself in multiphase codes for a simple reason. Just consider the case where you have a phase fraction 0 < alpha < 1, and you define

phi = (alpha*rho*U)_f

In this case, if you use a reconstruction procedure when correcting the velocity, you must divide by alpha_f*rho_f. This forces you to do something (typically not very clean) when alpha is small.

If you define phi as the velocity flux:

phi = (U)_f

and then consider alpha_f*rho_f*(U)_f where you need it, like in the pressure equation, your correction step is clean, and does not have any problem when alpha tends to zero. Of course if your fields are strongly non-uniforms (shocks), you need the first approach, but in incompressible solvers or slightly compressible ones, you probably won't see a significant difference.

Finally, b) does not enforce the conservation principle because you do not enforce it on U exactly, but on phi. Your pressure equation is

div(phi) = div(U_f) = 0

not

div(U_c) = 0

being U_c the cell-centered velocity. For the same reason, in transport equations you use div(phi, T).

Best,
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image.
OpenQBMM - An open-source implementation of quadrature-based moment methods
alberto is offline   Reply With Quote

Old   January 7, 2011, 14:44
Default
  #14
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 14
santiagomarquezd will become famous soon enough
Excellent!!
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   March 30, 2012, 06:28
Default variable porosity
  #15
New Member
 
Amir
Join Date: Feb 2012
Posts: 2
Rep Power: 0
amirhoshangm is on a distinguished road
Hi

I, simulating flow in a channel with variable porosity. so I need to change some thing in correction of flux. it means I have to multiply porosity to pressure. the flux in updated by

pEqn.flux()


could anyone let me know how I can use porosity*pressure instead of pressure in correction of flux?
amirhoshangm is offline   Reply With Quote

Old   April 5, 2013, 17:37
Default
  #16
Senior Member
 
Join Date: Nov 2012
Posts: 168
Rep Power: 4
hz283 is on a distinguished road
Hello,

I also had the problem about the pEqn.flux(). In these equations in your post, is pEqn.flux() returns the value [(1/a_P) (\nabla p)]_f&Sf?

In rhoPimpleFoam, there are two regimes for pressure: transonic and subsonic. In Transonic, the correcpsonding flux correction is:
phi == pEqn.flux()

while in subsonic regime, it is as follows
phi +=pEqn.flux

So from the above the sign is opposite. How can we justify this difference? I also check other solvers and found that some use "==" while some use "+=" . Does anybody know something about this?

Thanks.

Quote:
Originally Posted by santiagomarquezd View Post
Aha, it's true, phi is read if it is present, I had forgotten this point, and as you said if you delete this file and approximated phi is assembled to restarting.

Respect of using both methods while running it's true too, but the phi stored for the next time is what was calculated using method a). This method resembles Eqn. 144 from Jasak thesis, which comes from Eqn.140 & Sf. Eqn. 140 is the same of Eqn. 139 but interpolated at faces.

Checking the notation is worthy to note that Eqn. 140 is written as:

a) U_f=[H(U)/a_P]_f-(1/a_P)_f (\nabla p)_f

(the last has each factor previously interpolated)

and not

b) U_f=[H(U)/a_P]_f-[(1/a_P) (\nabla p)]_f

which one would obtain applying the interpolation the each term of the RHS in Eqn. 129. Checking the code in laplacianScheme.C we have:

Code:
00108     return fvmLaplacian(tinterpGammaScheme_().interpolate(gamma)(), vf);
so that code is implemented as in Jasak's thesis (option a). Assembling phi as

phi=U_f & S.f

leads to

phi=(U_f=[H(U)/a_P]_f-[(1/a_P) (\nabla p)]_f)&Sf

which uses b) and not a), this explains the differences we've found, but it's not clear for me and my advisor why a) is correct and not b).

Regards.
hz283 is offline   Reply With Quote

Old   May 25, 2013, 08:36
Default
  #17
New Member
 
Emad Tandis
Join Date: Sep 2010
Posts: 27
Rep Power: 6
EmadTandis is on a distinguished road
Hi
I have this problem, too. this function (pEqn.flux() ) works based on which relation in jasak's thesis?
thanks
EmadTandis is offline   Reply With Quote

Old   May 25, 2013, 08:59
Default
  #18
Cyp
Senior Member
 
Cyprien
Join Date: Feb 2010
Location: Stanford University
Posts: 227
Rep Power: 9
Cyp is on a distinguished road
Hi!

You can find in the following link a presentation I made regarding the details of icoFoam. I sure you will find all your answers there.
http://www.scribd.com/doc/143414962/...on-in-OpenFOAM

The slides are in French, but you could easily understand the mathematical and programming part

Enjoy!

Cyp
Cyp is offline   Reply With Quote

Old   May 26, 2013, 04:24
Question
  #19
New Member
 
Emad Tandis
Join Date: Sep 2010
Posts: 27
Rep Power: 6
EmadTandis is on a distinguished road
Quote:
Originally Posted by Cyp View Post
Hi!

You can find in the following link a presentation I made regarding the details of icoFoam. I sure you will find all your answers there.
http://www.scribd.com/doc/143414962/...on-in-OpenFOAM

The slides are in French, but you could easily understand the mathematical and programming part

Enjoy!

Cyp
Thanks
But there is no any note about pEqn.flux() overthere. I know this module return a surfaceScalarField type that must be subtracted form phi. But I don't Know Why?
phi -= pEqn.flux();

Last edited by EmadTandis; May 26, 2013 at 05:36.
EmadTandis is offline   Reply With Quote

Old   May 26, 2013, 05:05
Default
  #20
Senior Member
 
Pawel Sosnowski
Join Date: Mar 2009
Location: Trieste, Italy
Posts: 105
Rep Power: 9
psosnows is on a distinguished road
Hello,

Quote:
Originally Posted by EmadTandis View Post
But there is no any note about pEqn.flux() overthere.
Well, in fact there is. Check out the 5th slide (or 3rd from the end).

I think some time ago there was some explanation regarding this topic, but lets do it from the top.

A nice explanation of the way N-S is solved in OpenFOAM is presented in the OF-wiki about SIMPLE algo:
http://openfoamwiki.net/index.php/Th...hm_in_OpenFOAM
PISO algo is very similar:
http://openfoamwiki.net/index.php/Th...hm_in_OpenFOAM
Of course check also sub-topics.

Now for some notes. As we all know, OF uses Eulerial approach with non-staggered meshes. This can introduce some significant numerical errors (especially while using 2nd order central difference schemes). For that reason in OF we "do not" solve directly for velocity (U), but for the fluxes (phi).
If you check the divergence of U, you will almost always find that there is a not-negligible error. At the same time phi is guaranteed to be divergence free.

Now lets take a look into the PISO algo itself. I will not cite the code but it should be easy to follow.
1) create the UEqn using the last- known phi. Note that any "XEqn" is like a black box, with a void space for the unknown. The imporatant stuff are the coeffs.
2) if you like, solve momentum predictor (this is unnecessary and can be dropped for time saving).
3) extract the semi-central coeffs from the UEqn, reverse them and call rAU. This is the famous "operator splitting".
4) pressure loop:
4.1) recalculate velocity: U = rAU * H();
4.2) recalculate fluxes: phi = interpolate(U) * S; note that the flux field is NOT divergence free;
4.3) solve for pressure using the non divergence free flux
4.4) solve it several times...
5) now, we solved for pressure with non-div-free phi. But the literature (Jasak's thesis, Issa et al.) tells us, that we can correct the fluxes using the pressure field and this way ensure div-free condition. This is the famous phi -= pEqn.flux();
6) Finally, we correct the velocity field, acquiring a good approximation of the velocity field.

Note that the solution are the three fields: U, p and phi. And the imporatant one is in fact phi not U.

Hope I managed to clear the things a bit.

Best,
Pawel
psosnows is offline   Reply With Quote

Reply

Thread Tools
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
read scalar field phi, get flux through faces peterwy OpenFOAM Programming & Development 4 August 4, 2010 16:36
Turbulence Model phi vs phi_ doug OpenFOAM Running, Solving & CFD 4 November 10, 2009 05:33
Another phi question ehsan_vaghefi OpenFOAM Running, Solving & CFD 0 October 24, 2008 19:56
About phi in icoFoam kar OpenFOAM Running, Solving & CFD 3 February 20, 2008 06:20
PHI file structure Eugene Phoenics 9 November 2, 2001 23:00


All times are GMT -4. The time now is 12:14.