
[Sponsors] 
February 14, 2011, 16:17 
Temperature anomoly at pressure reference cell

#1 
New Member
Will Logie
Join Date: Sep 2010
Location: ANU, Canberra, Australia
Posts: 19
Rep Power: 7 
Hello Foamers,
I've been witnessing temperature divergence at the pressure reference cell in a buoyantBoussinesqPimpleFoam problem I've defined axisymmetrically; a cylindrical hot water tank with a helixcoiled heat exchanger became a 2° wedge of 1 cell thickness. The picture attached hopefully gives an idea of this. In it you can also see where I placed the pressure reference and how the temperature (shown here at 25 seconds) at this point is rising (inexplicably for my liking) and causing all sorts of mischief. Here my fvSolution: Code:
/** C++ **\  =========    \\ / F ield  OpenFOAM: The Open Source CFD Toolbox   \\ / O peration  Version: 1.7.0   \\ / A nd  Web: www.OpenFOAM.com   \\/ M anipulation   \**/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object fvSolution; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // solvers { p_rgh { solver PCG; preconditioner DIC; tolerance 1e8; relTol 0.01; } p_rghFinal { $p_rgh; relTol 0; } "(UTkepsilonomegaR)" { solver PBiCG; preconditioner DILU; tolerance 1e6; relTol 0.1; } "(UTkepsilonomegaR)Final" { $U; relTol 0; } } PIMPLE { momentumPredictor yes; nOuterCorrectors 1; nCorrectors 2; nNonOrthogonalCorrectors 0; //pRefPoint ( 0.25 0.13 0 ); pRefCell 684; pRefValue 0; } relaxationFactors { "(UTkepsilonomegaR)" 1; "(UTkepsilonomegaR)Final" 1; } // ************************************************************************* // Code:
/** C++ **\  =========    \\ / F ield  OpenFOAM: The Open Source CFD Toolbox   \\ / O peration  Version: 1.7.0   \\ / A nd  Web: www.OpenFOAM.com   \\/ M anipulation   \**/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object fvSchemes; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // ddtSchemes { default Euler; } gradSchemes { default Gauss linear; } divSchemes { default none; div(phi,U) Gauss upwind; div(phi,T) Gauss upwind; div(phi,k) Gauss upwind; div(phi,epsilon) Gauss upwind; div(phi,omega) Gauss upwind; div(phi,R) Gauss upwind; div(R) Gauss linear; div((nuEff*dev(grad(U).T()))) Gauss linear; } laplacianSchemes { default none; laplacian(nuEff,U) Gauss linear uncorrected; laplacian((1A(U)),p_rgh) Gauss linear uncorrected; laplacian(kappaEff,T) Gauss linear uncorrected; laplacian(DkEff,k) Gauss linear uncorrected; laplacian(DepsilonEff,epsilon) Gauss linear uncorrected; laplacian(DomegaEff,omega) Gauss linear uncorrected; laplacian(DREff,R) Gauss linear uncorrected; } interpolationSchemes { default linear; } snGradSchemes { default uncorrected; } fluxRequired { default no; p_rgh ; } // ************************************************************************* // Thanks, Will. 

February 16, 2011, 07:31 
Anomaly in the sense that...

#2 
New Member
Will Logie
Join Date: Sep 2010
Location: ANU, Canberra, Australia
Posts: 19
Rep Power: 7 
... the temperature at pRefCell oscillates in time from being the lowest (~265K) to the highest (~310K).
In previous cases where I had not yet included the axisymmetry of the problem (simple 2D extruded mesh) I was not witnessing this in any way. I am thus inclined to think that my error lies therein but am ill equipped to diagnose further. I'd really appreciate some thoughts from you experts out there! 

February 16, 2011, 09:33 
checkMesh failed?

