[MikeC]

I'm trying to setup a canonical supersonic cavity-in-duct flow simulation (2D for now) using rhoCentralFoam, but am having stability issues with the solvers when using anything but extremely small timesteps (CFL < 0.3). When using larger timesteps I get oscillations emanating from the top and bottom walls at the inlet on startup which eventually produce a floating point exception. I'm pretty new to openfoam and cfd in general, so I'm not sure if the issue is with my numerical schemes, my boundary conditions, or if I need to initialize better. Any suggestions are greatly appreciated!

My setup is:

BCs:
U - fixedValue at inlet, noSlip on walls, zeroGradient on outlet
p - zeroGradient at inlet and walls, waveTransmissive on outlet
T - fixedValue at inlet, zeroGradient on walls and outlet
k - fixedValue on inlet and walls, inletOutlet at outlet
alphat & nut - zeroGradient everywhere

Initialization:
Quiescent everywhere

Solvers:

{
    "(rho|rhoU|rhoE)"
    {
        solver          diagonal;
    }

    "(k|nut)"
    {
        solver          PBiCGStab;
        preconditioner	DILU;
        nSweeps         2;
        tolerance       1e-10;
        relTol          0;
    }

    U
    {
	solver		PCG;
	preconditioner	FDIC;
	tolerance	1e-10;
	relTol		0;	 
    }

    e
    {
        $U;
        tolerance       1e-10;
        relTol          0;
    }
}

Schemes:

fluxScheme      Kurganov;

ddtSchemes
{
    default         CrankNicolson 0.9;
}

gradSchemes
{
    default         cellLimited Gauss linear 1;
}

divSchemes
{
    default         none;
    div(tauMC)      Gauss linear;
    div(phi,k)      Gauss limitedLinear 1;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
    reconstruct(rho) vanLeer;
    reconstruct(U)  vanLeerV;
    reconstruct(T)  vanLeer;
}

snGradSchemes
{
    default         corrected;
}

[MikeC]

wait, is this paper correct? -> https://www.researchgate.net/publica...rk_of_OpenFOAM

is rhoCentralFoam actually an explicit solver, even though I'm using an implicit time-stepping scheme? If so, why am I having to invert matrices?

[chengdi]

Hi, MikeC

I encountered the same problem. I am also very confused.

I just tested the forwardStep example in OpenFOAM 201612+. It is true that CFL no. matters for this "implicit solver". I tested CFL=0.2, 0.3, 0.4 and 0.5, the oscillations behind shock becomes larger for larger CFL, and it goes to error at CFL = 0.5.

I am looking into the source code.

[MikeC]

Hi Chengdi,

I looked into it a little more, and if I'm understanding it correctly part of the system is being discretized explicitly (rhoe, rhoU, etc) and part implicitly (p, etc). The difference is in fvc versus fvm. I've been told that this method allows you to use a larger CFL than you could if you were using a purely explicit solver. Still not sure why the default isn't just a completely implicit solver.

[chengdi]

Yes, I have found the paper "Implementation of semi-discrete, non-staggered central schemes in a colocated, polyhedral, finite volume framework, for high-speed viscous flows".

Because OpenFOAM is a segregate solver which solve variables in cell volume center separately. So you can choose to solve one variable implicitly and other variables explicitly.

For rhoCentralFoam, it construct the KT or KNP numerical flux explicitly. So it must be bounded by CFL number like explicit scheme.

[jnanabrao]

I am replying on this old thread because I have struggled with similar questions. I believe a completely implicit solver will require a U, p and e to be strongly coupled. This would require an iterative semi-explicit method solver like SIMPLE. Usually the intermediate states of such solvers are not very stable. If you read Issa's pioneer paper on PISO (from 1985 or 1986), one of his motivations for the algorithm was to allow for large time steps. This is not stable for explicit solvers as they do not have any weights for the t+1 time steps and can be erroneous when the time variations are high. However, it is possible to avoid the errors of incorrect temperature at the cost of higher computation by increasing the outer corrector loops in rhoPimpleFoam and sonicFoam as they are based on combined PISO and SIMPLE (PIMPLE). In rhoCentralFoam, the only ways to increase stability are reducing the Courant number, using more stable (dissipative) schemes or better initialization, as far as I am aware.