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 : : #include "MultiSpecies/Mixture/Mixture.hpp"
28 : :
29 : : namespace inciter {
30 : :
31 : : //! AUSMMultiSpecies+up approximate Riemann solver
32 : : struct AUSMMultiSpecies {
33 : :
34 : : //! AUSM+up approximate Riemann solver flux function for multi-species flow
35 : : //! \param[in] fn Face/Surface normal
36 : : //! \param[in] u Left and right unknown/state vector
37 : : //! \return Riemann flux solution according to AUSM+up, appended by Riemann
38 : : //! velocities and volume-fractions.
39 : : //! \note The function signature must follow tk::RiemannFluxFn
40 : : static tk::RiemannFluxFn::result_type
41 : 1752300 : flux( const std::vector< EOS >& mat_blk,
42 : : const std::array< tk::real, 3 >& fn,
43 : : const std::array< std::vector< tk::real >, 2 >& u,
44 : : const std::vector< std::array< tk::real, 3 > >& = {} )
45 : : {
46 : 1752300 : auto nspec = g_inputdeck.get< tag::multispecies, tag::nspec >();
47 : :
48 : : // All-speed parameters
49 : : // These parameters control the amount of all-speed diffusion necessary for
50 : : // low-Mach flows. Setting k_u and k_p to zero does not add any all-speed
51 : : // diffusion, whereas setting k_u and k_p to 1 adds maximum recommended
52 : : // all-speed diffusion. See "Liou, M. S. (2006). A sequel to AUSM, Part II:
53 : : // AUSM+-up for all speeds. Journal of computational physics, 214(1),
54 : : // 137-170" for more mathematical explanation. k_u is the velocity diffusion
55 : : // term and k_p is the pressure diffusion term. These two terms reduce
56 : : // pressure-velocity decoupling (chequerboarding/odd-even oscillations).
57 : 1752300 : auto k_u = g_inputdeck.get< tag::lowspeed_ku >();
58 : 1752300 : auto k_p = g_inputdeck.get< tag::lowspeed_kp >();
59 : :
60 : 1752300 : auto ncomp = u[0].size()-1;
61 [ + - ]: 1752300 : std::vector< tk::real > flx( ncomp, 0 );
62 : :
63 : : tk::real rhol(0.0), rhor(0.0), pl(0.0), pr(0.0), hl(0.0), hr(0.0),
64 : : Tl(0.0), Tr(0.0), al(0.0), ar(0.0), a12(0.0), rho12(0.0);
65 : :
66 : : // initialize mixtures
67 [ + - ]: 1752300 : Mixture mixl(nspec, u[0], mat_blk);
68 [ + - ]: 1752300 : Mixture mixr(nspec, u[1], mat_blk);
69 : :
70 : : // Mixture densities
71 : 1752300 : rhol = mixl.get_mix_density();
72 : 1752300 : rhor = mixr.get_mix_density();
73 [ + - ]: 1752300 : Tl = u[0][ncomp+multispecies::temperatureIdx(nspec, 0)];
74 : 1752300 : Tr = u[1][ncomp+multispecies::temperatureIdx(nspec, 0)];
75 : :
76 : : // Velocities
77 [ + - ]: 1752300 : auto ul = u[0][multispecies::momentumIdx(nspec, 0)]/rhol;
78 : 1752300 : auto vl = u[0][multispecies::momentumIdx(nspec, 1)]/rhol;
79 : 1752300 : auto wl = u[0][multispecies::momentumIdx(nspec, 2)]/rhol;
80 : 1752300 : auto ur = u[1][multispecies::momentumIdx(nspec, 0)]/rhor;
81 : 1752300 : auto vr = u[1][multispecies::momentumIdx(nspec, 1)]/rhor;
82 : 1752300 : auto wr = u[1][multispecies::momentumIdx(nspec, 2)]/rhor;
83 : :
84 [ + - ]: 1752300 : pl = mixl.pressure( rhol, Tl );
85 [ + - ]: 1752300 : hl = u[0][multispecies::energyIdx(nspec, 0)] + pl;
86 [ + - ]: 1752300 : al = mixl.frozen_soundspeed( rhol, Tl, mat_blk );
87 : :
88 [ + - ]: 1752300 : pr = mixr.pressure( rhor, Tr );
89 [ + - ]: 1752300 : hr = u[1][multispecies::energyIdx(nspec, 0)] + pr;
90 [ + - ]: 1752300 : ar = mixr.frozen_soundspeed( rhor, Tr, mat_blk );
91 : :
92 : : // Average states for mixture speed of sound
93 : 1752300 : a12 = 0.5*(al+ar);
94 : 1752300 : rho12 = 0.5*(rhol+rhor);
95 : :
96 : : // Face-normal velocities from advective velocities
97 : 1752300 : auto vnl = ul*fn[0] + vl*fn[1] + wl*fn[2];
98 : 1752300 : auto vnr = ur*fn[0] + vr*fn[1] + wr*fn[2];
99 : :
100 : : // Mach numbers
101 : 1752300 : auto ml = vnl/a12;
102 : 1752300 : auto mr = vnr/a12;
103 : :
104 : : tk::real f_a(1.0);
105 : :
106 : : // Split Mach polynomials
107 : 1752300 : auto msl = splitmach_ausm( ml, f_a );
108 : 1752300 : auto msr = splitmach_ausm( mr, f_a );
109 : :
110 : : // Riemann Mach number
111 : 1752300 : auto m0 = 1.0 - (0.5*(vnl*vnl + vnr*vnr)/(a12*a12));
112 [ + + ]: 1752300 : auto mp = -k_p* std::max(m0, 0.0) * (pr-pl) / (f_a*rho12*a12*a12);
113 : 1752300 : auto m12 = msl[0] + msr[1] + mp;
114 : 1752300 : auto vriem = a12 * m12;
115 : :
116 : : // Riemann pressure
117 : 1752300 : auto pu = -k_u* msl[2] * msr[3] * f_a * rho12 * a12 * (vnr-vnl);
118 : 1752300 : auto p12 = msl[2]*pl + msr[3]*pr + pu;
119 : :
120 : : // Flux vector splitting
121 : 1752300 : auto l_plus = 0.5 * (vriem + std::fabs(vriem));
122 : 1752300 : auto l_minus = 0.5 * (vriem - std::fabs(vriem));
123 : :
124 : : // Conservative fluxes
125 [ + + ]: 4380750 : for (std::size_t k=0; k<nspec; ++k)
126 : : {
127 : 2628450 : flx[multispecies::densityIdx(nspec, k)] =
128 : 2628450 : l_plus *u[0][multispecies::densityIdx(nspec, k)] +
129 : 2628450 : l_minus*u[1][multispecies::densityIdx(nspec, k)];
130 : : }
131 : :
132 : 1752300 : flx[multispecies::energyIdx(nspec, 0)] = l_plus*hl + l_minus*hr;
133 : :
134 [ + + ]: 7009200 : for (std::size_t idir=0; idir<3; ++idir)
135 : : {
136 : 5256900 : flx[multispecies::momentumIdx(nspec, idir)] =
137 : 5256900 : l_plus *u[0][multispecies::momentumIdx(nspec, idir)] +
138 : 5256900 : l_minus*u[1][multispecies::momentumIdx(nspec, idir)] + p12*fn[idir];
139 : : }
140 : :
141 : 1752300 : return flx;
142 : : }
143 : :
144 : : //! Flux type accessor
145 : : //! \return Flux type
146 : : static ctr::FluxType type() noexcept { return ctr::FluxType::AUSM; }
147 : : };
148 : :
149 : : } // inciter::
150 : :
151 : : #endif // AUSMMultiSpecies_h
|
