Robust optimization of dense gas flows under uncertain operating conditions

Robust optimization of dense gas flows under uncertain operating conditions

Computers & Fluids 39 (2010) 1893–1908 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 ...

2MB Sizes 0 Downloads 27 Views

Computers & Fluids 39 (2010) 1893–1908

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 flows 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 quantification 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 profiles for transonic inviscid flows 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 defined 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 influential in the dynamic behavior of the fluid. 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 fluids characterized by particularly high molecular complexity have been theoretically predicted [7–10] to display a region of negative values of the fundamental derivative of fluid 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 flows of BZT type fluids, 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 fluids could potentially have a significant 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.compfluid.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 fluid, ORC turbines use an organic fluid such as hydrocarbons, silicon oils or other organic refrigerants. The low critical temperature, high density and elevated heat capacities of these fluids 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 fluids eliminates the risk of condensate at the turbine outlet, without heating the working fluid 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 configuration 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 efficiency. However, this small temperature jump (and therefore enthalpy jump) diminishes power output and first law efficiency, 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 fluids 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 fluids 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 efficiency of the turbine is increased at the cost of diminished global power output and overall cycle efficiency. 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 flows in ORC turbines, a compromise is suggested in Cinnella and Congedo [16]. The proposition is to allow the dense gas flow 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 efficiency. 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 flow 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 fluids currently considered, the inversion zone, where the turbine isentropic efficiency 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 flow 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 difficult 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 fluids. 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 final 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 fluid 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 final 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 difficult to achieve given fluctuating real world conditions. Uncertainty quantification consists of measuring the system response to random and unknown variations in input parameters using stochastic analysis. The AIAA [17] provides a definition of uncertainty as a ‘‘potential deficiency 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 confidence, reduced risks, improved safety and refinement of system operation range through reduction of over-optimization. Dense gas flows 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 flow 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 efficiency 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 fluid flow 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 modification 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 significantly 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 first and second order sensitivity derivatives, which can be difficult 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 flows. It is stated that PC methods can efficiently 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 significant modifications of the underlying deterministic flow solver. The intrusive approach may result in improved computational efficiency for highly dimensioned stochastic problems, however it dramatically decreases flexibility 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 flows, 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 simplified geometry, namely an airfoil placed into a dense gas stream with randomly varying thermodynamic conditions. The study is divided into two main sections. The first section considers uncertainty quantification applied to transonic dense gas flows 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 fluid flow in a porous medium, modeled by the Buckley–Leverett (BL) equation. This test model involves a non-convex flux 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 flow over half of a symmetrical NACA0012 airfoil at M1 = 0.95. The working fluid used is pf-perhydrofluorene (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 flows [27–29]. The robust optimization procedure generates a series of optimized 2D airfoils for dense gas flows, based on minimization of the mean and standard deviation of the drag coefficient. 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 profile optimization of dense gas flows [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 verified with both deterministic and stochastic simulations. Finally, a posteriori testing of the optimized individuals is carried out in the turbulent flow regime using the Baldwin–Lomax turbulence model.

2. Governing equations and flow solver 2.1. Physics of dense gas flows 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) fluids. In a specific 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 fluid density and s the entropy of the fluid. 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 specific heat ratio of the fluid) is strictly greater than 1 due to thermodynamic stability requirements. For the majority of classical working fluids, an isentropic compression leads to an increase in the speed of sound. For a BZT fluid, if C < 1, then ð@[email protected] qÞs < 0, and the flow 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 fluid property across the shock, v = 1/q is the specific 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 perfluorocarbon pf-perhydrofluorene, 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 fluid: the liquid/vapor coexistence curve exhibits a positive slope up to near-critical conditions. Numerical simulation of BZT flows over isolated airfoils and wings [30,31] has shown that with a negative upstream value of C, the flow 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 fluids in turbine design could lead to a potential increase in efficiency of approximately 3–8%, with advancement in fluid dynamic design and non-classical effect exploitation. 2.2. Dense gas model and flow solver Dense gas flows are governed by the equations for equilibrium for non-reacting flows. 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 flux density:

f ¼ ðqv ; qI; qvv; qv HÞT where v is the velocity vector, E the specific total energy, H = E + p/q the specific 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 specific 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 specific 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 specified. 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:



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 finite volume scheme for structured multi-block meshes of third-order accuracy, which allows the computation of flows 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 simplifies 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 fluxes 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 efficiently drive the solution to the steady state. For external flows, non-reflecting boundary conditions based on a multidimensional method of characteristics are applied at the far-field 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 quantification for dense gas flows

with Tc the critical temperature. The gas dependent coefficients 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 specific 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 quantification 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 final 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 field u subjected to random input variable fluctuations n(h), the decomposition is given by:

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

P GPC X

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

ð8Þ

rin correspond respectivelyptoffiffiffi 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 fluctuation 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 find 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 efficient 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 modifications to the classic method. The first modification 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 modification 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 modified 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

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

ð11Þ

Additionally, in the case of Gauss–Hermite quadrature, the followpffiffiffi 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 modified. 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 definitions 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 significant 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 flexibility. The PCM is therefore particularly suited to the study of both complex CFD simulations, without requiring modification 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) pffiffiffiare 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 coefficient or the mach number) for each grid point. The statistical moments of each output parameter distribution can finally 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 fluid flow in a porous medium and is often applied to oil reservoir simulations. Initially, an underground oil source is tapped and oil flows out under the existing high pressure. Once this initial flow 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 flux 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 flux 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 fluids. Oil flows less easily than water, so a < 1. The term q(x, t) is the fraction of fluid 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 coefficient 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 influences 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 fluid 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 defined 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 significantly 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 sufficient 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, defined 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 flow 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 flow of a BZT type retrograde dense gas fluid 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 fixed 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 flow over the airfoil is transonic and shock waves are formed at both airfoil surfaces. Due to the symmetrical nature of the flow, 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-reflecting boundary conditions are used at the far-field boundaries. Three structured half-C grids were created to discretize the fluid domain: a coarse grid with 50  16 cells, a medium grid with 100  32 cells and a fine grid with 200  64 cells. The far-field 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 flow characteristics. The results of the deterministic simulation on the medium and fine grids are visible in Fig. 6, where an oblique shock is visible close to the trailing edge of the profile. Good agreement was observed between the medium grid and fine grid solutions for the pressure coefficient Cp, Mach number M and fundamental derivative C on the surface of the profile. The oblique shock is sufficiently strong to approximately reestablish the far field pressure value (Cp = 0 by definition). The grid convergence of the deterministic solution was studied by observing the drag coefficient CD obtained from the three different grids. The drag coefficient 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 fine-grid computations while preserving a comparable accuracy on the computed aerodynamic coefficients, 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 fine 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 flow 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 coefficients 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 coefficient of variation CV of the drag coefficient 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 significant changes in physical flow behavior for the different collocation points. Several collocation points lie within or very close to the inversion zone, resulting in a reduction in drag coefficient due to the reduced magnitude of the compression shocks. As a result, the mean drag coefficient 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 coefficient 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 coefficient which tends towards the low drag side. This is confirmed 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 confirms 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 coefficient Cp over the NACA0012 profile for the PCM analysis with Np = 3 (therefore PPCM = 16). This figure 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 coefficient 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 fine grid solutions (200  64 cells).

Table 1 Drag coefficient 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 confirms this non-Gaussian behavior with probability density functions of the pressure coefficient 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 profile is carried out, a polynomial chaos order of Np = 2 is selected in order to reduce scalar computation time. In the near future, efficient 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 coefficient 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 coefficient 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 flows [27–29]. Evolutionary

