- Email: [email protected]

Contents lists available at ScienceDirect

Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

Robust optimization of dense gas ﬂows under uncertain operating conditions P. Cinnella *, S.J. Hercus Laboratoire DynFluid, Arts et Métiers ParisTech, 151 Bd de l’Hôpital, 75013 Paris, France

a r t i c l e

i n f o

Article history: Received 2 March 2010 Received in revised form 9 June 2010 Accepted 24 June 2010 Available online 30 June 2010 Keywords: Uncertainty quantiﬁcation Polynomial chaos Robust optimization Dense gas dynamics

a b s t r a c t A robust optimization procedure based on a multi-objective genetic algorithm (MOGA) is used to generate airfoil proﬁles for transonic inviscid ﬂows of dense gases, subject to uncertainties in the upstream thermodynamic conditions. The effect of the random variations on system response is evaluated using a non-intrusive polynomial chaos (PC) based method known as the probabilistic collocation method (PCM). After initial PCM simulations which showed that the dense gas system was highly sensitive to input parameter variation, a multi-objective genetic algorithm coupled to the PCM produced a Pareto front of optimized individual geometries which exhibited improvements in mean performance and/or stability over the baseline NACA0012 airfoil. This type of analysis is essential in improving the feasibility of organic Rankine cycle (ORC) turbines, which are typically designed to recover energy from variable sources such as waste heat from industrial processes. Ó 2010 Elsevier Ltd. All rights reserved.

1. Introduction Dense gases are deﬁned as single phase vapors, characterized by complex molecules and moderate to large molecular weights. Flows of dense gases are of considerable interest for many applications in energy production [1–3], refrigeration [4,5], and chemical processing [6]. For dense gases operating at temperatures and pressures of the same order of magnitude of their critical point, the ideal gas approximation is no longer valid and real gas effects become inﬂuential in the dynamic behavior of the ﬂuid. Real gas effects are particularly strong for highly complex molecules which display a non-monotonic variation of the sound speed with respect to pressure. Some ﬂuids characterized by particularly high molecular complexity have been theoretically predicted [7–10] to display a region of negative values of the fundamental derivative of ﬂuid dynamics C [11] in the vapor phase. This thermodynamic region is called the inversion zone [12] and the locus of thermodynamic states in the vapor phase such that C < 0 is referred to as the transition line. Fluids of this type are said to possess Bethe– Zel’dovich–Thompson (BZT) properties. For transonic ﬂows of BZT type ﬂuids, compression shock waves may be suppressed if the upstream thermodynamic conditions are selected within the inversion zone, since for C < 0 compression shocks violate the second principle of thermodynamics [12]. Harnessing the unusual properties of BZT type ﬂuids could potentially have a signiﬁcant impact on the performance of low temperature energy conversion

* Corresponding author. Tel.: +33 1 44246278; fax: +33 1 44246275. E-mail address: [email protected] (P. Cinnella). 0045-7930/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.compﬂuid.2010.06.020

cycles, such as in the case of organic Rankine cycle (ORC) turbines [2,3,13]. Whereas a traditional Rankine cycle operates with water vapor as the working ﬂuid, ORC turbines use an organic ﬂuid such as hydrocarbons, silicon oils or other organic refrigerants. The low critical temperature, high density and elevated heat capacities of these ﬂuids result in high suitability for low temperature operation, even as low as 80 °C [14]. Furthermore, the slope of the saturated vapor line for organic ﬂuids eliminates the risk of condensate at the turbine outlet, without heating the working ﬂuid into the superheated vapor region. ORCs represent a promising technology for the development of widely distributed, small yield (less than 1 MW) thermal energy conversion devices [15]. The heat source in this type of conﬁguration would typically be an external, low quality heat source such as solar thermal collectors or waste heat from industrial processes. The low temperature differential in this case would promote a high second-law efﬁciency. However, this small temperature jump (and therefore enthalpy jump) diminishes power output and ﬁrst law efﬁciency, which is one of the primary drawbacks to the development of real-world ORCs. A trade-off is therefore necessary between these requirements, taking into account the availability, quality (temperature) and cost of the heat source. Although residual heat from industrial processes or solar energy may be of low quality, the cost per kW of incoming energy is often very small or free. Practically speaking, low temperature operation favors widely distributed energy conversion due to smaller unit size and cost, the use of standard materials and improved safety considerations. Additional practical advantages are the self lubricating nature of organic ﬂuids and the absence of a superheating requirement. ORC turbines typically utilize a

1894

P. Cinnella, S.J. Hercus / Computers & Fluids 39 (2010) 1893–1908

single-stage turbine operating in the transonic/supersonic regime, where a major source of losses arises from the formation of shock waves and their interaction with boundary layers. The use of organic ﬂuids possessing BZT properties could potentially avoid shock formation and the consequent losses, if only turbine expansion occurred entirely within or very close to the inversion zone. To achieve this, a reduction in the temperature difference between the heater and condenser stages is generally required. As a consequence, the isentropic efﬁciency of the turbine is increased at the cost of diminished global power output and overall cycle efﬁciency. This important disadvantage has been an impediment to the development of ORC type turbines capable of fully exploiting BZT effects. In order to improve the feasibility of BZT ﬂows in ORC turbines, a compromise is suggested in Cinnella and Congedo [16]. The proposition is to allow the dense gas ﬂow to evolve partially outside the inversion zone, thus permitting loss inducing compression shocks and mixed waves. Such waves are expected to be relatively weak if the thermodynamic states on either side of the discontinuity are in the vicinity of the transition line [12]. These weak waves would result in associated losses which are lower than in the case of an ideal gas, leading to increased turbine efﬁciency. Operating the turbine cascade partially outside the inversion zone also allows for a larger difference between the maximum and minimum temperature of the thermodynamic cycle, boosting the overall power output. However, this approach of allowing variation in the upstream thermodynamic state of the dense gas ﬂow introduces a complex challenge at the modeling stage of development. ORC turbines are often designed to exploit waste heat from other processes, where inlet operating conditions are likely to be highly variable. For the BZT working ﬂuids currently considered, the inversion zone, where the turbine isentropic efﬁciency is expected to be the most favorable, is found to be relatively small and in close proximity to the saturated vapor curve. As a result, variations in upstream thermodynamic conditions may induce dramatic changes in the ﬂow physics and ultimately in the system performance. Given the relatively small size of the inversion zone and the variability of proposed energy sources, the physical behavior of the system is likely to be highly sensitive and difﬁcult to model using a classic approach. The variability of the operating conditions represents a crucial issue for all ORC turbines, not restricted to those using BZT working ﬂuids. Another important point is that the serialization of ORC turbogenerators leads to the development of a standard turbine geometry which is applied to a relatively wide range of situations. As a consequence, even when turbine admission conditions remain relatively stable for a given application, the ﬁnal performance may vary considerably according to the proximity of these conditions to the design point. Taking account of the variability of operating conditions is therefore critical to ensure optimal turbine performance and durability. Classic computational ﬂuid dynamics (CFD) simulations are generally deterministic: given a set of constant input parameters, the calculation returns a single result. Exponential growth in computing power has recently permitted the application of stochastic approaches to CFD. A stochastic study aims to incorporate the effects of input parameter variability on the ﬁnal solution, in order to determine a feasible optimal set point given a wide range of system operation scenarios. A deterministic only approach can lead to over-optimization of a mechanical system, which may occur when the theoretical optimal set point is difﬁcult to achieve given ﬂuctuating real world conditions. Uncertainty quantiﬁcation consists of measuring the system response to random and unknown variations in input parameters using stochastic analysis. The AIAA [17] provides a deﬁnition of uncertainty as a ‘‘potential deﬁciency in any phase or activity of the modeling process that is due to lack of knowledge”. Sources of uncertainty may be intrinsically related to boundary conditions, physical models, discretization and solu-

