# chtMultiRegionSimpleFoam fails to solve cases without gravity

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

January 4, 2016, 12:05
chtMultiRegionSimpleFoam fails to solve cases without gravity
#1
New Member

Karl Lindqvist
Join Date: Jul 2012
Posts: 20
Rep Power: 7
Dear all,
I am facing a rather peculiar problem with the chtMultiRegionSimpleFoam solver. I am interested in solving heat transfer from a solid to a fluid with a temperature dependent density, but would like to neglect buoyancy effects (for now). It seems, however, that the chtMultiRegionSimpleFoam solver is unable to solve the p_rgh equation correctly for cases where g = (0 0 0) and the (forced) velocity is low.

The attached case (adapted from the planeWall2D tutorial case at OpenFOAMWiki) demonstrates the issue in a simplistic case of forced convection over a heated floor (see figure below). A bottomAir region will also be present in the solution, but can be ignored (it is not solved for).

The conjugate heat transfer causes a temperature gradient through the floor, and a thermal boundary layer in the fluid. The case solves as expected with the attached parameters where the inlet velocity is 10 m/s, but changing this value to 1 m/s (in the corresponding changeDictionaryDict) causes the solver iterations for the p_rgh equation to reach its maximum limit after a short time (here at iteration 389):

Code:
```Time = 388

Solving for fluid region topAir
DILUPBiCG:  Solving for Ux, Initial residual = 0.0001642792187472528, Final residual = 2.086793044719327e-06, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 0.0003276496308532365, Final residual = 9.054436042068467e-06, No Iterations 1
DILUPBiCG:  Solving for h, Initial residual = 0.0001652516759873468, Final residual = 1.096169895436815e-05, No Iterations 1
Min/max T:299.9999249381647 485.7994890854322
GAMG:  Solving for p_rgh, Initial residual = 7.12858129068406e-05, Final residual = 7.03830111216478e-07, No Iterations 39
time step continuity errors : sum local = 4.644778702141491e-09, global = -3.757627749083547e-11, cumulative = 6.025277233104804e-08
Min/max rho:0.7154943606653349 1.158623823918645

Solving for solid region wall
DICPCG:  Solving for h, Initial residual = 0.006403279730421732, Final residual = 3.900945002001053e-05, No Iterations 1
Min/max T:min(T) [0 0 0 1 0 0 0] 483.9315996342299 max(T) [0 0 0 1 0 0 0] 500
ExecutionTime = 8.19 s  ClockTime = 8 s

Time = 389

Solving for fluid region topAir
DILUPBiCG:  Solving for Ux, Initial residual = 0.0001626781773516434, Final residual = 2.069908791677873e-06, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 0.0003204160935846439, Final residual = 9.270251346231881e-06, No Iterations 1
DILUPBiCG:  Solving for h, Initial residual = 0.0001638593213157633, Final residual = 1.087279304981332e-05, No Iterations 1
Min/max T:299.9999248933859 485.9034691736506
GAMG:  Solving for p_rgh, Initial residual = 6.801409757909533e-05, Final residual = 7.818732180590184e-07, No Iterations 1000
time step continuity errors : sum local = 4.87402540356367e-09, global = -3.229735324601451e-11, cumulative = 6.022047497780203e-08
Min/max rho:0.7153412496577295 1.158623824019448

Solving for solid region wall
DICPCG:  Solving for h, Initial residual = 0.006401407648684744, Final residual = 3.900699891191313e-05, No Iterations 1
Min/max T:min(T) [0 0 0 1 0 0 0] 484.0462054555338 max(T) [0 0 0 1 0 0 0] 500
ExecutionTime = 8.57 s  ClockTime = 9 s

Time = 390

Solving for fluid region topAir
DILUPBiCG:  Solving for Ux, Initial residual = 0.0001611182905095051, Final residual = 2.053164802211731e-06, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 0.0003136746051429199, Final residual = 9.491366256280922e-06, No Iterations 1
DILUPBiCG:  Solving for h, Initial residual = 0.0001624826897776338, Final residual = 1.078451248755692e-05, No Iterations 1
Min/max T:299.9999248494713 486.0066718410475
GAMG:  Solving for p_rgh, Initial residual = 6.507513719961976e-05, Final residual = 7.436642579933397e-07, No Iterations 1000
time step continuity errors : sum local = 4.695772893242314e-09, global = -3.700505525107244e-11, cumulative = 6.018346992255096e-08
Min/max rho:0.7151893481877345 1.158623824119854

Solving for solid region wall
DICPCG:  Solving for h, Initial residual = 0.006399522667559255, Final residual = 3.900445715954452e-05, No Iterations 1
Min/max T:min(T) [0 0 0 1 0 0 0] 484.1599792231602 max(T) [0 0 0 1 0 0 0] 500
ExecutionTime = 9.07 s  ClockTime = 9 s```
I have not been able to find any physical or algorithmic reason for why the pressure field cannot be solved. Continued iterations results in a physically sound pressure field, but the residuals flatten out at around 1e-6 and the high iteration count causes a long solution time. Imposing a stricter limit on the max. no. of iterations "solves" the latter problem, but obviously not the root cause.

