Multi-Level Multi-Objective Multi-Point Optimization System for Axial Flow Compressor 2D Blade Design

This paper introduces a multi-level framework to perform a multi-objective multi-point aerodynamic optimization of the axial compressor blade. This framework results in a considerable speed-up of the design process by reducing both the design parameters and the computational effort. To reduce the computational effort, optimization procedure is working on two levels of sophistication. Fast but approximate prediction methods has been used to find a near-optimum geometry at the firs-level, which is then further verified and refined by a more accurate but expensive Navier–Stokes solver. Surface curvature optimization was carried out in a first-level as a meta-function. Genetic algorithm and gradient-based optimization were used to optimize the parameters of first-level and second-level, respectively. This procedure considers both the aerodynamic and mechanical constraints. An initial blade has been optimized with three different design targets to highlight the ability of the design method and to develop design know-how. Leading-edge shape and curvature distributions of pressure and suction surface had major effects on the design philosophies of the blades. The result shows about −22.5 % reductions in pressure-loss coefficient at design condition and 23.6 % improvement in the allowable incidence-angle range at off-design conditions compared to the initial blade.

The growing market of turbomachinery products causes a reduction in the times and costs of design by considering alternative approaches, particularly in the aerodynamic design of blades. Many researchers have developed approaches for optimal design of blades aided by numerical algorithm. In this regard, one approach is to automate the conventional design process by coupling an optimization method, CAD-based blade generators, and a computational fluid dynamics (CFD) code. The earliest parametric design method was presented by Dunham [1]. In this method, the thickness distribution of blade was optimized around the camber line. However, the nature of flow phenomena in the pressure and suction sides of the blades promoted the next activity to design blade surfaces independently. Corral and Pastor [2] used the Bezier curve to generate turbomachinery blades' surface. Keshin et al. [3] developed a tool for automatic multi-objective optimization of blade geometry with respect to loss and compressor operating range. Bonaiuti and Zangeneh [4] proposed an optimization strategy, which enables the three-dimensional multi-point multi-objective aerodynamic optimization of turbomachinery blades. The design strategy was based on the coupling of three-dimensional inverse design, response surface method, multi-objective evolutionary algorithms, and computational fluid dynamics analyzes. Korakiantis [5] presented the hierarchical development of three direct blade-design methods of increasing utility for generating two-dimensional blade shape. These methods can be used for generating inputs to the direct or inverse blade-design sequences in subsonic or supersonic airfoils in compressors and turbines or isolated ones. These blade design methods can be applied for improving the aerodynamic and heattransfer performance of turbomachinery cascades, leading to perform airfoils with in very little iteration. Many studies in the turbomachinery field (Sonoda et al. [6], Shahpar and Radford [7], Kammerer et al. [8], Buche et al. [9], Benini and Tourlidakis [10], Bonaiuti and Pediroda [11]) have tried to exploit the time-saving potential of such a strategy and to integrate it into their design systems. Samad and Kim [12] performed a multi-objective optimization of an axial compressor rotor blade through genetic algorithm with total pressure and adiabatic efficiency as objective functions. By the optimization, maximum efficiency and total pressure are increased by 1.76 and 0.41 %, respectively. Ju and Zhang [13] established a multi-point and multi-objective optimization design method, particularly, aiming at widening the operating range while maintaining good performance at the acceptable expense of computational load. They used the artificial neural network and the genetic algorithm technique and back-propagation algorithm along with the computational fluid dynamics. Compared with their original design, optimized cascade decreases the total pressure loss coefficient by 1.54, 23.4, and 7.87 % at the incidences of 5 • , −9 • , and 13 • , respectively.
Pierret and Van den Braembussche [14] presented a multilevel method for the automatic design of more efficient turbine blades. An artificial neural network (ANN) has been used to construct an approximate model using a database containing Navier-Stokes solutions for all previous designs. This approximate model has been used for the optimization of the blade geometry by means of simulated annealing, which is then analyzed by a Navier-Stokes solver. This procedure resulted in a considerable speed-up of the design process by reducing the interventions of the operator and the computational effort.
The parametric description of the blade is related to the difficulties, which are conventionally performed by means of purely geometrical parameters [15]. A possible solution to this problem is to use a blade parameterization technique based on the parameters related to the fluid flow. Utilizing an optimized curvature distribution on blade surface prevents the occurrence of most of these drawbacks. This paper presents a framework to perform multi-objective multi-point aerodynamic optimization of axial compressor blade which considers all aerodynamic and mechanical constraints. A multi-level optimization is used to speed-up design time; therefore, the blade geometry parameters were classified into two categories and optimized in two loops with the distanced objective function. A combination of the gradient-based method and genetic algorithm is used to optimize the profiles. The main advantage of this approach is to use optimized surface curvature distribution to build up blade profile which was done in the first-level optimization. The aerodynamic performance of blade section in the design and off-design conditions is evaluated by a two-dimensional Navier-Stokes blade-to-blade flow solver. An initial blade is optimized with three different design philosophies to highlight the ability of design method in developing design know-how.