#3 
New Member
Will Logie
Join Date: Sep 2010
Location: ANU, Canberra, Australia
Posts: 19
Rep Power: 7 
The only thing I can see that seems remotely suspicious is a point in the front wedge slipping outside a threshold (1e8 meter?). I created the mesh in what I thought was a rather cleanly parameterised gmsh file:
Code:
/**\  =========    \\ / F ield  OpenFOAM: The Open Source CFD Toolbox   \\ / O peration  Version: 1.7.x   \\ / A nd  Web: www.OpenFOAM.com   \\/ M anipulation   \**/ Build : 1.7.xda2e560009e6 Exec : checkMesh Date : Feb 16 2011 Time : 14:25:58 Host : pcspf163 PID : 14620 Case : /home/will/OpenFOAM/will1.7.x/run/ihx/buoyantBoussinesqPimpleFoam/ihx35mmCourseWedge nProcs : 1 SigFpe : Enabling floating point exception trapping (FOAM_SIGFPE). // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Create time Create polyMesh for time = 0 Time = 0 Mesh stats points: 11948 internal points: 0 faces: 23470 internal faces: 11521 cells: 5846 boundary patches: 5 point zones: 0 face zones: 0 cell zones: 1 Overall number of cells of each type: hexahedra: 5761 prisms: 85 wedges: 0 pyramids: 0 tet wedges: 0 tetrahedra: 0 polyhedra: 0 Checking topology... Boundary definition OK. Point usage OK. Upper triangular ordering OK. Face vertices OK. Number of regions: 1 (OK). Checking patch topology for multiply connected surfaces ... Patch Faces Points Surface topology tesWall 209 418 ok (nonclosed singly connected) front 5846 6017 ok (nonclosed singly connected) pipeWall 48 96 ok (nonclosed singly connected) back 5846 6017 ok (nonclosed singly connected) defaultFaces 0 0 ok (empty) Checking geometry... Overall domain bounding box (0 0 0.00567203) (0.324951 0.3 0.00567203) Mesh (nonempty, nonwedge) directions (1 1 0) Mesh (nonempty) directions (1 1 1) Wedge front with angle 1.00004 degrees ***Wedge patch front not planar. Point (0.260332 0.0457038 0.00454412) is not in patch plane by 3.73728e08 meter. Boundary openness (2.60064e19 6.7127e20 6.3524e19) OK. Max cell openness = 3.04681e16 OK. Max aspect ratio = 9.45051 OK. Minumum face area = 1.83948e07. Maximum face area = 0.000137058. Face area magnitudes OK. Min volume = 1.54238e09. Max volume = 3.80042e07. Total volume = 0.000550207. Cell volumes OK. Mesh nonorthogonality Max: 41.1657 average: 8.53724 Nonorthogonality check OK. Face pyramids OK. Max skewness = 0.952198 OK. Failed 1 mesh checks. End 

February 18, 2011, 11:51 
... a bigger example

#4 
New Member
Will Logie
Join Date: Sep 2010
Location: ANU, Canberra, Australia
Posts: 19
Rep Power: 7 
... after 200 seconds.


April 5, 2011, 19:52 

#5 
Senior Member
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 12 
Well your checkMesh says the front plane of the wedge is nonplanar. If you still have the problem try replacing the wedge boundaries with symmetryPlane  if the problem goes away, then the issue is your mesh not being a perfect enough wedge.


April 6, 2011, 04:42 
ascii > binary

#6 
New Member
Will Logie
Join Date: Sep 2010
Location: ANU, Canberra, Australia
Posts: 19
Rep Power: 7 
If one looks inside the average controlDict file there are these two settings which I never really paid attention to before; writeFormat and writePrecision.
Changing writeFormat to binary worked a charm... (thank you Mattijs!) Checking geometry... Overall domain bounding box (0 0 0.00567203) (0.324951 0.5 0.00567203) Mesh (nonempty, nonwedge) directions (1 1 0) Mesh (nonempty) directions (1 1 1) Wedge front with angle 1 degrees Wedge back with angle 1 degrees All edges aligned with or perpendicular to nonempty directions. Boundary openness (1.06916e19 5.56146e20 8.77993e19) OK. Max cell openness = 3.05762e16 OK. Max aspect ratio = 4.96254 OK. Minumum face area = 1.22973e07. Maximum face area = 3.30134e05. Face area magnitudes OK. Min volume = 6.95689e10. Max volume = 7.61277e08. Total volume = 0.000910615. Cell volumes OK. Mesh nonorthogonality Max: 42.3953 average: 9.13043 Nonorthogonality check OK. Face pyramids OK. Max skewness = 0.529742 OK. Mesh OK. 

April 13, 2011, 13:06 

#7 
Member

