
[Sponsors] 
phi = pEqn.flux() vs. linearInterpolate(U) & mesh.Sf() 

LinkBack  Thread Tools  Search this Thread  Display Modes 
May 26, 2013, 11:42 

#21 
Member
Emad Tandis
Join Date: Sep 2010
Posts: 74
Rep Power: 11 
thanks
your points are very helpful, specially "U is not divergence free". but there is thing that I can not understand about openFoam operators: are these relations correct? a) fvc::laplacian(rUA,p)=fvc::div(rUA,fvc::grad(p)) ?? b) does pEqn.flux() return (fvc::interpolate(rUA*fvc::grad(p)) & mesh.Sf ( based on jasak' thesis) from my experience in openFOAM a) is not true but I know from CFD this should be true! about b): my experience from openFOAM says me this is not be true but from 3.144 jasak it should be true. 

May 26, 2013, 12:50 

#22 
Senior Member
Pawel Sosnowski
Join Date: Mar 2009
Location: Munich, Germany
Posts: 105
Rep Power: 13 
Regarding point a)
take a look into Programmers Guide. If you would like to expand fvc::laplacian, you would get someting like this: fvc::laplacian(A,p) = fvc::surfaceSum ( fvc::interpolate(A) * fvc::snGrad(p) ); but consider that in the code we have fvm, which is implicit and returns fvMatrix object. And also there is a difference between fvc::grad() and fvc::snGrad. The latter one returns surfaceTypeField, and the first is in fact an average of snGrad in the cell center (well, depending on the scheme). Regarding b) I think that this should be true. But I will not put my head on the line right now. I remember checking it in the past, but can not find my notes at the moment. I would strongly advice you to dig into the code and check the structure yourself. If you run into too much trouble, give a call Best, Pawel 

May 27, 2013, 03:07 

#23 
Member
Emad Tandis
Join Date: Sep 2010
Posts: 74
Rep Power: 11 
Actually,
fvc::laplacian(A,p)=fvc::surfaceIntegrate(fvc::int erpolate(A)*fvc::snGrad(p)*mag(mesh.Sf())) but I mean we can take from definition of laplacian in mathematics: laplacian(A,p)=napla.(A*napla(p)) therefore, in OpenFOAM these sentence should have same results (with some numerical error ): fvc::laplacian(A,p) fvc::div(A*fvc::grad(p)) but their results, based on my experiece in OF, are quite different. Why? 

May 27, 2013, 03:30 

#24 
Senior Member
Pawel Sosnowski
Join Date: Mar 2009
Location: Munich, Germany
Posts: 105
Rep Power: 13 
Hello,
yup, I missed the Sf() for some reason I thought it was hidden in snGrad(). Regarding the semantics of laplacian always take causion with which kind of field are you dealing with, and what is the input for the functions. The direct mathematical semantics makes sense, and gives the result which is another way of calculating the operator. At the same time, it performs few more surfacetovolume operations, thus smears the result. The corrected definition you gave is a shorter (and more precise) version of discrete representation of laplacian. What is does, is calculating the surface values of the operators inside divergence, and then using those surface fields (and Gauss theorem) it calculates the divergence at the cell center. In case of direct mathematical representation: First you calculate the gradient at the face centers, * then you average it to get the values in the cell centers, * then you interpolate the gradient again on the face centers, ** then you interpolate A to the surfaces, and finally you calculate the divergence in the CC's In case of OF representation, you drop the (*) points and instead of (**) you calc the gradient on the surface directly. Hope it clarified a bit the issue. Best, Paweł 

May 27, 2013, 04:10 

#25 
Member
Emad Tandis
Join Date: Sep 2010
Posts: 74
Rep Power: 11 
you are quite right !
Based on this points, this is numerical Error. regarding pEqn.flux() it returns: fvc::interpolate(rUA)*fvc::snGrad(p)*mag(mesh.Sf() ) NOT fvc::interpolate(rUA)*fvc::interpolate(fvc::Grad(p )) & mesh.Sf() that I said before. this difference is based on your points. Thanks Very Much 

May 30, 2013, 18:41 

