
[Sponsors] 
January 4, 2016, 12:05 
chtMultiRegionSimpleFoam fails to solve cases without gravity

#1 
New Member
Karl Lindqvist
Join Date: Jul 2012
Location: Trondheim, Norway
Posts: 18
Rep Power: 6 
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.086793044719327e06, No Iterations 1 DILUPBiCG: Solving for Uy, Initial residual = 0.0003276496308532365, Final residual = 9.054436042068467e06, No Iterations 1 DILUPBiCG: Solving for h, Initial residual = 0.0001652516759873468, Final residual = 1.096169895436815e05, No Iterations 1 Min/max T:299.9999249381647 485.7994890854322 GAMG: Solving for p_rgh, Initial residual = 7.12858129068406e05, Final residual = 7.03830111216478e07, No Iterations 39 time step continuity errors : sum local = 4.644778702141491e09, global = 3.757627749083547e11, cumulative = 6.025277233104804e08 Min/max rho:0.7154943606653349 1.158623823918645 Solving for solid region wall DICPCG: Solving for h, Initial residual = 0.006403279730421732, Final residual = 3.900945002001053e05, 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.069908791677873e06, No Iterations 1 DILUPBiCG: Solving for Uy, Initial residual = 0.0003204160935846439, Final residual = 9.270251346231881e06, No Iterations 1 DILUPBiCG: Solving for h, Initial residual = 0.0001638593213157633, Final residual = 1.087279304981332e05, No Iterations 1 Min/max T:299.9999248933859 485.9034691736506 GAMG: Solving for p_rgh, Initial residual = 6.801409757909533e05, Final residual = 7.818732180590184e07, No Iterations 1000 time step continuity errors : sum local = 4.87402540356367e09, global = 3.229735324601451e11, cumulative = 6.022047497780203e08 Min/max rho:0.7153412496577295 1.158623824019448 Solving for solid region wall DICPCG: Solving for h, Initial residual = 0.006401407648684744, Final residual = 3.900699891191313e05, 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.053164802211731e06, No Iterations 1 DILUPBiCG: Solving for Uy, Initial residual = 0.0003136746051429199, Final residual = 9.491366256280922e06, No Iterations 1 DILUPBiCG: Solving for h, Initial residual = 0.0001624826897776338, Final residual = 1.078451248755692e05, No Iterations 1 Min/max T:299.9999248494713 486.0066718410475 GAMG: Solving for p_rgh, Initial residual = 6.507513719961976e05, Final residual = 7.436642579933397e07, No Iterations 1000 time step continuity errors : sum local = 4.695772893242314e09, global = 3.700505525107244e11, cumulative = 6.018346992255096e08 Min/max rho:0.7151893481877345 1.158623824119854 Solving for solid region wall DICPCG: Solving for h, Initial residual = 0.006399522667559255, Final residual = 3.900445715954452e05, 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 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 blockstructured 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 

January 4, 2016, 16:25 

#2 
Senior Member
Alex
Join Date: Oct 2013
Posts: 334
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
Location: Trondheim, Norway
Posts: 18
Rep Power: 6 
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 xdirection rather than the ydirection (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
Location: Trondheim, Norway
Posts: 18
Rep Power: 6 
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 1e4.
Any suggestions? Regards, Karl 

January 9, 2016, 15:01 

#5 
Super Moderator
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 9,837
Blog Entries: 39
Rep Power: 108 
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:
Very well, let me briefly reexplain: 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.
If we take a look at the tutorial case "heatTransfer/chtMultiRegionFoam/multiRegionLiquidHeater" in OpenFOAM 2.2.x:
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 multiregion cases, at least last time I checked: http://www.cfdonline.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
__________________


January 15, 2016, 12:07 

#6 
New Member
Karl Lindqvist
Join Date: Jul 2012
Location: Trondheim, Norway
Posts: 18
Rep Power: 6 
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 1e10 at U = 0.1 m/s. Best regards, Karl 

Tags 
chtmultiregionfoam, chtmultiregionsimplefoam, gravity, p_rgh 
Thread Tools  
Display Modes  


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 