
[Sponsors] 
October 15, 2020, 02:16 
rhoSimpleFoam crash help

#1 
Member
ishan
Join Date: Oct 2017
Posts: 77
Rep Power: 8 
Hello.
I am simulating a compressible subsonic flow in a duct like geometry using rhoSimpleFoam. I am initializing this solution using simpleFoam solution. The rhoSimpleFoam solution keeps on crashing and I am am trying to understand why this is happening. I am monitoring maximum velocity and maximum temperature in the domain. The maximum velocity reaches nonphysical levels but the temperature also increases but upto to acceptable values. The maximum Mach number should be close to 0.7. Interestingly, the yplus in the domain always remains less than 1 through the simulation. With regards to the boundary conditions, I am imposing total pressure at the inlet and mass flow rate at the outlet. I have linked here the images of my residuals and monitor plots. Attached here are my fvSchemes and fvSolution files: Code:
/**/ FoamFile { version 2.0; format binary; class dictionary; location ""; object fvSchemes; } /**/ /**/ ddtSchemes { default steadyState; } gradSchemes { default Gauss linear; limited cellLimited Gauss linear 1; grad(U) $limited; grad(k) $limited; grad(omega) $limited; } divSchemes { default none; div(phi,U) Gauss linearUpwind limited; turbulence Gauss linearUpwind limited; energy Gauss linearUpwind limited; div(div(phi,U)) Gauss linear; div(phi,k) $turbulence; div(phi,omega) $turbulence; div(phi,e) $energy; div(phi,K) $energy; div(phi,Ekp) $energy; div(phiv,p) Gauss upwind; //for epsilon //div(phi,epsilon) bounded Gauss upwind; div(phid,p) Gauss upwind; div(phi,Ekp) Gauss upwind; //div((nuEff*dev2(T(grad(U))))) Gauss linear; //div(((rho*nuEff)*dev2(T(grad(U))))) Gauss upwind; div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear; //div((phiinterpolate(rho)),p) Gauss upwind; } laplacianSchemes { default Gauss linear corrected; } interpolationSchemes { default linear; } snGradSchemes { default corrected; } fluxRequired { default no; p; } wallDist { method meshWave; } Code:
FoamFile { version 2.0; format binary; class dictionary; location ""; object fvSolution; } /**/ /**/ FoamFile { version 2.0; format binary; class dictionary; location ""; object fvSolution; } /**/ /**/ /* solvers { Phi { solver GAMG; smoother DIC; tolerance 1e06; relTol 0.01; } p { $Phi; } } potentialFlow { nNonOrthogonalCorrectors 50; } */ solvers { "p.*" { solver GAMG; tolerance 1e10; relTol 0.05; smoother GaussSeidel; nPreSweeps 0; nPostSweeps 2; cacheAgglomeration on; agglomerator faceAreaPair; nCellsInCoarsestLevel 10; mergeLevels 1; maxIter 20; } /* p { solver PCG; preconditioner DIC; tolerance 1e06; relTol 0.05; } */ "U.*" { solver smoothSolver; smoother GaussSeidel; tolerance 1e10; relTol 0.1; nSweeps 1; } Phi { $p; } //limit the max iterations to solve energy equations //make a separate setup for solving e instead of clubbing it with turbulence "(kepsilonomeganuTildakFinalepsilonFinalomegaFinalerho)" { solver smoothSolver; smoother GaussSeidel; tolerance 1e10; relTol 0.1; nSweeps 1; maxIter 200; } } //copy setings from Ubend tutorial //check for transonic options SIMPLE { nNonOrthogonalCorrectors 5; //pMinFactor 0.1; //pMaxFactor 2; transonic yes; consistent yes; residualControl { p 1e5; U 1e4; k 1e4; omega 1e4; nuTilda 1e4; e 1e3; } } //change the relaxation factors potentialFlow { nNonOrthogonalCorrectors 50; } relaxationFactors { fields { p 0.3; rho 0.3; } equations { U 0.3; k 0.3; omega 0.3; epsilon 0.3; nuTilda 0.3; e 0.2; } } cache { grad(U); } 

October 15, 2020, 04:23 

#2 
New Member
Mathias Poulsen
Join Date: Feb 2018
Location: Denmark
Posts: 9
Rep Power: 8 
Hi Ishan.
Relaxing your pressure field when using SIMPLEC (consistent yes) can sometimes cause troubles. I recommend that you either turn consistent off, or that you don't relax your p field. hope it helps. Best Mathias 

October 15, 2020, 06:33 

#3 
Member
ishan
Join Date: Oct 2017
Posts: 77
Rep Power: 8 
Hello.
Thanks for your inputs. I switched off the consistent option but that still showed the same behavior. I took a look at the final iteration's crashed solution and noticed on the outlet patch that the maximum velocities of around 1000 were occurring. Interestingly, within the domain the values are all right. Perhaps the specified flow rate BC is not suited for my case ? I also tried using fixeValue BC on the outlet with patch values generated from the simpleFoam solution. It is running right now and it shows weird behavior. I will relook at my boundary conditions and also numerical schemes. My guess is that there are strong transonic flows occurring and this is causing problems? Not sure. 

October 16, 2020, 00:59 

#4 
Member
ishan
Join Date: Oct 2017
Posts: 77
Rep Power: 8 
Update: I have also tried running a sonicFoam simulation but it doesn't show any timedependent pattern on quantities of interest.
Now, I am thinking the issue seems to be more fundamental. Perhaps some incorrect boundary condition ? Unstable scheme for the velocity terms? Maybe some wrong form of data file I am using for internalField initialization from simpleFoam data ? Need to recheck. It is strange that only the U values are increasing whereas the T values remain within logically physical limits. I am unable to understand what is causing this nonphysical increase because the U residuals seem ok to me. 

October 16, 2020, 10:10 

#5 
Member
ishan
Join Date: Oct 2017
Posts: 77
Rep Power: 8 
Update: So I changed the boundary conditions from flowRateOutletVelocity to totalPressure at the outlet. Now, the model has total pressure values imposed both on inlet and outlet. I am using limiters and firstorder schemes. I have also reduced the relaxation factors with a maximum of 0.3.
The monitor points are oscillating/unstable and the pressure and velocity residuals look like they are on a trajectory to Mars. I will let it run for a few more iterations before solver crashes. I approximated the total pressure value at the outlet from simpleFoam solution. That was not a good idea. I followed the BC recommendations mentioned here: https://www.openfoam.com/documentati...binations.html Maybe someone with more experience/knowledge can help me overcome the uninformed trial and error approach. 

October 20, 2020, 04:27 

#6 
Senior Member
Pawan Ghildiyal
Join Date: Nov 2015
Posts: 135
Rep Power: 10 
try with relaxing rho to some small value . i.e may be 0.01 and see behaviour. Since it is steady state , it should not affect final solution but will help stabilize solution


October 21, 2020, 06:47 

#7 
Member
ishan
Join Date: Oct 2017
Posts: 77
Rep Power: 8 
Thanks for your input. I reduced the rho field relaxation to 0.01 and it still showed the same behavior.


October 21, 2020, 06:53 

#8 
Senior Member
Pawan Ghildiyal
Join Date: Nov 2015
Posts: 135
Rep Power: 10 
if you are using steady state solver, then you need to use bounded version of schemes,.
div(phi,U) bounded Gauss linearUpwind limited Check the tutorial. For transinet, bounded is not needed but for steady state it is needed 

October 21, 2020, 09:17 

#9 
Member
ishan
Join Date: Oct 2017
Posts: 77
Rep Power: 8 
I am using bounded schemes and now my setup looks like this:
Code:
FoamFile { version 2.0; format binary; class dictionary; location ""; object fvSchemes; } /**/ /**/ ddtSchemes { default steadyState; } gradSchemes { default Gauss linear; limited cellLimited Gauss linear 1; grad(U) $limited; grad(k) $limited; grad(omega) $limited; } divSchemes { default none; div(phi,U) bounded Gauss upwind limited; turbulence bounded Gauss linearUpwind limited; energy bounded Gauss linearUpwind limited; div(div(phi,U)) bounded Gauss linear; div(phi,k) $turbulence; div(phi,omega) $turbulence; div(phi,e) $energy; div(phi,K) $energy; div(phi,Ekp) $energy; div(phiv,p) Gauss upwind; //for epsilon //div(phi,epsilon) bounded Gauss upwind; div(phid,p) Gauss upwind; div(phi,Ekp) Gauss upwind; //div((nuEff*dev2(T(grad(U))))) Gauss linear; //div(((rho*nuEff)*dev2(T(grad(U))))) Gauss upwind; div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear; //div((phiinterpolate(rho)),p) Gauss upwind; // div((phiinterpolate(rho)),p) Gauss upwind; } laplacianSchemes { default Gauss linear corrected; } interpolationSchemes { default linear; } snGradSchemes { default corrected; } fluxRequired { default no; p; } wallDist { method meshWave; } EDIT I added bounded for div(phiv,p), div(phid,p) and also div(phi,Ekp) and starting it again. 

October 21, 2020, 17:00 

#10 
Senior Member
Join Date: Dec 2019
Posts: 215
Rep Power: 7 
Hi,
I am not familiar with rhoSimpleFoam, but I think the solver and schemes should be ok for your case. I think the error is most likely due to wrong BC at the outlet. To reduce the amount of possible mistakes I would suggest running sdome simulations with slip walls and without turbulence models. I dont know what you are modelling tho. 

October 22, 2020, 05:39 

#11 
Member
ishan
Join Date: Oct 2017
Posts: 77
Rep Power: 8 
I changed the setup with total pressure imposed at outlet as well. I have better results as compared to the case where I had imposed mass flow rate at the outlet.
I am monitoring maximum velocity and temperature and the results are oscillatory. Also, the pressure residuals are oscillating don't go below 1e2. The global time step continuity errors are also not decreasing. They are stuck around 0.09 levels. 

October 23, 2020, 03:08 

#12 
Member
ishan
Join Date: Oct 2017
Posts: 77
Rep Power: 8 
Update: I managed to get a stable solution on the case with total pressure imposed at inlet and outlet. I started by using firstorder scheme and low relaxation factors. But, the global continuity errors don't decrease and also pressure residuals are not going below 1e2 levels. Rest of the residuals are at least at 1e4 levels.
For the schemes I changed from upwind to linearUpwind limited for energy terms, and div(phi,U). They are bounded. I am wondering if now the issue is with the pressure matrix solver. I am using GAMG. 

November 6, 2020, 07:40 

#13 
Member
ishan
Join Date: Oct 2017
Posts: 77
Rep Power: 8 
Update: I rechecked the overall case setup and I have narrowed down the possible issues to these points:
With this in mind can someone please provide suggestions based on the following points:


March 8, 2021, 05:39 

#14 
Member
Arun subramanian
Join Date: Jun 2016
Location: Florence,Italy
Posts: 48
Rep Power: 10 
for the inlet, maybe you can try surfaceNormal velocity inlet condition.


Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
transsonic nozzle with rhoSimpleFoam  Unseen  OpenFOAM Running, Solving & CFD  8  July 1, 2022 06:54 
rhoSimpleFoam crash at first step  fxzf  OpenFOAM Running, Solving & CFD  2  August 2, 2017 04:13 
rhoSimpleFoam angledDuctExplicitFixedCoeff tutorial fails in parallel  donQi  OpenFOAM Running, Solving & CFD  1  February 22, 2016 19:47 
rhoSimpleFoam. patchField error.  123  OpenFOAM Running, Solving & CFD  4  June 6, 2014 15:22 
Problem with rhoSimpleFoam  mecbe2002  OpenFOAM  3  April 11, 2010 00:54 