On the Development of an Implicit Discontinuous Galerkin Solver for Turbulent Real Gas Flows

: The aim of this work is to describe an efficient implementation of cubic and multiparameter real gas models in an existing discontinuous Galerkin solver to extend its capabilities to the simulation of turbulent real gas flows. The adopted thermodynamic models are van der Waals, Peng–Robinson


Introduction
Nowadays, the increasing interest of industries towards highly accurate simulation tools motivates the implementation of complex physical models and numerical schemes to reproduce specific phenomena.In this context, non ideal compressible fluid dynamics (NICFD) is still quite a challenging task, mainly because the determination of correct thermophysical properties is crucial to obtain accurate and robust solvers, but also because non-classical behaviors may arise in these flows.Examples can be found in turbomachinery for organic Rankine cycles (ORC), carbon capture and storage (CCS), and refrigeration systems.
During the last decades, many real gas models have been proposed to overcome the limits of the perfect gas law.Actually, the most accurate models are the multi-parameter Helmholtz energy equations of state (MEoSs).Simpler models are available, generally written in terms of cubic polynomials of the density.These cubic equations of state (CEoSs) are widely used, and, sometimes, preferred to MEoSs for the easiness of implementation and use.CEoSs have simpler formulations, require a very limited number of fluid parameters, and their computational cost is an order of magnitude lower than MEoSs.However, for some particular problems or with highly accurate solvers, the adoption of ad hoc models is desirable to obtain accurate predictions.
The aim of the present work is the development of a highly accurate discontinuous Galerkin (dG) solver for the simulation of turbulent real gas flows, where the higher accuracy guaranteed by dG methods is coupled with reliable methods for the calculation of thermodynamic properties.In particular, different thermal EoSs are implemented: (i) the pressure-explicit van der Waals (vdW) [1] and Peng-Robinson (PR) [2] CEoSs, and (ii) the Helmholtz-explicit MEoS of Span-Wagner (SW) [3].
The test cases chosen to asses the new solver performances are (i) an unsteady shock tube problem [4], (ii) a stationary supersonic wedge-shaped channel [5], and (iii) an ORC turbine blade nozzle [6].Chung et al. model [7] is adopted to compute transport properties.
The implementation of a non ideal EoS in a dG solver requires also the modification of the algorithm to compute the convective numerical flux and the boundary conditions.The dG finite element method (FEM) provides by definition a solution that is discontinuous across elements interfaces in the grid, so a unique value for the convective flux must be determined to guarantee the conservation and the stability of the numerical scheme, just as in finite volume methods (FVM).Many procedures for the computation of the convective numerical flux are available, based on the exact or approximate solution of a Riemann problem, but they assume an ideal behavior of the flow.As a consequence, a thermodynamic generalization is required, and the extension to real gas flows of the Roe's linearization for the Riemann problem [8] proposed by Vinokur and Montagné [9] is adopted (the average speed of sound at the interfaces is computed according to Glaister [10]).Moreover, a generalized set of boundary conditions has to be determined, especially for inflow/outflow boundaries, which are normally based on the theory of the Riemann invariants.The extension of the Riemann invariants to real gas models is quite complex, and, for this reason, the linearization proposed by Colonna et al. [11] is employed in this work to solve the boundary problem in a consistent and generalized way.
The proposed implementation has been used to extend the prediction capability of the dG-FEM solver MIGALE [12][13][14], whose performance has been already assessed for turbulent flows with ideal behavior.The solver adopts an implicit time integration strategy, and, as a consequence, at each iteration the Jacobian matrix must be computed.In this work, the automatic differentiation (AD) tool Tapenade [15] is used to derive the exact Jacobian matrix to keep the solver able to reach the quadratic convergence speed on stationary problems, which is proper of the Newton-type method.AD has been also employed to derive some complex thermodynamic derivatives.
This paper is organized as follows.First, a brief description of the dG-FEM solver MIGALE is presented from the spatial and temporal points of view (Section 2).Then, all the details of implementation of the real gas models are discussed with a focus on the implementation of the auxiliary procedures (Section 3).After that, the results obtained from the validation test cases are discussed (Section 4) and, in the end, some conclusions are presented (Section 5).

