
[Sponsors] 
buoyantFoam crashes when solving mixed convection or using higher order schemes 

LinkBack  Thread Tools  Search this Thread  Display Modes 
June 27, 2024, 09:43 
buoyantFoam crashes when solving mixed convection or using higher order schemes

#1 
New Member
talop
Join Date: Feb 2023
Posts: 6
Rep Power: 3 
Hi all,
Long time reader, first time poster I have been trying to simulate convective heat losses from an isothermal, horizontal cylinder under cross flow in the mixed and natural convection regimes at steady state. This steady state solution will be used to initialize a transient case. Upon successfully simulating this reference case, I plan to move on to rotating bodies and detachededdy simulations. However, I was unable to get a solution, let alone comparing it with the literature. The problem geometry consits of a horizontal cylinder with a diameter of 1 meter placed in a 10x10 m square cross section duct. The cylinder is placed 6.5 m downstream of the inlet and the outlet is 12.5 m downstream from the cylinder. The mesh is generated in ANSYS Fluent using the polyhexcore method with quality above 0.5. Running checkMesh does not give any errors. Maximum skewness is 3.65 and nonorthogonality is maximum 59.85 with an average of 8.10. For comparison two more meshes were generated in ANSYS Meshing, one using a mapped mesh with hexahedral elements, other using tetrahedrons outside the boundary layer. So far, I have been able to run some cases with significant forced convection, using the upwind schemes only, with kepsilon and komega SST turbulence models. The air was flowing in at a rate 1 m/s. The inlet temperature was 288.15 K while hot cylinder was at 600 K (Re ~ 24000, Gr ~ 3.97e+09 and Ri ~ 6.9). 1) If a case dominated by natural convection is to be simulated (obtained by reducing the inlet velocity), solution crashes during p_rgh iterations after a few time steps. 2) Similarly, if the divergence scheme for the velociy term "div(phi,U)" is not "bounded Gauss upwind", the simulation crashes regardless of the inlet velocity at the same stage but at a different "time step" (anywhere between 5 to 500). 3) The same mesh performs well in ANSYS Fluent while using both segregated (SIMPLE) and coupled solver with second order upwind schemes except for the kequation with the same relaxation factors or higher default values. SIMPLE solver takes ~1000 iterations to reach a solution. I have tried to reduce the relaxation factors, use more correctors, use "type coupled" keyword for solvers and initialize the internal field with the inlet values. Using the coupled solver and reducing the relaxation factors helped slightly but the stability problem is still there. Using the blockMesh is not possible for me because I will move towards more complicated geometries, which will be difficult , if not impossible, to mesh with blockMesh. I have tried using the snappyHexMesh in a different setting but I was unable to obtain good boundary layer meshes around internal corners. Therefore I will stick with the ANSYS Mechanical or Fluent as a mesher for now. My guess is that I am using suboptimal settings at best, leading to divergent behavior for higher order schemes or highly transient caeses. I will be glad if more experienced users can suggest some better settings for running mixed convection simulations using twoequation models and/or DES. You can find the case files (except mesh because of size limitations) in the attached .zip file. Best Last edited by talop; July 4, 2024 at 12:22. 

June 28, 2024, 04:28 

#2 
New Member
talop
Join Date: Feb 2023
Posts: 6
Rep Power: 3 
Update: Using a hex mesh gave slightly better results but simulations for the natural and mixed convection cases still blow up even with upwind schemes.
Best, talop Last edited by talop; July 2, 2024 at 09:29. 

July 1, 2024, 03:49 

#3 
Senior Member
Join Date: Dec 2021
Posts: 230
Rep Power: 5 
Hey!
Natural convection has been a long time nemesis of mine too! From what I have read and tried, transient or pseudotransient are better at this kind of simulation. If you are using the ESI version, it is also possible to use hydrostaticPressure function object to initialize p and p_rgh with a hydrostatic distribution and a ph_rgh field. I found it can help. Mesh wise, hexaedral is indeed better from my experience. How low did you get with the relaxation factors? 

July 1, 2024, 04:02 

