Fully turbulent discrete adjoint solver for non-ideal compressible flow applications

Non-Ideal Compressible Fluid-Dynamics (NICFD) has recently been established as a sector of fluid mechanics dealing with the flows of dense vapors, supercritical fluids, and two-phase fluids, whose properties significantly depart from those of the ideal gas. The flow through an Organic Rankine Cycle (ORC) turbine is an exemplary application, as stators often operate in the supersonic and transonic regime, and are affected by NICFD effects. Other applications are turbomachinery using supercritical CO2 as working fluid or other fluids typical of the oil and gas industry, and components of air conditioning and refrigeration systems. Due to the comparably lower level of experience in the design of this fluid machinery, and the lack of experimental information on NICFD flows, the design of the main components of these processes (i.e., turbomachinery and nozzles) may benefit from adjoint-based automated fluid-dynamic shape optimization. Hence, this work is related to the development and testing of a fully-turbulent adjoint method capable of treating NICFD flows. The method was implemented within the SU2 open-source software infrastructure. The adjoint solver was obtained by linearizing the discretized flow equations and the fluid thermodynamic models by means of advanced Automatic Differentiation (AD) techniques. The new adjoint solver was tested on exemplary turbomachinery cases. Results demonstrate the method effectiveness in improving simulated fluid-dynamic performance, and underline the importance of accurately modeling non-ideal thermodynamic and viscous effects when optimizing internal flows influenced by NICFD phenomena.