Discontinuous Galerkin Solver
In this section, the main features of the dG-solver MIGALE [12][13][14] are outlined.In particular, the governing equations are presented in Section 2.1 and the discretization methods, in space (Section 2.2) and time (Section 2.3), are discussed, with a particular emphasis on the aspect that will be crucial for the implementation of the new thermodynamic models.

Governing Equations
The set of conservation laws to be solved for inviscid flows is given by the Euler equations, whereas for turbulent flows by the Reynolds-averaged Navier-Stokes (RANS) equations, here supplemented by the k − ω turbulence model [12][13][14].The RANS and k − ω model equations can be written using Einstein's notation as where u i is the flow velocity, p the pressure, ρ the density, and are the total mass-specific internal energy and enthalpy, the turbulent and the overall shear stress tensors, and the overall Fourier's conductive heat flux, calculated using both a molecular and a turbulent thermal conductivity λ and λ t = (µ t c p )/(Pr t ), where Pr t is the turbulent Prandtl number.The remaining quantities are the mean strain-rate tensor, i.e., the turbulent dynamic viscosity µ t = α * ρke − ω r , the limited mass-specific turbulent kinetic energy k = max(0, k) and the logarithm of the specific dissipation rate ω = log(ω).α, α * , β, β * , σ, σ * are the closure parameters [16].The production term of the energy equation and the destruction term of the k and ω equations are computed with the value ω r , which satisfies the realizability condition for the turbulent stresses [12].In this work, no additional terms are added in the turbulence model equations to account for compressibility effects.
In fact, this treatment should be considered just for hypersonic flows, i.e., M a > 5, with cold walls, which are not present in the proposed testcases.A brief review of the possible pressure corrections can be found in [17].
Equations ( 1)-( 5) can then be written in the following compact form as where w is the vector of the unknown variables, F c and F v are, respectively, the convective and the viscous flux, and s is the vector of the source terms.The matrix P(w) takes into account the change of variables from the conservative to the primitive set w = p, u 1 , u 2 , u 3 , T, k, ω T , where p = log(p), T = log(T), and ω = log(ω) are used to enhance the solver's robustness [12].