I have tried experimenting with different combinations of boundary conditions and initializations, without success. OpenFOAM 2.2, 2.4.x and 3.0.x all cause the same behavior, so does the transient solver version chtMultiRegionFoam. The mesh is block-structured and generated by blockMesh.

Comments and suggestions are highly appreciated. I also have some specific questions:
- Could this be related to the way the SIMPLE algorithm is constructed? Could the issue be solved by using e.g. the SIMPLER algorithm instead and avoid solving the pressure correction equation?
- My understanding of the OpenFOAM discretization is that p_rgh = p - rho*g*h_cell, such that setting g to zero should reduce the p_rgh equation to the p equation. Is this correct, or is the g*h term somehow still included in the equation system (through a momentum source term)?

Best regards,
Karl
Attached Images
 p_rgh_10.png (7.2 KB, 9 views)
Attached Files
 heatedFloor2D.tar.gz (5.9 KB, 9 views)

 January 4, 2016, 16:25 #2 Senior Member   Alex Join Date: Oct 2013 Posts: 336 Rep Power: 14 Hi Karl, Have you tried not to give the value of 0 to the gravity but a very low one different from 0? Regards, Alex __________________ Web site where I present my Master's Thesis: foamingtime.wordpress.com The case I talk about in this site was solved with chtMultiRegionSimpleFoam solver and involves radiation. Some basic tutorials are also resolved step by step in the web. If you are interested in these matters, you are invited to come in!

 January 5, 2016, 04:25 #3 New Member   Karl Lindqvist Join Date: Jul 2012 Posts: 20 Rep Power: 7 Alex, Thank you for your reply! Setting g to a small nonzero value seems to be a solution; U_in = 1 m/s requires g = (0 0.01 0) to solve and U_in = 0.1 m/s requires g = (0 0.1 0) to solve. Do you have any insight into why this is the case? I have, for instance, noticed that the case diverges if gravity is applied in the x-direction rather than the y-direction (irrespective of magnitude). Could this be related? Best Regards, Karl UPDATE: The proposed solution does not always work. I'll try to post an updated case (heated channel flow) within the next days for further discussion. Last edited by karlli; January 7, 2016 at 03:41.

January 8, 2016, 15:37
#4
New Member

Karl Lindqvist
Join Date: Jul 2012
Posts: 20
Rep Power: 7
As mentioned in the previous post, setting g to a small nonzero value is not a general solution. Consider, for example, the attached case where only the fluid region is considered but an unheated starting length has been added. All other parameters are the same as the case above, but the p_rgh equation breaks down after about 250 iterations reaching a residual of just above 1e-4.

Any suggestions?

Regards,
Karl
Attached Files
 partiallyHeatedFloor2D.tar.gz (141.8 KB, 4 views)

January 9, 2016, 15:01
#5
Super Moderator

Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 10,020
Blog Entries: 39
Rep Power: 109
Greetings to all!

I haven't taken the time to fully diagnose the problem, but I've found the main issue.

