A Unified Geometry Parametrization Method for Turbomachinery Blades

Turbomachinery design is increasingly carried out by means of automated workflows based on high-fidelity physical models and optimization algorithms. The parametrization of the blade geometry is an essential aspect of such workflows because it defines the design space in which an optimal solution can be found. Currently, parametrization methods used for this purpose are often tailored to one particular type of turbomachinery blade, do not provide shape derivatives required for gradient-based optimization, or are not suited to re-parametrize a baseline blade geometry defined by a set of scattered point coordinates in a systematic way. This paper thus presents a general blade parametrization method for axial, radial, and mixed flow blades based on typical turbomachinery design variables and NURBS curves and surfaces. The shape derivatives are computed by means of the complex-step method, allowing the integration of the parametrization into gradient-based shape optimization workflows. In addition, the method enables the re-parametrization of a blade geometry defined by a cloud of points by solving a two-step optimization problem. The capabilities of the method are demonstrated by replicating eight blade geometries in two and three dimensions with an accuracy comparable to the


Introduction
Driven by the ever-increasing requirements in performance, environmental impact, and life-time cost, turbomachinery design is increasingly carried out by means of automated workflows [1].These workflows integrate geometry parametrization tools, highfidelity physical models, and optimization algorithms to systematically explore the design space.The parametrization of the geometry is an essential aspect of the design chain because it defines the design space within which the optimization algorithm can find the optimal solution [2].Ideally, a parametrization method for turbomachinery blades should: 1. Support any type of blade configuration and contain the shapes that achieve the required design objectives.2. Allow the designer to impose geometric constraints due to mechanical or manufacturing requirements.3. Provide the sensitivity of the shape with respect to the design variables to enable gradient-based shape optimization [3].
4. Use conventional engineering parameters with an intuitive geometrical meaning. 5. Produce smooth geometries with continuous curvature (G 2  continuity) and continuous rate of change of curvature to avoid velocity spikes that may lead to flow separation [4].6. Retain compatibility with Computer-Aided Design (CAD) software for further analysis, geometry manipulation, and manufacturing.7. Be computationally cheap in terms of execution time and memory usage.
Shape parametrization methods can be classified into deformation and constructive methods.Deformation methods can be used to modify an existing geometry (a mesh or a CAD model) and are widely used in the context of turbomachinery shape optimization.These methods include mesh point displacement [5,6], CAD model control point displacement [7,8], superposition of shape functions such as Hicks-Henne bumps [9,10], and space morphing methods based on Free-Form Deformation (FFD) [11,12] or on Radial Basis Function (RBF) interpolation [13,14].Although these methods enable the exploration of rich design spaces, they are not suited for an effective handling of geometric constraints, making it difficult to obtain feasible shapes out of the optimization process.As a notable exception, the NSPCC method [8]  allows the designer to impose geometric and continuity constraints by evaluating these constraints at a finite number of test-points and using a projected gradient optimization algorithm to maintain feasibility.In contrast, constructive methods can be used to generate the geometry of a new blade from scratch, or possibly using design variable values obtained from preliminary design models such as mean-line [15,16] or throughflow models [17,18].In addition, they allow one to impose geometric constraints such as minimum blade thickness in a natural and non-intrusive way.Due to these strengths, constructive methods are widely used for turbomachinery blade parametrization and a large number of such methods have been developed over the years.Table 1 provides a comprehensive review of constructive blade parametrizations up to the present day.