Introduction
NICFD is a new branch of fluid-mechanics (NICFD, 2016) concerned with the flows of dense vapors, supercritical fluids, and two-phase fluids, in cases in which the ideal gas law does not apply.
In these flows, the isentropic variation of the speed of sound with density is different if compared to the flow of an ideal gas (Cramer, 1991); thus, the flow field is bound to be quantitatively (Harinck et al., 2009) or even qualitatively different (Zamfirescu et al., 2008).
NICFD internal flows occur in numerous heterogeneous industrial processes.The supersonic and transonic flow through an ORC turbine nozzle is an example (Colonna et al., 2008).Another case in the energy sector is the transonic flow occurring in the compressor of supercritical CO 2 power plant (Pecnik et al., 2012), or of CO 2 capture and sequestration plants (Baltadjiev et al., 2014).Similarly, flows of fluids in dense-vapor or two-phase conditions are relevant in throttling valves, compressors, and ejectors of refrigeration and heat pump systems (Bassi et al., 2000;Bartosiewicza et al., 2006).Turbomachinery and nozzles partly operating in the NICFD regime are common in the oil and gas industry (Boncinelli et al., 2004;Jassim et al., 2008;Pacheco et al., 2013), and these unconventional flows can also occur in pipelines for fuel distribution (Thorley and Tiley, 1987).Dense vapors made of heavy molecules can be used in supersonic wind tunnels instead of air to achieve higher Reynolds numbers, which can be varied almost independently from the Mach number (Anderson, 1991).Finally, supercritical CO 2 nozzle flows are used in the pharmaceutical industry to extract chemicals (Turk, 2000).
The technical and economic viability of these processes can be greatly enhanced if the performance of fluid flow components is improved.Fluid-dynamic shape optimization (FSO) arguably allows a quantum step progress in this respect (Mohammadi and Pironneau, 2004).FSO has played a crucial role in the development of more conventional technologies (Jameson and Reuther, 1994;Kim et al., 2004;Pierret et al., 2007), but it can be even more important in the case of technologies entailing NICFD flows, where design experience and experimental information are much more limited.Hence, concerted research efforts have been recently devoted to develop FSO techniques for NICFD applications, in particular for nozzles and turbomachinery blades (Harinck et al., 2013;Pini et al., 2014Pini et al., , 2015;;Fernandez and Persico, 2015;Persico, 2016).
The FSO of nozzles and of turbomachinery blades can be performed with either gradient-free (Lian and Liou, 2005) or gradient-based (Luo et al., 2014) methods.Gradient-free algorithms only demand the evaluation of the objective function (e.g., genetic algorithms), and they are often coupled with surrogate models to reduce the computational cost (Pierret et al., 2007;Samad and Kim, 2008;Samad et al., 2008).Nevertheless, the number of function evaluations necessary to converge to an optimum solution are comparatively large, and only few design variables can be concurrently optimized (Quagliarella, 1998).Instead, gradient-based methods can reach an optimal solution in far fewer iterations.However, these techniques require not only the computation of the objective function, but also the expensive estimation of its gradient with respect to the design variables.The use of the adjoint method makes the computational cost of the gradient evaluation of the same order of magnitude of that of the objective function, regardless of the number of design variables (Peter and Dwight, 2009).Thus, if the problem involves a large number of design variables and the estimation of the objective function is computationally expensive, adjoint-based methods are the only viable technique.
The already difficult task of linearizing the flow equations (Walther and Nadarajah, 2015) becomes even more challenging in the context of NICFD, where complex thermo-physical fluid models must be adopted to accurately estimate fluid properties, raising the need of specialized numerical methods (Colonna and Rebay, 2004;Rinaldi et al., 2014).Despite that, recent work on the subject (Pini et al., 2014(Pini et al., , 2015) ) has demonstrated the potential of adjoint-based method for the FSO of NICFD flows occurring in ORC turbine cascades.However, this approach was limited to inviscid flows, restricting the adjoint applicability to some supersonic flow cases where viscous effects are negligible.
In this work the adjoint method was extended to fully-turbulent NICFD flows; thus, the approach can be applied, without restrictions, to the FSO of any NICFD application.The adjoint solver was obtained by linearizing the discretized flow equations by means of AD.Nevertheless, the use of AD, if performed by differentiating individual subroutines like in (Pini et al., 2014), still requires additional error-prone steps whenever new numerical schemes or fluid models are added to the source-code.On the contrary, in this work a holistic linearization approach was adopted, whereby AD is applied in a black-box manner to the entire source code.This is accomplished with the help of modern meta-programming features (Sagebaum et al., 2017) in combination with a reformulation of the state constraint into a fixedpoint problem (Albring et al., 2016).The result is a fast and accurate discrete adjoint solver that includes all the flow solver features, such as arbitrary complex equations of state and turbulence models.
The new Reynolds-averaged Navier-Stokes equations (RANS) adjoint solver was developed by leveraging on the open-source software infrastructure of SU2 (Economon et al., 2015), a platform conceived for solving multi-physics Partial Differential Equation (PDE) and PDE-constrained optimization problems using general unstructured meshes.The SU2 flow solver, previously extended to model NICFD flows (Vitale et al., 2015), was adapted to simulate turbomachinery flows.
The developed adjoint solver was naturally integrated into the automatic constrained FSO framework already available in SU2.The optimizer uses the objective function and constraint sensitivities to accordingly re-shape the target geometry by moving the control points of a free-form deformation box.To avoid re-meshing at each design cycle, a linear-elasticity method is applied to propagate the surface perturbation to the entire mesh.
The capability of the new design tool was tested on two exemplary cases: a supersonic and a transonic ORC turbine cascade.The results demonstrate the importance of accurately modeling non-ideal thermodynamic and viscous effects for adjoint-based FSO applied to NICFD applications.This article is organized as follows.The next section describes the SU2 flow solver and its extension to accurately simulate and analyze turbomachinery NICFD flows.The Fluid Dynamic Design Chain section focuses on the derivation of the discrete adjoint solver, also taking into account the surface and volume mesh deformation.The Results and Discussion section reports the application of the new adjoint solver to the re-designing of two typical ORC blades, and discusses the different results obtained using a turbulent or an inviscid approach.Finally, in the last section, before the conclusions, the more accurate NICFD approach is compared with the more standard ideal-gas based method.

Spatial and time discretization
The RANS equations are discretized in SU2 using a finite volume method with a standard edge-based structure on a dual grid with control volumes that are constructed using a median-dual vertex-based scheme (Economon et al., 2015).The PDE semi-discretized integral form is where R(U) is the residual vector obtained by integrating the source term over the control volume Ω and summing up all the projected numerical convective and viscous fluxes associated with all the edges of Ω.
Ad-hoc methods must be used to compute numerical fluxes that are independent from the thermodynamic model of the fluid.For example, the Roe scheme (Roe, 1981) in the generalized formulation presented in (Montagne and Vinokur, 1990) was implemented in SU2 (Vitale et al., 2015).
The time integration is performed with an implicit Euler scheme resulting in the following linear system: where ΔU n := U n+1 −U n and Δt n is the (pseudo) time-step which may be made different in each cell by using the local time-stepping technique (Palacios et al., 2013).Equation (2) can be solved using different linear solvers implemented in the code framework (Economon et al., 2015).Furthermore, non-linear multi-grid acceleration (Borzi, 2005) is available.