Spatial Discretization
The weak formulation of the problem is obtained by multiplying Equation ( 10) by an arbitrary smooth test function v = {v 1 , . . ., v m } and integrating by parts over a physical domain Ω, with m being the total number of unknowns.Once a proper triangulation T h of the approximated domain Ω h in arbitrary shaped non-overlapping elements having the set of faces F h is given, the discrete weak problem is obtained by substituting the continuous solution w and the continuous test function v with their discrete finite element approximations w h and v h , each one belonging to the discrete polynomial space V h = [P q d (T h )] m expressed in physical coordinates.
The set of test and shape functions in every element K ∈ T h is chosen here as the set {φ} of N K do f orthonormal and hierarchical basis functions having compact support over K, defined from its principal inertial axes and with N K do f being the total number of degrees of freedom of the solution in K.Each component of the discrete elemental solutions w h,j , with j = 1, . . ., m, can then be expressed as the linear combination w h,j = φ l W j,l , with l = 1, . . ., N K do f and ∀K ∈ T h .The direct sum of all the discrete elemental solutions over T h represents the global discrete solution, which is the objective of the solver.The dG-FEM spatial discretization of the governing equations consists therefore in seek, at every time instant, the elements of the degrees of freedom's vector W ∈ R n e ×m×N K do f with n e being the total number of elements, such that for i = 1, . . ., N K do f and k = 1, . . ., n e , they represent the solution of the semi-discrete weak problem As the functional approximation is discontinuous, the sum of the convective and viscous flux functions F is not uniquely defined at each element's interface, so a numerical flux vector F is adopted.The convective part is based on the local solution of linearized Riemann problems, using the Roe solver [8] generalized to the case of an arbitrary gas model with the Vinokur-Montagné approach [9] and Glaister's [10] generalized average speed of sound.The viscous part is instead centered and discretized with the BR2 scheme [18], by employing either local and global lifting operators r F (•) and r(•) on the solution's componentwise jump [[w h ]] = w h|K + n F + + w h|K − n F − across mesh interfaces.In this sense, n F ± denotes the outward or inward pointing unit vector normal to the interface, whereas η F is the stability parameter of Brezzi et al. [19].
To avoid spurious oscillations of the solution, an artificial diffusion contribution is introduced inside each element using a shock sensor to detect discontinuities.The shockcapturing term SCT is added to the left-hand side of Equation (11), and, as reported in [14,20], is given by where b(w) is a unit vector representing the direction along which the dissipation is acting.In this work, b(w) is given by the logarithmic pressure gradient, where ε is a small value proportional to the machine precision.In Equation ( 12), p is the artificial diffusion coefficient defined in each element as where C is a user-defined constant, h K is a characteristic dimension of the element K, and s p w ± h , w h is a pressure-based shock sensor defined through the global lifting operator as s p w ± h , w h is always active in every element, but the numerical viscosity is introduced only in regions where unphysical oscillations are present.The remaining terms d p (w h ) and f p (w h ) in Equation ( 14) introduce also the dependence of the numerical viscosity from the magnitude of the divergence of the convective flux and from the polynomial degree of discretization, respectively.

Temporal Discretization
Assembling the elemental contribution of Equation ( 11), the following system of ordinary differential Equations (ODEs) in time is obtained where R(W) is the vector of the global residuals and M P (W) is the global block diagonal mass matrix arising from the calculation of the first integral in Equation ( 11).The Linearized Backward Euler (LBE) scheme with a pseudo-transient continuation strategy for stationary problems is adopted to solve Equation ( 16) [21], which can be written as When the final solution is steady, an exponential CFL law, the function of the residuals norms, enables the usage of progressively higher values of ∆t that reduce Equation ( 17) to a Newton-Rhapson method, which guarantees quadratic convergence rates, once an exact Jacobian matrix ∂R(W n )/∂W is provided at every timestep n.The algebraic system described by Equation ( 17) is nonlinear, and an iterative solver is required ∀n.In this work, a restarted version of the generalized minimal residual (GMRES) Krylov's subspace-type method is used, as available in the PETSc library [22].These kinds of methods have been extensively used and developed during the last decades [23] for their generality and robustness.They are still the subject of intense research activity to improve their convergence speeds through techniques such as globalization [24], efficient preconditioning [25], and Jacobian approximation [26].Here, the GMRES convergence is enhanced by system preconditioning; MIGALE allows us to choose the block Jacobi method with one block per process, each of which is solved with ILU(0), or the additive Schwarz method (ASM), as available in the PETSc library.The ASM [27] is used for the simulation presented in this work.
Devising an effective and robust strategy to increase the CFL number as the residual decreases is far from trivial, especially for transitional or turbulent simulations, an empirically determined "CFL law" is here used to speed up convergence.It is based on the L ∞ and L 2 norms of the residual and depends on three user-defined parameters.The first and second ones are CFL min and CFL max to set the minimum and maximum limits of the CFL number during the simulation.The third one is an exponent α governing the growth rate of the CFL number, where typically α ≤ 1.The "CFL law" is where CFL exp = min(1/(2q + 1), CFL min ) is the minimum value between the maximum CFL number proper of an explicit scheme and the user-defined minimum value, with q being the polynomial degree of discretization of the solution.The remaining terms are β = CFL min − CFL exp and ξ, which is defined as where || • || 2 and || • || ∞ are the L 2 and L ∞ norms of the residual vector of the i-th equation of the system R i and R i 0 is the corresponding residual at the first iteration.

