Nektar++
AUSM3Solver.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: AUSM3Solver.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: AUSM3 Riemann solver.
32//
33///////////////////////////////////////////////////////////////////////////////
34
36
37namespace Nektar
38{
39
40std::string AUSM3Solver::solverName =
42 "AUSM3", AUSM3Solver::create, "AUSM3 Riemann solver");
43
45 : AUSM0Solver(pSession)
46{
47 pSession->LoadParameter("Mco", m_Mco, 0.01);
48}
49
50/**
51 * @brief AUSM3 Riemann solver
52 *
53 * @param rhoL Density left state.
54 * @param rhoR Density right state.
55 * @param rhouL x-momentum component left state.
56 * @param rhouR x-momentum component right state.
57 * @param rhovL y-momentum component left state.
58 * @param rhovR y-momentum component right state.
59 * @param rhowL z-momentum component left state.
60 * @param rhowR z-momentum component right state.
61 * @param EL Energy left state.
62 * @param ER Energy right state.
63 * @param rhof Computed Riemann flux for density.
64 * @param rhouf Computed Riemann flux for x-momentum component
65 * @param rhovf Computed Riemann flux for y-momentum component
66 * @param rhowf Computed Riemann flux for z-momentum component
67 * @param Ef Computed Riemann flux for energy.
68 */
69void AUSM3Solver::v_PointSolve(double rhoL, double rhouL, double rhovL,
70 double rhowL, double EL, double rhoR,
71 double rhouR, double rhovR, double rhowR,
72 double ER, double &rhof, double &rhouf,
73 double &rhovf, double &rhowf, double &Ef)
74{
75 // Left and Right velocities
76 NekDouble uL = rhouL / rhoL;
77 NekDouble vL = rhovL / rhoL;
78 NekDouble wL = rhowL / rhoL;
79 NekDouble uR = rhouR / rhoR;
80 NekDouble vR = rhovR / rhoR;
81 NekDouble wR = rhowR / rhoR;
82
83 // Internal energy (per unit mass)
84 NekDouble eL = (EL - 0.5 * (rhouL * uL + rhovL * vL + rhowL * wL)) / rhoL;
85 NekDouble eR = (ER - 0.5 * (rhouR * uR + rhovR * vR + rhowR * wR)) / rhoR;
86 // Pressure
87 NekDouble pL = m_eos->GetPressure(rhoL, eL);
88 NekDouble pR = m_eos->GetPressure(rhoR, eR);
89 // Speed of sound
90 NekDouble cL = m_eos->GetSoundSpeed(rhoL, eL);
91 NekDouble cR = m_eos->GetSoundSpeed(rhoR, eR);
92
93 // Average speeds of sound
94 NekDouble cA = 0.5 * (cL + cR);
95
96 // Local Mach numbers
97 NekDouble ML = uL / cA;
98 NekDouble MR = uR / cA;
99
100 // Parameters for specify the upwinding
101 // Note: if fa = 1 then AUSM3 = AUSM3
102 NekDouble Mtilde = 0.5 * (ML * ML + MR * MR);
103 NekDouble Mo = std::sqrt(std::min(1.0, std::max(Mtilde, m_Mco * m_Mco)));
104 NekDouble fa = Mo * (2.0 - Mo);
105 NekDouble beta = 0.125;
106 NekDouble alpha = 0.1875;
107 NekDouble sigma = 1.0;
108 NekDouble Kp = 0.25;
109 NekDouble Ku = 0.75;
110 NekDouble rhoA = 0.5 * (rhoL + rhoR);
111 NekDouble Mp = -(Kp / fa) * ((pR - pL) / (rhoA * cA * cA)) *
112 std::max(1.0 - sigma * Mtilde, 0.0);
113
114 NekDouble Mbar = M4Function(0, beta, ML) + M4Function(1, beta, MR) + Mp;
115
116 NekDouble pu = -2.0 * Ku * rhoA * cA * cA * (MR - ML) *
117 P5Function(0, alpha, ML) * P5Function(1, alpha, MR);
118
119 NekDouble pbar =
120 pL * P5Function(0, alpha, ML) + pR * P5Function(1, alpha, MR) + pu;
121
122 if (Mbar >= 0.0)
123 {
124 rhof = cA * Mbar * rhoL;
125 rhouf = cA * Mbar * rhoL * uL + pbar;
126 rhovf = cA * Mbar * rhoL * vL;
127 rhowf = cA * Mbar * rhoL * wL;
128 Ef = cA * Mbar * (EL + pL);
129 }
130 else
131 {
132 rhof = cA * Mbar * rhoR;
133 rhouf = cA * Mbar * rhoR * uR + pbar;
134 rhovf = cA * Mbar * rhoR * vR;
135 rhowf = cA * Mbar * rhoR * wR;
136 Ef = cA * Mbar * (ER + pR);
137 }
138}
139
140} // namespace Nektar
double P5Function(int A, double alpha, double M)
double M4Function(int A, double beta, double M)
static RiemannSolverSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession)
Definition: AUSM3Solver.h:46
AUSM3Solver(const LibUtilities::SessionReaderSharedPtr &pSession)
Definition: AUSM3Solver.cpp:44
static std::string solverName
Definition: AUSM3Solver.h:52
void v_PointSolve(double rhoL, double rhouL, double rhovL, double rhowL, double EL, double rhoR, double rhouR, double rhovR, double rhowR, double ER, double &rhof, double &rhouf, double &rhovf, double &rhowf, double &Ef) override
AUSM3 Riemann solver.
Definition: AUSM3Solver.cpp:69
EquationOfStateSharedPtr m_eos
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
@ beta
Gauss Radau pinned at x=-1,.
Definition: PointsType.h:59
RiemannSolverFactory & GetRiemannSolverFactory()
double NekDouble
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285