#4 
New Member
talop
Join Date: Feb 2023
Posts: 6
Rep Power: 3 
Hi, Alczem
I have tried initializing the flow field with some velocity to designate the flow direction and it did work partially. At least the simulations do not crash now. But it did not solve the negative p_rgh and back flow problem completely. I have seen that the fluid "collapses" under gravity and this leads to a backflow. To solve this, I tried initializing the pressure field (0/p) with overpressure but this did not solve the back flow either. I will try the hydrostaticPressure to see if it helps. The fvSchemes uses "stadyState" for the ddtSchemes and therefore it runs in the SIMPLE mode even if it reads PIMPLE in the fvSolution file below. Code:
solvers { p_rgh { // type coupled; // segregated; solver GAMG; // PBiCICG; // GAMG; tolerance 1e7; relTol 0.01; smoother DICGaussSeidel; preconditioner DILU; } "(Uhkepsilonomega)" { // type coupled; // segregated; solver PBiCGStab; // PBiCICG; // PBiCGStab; preconditioner DILU; tolerance 1e8; relTol 0.1; } } PIMPLE { momentumPredictor no; nCorrectors 5; // 3; nNonOrthogonalCorrectors 5; // 3; nOuterCorrectors 2; // 1; pRefCell 0; pRefValue 0; residualControl { p_rgh 1e4; U 1e4; h 1e4; // possibly check turbulence fields "(kepsilonomega)" 1e3; } } relaxationFactors { fields { rho 0.85; // 1.0; p_rgh 0.5; // 0.7; } equations { U 0.5; // 0.3; h 0.3; "(kepsilonomega)" 0.4; // 0.7; } } // ************************************************************************* // 

July 1, 2024, 05:39 

#5 
Senior Member
Join Date: Dec 2021
Posts: 230
Rep Power: 5 
Hey,
Regarding your fvSchemes, I would use limited 0.33 instead of orthogonal since your mesh is not perfectly orthogonal. You could also try to add a limiter to the gradient calculation. Why does your p_rgh values are set to 0 instead of something like 101325? I noticed you set pRefValue to 0 in fvSolution and there is a pRef file with 101325 in constant. I am not familiar with OF10, but are you sure the pressure is properly initialized? Again, I don't know exactly how it should be configured. Cheers! 

July 4, 2024, 12:02 

#6 
New Member
talop
Join Date: Feb 2023
Posts: 6
Rep Power: 3 
Hi,
I have solved some of the issues and implemented better schemes. Getting there slowly The p_rgh = p  rho*g*(h  h_ref) is the pressure without hydrostatic effects w.r.t. a reference point and gravity vector. As far as I know, it is defined as gauge pressure referencing the pRef file. In the tutorial cases this definition was used in the tutorial cases. After tinkering with the case more, I have identified two problems: 1) constant/physicalProperties: equationOfState perfectGas; The perfectGas uses variable pressure in rho = P/RT. FOr this case, pressure variation is small so incompressiblePerfectGas can be used with a reference pressure of 101325 Pa. 2) system/fvSolution: nOuterCorrectors 3; By increasing the number of outer correctors to 1015, I have been able to use higher order schemes for the divergence terms (limitedLinear, vanLeer ... etc.) 

Tags 
buoyantfoam, external mesh, mixed convection, unstable solver 
Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
chtMultiRegionSimpleFoam: maximum number of iterations excedeed.  Nkl  OpenFOAM Running, Solving & CFD  19  October 10, 2019 02:42 
conjugate heat transfer in OpenFOAM  skuznet  OpenFOAM Running, Solving & CFD  99  March 16, 2017 05:07 
simpleFoam error  "Floating point exception"  mbcx4jc2  OpenFOAM Running, Solving & CFD  12  August 4, 2015 02:20 
Moving mesh  Niklas Wikstrom (Wikstrom)  OpenFOAM Running, Solving & CFD  122  June 15, 2014 06:20 
Unstabil Simulation with chtMultiRegionFoam  mbay101  OpenFOAM Running, Solving & CFD  13  December 28, 2013 13:12 