Problem Description
High aero-performance is not the only objective of an optimization. A good design should perform well at off-design  Figure 1 schematically shows the loss variations versus incidence angles for a sample blade. The value and location of minimum loss and incidence angle range between blade's stall and choke can be observed in this figure. Ideally, design is done with the aim of increasing α and reducing ω min at specific ∝ min . However, boosting the blade loading will demand a decrease in α and an increase in ω min . When incidence variation at the off-design is not high, a blade can be designed with the aim of reducing the minimum loss, and when incidence variations at the off-design is high, a blade should be designed with the aim of increasing α.
Performance of a blade is only one of the many considerations among the design objectives. In fact, the mechanical constraints and the imposed flow angle must also be satisfied. The general approach to this problem is to build an objective function, which is the summation of penalty terms, to limit the violations of the constraints. Equation (1) shows the proposed objective function in this paper for blade design where ω min.ref , α min.ref and α ref are the reference values calculated in the through-flow design step. w 1 , w 2 , and w 3 are the weights of each term at the objective function. These weights were used to dictate design philosophies to the design system [16].
where P T is the penalty term which is used to consider optimization constraints. Computational cost of mentioned objective function is high, because it needs a calculation of cascade flow filed at several operational points (multi-point optimization). Increasing the number of parameters improves blade geometry generation flexibility, but time and cost of blade design would increase. If blade parameters have a closer relationship with the parameters influencing the flow field, the speed of design process may be increased. Boundary layers on the blade surfaces and blade losses are affected by both curvature value and curvature gradient of the blade surfaces. The terms 1/r and 1/r 2 in compressible Navier-Stokes equations in polar coordinate show strong dependence of velocity and pressure on the curvature. Small slope-ofcurvature discontinuities in the blade surfaces result in Mach number or pressure coefficient surface-distribution spikes or humps, which may introduce regions of local acceleration and deceleration, separation, or undesirable loading distributions along the length of the blades [15]. In this paper, velocity diagram and its parameters were predefined from axisymmetric through-flow design. The following non-dimensional parameters were also specified: the tangential spacing (s/c) between blades, opening diameter o/c, the trailing-edge thickness and the blade maximum thickness (the minimum allowable value of blade's maximum thickness was specified from mechanical constraints).

Airfoil-Design System
Typically, the blade-design system developed in this work includes a geometry generation, a cascade flow analysis tool, and an optimization algorithm. Because of high computational cost of the objective function, the optimization and geometry tools are modified to optimize the airfoil with a minimum iteration. However, fast and reliable flow solver is also necessary.

Geometry Tools
As shown in Fig. 2, the cascade geometry is specified from three curves: leading edge, suction surface and pressure surface. Each part of the blade surface has an analytical polynomial, which is evaluated by mapping a described curvature distribution. System of differential equations with boundary value [Eq. (2)] is solved to calculate surface coordinates using the curvature distribution.
where a and b are boundary points and C(x) is the curvature distribution on surfaces.
The suction surface is made of two polynomials: one from p f (leading edge) to p o (indicating the minimum distance between suction surface and stagnation point, evaluated by opening and the ϕ); another, from p o to the trailing edge. Each part of curvature distribution is made of a curve with four points. The first part of suction surface is evaluated by two conditions: locations of p f and p o , and curvature distribution specified by the curvature of these points plus two intermediate points. The second part is evaluated by four conditions: location of p o , slope of suction surface at p o , location of p t (trailing edge point), and slope of suction surface at the trailing edge (β s−out ). The second part of the suction surface is defined by four boundary conditions and by Bezier curve with four points. This is because the curvature value of two points of Bezier is a parameter.
Pressure surface consists of two curves: one from leading edge to maximum thickness location and another from maximum thickness to trailing edge. Pressure surface should be defined to satisfy the maximum thickness location and maximum thickness value. First and second derivative of thickness distribution should be zero and negative, respectively. Therefore, suction surface coordinate is calculated by Eq. (3).
There are three conditions for calculating the thickness distribution of the first part: thickness value at the leading edge, maximum thickness value and thickness slope at x max . Curvature distribution of pressure surface is also determined by the four-point Bezier curve. Equation (4) is solved to calculate the thickness distribution of each part of the suction surface.
There are three and four boundary conditions for the first and second sections, respectively [value and slope of the thickness curve at x max , and the trailing edge where (t = y p − y s )]. Pressure surface curvature distribution in each section is defined by the four-point Bezier curve. Therefore, curvature value of one point at the first section and two points at the second section should be a parameter at each Bezier curve.
A leading-edge bluntness mode is implemented by fitting a profile to the nose of the airfoil section. The profile is defined by Eq. (5).
is a concave-down parabola centered at x = 1/2 with maximum value of 1. As b increases, so does the bluntness of the leading edge. Figure 3 shows the nose of an airfoil to which the five different profiles with b values of 1-5 have been fitted.
Trailing edge is defined by a semicircle. The semicircle diameter is a manufacturing limit which should not be less than the allowed minimum in the optimization process.

Flow Solver
Accurate calculations of the cascade loss in different conditions are necessary to determine an objective function. Thus, a flow analysis tool including mesh-generation tool and computational fluid dynamics solver are used to calculate cascade losses. The flow field is automatically divided into several zones in a way that one structured mesh is used for each zone which is able to make a fine mesh near the walls. A two-dimensional compressible viscous flow code is used to simulate the cascade fluid flow. This code is developed based on Roe scheme, so that it could be used to design transonic blades also. A method has been presented by Kermani and Plett [17] is used to solve the fluid equations by the Roe scheme in the computational domain and to give a formula for the Roe's numerical fluxes in generalized coordinates. The fluid governing equations for the viscous, unsteady and compressible flow in generalized coordinates with no body force could be shown as follows: where Q 1 is the conservative vector, F 1 and G 1 are the inviscid flux vectors, and G 1VT is the viscous flux vector. Because of high-speed flow in cascade, all the viscous derivatives along the mainstream of the flow are neglected (thin-layer approximation). Therefore, the governing equation can be discretized as follows.
For the time discretization, the two-step explicit scheme, from the Lax-Wendroff family of predictor-correctors, is used. The inviscid numerical flux F 1E based on the Roe scheme is written in generalized coordinates, according to Ref. [17]. Other flux term can be written the same as F 1E .
To obtain the left (L) and right (R) flow conditions, a thirdorder upwind-based algorithm with the MUSCL extrapolation strategy [18] is applied to the primitive variables (pressure, velocity components and temperature). For example, at the east cell face of the control volume, E, the L and R flow conditions are determined as follows (k = 1/3 in this study).
For successful modeling of the turbulent effects on the flow, the K − ω (SST) method has been utilized which incorporates modifications for low-Reynolds number effect, compressibility and shear flow spreading. k − ω (SST) model consists of a transformation of the k − e model in the outer region to a k − ω formulation near the surface by a blending function F1. Equations (10,11) are k and ω equations of the model, respectively.  In this model, eddy-viscosity formulation by a limiter is used to consider the transport of the turbulent shear stress.
F 2 is a blending function which restricts the limiter to the wall boundary layer. S is an invariant measure of the strain rate. The blending functions are defined based on the distance to the nearest wall and on the flow variables by Eqs. (13,14). where: All coefficients of model are listed in Table 1.
To achieve accurate and reliable results, the tools could make the mesh coarser or finer considering the input condition and the turbulence model. Details of domain mesh have been presented in Table 2. Boundary layer mesh is used to satisfy y + criteria. Domain mesh has been shown in Fig. 4.
The surface pressure coefficient [Eq. (15)] distributions at design point conditions resulting from the experimental cascade tests [19] and Navier-Stokes solver calculations for a controlled diffusion blade are shown in Fig. 5. A good agreement between experimental data and simulated by presented flow solver was observed. This agreement was pronounced along the whole blade surfaces. Therefore, the comparison indicated that the present flow solver has enough accuracy to be used in this work.

Algorithm and Optimization Tool for Blade Design
High computational cost of objective function and large number of parameters lead to high design computational cost. To reduce optimization time, design variable parameters are divided into two categories. For each category, an optimization level is considered. The first category is related to the curvature and the curvature slope of blade surfaces which directly related to the both Mach and the pressure distribution. The second category is used to define the overall shape of blades (maximum thickness location, thickness and shape of the leading edge). These parameters are optimized using the objective function in Eq. (1). Figure 6 depicts blade-design system components. The input blade has been re-parameterized to create initial guess of optimization process. The parameters are calculated with least deviations between input blade and resulting geometry of these parameters. In addition to coordinates of points, slope and curvature of the surfaces are compared to calculate these deviations.
The whole shape of blade including the leading-edge and trailing-edge thickness to the maximum thickness, the location of maximum thickness, and the leading-edge shape factor is optimized at an outer loop with Eq. (1) as objective function. Table 3 define optimization problem for both optimization levels. As shown in Fig. 6, in each iteration of outer loop, first-level variables need  In the first-level optimization, the surfaces curvature distribution optimization is carried out at an inner loop using the integration of normalized curvature slope as its objective function which summarized in Table 4. Considering that the cost of optimization based on the first-level objective function is low, therefore, optimization of these parameters has little impact on the total time of the blade-design process.
A gradient-based optimization method is used to optimize the second-category parameters on an outer loop. Because of the design is established aiming at minor changes on the overall shape of the blade. Objectivity, computational cost of gradient-based methods is low. A genetic algorithm is also used to optimize the parameters of first category on an interior loop. The genetic algorithm is applied to prevent the optimization process from falling into local optimal points. Since the computational cost of an inner-loop objective function is low, the use of genetic algorithms for an inner loop is justifiable. Finally, the output geometry is determined based on the criteria of objective function in Eq. (1). Interior optimization loop are used to reduce the number of influencing parameters and optimization time.

Results and Discussion
Ability of each level of multi-level optimization system is invested in this section. An initial profile is optimized with nine different objects by first-level optimization process to evaluate effect of curvature distribution. Three blades with different objective function have been optimized with twolevel optimization. In the end, performance of optimized profiles is compared with initial profile and a well-known NACA-65 profile.

Firs-Level Optimization by Curvature Distribution Assessment
The objective of airfoil design is to curve flow to create lift force with minimal stagnation pressure drop which is usually due to: (1) shear forces between fluid and airfoil, and (2) an imbalance input and output pressures; both of which are affected by the side-wall boundary layer. The blade design affects boundary-layer development in three ways: first, the positive pressure gradient on the blade suction surface that causes flow separation; second, drastic changes in surface curvature of the leading edges that make over-speed zones and, consequently, the local separation, and, third, diffusing regions created on the pressure surface. The designer should define the surface curvature avoiding drastic changes on the trailing-edge curvature and optimizing surface curvature distribution which controls the boundary layer at pressure and suction surfaces. Blade curvature distribution can be changed using the weight function, H (x), which has been defined as a linear combination of three sinusoidal-shape function illustrates by Eq. (16). The weighing coefficients, y i , are the outer-loop variables which are determined through multi-level optimization process. Actually, the optimization loops are related to each other by these weighting coefficients.  , x fo1·s , c fo2·s , c fo·s , x ot1·s , c ft1·p , x ft1·p , c ft2·p , x ft2·p , c tt1·p , x tt1·p , c fp , c max·p ) Subject to: Here, xk i values are the locations of maximum height of corresponding shape functions. Figure 7 shows the effect of y i variation on H (x) where y 1 + y 2 + y 3 = 1, y 2 = 0.1 : 0.5, y 3 = 0 : 0.3. For example, if y 1 become larger, location of the maximum height move toward the front of the blade, consequently front curvature slope of blade became larger (front loaded). This blade had relatively long surface lengths, with relatively less adverse pressure gradients toward the pressure and suction surface.
The surfaces curvature optimization is carried out at the first level. Three groups of curvature distribution are investigated. In the first one, weight coefficients of H (x) are of same order (y 1 ∼ y 2 ∼ y 3 ), the second group is mid-loaded (y 3 > y 2 > y 1 ), and third group is also front-loaded (y 1 > y 2 > y 3 ). An initial blade with inlet metal angles of 50 • , outlet metal angle of 30 • and t max /c = 0.2 has been used as initial guess. Design variables which are presented in Table 2 are summarized in the three mentioned groups and the other geometric parameters are kept fixed. Figure 8 shows loss coefficient distribution and curvature distribution comparison between the reference and optimized blades. In the first group, the maximum curvature is located near leading edge which is reduced with almost uniform slope. Allowable incidence rang has been increased and minimum loss has been decreased in this case. In the second group, the maximum curvature is located in the middle of the blade surfaces and only the minimum loss is reduced. In the third one, only allowable incidence angles range has increased. In the best cases, about 16 % decreasing in minimum loss and 20 % increment in allowable incidence angle rang has been achieved. Similar changes on loss curve of optimized blade in each group result in a very close relationship between design philosophy and curvature distributions of surfaces. A blade with an expected performance can be designed by setting the proper weights of H (x) which are design variables of the outer optimization loop.