Everything can look very innocent until we start looking at the details .
In this case:
1. When using 10m/s at the inlet, we get the first image in attachment: "10m_s-absolute_pressure.png"
• The pressure range in the "topAir" section is only of about 2 Pa.
• The wall at the bottom is showing the temperature, so that there is a gradient being viewed
2. For the same speed, if we use the "Calculator" filter to do "p-10000", we get what's in the second image: "10m_s-relative_pressure.jpg"
• The actual pressure range is of 1.95 Pa. Still looks innocent enough, doesn't it?
3. When using 1m/s at the inlet, we get the 3rd image, when seeing the "p" field: "1m_s-absolute_pressure.png"
• "1e+5 = 100000" - Which means that ParaView seems to be having a rendering problem. Well, not exactly, since there is a pressure gradient in the "topAir" region.
4. If we do the same calculation as before, we get the 4th image: "1m_s-relative_pressure.png"
• The pressure range is only of 0.102 Pa! This is the reason why the solver seems to hint that there is nothing more to solve!
• This is because the numerical precision needed requires "100000.102" = 9 digits, while the previous one only needs around 6 digits. Well, I'm oversimplifying the digits needed, but this is the ballpark numbers.
So, what now? Let me see if I can find a few posts on the topic of why simpleFoam doesn't use absolute pressure values... OK, I'm having trouble finding threads on this topic
Very well, let me briefly re-explain: simpleFoam solves the "kinematic pressure", instead of the actual pressure, because the density "rho" is constant. For air, this might seem innocuous, but for water with a "rho=1000", this means that there is somewhat of an inconsistency in the scales of values at play here. For example, rho=1000 versus U=1 can potentially lead to numerical errors a bit larger than if rho=1 versus U=1.
A similar situation occurs in this scenario with "U=1 m/s", because the pressure field varies very little, therefore using the absolute pressure is counter-productive.

If we take a look at the tutorial case "heatTransfer/chtMultiRegionFoam/multiRegionLiquidHeater" in OpenFOAM 2.2.x:
• The region "bottomWater" uses for the file "constant/bottomWater/thermophysicalProperties" the following settings:
Code:
```thermoType
{
type            heRhoThermo;
mixture         pureMixture;
transport       const;
thermo          hConst;
equationOfState rhoConst;
specie          specie;
energy          sensibleEnthalpy;
}

mixture
{
specie
{
nMoles          1;
molWeight       18;
}
equationOfState
{
rho             1000;
}
thermodynamics
{
Cp              4181;
Hf              0;
}
transport
{
mu              959e-6;
Pr              6.62;
}
}

// ************************************************************************* //```
In other words, it uses fixed density.
• If we take a look at the file "system/bottomWater/changeDictionaryDict", the pressure values are set to "0", instead of "100000". If I remember correctly, this is only possible/practical because of the thermodynamic properties defined in the previous file.
If you do the same for your case without gravity, i.e. if you define the fluid properties in a way that makes the air pretty much incompressible, then you should be able to no longer have these issues.

The other possibility is to assume that when the pressure equation reaches this numerical precision, it could be considered as converged... although it doesn't feel like it in this case. And since residual controls don't work in multi-region cases, at least last time I checked: http://www.cfd-online.com/Forums/ope...implefoam.html - then this is not exactly the best option.

There is another possibility, which is to rely on reducing the tolerance requested from "fvSolution" after several iterations. This can be achieved with the function object "timeActivatedFileUpdate": https://github.com/OpenFOAM/OpenFOAM...te/controlDict

Best regards,
Bruno
Attached Images
 10m_s-absolute_pressure.png (20.0 KB, 17 views) 10m_s-relative_pressure.jpg (35.2 KB, 12 views) 1m_s-absolute_pressure.png (18.8 KB, 11 views) 1m_s-relative_pressure.png (24.1 KB, 10 views)
__________________

January 15, 2016, 12:07
#6
New Member

Karl Lindqvist
Join Date: Jul 2012
Posts: 20
Rep Power: 7
Dear Bruno,
Thank you for your helpful diagnosis! That is indeed quite a big difference in magnitude between the pressure gradient and absolute value of pressure!

Setting rho = constant, as in the tutorial you mention, is not an option in my application. Fortunately, OpenFOAM provides an "incompressible" version of the ideal gas law - incompressiblePerfectGas.

I was not aware of this option since it's not used in any tutorial, but it allows the definition of a reference pressure to calculate density which is different from the pressure used in calculations.

For completeness, I have attached a case that happily solves down to a residual of 1e-10 at U = 0.1 m/s.

Best regards,
Karl
Attached Files
 partiallyHeatedFloor2D.tar.gz (7.5 KB, 32 views)

 Tags chtmultiregionfoam, chtmultiregionsimplefoam, gravity, p_rgh

 Thread Tools Display Modes Linear Mode

 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 OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post salazardetroya OpenFOAM Running, Solving & CFD 1 July 29, 2015 22:34 incompressible OpenFOAM Running, Solving & CFD 1 July 5, 2015 11:09 virendra_p FLUENT 0 May 1, 2014 08:42 7islands OpenFOAM 2 August 26, 2008 05:25 edvardsenpriv Open Source Meshers: Gmsh, Netgen, CGNS, ... 2 December 6, 2005 06:23

All times are GMT -4. The time now is 10:02.