After run makeAxialMesh and collapseEdges.
I found a similar problem. however ascii to binary doesn't work for me. My grid looks fine and the checkMesh before the makeAxialMesh and collapseEdges reports no errors. Any suggestions? Thanks Carlos Code:
Create time Create polyMesh for time = 0 Time = 0 Mesh stats points: 2900 internal points: 0 faces: 5558 internal faces: 2659 cells: 1386 boundary patches: 7 point zones: 0 face zones: 0 cell zones: 0 Overall number of cells of each type: hexahedra: 1287 prisms: 99 wedges: 0 pyramids: 0 tet wedges: 0 tetrahedra: 0 polyhedra: 0 Checking topology... Boundary definition OK. Point usage OK. Upper triangular ordering OK. Face vertices OK. Number of regions: 1 (OK). Checking patch topology for multiply connected surfaces ... Patch Faces Points Surface topology frontback 0 0 ok (empty) inlet 14 29 ok (nonclosed singly connected) outlet 14 29 ok (nonclosed singly connected) simetry 0 0 ok (empty) wall 99 200 ok (nonclosed singly connected) frontback_pos 1386 1500 ok (nonclosed singly connected) frontback_neg 1386 1500 ok (nonclosed singly connected) Checking geometry... Overall domain bounding box (0 1 0.360028) (10 3.035 1.36003) Mesh (nonempty, nonwedge) directions (1 1 0) Mesh (nonempty) directions (1 1 1) Wedge frontback_pos with angle 22.91 degrees ***Wedge patch frontback_pos not planar. Point (0.10101 1 0.5) is not in patch plane by 1.28966e06 meter. Boundary openness (5.88638e18 1.0902e15 5.47392e16) OK. Max cell openness = 2.36954e16 OK. Max aspect ratio = 2.83371 OK. Minumum face area = 0.00215678. Maximum face area = 0.241098. Face area magnitudes OK. Min volume = 0.000218165. Max volume = 0.0227451. Total volume = 8.04848. Cell volumes OK. Mesh nonorthogonality Max: 21.5012 average: 8.88707 Nonorthogonality check OK. Face pyramids OK. Max skewness = 1.05346 OK. Failed 1 mesh checks. End 

April 14, 2011, 04:55 

#8  
New Member
Will Logie
Join Date: Sep 2010
Location: ANU, Canberra, Australia
Posts: 19
Rep Power: 7 
Quote:
23° is a pretty big wedge too. Maybe try <5°...? Which settings did you use for collapseEdges? 

April 14, 2011, 05:03 

#9 
Member

Hello Will and thank you for your reply.
I have already tried with 5 degrees and I get the same result. For the collapseEdges I used 1e8 and 180 degrees Code:
Create polyMesh for time = 0 Time = 0 Mesh stats points: 2900 internal points: 0 faces: 5558 internal faces: 2659 cells: 1386 boundary patches: 7 point zones: 0 face zones: 0 cell zones: 0 Overall number of cells of each type: hexahedra: 1287 prisms: 99 wedges: 0 pyramids: 0 tet wedges: 0 tetrahedra: 0 polyhedra: 0 Checking topology... Boundary definition OK. Point usage OK. Upper triangular ordering OK. Face vertices OK. Number of regions: 1 (OK). Checking patch topology for multiply connected surfaces ... Patch Faces Points Surface topology frontback 0 0 ok (empty) inlet 14 29 ok (nonclosed singly connected) outlet 14 29 ok (nonclosed singly connected) simetry 0 0 ok (empty) wall 99 200 ok (nonclosed singly connected) frontback_pos 1386 1500 ok (nonclosed singly connected) frontback_neg 1386 1500 ok (nonclosed singly connected) Checking geometry... Overall domain bounding box (0 1 0.411235) (10 3.035 0.588765) Mesh (nonempty, nonwedge) directions (1 1 0) Mesh (nonempty) directions (1 1 1) Wedge frontback_pos with angle 2.49762 degrees ***Wedge patch frontback_pos not planar. Point (0.10101 1 0.5) is not in patch plane by 5.03641e07 meter. Boundary openness (6.12585e19 3.61958e17 8.53818e16) OK. Max cell openness = 2.13747e16 OK. Max aspect ratio = 2.83371 OK. Minumum face area = 0.000222607. Maximum face area = 0.0248842. Face area magnitudes OK. Min volume = 2.25174e05. Max volume = 0.00234757. Total volume = 0.830702. Cell volumes OK. Mesh nonorthogonality Max: 21.5012 average: 8.88707 Nonorthogonality check OK. Face pyramids OK. Max skewness = 0.516369 OK. Failed 1 mesh checks. 

April 14, 2011, 05:30 
Hmmm...

#10  
New Member
Will Logie
Join Date: Sep 2010
Location: ANU, Canberra, Australia
Posts: 19
Rep Power: 7 
Hey Carlos,
I've just had a quick look at my makeAxialMesh example and found that this problem you mention persists for me too when I use makeAxialMesh, regardless of binary and writePrecision. The example where ascii > binary worked for me (post from April 6th) does not use makeAxialMesh. There I create the wedge in gmsh (translate mesh 1° about symmetryPlane and extrude/revolve 2° about symmetryPlane) and convert that to foam. So I'm not sure I can help you much further with makeAxialMesh  are you using gmsh? I'm also not sure if you need to straddle the xyplane equally or not. Maybe someone else can comment on the interpretation from the user manual: Quote:


April 14, 2011, 05:34 

#11 
Member

I'm using Pointwise for the grid.
Thanks for your help anyway. What is curious is that the test case that came with makeAxial reports the same error. 

April 2, 2012, 16:39 

#12 
New Member
Join Date: Jan 2010
Posts: 23
Rep Power: 7 
I'm seeing the same issue, except the mesh isn't axisymmetric. The mesh is a cube of hexahedral cells, looking at a buoyancy flows (solving for piezometric pressure "pd" instead of "p = pd + rhogh"). "pd" is set to "zeroGradient" on all boundaries, and the reference value for pdEqn is supplied by calculating the total pressure:
if (pd.needReference()) { p += dimensionedScalar ( "p", p.dimensions(), pRefValue  getRefCellValue(p, pRefCell) ); pd = p  rhokamb*gh; } Everything runs fine when you specify the temperature of the top/bottom boundaries (e.g. "fixedValue")  but if you change one of the boundaries (say to "zeroGradient"), the temperature in the cell that corresponds to "pRefCell" increases, like will.logie showed. I've tried moving the reference cell around, but wherever you move it, that's where you see the temperature anomaly. It's strange because one of the temperature boundaries is still fixed, so a check to something like: if (T.needReference()) { Info << "Applying reference temperature " << getRefCellValue(T,pRefCell) << " at point " << mesh.C()[pRefCell] << endl; } doesn't spit out anything, so I can only assume that the TEqn doesn't need a reference cell. Of course, if I set all the temperature boundaries to "zeroGradient" this loop spits out info, and you can see the temperature in that cell getting changed... How/where could "pRefCell" be entering into the TEqn? Any thoughts? Thanks. 

April 5, 2012, 16:36 

#13 
New Member
Join Date: Jan 2010
Posts: 23
Rep Power: 7 
To try to offer up some more information...
The error comes during the TEqn.solve() step. If you look at the temperature at the cell "pRefCell" you can see it changing wildly after the solution. Changed the top boundary type on T to "fixedGradient" (instead of trying "zeroGradient"), setting the gradient equal to the temperature gradient for the stratification (linear stratification). Bottom boundary still constant value. This appears to fix the problem running in serial, but in parallel the error occurs  the temperature in the cell that corresponds to "pRefCell" becomes unphysical. Attached is a snapshot of what the parallel solution looks like after about 50 iterations. Any thoughts on how pRefCell could be interfering with the TEqn, even though one of the temperature equation boundaries is a "fixedValue"? This happens in both OF2.1.0 and OF1.6ext (tried both). Tried different solvers for T as well (PCBiCG and GAMG), same deal. Also, when I use a "fixedGradient" at the top, why would it work in serial but not in parallel? 

April 10, 2012, 08:50 

#14 
New Member
Join Date: Jan 2010
Posts: 23
Rep Power: 7 
Bump. To try and provide some more info, wrote out the location of the pressure reference cell, and then applying:
TEqn.setReference(pRefCell,T.internalField()[pRefCell]); to see if it's actually setting a temperature reference. Here it would basically set the reference temperature to whatever it was before the TEqn.solve(). Here is an example output  notice how T[pRefCell] is increasing and not equal to what it was before the solve.... Code:
Time = 0.01 Courant Number mean: 0.00783788 max: 0.346089 DILUPBiCG: Solving for Ux, Initial residual = 0.00193302, Final residual = 4.01574e10, No Iterations 3 DILUPBiCG: Solving for Uy, Initial residual = 0.00193302, Final residual = 4.01574e10, No Iterations 3 DILUPBiCG: Solving for Uz, Initial residual = 0.00298749, Final residual = 7.85053e10, No Iterations 3 DICPCG: Solving for pd, Initial residual = 1, Final residual = 9.99605e09, No Iterations 313 time step continuity errors : sum local = 2.48078e09, global = 2.47807e09, cumulative = 2.47807e09 DICPCG: Solving for pd, Initial residual = 0.00911823, Final residual = 9.80415e09, No Iterations 283 time step continuity errors : sum local = 2.48337e09, global = 2.47807e09, cumulative = 4.95615e09 pRefCell location = (7.5 8.9 19.3) Setting reference temperature at (7.5 8.9 19.3) to 299.397 T[pRefCell] = 299.397 DILUPBiCG: Solving for T, Initial residual = 8.09426e06, Final residual = 8.37418e09, No Iterations 1 T[pRefCell] = 299.768 ExecutionTime = 50.21 s ClockTime = 51 s Time = 0.02 Courant Number mean: 0.000653823 max: 0.393349 DILUPBiCG: Solving for Ux, Initial residual = 0.952681, Final residual = 1.70106e10, No Iterations 3 DILUPBiCG: Solving for Uy, Initial residual = 0.952681, Final residual = 1.70104e10, No Iterations 3 DILUPBiCG: Solving for Uz, Initial residual = 0.92743, Final residual = 3.8282e10, No Iterations 3 DICPCG: Solving for pd, Initial residual = 0.955734, Final residual = 9.13164e09, No Iterations 307 time step continuity errors : sum local = 1.65392e09, global = 1.65143e09, cumulative = 6.60758e09 DICPCG: Solving for pd, Initial residual = 0.0156603, Final residual = 9.83983e09, No Iterations 288 time step continuity errors : sum local = 1.65325e09, global = 1.65143e09, cumulative = 8.25901e09 pRefCell location = (7.5 8.9 19.3) Setting reference temperature at (7.5 8.9 19.3) to 299.768 T[pRefCell] = 299.768 DILUPBiCG: Solving for T, Initial residual = 7.12715e06, Final residual = 2.1977e09, No Iterations 1 T[pRefCell] = 299.995 ExecutionTime = 85.37 s ClockTime = 86 s 