tion errors and geometrical uncertainty due to machining tolerances. This type of analysis leads to an approach known as riskbased design, where the distribution of possible solutions is used to integrate the concept of risk as a fundamental design constraint into the development process. Improved knowledge of the effect of input uncertainties will improve the reliability of a risk-based design approach, leading to increased design conﬁdence, reduced risks, improved safety and reﬁnement of system operation range through reduction of over-optimization. Dense gas ﬂows are particularly suited to stochastic analysis due to their high sensitivity to variations in upstream thermodynamic conditions. A recent study [18] has already applied a stochastic analysis to a BZT type dense gas ﬂow simulation over an airfoil in order to quantify the effect of uncertainties of parameters in different models of the equation of state (EOS). The most basic type of stochastic analysis is the Monte Carlo (MC) method, originally conceived at the end of WWII during the development of the atomic bomb. Essentially a brute force approach, the MC method consists of randomly sampling N values of an input variable with a known or supposed input distribution, and then calculating the deterministic solution for each input value. The resulting output distribution can then be used to calculate the statistical moments of the system. Theoretically, the MC method will converge to the exact stochastic solution when the number of samples N ? 1. In practice, several thousand samples are required to obtain an acceptable level of accuracy. Although efﬁciency improvement methods such as variance reduction or latin hypercube sampling [19] can reduce the number of deterministic simulations required, the MC method is very computationally intensive when applied to ﬂuid ﬂow simulations. Nevertheless, the method proves to be a useful benchmarking tool for the stochastic analysis of simple CFD problems. An additional advantage of the method is that it is non-intrusive, where the simulation is considered an independent ‘‘black box” and direct modiﬁcation of the CFD code is not required. A second type of stochastic methods, known as Moment Methods, are able to predict the statistical moments with a similar precision to the MC method, but with a signiﬁcantly smaller computational effort. In the case of a 1D viscous Burgers equation with random input viscosity, Walters and Huyse [20] observed that N = 1000 MC simulations produced the same error as a second order Moment Method approximation, but with a computational cost roughly 300 times superior. Moment Methods calculate an approximation of the statistical moments via truncated Taylor series expansions about the expected value of the random input parameters. In this case only the statistical moments are determined, whereas the MC method delivers a complete output distribution. An additional disadvantage of the Moment Methods are that the Taylor expansions require the calculation of ﬁrst and second order sensitivity derivatives, which can be difﬁcult to determine for complex CFD problems, often requiring dedicated automated differentiation tools such as INRIA’ s TAPENADE software [21]. A more sophisticated class of stochastic analysis methods are known as the polynomial chaos (PC) methods. PC methods characterize the behavior of a random process via series expansions of orthogonal chaos polynomials (Lucor et. al [22]). Essentially a spectral approach, PC methods can provide detailed statistical information of system response to input parameter variations at a fraction of the computational cost of the MC method. Knio and LeMaître [23] highlight several advantages of PC methods which justify their implementation in the study of dense gas ﬂows. It is stated that PC methods can efﬁciently obtain accurate estimates of uncertainty for models which describe physical phenomena in terms of partial differential equations. Additionally, they have the capacity to deal with situations exhibiting strong non-linear dependence on random input variations. The classic PC approach is an intrusive meth-

P. Cinnella, S.J. Hercus / Computers & Fluids 39 (2010) 1893–1908

od, which requires signiﬁcant modiﬁcations of the underlying deterministic ﬂow solver. The intrusive approach may result in improved computational efﬁciency for highly dimensioned stochastic problems, however it dramatically decreases ﬂexibility and ease of implementation. For this reason, several PC based non-intrusive methods have been proposed in the literature [19,24,25]. In the present study, the probabilistic collocation method (PCM) proposed by Loeven et al. [19,24] is selected due to its exponential convergence and non-intrusive properties. The aim of the present work is to demonstrate the feasibility and usefulness of robust design techniques for dense gas ﬂows, given uncertainties applied to the upstream thermodynamic conditions. The robust design procedure involves the use of a PCM stochastic solver coupled with a multi-objective genetic algorithm (MOGA). At this early stage of the research work, we focus our attention on a simpliﬁed geometry, namely an airfoil placed into a dense gas stream with randomly varying thermodynamic conditions. The study is divided into two main sections. The ﬁrst section considers uncertainty quantiﬁcation applied to transonic dense gas ﬂows using a polynomial chaos (PC) based stochastic approach known as the probabilistic collocation method (PCM). Initially we validate the PCM developed by Loeven et al. [19,24] with the known test case of diphasic ﬂuid ﬂow in a porous medium, modeled by the Buckley–Leverett (BL) equation. This test model involves a non-convex ﬂux function, which may permit non-classical shock wave behavior such as combined shock/fan waves. These complex waves can also arise in the case of dense gas models, thus the BL system represents a pertinent test case with a known analytical solution for comparison. The validated PCM algorithm is then applied to an existing solver [26] of a inviscid transonic dense gas ﬂow over half of a symmetrical NACA0012 airfoil at M1 = 0.95. The working ﬂuid used is pf-perhydroﬂuorene (C13F22, commercially known as PP10) with the thermodynamic properties modeled by the Martin–Hou equation of state. The aerodynamic performance of the airfoil at 0° incidence is examined with Gaussian random variations in the upstream pressure and temperature. In the second section of this study, a PCM stochastic analysis is coupled with an existing multi-objective Pareto-based genetic algorithm (MOGA), previously applied to transonic dense gas ﬂows [27–29]. The robust optimization procedure generates a series of optimized 2D airfoils for dense gas ﬂows, based on minimization of the mean and standard deviation of the drag coefﬁcient. The high computational cost of the genetic algorithm is mitigated by the use of the Richardson extrapolation method (REX) on each individual calculated. This mesh extrapolation method, already applied to robust proﬁle optimization of dense gas ﬂows [29], accelerates convergence and increases the solution accuracy for each individual, which can improve MOGA convergence. Optimal individual geometries from the MOGA optimization process are selected and their aerodynamic performance analyzed and veriﬁed with both deterministic and stochastic simulations. Finally, a posteriori testing of the optimized individuals is carried out in the turbulent ﬂow regime using the Baldwin–Lomax turbulence model.

2. Governing equations and ﬂow solver 2.1. Physics of dense gas ﬂows In this study, we focus our attention on a class of dense gases of the retrograde type (gases which superheat when expanded) known as Bethe–Zel’dovich–Thompson (BZT) ﬂuids. In a speciﬁc range of thermodynamic conditions in the vapor region, BZT gases exhibit non-classical dynamic behavior, such as expansion shock waves, mixed shock/fan waves and splitting shocks. This unusual

1895

behavior arises when the value of C, the fundamental derivative of gas dynamics [11], becomes negative:

C¼1þ

q @a a @q

ð1Þ s

where a is the speed of sound, q the ﬂuid density and s the entropy of the ﬂuid. The fundamental derivative measures the rate of change of the sound speed in isentropic perturbations. For perfect gases, C = (c + 1)/2 > 1, where c (the speciﬁc heat ratio of the ﬂuid) is strictly greater than 1 due to thermodynamic stability requirements. For the majority of classical working ﬂuids, an isentropic compression leads to an increase in the speed of sound. For a BZT ﬂuid, if C < 1, then ð@[email protected] qÞs < 0, and the ﬂow exhibits a reversed variation in the speed of sound, where a grows in isentropic expansions and falls in isentropic compressions. An important consequence of this phenomenon is observed when considering a weak shock wave. It can be shown [12] that the entropy difference Ds across the shock is given by the expression:

Ds ¼

a2 C ðDv Þ3 þ OððDv Þ4 Þ v 3 6T

ð2Þ

where D is the change in a ﬂuid property across the shock, v = 1/q is the speciﬁc volume, and T is the absolute temperature. It follows that in regions where C < 0, compression shocks cannot form as a consequence of the entropy inequality, whereas expansion shocks are theoretically admissible. The region where C < 0 is referred to as the inversion zone, and the C = 0 contour is known as the transition line. These features can be observed on the state diagrams shown in Fig. 1 for the heavy perﬂuorocarbon pf-perhydroﬂuorene, commercially known as PP10. The thermodynamic properties of this gas are modeled using the realistic Martin–Hou equation of state. The T s diagram displayed in Fig. 1 shows the retrograde behavior of the ﬂuid: the liquid/vapor coexistence curve exhibits a positive slope up to near-critical conditions. Numerical simulation of BZT ﬂows over isolated airfoils and wings [30,31] has shown that with a negative upstream value of C, the ﬂow remains subsonic throughout the entire domain, whereas a shock is observed for a perfect gas at the same operating conditions. BZT effects can essentially delay shock wave formation and the related losses, increasing the critical Mach number to near-sonic speeds. This is in clear contrast to the classical behavior of regular gases such as air, oxygen, nitrogen and steam that exhibit much lower critical Mach numbers (around 0.8 or even less). As a result of further numerical simulations conducted by Brown and Argrow [2], Cinnella et al. [32], Congedo et al. [33], the conclusion has been made that the use of BZT ﬂuids in turbine design could lead to a potential increase in efﬁciency of approximately 3–8%, with advancement in ﬂuid dynamic design and non-classical effect exploitation. 2.2. Dense gas model and ﬂow solver Dense gas ﬂows are governed by the equations for equilibrium for non-reacting ﬂows. In the present study we consider the Euler equations, written in integral form for a control volume X with boundary @ X:

d dt

Z X

w dX þ

Z

f n dS ¼ 0

ð3Þ

@X

In Eq. (3), w is the conservative variable vector, where

w ¼ ðq; qv ; qEÞT n is the outer normal to @ X, and f, is the ﬂux density:

f ¼ ðqv ; qI; qvv; qv HÞT where v is the velocity vector, E the speciﬁc total energy, H = E + p/q the speciﬁc total enthalpy, p is the pressure and I is the unit tensor.

1896

P. Cinnella, S.J. Hercus / Computers & Fluids 39 (2010) 1893–1908

Fig. 1. (left) Pressure-volume diagram and (right) temperature-entropy diagram for PP10 (C13F22) using the Martin–Hou equation of state. The subscript c represents a critical value.

The preceding equations are completed by a thermal equation of state:

p ¼ pðqðwÞ; TðwÞÞ

ð4Þ

and by a caloric equation of state for the speciﬁc internal energy e, which must satisfy the compatibility relation:

e ¼ eðqðwÞ; TðwÞÞ Z T Z 0 ¼ er þ cv ;1 ðT 0 Þ dT

q

qr

Tr