Non-reflecting boundary conditions
Non-reflecting boundary conditions (NRBC) were implemented according to the method proposed in (Giles, 1990).Any incoming boundary characteristic δc i BC can be seen as a contribution of two different components: The harmonic boundary solution δ c ^is calculated using the 2D non-reflecting theory, and this is the component that prevents the formation of non-physical boundary reflections.The average component δ c (also known as zero-th fourier mode) is computed according to the standard 1D characteristic- based approach, and allows the user to specify quantities at the boundary.While the calculation of δ c ŵas implemented following the work of Giles, the computation of δ c was reformulated for any arbitrary fluid model based on an Equation of State (EoS).
When NRBCs are imposed at the inlet boundary, the computed average entropy, stagnation enthalpy, and flow direction can differ from the user-specified values (Giles, 1991).Therefore, to ensure that a correct solution is found, the average incoming characteristics are computed by driving the difference between the computed average quantities and the user-specified quantities, i.e., to zero at each time step.The resulted non-linear system is solved via a Newton-Raphson iteration, where the Jacobian of the residuals with respect to the characteristic variables for an arbitrary thermodynamic model can be written as The expression of the enthalpy derivatives and the entropy derivatives are available in SI (i).
The same problem occurs at the outflow boundary; thus, the average incoming characteristic is determined in such a way the exit pressure at convergence is the same as that specified by the user.Since the pressure is directly related to the variation of the incoming characteristic, the average change on the fourth characteristic can be directly computed as Equation ( 7) is also used in case of supersonic outflow conditions whereby the normal component of the outflow Mach number is subsonic.
Once the average and the harmonic component of each incoming characteristic, and the local outgoing characteristic have been computed, the characteristic change at the boundary nodes can be converted back into the change of the primitive variables.This change is summed to the equivalent average values, around which the flow equations are linearized at the boundary.The primitive variables are averaged using the mixed-out procedure, which is the only physically consistent method for averaging flow quantities (Saxer, 1992).The updated values of the primitive variables are used to compute the boundary conservative variables, which, in turn, are used to calculate the convective and the viscous numerical fluxes for the residual in Equation ( 1).

Performance parameters
The performance parameters are computed with the same averaged primitive values used for the NRBCs.Those implemented are the entropy-generation rate, the total-pressure, and the kinetic loss coefficients, i.e., J can be any of the parameters in Equation ( 8), while D is the surface variation.A change in D causes a variation in the surface coordinates, X surf , which in turn requires a continuous deformation of the volume mesh X.The deformed mesh is then used as an input to the flow solver to compute the performance parameter.
The use of the surface nodes as design variables (i.e., D = X surf ) may lead to discontinuous solutions.Thus, the control points (CP) of a Free-Form Deformation (FFD) box were selected as a design variables (Sederberg and Parry, 1986).The surface variation imposed by a perturbation of the FFD CPs is propagated through the volume mesh using the linear elasticity theory as proposed of (Dwight et al., 2009).

Discrete adjoint solver
The discrete adjoint solver was implemented following the approach in (Albring et al., 2016).Based on the design chain described in Figure 1, the optimization problem is formulated as (11) Following the fixed point iteration approach, the constraint on the solution, Equation (10), is the method to solve the flow equations itself described in Equation (2).M(D), instead, is a linear function that formally contains the surface and the volume mesh deformation as shown in Figure 1.
The Lagrangian associated to this problem is T T (13) where the shifted Lagrangian N is and X , U are arbitrary Lagrangian multipliers.Differentiating L with respect to D, and by choosing X and U in such a way that the terms can be eliminated, leads to the adjoint and the mesh sensitivity equations: T T (16) Similarly to the iterative method used for the solution of the flow equations, Equation ( 10), the adjoint equation, Equation ( 15), is solved iteratively with the fixed-point iteration: , ), where U* is a numerical solution of the flow equations.Once the adjoint solution, U , has been found, the mesh node sensitivity, X , is computed with Equation ( 16) and the total derivative of J with respect to the design variables, which is also the total derivative of the Lagrangian, is finally given by