Thermodynamic Models
The MIGALE solver's predicting capabilities are extended with three real gas models: the pressure-explicit CEoSs of van der Waals [1] and Peng-Robinson [2] (Section 3.1) and the Helmholtz-explicit MEoS of Span-Wagner [3] (Section 3.2).MEoSs generally require more coefficients with respect to CEoSs and their computational cost is, therefore, much higher.On the other hand, they also guarantee superior accuracy for thermodynamic quantities whose values are crucial during a fluid dynamic simulation, such as the speed of sound.The implementation of all models is discussed, also describing modifications needed by other algorithms of the solver, such as the numerical flux computation and boundary conditions (Sections 3.3 and 3.4).

Peng-Robinson and van der Waals Models
The pressure-explicit CEoSs of van der Waals [1] and Peng-Robinson [2] can be obtained from the general formulation [28] as where p is the fluid pressure, ρ the density, T the temperature, R * = R/m M the mass-specific gas constant, R = 8314.463J/(kmolK) the universal gas constant, and m M the fluid's molecular weight.Equation (20) shows also the term A(T), which accounts for intermolecular attractions, and the terms B, C, and D that account for molecular volume.The term A is usually written as A = aα 2 (T), where the function α(T) (if not null) contains the dependence of A(T) from the molecular shape, whereas a is a constant.For all models that can be obtained from Equation (20), A, B, C, and D assume different values, depending on the working fluid.In fact, they depend from some input parameters, which are the critical pressure p cr and temperature T cr , the molecular weight, and the acentric factor ω, which is an estimation of the non-sphericity of the molecules defined as ω = (− log 10 (p sat r ) − 1)| T r =0.7 , where p sat r = p sat /p cr , T r = T/T cr , and p sat (T) is the saturation pressure.Table 1 summarizes the expressions that must be used for A, B, C, and D to obtain the van der Waals and the Peng-Robinson gas models, which are given by Starting from these equations, a complete characterization of a pure single-phase substance comes from the determination of at least one caloric EoS for each model [29].A general procedure for any thermal pressure-explicit EoS like the one in Equation ( 20) is given by Reynolds [30].The expression for the mass-specific internal energy takes the form whereas the mass-specific entropy is where ξ and η are used as symbolic substitutes of ρ and T in the integral functions and (ρ 0 , T 0 ) identifies an arbitrary reference state.The last terms in both Equations ( 22) and ( 23) represent departure functions from the non-polytropic ideal gas behavior since they vanish for sufficiently rarefied thermodynamic states, i.e., ρ → 0. The remaining two integrals require instead an expression for the ideal gas contribution to the isochoric specific heat c 0 v (T), which is by definition the limit of c v (ρ, T) as ρ → 0. In this work, a polynomial function of the absolute temperature in the form employed for each considered fluid, where c 0 p (T) is the ideal gas contribution to the isobaric specific heat.Coefficients c i for i = 0, . . ., 3 can be determined theoretically from chemical group contribution methods such as the one in [31], or from given polynomial fittings of experimental data available in the literature.
Table 1.Expressions for all the quantities involved in Equation (20).Once the expressions for p(ρ, T) and e(ρ, T) are known, all the other relevant thermodynamic properties can be determined using a combination of them and of their derivatives.For example, by definition the mass-specific enthalpy and real gas isochoric specific heat are obtained as