" # @p dq0 T p 02 @T q q

ð5Þ

In Eq. (5), cv,1 is the ideal gas speciﬁc heat at constant volume, quantities with a prime superscript are auxiliary integration variables, and subscript r indicates a reference state. The caloric equation of state is completely determined once a variation law for cv,1 has been speciﬁed. In the present work, the gas response is modeled through the comprehensive thermal equation of state of Martin and Hou [34], which provides a realistic description of the gas behavior close to saturation conditions [35,36] and requires just a minimum amount of thermodynamic input data:

p¼

5 X RT fi ðTÞ þ v b i¼2 ðv bÞi

ð6Þ

where R is the gas constant, b is the molecular covolume, and the functions fi(T) are of the form: T

fi ðTÞ ¼ Ai þ Bi T þ C i ekT c

n T Tc

The governing equations are discretized using a cell-centered ﬁnite volume scheme for structured multi-block meshes of third-order accuracy, which allows the computation of ﬂows governed by an arbitrary equation of state [26]. The scheme is constructed by correcting the dispersive error term of the second-order-accurate Jamesons scheme [37]. The use of a scalar dissipation term simpliﬁes the scheme implementation in the case of highly complex equations of state and greatly reduces computational costs. In order to preserve the high accuracy of the scheme on non-Cartesian grids, the numerical ﬂuxes are evaluated using weighted discretization formulas, which take into account the stretching and the skewness of the mesh: this ensures true third-order accuracy on moderately deformed meshes and at least second-order accuracy on highly distorted meshes (see [38] for details). The equations are then integrated in time using a four-stage Runge–Kutta scheme [37]. Local time stepping, implicit residual smoothing and multigrid acceleration are used to efﬁciently drive the solution to the steady state. For external ﬂows, non-reﬂecting boundary conditions based on a multidimensional method of characteristics are applied at the far-ﬁeld boundaries and an adiabatic wall condition is imposed at solid boundaries. The accuracy of the numerical solver has been demonstrated in previous works [16,26], and will not be discussed further.

3. Uncertainty quantiﬁcation for dense gas ﬂows

with Tc the critical temperature. The gas dependent coefﬁcients Ai, Bi and Ci can be expressed in terms of the critical temperature and pressure, the critical compressibility factor, the Boyle temperature and one point on the vapor pressure curve. A power-law of the form:

cv ;1 ðTÞ ¼ cv ;1 ðT c Þ

2.3. Flow solver

ð7Þ

is used to model variations of the low density speciﬁc heat with temperature, where n is a material dependent parameter. Eqs. (6) and (7) completely determine the gas thermodynamic behavior and are referred to as the Martin–Hou model.

3.1. Polynomial chaos methods for uncertainty quantiﬁcation The most general form of the PC method is known as the generalized polynomial chaos (GPC) method [39]. The broad idea behind the GPC method is to consider an uncertainty as a generator of a new dimension upon which the ﬁnal solution depends. Mutually orthogonal functions are used to create a convergent expression which describes this new dimension. The type of orthogonal polynomial used in the chaos decomposition depends on the form of the distribution of the random input variable. Each input distribution is associated with a corresponding orthogonal polynomial according to the Askey scheme (summarized in Xiu and Karniadakis [39]). For example, in

P. Cinnella, S.J. Hercus / Computers & Fluids 39 (2010) 1893–1908

the case of a normal or Gaussian input distribution, the corresponding orthogonal polynomial is the Hermite polynomial. The implementation of the GPC method involves the decomposition of a random function into deterministic and stochastic components. In the case of a velocity ﬁeld u subjected to random input variable ﬂuctuations n(h), the decomposition is given by:

uðx; t; nðhÞÞ ¼

P GPC X

ui ðx; tÞ wi ðnðhÞÞ

ð8Þ

rin correspond respectivelyptoﬃﬃﬃ the mean and standard deviation of

the input distribution. The 2 factor arises from the difference between the classic expression used to calculate the Gaussian quadrature weights (Eq. (11)) and the form of the Gaussian probability distribution function. Finally, the output mean lu and variance r2u can now be calculated:

lu ¼

where the chaos polynomial wi is interpreted as the stochastic part of the decomposition. It is the random basis function corresponding to the ith ﬂuctuation mode, and the deterministic solution ui(x, t), is its amplitude. The total number of terms in the sum is P GPC ¼ ðnþpÞ! , n!p! which depends on the order p of the polynomial chaos and the number of random input variables n. Eq. (8) can be substituted into the differential equations describing the system and then solved using a Galerkin method to ﬁnd the modes ui of the solution. These modes provide statistical information to characterize the variation of the output solution. Recent research by Loeven et al. [19,24], has developed an efﬁcient non-intrusive variation on the standard generalized polynomial chaos (GPC) method. Based on the idea of a standard chaos transformation, the probabilistic collocation method (PCM) approach consists of two important modiﬁcations to the classic method. The ﬁrst modiﬁcation is that a chaos version of Lagrange interpolation is used to approximate the chaos polynomial wi(n) even with a minimum of two collocation points. The second modiﬁcation is to use Gaussian quadrature to compute the Galerkin projection and the integration of the distribution function approximation. In practical terms, the probabilistic collocation method consists of calculating the deterministic solution at selected points (nodes) in the input distribution, then multiplying the solutions by a weighting function in order to compute output statistical information. Here the PCM is presented considering the case of a Gaussian random input distribution, therefore a Hermite polynomial is used and Gauss–Hermite chaos quadrature is employed to compute the Galerkin projection. In the case of a non-Gaussian distribution, the approach is modiﬁed in accordance with the Askey framework in terms of the orthogonal polynomial selected [39]. The decomposition of the solution into deterministic: ui(x, t), and stochastic: hi (n(h)), parts is presented below: P PCM X

ui ðx; tÞ hi ðnðhÞÞ

ð9Þ

i¼1

hi ðnðhÞÞ ¼

Np Y nðhÞ nðhk Þ nðh i Þ nðhk Þ k¼1

ð10Þ

k–i

where ui(x, t) is the deterministic solution at the collocation point n(hi). In the PCM, p is the order of the quadrature polynomial and the number of collocation points is given by PPCM = pn, where n represents the number of random input variables. The term hi is the Lagrange interpolating polynomial chaos of order Np = p 1 that passes through the PPCM collocation points, with hi(n(hk)) = dik. The collocation points hi are chosen such that they correspond to the Gaussian quadrature points, which are simply the roots of the quadrature polynomial. Finally, in order to compute the statistics of the output, the solution is integrated using Gauss–Hermite chaos quadrature. In this case, the quadrature nodes correspond to the collocation points hi, and the weights wi are computed by:

wi ¼

P PCM X

ui ðx; tÞ wi

ð12Þ

ðui ðx; tÞÞ2 wi ðlu Þ2

ð13Þ

i¼1

i¼1

uðx; t; nðhÞÞ ¼

1897

pﬃﬃﬃﬃ 2n1 n! p n2 ½Hn1 ðhi Þ2

ð11Þ

Additionally, in the case of Gauss–Hermite quadrature, the followpﬃﬃﬃ ing transformation is required: nðhÞ ¼ lin þ 2rin hi , where lin and

r2u ¼

P PCM X i¼1

In the case of multiple input variables, the PCM is slightly modiﬁed. The transformed input vector becomes n(h) = {ns=1(h), ns=2(h), . . ., ns=n(h)}, and Eq. (10) can be rewritten in vector form:

2 3 Np n Y 6Y ns ðhÞ ns ðhk Þ 7 hi ðnðhÞÞ ¼ 4 5 n ðh Þ ns ðhk Þ s¼1 k¼1 s i

ð14Þ

k–i

In order to determine output statistics given multiple uncertain input parameters, Eqs. (12) and (13) must be analytically determined with the multi-variable PC expansion (c.f. Eq. (9), Xiu and Karniadakis [39]) and the standard deﬁnitions of the statistical moments. This procedure is relatively time consuming and complex, especially with several uncertain input variables (i.e. high values of n). A more practical approach to determine the solution statistics is to reconstruct the problem using a simple MC method on the Lagrange interpolation equation, denoted the reconstructed Monte Carlo method (RMC). Once the deterministic solution at each collocation point, ui(x, t), has been determined, a MC analysis is used to generate a large number, M, of values for ns(h), thereby constructing M possible variations of hi(n(h)). Since Eq. (9) is linear in the terms ui(x, t), the RMC can be carried out inexpensively for large values of M. The result is a complete set of output solutions, from which the statistical moments can be easily calculated. In terms of calculation cost, both PC methods show signiﬁcant improvements over the MC analysis and Moment Methods. Walters and Huyse [20] demonstrate the advantage of the GPC method by applying it to a 1D viscous Burgers equation with a random variation on the input viscosity. Even at small values of p, the order of the quadrature polynomial, the GPC method obtains accurate estimates of output statistical information. With p = 1, the variance of the output solution is found to be of the same order as a second order Moment Method. With p = 2, the error on the variance is one order of magnitude smaller than in the case of 1000 MC simulations. As it is critical to reduce the number of deterministic simulations to be carried out for computationally intensive CFD simulations, both PC methods are therefore advantageous. In order to consider potential performance differences between the two PC methods, Loeven et al. [24] conducted a comparison between the PCM and GPC method. While both methods demonstrate exponential convergence with respect to the order of the polynomial, it is noted that for increasing values of n, PPCM > PGPC. Although this would suggest that fewer simulations are required for the GPC method to obtain convergence, the non-intrusive nature of the PCM provides a substantial increase in ﬂexibility. The PCM is therefore particularly suited to the study of both complex CFD simulations, without requiring modiﬁcation to the CFD code framework. The PCM is therefore the stochastic method used in the present study. To illustrate in detail the non-intrusive PCM implementation procedure, an uncertain differential equation system is considered:

Lðx; t; nðhÞÞ ¼ f ðx; t; nðhÞÞ where L is a differential operator, f is a source term, and n(h) is a vector of uncertain input variables as shown earlier. We consider

1898

P. Cinnella, S.J. Hercus / Computers & Fluids 39 (2010) 1893–1908

the case of two Gaussian random input variables (therefore n = 2), and a Lagrange interpolating polynomial of order Np = 1 (therefore p = 2, the order of the quadrature polynomial). As a result, the number of decoupled deterministic simulations required is PPCM = pn = 4. The PCM procedure is outlined as follows: According to the Askey scheme (see [39]), the Gaussian random input distributions require the use of a Hermite chaos polynomial. The roots of the Hermite polynomial correspond to collocation points hi for the uncertain parameter h, where i = 1, . . ., p. For each uncertain input parameter s = 1, 2, the p = 2 collocation points ns(hi) pﬃﬃﬃare determined from the transformation ns ðhÞ ¼ lin;s þ 2rin;s hi (Gaussian distribution case). The collocation points are now represented in terms of the parameters of the input distribution: ns ðhÞ Nðlin;s ; r2in;s Þ. A decoupled deterministic CFD simulation is carried out at the p collocation points for each uncertain parameter (n1, n2), permuted by a pn full factorial design:

3 n1 ðh1 Þ n2 ðh1 Þ 6 n ðh Þ n ðh Þ 7 2 1 7 6 1 2 D¼6 7 4 n1 ðh1 Þ n2 ðh2 Þ 5 2

n1 ðh2 Þ n2 ðh2 Þ where each row of the matrix D represents the input parameters for a single deterministic simulation. In this example, a total of PPCM = 4 deterministic simulations are carried out. Once all of the deterministic simulations have been completed, the output statistics can be calculated using the RMC method. A large number (M 10,000) of values of ns(h) are selected from the Gaussian distribution of each input parameter. These values are substituted into the multi-variable form of the Lagrange interpolation expression shown in Eq. (14). This linear equation rapidly returns M values of hi, which are used in Eq. (9) to reconstruct a distribution of M values for each output parameter (for example the pressure coefﬁcient or the mach number) for each grid point. The statistical moments of each output parameter distribution can ﬁnally be calculated arithmetically. 3.2. Preliminary validation of the probabilistic collocation method In order to conduct a preliminary evaluation of the PCM, the Buckley–Leverett (BL) equation, as presented by Leveque [40], is considered as a test case. This time-dependent model describes diphase ﬂuid ﬂow in a porous medium and is often applied to oil reservoir simulations. Initially, an underground oil source is tapped and oil ﬂows out under the existing high pressure. Once this initial ﬂow diminishes, water is pumped into the reservoir in order to recover the remaining oil. This process is known as secondary recovery. The BL system represents an appropriate test case for the stochastic PCM analysis, as the non-convex ﬂux may permit non-classic shock wave behavior which can also arise in the context of dense gases. The 1D oil–water interface is modeled by the scalar hyperbolic BL equation with a non-convex ﬂux function f(q):

f ðqÞ ¼

q2 q2 þ að1 qÞ2

ð15Þ

where a is a constant which relates the difference in properties between the two ﬂuids. Oil ﬂows less easily than water, so a < 1. The term q(x, t) is the fraction of ﬂuid that is water (the water saturation), so that 0 6 q 6 1 and 1 q(x, t) is the fraction of oil. Hence in the case of pure water q = 1 and for pure oil q = 0. The initial conditions are:

qðx; 0Þ ¼

1

if x ¼ 0

0:05 if x ¼ 1

ð16Þ

The solution of the Buckley–Leverett system is determined via an exact Riemann solver on a grid with cell size dx = 0.02. In this study, the solution of the time dependent problem is analyzed after an elapsed time of t = 0.7. The input parameter a is varied with a Gaussian input distribution with the statistics: al,in = 0.5, CV = r/ al,in = 1.5% (where CV represents the coefﬁcient of variation, a measure of the variance). Note that with these input parameters, the probability of a falling outside of the x 2 [0,1] range is negligible for this application. The stochastic analysis was implemented on the BL system using the PCM method, which required the use of a Hermite chaos polynomial according to the Askey scheme. The solution to the BL equation involves both a shock and a rarefaction wave, referred to as a compound wave [40]. The variation of the parameter a inﬂuences the position and size of the compound wave. Mean solutions obtained with increasing orders Np = 0 of the Lagrange polynomial chaos using the PCM method are compared in Fig. 2 to a Monte Carlo analysis with 10,000 deterministic simulations. Also included for comparison is the case where Np = 0, which represents the deterministic solution calculated at the input mean, denoted the DSIM. This is the solution obtained by classic CFD where random variations in input parameters are not taken into account. Outside the shock zone, the mean solutions provided by the PCM are in excellent agreement with the mean of the MC analysis. However, the mean shock value predicted by the PCM is spread abruptly over several cells, which generates a non-physical staircase appearance. This phenomenon arises due to the nature in which the mean solution is calculated in Eq. (12), where discrete deterministic solutions are multiplied by a weighting function. When this is applied to a large discontinuity with discrete displacement over several cells, the non-physical staircase form appears. Note that the mean solution obtained for Np = 8 is closer to the MC mean solution than the Np = 2 case, due to the fact that increasing the polynomial chaos order results in more deterministic solutions used in the calculation. For strong discontinuities, the MC stochastic analysis is locally superior to the PCM. However, the elevated computational cost of this method restricts its use to basic ﬂuid dynamics problems. In order to examine the convergence of the PCM, the mean l(xshock) and standard deviation r(xshock) of the shock position is calculated for increasing orders of the polynomial chaos (Fig. 3) using the exact expressions deﬁned by Eqs. (12) and (13). These values are compared to the values of the mean and standard deviation of xshock obtained by N = 10,000 MC simulations. Good agreement is observed between the output statistics provided by the PCM with a signiﬁcantly reduced number of deterministic simulations carried out. Although convergence of the PCM is not obtained to machine precision, the accuracy of the output statistics is sufﬁcient for use in a practical engineering context. In the case of complex CFD simulations with long calculation times, a good prediction of the response of the system to input parameter variation can be obtained with Np = 2, equivalent to PPCM = 3 deterministic simulations required at the collocation points for the BL test case with a single uncertain parameter. The reconstructed Monte Carlo (RMC) method is applied to the PCM results for Np = 2 in order to obtain a more complete image of the BL system response to the variation of the parameter a. Using M = 10,000 linear simulations on the Lagrange interpolating polynomial, the RMC is an inexpensive way to obtain additional statistical parameters. The method is shown to be reliable, with an excellent agreement (less than 0.1% difference) obtained between the mean value of xshock using Eq. (12) and the RMC with M = 10,000 simulations. Good agreement (approximately 2% difference) is observed between the two methods for the standard deviation of xshock. The RMC method becomes essential when more than one input parameter is varied, as the analytical expressions for the mean and variance of the output solutions become increas-

1899

P. Cinnella, S.J. Hercus / Computers & Fluids 39 (2010) 1893–1908

1.1 0.7

1.0 0.9

0.6

0.8

0.5

0.7

0.4

0.6

q

0.3

0.5 0.2

0.4 0.1

DSIM (Np= 0) Np= 2 Np= 8 MC10k

0.3 0.2

0.0 0.85

0.90

0.95

1.00

1.05

1.10

1.15

0.1 0.0 -0.4

-0.2

0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

x Fig. 2. Mean solutions of the Buckley–Leverett system, PCM and MC method comparison, with zoom of the shock region.

0.050

1.055

0.045 1.050

0.040 0.035 0.030

σ(xshock) PCM σ(xshock) MC10k μ(xshock) PCM μ(xshock) MC10k

1.040

1.035

0.025

σ(xshock)

μ(xshock)

1.045

0.020 0.015

1.030

0.010 1.025

0.005 1.020

0

1

2

3

4

5

6

7

0.000

Np Fig. 3. Convergence of the PCM analysis for the mean l and standard deviation r of the shock position xshock of the BL system.

ingly complex. The RMC is particularly useful if a complete set of output solutions is desired, from which higher order moments and other statistical information can be determined. Fig. 4 shows an envelope of possible output solutions of the BL system given a random variation of the parameter a, deﬁned by the 10th and 90th percentile solutions. Additionally, the mean and median solutions are compared to the DSIM in Fig. 4. Note that the staircase effect is still visible, as the reconstruction stage of the RMC method is carried out using a base of a discrete number of exact deterministic solutions obtained from the original PCM analysis. 3.3. Stochastic analysis of a transonic dense gas ﬂow over the NACA0012 airfoil The non-intrusive probabilistic collocation method (PCM) stochastic analysis was coupled to an existing dense gas code developed by Cinnella and Congedo [26] to model the ﬂow of a BZT type retrograde dense gas ﬂuid over one half of a symmetrical NACA0012 airfoil at 0° incidence. The upstream Mach number