Gradients evaluation with algorithmic differentiation
Algorithmic Differentiation, also known as Automatic Differentiation, is a method to calculate the derivative of a programmed function by manipulating the source-code (Griewank and Walther, 2008).
As an example, consider a generic function y = f(x).It can be demonstrated that the product where y is an arbitrary seed vector, can be obtained by using the reverse mode of AD.Equation ( 19) resembles the second term of the right hand side of the adjoint equation, Equation (15), where f , y, and y denote the fixed-point iteration associated with the flow equations G, the state vector U, and the adjoint vector U , respectively.
The open-source AD tool CoDiPack (Sagebaum et al., 2017), which had been already successfully applied to SU2 (Albring et al., 2016), was also selected for this work.Compared to other AD approaches, CoDiPack exploits the Expression Templates feature of C++.This approach introduces only a small overhead in terms of statement run-time.Thus, the black-box application of AD to complicated non-linear iterative functions becomes feasible.Furthermore, this method is computationally efficient whenever there are code parts containing statements that are not directly involved in the calculation of derivatives, which avoids to perform the extensive derivative-dependency analysis that is necessary in most applications of AD (Bischof et al., 2008).

Results and discussion
The capabilities of the FSO method described in the previous section were demonstrated by redesigning a supersonic and a transonic cascade that are representative of typical cascades adopted in single-stage and multi-stage ORC turbines (Colonna et al., 2015).The illustration of the two test cases follows the same structure.First, the gradient validation is reported, in which the adjoint sensitivities of the objective function and constraint were compared with their finite-difference (FD) equivalent.Second, the results of the optimization are documented.Lastly, the optimization results obtained with the inviscid flow and adjoint solver are also discussed and compared with the ones obtained with the RANS solvers.
In both test-cases the chosen working medium was the siloxane MDM, modeled with the polytropic Peng-Robinson (PR) EoS.The parameters of the model are listed in Table 1.A detailed description of the implementation of the PR EoS in SU2 is available in (Vitale et al., 2015).Turbulent computations were carried out with the k-ω SST model (Menter, 1993), ensuring wall y+ below unity along the blade surface.The optimizer is the minimize routine of the python library scipy.optimize.This routine implements the Sequential Least SQuares Programming (SLSQP) algorithm introduced in (Kraft, 1988).

Supersonic cascade
The supersonic cascade considered in this work was previously investigated and documented in (Colonna et al., 2008;Pini et al., 2014Pini et al., , 2015)).The simulation of the flow around the baseline geometry shows that significant fluid-dynamic penalties would be present because of a strong shockwave forming on the rear suction side of the blade.Therefore, to improve the cascade performance, the entropy-generation rate, s gen , was minimized under the constraint of preserving the baseline mass-flow rate = m m ( ˙˙) b .
The boundary conditions for the simulation are summarized in Table 2.Although total inlet temperature and pressure quantities were given as input for the inflow conditions, these are internally converted to total enthalpy and entropy, which are eventually imposed at the boundary (cf.Equation [4]).
According to the convergence study on the baseline geometry, a maximum number of 1,200 iterations was set for both the flow and the adjoint solver, and the solutions were obtained with an implicit timemarching Euler scheme with a Courant-Friedrichs-Lewy (CFL) number of 40 on 44,000 grid points mesh.
The gradient of the objective function and constraint were first validated against FD.To ease the process, a simple FFD box of 9 CPs was selected (Figure 2a), and the validation was conducted by using a FFD step-size equal to 1E-05.Figure 2b shows that the objective function and constraint gradient provided by the NICFD adjoint very accurately correlate with the correspondent FD values.This result is also a confirmation that the selected convergence criterion guarantees an accurate estimation of the gradient.
Differently from the validation, in which the process can be considered FFD-degree independent, the FFD box for the optimization is composed by 121 CPs, as can be seen in Figure 3a.This choice was made to ensure a high level of design flexibility, which is of primary concern for cascades operating in the supersonic regime, for which slight geometrical modifications can lead to largely different performance.In addition, to prevent unfeasible designs, the sensitivity around the trailing-edge was nullified at each optimization iteration.This has a twofold effect.First, it ensures that a minimum acceptable value of the blade thickness is maintained, without directly introducing a geometrical constraint on the thickness distribution.Second, it alleviates mesh deformation issues at the trailingedge.To guarantee a smooth convergence, the step size of the SLSPQ optimizer was under-relaxed with a value equal to 1E-05 for both the objective function sensitivity and the constraint sensitivity.
The normalized optimization history in Figure 3b shows that, within 5 iterations, the blade entropygeneration rate is reduced by as much as 45%.A similar reduction is found for the other two performance parameters (44.6% for z p,tot and 47.1 % for z kin ), suggesting that the optimal shape is independent from the type of performance parameter selected.The equality constraint on the massflow rate is satisfied with an negligible error of 0.01% with respect to the prescribed value.
The Mach contours of both the baseline and the optimal solution are depicted in Figures 4a and 4b, and the normalized pressure distributions in Figure 4c.In the baseline configuration the simulated flow overaccelerates in the semi-blade region after the outlet divergent section, and this generates an oblique  shock on the rear suction side.The over-acceleration is promoted by the continuing increase of flow passage area in the semi-blade region.On the contrary, in the optimized nozzle, the flow acceleration is more pronounced in the divergent channel, after which the flow keeps smoothly expanding through the re-designed straight semi-bladed channel.This results in the complete removal of the shock-wave from the suction-side, and, consequently, in a much more uniform flow at the outlet.
To gain knowledge of the influence of simulated viscous effects on the design of supersonic ORC cascades, the optimization was repeated using the inviscid adjoint solver.As expected, the use of the inviscid model results into an underestimation of the actual value of the entropy-generation rate, which in the case of the baseline geometry is equal to 129% compared to around 177% predicted by the turbulent solver.Nevertheless, the inviscid and the turbulent optimizations converge to the same optimal geometry, suggesting that the solution of the design problem is driven by the inviscid flow  phenomena (i.e., shock-waves).In addition, the computed boundary-layer viscous losses are very similar in the two studied cases.Further numerical experiments demonstrated that this is only valid if the trailing-edge thickness is constrained to its initial value; since, if the thickness is not constrained, the optimization with the turbulent adjoint would lead to a sharp trailing edge, as expected.