van der Waals Peng-Robinson
As reported in [32], the real gas isobaric specific heat and the speed of sound are obtained from Another important quantity that must be determined is the fundamental derivative of gas dynamics Γ, that following the work of Cramer [33] can be again expressed as a function of temperature and density only, as This derivative is crucial in real gas dynamics, since with negative values of Γ some non-classical phenomena may arise, such as expansion shocks or compression fans [34].
Lastly, since the solver works with p and T as independent variables, the computation of the inverse problem is needed, as ρ = ρ(p, T), T = T(ρ, p), T = T(ρ, e). ( The fluid density is determined from the equation for the pressure.When the models are derived by Equation ( 20), the thermal EoS can be reformulated as a third-degree polynomial in the density, whose coefficients are a function of temperature and pressure, i.e., d 0 + d 1 ρ + d 2 ρ 2 + d 3 ρ 3 = 0, with d i = d i (p, T) for i = 0, . . ., 3. The analytical resolution method of Cardano is employed in this work, whereas the physical meaning and validity of each root have been determined using the considerations in [35].For temperature, some Newton's iterations are employed on the functions p(ρ, T) and e(ρ, T), since their derivatives are known and the resulting formulation is more complicated.Initial guesses are calculated using the polytropic ideal gas model with γ = c 0 p (T 0 )/c 0 v (T 0 ).

Span-Wagner Model
The Helmholtz-explicit MEoS of Span-Wagner [3] is formulated in terms of an optimized functional fit of experimental measurements, which can be derived for any fluid having a sufficiently wide and precise range of data [36].The derived EoS is formulated for the free Helmholtz energy state function a(ρ, T) = e(ρ, T) − Ts(ρ, T), described in a non-dimensional form with the summation of an ideal gas contribution and a real gas residual as where δ = ρ/ρ cr is the reduced density and τ = T cr /T is the inverse of the reduced temperature.In Equation ( 28), the dimensional ideal gas part is defined as since for the ideal gas p/ρ = R * T. So, once a suitable approximation of c 0 p (T) is provided, a 0 (ρ, T) can be completely determined by computing two integrals.In this work, four different functional forms can be activated by the user, since c 0 p (T) is implemented as where each term represents an approximation of a statistical mechanical behavior of the ideal gas heat capacity as suggested by Aly and Lee [37].In Equation (30), coefficients (c 1,i , c 2,i ) are considered as user parameter, since many functional fittings can be found in the literature.The non-dimensional residual part of Equation ( 28) is similarly provided as a summation of various activatable terms as with 6,i and where the last Gaussian bell-shaped sums are generally used to improve the fluid description near the critical point [36].
Thanks to the Helmholtz energy definition, all the other relevant thermodynamic properties can be computed with Maxwell's relations, such as whereas Equations ( 24)-( 26) still hold for the calculation of the enthalpy, specific heats, speed of sound, and fundamental derivative.The inverse problem of Equation ( 27) is here treated with Newton's iterations also for the density, but since the number of roots may be higher than the CEoS case, some efficient initial guesses are chosen as suggested by Span [36].In particular, the initial guess for the density is provided by the Peng-Robinson model, whose coefficients are calculated and stored once.For the temperature, a simplified version of the van der Waals model with a power law ideal gas-specific heat is analytically inverted.The adopted expression is ]/ log(T 2 /T 1 ) and T 1 < T 0 < T 2 as suggested by [33].Furthermore, the van der Waals coefficients are calculated and stored before computations.

Derivatives
The first and second derivatives of the thermodynamic properties are needed for the Jacobian matrix of the implicit time integration scheme, the shock-capturing term, the permutation matrix, and the convective fluxes.In particular, the following derivatives must be provided: where x(p, T) can represent e, h, c, s, c v , c p .Since all the properties are formulated as functions of ρ and T, the exact expressions of their first and second derivatives with respect to these variables are obtained with the AD tool Tapenade [15].Then, using the relations from [38], which involve just the derivatives of p = p(ρ, T), the values of are calculated.Thanks to the chain rule on x = x[ρ(p, T), T], and considering that ∂/∂ ỹ = (∂/∂y)(∂ ỹ/∂y) −1 = y(∂/∂y), where y can be either p or T, the last five derivatives in Equation ( 33) can be rewritten as The first two derivatives in Equation ( 33) are determined as suggested by Cinnella [39].However, despite Equations ( 35)-( 39) being valid for all the chosen EoSs, the Span-Wagner model requires a further step, i.e., the computation of all the pure and mixed derivatives, from the first to the third order of the non-dimensional Helmholtz energy state function.This task is here performed with the AD tool Tapenade [15].

