[karlli]

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):

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

[zfaraday]

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

[karlli]

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.

[karlli]

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

[wyldckat]

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.

• Here's a similar explanation I wrote almost a year ago on numerical precision: http://www.cfd-online.com/Forums/ope...tml#post529800 - post #4

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:
  

  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

[karlli]

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