Transonic cascade
The second exemplary test case is representative of the fluid dynamic design of transonic cascades commonly found in multi-stage ORC turbines.The operating conditions of the cascade are listed in Table 3.The inlet total pressure and temperature correspond to the inlet turbine conditions of a superheated MDM thermodynamic cycle; the outlet back pressure was chosen such to obtain a sonic flow at the outlet.
Under the above conditions, the simulated suction-side flow reaches a maximum Mach number of 1.12 at about 90% of the axial-chord length.Afterwards, a shock-wave is generated allowing for a pressure matching downstream of the trailing-edge.The shock causes a sudden increase of the boundary-layer thickness, which eventually results in a more pronounced wake, and, correspondingly, to significant profile losses.Hence, to improve the performance, the cascade was optimized by minimizing its total-pressure loss coefficient.Nevertheless, since the profile-losses are proportional to the flow deflection (Denton, 1993), an unconstrained optimization would ideally converge to a zerodeflection blade.To overcome this issue, the optimization problem was constrained by imposing a minimum tolerable value of the outlet flow angle (β out > 74.0°).
The optimization was performed with a FFD box of degree four in both directions (Figure 5a).However, only 20 of the 25 FFD CPs were chosen as a design variables, and the remaining five, belonging to the short-side of the box next to the trailing-edge were kept fixed to prevent unfeasible designs.The flow and adjoint simulations were marched in time for 2,000 iterations on a grid of 17,500 elements using an implicit Euler scheme with a CFL equal to 30.
In order to verify the accuracy of the adjoint for predominantly viscous flows, the gradient values was validated against the FD ones using a step-size of 5.0E-06.Figure 5b shows that the norm of the two gradients correlates well, confirming the correct linearization of the RANS flow solver.
As Figure 5c shows, the optimization process converges in 16 iterations.The total-pressure loss coefficient is reduced of about 20% with respect to the initial value, and, as expected, similar reductions are also obtained for the other performance parameters.Furthermore, the optimized flow angle is set to the lower bound of 74°so that, given the constant inflow angle, the flow deflection is reduced to its minimum allowable value.
The baseline and the optimal geometry are compared in terms of normalized pressure distribution in Figure 6a.It can be observed that the suction-side velocity peak occurs much more upstream in the optimized blade than in the baseline blade, leading to a smoother flow deceleration on the rear suction side, thus weakening the strength of the shock-wave.This is more evident if the Mach number contours depicted in Figures 6b and 6c are considered.The optimal blade features a reduced outgoing boundary-layer thickness and wake with respect to the baseline.On the contrary, the simulated flow around both geometries display similar characteristics on the pressure side, confirming that most of the losses in transonic ORC cascades are caused by adverse pressure gradients on the rear suction side.
The optimization of the transonic cascade was repeated by using the inviscid flow solver in order to gain insight whether the use of a simpler model can still be advantageous for the design of transonic ORC blades.Interestingly, though the optimization process leads to a reduction of the (inviscid) totalpressure loss coefficient of about 13%, the optimal geometry is practically coincident with the baseline, and the performance improvement predicted by the turbulent solver is negligible (around 1% w.r.t. the initial turbulent value).Apparently, since the main cause of irreversibility are the viscous losses in  the boundary-layer, the inviscid adjoint is inadequate for such cases and only a fully viscous adjoint method is suitable for the optimization of transonic ORC cascades.

