Nektar++
AUSM0Solver.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: AUSM0Solver.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: AUSM0 Riemann solver.
32//
33///////////////////////////////////////////////////////////////////////////////
34
36
37namespace Nektar
38{
39
40std::string AUSM0Solver::solverName =
42 "AUSM0", AUSM0Solver::create, "AUSM0 Riemann solver");
43
45 : CompressibleSolver(pSession)
46{
47}
48
49/**
50 * @brief AUSM0 Riemann solver
51 *
52 * @param rhoL Density left state.
53 * @param rhoR Density right state.
54 * @param rhouL x-momentum component left state.
55 * @param rhouR x-momentum component right state.
56 * @param rhovL y-momentum component left state.
57 * @param rhovR y-momentum component right state.
58 * @param rhowL z-momentum component left state.
59 * @param rhowR z-momentum component right state.
60 * @param EL Energy left state.
61 * @param ER Energy right state.
62 * @param rhof Computed Riemann flux for density.
63 * @param rhouf Computed Riemann flux for x-momentum component
64 * @param rhovf Computed Riemann flux for y-momentum component
65 * @param rhowf Computed Riemann flux for z-momentum component
66 * @param Ef Computed Riemann flux for energy.
67 */
68void AUSM0Solver::v_PointSolve(double rhoL, double rhouL, double rhovL,
69 double rhowL, double EL, double rhoR,
70 double rhouR, double rhovR, double rhowR,
71 double ER, double &rhof, double &rhouf,
72 double &rhovf, double &rhowf, double &Ef)
73{
74 // Left and Right velocities
75 NekDouble uL = rhouL / rhoL;
76 NekDouble vL = rhovL / rhoL;
77 NekDouble wL = rhowL / rhoL;
78 NekDouble uR = rhouR / rhoR;
79 NekDouble vR = rhovR / rhoR;
80 NekDouble wR = rhowR / rhoR;
81
82 // Internal energy (per unit mass)
83 NekDouble eL = (EL - 0.5 * (rhouL * uL + rhovL * vL + rhowL * wL)) / rhoL;
84 NekDouble eR = (ER - 0.5 * (rhouR * uR + rhovR * vR + rhowR * wR)) / rhoR;
85 // Pressure
86 NekDouble pL = m_eos->GetPressure(rhoL, eL);
87 NekDouble pR = m_eos->GetPressure(rhoR, eR);
88 // Speed of sound
89 NekDouble cL = m_eos->GetSoundSpeed(rhoL, eL);
90 NekDouble cR = m_eos->GetSoundSpeed(rhoR, eR);
91
92 // Average speeds of sound
93 NekDouble cA = 0.5 * (cL + cR);
94
95 // Local Mach numbers
96 NekDouble ML = uL / cA;
97 NekDouble MR = uR / cA;
98
99 // Parameters for specify the upwinding
100 NekDouble beta = 0.0;
101 NekDouble alpha = 0.0;
102 NekDouble Mbar = M4Function(0, beta, ML) + M4Function(1, beta, MR);
103 NekDouble pbar =
104 pL * P5Function(0, alpha, ML) + pR * P5Function(1, alpha, MR);
105
106 if (Mbar >= 0.0)
107 {
108 rhof = cA * Mbar * rhoL;
109 rhouf = cA * Mbar * rhoL * uL + pbar;
110 rhovf = cA * Mbar * rhoL * vL;
111 rhowf = cA * Mbar * rhoL * wL;
112 Ef = cA * Mbar * (EL + pL);
113 }
114 else
115 {
116 rhof = cA * Mbar * rhoR;
117 rhouf = cA * Mbar * rhoR * uR + pbar;
118 rhovf = cA * Mbar * rhoR * vR;
119 rhowf = cA * Mbar * rhoR * wR;
120 Ef = cA * Mbar * (ER + pR);
121 }
122}
123
124/**
125 *
126 */
127double AUSM0Solver::M1Function(int A, double M)
128{
129 double out;
130
131 if (A == 0)
132 {
133 out = 0.5 * (M + fabs(M));
134 }
135 else
136 {
137 out = 0.5 * (M - fabs(M));
138 }
139
140 return out;
141}
142
143double AUSM0Solver::M2Function(int A, double M)
144{
145 double out;
146
147 if (A == 0)
148 {
149 out = 0.25 * (M + 1.0) * (M + 1.0);
150 }
151 else
152 {
153 out = -0.25 * (M - 1.0) * (M - 1.0);
154 }
155
156 return out;
157}
158
159/**
160 *
161 */
162double AUSM0Solver::M4Function(int A, double beta, double M)
163{
164 double out;
165
166 if (fabs(M) >= 1.0)
167 {
168 out = M1Function(A, M);
169 }
170 else
171 {
172 out = M2Function(A, M);
173
174 if (A == 0)
175 {
176 out *= 1.0 - 16.0 * beta * M2Function(1, M);
177 }
178 else
179 {
180 out *= 1.0 + 16.0 * beta * M2Function(0, M);
181 }
182 }
183
184 return out;
185}
186
187/**
188 *
189 */
190double AUSM0Solver::P5Function(int A, double alpha, double M)
191{
192 double out;
193
194 if (fabs(M) >= 1.0)
195 {
196 out = (1.0 / M) * M1Function(A, M);
197 }
198 else
199 {
200 out = M2Function(A, M);
201
202 if (A == 0)
203 {
204 out *= (2.0 - M) - 16.0 * alpha * M * M2Function(1, M);
205 }
206 else
207 {
208 out *= (-2.0 - M) + 16.0 * alpha * M * M2Function(0, M);
209 }
210 }
211
212 return out;
213}
214
215} // namespace Nektar
double M2Function(int A, double M)
static RiemannSolverSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession)
Definition: AUSM0Solver.h:46
static std::string solverName
Definition: AUSM0Solver.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
AUSM0 Riemann solver.
Definition: AUSM0Solver.cpp:68
double P5Function(int A, double alpha, double M)
double M1Function(int A, double M)
double M4Function(int A, double beta, double M)
AUSM0Solver(const LibUtilities::SessionReaderSharedPtr &pSession)
Definition: AUSM0Solver.cpp:44
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