Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/Riemann/AUSMMultiSpecies.hpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.
7 : : All rights reserved. See the LICENSE file for details.
8 : : \brief Advection Upstream Splitting Method (AUSM+) Riemann flux function
9 : : for multi-species fluid dynamics.
10 : : \details This file implements the Advection Upstream Splitting Method (AUSM)
11 : : Riemann solver, with the all-speed corrections.
12 : : Ref. Liou, M. S. (2006). A sequel to AUSM, Part II: AUSM+-up for
13 : : all speeds. Journal of computational physics, 214(1), 137-170.
14 : : */
15 : : // *****************************************************************************
16 : : #ifndef AUSMMultiSpecies_h
17 : : #define AUSMMultiSpecies_h
18 : :
19 : : #include <vector>
20 : :
21 : : #include "Fields.hpp"
22 : : #include "SplitMachFns.hpp"
23 : : #include "FunctionPrototypes.hpp"
24 : : #include "Inciter/Options/Flux.hpp"
25 : : #include "EoS/EOS.hpp"
26 : : #include "MultiSpecies/MultiSpeciesIndexing.hpp"
27 : :
28 : : namespace inciter {
29 : :
30 : : //! AUSMMultiSpecies+up approximate Riemann solver
31 : : struct AUSMMultiSpecies {
32 : :
33 : : //! AUSM+up approximate Riemann solver flux function for multi-species flow
34 : : //! \param[in] fn Face/Surface normal
35 : : //! \param[in] u Left and right unknown/state vector
36 : : //! \return Riemann flux solution according to AUSM+up, appended by Riemann
37 : : //! velocities and volume-fractions.
38 : : //! \note The function signature must follow tk::RiemannFluxFn
39 : : static tk::RiemannFluxFn::result_type
40 : 0 : flux( const std::vector< EOS >& mat_blk,
41 : : const std::array< tk::real, 3 >& fn,
42 : : const std::array< std::vector< tk::real >, 2 >& u,
43 : : const std::vector< std::array< tk::real, 3 > >& = {} )
44 : : {
45 : 0 : auto nspec = g_inputdeck.get< tag::multispecies, tag::nspec >();
46 : 0 : auto k_p = g_inputdeck.get< tag::lowspeed_kp >();
47 : :
48 : 0 : auto ncomp = u[0].size();
49 [ - - ]: 0 : std::vector< tk::real > flx( ncomp, 0 );
50 : :
51 : 0 : tk::real rhol(0.0), rhor(0.0), pl(0.0), pr(0.0), hl(0.0), hr(0.0),
52 : 0 : al(0.0), ar(0.0), a12(0.0), rho12(0.0);
53 : :
54 : : // Mixture densities
55 [ - - ]: 0 : for (std::size_t k=0; k<nspec; ++k)
56 : : {
57 : 0 : rhol += u[0][multispecies::densityIdx(nspec, k)];
58 : 0 : rhor += u[1][multispecies::densityIdx(nspec, k)];
59 : : }
60 : :
61 : : // Velocities
62 : 0 : auto ul = u[0][multispecies::momentumIdx(nspec, 0)]/rhol;
63 : 0 : auto vl = u[0][multispecies::momentumIdx(nspec, 1)]/rhol;
64 : 0 : auto wl = u[0][multispecies::momentumIdx(nspec, 2)]/rhol;
65 : 0 : auto ur = u[1][multispecies::momentumIdx(nspec, 0)]/rhor;
66 : 0 : auto vr = u[1][multispecies::momentumIdx(nspec, 1)]/rhor;
67 : 0 : auto wr = u[1][multispecies::momentumIdx(nspec, 2)]/rhor;
68 : :
69 [ - - ]: 0 : pl = mat_blk[0].compute< EOS::pressure >( rhol, ul, vl, wl,
70 : 0 : u[0][multispecies::energyIdx(nspec, 0)] );
71 : 0 : hl = u[0][multispecies::energyIdx(nspec, 0)] + pl;
72 [ - - ]: 0 : al = mat_blk[0].compute< EOS::soundspeed >( rhol, pl );
73 : :
74 [ - - ]: 0 : pr = mat_blk[0].compute< EOS::pressure >( rhor, ur, vr, wr,
75 : 0 : u[1][multispecies::energyIdx(nspec, 0)] );
76 : 0 : hr = u[1][multispecies::energyIdx(nspec, 0)] + pr;
77 [ - - ]: 0 : ar = mat_blk[0].compute< EOS::soundspeed >( rhor, pr );
78 : :
79 : : // Average states for mixture speed of sound
80 : 0 : a12 = 0.5*(al+ar);
81 : 0 : rho12 = 0.5*(rhol+rhor);
82 : :
83 : : // Face-normal velocities from advective velocities
84 : 0 : auto vnl = ul*fn[0] + vl*fn[1] + wl*fn[2];
85 : 0 : auto vnr = ur*fn[0] + vr*fn[1] + wr*fn[2];
86 : :
87 : : // Mach numbers
88 : 0 : auto ml = vnl/a12;
89 : 0 : auto mr = vnr/a12;
90 : :
91 : : // All-speed parameters
92 : : // These parameters control the amount of all-speed diffusion necessary for
93 : : // low-Mach flows. Setting k_u and k_p to zero does not add any all-speed
94 : : // diffusion, whereas setting k_u and k_p to 1 adds maximum recommended
95 : : // all-speed diffusion. See "Liou, M. S. (2006). A sequel to AUSM, Part II:
96 : : // AUSM+-up for all speeds. Journal of computational physics, 214(1),
97 : : // 137-170" for more mathematical explanation. k_u is the velocity diffusion
98 : : // term and k_p is the pressure diffusion term. These two terms reduce
99 : : // pressure-velocity decoupling (chequerboarding/odd-even oscillations).
100 : 0 : tk::real k_u(1.0), f_a(1.0);
101 : :
102 : : // Split Mach polynomials
103 : 0 : auto msl = splitmach_ausm( f_a, ml );
104 : 0 : auto msr = splitmach_ausm( f_a, mr );
105 : :
106 : : // Riemann Mach number
107 : 0 : auto m0 = 1.0 - (0.5*(vnl*vnl + vnr*vnr)/(a12*a12));
108 : 0 : auto mp = -k_p* std::max(m0, 0.0) * (pr-pl) / (f_a*rho12*a12*a12);
109 : 0 : auto m12 = msl[0] + msr[1] + mp;
110 : 0 : auto vriem = a12 * m12;
111 : :
112 : : // Riemann pressure
113 : 0 : auto pu = -k_u* msl[2] * msr[3] * f_a * rho12 * a12 * (vnr-vnl);
114 : 0 : auto p12 = msl[2]*pl + msr[3]*pr + pu;
115 : :
116 : : // Flux vector splitting
117 : 0 : auto l_plus = 0.5 * (vriem + std::fabs(vriem));
118 : 0 : auto l_minus = 0.5 * (vriem - std::fabs(vriem));
119 : :
120 : : // Conservative fluxes
121 [ - - ]: 0 : for (std::size_t k=0; k<nspec; ++k)
122 : : {
123 : 0 : flx[multispecies::densityIdx(nspec, k)] =
124 : 0 : l_plus *u[0][multispecies::densityIdx(nspec, k)] +
125 : 0 : l_minus*u[1][multispecies::densityIdx(nspec, k)];
126 : : }
127 : :
128 : 0 : flx[multispecies::energyIdx(nspec, 0)] = l_plus*hl + l_minus*hr;
129 : :
130 [ - - ]: 0 : for (std::size_t idir=0; idir<3; ++idir)
131 : : {
132 : 0 : flx[multispecies::momentumIdx(nspec, idir)] =
133 : 0 : l_plus *u[0][multispecies::momentumIdx(nspec, idir)] +
134 : 0 : l_minus*u[1][multispecies::momentumIdx(nspec, idir)] + p12*fn[idir];
135 : : }
136 : :
137 : 0 : return flx;
138 : : }
139 : :
140 : : //! Flux type accessor
141 : : //! \return Flux type
142 : : static ctr::FluxType type() noexcept { return ctr::FluxType::AUSM; }
143 : : };
144 : :
145 : : } // inciter::
146 : :
147 : : #endif // AUSMMultiSpecies_h
|