Nomenclature
The early constructive parametrization methods used circular arcs and polynomials in monomial-basis form (that is, polynomials in the form ∑ n i=0 a i x i ) to define the geometry of the blades [19][20][21][22][23][24].This type of parametrization gained significant popularity among industry practitioners, but it has severe limitations arising from the use of a monomial basis.Specifically: (1) the polynomial coefficients convey little insight about the shape of the blade, (2) ensuring geometric continuity at the connecting points between segments requires the solution of a linear system that may not have a unique solution, (3) the surface of the blade is prone to undesirable inflection points, and (4) the resulting shapes are not compatible with the geometric representation used by modern CAD systems.
To overcome these shortcomings, several authors proposed new constructive parametrizations based on Bézier [25][26][27][28][29][30][31], Bspline [32][33][34][35][36][37], and NURBS [38][39][40][41][42] curves and surfaces.These mathematical functions have become the standard to represent geometric objects in modern CAD packages due to their favorable mathematical properties and the availability of a wide range of algorithms to define and manipulate curves and surfaces [43,44].Currently, most of the constructive CAD-based parametrizations for turbomachinery blades described in the open literature are not suitable for automated design workflows.This is because they do not offer a robust way to handle trimming and intersection operations [45,46] or do not provide sensitivity information required by gradient-based optimization algorithms [3].In addition, to optimize an existing blade, it is essential to find a parametric representation of the baseline geometry, available, for instance, in the form of a large set of points in the Cartesian space.Solving this reverse engineering problem by trial and error is doable for simple cases [36,42], but it becomes impractical for complex blade geometries.Despite the practical relevance of this problem, a robust and automatic method to re-parametrize the geometry of a blade defined by a scattered set of points is still lacking.
In response to the limitations of the existing methods, this paper presents a general constructive parametrization method for axial, radial, and mixed-flow turbomachinery blades.The method exploits conventional engineering design variables (leading/trailing edge radius, metal angles, blade thickness, etc.) and NURBS curves and surfaces to represent the blade geometry.The method is formulated in an explicit way that avoids the use of intersection and trimming operations to define the geometry of the blade and flow domain and produces blades satisfying G 2 continuity by construction.The sensitivity of the geometry with respect to the design variables is computed with machine accuracy by means of the complex-step method [47][48][49].In addition, the method is also adapted to re-parametrize the geometry of an existing blade defined by a scattered set of point coordinates.This problem, often referred to as blade matching, is formulated as a two-step optimization problem and it allows one to find the set of design variable values that best approximates the prescribed geometry in a systematic way.The flexibility and accuracy of the proposed method is demonstrated by replicating the geometry of eight turbomachinery blades in two and three dimensions.The rest of the paper is organized as follows.Section 2 documents the definition and properties of NURBS curves and surfaces.The blade parametrization in two and three dimensions is described in Sections 3 and 4, respectively, and the computation of the geometry sensitivity using the complex-step method is introduced and verified in Section 5.The blade matching method is presented and applied to replicate a wide range of blade geometries in Section 6.Finally, the software implementation of the method is described in Section 7 and the conclusions are summarized in Section 8.

Background on NURBS curves and surfaces
The origin of Non-Uniform Rational Basis Spline (NURBS) curves and surfaces can be traced back to the research efforts in computer-aided geometric design in the late 60s and early 70s [50].Since then, NURBS curves and surfaces have been universally used for geometrical modeling thanks to their intuitive geometrical interpretation, favorable mathematical properties, and efficient computational algorithms.A NURBS curve, see Fig. 1(a), is a parametric curve defined by where p is the degree of the curve, the coefficients P i and w i are the coordinates and weights of the n + 1 control points, and N i,p are B-spline basis functions defined on the non-decreasing, clamped knot vector with r = n + p + 1.The B-spline basis functions are given by the recursive relation Similarly, a NURBS surface, see Fig. 1(b), is a parametric surface defined by where p and q are the degrees of the surface in the u-and v-directions, the coefficients P i,j and w i,j are bidirectional nets containing the coordinates and weights of the (n + 1) × (m + 1) control points, and N i,p (u)N j,q (v) are the product of univariate Bspline basis functions defined on the non-decreasing, clamped knot vectors with r = n + p + 1 and s = m + q + 1.The u-direction basis functions N i,p (u) are given by Eqs. ( 3) and ( 4), whereas the vdirection basis functions N i,q (v) are defined in an analogous way replacing the variable u by v and the indices i and p by j and q, respectively.NURBS curves and surfaces have the following mathematical properties that make them particularly suited for geometric modeling [43, pp.117-139]: Fig. 1. Construction of a NURBS curve (left) and surface (right).Note that the NURBS curve interpolates its endpoints and it is tangent to the control polygon at its ends.The control net of the NURBS surface interpolates its four corner points and it was represented at an offset distance in the x-direction for clarity.

