CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Running, Solving & CFD

chtMultiRegionSimpleFoam fails to solve cases without gravity

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

Like Tree2Likes
  • 1 Post By wyldckat
  • 1 Post By karlli

Reply
 
LinkBack Thread Tools Display Modes
Old   January 4, 2016, 12:05
Question chtMultiRegionSimpleFoam fails to solve cases without gravity
  #1
New Member
 
Karl Lindqvist
Join Date: Jul 2012
Location: Trondheim, Norway
Posts: 18
Rep Power: 6
karlli is on a distinguished road
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
File Type: png p_rgh_10.png (7.2 KB, 7 views)
Attached Files
File Type: gz heatedFloor2D.tar.gz (5.9 KB, 8 views)
karlli is offline   Reply With Quote

Old   January 4, 2016, 16:25
Default
  #2
Senior Member
 
Alex
Join Date: Oct 2013
Posts: 332
Rep Power: 14
zfaraday will become famous soon enough
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!
zfaraday is offline   Reply With Quote

Old   January 5, 2016, 04:25
Default
  #3
New Member
 
Karl Lindqvist
Join Date: Jul 2012
Location: Trondheim, Norway
Posts: 18
Rep Power: 6
karlli is on a distinguished road
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.
karlli is offline   Reply With Quote

Old   January 8, 2016, 15:37
Default
  #4
New Member
 
Karl Lindqvist
Join Date: Jul 2012
Location: Trondheim, Norway
Posts: 18
Rep Power: 6
karlli is on a distinguished road
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
File Type: gz partiallyHeatedFloor2D.tar.gz (141.8 KB, 2 views)
karlli is offline   Reply With Quote

Old   January 9, 2016, 15:01
Default
  #5
Super Moderator
 
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 9,748
Blog Entries: 39
Rep Power: 103
wyldckat is a glorious beacon of lightwyldckat is a glorious beacon of lightwyldckat is a glorious beacon of lightwyldckat is a glorious beacon of lightwyldckat is a glorious beacon of light
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
File Type: png 10m_s-absolute_pressure.png (20.0 KB, 15 views)
File Type: jpg 10m_s-relative_pressure.jpg (35.2 KB, 10 views)
File Type: png 1m_s-absolute_pressure.png (18.8 KB, 9 views)
File Type: png 1m_s-relative_pressure.png (24.1 KB, 8 views)
arvindpj likes this.
__________________
wyldckat is offline   Reply With Quote

Old   January 15, 2016, 12:07
Default
  #6
New Member
 
Karl Lindqvist
Join Date: Jul 2012
Location: Trondheim, Norway
Posts: 18
Rep Power: 6
karlli is on a distinguished road
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
File Type: gz partiallyHeatedFloor2D.tar.gz (7.5 KB, 21 views)
wyldckat likes this.
karlli is offline   Reply With Quote

Reply

Tags
chtmultiregionfoam, chtmultiregionsimplefoam, gravity, p_rgh

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
Can OpenFoam solve this problem? salazardetroya OpenFOAM Running, Solving & CFD 1 July 29, 2015 22:34
Additional 2D solve at an inlet patch during each iteration of a usual 3D solve incompressible OpenFOAM Running, Solving & CFD 1 July 5, 2015 11:09
The issue with gravity!!! virendra_p FLUENT 0 May 1, 2014 08:42
Directory constructions of complex cases 7islands OpenFOAM 2 August 26, 2008 05:25
Importing mesh for complex 2D cases edvardsenpriv Open Source Meshers: Gmsh, Netgen, CGNS, ... 2 December 6, 2005 06:23


All times are GMT -4. The time now is 07:11.