Dear Aidyn,

Thank you for raising this issue.

I am not too familiar with the AUSM Riemann solvers so I will need to look into it a bit more carefully, but the reference you suggest seems the correct one and I agree with your corrections.

I tried quickly to test a 2D shock tube with your corrections but the solution diverges after 11 iterations, there might be other errors.

I opened an issue on the Nektar++ gitlab so that we can sort this out.

One question, does your subsonic case run with the modifications you suggested?

Cheers,

Giacomo

On Thu, Jan 10, 2019 at 4:23 PM Aidyn Aitzhan <aidyn.b.aitzhan@gmail.com> wrote:

Dear Nektar++ Developers,

I was working on CompressibleFlowSolver and noticed, that AUSM3 Riemann solver behaves unstable attempting to solve subsonic flows.

I was getting a NaN error every time I tried to use it at low Mach flows. Then I figured out that there is a bug in the code. As I understood AUSM3 is

AUSM+up solver from paper by Meng-Sing Liou, "A sequel to AUSM, Part II: AUSM+-up for all speeds," Journal of Computational Physics, Vol. 214, No. 1, 2006, pp. 137-170.

Comparing the code and section "3.3. Algorithm: AUSM+-up for all speeds" in the paper I made several corrections.

        // Parameters for specify the upwinding
        // Note: if fa = 1 then AUSM3 = AUSM2 // AUSM3 = AUSM3
        NekDouble Mco    = 0.01;
        NekDouble Mtilde = 0.5 * (ML * ML + MR * MR);
        NekDouble Mo     = sqrt(std::min(1.0, std::max(Mtilde, Mco*Mco))); //std::min(1.0, std::max(Mtilde, Mco*Mco));
        NekDouble fa       = Mo * (2.0 - Mo);
        NekDouble beta    = 0.125;
        NekDouble alpha  = 0.1875*(5.0*fa*fa - 4.0); //0.1875;
        NekDouble sigma = 1.0;
        NekDouble Kp      = 0.25;
        NekDouble Ku      = 0.75;
        NekDouble rhoA   = 0.5 * (rhoL + rhoR);
        NekDouble Mp     = -(Kp / fa) * ((pR - pL) / (rhoA * cA * cA)) *
                            std::max(1.0 - sigma * Mtilde, 0.0);
       
        NekDouble Mbar   = M4Function(0, beta, ML) +
                           M4Function(1, beta, MR) + Mp;
       
        NekDouble pu     = -2.0 * Ku * rhoA * cA * fa * (uR - uL) *  // -2.0 * Ku * rhoA * cA * cA * (MR - ML) *
                           P5Function(0, alpha, ML) * P5Function(1, alpha, MR);

An AUSM3Solver.cpp file is also attached.

If you are agree with me please update the code in repository.

Best regards,

Aidyn

_______________________________________________
Nektar-users mailing list
Nektar-users@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/nektar-users