Numerical Fluxes and Boundary Conditions
The first thermodynamic generalization required by the solver is the adoption of a consistent numerical flux for real gas computations.In this work, the generalization of the approximate Riemann solver of Roe [8] proposed by Vinokur-Montagné [9] is used for the convective part.This procedure differs from the original Roe version since in the real gas regime the description of the Roe averaged state must be enriched with the definition of averaged values of the pressure derivatives χ = (∂p/∂ρ) e and κ = (∂p/∂e) ρ between the two sides of every mesh interface.These values are here obtained following the procedure proposed by Glaister [10] and then used to generalize the Roe averaged a speed of sound for the determination of convective eigenvalues.
For the viscous part, the generalized multiparameter correlation of Chung et al. [7] for the determination of transport properties in real gas regime is applied.In particular, the procedure allows for estimating reliable values of the molecular dynamic viscosity and thermal conductivity of polar and non-polar fluids as functions of ρ and T. The required additional input data are the critical density ρ cr , the dipole moment of the fluid molecules, and the equilibrium dissociation constant of the substance.
The contribution of the new flux to the global Jacobian matrix is computed with the AD tool Tapenade [15], an open source algorithm developed by the Institut National de Recherche en Sciences et Technologies du Numérique (INRIA).AD guarantees that every derivative will be mathematically exact and will not suffer any truncation error, which is typical of the finite differences (FD) approach [40].In fact, every derivative is obtained with a symbolic optimized differentiation of all the lines of a source code, to generate a new program that will contain the calculations for both the original outputs and their derivatives.This is made possible by an iterative application of the chain rule of differentiation since the whole source code is interpreted as a composite function of all its lines.The chain can be traveled from top to bottom with the tangent (or direct) differentiation mode or from bottom to top with the adjoint (or reverse) mode (see [15] for further details).In this work, the tangent mode has been used, since it is best suited for large amounts of inputs and is easy to use.In particular, the focus is on the term F = F(w ± , ∇w ± ) in Equation ( 11), where w − and w + are the unknown variables at the inner and outer side of an element face.The Jacobian matrix of F is generated column by column, differentiating F one time in tangent mode for every component of w ± and ∇w ± .Every column is then wrapped with the others to assemble the Jacobian matrix.This often results in an increment of the computational cost with respect to manually derived analytical procedures, which are often difficult to obtain.In this work, an ad hoc automated strategy for the use of Tapenade has been derived, that is able to scan and modify the generated routines to avoid or regroup redundant computations.The new Jacobian matrix is thus characterized by a lower computational cost with respect to FD, especially when the thermodynamic or physical complexity is high.Table 2 reports the time required to perform 10 5 calls to the routines to build the Jacobian matrix of an inviscid two-dimensional convective flux and a three-dimensional turbulent diffusive flux.AD is always less expansive than the FD counterpart, and shows a maximum reduction in the computing time ≈ 60%.The second generalization concerns inflow and outflow boundary conditions, which are implemented following the work of Colonna and Rebay [11].The approach relies on the determination of a linearized form of the Riemann invariants, that allows the imposition of the proper set of physical quantities at every boundary face in both subsonic and supersonic regimes, for both incoming and outgoing flows.The contribution to the global Jacobian matrix of the residual is also here derived with Tapenade [15], following the same approach described for the convective and viscous fluxes.

Results
In this section, the results obtained with the new solver are discussed and the predicted solutions are compared with available experimental and numerical reference data.

