[ishan_ae]

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:

/*---------------------------------------------------------------------------*/

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((phi|interpolate(rho)),p)  Gauss upwind;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}


fluxRequired
{
    default         no;
    p;
}


wallDist
{
    method meshWave;
}

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       1e-06;
        relTol          0.01;
    }

    p
    {
        $Phi;
    }
}

potentialFlow
{
    nNonOrthogonalCorrectors 50;
}
*/

solvers
{

"p.*" 
{
        solver           GAMG;
        tolerance        1e-10;
        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       1e-06;
	relTol          0.05;
}
*/
"U.*"
{
        solver           smoothSolver;
        smoother         GaussSeidel;
        tolerance        1e-10;
        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 
"(k|epsilon|omega|nuTilda|kFinal|epsilonFinal|omegaFinal|e|rho)"
{
        solver           smoothSolver;
        smoother         GaussSeidel;
        tolerance        1e-10;
        relTol           0.1;
        nSweeps          1;
	maxIter		200;
}


}

//copy setings from U-bend tutorial 
//check for transonic options 
SIMPLE
{
    nNonOrthogonalCorrectors 5;
    //pMinFactor      0.1;
    //pMaxFactor      2;
    transonic       yes;
    consistent      yes;

    residualControl
    {
        p               1e-5;
        U               1e-4;
        k               1e-4;
        omega           1e-4;
	nuTilda         1e-4;
	e		1e-3;
    }
}


//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);
}

Can someone please suggest some troubleshooting points on this ?

[SvenBent]

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

[ishan_ae]

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 re-look 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.

[ishan_ae]

Update: I have also tried running a sonicFoam simulation but it doesn't show any time-dependent 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 re-check.

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.

[ishan_ae]

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 first-order 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 un-informed trial and error approach.

[pgh]

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

[ishan_ae]

Thanks for your input. I reduced the rho field relaxation to 0.01 and it still showed the same behavior.

[pgh]

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

[ishan_ae]

I am using bounded schemes and now my setup looks like this:

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((phi|interpolate(rho)),p)  Gauss upwind;
   // div((phi|interpolate(rho)),p)  Gauss upwind;

}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}


fluxRequired
{
    default         no;
    p;
}


wallDist
{
    method meshWave;
}

It still shows the same behavior. I have observed that after a while the maximum velocities that I am monitoring start accumulating on the flowRateOutletVelocity patch. I have been wondering from the start whether this issue is with the schemes or the BC.


EDIT- I added bounded for div(phiv,p), div(phid,p) and also div(phi,Ekp) and starting it again.

[shock77]

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.

[ishan_ae]

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 1e-2. The global time step continuity errors are also not decreasing. They are stuck around 0.09 levels.

[ishan_ae]

Update: I managed to get a stable solution on the case with total pressure imposed at inlet and outlet. I started by using first-order scheme and low relaxation factors. But, the global continuity errors don't decrease and also pressure residuals are not going below 1e-2 levels. Rest of the residuals are at least at 1e-4 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.

[ishan_ae]

Update: I rechecked the overall case setup and I have narrowed down the possible issues to these points:

1. At the outlet patch where mass flow rate is imposed using flowRateOutletVelocity I have noticed that there is an nonphysical increase in p,U, T and rho quantities over the iterations where as inside the domain the quantities seem normal and follow an expected pattern. Only a few iterations before the crash the quantities inside domain increase to nonphysical values. This makes me think that this BC is not suited for my case or maybe the flow is reversing ?
2. I also found that I had used incorrect pressure values at the inlet forsimpleFoam setup and maybe these results are causing issues. I am initializing with simpleFoam values. To be specific, I am extracting only internalField values from simpleFoam converged solution.

Also for schemes I am using bounded first order initially and for solution controls I am using 0.2 and 0.01 for p and rho respectively.

With this in mind can someone please provide suggestions based on the following points:

1. Which BC apart from flowRateOutletVelocity can be used in my case ? I have already tried imposing total pressure at the outlet but the system goes unstable after a while.
2. How can I correctly setup in-viscid compressible flow simulations for my case? I think the initial values are definitely playing a major role here and intend to use in-viscid solution as initial.
3. Is lowering of inlet total pressure initially recommended ? If yes, how should I go about using this method? I have already tried using uniformTotalPressure on a clean case with values starting from 60,000 Pa. That did not work.
4. Finally, can someone please check my case setup with mass flow rate imposed at the outlet and pressure at inlet? I have attached a RAR file here.

mySetup.tar.gz

[arun1994]

for the inlet, maybe you can try surfaceNormal velocity inlet condition.