was set to M1 = 0.95. The upstream thermodynamic conditions were ﬁxed in terms of critical point values at p/pc = 0.91 and T/ Tc = 1.02. In these conditions C = 0.43, which represents an operating point just outside the inversion zone. For this choice of the operating conditions, the ﬂow over the airfoil is transonic and shock waves are formed at both airfoil surfaces. Due to the symmetrical nature of the ﬂow, only half of the domain is simulated by explicitly enforcing a symmetry condition. A slip boundary condition was imposed on the surface of the airfoil and non-reﬂecting boundary conditions are used at the far-ﬁeld boundaries. Three structured half-C grids were created to discretize the ﬂuid domain: a coarse grid with 50 16 cells, a medium grid with 100 32 cells and a ﬁne grid with 200 64 cells. The far-ﬁeld boundary is located approximately 10 chords away from the airfoil. The medium grid is shown in Fig. 5. Initially, a set of deterministic simulations were carried out in order to evaluate the grid convergence and general ﬂow characteristics. The results of the deterministic simulation on the medium and ﬁne grids are visible in Fig. 6, where an oblique shock is visible close to the trailing edge of the proﬁle. Good agreement was observed between the medium grid and ﬁne grid solutions for the pressure coefﬁcient Cp, Mach number M and fundamental derivative C on the surface of the proﬁle. The oblique shock is sufﬁciently strong to approximately reestablish the far ﬁeld pressure value (Cp = 0 by deﬁnition). The grid convergence of the deterministic solution was studied by observing the drag coefﬁcient CD obtained from the three different grids. The drag coefﬁcient results are shown in Table 1, where the calculation times are obtained with an IntelÓ XeonÓ W3503 CPU with a 2.40 GHz clock speed. To reduce the computational cost with respect to ﬁne-grid computations while preserving a comparable accuracy on the computed aerodynamic coefﬁcients, the Richardson extrapolation (REX) method was applied to results from the coarse and medium grids. The details of this method are presented in Cinnella and Congedo [29] and will not be outlined here. The advantage of the REX method was a reduced calculation time with the same level of accuracy as the ﬁne grid solution. For this reason, the following computations are carried out using the REX method applied, using the coarse and medium grids. A stochastic analysis was carried out in order to examine operating point variability on the dense gas ﬂow simulation. The

1900

P. Cinnella, S.J. Hercus / Computers & Fluids 39 (2010) 1893–1908

1.1 0.7 1.0 0.9

0.6

0.8

0.5

0.7 0.4

q

0.6 0.3

0.5 0.4

0.2

DSIM Mean Solution Median 10th, 90th percentiles

0.3 0.2

0.1 0.0 0.8

0.1 0.0 -0.4

-0.2

0.0

0.2

0.4

0.6

0.8

1.0

1.2

0.9

0.9

1.0

1.0

1.1

1.1

1.2

1.4

x Fig. 4. Solution percentiles for the PCM analysis of the BL system, Np = 2, with zoom of the shock region.

Fig. 5. Structured C-grid for the upper half of a NACA0012 airfoil. Medium grid 100 32 cells, with close-up of the grid around the airfoil surface (right).

non-intrusive PCM approach was used to study the effects of Gaussian random variations about the upstream states of p/pc = 0.91 and T/Tc = 1.02, with the input coefﬁcients of variation (CV) set to 1.0% and 2.3% respectively. The Hermite polynomial was used in accordance with the Askey scheme and a reconstructed Monte Carlo (RMC) method with M = 10,000 linear simulations was used on the Lagrange interpolating polynomial in order to generate a complete set of output solutions. The PCM analysis was implemented with chaos polynomials of order Np = 2 and Np = 3 to enable a study of the stochastic method convergence properties. The location of the thermodynamic operating conditions of the collocation points can be observed on the p v diagram in Fig. 7. Results for the mean l, standard deviation r, and the coefﬁcient of variation CV of the drag coefﬁcient CD obtained by the PCM for Np = 2 and Np = 3 are shown in Table 2. As in the BL test case, good agreement is observed between the results of increasing polynomial order. In both cases the stochastic solution exhibits a high sensitivity to input parameter variation, due to the relatively signiﬁcant changes in physical ﬂow behavior for the different collocation points. Several collocation points lie within or very close to the inversion zone, resulting in a reduction in drag coefﬁcient due to the reduced magnitude of the compression shocks. As a result, the mean drag coefﬁcient obtained by the stochastic simulation is lower than in the case of the deterministic simulation at the input mean (DSIM). Using the complete set of output solutions obtained by the RMC method on the PCM solutions, a probability

density function of the drag coefﬁcient CD can be generated (Fig. 8). The distribution of the output CD is clearly non-Gaussian, due to the non-linearity of the dense gas system. In physical terms, Fig. 8 shows that a random Gaussian variation in an operating point close to the saturation curve results in a non-Gaussian distribution of the drag coefﬁcient which tends towards the low drag side. This is conﬁrmed with a skewness value of 2.2. Further evidence of the high sensitivity of the dense gas simulation to operating point variability is visible by the large size of the error bars in Fig. 9. The largest contributor of variability the position and magnitude of the oblique shock, which is directly affected by the reduction of the fundamental derivative C for operating points close to the inversion zone. Fig. 9 conﬁrms the base of the shock as a principal source of solution variation, with the largest values in standard deviation of the Mach number observed in this region. As in the BL test case, the uneven representation of the shock waves in Fig. 9 occurs due to the limited number of deterministic solutions used to model the discontinuity. These discrete individual deterministic solutions at the collocation points can be seen in Fig. 10 for the pressure coefﬁcient Cp over the NACA0012 proﬁle for the PCM analysis with Np = 3 (therefore PPCM = 16). This ﬁgure clearly shows that compression shocks can be avoided if the upstream operating point is located within the inversion zone. Using the complete set of output solutions generated by the RMC method with M = 10,000 linear simulations, it is possible to determine the percentiles of the output solution

1901

P. Cinnella, S.J. Hercus / Computers & Fluids 39 (2010) 1893–1908

(b) -2.0

(a)

-1.5

-1.0

0.0 1.

10

1.05

1.00

0.95

Cp

-0.5

1. 1

5

0.5

0.90

90 0.

0. 8 5

1.0 0 0. 8

1.5 0.0

(c)

(d)

0.2

0.4

0.2

0.4

x/c

0.6

0.8

1.0

0.6

0.8

1.0

0.8 0.7

1.2

0.7 0.6 1.0

Γ

M

0.6 0.5 0.5

0.8 0.4 0.4 0.6 0.0

0.2

0.4

x/c

0.6

0.8

1.0

0.3 0.0

x/c

Fig. 6. Deterministic simulation results for p/pc = 0.91 and T/Tc = 1.02. (a) Contours of the Mach number, M. (b) Pressure coefﬁcient Cp, (c) Mach number M, and (d) fundamental derivative C on the surface of the airfoil. Solid lines represent the medium grid solutions (100 32 cells) and dotted lines represent the ﬁne grid solutions (200 64 cells).

Table 1 Drag coefﬁcient CD results for the deterministic dense gas simulations. Grid

Number of cells

CD

Iterations

Residual (10x)

Calculation time (MM:SS)

Coarse Medium Fine Richardson extrapolation (REX)

50 16 100 32 200 64 Coarse + medium

0.0728 0.0724 0.0723 0.0723

5000 7000 10,000 4000 + 5000

6.0 9.0 5.6 10.2

01:22 08:18 46:25 07:11

(Fig. 10). Note that as the percentile calculation is carried out cell by cell, it does not produce a range of possible physical solutions, instead an envelope of possible Cp values is obtained at each point. Note that the median and the mean solutions are not coincident, suggesting the presence of non-Gaussian effects. Fig. 11 conﬁrms this non-Gaussian behavior with probability density functions of the pressure coefﬁcient at two points on the surface of the airfoil. Good overall convergence of the stochastic method is observed in Fig. 11 with the close accordance of the local solution distributions obtained by the Np = 2 and Np = 3 simulations. However, the increased solution variability close to the shock region (x/c = 0.9)

slightly reduces the quality of this convergence. In the following section of this work where a robust optimization of the airfoil proﬁle is carried out, a polynomial chaos order of Np = 2 is selected in order to reduce scalar computation time. In the near future, efﬁcient parallelization of the robust optimization procedure will permit the use of higher order polynomial chaos. Despite the limited number of deterministic solutions used compared with the classic MC method, the stochastic PCM analysis coupled with the RMC method can be used to provide critical uncertainty information for use in an engineering application with a reduced computational effort.

1902

P. Cinnella, S.J. Hercus / Computers & Fluids 39 (2010) 1893–1908

1.5

Saturation Curve Critical Isotherm Isotherms Γ= 0 Γ= 1 Np = 2 Np = 3