Unsteady Shock Tube
The first case consists of the Euler solution of a one-dimensional Riemann problem proposed by Guardone et al. [4] for the fluid PP10 (C 13 F 22 ).The aim of the original setup is to reproduce and capture a non-classical expansion shock at a precise time instant and spatial location.The length of the tube is 5 m, which is divided into 400 uniform elements.At x = 3 m, a diaphragm separates two regions, where the fluid is at different densities and pressure, as reported in Table 3.At time t = 0 s the diaphragm is removed and a rightrunning compression shock wave, a contact discontinuity, and a left-running rarefaction shock wave start to travel along the tube.Peng-Robinson EoS is used against the Martin-Hou (MH) model employed by Guardone et al. [4] and a second-order approximation is adopted for the spatial discretization.
The simulation is stopped at t = 29.46 × 10 −3 s and Figure 1 shows the density and pressure profiles.The initial condition resulted in an almost isothermal domain and onedimensional wave propagation is well captured by the solver.In particular, a rarefaction shock is observed, as expected, and oscillations of the solution in the neighborhood of discontinuities are kept small thanks to the shock-capturing term.In this case, the polynomial degrees and the number of mesh elements do not influence the predicted result, whereas a sufficiently low value for the maximum "CFL" number is mandatory to achieve a satisfactory accuracy in time, since the adopted time integration scheme, i.e., the LBE scheme, has a higher truncation error.

Supersonic Wedge
The second case consists of a supersonic inviscid flow of supercritical gaseous MDM (C 8 H 24 O 2 Si 3 ) in a wedge-shaped channel.The free stream Mach number at the entry of the domain is 1.7, whereas pressure and density are 15 bar and 202.888 kg/m 3 , respectively.The temperature is 571.72 K, calculated in this condition from Equation (21).The original problem is used by Pini et al. [5] to assess the performances of the code SU2 using the van der Waals EoS, which allows for predicting a negative Γ zone for MDM where nonclassical phenomena are possible.Different simulations have been performed for different polynomial degrees and a number of mesh elements: (i) 42 400 elements and P 0 solution approximation, (ii) 2 650 elements and P 4 solution approximation, and (iii) 42 400 elements and P 2 solution approximation.The total number of degrees of freedom for every solution can be calculated in this case as where q is the polynomial degree of the discretization, m the number of unknowns, and n e the number of elements.The angle formed by the rarefaction shock is reported together with the solution parameters in Table 4, whereas the value of the same angle from the theoretical relation reported by Pini et al. [5] is where subscripts u and d refer to the quantities upstream and downstream from the shock, respectively.Furthermore, Figure 2 shows the Mach number contours obtained for each solution.As expected, a reasonable accuracy is obtained even on very coarse grids, thanks to the high-order dG spatial discretization.Figure 2 shows as increasing the spatial resolution, both in terms of mesh density and polynomial approximation, the discontinuity given by the shock can be confined in a narrower stripe of elements.Furthermore, for this case, the non-classical behavior is well captured.