Influence of the thermodynamic model on the optimal solution
The nozzle of ORC turbines commonly operate with the fluid in thermodynamic states far from ideal gas conditions.As demonstrated in (Colonna et al., 2008), the use of over-simplified thermodynamic models may lead to inaccurate predictions of the main flow features.Table 4 shows, for example, that there are some discrepancies between supersonic cascade performance values calculated with the Ideal Gas (IG) EoS and those calculated with the more accurate PR EoS.Correspondingly, neglecting nonideal effects in the gradient calculation may lead to sub-optimal blade shape configurations.Nevertheless, as the majority of the available turbomachinery adjoint solvers only implement the polytropic ideal gas model, it is of practical interest to evaluate when the adoption of such model may lead to satisfactory results.
To this purpose, the blade geometries documented in the previous section were re-optimized using the polytropic ideal gas model, and the new optimal solutions were compared with the previous ones.The optimal solution calculated using the polytropic ideal gas model is labeled "Optimal-IG," while that calculated with the PR model maintains the label "Optimal." In order to correctly evaluate the performance of the optimized geometry obtained with the flow and adjoint solver linked to the polytropic ideal gas model, the flow solution around such geometry was computed also with the PR model.This solution was then compared with the one illustrated in the previous section.Except the thermodynamic model, all the other test-case parameters (optimization problems, boundary conditions, meshes, FFD box, etc.) were kept the same.
For the supersonic nozzle test-case, an optimal solution is found within 9 iterations.The predicted entropy generation reduction is around 37%, and the constant mass-flow rate constraint is tightly satisfied.As shown in Figure 7a, the optimizer is able to calculate a solution for which no shock wave occurs at the rear suction-side of the blade, even if the ideal gas model is used for the computation of fluid properties.The obtained geometry is very similar to the one obtained using the PR model, except for a slight offset, see Figure 7b.Even if the geometries are very similar, the calculated performance is different.The performance improvement calculated for the Optimal-IG geometry is 5% lower if compared to that calculated for the Optimal geometry.
Nonetheless, even the use of the inaccurate ideal gas model allows to calculate a substantial improvement with respect to the baseline, though suboptimal.This result can be explained by observing that the geometry optimization occurs mainly after the throat of the nozzle where the computed compressibility factor is close to unity, see Figure 7c.
On the contrary, if the cascade is transonic and operates with the fluid completely in the non-ideal thermodynamic regime, the use of the more accurate thermodynamic model is mandatory (Figure 8a).The simulation corresponding to the Optimal-IG geometry leads to an estimation of the total-pressure loss coefficient that is higher than that of the baseline (6.49% against 5.27%).As shown in Figure 8b, this is due to the fact that the peak velocity is higher, and the normal shock wave on the suction side is stronger.