1.4 1.3 1. 0

p/pc

1.2

6

1.0 4

1.1

1.02

1.0

1.00

Input mean

0.9 0.8 0.7 0.6 0.5

1.0

1.5

2.0

2.5

3.0

v/vc Fig. 7. Collocation points corresponding to the input conditions for the stochastic PCM dense gas analysis.

Table 2 Results of the drag coefﬁcient CD, calculated using the REX method, for the PCM stochastic analysis. Np

Number of deterministic solutions

l(CD)

r(CD)

CV(CD) (%)

2 3 DSIM

9 16 1

0.0682 0.0685 0.0723

0.0120 0.0123 –

17.5 17.9 –

100

median (Np = 2)

median (Np = 3)

80

Np = 2 Np = 3

pdf(C D )

60

40

20

0 0.00

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.10

CD Fig. 8. Probability density function of the drag coefﬁcient CD, calculated using the REX method, with mean and median values for the PCM with Np = 2 and Np = 3.

4. Robust optimization under uncertain operating conditions 4.1. Robust optimization with a multi-objective genetic algorithm (MOGA) In spite of their relatively high computational cost, genetic algorithms have been successfully applied to aerodynamic shape optimization to a transonic dense ﬂows [27–29]. Evolutionary

optimization strategies employed in Pareto-based genetic algorithms are a ﬂexible and robust means of determining global optima of multi-point problems. Implementation is relatively straightforward without signiﬁcant code modiﬁcation, as only the evaluation of selected objective functions for each individual is required. In this study, the PCM stochastic analysis of transonic dense gas ﬂows presented earlier is coupled to an existing multiobjective genetic algorithm (MOGA), based on the non-dominated sorting algorithm proposed by Deb [41], Srinivas and Deb [42]. The problem is formulated as follows: ﬁnd a symmetric airfoil shape R such that the associated mean drag coefﬁcient l(CD) and its standard deviation r(CD) are minimized. The airfoil shape is constrained to satisfy the following conditions: (i) the coordinates of the leading edge and the trailing edge normalized by the airfoil chord are respectively (0, 0) and (1, 0); (ii) the thickness-to-chord ratio is ﬁxed and equal to 12%. Due to the ﬂow symmetry, only the upper surface of the airfoil is considered. The surface of the airfoil is parameterized with a Bezier curve, determined by specifying the coordinates of eight Bezier control points: the ﬁxed leading edge P0 and trailing edge P7; ﬁve control points Pk(xk, yk) regularly spaced along the chord; and a control point P1(0, y1) which ensures that the leading edge of the airfoil is tangent to the y-axis. These parameters are subsequently normalized to maintain a thicknessto-chord ratio of 12%. Using 11 distinct design parameters which vary between prescribed extrema, a family of airfoil shapes can be generated. In this study, the MOGA was used to generate a series of optimized symmetrical 2D airfoils for dense gas ﬂows, based on minimization of the mean and standard deviation of the drag coefﬁcient. 28 generations of 36 geometries were obtained, resulting in over 800 stochastic PCM analyses carried out (duplicated geometries were not recalculated). Each individual PCM analysis was calculated with a chaos polynomial of order Np = 2 (requiring nine deterministic simulations), which resulted in an overall calculation time of over 40 days on an IntelÓ XeonÓ W3503 2.4 GHz processor. The high computational cost of the genetic algorithm was mitigated by the use of the Richardson extrapolation method (REX) on each individual calculated. This solution extrapolation method, already applied to robust proﬁle optimization of dense gas ﬂows [29], accelerates convergence and increases the solution accuracy for each individual. Furthermore, Cinnella and Congedo [29] show that an increase in solution accuracy can improve MOGA convergence. The structure of the code for the robust optimization algorithm coupled with the stochastic dense gas analysis is strongly modular and is therefore very well suited for parallel operation. This is a critical point when considering potential industrial applications such as the design of ORC turbines, where more complex and costly CFD simulations would be required. A complete parallelization of the MOGA and of the PCM solver are warranted as forthcoming work. For comparison, a deterministic single parameter optimization was also carried out for 28 generations of 36 individuals, with only the drag coefﬁcient as the minimization function. With only a single deterministic simulation to be carried out for each individual at p/pc = 0.91 and T/Tc = 1.02, the overall scalar calculation time was reduced by a factor of 9, however this left the variation of the solution uncontrolled. Finally, several optimized proﬁles were selected from both the deterministic and stochastic optimization procedures and their performance veriﬁed on the medium grid (100 32 cells) with a higher polynomial chaos order (Np = 3). Since optimization based on an inviscid ﬂow model may lead to airfoil shapes with low wave drag but with poor viscous behavior [28], the airfoil performance was veriﬁed by applying the dense gas solver to the solution of the Reynolds-averaged Navier–Stokes (RANS) equations completed by the algebraic Baldwin–Lomax turbulence model. We refer to Cinnella and Congedo [26] for details about the viscous ﬂow solver. The turbulent re-

1903

P. Cinnella, S.J. Hercus / Computers & Fluids 39 (2010) 1893–1908

(b) -2.0

(a)

-1.5

Contours of μ(M)

-1.0

-0.5

0.0

1 1.05

0.95

0.20 0.18 0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02

Cp

0.95

1

σ(M)

0.9

1.1 5

0.9

1. 1

0.5

0.85

1.0

1.5 0.0

(c)1.4

(d)

1.3

0.2

0.4

x/c

0.6

0.8

1.0

0.9

0.8

1.2 0.7 1.1

Γ

M

0.6 1.0

0.5 0.9 0.4 0.8 0.3

0.7

0.6 0.0

0.2

0.4

x/c

0.6

0.8

1.0

0.2 0.0

0.2

0.4

x/c

0.6

0.8

1.0

Fig. 9. Stochastic PCM results on the medium grid (100 32) with Np = 3. (a) Contours of the mean Mach number shown by lines and shaded areas representing contours of standard deviation of the Mach number. (b) Pressure coefﬁcient Cp, (c) Mach number and (d) the fundamental derivative C on the surface of the airfoil x/c 2 [0, 1]. Solid lines represent the mean solution, dotted lines represent the DSIM solution. The error bars represent ±r for each parameter.