#26 
Member
Thomas Boucheres
Join Date: May 2013
Posts: 41
Rep Power: 8 
Hi,
ok there has been many interesting questions/answers. I just give complements because I think the basic knowledge of different peoples is very different in this post and maybe it is usefull to recall the cornerstone of the problem and its basis. OF uses collocated scheme for solving NS equation. Means velocity and pressure are both stored in cell centers. Since incompressible NS is not a classical evolution equation (because the presence of the constraint div(U) which in turns implies that one don't have equation for pressure), this collocated approach is numerically instable in the finite volume framework. Long to explain why but the interested reader should search and study the following keymords: collocated scheme versus staggered scheme pressure checkboard instability Rhie Chow stabilisation infsup condition for mixed problems in finite element context All the core reasons why one cannot use an interpolation scheme (whatever it is precise type, linear or ...) for phi field come from those notions. Not so difficult but if the reader has strictly no knowledge on numerical analysis, it will be to hard to understand in deep. Good luck. 

October 18, 2013, 12:43 

#27  
Senior Member
Vesselin Krastev
Join Date: Jan 2010
Location: University of Tor Vergata, Rome
Posts: 368
Rep Power: 14 
Quote:
after reading your post and most of the previous ones I think that now pressurevelocity coupling in OpenFOAM is getting much clearer to me, so thank you all (also Alberto, Santiago and the others) for the interesting explanations. I have just one little doubt left about the momentum predictor stage: you said above that it is, in principle, unnecessary for the successful application of the correction procedure; however, both in H.Jasak's and E. De Villiers Theses, I found passages claiming that the predicted velocities are used to update the H(U) vector before calculating the "pseudovelocities" U=rAU*H(U) and interpolating them for the nonnull divergence face fluxes evaluation (the points 4.1 and 4.2 of your post). These affirmations are apparently supported by the fact that in OpenFOAM H(U) appears to be calculated right after the momentum predictor is solved (see the pisoFoam.C source below). Code:
if (momentumPredictor) { solve(UEqn == fvc::grad(p)); } //  PISO loop for (int corr=0; corr<nCorr; corr++) { volScalarField rAU(1.0/UEqn.A()); volVectorField HbyA("HbyA", U); HbyA = rAU*UEqn.H(); surfaceScalarField phiHbyA ( "phiHbyA", (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, U, phi) ); So, my questions are: 1) how and where is actually used the velocity field obtained from the momentum predictor? Is it, for instance, only used for the neighbour velocities update inside H(U), whereas face fluxes phi, which are "hidden" in the convective discretization coefficients, are left "frozen" until the pressure equation is solved? 2) depending on the answer to question 1), is it really negligible the predictor stage influence on the correction procedure (or, alternatively, am I missing something in my considerations)? Can you please comment on this? Of course comments from other sources are also appreciated Thanks in advance V. 

October 18, 2013, 16:29 

#28 
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 445
Rep Power: 19 
Hi Vesselin, each time you do H(U) the last calculated U is used. In the first PISO loop is the velocity obtained in the momentum predictor, in case if it is not performed is the velocity of the previous timestep. From the second corrector and later the last reconstructed velocity is used. In this case the coefficients of H() are frozen.
If you use pimple solvers you can do outer loops which will renew the H() operator since a new UEqn is assembled in each outer loop Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Research Scientist Research Center for Computational Methods (CIMEC)  CONICET/UNL Tel: 543424511594 Int. 7032 Colectora Ruta Nac. 168 / Paraje El Pozo (3000) Santa Fe  Argentina. http://www.cimec.org.ar 

October 19, 2013, 05:07 

#29  
Senior Member
Vesselin Krastev
Join Date: Jan 2010
Location: University of Tor Vergata, Rome
Posts: 368
Rep Power: 14 
Hi Santiago,
and thanks a lot for your fast reply. However, just for a better comprehension, I would like to further address some the points of your response: Ok, but is it the same for phi (which is included inside convection discretization coefficients, that contribute to build H)? If it is so, in the first PISO loop the phi values (and thus the discretization coefficients) should be the ones created by "createPhi.H" or coming from the previous time step/correction loop, and should remain the same from the momentum predictor (if present) to the pressure equation solution, after which they are corrected enforcing mass conservation. Thus, the second PISO loop will start with H(U) updated with the last corrected cell center velocity field (explicit correction after the pressure equation solution), but also with the last (at this stage, conservative) face fluxes phi. And so on for the next (eventual) PISO sequences. Is it right? Quote:
Thanks once again V. 

October 19, 2013, 07:57 

#30  
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 445
Rep Power: 19 
Quote:
Quote:
In the case of steady solvers like simpleFoam the trick to avoid raw approximations at the begginning is to low the relaxation factors. Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Research Scientist Research Center for Computational Methods (CIMEC)  CONICET/UNL Tel: 543424511594 Int. 7032 Colectora Ruta Nac. 168 / Paraje El Pozo (3000) Santa Fe  Argentina. http://www.cimec.org.ar 

October 19, 2013, 12:25 

#31  
Senior Member
Vesselin Krastev
Join Date: Jan 2010
Location: University of Tor Vergata, Rome
Posts: 368
Rep Power: 14 
Quote:
So, to have a completely renewed H(U), one has to move to the next time step (or outer iteration, if a PIMPLE solver is used), at the beginning of which UEqn is ALWAYS assembled (but solved as a momentum predictor only if the option is on). Ok, I think that now things are all in place, sorry if I made it longer than strictly necessary and thanks for your patience (and of course please correct me if you find some other inconsistencies). Regards V. 

April 9, 2014, 08:28 

#32 
Member
Fabian E.
Join Date: Nov 2009
Posts: 38
Rep Power: 11 
Hello,
it is hard for me to follow your interesting discussion, so I wanted to have a practical approach on the phi variable. What I thought that right after calling phi = pEqn.flux(); mass conversion should be valid. I tried to check this for a simpleFoam and a rhoPimpleFoam tutorial case. Since phi is given as (kg s1) for all faces, the sum of phi for all faces of one cell should be zero. I took into account the normal vecs orientation: if (mesh.faceOwner()[mesh.cells()[cell_label][face_label]]) then the normal vec is pointing inside the cell. Else I multipy phi with 1, since phi is pointing towards the neighbouring cell (same in case of boundary face). I never see a sum 0, but rather a sum of ca. 1020% of the sum of the absolute values. So how do I calculate the correct mass flow in and out of one arbitary cell? 

June 12, 2014, 01:50 

#33 
Member
Fabian E.
Join Date: Nov 2009
Posts: 38
Rep Power: 11 
I finally solved it, the face orientation was the problem.
May this helps anyone in the future: Code:
volVectorField phi_V ( IOobject ( "phi_V", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), (rho*U) ); surfaceVectorField phi_V_interpolated ( IOobject ( "phi_V_interpolated", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), linearInterpolate(phi_V) ); forAll(U, cell_label) { double flux_sum_neg = 0; double flux_sum_pos = 0; forAll(mesh.cells()[cell_label], face_label) { double sign = 0; double sign_phi = 0; if (!mesh.isInternalFace(mesh.cells()[cell_label][face_label])) { sign = 1; label patch_index = mesh.boundaryMesh().whichPatch(mesh.cells()[cell_label][face_label]); label face_index = mesh.cells()[cell_label][face_label]  mesh.boundaryMesh()[patch_index].start(); vector normal = sign * mesh.Sf().boundaryField()[patch_index][face_index] / mesh.magSf().boundaryField()[patch_index][face_index] ; sign_phi = (((phi_V_interpolated.boundaryField()[patch_index][face_index] & normal) < 0) ? 1 : 1); if (sign_phi < 0) flux_sum_neg += std::abs(phi.boundaryField()[patch_index][face_index]); if (sign_phi > 0) flux_sum_pos += std::abs(phi.boundaryField()[patch_index][face_index]); } else { sign = ((mesh.faceNeighbour()[mesh.cells()[cell_label][face_label]] == cell_label) ? 1 : 1); vector normal = sign * mesh.Sf()[mesh.cells()[cell_label][face_label]] / mesh.magSf()[mesh.cells()[cell_label][face_label]]; sign_phi = (((phi_V_interpolated[mesh.cells()[cell_label][face_label]] & normal) < 0) ? 1 : 1); if (sign_phi < 0) flux_sum_neg += std::abs(phi[mesh.cells()[cell_label][face_label]]); if (sign_phi > 0) flux_sum_pos += std::abs(phi[mesh.cells()[cell_label][face_label]]); } } double error = (flux_sum_pos  flux_sum_neg) / flux_sum_pos * 100; //in % if (std::abs(error) > 1) { .. } else { Info << "massflow into " << cell_label << " = " << flux_sum_pos * mesh.cellVolumes()[cell_label] << " kg/s" << endl; } } 

Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
read scalar field phi, get flux through faces  peterwy  OpenFOAM Programming & Development  5  April 28, 2017 19:40 
Turbulence Model phi vs phi_  doug  OpenFOAM Running, Solving & CFD  4  November 10, 2009 04: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 05:20 
PHI file structure  Eugene  Phoenics  9  November 2, 2001 22:00 