ORC Turbine Nozzle
The last problem investigated is the turbulent flow through an ORC turbine nozzle (the geometry has been provided by Turboden) with the Siloxane MDM as working fluid.The operating condition is characterized by a total inlet pressure of 8 bar, a total inlet temperature of 270.5 • C, and a total to static pressure ratio of 6.An inviscid two-dimensional solver was used by Colonna et al. [6], with a Span-Wagner-type MEoS given in the form of Equation ( 28), whose coefficients are given in [41].The same thermodynamic model is used in this work, but the RANS equation coupled with the kω turbulence model are solved on a quasi 3D domain, as the mesh is two-dimensional and it is extruded on one element in the third direction.The grid is unstructured with 5305 elements, as shown in Figure 3.At the inflow, the total temperature, the total pressure, the flow angle α 1 = 0 • , and the turbulence intensity Tu 1 = 4.0% are prescribed, whereas at the outflow, the static pressure is set.Blade wall is considered adiabatic.
All computations are performed for a P 2 solution approximation and a convergence tolerance of 10 −10 on the residuals norm is reached for every equation.All the computations have been run in parallel, initialising the P 0 solution from the uniform flow at outflow conditions and the higher order solutions from the lower order ones.Figure 4 shows the convergence history of the simulation with PR EoS.  Figure 5 shows the Mach number and compressibility factor z = p/(ρR * T) contours of a mesh section for the solution employing SW model.Some spurious reflections appear at the outflow boundary and pollute a bit the solution for the lack of non-reflective boundary conditions.However, the main structures of the flow field are captured and well-represented by the solver.Figure 6 (left) compares the predicted pressure coefficient distribution along the blade with the reference numerical data by Colonna et al. [6].Curves obtained with Peng-Robinson ideal gas EoS are also added to highlight the effect of different thermodynamic models.The matching with the available reference data is satisfactory and the small differences can be ascribed to the different set of Equations (Euler for reference and RANS for present computations).Figure 6 (right) shows also the value of the non-dimensional wall shear stress along the blade for the three models.Both figures demonstrate an almost perfect matching between PR and SW curves, whereas some differences are evident when using the ideal gas law.More differences between SW/PR models and ideal gas EoS can be seen from the pitch-wise distribution at the outflow section for the flow angle, the Mach number, and the total pressure loss coefficient ξ = (p 01 − p 02 )/(p 02 − p 2 ), as shown in Figure 7.The computations reveal a very similar behavior between Peng-Robinson and Span-Wagner models, whereas the ideal gas EoS shows a very different distribution for the Mach number and the total pressure loss coefficient distributions.In particular, the ideal gas EoS predicts higher peaks in the Mach number distribution, whereas the loss coefficient is lower.These trends are confirmed also by the mixed-out quantities, as reported in Table 5.The mixed-out values for these quantities are calculated for the generic property x with the mass flow average x MO = ( A ρVxdA)/( A ρVdA), where A is the outflow section.

Conclusions and Future Works
A dG-FEM solver with complex thermodynamic models is developed and assessed with reference literature problems characterized by classical and non-classical gas dynamics phenomena.Cubic and multiparameter equations of state are implemented to achieve the best possible accuracy in the determination of thermophysical and transport properties in a RANS framework.A good agreement of results with respect to the reference is obtained and non-classical real gas dynamic phenomena are well captured by the solver.
Future works will cover the implementation of ad hoc numerical procedures, with the aim of achieving a systematical use of the solver for the design process of ORC turbomachinery, e.g., non-reflective boundary conditions and mixing-planes.The generation and use of efficient look up tables to speed up computations with heavy gas models will also be evaluated.

Figure 2 .
Figure 2. Wedge.Mach number contours predicted with the following set of mesh elements and polynomial order: 42,400 and P 0 (left), 2650 and P 4 (center), and 42,400 and P 2 (right).

Figure 3 .
Figure 3. ORC nozzle.Mesh of the blade channel, 5305 hybrid elements (hexahedral in the boundary layer and prisms outside).The geometry is distorted because the blade design is confidential property of the manufacturer.

Figure 4 .
Figure 4. ORC nozzle.Convergence history of the simulation with PR EoS.

Figure 5 .Figure 6 .
Figure 5. ORC nozzle.Mach number (top) and compressibility factor (bottom) contours, P 2 solution.Distorted geometries are depicted because the blade design is confidential property of the manufacturer.

Table 2 .
Time required for the evaluation of the Jacobian matrix of some routines through AD and FDs with ideal gas law.

Table 4 .
Wedge.Angles for the rarefaction shock in the wedge-shaped channel case.

Table 5 .
ORC nozzle.Mixed-out values of the flow angle, Mach number, and total pressure loss coefficient with various thermodynamic models, P 2 solution approximation.