gime simulations were performed on a 160 64 half C-grid with a reﬁned boundary layer. 4.2. Results of the robust optimization procedure The mean and standard deviation of the drag coefﬁcients CD of all the solutions obtained from the stochastic robust optimization procedure are shown in Fig. 12, including the ﬁnal Pareto front and the stochastic results for the NACA0012 airfoil and the single parameter optimization (#753d). The stochastic robust optimization procedure provided a Pareto front in which several individuals exhibited improvements in the mean and standard deviation of the drag coefﬁcient over the standard NACA0012 airfoil. Very good agreement was observed between the Pareto front after 23 generations and the ﬁnal 28 generations, indicating that convergence of the MOGA procedure was achieved. The ﬁnal Pareto front was found to be non-convex, due to the non-linearity of the dense gas system. Three individual geometries were selected from the

front to be examined in further detail (Fig. 12). Individual #258 was selected as it was the individual with the lowest mean drag coefﬁcient and it exhibited a reduction of the standard deviation of the CD over the NACA0012 airfoil. Individual #767 was chosen as the individual with a similar mean drag coefﬁcient as the NACA0012, but with a lower value of r(CD), which represents an important improvement in stability. Finally, individual #791 was chosen as the most stable airfoil in the Pareto front, despite its poorer performance in terms of the mean drag coefﬁcient value. Note that for some industrial applications, a certain loss in mean performance may be acceptable if a highly stable solution is obtained. The ﬁnal compromise between performance and stability would require consideration of the eventual application of the ORC turbine, including a study of the availability, quality and cost of the initial energy source. Finally, using a single objective deterministic optimization, the lowest drag individual #753d (also shown on Fig. 12) is selected. This individual was optimized solely on the basis of a low deterministic drag coefﬁcient. The variability

1904

P. Cinnella, S.J. Hercus / Computers & Fluids 39 (2010) 1893–1908

-2.0

-1.5

-1.5

-1.0

-1.0

-0.5

-0.5

x/c = 0.9 x/c = 0.3

Cp

Cp

-2.0

0.0

0.0

0.5

0.5

Mean solution Deterministic solutions at collocation points

1.0

1.5 0.0

0.2

0.4

0.6

Mean th Median (50 percentile) th th 25 , 75 percentiles 10th, 90 th percentiles

1.0

0.8

1.5 0.0

1.0

x/c

0.2

0.4

0.6

0.8

1.0

x/c

Fig. 10. (left) Pressure coefﬁcient Cp over the NACA0012 proﬁle for each deterministic solution at the 16 collocation points with Np = 3. (right) Percentiles of the pressure coefﬁcient over the proﬁle for Np = 3.

6.0

μ (Np = 2)

μ (Np = 3)

1.8

μ (Np = 2)

μ (Np = 3)

5.5 1.6 1.4

4.5

Np = 2

4.0

Np = 3

Np = 2 Np = 3

1.2

3.5

pdf(Cp )

pdf(C p )

5.0

3.0 2.5 2.0

1.0 0.8 0.6

1.5

0.4

1.0 0.2

0.5 0.0 -1.8

-1.6

-1.4

-1.2

-1.0

-0.8

-0.6

C p (x/c = 0.3)

0.0 -4.0

-2.0

0.0

2.0

4.0

6.0

C p (x/c = 0.9)

Fig. 11. Probability distribution functions for the pressure coefﬁcient Cp at two points on the airfoil surface, for the Np = 2 and Np = 3 cases. The mean Cp values are shown by the vertical lines.

of its performance was not taken into account during the optimization procedure. The optimized geometries and the NACA0012 airfoil are shown in Fig. 13. In the case of the lowest drag airfoils (#258, #753d), the point of maximum thickness is located downstream (x/c 0.55) of the NACA0012 maximum thickness point. The more stable airfoil proﬁles (#767, #791) have a maximum thickness in approximately the same region as the NACA0012 airfoil (x/c 0.30) but exhibit a more shallow trailing edge gradient compared to the other proﬁles. The only source of drag in the non-viscous simulations is the trailing edge compression shock, so any modiﬁcation of the location or magnitude of this shock affects the overall drag coefﬁcient for the airfoil. A veriﬁcation of the stochastic performance for Np = 3 of the dense gas ﬂow over the optimized airfoils is shown in Fig. 14. Evidence of the instability of individuals #258 and #753d is observed

by the small region of elevated values of r(M) at the base of the shock region. The the steep trailing edge gradient of these individuals induces an oblique and unstable trailing edge shock. In comparison, the most stable airfoil (#791) exhibits smaller values of r(M) spread over a larger region at the base of the shock. The shallow trailing edge gradients of these more stable individuals (#767 and #791) promote the formation of stronger shock waves. Although this induces a larger drag coefﬁcient, the stability of the shock is improved and the variation of the drag coefﬁcient is reduced. Good convergence of the stochastic PCM analysis was observed in the comparison of the drag coefﬁcient data for Np = 2 and Np = 3 (Table 3). The hierarchy of stability was preserved in the higher order case, with individuals #258 and #753d exhibiting the highest CV values. The stochastic convergence of the mean and standard deviation of the drag coefﬁcient shown in Table 3 val-

P. Cinnella, S.J. Hercus / Computers & Fluids 39 (2010) 1893–1908

0.018

MOGA airfoils, 28 Generations Pareto Front Individuals NACA0012

0.017

0.015

#753d

0.016

0.014 #753d

σ(CD )

0.013 0.012

NACA0012

0.011 0.010 0.009 #258

0.008 #767

0.007 0.006 0.005 0.050

#791

0.055

0.060

0.065

0.070

μ(C D )

0.075

0.080

0.085

Fig. 12. Mean l and standard deviation r values of the drag coefﬁcient CD of the ﬁnal population generated by the two parameter PCM stochastic robust optimization problem (REX grid method, stochastic polynomial order Np = 2). The Pareto front individuals after 28 generations of 36 individuals are encircled. Stochastic results for the NACA0012 and the single parameter optimal individual (#753d) are also shown.

0.07 0.06 0.05

y/c

0.04 0.03

NACA0012 #258 #767 #791 #753d

0.02 0.01 0.00 0.0

0.2

0.4

x/c

0.6

0.8

1.0

Fig. 13. Geometries of selected optimized airfoils from the Pareto front. Also shown is the geometry obtained from the deterministic optimization (#753d) and the NACA0012 airfoil. Aspect ratio not preserved.

idates the selection of Np = 2 in the MOGA to reduce the overall scalar computational effort. As a ﬁnal validation of the performance of the optimized individuals, a series of simulations in the turbulent regime was carried out using the Baldwin–Lomax turbulence model on a 160 64 cell grid with a reﬁned boundary layer. The performance of the optimized airfoils is summarized in Table 4, where the mean, standard deviation and DSIM drag coefﬁcients are presented. The total viscous drag coefﬁcient included both wave drag and viscous drag components. However, the viscous drag component was relatively minor, typically contributing less than 7% to the total drag coefﬁcient. Consequently, the results obtained by simulations in the turbulent regime were relatively similar to the inviscid case. As in the inviscid case for the NACA0012 airfoil examined earlier, all of the mean drag coefﬁcients obtained by the stochastic simulations are smaller than the DSIM results for all of the optimized proﬁles considered. All of the airfoils show a high sensitivity to input parameter variation, with the most stable airfoil #791 only achieving minimum CV(CD) of 7.7% in the inviscid case and 7.9% in the turbu-

1905

lent regime. Once again, the hierarchy of stability is preserved, with individuals #258 and #753d exhibiting the highest levels of instability. The results of stochastic simulations of dense gas ﬂow over the optimized airfoils presented in Tables 3 and 4 show the trend that a highly stable airfoil can be obtained with a shallow trailing edge gradient. The price of this improved stability is an increased shock wave intensity and a subsequent increase in the drag coefﬁcient. The characteristics of the turbulent dense gas ﬂow over the most stable (#791), most unstable (#753d), and the NACA 0012 airfoils are presented in Fig. 15, via contours of the mean and standard deviation of the Mach numbers. Superimposed on the contour plot of the mean Mach number values is a streamline representing the approximate extent of the separated recirculation bubble (when this is present). In the case of the unstable but low drag #753d airfoil, a small region of intense variation of the Mach number is observed at the base of the trailing edge shock. Additionally, the steep trailing edge gradient of this airfoil promotes the formation of a large region of separated ﬂow. These two effects are largely responsible for the poor stability performance of this airfoil. On the other hand, the stabilizing effects of a shallow trailing edge gradient are clear for the #791 airfoil. The ﬂow does not separate at the trailing edge, and the values of the standard deviation of the Mach number are lower than in the #753d airfoil case and are spread over a larger region at the trailing edge. However, the cost of this increased stability is an increase in drag coefﬁcient due to the formation of a stronger shock wave. The elevated cost of the two parameter stochastic optimization is justiﬁed when a compromise between a low mean drag coefﬁcient and a small solution variation is desired. The airfoil #753d obtained from the single parameter optimization exhibits very high CV values of 24.5% (inviscid) and 19.9% (turbulent regime) for the same order of polynomial chaos, as the variation of the solution is not taken into account during the optimization procedure. Note that for all of the simulations carried out in the present work, only the upper surface of the airfoils were considered. In this conﬁguration, it is possible for a detached ﬂuid ﬂow to reattach to the virtual boundary downstream of the airfoil proﬁle. In the full proﬁle case, this is not possible and an unstable recirculation zone may be created. As a result, the airfoils with the lowest drag coefﬁcients (#258, #753d) may eventually exhibit poorer drag performance in an unsteady simulation of the full proﬁle. Timedependent simulations in both two and three dimensions are required to conﬁrm this phenomenon. The robust optimization procedure based on a genetic algorithm coupled to a stochastic uncertainty quantiﬁcation method proves to be a powerful design tool. The robust optimization provides detailed engineering information of system response to input parameter variation, which will allow increased design conﬁdence. In the case of an ORC turbine where the input energy source is highly variable, predictable knowledge of the variation of system response is critical for evaluating system feasibility.

5. Conclusions In the present work, a series of optimized airfoils was developed for dense gas ﬂows under uncertain operating conditions using a robust optimization procedure based on a multi-objective genetic algorithm (MOGA). This type of analysis is essential in improving the feasibility of organic Rankine cycle (ORC) turbines, which are typically designed to recover energy from variable sources such as waste heat from industrial processes. Initially, the dense gas system response to input parameter variation was quantiﬁed using a stochastic polynomial chaos (PC) approach developed by Loeven et al. [19,24] known as the

1906

P. Cinnella, S.J. Hercus / Computers & Fluids 39 (2010) 1893–1908

(a)

(b) #258

#767

0.8 5

0.95 1 1.05 1.1 1.1 5

0. 9

5 0.8

1.05 1 1. 1 . 1 5

0.9

0. 9

0.8 5

0.95

1. 05

C ontours of μ (M)

0.9

1

C ontou urs ooff μ (M)

0.95

0.95

0.20 0.18 0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02

1

σ (M) σ (M (M) 0.20 0.18 0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02

2 1.

(d)

(c)

#75 7 3d

0.20 0.18 0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02

1

0.95

C on ntou urs ooff μ (M)

0.9

0.9

1.05

1.1 1.1 5

2

0.8 5

0.8 5

0.20 0.18 0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02

0. 9

1 1.1 .1 5

1 1 . 05

0. 9

1.

σ (M)

0.95

0.95

0. 9

C ontoour s ooff μ (M)

σ (M)

1 1.05

0.95

1

#7 1 #79

0. 8 5

Fig. 14. Contours of the mean l and standard deviation r of the Mach number M for the optimized airfoils. Stochastic study with Np = 3. Medium grid 100 32. Airfoils (a) #258. (b) #767 (c) #791 and (d) #753d.

Table 3 Drag coefﬁcients (REX grid) for optimized airfoils, Np = 2 and Np = 3. Airfoil geometry

NACA0012 #258 #767 #791 #753d

CD(DSIM)

Np = 2

0.0722 0.0595 0.0717 0.0803 0.0599

Np = 3

l(CD)

r(CD)

CV(CD) (%)

l(CD)

r(CD)

CV(CD) (%)

0.0682 0.0573 0.0687 0.0782 0.0565

0.0120 0.0104 0.0083 0.0060 0.0138

17.5 18.7 12.1 7.7 24.5

0.0685 0.0566 0.0690 0.0782 0.0567

0.0123 0.0116 0.0095 0.0079 0.0130

17.9 20.5 13.8 10.1 22.9

Table 4 Viscous drag coefﬁcients (turbulent grid, 160 64 cells) for optimized airfoils, Np = 2. Airfoil geometry

CD(DSIM)

l(CD)

r(CD)

CV(CD) (%)

NACA0012 #258 #767 #791 #753d

0.0704 0.0544 0.0730 0.0821 0.0467

0.0669 0.0568 0.0707 0.0804 0.0526

0.0108 0.0099 0.0085 0.0064 0.0104

16.2 17.5 12.0 7.9 19.9

probabilistic collocation method (PCM). The PCM was initially validated with the known test case of diphasic ﬂuid ﬂow in a porous medium, modeled by the Buckley–Leverett (BL) equation. This preliminary test case using the exact solution of the BL equation allowed a comparison of the PCM method against a classic Monte Carlo (MC) approach, when applied to a system exhibiting nonclassical shock wave behavior. The PCM is able to provide statistical output solutions which are in excellent agreement with the MC analysis at a fraction of the computational cost. Convergence of the

1907

P. Cinnella, S.J. Hercus / Computers & Fluids 39 (2010) 1893–1908

(a)

N CA NA C 0012

N CA NA C 0012 μ (M (M): 0.00 0.25 0.50 0.75 1.00 1.25

(b)

σ (M (M): 0.02 0.10 0.18 0.26 0.34

#753 7 d

#753d

σ (M (M): 0.02 0.10 0.18 0.26 0.34

μ (M (M): 0.00 0.25 0.50 0.75 1.00 1.25

(c)

#791

#791 μ (M (M): 0.00 0.25 0.50 0.75 1.00 1.25

σ (M ( ): 0.02 0.10 0.18 0.26 0.34

Fig. 15. Contours of the mean (left) and standard deviation (right) of the Mach number for the airfoils (a) NACA 0012, (b) #753d [most unstable], (c) #791 [most stable]. The streamline shows the approximate extent of the separated ﬂow. NB: aspect ratio not preserved.

PCM to an appropriate engineering accuracy is obtained rapidly with a polynomial chaos order Np of less than 7. Using an inexpensive reconstructed MC analysis (RMC) coupled to the PCM, a full set of output solutions can be generated which enable the calculation of higher order moments and other statistical information. The validated PCM algorithm with a subsequent RMC is then applied to an existing solver Cinnella and Congedo [26] of an inviscid transonic dense gas ﬂow over half of a symmetrical NACA0012 airfoil at M1 = 0.95. The working ﬂuid used is pf-perhydroﬂuorene (C13F22, commercially known as PP10) with the thermodynamic properties modeled by the Martin–Hou equation of state. The aerodynamic performance of the airfoil at 0° incidence is examined with Gaussian random variations in the upstream pressure and temperature. Using a small number of deterministic solutions, the PCM determined that the dense gas system was highly sensitive to input parameter variation, due to the large changes in physical ﬂow behavior close the the inversion zone of PP10. In this zone where the fundamental derivative of gas dynamics C < 0, compression shock waves can be suppressed as a consequence of the entropy inequality. Furthermore, the set of statistical output solutions was found to be non-Gaussian which tends towards the low drag side, due to the non-linearity of the dense gas system. Finally, the PCM analysis of the dense gas system with Np = 2 was coupled with an existing MOGA based robust optimization procedure, which generated a series of optimized 2D airfoils for dense gas ﬂows based on minimization of the mean and standard deviation of the drag coefﬁcient. Several individual optimized geometries were selected from the ﬁnal Pareto front which exhibited improvements in mean performance and/or stability over the NACA0012 airfoil. A posteriori testing of the the stochastic perfor-

mance of these individuals was veriﬁed for a polynomial chaos order Np = 3 on a standard medium grid (100 32 cells) grid in the inviscid case, and for Np = 2 on a grid with a reﬁned boundary layer (160 64 cells) in the turbulent regime simulation. The viscous drag component was found to be a relatively minor contributor to the overall drag coefﬁcient, whereas the wave drag of the trailing edge shock proved to be the dominant factor. The position, intensity and sensitivity of the trailing edge shock to ﬂuctuations in the operating conditions were modiﬁed by the geometry of the optimized individuals, which in turn strongly inﬂuenced the drag coefﬁcient. In particular, airfoils with a shallow trailing edge gradient were found to exhibit increased stability at the cost of an increased mean drag coefﬁcient. The two parameter stochastic robust optimization was compared to a standard deterministic optimization procedure where only the mean drag coefﬁcient was optimized. The elevated cost of the two parameter stochastic optimization is justiﬁed when both a low mean drag coefﬁcient and a small solution variation are desired. Several methods exist to mitigate the elevated computational cost of a multi-parameter optimization procedure. Large scale parallelization of the PCM analysis and the MOGA codes could be implemented with relative ease, as each collocation point of each individual in a generation could be calculated concurrently. Other methods to reduce calculation cost include surrogate models such as the artiﬁcial neural network (ANN), already employed to shape optimization of dense gas ﬂows [28]. Following the implementation of these computational cost reduction methods, multiobjective robust optimization procedures could help to develop the next generation of highly efﬁcient turbomachinery. The present study has shown that a robust optimization based on a genetic

1908

P. Cinnella, S.J. Hercus / Computers & Fluids 39 (2010) 1893–1908

algorithm coupled to a stochastic analysis proves to be a powerful design tool. The robust optimization provides detailed engineering information of system response to operating point variation, which could ultimately lead to reﬁnements in the robust design of ORC turbines. References [1] Horen J, Talonpoika T, Larjola J, Siikonen T. Numerical simulation of real-gas ﬂow in a supersonic turbine nozzle ring. J Eng Gas Turb Power 2002;124(2):395–403. [2] Brown B, Argrow B. Application of Bethe–Zel’dovich–Thompson ﬂuids in organic Rankine cycle engines. J Propul Power 2000;16(6):1118–23. [3] Monaco J, Cramer M, Watson L. Supersonic ﬂows of dense gases in cascade conﬁgurations. J Fluid Mech 1997;330:31–59. [4] Angelino G, Invernizzi C. Supercritical heat pump cycles. Int J Refrig 1994;17(8):543–54. [5] Zamﬁrescu C, Dincer I. Performance investigation of high-temperature heat pumps with various BZT working ﬂuids. Thermochim Acta 2009;488:66–77. [6] Kirillov N. Analysis of modern natural gas liquefaction technologies. Chem Petrol Eng 2004;40(7–8):401–6. [7] Lambrakis K, Thompson P. Existence of real ﬂuids with a negative fundamental derivative C. Phys Fluids 1972;15(5):933–5. [8] Thompson P, Lambrakis K. Negative shock waves. J Fluid Mech 1973;60:187–208. [9] Cramer M. Negative nonlinearity in selected ﬂuorocarbons. Phys Fluids 1989;1(11):1894–7. [10] Colonna P, Nannan N, Guardone A, van der Stelt T. On the computation of the fundamental derivative of gas dynamics using equations of state. Fluid Phase Equilibr 2009;286(1):43–54. [11] Thompson P. A fundamental derivative in gas dynamics. Phys Fluids 1971;14:1843–9. [12] Cramer M, Kluwick A. On the propagation of waves exhibiting both positive and negative nonlinearity. J Fluid Mech 1984;142:9–37. [13] Cramer M, Tarkenton G. Transonic ﬂows of Bethe–Zel’dovich–Thompson Fluids. J Fluid Mech 1992;240:197–228. [14] Maizza V, Maizza A. Unconventional working ﬂuids in organic rankine-cycles for waste energy recovery systems. Appl Therm Eng 2001;21:381–90. [15] Hung T, Shai T, Wang S. A review of organic rankine cycles (ORCs) for the recovery of low-grade waste heat. Energy 1997;22(7):661–7. [16] Cinnella P, Congedo P. Aerodynamic performance of transonic Bethe– Zel’dovich–Thompson ﬂows past an airfoil. AIAA J 2005;43:370–8. [17] AIAA. Guide for the veriﬁcation and validation of computational ﬂuid dynamics simulations. G-077-1998, American Institute of Aeronautics & Astronautics; 1998. ISBN 1563472856. [18] Cinnella P, Congedo P, Pediroda L. Quantiﬁcation of thermodynamic uncertainties in real-gas ﬂows. Int J Eng Syst Modell Simul 2009;2(1-2):12–24. [19] Loeven A, Witteveen J, Bijl H. Efﬁcient uncertainty quantiﬁcation using a twostep approach with chaos collocation. In: ECCOMAS CFD, TU Delft, Netherlands; 2006.

[20] Walters R, Huyse L. Uncertainty analysis for ﬂuid mechanics with applications. ICASE Report No. 2002-1 (NASA/CR-2002-211449); 2002. [21] INRIA. TAPENADE: on-line automatic differentiation engine; 2002.