Conclusion
The technical and economic viability of energy conversion technologies, such as organic Rankine or supercritical CO 2 cycle power systems, often depend on the performance of their turbomachinery components, operating partly in the so-called non-ideal compressible fluid dynamic regime.Due to the limited design experience and to the scarce experimental information on NICFD flows, the design of  these components greatly benefits from automated fluid dynamic shape optimization techniques.Since the achievement of satisfactory performance of the turbomachinery dictates the use of a large number of design variables, adjoint-based FSO is the only viable technique.
A fully-turbulent adjoint method for the design of nozzles and of turbomachinery operating in the NICFD regime was, therefore, devised.The implementation of the adjoint solver and its application to the design of two exemplary test-cases were presented.The following conclusions can be drawn.
1. Despite the high level of additional complexity due to the need of treating accurate fluid thermodynamic models, it was possible to obtain an exact fully-turbulent adjoint.The linearization of the flow solver was performed with an opensource operator-overloading AD tool (CoDiPack).
Given its small overhead in terms of run-time, the AD was applied to the entire flow solver iteration in a black-box manner, thus avoiding a posteriori and ad hoc intervention, which is errorprone.
2. The capability of the new tool was demonstrated by redesigning two typical 2D ORC turbine cascades, in which fluid properties are calculated with the polytropic PR model.In both cases, the optimization performed with the RANS flow and adjoint solver substantially improved the simulated cascades performance, while satisfying the imposed constraints.The optimization process was also repeated assuming inviscid flows.The results emphasize the importance of incorporating the viscous and turbulent gradient contributions, especially for the design of transonic ORC cascades.
3. The potential of the non-ideal compressible RANS adjoint was also shown by comparison with a more conventional ideal gas adjoint method, using the same test cases.The results demonstrate that the simplified approach provides physically inaccurate gradient information leading to suboptimal cascade configurations, if NICFD effects are moderate.Instead, in case of strong NICFD effects, the calculated performance can be even worse than that of the starting geometry.
Future work will be devoted to the extension of the proposed approach to deal with fluid thermodynamic libraries that are external with respect to the flow solver.In this way, the method can leverage the large number of available software.In addition, the capability of optimizing 3D and multi-stage geometries is being implemented.

Supporting information
Thermodynamic derivatives for the NRBCs The enthalpy derivatives are computed as The explicit expression of the above partial derivative for different EoS can be found in (Vitale et al., 2015).

Figure 1
Figure 1 schematically shows the dependence of the objective function, ) = D U D X D J J ( ) : ( ( ), ( )), from the design variables D. In the implementation described here, )J can be any of the parameters in Equation (8), while D is the surface variation.A change in D causes a variation in the surface coordinates, X surf , which in turn requires a continuous deformation of the volume mesh X.The deformed mesh is then used as an input to the flow solver to compute the performance parameter.

Figure 1 .
Figure 1.Formal representation of the evaluation of the performance parameter.

Figure 2 .
Figure 2. Objective function and constraint gradient validation for the supersonic cascade test-case.a) 2D FFD box of degree two on both directions (9 CPs) used for the validation; b) comparison between the FD and adjoint gradient values.

Figure 3 .
Figure 3. Supersonic cascade optimization.a) 2D FFD box of degree 10 on both directions (121 CPs) used for the optimization; b) optimization convergence history.

Figure 4 .
Figure 4. Flow features of the baseline and the optimized geometry of the supersonic cascade test-case.a) mach contour of the baseline geometry; b) mach contour of the optimal geometry; c) normalized pressure distribution around the baseline and the optimal geometry.

Figure 5 .
Figure 5. Objective function and constraint gradient validation, and optimization history for the transonic cascade test-case.a) 2D FFD box of degree four on both directions (25 CPs) used for the validation and the optimization; b) comparison between the FD and adjoint gradient values; c) optimization convergence history.

Figure 6 .
Figure 6.Flow features of the baseline and the optimized geometry of the transonic cascade test-case.a) normalized pressure distribution around the baseline and the optimal geometry; b) mach contour of the baseline geometry; c) mach contour of the optimal geometry.

Figure 7 .
Figure 7. Optimization of the supersonic nozzle geometry using the IG EoS.a) mach contour of the Optimal-IG geometry; b) comparison between the Optimal and the Optimal-IG geometry; c) compressibility factor for the baseline geometry.

Figure 8 .
Figure 8. Optimization of the transonic cascade geometry using the IG EoS.a) compressibility factor around the baseline geometry; b) normalized pressure distribution around the baseline, the Optimal-IG, and the Optimal geometry.

2
The entropy derivatives can be written in terms of the partial derivatives of pressure with respect to temperature Abbreviations NICFD Non-Ideal Compressible Fluid-Dynamics ORC Organic Rankine Cycle AD Algorithmic Differentiation NRBC Non Reflecting Boundary Condition FFD Free-Form Deformation FD Finite Difference PR Peng-Robinson EoS Equation of State SLSPQ Sequential Least SQuares Programming IG Ideal-Gas

Table 1 .
Peng-Robinson EoS parameters for the MDM organic fluid.

Table 2 .
Inlet and outlet boundary conditions values for the supersonic cascade test-case.

Table 3 .
Inlet and outlet boundary conditions values for the transonic cascade test-case.

Table 4 .
Summarized results of the supersonic-nozzle test-case for different EoSs.