optimization strategies employed in Pareto-based genetic algorithms are a flexible and robust means of determining global optima of multi-point problems. Implementation is relatively straightforward without significant code modification, as only the evaluation of selected objective functions for each individual is required. In this study, the PCM stochastic analysis of transonic dense gas flows 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: find a symmetric airfoil shape R such that the associated mean drag coefficient 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 fixed and equal to 12%. Due to the flow 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 fixed leading edge P0 and trailing edge P7; five 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 flows, based on minimization of the mean and standard deviation of the drag coefficient. 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 profile optimization of dense gas flows [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 coefficient 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 profiles were selected from both the deterministic and stochastic optimization procedures and their performance verified on the medium grid (100  32 cells) with a higher polynomial chaos order (Np = 3). Since optimization based on an inviscid flow model may lead to airfoil shapes with low wave drag but with poor viscous behavior [28], the airfoil performance was verified 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 flow 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 coefficient 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 refined boundary layer. 4.2. Results of the robust optimization procedure The mean and standard deviation of the drag coefficients CD of all the solutions obtained from the stochastic robust optimization procedure are shown in Fig. 12, including the final 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 coefficient over the standard NACA0012 airfoil. Very good agreement was observed between the Pareto front after 23 generations and the final 28 generations, indicating that convergence of the MOGA procedure was achieved. The final 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 coefficient 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 coefficient 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 coefficient value. Note that for some industrial applications, a certain loss in mean performance may be acceptable if a highly stable solution is obtained. The final 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 coefficient. 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 coefficient Cp over the NACA0012 profile for each deterministic solution at the 16 collocation points with Np = 3. (right) Percentiles of the pressure coefficient over the profile 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 coefficient 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 profiles (#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 profiles. The only source of drag in the non-viscous simulations is the trailing edge compression shock, so any modification of the location or magnitude of this shock affects the overall drag coefficient for the airfoil. A verification of the stochastic performance for Np = 3 of the dense gas flow 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 coefficient, the stability of the shock is improved and the variation of the drag coefficient is reduced. Good convergence of the stochastic PCM analysis was observed in the comparison of the drag coefficient 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 coefficient 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 coefficient CD of the final 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 final 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 refined boundary layer. The performance of the optimized airfoils is summarized in Table 4, where the mean, standard deviation and DSIM drag coefficients are presented. The total viscous drag coefficient 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 coefficient. 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 coefficients obtained by the stochastic simulations are smaller than the DSIM results for all of the optimized profiles 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 flow 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 coefficient. The characteristics of the turbulent dense gas flow 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 flow. 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 flow 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 coefficient due to the formation of a stronger shock wave. The elevated cost of the two parameter stochastic optimization is justified when a compromise between a low mean drag coefficient 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 configuration, it is possible for a detached fluid flow to reattach to the virtual boundary downstream of the airfoil profile. In the full profile case, this is not possible and an unstable recirculation zone may be created. As a result, the airfoils with the lowest drag coefficients (#258, #753d) may eventually exhibit poorer drag performance in an unsteady simulation of the full profile. Timedependent simulations in both two and three dimensions are required to confirm this phenomenon. The robust optimization procedure based on a genetic algorithm coupled to a stochastic uncertainty quantification 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 confidence. 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 flows 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 quantified 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 coefficients (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 coefficients (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 fluid flow 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 flow. 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 flow over half of a symmetrical NACA0012 airfoil at M1 = 0.95. The working fluid used is pf-perhydrofluorene (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 flow 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 flows based on minimization of the mean and standard deviation of the drag coefficient. Several individual optimized geometries were selected from the final 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 verified 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 refined 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 coefficient, 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 fluctuations in the operating conditions were modified by the geometry of the optimized individuals, which in turn strongly influenced the drag coefficient. In particular, airfoils with a shallow trailing edge gradient were found to exhibit increased stability at the cost of an increased mean drag coefficient. The two parameter stochastic robust optimization was compared to a standard deterministic optimization procedure where only the mean drag coefficient was optimized. The elevated cost of the two parameter stochastic optimization is justified when both a low mean drag coefficient 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 artificial neural network (ANN), already employed to shape optimization of dense gas flows [28]. Following the implementation of these computational cost reduction methods, multiobjective robust optimization procedures could help to develop the next generation of highly efficient 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 refinements in the robust design of ORC turbines. References [1] Horen J, Talonpoika T, Larjola J, Siikonen T. Numerical simulation of real-gas flow 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 fluids in organic Rankine cycle engines. J Propul Power 2000;16(6):1118–23. [3] Monaco J, Cramer M, Watson L. Supersonic flows of dense gases in cascade configurations. 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] Zamfirescu C, Dincer I. Performance investigation of high-temperature heat pumps with various BZT working fluids. 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 fluids 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 fluorocarbons. 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 flows of Bethe–Zel’dovich–Thompson Fluids. J Fluid Mech 1992;240:197–228. [14] Maizza V, Maizza A. Unconventional working fluids 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 flows past an airfoil. AIAA J 2005;43:370–8. [17] AIAA. Guide for the verification and validation of computational fluid dynamics simulations. G-077-1998, American Institute of Aeronautics & Astronautics; 1998. ISBN 1563472856. [18] Cinnella P, Congedo P, Pediroda L. Quantification of thermodynamic uncertainties in real-gas flows. Int J Eng Syst Modell Simul 2009;2(1-2):12–24. [19] Loeven A, Witteveen J, Bijl H. Efficient uncertainty quantification using a twostep approach with chaos collocation. In: ECCOMAS CFD, TU Delft, Netherlands; 2006.

[20] Walters R, Huyse L. Uncertainty analysis for fluid mechanics with applications. ICASE Report No. 2002-1 (NASA/CR-2002-211449); 2002. [21] INRIA. TAPENADE: on-line automatic differentiation engine; 2002. . [22] Lucor D, Xiu D, Su C, Karniadakis G. Predictability and uncertainty in CFD. Int J Numer Meth Fl 2003;43:483–505. [23] Knio O, LeMaître O. Uncertainty propagation in CFD using polynomial chaos decomposition. Fluid Dyn Res 2006;38:616–40. [24] Loeven A, Witteveen J, Bijl H. Probabilistic collocation: an efficient nonintrusive approach for arbitrarily distributed parametric uncertainties. In: 45th AIAA aerospace sciences meeting and exhibit. AAIA Paper 2007-317, Reno, Nevada, USA; 2007. [25] Abgrall R. A simple, flexible and generic deterministic approach to uncertainty quantification in non-linear problems: application to fluid flow problems. Tech. Rep., INRIA; 2008. . [26] Cinnella P, Congedo P. Numerical solver for dense gas flows. AIAA J 2005;43:2458–61. [27] Congedo P, Corre C, Cinnella P. Airfoil shape optimization for transonic flows of Bethe–Zel’dovich–Thompson fluids. AIAA J 2007;45:1303–16. [28] Cinnella P, Congedo P. Optimal airfoil shapes for viscous transonic flows of Bethe–Zel’dovich–Thompson fluids. Comput Fluids 2008;37:250–64. [29] Cinnella P, Congedo P. GA-hardness of aerodynamic optimization problems: analysis and proposed cures. In: 18th AIAA computational fluid dynamics conference. AAIA Paper 2007-3828, Miami, Florida, USA; 2007. [30] Cinnella P, Congedo P. Inviscid and viscous aerodynamics of dense gases. J Fluid Mech 2007;580:179–217. [31] Cinnella P. Transonic flows of dense gases over finite wings. Phys Fluids 2008;20:1–17. 046103. [32] Cinnella P, Congedo P, Laforgia D. Transonic flows of BZT fluids through turbine cascades. In: Proceedings of the 3rd ICCFD conference. Toronto, Canada, Berlin: Springer-Verlag; 2004, ISBN 3-540-31800-3. p. 227–32. [33] Congedo P, Cinnella P, Corre C. Shape optimization for dense gas flows in turbine cascades. In: Proceedings of ICCFD 4, Ghent, Belgium; 2006. [34] Martin J, Hou Y. Development of an equation of state for gases. AIChE J 1955;1:142–51. [35] Emanuel G. Assessment of the MartinHou equation for modeling a nonclassical fluid. J Fluids Eng 1996;116:883–4. [36] Guardone A, Vigevano L, Argrow B. Assessment of thermodynamic models for dense gas dynamics. Phys Fluids 2004;16:3878–87. [37] Jameson A, Schmidt W, Turkel E. Solutions of the euler equations by finite volume methods using Runge-Kutta time-stepping schemes. AIAA Paper (8ll259); 1981. [38] Rezgui A, Cinnella P, Lerat A. Third-order finite volume schemes for euler computations on curvilinear meshes. Comput Fluids 2001;30:875–901. [39] Xiu D, Karniadakis G. Modeling uncertainty in flow simulations via generalized polynomial chaos. J Comput Phys 2003;187:137–67. [40] Leveque R. Finite volume methods for hyperbolic problems. Cambridge University Press; 2002 [p. 350–3]. [41] Deb K. Multi-objective optimization using evolutionary algorithms. Wiley; 2001. [42] Srinivas N, Deb K. Multi-objective optimization using nondominated sorting genetic algorithms. Evol Comput 1994;2:221–48.