Table 2
Two-dimensional design variables.Each design variable is provided as a scalar value, except for the upper and lower thickness that are given as sets of control points.

Variable name Symbol
Spacing s Leading edge abscissa and ordinate • Convex hull.NURBS curves and surfaces are within the convex hull of their control points.When the control points are contained in a certain region of space, this property guarantees that the curve or surface will not blow up arbitrarily far away from this region.
• Endpoint interpolation NURBS curves and surfaces coincide with the polytope formed by the control points at the endpoints.
• Endpoint tangency.NURBS curves and surfaces are tangent to the polytope formed by the control points at the endpoints.
• Generalization.Bézier curves and surfaces are a special case of NURBS when p = n and q = m.In addition, B-spline curves and surfaces are an special case of NURBS when all the weights have the same value.
Most of the curves and surfaces used in the proposed blade parametrization method are B-splines.However, the parametrization is formulated in a general way using NURBS so that the user can include the control point weights as design variables to gain more control over the resulting geometry.

Blade parametrization in two dimensions
The proposed two-dimensional blade parametrization is based on typical blade design variables which are listed in Table 2.The geometry of the blade is generated by defining a camber line and subsequently imposing on it two independent thickness distributions in a way that ensures G 2 continuity at the junction between the upper and the lower sides.
The camber line C c (u) is a cubic B-spline curve defined by four control points as shown in Fig. 2(a).The coordinates of the control points are given by where ξ is the stagger angle, c ax = c cos(ξ ) is the axial chord length, θ in and θ out are the inlet and outlet metal angles, and d in and d out are the inlet and outlet tangent proportions.This construction of the camber line ensures that the blade has the specified axial chord length and that the slope at the leading and trailing edges agrees with the input metal angles thanks to the endpoint tangency property of B-spline curves [43, p. 97].
The upper and lower sides of the blade, C l (u) and C u (u), are defined as B-spline curves of degree four as it is the lowest degree that guarantees continuous rate of change of curvature at the spline knots.The coordinates of the control points {P l i } and {P u i }, see Fig. 2(b), are computed according to and The sampling values ûi are given by ûi = The upper and lower thickness distributions, t u (u) and t l (u), are given by B-spline polynomials of degree three with an arbitrary number of control points, {t u i } and {t l j }, specified by the user, see Fig. 2(c).The unitary vectors normal to the camber line n(u) are computed from the unitary tangent vector τ(u) according to where Ċc (u) is computed using analytical derivative formulas for B-spline curves [43, pp. 91-100].The functions f (r) and g(r) appearing in Eqs. ( 12) and ( 13) are used to impose the radii of curvature r in and r out at the leading and trailing edges, ensuring that the parametrization satisfies G 2 continuity by construction, see Fig. 2(d).This feature is important for the aerodynamic design of turbomachinery blades because a sudden change in curvature could cause a spike in the surface pressure distribution or even a local separation bubble [4].The functions f (r) and g(r) are based on the end point curvature formulas for NURBS curves and their derivation is detailed in the Appendix.Once that the upper and lower sides are defined, they can be combined into a single When performing the assessment of the fluid-dynamic performance of the blades via computational fluid dynamics, it is necessary to define the geometry of the flow domain around the blade.For the majority of turbomachinery flow problems one can resort to the periodicity of the flow to reduce the size of the computational domain.In this case, it is therefore sufficient to describe the flow domain around a single blade, which is characterized by the inflow, outflow, and periodic boundaries, as illustrated in Fig. 3.The periodic boundaries are given by two cubic B-spline curves defined by extending the camber line while keeping zero slope at the inlet and outlet.The periodic boundaries are located at an offset distance of half of the blade spacing, s, with respect to the blade camber line.Finally, the inflow and outflow boundaries are defined as two straight lines connecting the upper and lower periodic boundaries.
The proposed parametrization produces blade profiles that have continuous curvature and rate of change of curvature, therefore reducing the possibility of flow separation [4].This contrasts with most of the two-dimensional methods available in the open literature, which produce blades with discontinuous curvature [19,22,28,29,33], or rate of change of curvature [21,23,24,30].As a notable exception, the second and third methods proposed by Korakianitis [23], see also [51,52], produce blade profiles with continuous curvature and slope-of-curvature.However, the methods proposed by Korakianitis involve the solution of systems of equations, are not compatible with CAD representations, and are not easily extended from two to three dimensions.In addition, to the knowledge of the authors, it is the first time that the endpoint curvature formulas for NURBS curves are used to impose the curvature of turbomachinery blades at the leading and trailing edges.This is different than what is documented in previous publications [34,37,42], where all the reported methods used the endpoint curvature formulas for Bézier curves to ensure G 2 continuity, with the limitation that the curvature is not imposed exactly when the blades are described by B-spline or NURBS curves.