Multi-Level Optimization of Blade Profile
To investigate design tool abilities, three blades with different design philosophies were optimized. First case is aimed to reduce the value of loss at the design point (w 1 = 0.8, w 2 = 0.2, w 3 = 0.0). In the second case, the goal is to widen the operating range for off-design conditions (w 1  The first-case airfoil Comparison of the airfoils geometry is presented in Fig. 9a, the optimized profile shows almost same stagger angles and thinner leading edge than the initial airfoils. Similar to initial profile, the new airfoils show more camber in the mid-part and less camber in the front and rearpart (Fig. 9b). In as a starting profile, the design losses were reduced by more than 22 %; the incidence range is reduced by more than 3 • . As shown in Fig. 9c, design incidence angle has not changed more. The corresponding design and off-design Mach number distributions are shown in Fig. 9d. At design condition, first case is characterized by a smaller Mach number peak on leading edge than initial profile, flow acceleration at beginning of the suction surface, boundary-layer transition at about 20 % chord, and smaller deceleration gradients in the mid-part and rear-part of the airfoils which result in lower pressure losses than initial profile. At positive incidence angle (α = +5 • ), Mach number peak on the leading edge is smaller than initial profile but boundary-layer separation and re-attachment occur at pressure surface; therefore, pressure losses is not significantly improved. Compared to initial blade, optimized blade pressure loss has been slightly improved at negative off-design incidence angle (α = −15 • ) because of the following flow characteristics: Mach number declaration gradient is less on the suction surface, flow separation on the pressure surface is occurred at x/c = 0.7 than x/c = 0.3, and a smaller Mach number peak at the leading edge.
The second-case airfoil Figure 10a, b compares geometry and surface curvature distributions of the second case and initial profile. In the second surface case, curvature peak is located near the leading edge same as second airfoil group in last section. As shown in Fig. 10c, the considerable improvement in the off-design behavior of the optimized profile was carried out where the incidence range increased and the design losses decreased. Comparison of Mach number distributions is shown in Fig. 10d; the most significant difference is visible in the upstream propagation of the peak suction side Mach number for the optimized airfoil. Suction side diffusion starts shortly after the Mach number peak location, when the turbulent boundary layer is still thin. On the other hand, this leads to a smaller deceleration gradient for the optimized airfoil in the downstream of the peak Mach number. On the one hand, this means that pressure loss at design incidence slightly improved by optimization.
The most significant improvement is observed at offdesign negative incidence angle. Compared to the initial profile, the second-case airfoil is characterized by leading-edge shape treatment (the lower suction side peak Mach number at the optimized airfoil), a flattened mid-part, a front-loaded Mach number distribution, and removing the flow separation at pressure side result in lower total pressure losses and considerably higher operating ranges at negative incidence angle.
The higher local deceleration gradient results in the higher local momentum thickness growth rate. At positive incidence angle (α = +5 • ), between 5 and 80 % of chord (Fig. 10d), the diffusion on the optimized geometry is significantly higher and the boundary-layer momentum thickness grows with a higher gradient; therefore, the total pressure losses for optimized airfoil has been deteriorated at positive incidence angle.
The third-case airfoil An equal weight was used for the first and third terms of the objective function. Geometry and surface curvature distributions of third case and initial one have been compared in Fig. 11a and b, respectively. In the second surface case, curvature peak is located near leading edge as same as the second airfoil group in last section. It was observed that the minimum loss reduced about 6.6 % and the range of incidence angle improved by 9.6 % than the initial profile (Fig. 11c). To explain the reasons for the considerable increase in operating range, the Mach number distributions at +5 • and −15 • incidence for the third cascade are compared in Fig. 11d. The corresponding incidence flow angles are marked in the total pressure-loss diagram in Fig. 11c. The importance of leading-edge treatment is demonstrated by the off-design behavior of optimized profile, while the peak Mach numbers of initial blade is about 0.8, the new airfoil avoids these peaks at −15 • incidence and stays at a lower Mach number level. This peak reduction leads to a significant boundary layer unloading in the vicinity of the leading edge and finally results in an improvement of the separation behavior. The more camber in the front, the almost flattened mid-part and the lower curvature of rear-part of the profile lead to the reduced design point losses.
The NACA-65 profile minimum loss is greater than all other cases because of great much number at leading edge and great curvature slope near leading edge. But off-design behavior of this blade is better than that of initial profile and the first case. To explain the reasons, there is no boundarylayer separation at negative off-design condition (Fig. 11d). In positive incidence angles, separation is occurred near leading edge rather than initial profile and first case, therefore, its loss at α = +5 • is somewhat greater.   Table 5 summarized run time and aerodynamic performance of optimized blade. Optimized profile results are compared with initial one and a well-known NACA-65 profile which has the same inlet and outlet metal angle as initial profile. As expected, Maximum reduction of design point loss coefficient has been occurred in case 1 by −22.5 % and the greatest increase in allowed angles of attack ranges were being entailed to the second case by 23.6 %. Calculations have been done by a computer with Intel core(TM) 2 quads CPU 2.8 GHz. For typical problem, presented multi-level optimization is about 2.5 times faster than single-level optimization via genetic algorithm because of numbers of parameters that are optimized using Navier-Stokes analysis in multi-level procedure (7 parameters) are lesser than that in single-level optimization (17 parameters). The main deficiency of using multi-level optimization is the risk that the discrepancies between the predictions by the meta-function and the Navier-Stokes results drive the optimizer to a false optimum. Presented meta-function for first-level optimization is fast but inaccurate and using it may drive optimization process to a non-optimum combination of design parameters. Any further control by optimization of y i using an accurate Navier-Stokes solver will diminish the inherent inaccuracy of the fast calculation method but there is no mechanism to eliminate it.
The method presented in this paper can consider different design philosophies using different weights in the objective function (multi-objective ability). Weights in objective function can be determined based on the expected application of the blades and the designer considerations. For instance, the blade designed in the first case is suitable for the middle stage of compressors, because the range variation of their angle of attack is not very high at the off-design condition; where the compressor performance can be optimized by reducing the loss at the design condition. The third case is appropriate for the rear and front stages of the compressors where blade incidence variations are high. In optimum conditions, the spectrum of different weights may be considered by the designer for different sections of compressor blades.

Conclusions
In this paper, a multi-level airfoil design tool was developed aiming to decrease loss in design and off-design conditions. To verify the ability of the first-level optimization (metafunction) in the blade design process, an initial blade was optimized with nine different design objectives. In best cases, optimization for off-design condition resulted in an increase of 20 % in allowed incidence angle range. Design incidence angle loss has been reduced by 16 % when optimization for design condition was desired.
Three different blade design philosophies are considered to optimize an initial blade with multi-level optimization; the first one is to reduce blade loss at the design condition; the second one aims to increase the range of incidence angle at the off-design, and the third one is to reduce the minimum loss and increase the incidence angle range. The highest performance improvement at design condition in the first case is −22.5 % reductions in loss coefficient at design incidence angle. In the second case, 23.6 % increase in allowable incidence angle range, improves blade's performance at off-design conditions. It is also found that the design philosophy affects the leading-edge shape, maximum thickness location and maximum curvature location and value. Presented multi-level optimization is about 2.5 times faster than single-level optimization via genetic algorithm with same geometry parameters and objective function. Therefore, the presented system is successful in both time-saving and optimization with various targets.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.