DOI: 10.3303/CET24114106 Paper Received: 15 June 2024; Revised: 27 August 2024; Accepted: 6 November 2024 Please cite this article as: Dittmann F.N., Nemer M., Bouallou C., 2024, Discrete Adjoint Shape Optimisation of Supersonic Ejectors for Heat Recovery Refrigeration Cycles, Chemical Engineering Transactions, 114, 631-636 DOI:10.3303/CET24114106 CHEMICAL ENGINEERING TRANSACTIONS VOL. 114, 2024 A publication of The Italian Association of Chemical Engineering Online at www.cetjournal.it Guest Editors: Petar S. Varbanov, Min Zeng, Yee Van Fan, Xuechao Wang Copyright © 2024, AIDIC Servizi S.r.l. ISBN 979-12-81206-12-0; ISSN 2283-9216 Discrete Adjoint Shape Optimisation of Supersonic Ejectors for Heat Recovery Refrigeration Cycles Florian N Dittmann*, Maroun Nemer, Chakib Bouallou CEEP Mines Paris, PSL University, 5 rue Léon Blum, 91120 Palaiseau, France florian.dittmann@minesparis.psl.eu The performance of supersonic ejectors for heat recovery refrigeration cycles is optimised by means of fluid dynamic shape optimisation with a discrete adjoint method to efficiently evaluate gradients of an objective function with respect to an arbitrary number of design variables. The Reynolds-averaged Navier-Stokes equations and the discrete adjoint equations are solved using the finite volume solver SU2. The model is validated on cases from the literature. It is conclusively demonstrated that the method employed can generate a well performing ejector shape from a failed design despite concerns about the well-posedness of the problem given the apparent discontinuity of the entrainment ratio objective function at the critical point in the parameter space. This is achieved in only 10-20 iterations with a cost of one or two direct solutions each. The entrainment ratios predicted with the optimal shape exceed those announced by a leading manufacturer for identical cycle conditions by around 15 %. More importantly, discrete adjoint shape optimisation is here applied successfully for the first time to transonic ejector flows and is shown to be reliable and robust, paving the way for the application of topology optimisation methods to ejectors. 1. Introduction An ejector cycle may be thought of as a fusion of a Rankine cycle and a refrigeration cycle, where the expanding stream of the former and the compressed stream of the latter are brought into direct contact. The device in which this occurs, the ejector, has no moving parts and is designed to accelerate the flow to high velocities to exchange sufficient kinetic energy. Its simplicity makes it reliable and inexpensive compared to a turbine- compressor couple or an ad/absorption cycle, but it is also inherently less energy efficient due to the dominance of dissipative shear forces as opposed to normal forces exerted by blades. 1.1 Working Principle of Ejector Refrigeration Cycles The refrigerant enters the ejector and expands through a nozzle (see Figure 1) which it exits at a high, usually supersonic, velocity to attain a pressure as low as that of the evaporator from where it draws the cold flow, usually in a converging section. The combined flow then decelerates in a diffuser with a diverging outlet section, increasing in pressure as it reaches the condenser. The cycle considered in this work operates at around 8.8 bar condenser pressure and 3.4 bar evaporator pressure close to the saturation curve of the refrigerant R134a (1,1,1,2-Tetrafluoroethane). Its nozzle inlet flow rate is 1.6 kg/s. Figure 1: Lateral section of a generic ejector – illustration of entrainment Ejector designs for refrigeration cycles are constrained by the evaporator and condenser pressures at the intake and the outlet respectively. Those are determined by the cooling needs and the saturation curve of the 631 refrigerant. The ratio of entrained mass flow to nozzle mass directly determines the flow rate through the evaporator and therefore the coefficient of performance (COP) of the cycle. Thus, the objective of any ejector design is to obtain the highest entrainment ratio for given conditions. Due to choked supersonic flow the flow rate does not increase with decreasing outlet pressure below what is called the critical pressure. 1.2 State of the Art of Ejector Optimisation Two-dimensional (axisymmetric) Reynolds-averaged Navier-Stokes (RANS) equations with an eddy-viscosity turbulence closure (such as k-ω SST) are commonly solved to model ejector (Bartosiewicz et al., 2005). Accurate predictions of the entrainment ratio can be expected in the choked flow regime, but results become unreliable above the critical pressure (Del Valle et al., 2015). For refrigeration applications the resolution is further complicated by the need for non-ideal gas thermodynamic models in the form of cubic equations of state or table interpolations. Given the high computational cost of solving these models, a manual parametric analysis only allows to optimise a limited number of interdependent parameters. Gradient-based optimisation uses approximations of derivatives of a global objective, such as maximising entrainment, with respect to variations in all design variables to improve the design. The conventional way to approximate this gradient, the finite difference method, requires at least as many objective function evaluations as there are design variables for each shape modification (iteration), making the computation prohibitively expensive even for a sparse parameterisation of an ejector shape. A common approach found in the recent literature on ejectors to estimate the value of the objective function is to use dimensionally reduced surrogate models (Carrillo et al., 2018), which are only valid for a very limited range of operating conditions and values of design variables of conventional ejector shapes and need to be developed by experts on a case-by-case basis. Unconventional topologies containing multiple inlets or intakes as well as unsteady flow ejectors today make up an active area of research (Han et. al, 2021) including in refrigeration (Van den Berghe et al., 2023). The highly dissipative nature of energy exchanges in ejector flows and their poor adaptability to changing operating conditions is attempted to be mitigated with innovative shapes. Algorithmic exploration of the design space of such devices, in particular by topology optimisation which requires hundreds, if not thousands, of parameters, has been held back by the computational cost and lack of generality of established modelling and optimisation methods. The approach presented in this work can be applied to any shape with any parameterisation and drastically decreases the computational cost by an evaluation of the gradient of the objective function independent of the number of design variables by the discrete adjoint method. We are here taking an important step by applying this method for the first time to a transonic ejector flow of a refrigerant with a classical shape parameterisation (free form deformation). This addresses concerns about the application of the method in cases where the objective function shows an apparent discontinuity, here at the critical point in the parameter space. It also demonstrates the flexibility of the implementation of the discrete adjoint method, which allows to include arbitrary thermodynamic models necessary to model refrigerants under high pressure, while leveraging state of the art methods of rapidly solving systems of partial differential equations with high resolution schemes using the open-source finite volume solver SU2 (Economon et al., 2016). 2. Computational Model and Optimisation Method The Reynolds-averaged Navier-Stokes (RANS) equations for compressible 2D-axisymmetric flows, the k-ω SST turbulence model and a thermodynamic model for a general gas are coupled to describe the ejector flow. Gradient-based shape optimisation is applied to drive design variables towards the maximisation of the entrainment ratio objective function. The gradient of the objective function with respect to changes in design variables is evaluated using the discrete adjoint method. 2.1 Computational Fluid Dynamics Model The conservation of mass, momentum and energy are expressed in terms of time averaged flow variables density ρ, velocity u and energy e (Eq(1) to Eq(3)) and the system is closed with two conservation equations for turbulence kinetic energy kt and the rate of dissipation of turbulence kinetic energy ω (not given here). 𝜕𝜌 𝜕𝑡 + 𝛻 ⋅ (𝜌𝒖) = 0 𝜕𝜌𝒖 𝜕𝑡 + 𝛻 ⋅ (𝜌𝒖𝒖) = 𝛻𝑝 + 𝛻 ⋅ (𝝉 + 𝝉𝒕) 𝜕𝜌𝑒 𝜕𝑡 + 𝛻 ⋅ (𝜌𝒖𝑒) = 𝛻(𝑝𝒖) + 𝛻 ⋅ (𝝉 ⋅ 𝒖 + 𝝉𝒕 ⋅ 𝒖) − 𝛻 ⋅ (𝒒 + 𝒒𝒕) (1) (2) (3) where the viscous stress tensor (assuming no bulk viscosity) and the Reynolds stress tensor are given by 632 𝝉 = 𝜇 (𝛻𝒖 + 𝛻𝒖𝑇 − 2 3 (𝛻 ⋅ 𝒖)𝑰) 𝝉𝒕 = 𝜇𝒕 (𝛻𝒖 + 𝛻𝒖𝑇 − 2 3 (𝛻 ⋅ 𝒖)𝐼) − 2 3 𝜌𝑘𝑡𝑰 (4) (5) the sum of specific internal and kinetic energy (mean flow and turbulence) is given by 𝑒 = 𝑒𝑖 + 1 2 |𝒖|2 + 𝑘𝑡𝑰 (6) and the heat flux and turbulent heat flux are given by 𝒒 = 𝜆𝛻𝑇 𝒒𝒕 = 𝜆𝒕𝛻𝑇 (7) (8) the eddy viscosity μt , turbulence conductivity λt and kinetic energy kt are related by the k-ω SST turbulence model which was implemented in SU2 in its original form (Menter, 1994) at the time of this work. To allow 2D modelling of axisymmetric turbulent flow of a general gas we extended the calculations of the numerical residual in SU2 with the appropriate source terms (convective and viscous contributions due to fluxes in 3D regarded as sources in 2D cylindrical coordinates) in all the conservation equations including the turbulence model. The system of equations is closed by an equation of state (EOS) to determine the thermodynamic state of the refrigerant R134a. 𝑝 = 𝐸𝑂𝑆(𝜌, 𝑒) (9) A lookup table-based EOS was implemented and compared to the cubic Peng-Robinson EOS already implemented in SU2 (Re and Guardone, 2019). Good agreement was obtained especially in the prediction of the entrainment ratio. The latter is therefore used to compute the results presented hereafter due to its slightly lower computational cost. It depends on the following empirical parameters (Table 1). Table 1: Peng-Robinson equation of state parameters for R134a critical temperature critical pressure acentric factor 374.23 K 4.0603 MPa 0.32684 Total pressure conditions were imposed at all inlet and outlet boundaries. This stabilises calculations as a different pressure value at every cell reduced by the rate of change of momentum allows for a realistic velocity profile. Temperatures slightly above the saturation point were imposed at the inlet and intake. Walls were modelled as adiabatic and with a no-slip condition for the velocity and turbulence kinetic energy. 2.2 Discretisation and Resolution For the space discretisation a structured boundary layer, with values of y+<1 almost everywhere, is combined with an unstructured quad mesh generated with a Frontal-Delaunay algorithm. The grid of around 80,000 cells is highly refined at the nozzle exit. It was produced using the open-source tool Gmsh (Geuzaine and Remacle, 2009) with its own scripting language allowing to create a mesh generator producing high quality grids for a wide range of possible ejector shapes including sharp corners and flexible local refinement. Figure 2: Computational Grid of the ejector nozzle and intake The generalised approximate Riemann solver of Vinokur and Montagné (1990) based on a Roe averaged flux, implemented in SU2 by Vitale et al. (2015), is employed for the upwind discretisation of convective terms with MUSCLE reconstruction using the Venkatakrishnan-Wang limiter to obtain second order accuracy. Convective terms in the turbulence equations are discretised with a first order upwind scheme. Central schemes discretise all diffusive terms with second order accuracy. First order implicit time integration is used with locally varying time steps determined by adaptive CFL numbers between 10 and 100 to bring residuals down by 8 orders of magnitude within 10,000 – 20,000 iterations. SU2 solves the system of flow equations and of the turbulence model each in a coupled manner. An incomplete LU decomposition (ILU) preconditioner as well as a three-level multigrid method and a hybrid parallelisation (Blühdorn et al., 2024) are used to accelerate convergence. 633 2.3 Objective Function and Shape Parameterisation A mass flow (entrainment) ratio objective function was implemented within SU2 with mass flow rates at inflow and outflow boundaries of the domain are determined as the sum of flow rates at each grid vertex ν which are calculated using the normal vector n whose magnitude is the area of the element (edge length in 2D). The ratio of the flow rate at the intake to that at the inlet then gives the objective function to be maximised. Free Form Deformation uses boxes shown in Figure 3 to parameterise the entire shape. So-called sensitivities of the objective function are determined with respect to changes in the location of the vertices of these boxes which are mapped to boundary deformation by a linear elasticity model as a function of stiffness which depends on the distance from the wall and the volume of the cell. Vertices on the axis are fixed. Figure 3: Free Form Deformation boxes 2.4 Discrete Adjoint Gradient Computation and Optimisation Algorithm Given a design variable vector α the optimisation problem can be formulated as follows. 𝑚𝑖𝑛 𝜶 𝐽(𝑼(𝜶), 𝑋(𝜶)) 𝑠𝑢𝑏𝑗𝑒𝑐𝑡 𝑡𝑜 𝑹(𝑼(𝜶), 𝑋(𝜶)) = 0 𝑋(𝜶) ≡ 𝑀(𝜶) (12) (13) (14) where U represents the mean flow and turbulence variables, X the locations of the grid points and M the shape deformation function. A solution U of the flow equations is found when the residual R approaches zero for a given X. The discrete adjoint method provides a way to determine the gradient of the objective function J with respect to an unlimited number of design variables α by a single solution of the adjoint equations with the same computational cost as solving the flow equations. This is achieved by explicitly including the constraints in the objective function using Lagrange multipliers which represent the relative magnitude of the gradient with respect to the constraints. Albring et. al (2016) describes the specific implementation of the method to derive, from the implicit numerical scheme (space-time discretisation of the flow equations), the adjoint equations and a duality preserving iteration to solve them for the adjoint variables by exploiting the fixed-point structure of the flow solver. Each iteration as well as the subsequent direct calculation of the sensitivities essentially amount to the evaluation of gradients with respect to changes in the flow variables. To overcome the problem of excessive the working memory requirements of the discrete adjoint method this is done by algorithmic differentiation at the statement level using expression templates and operator overloading, features of the C++ programming language, for a derived type which stores the value and the gradient of a variable (Sagebaum et al., 2019). Finally, by differentiating the transformation of �̃� by M, the gradient of the objective function with respect to the design variables is obtained. 𝑑𝐽𝑇 𝑑𝜶 = 𝑑(𝑀𝑇�̃�) 𝑑𝜶 (15) The python library Scipy provides the required optimisation algorithm which executes objective function and gradient evaluations to drive shape deformations towards an optimum by sequential quadratic least squares programming. The gradient is not required at every iteration of the algorithm and only computed when it is. 3. Results and Discussion The model was first validated on cases by Del Valle et al. (2015). Entrainment ratios were overestimated by 9- 12 % compared to the experimental measurements. However, our predictions were almost identical to those obtained by the authors when the k-ω SST turbulence model was used. This suggests that some sources of dissipation are neglected in the common RANS model in general, possibly viscous dissipation near rough walls and turbulent heat transfer modelled by an inadequate Reynolds analogy. It would be interesting to investigate a systematic inadequacy of eddy-viscosity turbulence modelling in ejectors since they exhibit an unusual combination of free jet flow and wall-bounded adverse pressure gradient flow which are distinctly treated in turbulence model calibration and could explain the typical underestimation of dissipation. Further refining the mesh resulted in identical entrainment ratios in our case while a significantly lower mesh resolution decreases the predicted performance due to artificial dissipation as do inadequate numerical schemes when applied to transonic flow with shocks. An initial condition was created to challenge the optimisation method, a deliberately failed design causing backflow at the intake, given a boiler pressure of 25.6 bar, an evaporator pressure of 3.4 634 bar (86° and 8°C respectively) and an outlet pressure of 8.8 °C. Its diffuser throat is too narrow and its intake section too large, causing a recirculation which dissipates so much energy that the flow cannot overcome the adverse pressure gradient (Figure 4 bottom). This initial condition could typically result from a one-dimensional design model since the too large intake section does not violate the equations and empirical rules of common models. It is not too far removed from the final shape to avoid problems of low mesh quality after free form deformation while it challenges the optimisation algorithm increase and decrease different sections in a range of high performance sensitivity. Only the diffuser optimisation is shown in this publication. Figure 4: Mach number of the optimised design (top) and the failed design (bottom) Figure 5 (left) shows the convergence of the optimisation which increases the entrainment ratio from a negative value to above 0.3 in only 16 iterations. Figure 4 (top) shows the optimised design with supersonic flow across the diffuser. The weakness of the shock suggests that the outlet pressure is just below the critical point where the ejector operates most efficiently. Lowering it decreases the condenser temperature without increasing the evaporator flow rate, due to choked flow, and therefore results in a lower cycle COP. If the outlet pressure is increased the entrainment ratio drops drastically as shown in figure 5 (right). Figure 5: Convergence from a failed design and entrainment ratios as a function of outlet pressure Increasing the evaporator pressure results in a higher entrainment ratio and critical pressure while increasing the boiler pressure results in a much higher critical pressure but a lower entrainment. This is because the flow reaches the exit of the nozzle at a higher pressure which causes a stronger fanning out of the jet and leaves less room for the choked entrained flow (Lamberts et al. 2018). Since less kinetic energy is lost in the entrainment more is left to be converted to pressure at the outlet. Applying the optimisation only to the diffuser under these conditions would result in a larger intake section to compensate this but free form deformation is capable of also increasing the nozzle outlet section which would result in an overall more efficient ejector. Table 2: Diffuser shape coordinates w.r.t axis origin at nozzle inlet (80 is at the nozzle exit) x axis (mm) 80 100 120 240 340 480 520 685 700 y initial (mm) 22.75 18.59 16.33 14.24 13.20 13.31 16.52 41.94 42.00 y opt (mm) 21.76 17.58 15.35 13.87 13.73 14.33 17.63 42.71 42.75 The entrainment ratio just above 0.3 is around 15 % higher than the highest quote of 0.26 from the manufacturer Schutte and Koerting which already significantly exceeded quotes from the manufacturers GEA and Koerting. for identical conditions. However, as our model validations have shown it is often around 10 % overestimated. These gains are modest, but the proposed method does not require an expert and can be applied to arbitrary shapes and operating conditions where a conventional design model is insufficient. The steep drops in the entrainment above the critical outlet pressure suggest that the optimisation problem is ill-posed from a theoretical point of view. While there exists a solution its behaviour does not change continuously with initial conditions. The optimum does not appear to constitute a plateau in the parameter space but rather an edge. As this point is approached by design changes at constant outlet pressure a small change may lead to a large change away from the objective due to the abrupt transition from supersonic flow across the entire section to subsonic flow around the jet where pressure waves are suddenly able to travel upstream. Despite this the applied method has proven remarkably robust which is a promising result for its potential application in topology optimisation of domains with transonic flow. 635 4. Conclusion Fluid dynamic shape optimisation with a discrete adjoint method was successfully applied for the first time to supersonic converting a failed design to a highly efficient one. The approach completely automates the design and can be employed by non-experts using only freely available tools and desktop computing resources due to the low cost of gradient evaluation with respect to an arbitrary amount of design variables. It could be used to design ejectors with unsteady flow and in topology optimisation or the level-set method. Increasing the efficiency of ejectors has the potential to save large amounts of electricity by recovering heat to power refrigeration cycles with a reliable inexpensive device. Nomenclature t – time, s ρ – density, kg/m³ u – velocity vector, m/s e – total specific energy, m²/s² ei – specific internal energy, m²/s² p – pressure, kg/(m s²) τ – viscous stress tensor, kg/(m s²) q – heat flux vector, kg/s³ μ – dynamic viscosity, kg/(m s) T – temperature, K λ – conductivity, kg/(s m K) t – turbulence (subscript) I – identity matrix U – flow solution vector J – objective function X – mesh coordinates vector α – design variables vector R – numerical residual vector M – mesh deformation function n – normal vector References Albring T.A., Sagebaum M., Gauger N.R., 2016, Efficient aerodynamic design using the discrete adjoint method in SU2. 17th AIAA/ISSMO multidisciplinary analysis and optimization conference, DOI: 10.2514/6.2016- 3518. Bartosiewicz Y., Aidoun Z., Desevaux P., Mercadier Y., 2005, Numerical and experimental investigations on supersonic ejectors. International Journal of Heat and Fluid Flow, 26(1), 56-70. Blühdorn J., Gomes P., Aehle M., Gauger N.R., 2024, Hybrid Parallel Discrete Adjoints in SU2. arXiv preprint ,, accessed 09.05.2024. Carrillo J.A.E., de La Flor F.J.S., Lissén J.M.S., 2018, Single-phase ejector geometry optimisation by means of a multi-objective evolutionary algorithm and a surrogate CFD model. Energy, 164, 46- 64. Del Valle J.G., Sierra-Pallares J., Carrascal P.G., Ruiz F.C., 2015, An experimental and computational study of the flow pattern in a refrigerant ejector. Validation of turbulence models and real-gas effects. Applied Thermal Engineering, 89, 795-811. Economon T.D., Palacios F., Copeland S.R., Lukaczyk T.W., Alonso J.J., 2016, SU2: An open-source suite for multiphysics simulation and design. AIAA Journal, 54(3), 28-846. Geuzaine C., Remacle J.F., 2009, Gmsh: A 3‐D finite element mesh generator with built‐in pre‐and post‐ processing facilities. International Journal for Numerical Methods in Engineering, 79(11), 1309-1331. Han J., Feng J., Hou T., Peng X., 2021, Performance investigation of a multi‐nozzle ejector for proton exchange membrane fuel cell system. International Journal of Energy Research, 45(2), 3031-3048. Lamberts O., Chatelain P., Bartosiewicz Y., 2018, Numerical and experimental evidence of the Fabri-choking in a supersonic ejector. International Journal of Heat and Fluid Flow, 69, 194-209. Menter F., 1994, Two-equation eddy-viscosity turbulence models for engineering applications. AIAA Journal, 32, 269-289. Papapetrou M., Kosmadakis G., Cipollina A., La Commare U., Micale G., 2018, Industrial waste heat: Estimation of the technically available resource in the EU per industrial sector, temperature level and country. Applied Thermal Engineering, 138, 207-216. Re B., Guardone A., 2019, An adaptive ALE scheme for non-ideal compressible fluid dynamics over dynamic unstructured meshes. Shock Waves, 29(1), 73-99. Sagebaum M., Albring T., Gauger N.R., 2019, High-performance derivative computations using codipack. ACM Transactions on Mathematical Software (TOMS), 45(4), 1-26. Van den Berghe J., Dias B.R., Bartosiewicz Y., Mendez M.A., 2023, A 1D model for the unsteady gas dynamics of ejectors. Energy, 267, 126551. Vinokur M., Montagné J.L., 1990, Generalized flux-vector splitting and Roe average for an equilibrium real gas. Journal of Computational Physics, 89(2), 276-300. Vitale S., Gori G., Pini M., Guardone A., Economon T.D., Palacios F., Alonso J.J., Colonna P., 2015, Extension of the SU2 open source CFD code to the simulation of turbulent flows of fluids modelled with complex thermophysical laws. In: 22nd AIAA computational fluid dynamics conference, DOI: 10.2514/6.2015-2760. 636