Application
The flexibility of the proposed two-dimensional blade parametrization method is demonstrated by reconstructing the four blade profiles illustrated in Fig. 4.Each blade profile was defined using 6 control points for each thickness distribution, resulting in a total of 22 design variables.The LS89 [53,54] and T106A [55] are representative of high-pressure and low-pressure axial gas turbine blade rows, respectively.In addition, the SIRT profile is typical of a supersonic impulse turbine rotor [56] and the STD10 profile is representative of an axial compressor blade derived from a NACA 0006 airfoil profile [57].It can be observed that the parametrization method produces blades with smooth curvature variations, which is essential to avoid spikes and dips in the surface-pressure distribution.The numerical values of the design variables used to produce the blade profiles were computed from a set of scattered point coordinates using the method described in Section 6.From top to bottom: LS89 [53,54], T106A [55], SIRT [56], and STD10 [57].The abscissa of the curvature distribution is the normalized axial length.

Table 3
Three-dimensional design variables.Each design variable is provided as a set of control points that defines a continuous variation, except for the number of blades that is a single integer value.

Variable name Symbol
Number

Blade parametrization in three dimensions
The proposed three-dimensional parametrization is formulated as an extension of the two-dimensional parametrization and uses the design variables listed in Table 3.Similar to the twodimensional case, the parametrization starts by defining a camber surface and subsequently imposing two independent thickness distributions perpendicular to the camber surface in a way that ensures G 2 continuity.
The camber surface is determined by the shape of the blade in the meridional plane and the spanwise variation of the design variables.The shape of the blade in the meridional plane is described by four curves, namely, leading edge, trailing edge, hub, and shroud, as illustrated in Fig. 5.In contrast with other parametrization methods that are limited to axial turbomachines [20,25,27,31,32], the proposed method is suited to describe any kind of turbomachinery configuration, including axial, radial, and mixed-flow machines.The number of control points required to describe the shape of the blade in the meridional plane depends on the complexity of the geometry.For instance, it is possible to define a purely axial turbine using only four control points, but it may be necessary to use 10-20 control points to describe the shape of a mixed-flow machine such as a centrifugal compressor.
The spanwise variation of the some design variables α(v), see Table 3 footnote, is defined as law of evolution through a B-spline of, at most, degree three with an arbitrary number of control points as illustrated in Fig. 7.The number of control points used for each design variable is specified by the user, and its selection is based on the complexity of the blade geometry.As an example, it is sufficient to use a single constant value to define a prismatic blade, but it might be necessary to use 3-6 control points to describe the geometry of a blade with large twist from the root to the tip.
As illustrated in Fig. 9(a), the camber surface S c (u, v) is defined as a bi-quartic B-spline surface with control points P c i,j = [x c i,j , y c i,j , z c i,j ].The coordinates of the control points are computed using the shape of the blade in the meridional plane and the spanwise evolution of the design variables.More specifically, the (x, z) coordinates of the camber surface control points are computed by transfinite interpolation [58] of the four curves that define the meridional plane, see Fig. 6, and are given by In addition, the y coordinates of the camber surface control points at each spanwise location v are given by a third order B-spline curve y c (u, v) with control points {y c 0 , y c 1 , y c 2 , y c 3 } that are computed according to This formulation ensures that the metal angles at the leading and trailing edges, θ in and θ out , are respected, as illustrated in Fig. 8.The arc length of the blade meridional plane at each span location L(v) is defined as and it is computed using 8th order Gaussian quadrature [59], which provides a satisfactory trade-off between computational speed and accuracy.
The upper and lower sides of the blade, S l (u, v) and S u (u, v), are defined as B-spline surfaces of degree four as it is the lowest degree that guarantees continuous rate of change of curvature at the spline knots.The coordinates of the control points {P l i,j } and {P u i,j }, see Fig. 9(b), are computed according to 6. Construction of the four B-splines that define the shape of the blade in the meridional plane (left) and point evaluation by transfinite interpolation (right).Note that the corner control points of the B-splines are shared.and where the sampling values ( ûi , vj ) are given by ûi = The upper and lower thickness distributions, t u (u, v) and t l (u, v), are given by bi-variate B-spline polynomials of degree three with an arbitrary number of control points {t u i,j } and {t l i,j }, specified by the user.The unitary vectors normal to the camber surface n(u, v) are computed according to where the tangent vectors τ u and τ v are given by The partial derivatives of the camber surface with respect to u and v are computed analytically using B-spline surface derivative  formulas [43, pp. 110-115].In addition, the functions f (r) and g(r) appearing in Eqs. ( 22) and ( 23) are used to impose the radius of curvature at the leading and trailing edges, ensuring that the upper and lower surfaces of the blade are smoothly joined with G 2 continuity.The derivation of the functions f (r) and g(r) is detailed in the Appendix.Once that the upper and lower sides are defined, they can combined into a single B-spline surface S b (u, v) = S l ∪ S u that represents the entire blade.
The parametrization just described is suitable to model linear cascades, which are commonly used for wind tunnel tests.However, in actual turbomachines, the blades are arranged in an axisymmetric way forming an annular cascade.In order to achieve this, the linear blade configuration is transformed into an annular one with the mapping H : R 3 → R 3 given by The rationale behind this transformation is to associate the Cartesian coordinates (x, y, z) of a linear cascade with the cylindrical coordinates (x, rθ , r) of an annular cascade and then convert from cylindrical to Cartesian coordinates.
The flow domain around a blade is characterized by the hub, shroud, inlet, outlet, and periodic boundaries as illustrated in Fig. 10.The hub boundary consists of two surfaces that conform with the blade at its root, see Fig. 10(a).Each of these surfaces is defined as a Coons patch [43, pp. 456-507] that is characterized by four edges.The blade edge is given by a B-spline curve formed by extending the lower side of the blade into the upstream and downstream directions following the slope of the camber line at the leading and trailing edges, respectively.The periodic edge is formed by extending the camber line in a similar way and rotating the resulting B-spline about the x-axis through an angle θ b /2, where θ b = 2π /N b .Finally, the inlet and outlet edges are defined as NURBS circular arcs that connect the periodic edge with the blade edge.The shroud surface, see Fig. 10(b), is defined in an analogous way and, for the case of rotor blades, it is possible specify a clearance between the tip of the blade and the shroud.Once that the hub and shroud surfaces are defined, it is straightforward to construct the inlet, outlet and periodic surfaces as ruled surfaces [43, pp. 337-340] that connect the limits of the hub and shroud surfaces as illustrated in Figs.10(c) and 10(d).Note that the parametrization of the blade and flow domain is watertight by construction and it does not rely on intersection and trimming operations.This contrasts with other blade parametrization methods [34,35,[37][38][39][40][41] that rely on intersection operations between the blade surface and the hub/shroud surfaces and produce trimmed NURBS patches that have to be specially treated during a shape optimization workflow [45,46].

Application
The flexibility of the proposed three-dimensional blade parametrization method is demonstrated by reconstructing the four blade geometries shown in Fig. 11.The first example, the AACHEN case, is a prismatic axial turbine stator blade [60].The meridional plane is defined by 4 control points and the design variables are constant in the spanwise direction (1 control point), resulting in a total of 26 design variables.The second case, NASA R67, is an axial compressor rotor blade [61,62].The blade is highly twisted due to the change in blade speed from root to tip and it was necessary to use 4 control points to describe the spanwise variation of the design variables, resulting in a total of 111 design variables.Similarly the XPROP case illustrates the geometry of an aircraft propeller blade [63].In this case it was necessary to use 5 control points to describe the twist of the blades, resulting in 113 design variables.Finally, the APU blade is the rotor of a mixed-flow turbine (radial-inflow, axial-outflow) [64,65].The complex shape of the blade in the meridional plane was described using 14 control points and the spanwise variation of the blade sections was described using 3 control points per design variable, giving rise to 86 design variables.The numerical values of the design variables used to produce the blades were computed from a set of scattered point coordinates using the method described in Section 6.

Sensitivity computation and verification
One simple way to approximate partial derivatives of a function is by using a finite difference approximation such as forward finite differences given as or central finite differences given as where F (α) can be identified with C b (u, α) in two dimensions or S b (u, v, α) in three dimensions and h is the step size used for finite difference computation.Finite difference approximations are susceptible to cancellation error when the step size is small because of the subtraction of very similar numbers in the numerator [66, pp 229-232].As a result, one is faced with the dilemma of selecting a small step size that minimizes the truncation error but does not lead to a large cancellation error.An alternative method that avoids the occurrence of cancellation error is the complex-step method [47][48][49].This method uses a complex argument to compute the first derivative of a realvalued function.Indeed, the Taylor series expansion of F (α) in the imaginary axis gives The sensitivity of the 3D case (right) is represented as a colormap of the scalar field given by the dot product of the sensitivity and the unitary vector normal to the blade surface.
Re-arranging the imaginary part of the equation leads to which is the complex-step method formula.In contrast to finite difference approximations, this method is not susceptible to subtraction error, allowing one to compute first derivatives accurate to the round-off precision of floating point arithmetic by using an arbitrarily small step size.Algorithmic Differentiation (AD) provides yet another alternative to compute the derivatives of a function with machine precision [67].AD is a set of techniques to numerically evaluate the derivatives of a function specified as a computer program by applying the chain rule of differentiation to each arithmetic operation of the program.This method offers more functionality and computational efficiency (first and higher order derivatives, forward and reverse modes) than the complex-step method (first derivatives and forward mode only), but it is also more difficult to implement [49].
In this work, the complex-step method was adopted to compute the sensitivity of the surface coordinates with respect to the design variables due to its accuracy, simplicity, and ease of implementation.Fig. 12 illustrates the sensitivity of the blade surface with respect to one thickness distribution control point in two and three dimensions.It can be observed that the sensitivity of the blade changes from point to point and that there may be regions where the sensitivity is zero as a result of the compact-support property of NURBS basis functions [43, p.118].
To verify the correctness of the sensitivity computation the authors performed a convergence study comparing the sensitivities computed using forward finite differences, central finite differences, and the complex-step method for the NASA R67 test case [61,62].The geometry of the NASA R67 rotor, see Fig. 11, was sampled with 10000 uniformly spaced points within the box (u, v) ∈ [0, 1] × [0, 1] and the sensitivity was computed with respect to one design variable (the stagger angle at the hub) for different step sizes in the interval h ∈ [10 −1 , 10 −15  ].The error of the sensitivity computation was defined as the mean-squareroot deviation between the exact and the estimated sensitivities.
The exact sensitivity was assumed to be that computed with the complex-step method using an step size h = 2.22 • 10 −16 , which corresponds to the machine epsilon of double-precision arithmetic [66, pp 8-11].
The results of the convergence study are shown in Fig. 13.For the case of the complex-step method (CS), reducing the step size decreases the error until the trend flattens to a value close to the machine precision.In contrast, the forward finite difference (F-FD) and central finite difference (C-FD) errors decrease as the step size decreases down to a minimum value and then increase because the cancellation error becomes more prominent than the a Defined as the arithmetic mean deviation between the prescribed and the matched blades.b Defined as the quotient of the mean error and the arc length of blade camber line (three-dimensional cases use the camber line at the hub).
Fig. 13.Sensitivity error at different step sizes for the complex-step method (CS), forward finite differences (F-FD), and central finite differences (C-FD).
truncation error.In addition, it can be observed that the complexstep method and the central finite differences agree in interval when the truncation error dominates (h ≲ 10 −6 ).This verifies that the implementation of the complex-step method is correct.
Although not shown here for brevity, the authors performed similar convergence studies for all the design variables of each of the test cases summarized in Table 4 and obtained similar results.

Blade matching methodology
In order to optimize the performance of an existing turbomachinery blade, it is essential to find a parametric representation of its geometry, which is usually available as a set of scattered points coordinates Q i , with i = 1, 2, . . ., N Q obtained from a mesh, from sampling a CAD model, or from laser scan measurements.This section proposes a systematic method to find the set of design variables that best approximates the shape of a prescribed blade geometry.The method can be divided in two phases: (1) the point projection phase and (2) the geometry update phase.It is assumed that the designer starts from an initial guess for the design variables that roughly approximates the existing geometry, see Fig. 14(a).
In the point projection phase, the goal is to find the parametric values u i , in two dimensions, or (u i , v i ) in three dimensions, that minimize the distance with respect to each prescribed point Q i .The two-dimensional point projection problem can be formulated as where J is the distance between the prescribed and the parametrized point.The gradient of the objective function J can be computed analytically as Similarly, in three dimensions, the point projection problem is given by minimize subject to 0 ≤ u, v ≤ 1 and the gradient of the objective function can be computed according to Note that the geometry of the parametrized blade does not change during the point projection phase.One common pitfall when solving the point projection problem is that the optimization may converge to a local minimum as illustrated in Fig. 14(b).This limitation can be addressed by solving the point projection problem from different starting points and then selecting the global minimum among the various solutions or by sampling the parametrized blade at several locations and then starting the optimization from the test point that is closest to Q i [43, pp. 229-234].
In the geometry update phase, the goal is to find the set of design variables α that minimizes the deviation between the parametrized and the prescribed blades.This can be formulated as an unconstrained minimization problem where the objective function is the sum of the distances between each projected point and the corresponding prescribed point.This optimization problem is given by minimize in two dimensions and by minimize in three dimensions, where α are the design variables listed in Tables 2 and 3, respectively.The gradient of these objective functions is computed using the complex-step method as described in Section 5.In contrast to the point projection phase, the geometry of the parametrized blade is updated until the deviation with respect to the prescribed geometry is minimized, see Fig. 14(c).In order to improve the matching of (u, v) and α, the point projection and geometry update problems are solved alternatively until the relative deviation between the prescribed and the parametrized blades does not change more than a small tolerance, e.g. 10 −8 .To demonstrate its flexibility and accuracy, the blade matching method was applied to replicate the geometry of eight exemplary blades.The set of test cases is summarized in Table 4 and it   was conceived to cover a wide range of turbomachinery blade geometries in two (Fig. 4) and three dimensions (Fig. 11).The results of the blade matching in terms of absolute and relative error are summarized in Table 4.
The relative matching error is below 0.38% for all cases and the highest absolute error is 0.127 mm for the XPROP test case, which is of the same order of magnitude as the tolerances used to manufacture axial gas turbine blades (≈ 0.05 mm) [68].In addition, the deviation could be further reduced by increasing the number of control points used to parametrize the blade.This is illustrated in the convergence study shown in Fig. 15, where the number of control points used to describe the thickness distribution of the T106A case is increased from 3 to 10 points.Specifically, the mean deviation is reduced from 0.057 mm to 0.031 mm when the number of control points is increased from 6 to 10.In addition, Fig. 16 shows the curvature distribution for the T106A blade described using 6 and 10 thickness distribution control points.It can be observed that the curvature variation is smooth for both cases and that increasing the number of control points does not introduce high-frequency undulations that would deteriorate the fluid dynamic performance of the blade.These results indicate that the parametrization and matching methodologies proposed in this work enable the replication of a wide range of geometries with an accuracy comparable to the typical tolerances of modern manufacturing techniques and that the reparametrization accuracy can be increased by refining the design space.

Software structure
The blade parametrization method proposed in this work was implemented in the Python programming language [69] and released under a permissive open source license as the Para-Blade software package [70].The code was written using objectoriented programming principles and the structure of the package is subdivided in the classes shown in Fig. 17.The implementation integrates the in-house NurbsPy package [71] to define and manipulate NURBS curves and surfaces.In addition, ParaBlade relies on the NumPy library [72] for array operations and on the Pagmo/Pygmo optimization library [73] to solve the blade matching optimization problems by means of the limited-memory BFGS algorithm [74,75].

Conclusions
This paper presented a general constructive parametrization method for turbomachinery blades.The method uses typical turbomachinery design variables and NURBS curves and surfaces to produce blade geometries with continuous curvature and rate of change of curvature.In contrast with existing methods, the flow domain parametrization was formulated in an explicit way that avoids intersection and trimming operations and the sensitivity of the geometry is computed by means of the complex-step method, allowing the integration of the parametrization into automated, gradient-based shape optimization workflows.
In addition, the method enables the re-parametrization of a baseline blade geometry defined by a set of scattered point coordinates in a systematic way by solving a two-step optimization problem.To demonstrate its capabilities, the re-parametrization method was applied to replicate the geometry of eight exemplary blades, showing that the proposed parametrization can replicate the geometry of a wide range of turbomachinery blades with an accuracy comparable to the tolerances of current manufacturing techniques for axial gas turbine profiles.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

ξ
AbbreviationsADAlgorithmic Differentiation B-spline Basis Spline BFGS Broyden-Fletcher-Goldfarb-Shanno C-FD Central Finite Differences CAD Computed-Aided Design CS Complex-Step F-FD Forward Finite Differences FD Finite Differences FFD Free Form Deformation NSPCC NURBS-based Parametrization with Complex Constraints NURBS Non Uniform Rational Basis Spline RBF Radial Basis Functions Subscripts and superscripts b Blade c Camber in, out Inlet and outlet l, u Lower and upper sides m Meridional u, v u and v parametric directions 1, 2, 3, 4 Meridional channel edges

x
in , y in Axial chord length c ax Stagger angle ξ Inlet and exit metal angles θ in , θ out Inlet and exit tangent proportions d in , d out Inlet and exit radii of curvature r in , r out Upper and lower thickness distributions t u , t l • Affine invariance.It is possible to apply affine transformations such as rotations, displacements, and scalings to NURBS curves and surfaces by applying the transformation to their control points.

Fig. 2 .
Fig.2.Construction of the blade geometry in two dimensions.The upper and lower thickness distributions (bottom-left) are imposed in the direction normal to the camber line (top-left) to compute the location of the blade control points (top-right).The second and second-to-last control points are computed in a special way to impose the radii of curvature at the leading and trailing edges and to ensure that the blade profile is G 2 continuous (bottom-right).

Fig. 3 .
Fig. 3. Construction of the blade flow domain in two dimensions.The fluiddynamic performance of the blade cascade (left) can be evaluated analyzing the flow around a single blade.The flow domain is defined by four boundaries: inflow, outflow, lower periodic, and upper periodic (right).

Fig. 5 .
Fig. 5. Geometry of the blade in the meridional plane.

Fig. 10 .
Fig. 10.Construction of the blade flow domain in three dimensions.

Fig. 12 .
Fig. 12. Sensitivity of the blade geometry with respect to a thickness distribution control point.The sensitivity of the 2D case (left) is represented as a quiver plot.The sensitivity of the 3D case (right) is represented as a colormap of the scalar field given by the dot product of the sensitivity and the unitary vector normal to the blade surface.

Fig. 14 .
Fig. 14.Illustration of the blade matching problem in two dimensions.The deviation between the prescribed and parametrized blades after the geometry update phase was exaggerated to improve visibility.

Fig. 15 .
Fig. 15.Matching error as a function of the number of thickness distribution control points for the T106A test case.

Fig. 16 .
Fig.16.Curvature distribution of the T106A profile when the matching is performed using 6 and 10 control points.

Fig. 17 .
Fig. 17.Class diagram of the ParaBlade and NurbsPy packages.Each class is represented by a box with three compartments containing its name and main attributes/methods.The diamond-ended arrows (♦) represent a composition relation between two classes.This means that the class next to the diamond symbol contains one or more instances of the class at the other end of the line.For instance, the BladeMatch3D contains one instance of the BladeGeom3D class, which in turn contains several NurbsCurve and NurbsSurface objects.

Table 1
Review of constructive blade parametrization methods documented in the open literature.
cThe parametrization satisfies curvature continuity everywhere except at the trailing edge.

Table 4
Summary of the test cases and matching results.