August 29, 2012, 09:39 

#15 
New Member
Join Date: Jan 2010
Posts: 23
Rep Power: 7 
So I think I found a solution. The strange temperature behavior at pRefCell was not coming from the temperature equation, but from the pressure equation. During the pressure solution, the value of the density perturbation (rhok) at pRefCell would fluctuate. I'm not 100% certain yet, but it's related to the fact that the flux needs to be adjusted AFTER the inclusion of buoyancy. That is, in pEqn.H:
Code:
surfaceScalarField phiU ( (fvc::interpolate(U) & mesh.Sf()) ); //adjustPhi(phiU, U, pd); phi = phiU + (rAUf * ((fvc::interpolate(rhok)) * (g & mesh.Sf()))); adjustPhi(phi, U, pd); 

June 5, 2013, 18:55 
pRefCell on Parallel Simuations

#16 
New Member
Join Date: Apr 2012
Posts: 21
Rep Power: 5 
Hi
I have developed a solver which works fine on single processor, but diverges when I am running it on multiple processors. For the later, by changing the cell at which pressure is set (pRefCell), my results change. i.e. pRefCell is affecting the results on parallel simulations! Any hints? 

August 30, 2013, 16:58 

#17 
Senior Member

Dear all,
For those who have the problem about wedge patch which was mentioned above, I want to inform that I solved the problem using Bernhard Gschaider's post in here: How to make a pipe using makeAxialMesh
Mojtaba
__________________
Complex Heat & Flow Simulation Research Group If you can't explain it simply, you don't understand it well enough. "Richard Feynman" 

September 1, 2013, 11:21 

#18 
Senior Member
Ehsan
Join Date: Oct 2012
Location: Iran
Posts: 2,205
Rep Power: 18 
Hi Mojtaba
how writePrecision can help?
__________________
Injustice Anywhere is a Threat for Justice Everywhere.Martin Luther King. To Be or Not To Be,Thats the Question! The Only Stupid Question Is the One that Goes Unasked. 

September 1, 2013, 14:29 

#19 
Senior Member

Hi Ehsan,
Well more values of writePrecision will make OF to write values with more accuracy. In wedge problems, having a symmetry geometry is one of the main concerns, therefore we are trying to obtain more accuracy in order to have a much better symmetry in our geometry. Changing writePrecision will help this.
__________________
Complex Heat & Flow Simulation Research Group If you can't explain it simply, you don't understand it well enough. "Richard Feynman" 

March 19, 2014, 09:11 

#20  
Senior Member

Quote:
Code:
fvScalarMatrix p_rghEqn ( fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA) ); p_rgh = prhok*gh; p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
__________________
Training Course on OpenFOAM at (http://www.isme.ir/) My Weblog (http://openfoam.blogfa.com/) 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Installing OpenFOAM1.5dev on a cluster  ZKM  OpenFOAM Installation  4  December 25, 2010 16:59 
Neumann pressure BC and velocity field  Antech  Main CFD Forum  0  April 25, 2006 02:15 
what the result is negatif pressure at inlet  chong chee nan  FLUENT  0  December 29, 2001 06:13 
Hydrostatic pressure in 2phase flow modeling (CFX4.2)  HB &DS  CFX  0  January 9, 2000 14:19 
Hydrostatic pressure in 2phase flow modeling (long)  DS & HB  Main CFD Forum  0  January 8, 2000 16:00 