key: cord-0045658-x530hkqg authors: Raul, Vishal; Leifsson, Leifur; Koziel, Slawomir title: Aerodynamic Shape Optimization for Delaying Dynamic Stall of Airfoils by Regression Kriging date: 2020-05-25 journal: Computational Science - ICCS 2020 DOI: 10.1007/978-3-030-50426-7_5 sha: 56c5ad79b78bf204f41391880e9cd6cc747ae394 doc_id: 45658 cord_uid: x530hkqg The phenomenon of dynamic stall produce adverse aerodynamic loading which can adversely affect the structural strength and life of aerodynamic systems. Aerodynamic shape optimization (ASO) provides an effective approach for delaying and mitigating dynamic stall characteristics without the addition of auxiliary system. ASO, however, requires multiple evaluations time-consuming computational fluid dynamics models. Metamodel-based optimization (MBO) provides an efficient approach to alleviate the computational burden. In this study, the MBO approach is utilized for the mitigation of dynamic stall characteristics while delaying dynamic stall angle of the flow past wind turbine airfoils. The regression Kriging metamodeling technique is used to approximate the objective and constrained functions. The airfoil shape design variables are described with six PARSEC parameters. A total of 60 initial samples are used to construct the metamodel, which is further refined with 20 infill points using expected improvement. The metamodel is validated with the normalized root mean square error based on 20 test data samples. The refined metamodel is used to search for the optimal design using a multi-start gradient-based method. The results show that an optimal design with a [Formula: see text] delay in dynamic stall angle as well a reduction in the severity of pitching moment coefficients can be obtained. The dynamic stall phenomenon was first observed on retreating blades of helicopter rotor [6] . Horizontal and vertical axis wind turbines are prone to dynamic stall. Wind turbines are subjected to dynamic loading from multiple sources, such as wind shear, turbulence, yaw angles, upwind turbine wake, and tower shadow, that cause unsteady inflow to the turbine rotor which results in dynamic stall. In vertical axis wind turbines (VAWT), dynamic stall arises from rapid changes in angle of attack on each blade in every rotation cycle [2, 25] . The dynamic loading in wind turbines generates adverse loading conditions, significantly impacting the blade, hub, tower structure, performance and turbine life. Significant research has been conducted to mitigate or control dynamic stall via active and passive control systems [10, 15, 27, 28] . The addition of structures and control systems to the wind turbines increases their mass as well as their cost and complexity. Mitigating the adverse dynamic stall characteristics passively through aerodynamic shape optimization (ASO) has recently received interest from multiple researchers offering promising improvement in airfoil performance [12, 16, 23, 24, 26] . ASO studies for dynamic stall mitigation are typically done with adjoint-based computational fluid dynamics (CFD) simulations [4, 12, 16, 26] and have shown promising results for multiple dynamic stall optimization cases. Adjoint-based CFD simulations is a modern approach to solve ASO problems using gradient-based optimization (GBO) algorithms [8] . The advantage of the adjoint method is the ability to estimate gradient information cheaply. The GBO approach, however, can get easily get stuck in local minima, especially if the CFD data is noisy. Wang et al. [23, 24] used sequential quadratic programming (SQP) to alleviate aerodynamic loads during dynamic stall cycle on rotor airfoils. Genetic algorithms have the ability to search the design space globally, but they require multiple design evaluations and can be impractical to use for high dimensional design problems. Ma et al. [11] used a multi-island genetic algorithm, which is a global search method, for VAWT performance improvement. Metamodel-based optimization (MBO) (also called surrogate-based optimization) [22] is an approach to alleviate the computational burden of costly simulation-based design problems. In MBO, a metamodel (also called a surrogate) of the objective function is constructed using a limited number of the time-consuming simulations. The surrogate model is fast to evaluate and can be used within GBO or with genetica algorithms to search for the optimal design. To the best of our knowledge, MBO has not yet been utilized for ASO to mitigate dynamic stall characteristics of airfoils. In this work, MBO is used for ASO of wind turbine airfoils to delay stall. The surrogate is constructed using regression Kriging [7] and is sequentially refinement using expected improvement infill criteria. The PARSEC airfoil parameterization technique [20] with six design variables is used for generating the airfoil shapes. The surrogate model is searched using a multi-start gradient-based optimizer. The next section presents the problem statement for dynamic stall mitigation and the setup of the computational model. The following section describes the MBO approach. Results of numerical experiments are presented for the ASO. Conclusions and suggestions of future work are then described. This section describes the problem formulation and the airfoil parameterization method used for the current study, as well as the CFD modeling and validation. The dynamic stall phenomenon is generally studied with sinusoidal oscillating airfoil in a uniform free-stream flow. The pitching motion of the airfoil is described using the angle of attack as a function of time t given as where α m , A and ω represent the mean angle of attack, amplitude of oscillation, and rotational rate, respectively. The reduced frequency, k, is another important parameter and is defined as where c is the airfoil chord length, and U is the free-stream speed. In this work, a deep dynamic stall case from Lee et al. [9] is used. The parameters defining the case are: α m = 10 • , A = 15 • , k = 0.05, and a Reynolds number of Re = 135, 000. The objective of the study is to produce an optimum airfoil shape which mitigates the dynamic stall adverse loading by delaying the dynamic stall angle. This objective is achieved by delaying the formation of the dynamic stall vortex responsible for sudden divergence in the drag and pitching moment coefficients. The optimization problem is formulated as: Here, x is the design variable vector. x l and x u are the lower and upper bounds of x, respectively. The parameters c di , c mi , α ds represent the time variant drag coefficient, pitching moment coefficient at the i th timestep and dynamic stall angle of the airfoil. The subscript '0' represents the baseline airfoil shape, which is the NACA0012 airfoil. Δα denotes the minimum delay in the dynamic stall angle expected in the optimum design, which is set to Δα = 3 • in this work. N denotes the number of time steps in each pitching cycle. For this study, we will only consider the upstroke part of the pitching cycle, which is predominantly affected by formation of dynamic stall vortex. In this work, the PARSEC [20] parameterization technique is used for describing the airfoil shapes. In PARSEC, there are 12 parameters defining the airfoil shape of unit chord. The parameters affecting only the upper surface of the airfoil are considered in this study. The trailing edge offset and thickness are set to zero, which generates a sharp trailing edge airfoil. For this study, we have selected six parameters (see Table 1 ). The current study is performed with the Stanford University Unstructured (SU 2 ) unsteady compressible Navier-Stokes (URANS) solver [17] . The dynamic stall simulations are performed using dual time stepping strategy, rigid grid motion and Menter's shear stress transport (SST) turbulence model [14] . The convective fluxes calculated using second-order Jameson-Schmidt-Turkel (JST) scheme [17] and time discretization is done by the Euler implicit scheme [17] with maximum Courant-Friedrichs-Lewy (CFL) number selected as 4. The two-level multigrid W-cycle method [17] is also used for convergence acceleration. The Cauchy convergence criteria [1] is applied with Cauchy epsilon as 10 −6 over last 100 iterations. No-slip boundary condition is used on airfoil surface with farfield condition on external boundary with Reynolds number of 135,000 and Mach number of 0.1. The c-grid mesh is set up an with outer boundary at 55c from airfoil is generated using blockmesh utility provided by OpenFoam [3] . The mesh is refined near the airfoil surface with first layer thickness to obtain y + ≤ 0.5 and growth ratio of 1.05, which is necessary to accurately capture the onset of the dynamic stall vortex. Figure 1 show a coarse version of the mesh. The grid and time independence study is done in two steps. Initially, the spatial resolution of the mesh is obtained by grid study. This mesh is then used to conduct time study to attain accurate physical time step. The flow and motion parameters are selected from study done by Lee et al. [9] as mentioned in Sect. 2.1. The grid study is done at Re=135,000, angle of attack α = 4 • and turbulence intensity T I = 0.08%. The details of grid study are shown in Table 2 . Meshes 2, 3 and 4 show minimal change in lift coefficient Δc l ≤ 0.003 with the drag counts variation within 4 counts. Considering the simulation time requirement and accuracy of the results, mesh 2 with 387,000 cells is selected for the study. After selecting the spatial resolution, a time independent study is conducted with multiple time steps of an airfoil in a sinusoidal pitching cycle in order to select the temporal resolution. This is done using the generalized Richardson extrapolation method (REM) [18] with the use of average drag coefficient per oscillation cycle c davg as a lower order value to an estimation parameter. The REM estimate c dEst represents the average drag coefficient per cycle at a zero time step, which is calculated as c dEst = 2, 108 counts. Table 3 summarizes the results. The simulation time and estimated error Est err are then considered to select time step of 0.0015 for all further investigations. This section describes the MBO algorithm and the mathematical details of the metamodeling. In particular, the details of the workflow, sampling plan, regression Kriging, infill criteria, and validation are described. A flowchart of the MBO algorithm is shown in Fig. 2 . The presented algorithm consist of an automated loop which sequentially improves the metamodel accuracy. The optimization algorithm starts with a sampling plan where the design space is sampled for initial samples. The initial samples are then evaluated with the CFD model. The regression Kriging metamodel is then constructed for the objective and constraint functions from the initial samples. The constructed metamodel is validated against a test data set. If the model does not pass termination criteria, then an infill strategy is used to refine the metamodel and the above steps are repeated until the metamodel accuracy satisfies the termination criteria. Finally, an optimum design is found by optimizing the metamodel. The accurate construction of metamodel requires an appropriate sampling plan which captures the trend of objective function throughout design space. In this study, Latin hypercube sampling (LHS) [5, 13] is used to generate initial and test data samples. For this study, an initial sample size is considered as ten times the number of design variables. Kriging [19] is a Gaussian based interpolation method widely used in surrogate based optimization [19] . It mainly takes the training point as the realization of the unknown process and approximates as a combination of global trend function plus a localised departure as where x is any sample is the unknown function, G(x) is a known polynomial function and Z(x) is a normally distributed Gaussian process with a zero mean, variance σ 2 , and non-zero covariance [19] providing localised deviation to global trend function. The training samples (x 1 , x 2 , ..., x ns ) in design domain D are correlated with each other through covariance matrix of function Z(x) given by where R is (n s , n s ) symmetric correlation matrix with R ij = R(x i , x j ) a correlation function between any two sample points x i and x j . In this work, we have used the Gaussian spatial correlation function where θ p denotes pth component of vector θ = [θ 1 θ 2 ... θ P ] T , a vector of unknown hyper-parameters to be tuned. The Kriging predictor is given by [19] y where y is the column vector (n s ,1) containing response at sample points, G is a column vector (n s ,1) and filled with ones when G(x) is considered constant. The vector r T (x) = [R(x, x 1 ), R(x, x 2 ), ..., R(x, x ns )] is the correlation vector between known observed points (x 1 , x 2 , ..., x ns ) and the new sample points x. The vectorβ in (9) can be evaluated aŝ The Kriging model is trained over sample data by tuning hyperparameters θ to maximize concentrated likelihood function [5] given by where estimated variance of Kriging modelσ 2 is computed aŝ The Kriging method assumes that the sampled responses are true and do not contain any errors. Typically, most of the engineering functions does have some inherent errors due to involved evaluation process. The objective function in this study could involve errors from the CFD simulation of separated flow region in the dynamic stall cycle. This would produce error in the Kriging approximation when more points are added in close proximity to each other during the optimization process. This problem can be alleviated by using regression Kriging [7] , which allows the Kriging model to do a regression over the sampled data [7] . This is achieved by the addition of a regularization parameter λ to the diagonal terms of the Kriging correlation matrix R, making it R + λI for regression Kriging method where I is an identity matrix. The regularization parameter λ is evaluated by maximizing likelihood function along with θ hyperparameters. The regression Kriging predictor is now given as [5] and varianceσ r 2 of regression Kriging model is computed bŷ where the subscript r denotes regression. The metamodel constructed with regression Kriging using the initial sample data is an approximation of true objective function. The search of the optimal design depends on the accuracy of the metamodel. Although, higher number of initial samples will improve model accuracy it is wise to add infill points strategically in the design space where further improvements in the metamodel are possible. For this study, we will use the expected improvement (EI) infill criteria to provide a balanced exploration and exploitation of the objective function. The EI for regression Kriging is an extension of the EI for Kriging, which uses a re-interpolation technique to make sure resampling of points are avoided. The infill points are obtained by the maximizing EI function, which is written as where y min is the current minimum response, Φ() and φ() are normal cumulative distribution and probability density functions, respectively. The mean square error of the regression Kriging metamodel is given bŷ whereσ 2 ri is the variance of the metamodel with re-interpolation technique given asσ The EI method for regression Kriging with the re-interpolation technique is described in detail by Forrester et al. [7] . In this work, the global accuracy of the metamodel is validated using the normalized root mean squared error (NRMSE) defined as where y i T est andŷ i T est represent responses from the CFD evaluation and metamodel prediction at i th test samples, respectively. The response value y could be an objective function f (x) or constraint function g 1 (x) values for their respective error estimation. The n T indicates the number of test data samples. The denominator of (y max − y min ) I represents maximum and minimum of response values of initial sample I data. In this work, NRMSE ≤ 10% and a fixed budget of 20 infill samples are considered as acceptable criteria for accurate global metamodel. Once an accurate metamodel is obtained it is used by the optimizer to find an optimal design for given problem. For this study, we use a multi-start gradientbased search algorithm to find the optimal design. The sequential least squares programming (SLSQP) algorithm offered by Scipy [21] python package is utilized in this work. A total 240 starting points are used in this study. These start points are distributed over the design space by using the LHS technique. The best obtained result is reported as optimal design. This section presents the results of the metamodel generation and the validation study for the dynamic stall mitigation problem. The optimization results are discussed. As discussed earlier, the optimization algorithm generates the metamodel and sequentially refines it. Initially, the design space is sampled using LHS. A total of 60 design samples (10× number of design variables) are generated. Each design sample is then evaluated with the CFD module to generate the objective and constraint function values. Note that in this study we only simulated the upstroke of the pitching cycle where dynamic stall vortex formation occurs. The obtained observations are used to construct two separate metamodels, one for the objective and another for the constrained function. Both these metamodels are validated with 20 test data points (one third of initial samples). The test data points are also generated using LHS technique separately and evaluated with the CFD module. The global accuracy of the metamodel is tested using the NRMSE metric. If the accuracy of the model satisfies the termination criteria then it is passed to the optimizer, else an infill point is evaluated and added to the initial sampling plan to construct a new metamodel. This process is iterated until the metamodel satisfies the termination criteria of NRMSE ≤ 10% and fixed budget of 20 infill points. Figure 3 shows a plot of the NRMSE for the objective and constraints functions every 5 infill points. It can be seen that both the metamodels satisfy global accuracy error criteria well before infill points reach the fixed budget criteria. The constraint function metamodel shows a higher accuracy than the objective function metamodel reaching 2.4% and 8.8%, respectively, by total 80 sample points (60 initial samples plus 20 infill points). Figure 4 shows the baseline and optimum airfoil results. Table 4 gives the aerodynamic characteristics of the airfoils. There are major shape variations between the baseline (NACA0012) and the optimized airfoil. The optimized airfoil has a higher maximum thickness (t/c max = 0.146) with a maximum camber (M ) = 1.89% located at x/c = 0.62. The optimum design is able to delay the dynamic stall angle (α ds ) by more than 3 • , whereas the moment stall angle α ms is delayed to 20.26 • . The α ms indicates formation of dynamic stall vortex which is responsible for sudden divergence in drag and pitching moment coefficients. The delay in dynamic stall vortex formation provides an increase in operational range without adverse loading on the airfoil. Moreover, optimum shape also shows the reduction in severity of pitching moment (Fig. 4d) . Figure 5 shows z-vorticity contour plots of baseline and optimum airfoil near moment stall and dynamic stall angles. It can be seen that near the moment stall and dynamic stall point of baseline airfoil, the optimal shape does not show any signs of dynamics stall vortex formation which verify details given in Table 4 . In this work, efficient aerodynamic shape optimization using regression Kriging metamodeling is used for mitigating the adverse effects of dynamic stall on wind turbine airfoil shapes. The optimal airfoil shape shows a significant delay in the dynamic stall angle when compared to a baseline airfoil. It was found that the optimal shape has a higher maximum thickness and maximum camber compared to the baseline airfoil. Future work will consider global sensitivity analysis to provide the sensitivities of the individual variables with respect to objective and constraint functions, and to explore the interaction effects of variables. This will reveal how the airfoil aerodynamics affects dynamic stall response. Understanding Analysis Dynamic stall in vertical axis wind turbines: comparing experiments and computations Openfoam for computational fluid dynamics Unsteady aerodynamic design on unstructured meshes with sliding interfaces Engineering Design via Surrogate Modelling: A Practical Guide Blade stall half fact, half fiction Design and analysis of "Noisy" computer experiments Comparison of gradient and response surface based optimization frameworks using adjoint method Investigation of flow over an oscillating airfoil Dynamic stall flow control via a trailing-edge flap Airfoil optimization to improve power performance of a high-solidity vertical axis wind turbine at a moderate tip speed ratio Adjoint-based unsteady airfoil design optimization with application to dynamic stall A comparison of three methods for selecting values of input variables in the analysis of output from a computer code Two-equation eddy-viscosity turbulence models for engineering applications Dynamic stall control via adaptive blowing Optimum shape design for unsteady flows with time-accurate continuous and discrete adjoint method Stanford university unstructured (SU 2): an open-source integrated computational environment for multi-physics simulation and design Grid convergence error analysis for mixed-order numerical schemes Metamodels for computerbased engineering design: survey and recommendations Parametric airfoils and wings SciPy 1.0-fundamental algorithms for scientific computing in python Review of metamodeling techniques in support of engineering design optimization Rotor airfoil profile optimization for alleviating dynamic stall characteristics Aerodynamic shape optimization for alleviating dynamic stall characteristics of helicopter rotor airfoil Numerical investigations on dynamic stall of low reynolds number flow around oscillating airfoils Investigation of effect of dynamic stall and its alleviation on helicopter performance and loads Dynamic stall control for advanced rotorcraft application Dynamic stall control optimization of rotor airfoil via variable droop leading-edge Acknowledgements. The second and third authors were supported in part by RAN-NIS grant number 174573.