Open-Source Implementation and Validation of a 3D Inverse Design Method for Francis Turbine Runners

The hydraulic design of Francis turbines and pump-turbines is an expensive project-specific engineering effort that typically involves a direct iterative exploration of the design space. An inverse design method for turbomachinery has been previously introduced in the literature, and several recent applications have demonstrated its advantages; however, only a commercial implementation of the method is currently available. In this work, an open-source implementation of the inverse design method is introduced. First, the governing equations in cylindrical and curvilinear coordinate systems are derived, consolidating the somewhat inconsistent formulations that are available in the literature. Then, a convergence analysis of the method is performed in order to characterize the behavior of the discretization error and deduce the mesh resolution requirements. A validation of the method output with respect to high-fidelity computational fluid dynamics simulations is then presented; it is demonstrated that the velocity fields are well predicted, the pressure distribution on the blades is reasonably well approximated, and the flow angular momentum extraction is achieved in the prescribed manner. Possible improvements to the open-source implementation of the method are discussed.


Introduction
Hydropower accounted for 62% of the global renewable electricity production of 2017 [1], and its global importance is bound to increase for two fundamental reasons. On the one hand, a rapid electrification of transportation is already underway, implying an ever-growing demand for inexpensive electricity with low lifecycle greenhouse gas emission intensity. On the other hand, the exponential increase of the installed capacity of intermittent renewable energy sources such as wind and solar urges a similar increase of the energy storage capacity, 96% of which is currently supplied by pumped storage [1].
The hydraulic design of Francis turbines is a project-specific engineering effort that involves an iterative exploration of the design space. In spite of modern optimization techniques implemented in the framework of computational fluid dynamics (CFD), the hydraulic design of Francis turbines still demands considerable engineering input and know-how. Design methodologies currently used in industry are direct: they require the specification of the blade geometry in terms of a distribution of angles at the leading and trailing edges, as well as a distribution of the blade angle variation along the chord. The design is then evaluated by means of CFD and modified iteratively until a satisfactory compromise among the design objectives is achieved, after which expensive experimental verification is undertaken. Improvement of the design methodology may result in greater competitiveness and faster project development.
An inverse design method involves calculating the set of system parameters, e.g., the turbine geometry, that will result in the desired system behavior, e.g., the runner efficiency or flow field characteristics. One of the earliest 3D inverse design methods for turbomachinery was proposed by Tan et al. [2] for highly loaded, infinitely thin blades in an annular cascade configuration, assuming inviscid and incompressible flow as well as a simplified meridional channel geometry. This seminal work was later extended by Borges [3,4] to account for meridional channels of arbitrary geometry, followed by Zangeneh [5], who further generalized the method to describe compressible flows and blades of finite thickness. Soon thereafter, Zangeneh co-founded Advanced Design Technology, a software company that commercializes an implementation of the inverse design method under the brand name TURBOdesign, among other related products.
Rather than directly establishing a trail blade geometry, in the inverse design method the user specifies the desired blade loading distribution, which is closely related to the resulting pressure distribution. The method therefore has the potential to render the design of Francis turbines and pump-turbines more physically intuitive: instead of defining blade angle distributions that are only indirectly related to the desired flow fields, the user specifies the work distribution to be achieved by the runner and the method calculates the blade angle distributions that accomplish the task. In doing so, the inverse design method has the potential to reduce the number of expensive CFD iterations required to achieve a satisfactory design, allowing for a considerable reduction of the time-to-solution.
The commercial implementation of the inverse design method has been successfully used in recent investigations focused on the design optimization of pump-turbines and Francis turbines. Zhu et al. [6,7] performed a multiobjective optimization of a reversible pump-turbine and investigated the effect of the blade lean angle on the hydraulic efficiency and on the magnitude of the pressure fluctuations. Daneshkah and Zangeneh [8] parametrically studied the influence of the blade loading distribution and stacking condition on the efficiency and cavitation performance of a Francis runner, leading to some design guidelines and considerable insight. Wang et al. [9] built on the aforementioned publications by Zhu et al. in order to investigate the trade-offs present in the optimum design of a pump-turbine; it was found that minimizing secondary flow losses improves the pump efficiency, whereas minimizing the profile losses increases the turbine efficiency, resulting in a set of Pareto-optimal designs between the two objectives. For a comprehensive review of the inverse design method and many of its recent applications in hydraulic turbomachinery, see Yang et al. [10].
In this context, the objective of the present work is to provide an open-source implementation of the inverse design method (see Supplementary Materials) that may contribute to further research into the optimal design of Francis turbines and pump-turbines. Indeed, although the method has been available in the literature for many years, no open-source implementation is currently available and not every researcher may afford a commercial license. More importantly, there are some mathematical inconsistencies between the articles of Borges [3] and Zangeneh [5], specifically in the final form of the equations in the curvilinear coordinate system which hinder independent implementations of the method. Apart from a detailed derivation of the final equations, this article presents a validation of the method with respect to high-fidelity CFD simulation results, and a thorough characterization of the method convergence.
The article is structured as follows: Section 2 presents an overview of the inverse design method, the governing equations in cylindrical and curvilinear coordinate systems, and several details about the solution algorithm; Section 3 describes the convergence analysis and validation of the method; Section 4 contains a discussion and the outlook.

Method Overview
The inverse design method assumes that the flow is steady, inviscid and incompressible, so it can be described by potential flow theory. The blades are represented by sheets of bound vorticity, i.e., the vorticity is not shed downstream. Furthermore, the flow is assumed to be irrotational at the inlet, with uniform radial and tangential velocity components, whereas at the outlet a uniform axial velocity profile is considered, in line with the requirement of no residual angular momentum downstream of the runner. These assumptions imply that the blades extract constant work along each streamline. A correction to account for blockage effects due to the blade thickness is introduced in the mean flow equations, although the blades are otherwise assumed to be infinitely thin.
The user input to the inverse design method includes: • The available head H, discharge Q, rotational speed ω, and number of blades B; • The meridional channel geometry: hub, shroud, leading edge and trailing edge curves; • The camber surface angle f , also termed blade wrap angle, along the blade leading edge; • The blade thickness distribution t n (r, z) normal to the camber surface; • The blade loading distribution λ(r, z).
The method outputs the blade shape that satisfies the inputs, including the required loading distribution, and the associated flow fields. The blade geometry is defined by its camber surface, itself determined by a distribution of the wrap angle f (r, z), and by the blade thickness distribution specified by the user. An example camber surface is presented in Figure 1 on the cylindrical coordinate system (r, θ, z) used in this work. Blade camber surface example on the cylindrical coordinate system (r, θ, z). The planes defined by θ = [0, π 2 ] illustrate the meridional projection of the blade leading edge (1) and trailing edge (2), as well as the curves that define the axisymmetric shroud (3) and hub (4) surfaces. The blade camber surface is defined by a set of coordinates (r, θ, z), with θ = f (r, z). More generally, the B camber surfaces are defined by θ = if (r, z), with i = 1, 2, ..., B.
In the current implementation, the meridional channel geometry is automatically calculated based on the runner specific speed and the hydraulic conditions, as detailed in Section 2.4.2. The normal blade thickness t n and blade loading λ are specified at the hub and shroud only, and then interpolated in the spanwise direction to compute their distribution on the blade.
The blade loading is defined as where r is the radius,m is the streamwise meridional coordinate, and C θ is the average tangential velocity, defined as The blade loading is therefore directly proportional to the variation of the flow angular momentum along a streamline. From Euler's equation, we have that ω ∂rC θ ∂m dm = gH where E is the transferred specific energy; consequently, in the current implementation of the method the user specifies a normalized blade loading distribution that is then scaled according to the available head. Similarly, the normal blade thickness is specified relative to the average blade chord c. Example distributions of blade loading and normal blade thickness are presented in Figure 2a,b, respectively. Both of these are specified in terms of the control points of non-uniform rational basis splines (NURBS).  The fluid velocity C (r, θ, z) is decomposed into an axisymmetric mean velocity C (r, z) and a tangentially-periodic velocity c (r, θ, z). The blade camber surface, defined by the wrap angle f (r, z), is updated iteratively, together with the flow field, until convergence is achieved. To begin the iterative procedure, an initial guess of f is computed based on the rC θ distribution, together with the assumption of uniform meridional velocity C m , itself computed from the discharge. As seen in Figure 2a, the Kutta-Joukowski condition is satisfied by setting ∂rC θ ∂m = 0 at the trailing edge, such that the pressure difference there is null, whereas the well-aligned flow condition is satisfied by setting ∂rC θ ∂m = 0 at the leading edge, meaning that the incidence angle is defined to be zero.

Governing Equations in Cylindrical Coordinates
The governing equations are hereafter presented in the form they take on the cylindrical coordinate system illustrated in Figure 1.

Mean Flow
The mean flow is described by Stokes' axisymmetric stream function Ψ (r, z) [11], which is found by solving the following elliptic equation where the right-hand side is equal to zero outside the blade region. The blockage factor B f , which accounts for the accelerating effect of the blade thickness on the mean flow, is defined as where t θ represents the blade thickness in the tangential direction, which is proportional to the user-defined normal blade thickness t n ; it is calculated as The mean radial and axial velocity components are defined in terms of the stream function as and Equation (3) is subject to the following boundary conditions. At the inlet and outlet, uniform radial and axial velocities are imposed, respectively, resulting in Neumann boundary conditions on Ψ which depend on the discharge and the inlet and outlet areas. At the walls the normal velocity is zero, which implies that the hub and shroud curves are streamlines; consequently they are defined by constant Ψ values. By convention Ψ = 0 at the hub, whereas at the shroud Ψ = Q 2π .

Periodic Flow
The periodic velocity is expressed using the Clebsch formulation [2] as where Φ (r, θ, z) is the potential function of the periodic flow and S (θ − f ) is the periodic sawtooth function with zero mean value, which introduces a velocity jump across the blade camber surface. A Fourier series in the tangential direction is used to express these functions as and respectively. Then, the governing equation for the nth harmonic of the potential function of the periodic flow, Φ n (r, z), can be shown to be where the right-hand side is equal to zero outside the blade region. These Helmholtz-type equations are subject to the following boundary conditions. At the inlet and outlet, where the flow is assumed to be uniform, i.e., the periodic velocity is null, a Dirichlet condition of the form Φ n = 0 is imposed. At the hub and shroud, where ∂c ∂n = 0, the corresponding condition ∂Φ n ∂n = 0 is imposed.

Blade Camber Surface
Having determined the mean and periodic flow components, the blade shape is calculated by enforcing that the camber surface be aligned with the flow field. This condition reads where the periodic velocity at the blade is defined as c bl = 1 2 (c + + c − ), i.e., taking the average between the periodic velocities at the pressure and suction sides of the blade, c ± .
This hyperbolic equation is solved using the method of characteristics, which amounts to integrating along the meridional projection of the flow streamlines on the blade surface. An initial condition is necessary, namely the user-defined camber surface angle f along the leading edge, sometimes termed the blade stacking condition.

Calculation of the Pressure
The pressure difference across the blade can be calculated as where ρ is the fluid density and W bl is the relative flow velocity at the blade, averaged between the pressure and suction sides. Note that the mean velocity increment caused by the blockage effect, which affects both blade sides, should have no impact on the pressure difference; for this reason the mean velocities C r,z disregarding the blockage effect are used, i.e., they are multiplied by B f . Indeed, only when disregarding the blockage effect is the torque calculated from the integration of the pressure difference over the blade equal to the torque calculated from the conservation of angular momentum. The momentum balance described by Euler's equation can be used to estimate the pressure gradient based on the mean flow velocity C as where g r,z accounts for the gravitational acceleration as well as the acceleration associated with the force exerted by the blades on the fluid. By integrating the pressure gradient, the pressure difference between any two points can be calculated. Taking the inlet midspan as the reference datum, and knowing that the pressure there is equal to p • = ρgH − 1 2 ρC 2 , the pressure at any point can be calculated as p (r, z) = p • + l 0 ∇p δl.

Governing Equations in Non-Orthogonal Curvilinear Coordinates
A finite difference scheme is used to solve the governing equations. A discretized meridional channel example is illustrated in Figure 3. The domain is meshed along curvilinear coordinates (ξ, η) that conform to the meridional channel geometry, including the blade leading and trailing edges, resulting in non-orthogonality. The governing equations on the non-orthogonal curvilinear coordinate system (ξ, θ, η) are hereafter derived. The transformation between cylindrical coordinates (r, z) and curvilinear coordinates (ξ, η) is represented by the Jacobian defined at each grid point. Its determinant is The Jacobian components can be used to predefine other transformation parameters that are going to appear in the governing equations, namely where the right-hand side is equal to zero outside the blade region; µ and ν are defined as and The blockage factor B f depends on the blade thickness in the tangential direction, that now reads The mean radial and axial velocity components can be respectively calculated as and The camber surface equation, presented further down, includes the mean streamwise and spanwise velocity components, C ξ and C η , respectively. These are calculated as and

Periodic Flow
Following the same procedure, i.e., application of the chain rule and extensive algebraic manipulation, the governing equation for the nth harmonic of the potential function of the periodic flow, Φ n (ξ, η), can be expressed as where ∇ 2 rC θ is calculated as whereasμ,ν and ζ are defined asμ and

Blade Camber Surface
The coordinate transformation results in the following form of the camber surface alignment condition used to calculate the blade shape: where the streamwise, spanwise and tangential components of the periodic flow velocity at the blade are computed, respectively, as and where N is the number of harmonics solved for, and Φ c,s n corresponds to the trigonometric decomposition of the imaginary Fourier coefficients Φ n .
These three equations are derived from Equation (8), noting that the sawtooth function cancels out when averaging between the pressure and suction sides of the blade, and using Euler's formula to express the imaginary Fourier coefficients in terms of their corresponding trigonometric counterparts.

Additional Implementation Details
The method has been implemented in MATLAB R2017b, although translating the code to other languages should be relatively straightforward. Three open-source scripts, related to the NURBS definition [12], logarithmic spacing [13], and plotting capabilities [14], have been used within the code. Hereafter some details of the implementation are presented.

Algorithm Overview
As illustrated in Algorithm 1, the code comprises an initialization stage, where the meridional channel geometry is generated and discretized, and an iterative solution stage, which includes three main steps: the solution of Equation (22) for the mean flow, the solution of Equations (30) for the periodic flow, and the integration of Equation (35) for the camber surface.

Meridional Channel Geometry Generation
The inverse design method requires the input of the meridional channel geometry, including the hub, shroud, blade leading edge and blade trailing edge meridional projections. The current implementation includes a function that calculates a trail meridional channel geometry that serves as a good starting point; the code can be extended to allow for the input of an arbitrarily parametrized geometry, although this is left for future work.
Algorithm 1: Implemented inverse design method. read user inputs generate meridional channel geometry generate finite difference mesh initialize fields and constant integration parameters: f , t n , rC θ , ∇rC θ , ∇ 2 rC θ ,μ,ν. while camber surface angle f has not converged do initialize variable integration parameters: Function ComputeMeanFlow() assemble and solve system matrix for Ψ calculate mean flow velocity C End Function ComputePeriodicFlow() foreach n ← 1 to N do assemble and solve system matrix for Φ n end calculate periodic flow velocity c End integrate blade camber surface equation to calculate update: f ← f + δf check whether camber surface angle f has converged end output runner geometry and flow fields The trail meridional channel geometry is calculated as a function of the available discharge Q, rotational speed ω, and the runner specific speed, defined as Bovet [15,16] proposed expressions that define the shape of the meridional channel as a function of n q ; the shape is then scaled according to a sizing parameter that depends on Q, ω and n q . The expressions that determine the meridional channel geometry contain 10 parameters that were fitted by Bovet to a database of Francis turbines designed by several manufacturers for a wide variety of operating conditions. In the current implementation, these parameters were updated somewhat by incorporating the more recent Francis turbine and pump-turbine designs that are reported by Henry [17]. Eight meridional channel geometry examples generated using the present approach are presented in Figure 4.

Discretization and Numerical Methods
As illustrated in Figure 3, the meridional channel domain is discretized with a mesh that conforms to the blade leading and trailing edges. Finite difference approximations are used to discretize the differential operators: ∂ ∂ξ , ∂ ∂η , ∂ 2 ∂ξ 2 and ∂ 2 ∂η 2 . Second-order accurate central differences are used for interior nodes, whereas second-order one-sided approximations are used for the boundary nodes [18]. Equivalent finite difference approximations are used to calculate the components of the Jacobian J, needed to compute the mesh metric coefficients, i.e., α, β, γ, J, τ , and σ.
The implicit algorithm used to solve the mean and periodic flow equations entails the solution of a large sparse linear systems assembled from the finite difference approximations of each term present in the respective governing equation. The iterative biconjugate gradient stabilized method (BiCSTAB) is used, together with an incomplete lower-upper (ILU) factorization of the system matrix as a preconditioner.
The blade camber surface equation is solved by integrating along the flow streamlines on the blade surface, which is to say along the characteristic curves of the equation. The camber surface angle f at the blade leading edge, i.e., the stacking condition, is used as initial condition for this integration. As evidenced in Figure 5a, the mean flow streamlines tend to recede from the hub and concentrate towards the shroud. Furthermore, during the initial iterations the periodic flow velocity is of considerable magnitude and may further disrupt the streamlines, making their distribution much less homogeneous. This implies that, when integrating equispaced streamlines starting at the blade leading edge, some blade regions may end up sparsely populated, resulting in significant interpolation error when calculating f at nodes that are relatively far from the closest integrated streamline.
To address this problem, the integration is performed in steps. For each of the nodes in a given grid line defined by constant η-coordinate, the streamline is integrated backwards using the second-order accurate Crank-Nicolson scheme until the upstream grid line is intersected; the calculated δf is used, together with the f value interpolated at the intersection with the upstream grid line, to compute the new f value at the node from which the streamline originated. This process is repeated from leading edge to trailing edge, as illustrated in Figure 5b, ensuring a uniform distribution of streamline sections and therefore avoiding the aforementioned interpolation error. Once the camber surface angle f new that respects the current flow field is calculated, a relaxation factor φ f is used to compute the updated camber surface angle as f = f new φ f + f old 1 − φ f . The relaxation factor ensures stability of the algorithm by damping the variation of the camber surface angle; φ f is adaptively selected according to where |∆f | l 2 and |∆f | l ∞ are respectively the l 2 -norm and the l ∞ -norm, with |x| l n = n ∑ i |x i | n , calculated over all the grid nodes on the blade, of the variation ∆f = f new − f old . This expression simply states that the relaxation factor is inversely proportional to the magnitude of the blade camber surface angle update, measured in terms of its average and maximum values over the blade; the parameters were fine-tuned based on extensive testing over several geometries. The relaxation factor φ f is typically close to 0.5. An additional stability measure that is implemented is to scale the periodic velocity by a factor φ c , which is smoothly relaxed from a value of φ c = 0.5 at the first iteration to a value of φ c = 1 from the tenth iteration onwards.
The solution convergence is assessed based on |∆f | l 2 over the blade grid nodes and | ∆C C | l 2 over all nodes. In other words, for a solution to be considered converged, the variation of both the blade camber surface angle and the velocity field between successive iterations must fall below an arbitrary threshold value, typically 0.1 • and 0.1%, respectively.

Convergence Analyses
In this section the convergence characteristics of the implemented method, both in terms of iterations and discretization resolution, are investigated with the objective of verifying the stability of the method and establishing the grid resolution requirements.
The iterative convergence behavior of the inverse design method is illustrated in Figure 6 for a representative Francis turbine case, using a mesh resolution that guarantees negligible discretization error. Figure 6a presents the l 2 -norm of the variations ∆f and ∆C C between successive iterations as a function of the iteration number. The convergence is mostly monotonic, except for a small bump close to the tenth iteration after which the full periodic velocity is considered by having φ c = 1. About 25 iterations are required to achieve a converged solution; a faster convergence is possible for this particular case by setting the relaxation factor φ f ≈ 1.0, although the more conservative adaptive approach is preferred in order to ensure convergence for all cases. Figure 6b evidences the monotonic convergence of the camber surface angle f at five points; the results are normalized by the respective angle value at the last iteration, f 35 . As expected, the points closer to the leading edge (m = 0) change significantly less than the points further down, in relation to their initialization value f 1 , given that the camber surface angle at the leading edge is fixed by the stacking condition selected. The convergence of the method with respect to the mesh resolution has been studied in detail. The grid resolution R is defined such that the number of nodes in the spanwise direction is 2 R + 1. For reference, the mesh illustrated in Figure 3 has R = 6, implying 65 nodes across the span.
A representative Francis turbine case is used for the grid convergence analysis. The pressure difference across the blade p + − p − . = ∆p as well as the blade angle β b = tan −1 r ∂f ∂m are selected to assess the convergence, given that these variables depend on the mean and periodic flow velocity components, as well as on the camber surface angle f . A total of 25 points uniformly distributed over the blade and four different resolutions R = (4, 5, 6, 7) are considered. Figure 7a presents the relative error for τ Ri = (∆p, β b ) Ri with respect to the corresponding values computed on the finest grid, τ R7 . The relative error for each of the 25 points analyzed is illustrated by thin lines, whereas the bold lines represent the average relative error, with error bars computed from the standard deviation. A clear convergence of the solution is evidenced, with the relative error monotonically decreasing with the mesh resolution. For the coarsest mesh the average relative error for the pressure difference is about 2%, whereas for the blade angle it is less than 1%.
Richardson extrapolation is used to compute the average discretization error. Any given solution metric can be described by the general model τ = k 1 + k 2 ∆x k 3 , where ∆x is the discretization size, k 3 is the convergence rate, k 2 is a constant and k 1 is the converged metric [18]; as such, the discretization error is equal to k 2 ∆x k 3 . This error model has been independently fitted to the 25 points analyzed, each of which includes the resulting ∆p and β b at four resolutions. Figure 7b presents the averaged results of the Richardson extrapolation analysis. The average convergence rate computed is k 3 = 1.79, whereas the average discretization error for the coarsest resolution is 2%. Note that the error presented in Figure 7b is relative to the finest resolution, which is not necessarily converged; on the contrary, the error presented in Figure 7b is with respect to the converged solution with ∆x → 0, implying that it truly is the magnitude of the discretization error. Even though all the operators in the present implementation are discretized with second-order accurate approximations, the effective convergence rate is ≈ 1.8, below the formal 2.0. This behavior is probably due to the fact that the number of harmonics N that can be resolved is proportional to the mesh resolution: the source terms of the governing equation for the potential function of the periodic flow, Equation (11), are imaginary exponential functions of frequency proportional to nB, with 1 ≤ n ≤ N , and therefore the highest order harmonic that can be resolved is subject to the Nyquist-Shannon sampling criterion. If higher order harmonics are used, the source terms suffer from aliasing and inject spurious information into the periodic velocity solution. In short, since the number of harmonics considered increases with the mesh resolution, the convergence rate is slightly below 2.0.
The time to solution, averaged over five runs, for each of the four resolution levels is presented in Table 1. The computations are performed on a single i7-4810MQ CPU core running at 3.8 GHz on Ubuntu Linux 16.04 LTS. It is evidenced that the time to solution increases very significantly with R, given that the number of equations being solved increases with the resolution, as explained above. Although these data are representative of the code performance, the time to solution does depend on other factors such as the turbine geometry: low specific speed runners with long blades tend to have more grid nodes for a given resolution level, since R only depends on the number of nodes along the span, and are therefore relatively more expensive to compute. Based on these results, it can be concluded that a resolution R = 5 can be used to efficiently explore the design space at a computational cost of about 10 s per simulation, whereas a resolution R = 6 can be used to compute the definite mesh-independent results in just over one minute. A considerable speed-up is expected if the code is ported to C++ or another high-performance programming language.

Implementation Validation
The inverse design method implementation is validated on a Francis turbine case characterized by a specific speed n q = 0.338 and a hydraulic power of 156.6 MW. The blade loading and blade thickness distributions specified are similar to the ones presented in Figure 2, whereas the mesh resolution is selected as R = 6 to ensure converged results; the corresponding meridional channel discretization is illustrated in Figure 3. The number of blades is set to 16, and a quadratic camber surface angle distribution is specified at the leading edge as the stacking condition, with a ∆f of 8.2 • at the shroud with respect to the hub. This stacking condition corresponds to an average leading edge slant of where r LE is the leading edge radius at the shroud and z LE is the leading edge height. A render of the resulting runner is shown in Figure 8, whilst Figure 9 illustrates the distributions of the blade camber surface wrap angle f and the blade angle β b = tan −1 r ∂f ∂m . The runner geometry output by the method is evaluated using the commercial CFD software ANSYS CFX 19.2, where the Navier-Stokes equations coupled with the k − ω SST turbulence model are solved in steady-state. Only one blade-to-blade channel is simulated. A mesh with 8.54 million elements and an average dimensionless wall distance over the blade of y + = 15.1 is employed. The velocity orientation and magnitude is prescribed at the inlet, and an average static pressure of 0 Pa is imposed at the outlet.
This simulation configuration is commonplace in the field of hydraulic turbomachines, both in academia and in industry, and has been shown to provide numerical results that are in good agreement with experiments [19][20][21][22][23], at least at the best efficiency operating conditions. Using this well-established CFD methodology is adequate to validate the correctness and physical soundness of the solution obtained with the implemented inverse design method, although it is acknowledged that both numerical solutions will have some degree of discrepancy with respect to physical reality.

Qualitative Validation
First the velocity fields are used to qualitatively validate the output of the inverse design method implementation. On the following figures, even though the color scale is not identical in MATLAB and CFX, the same variable range is enforced and both color bars are presented. To further aid in the assessment of the results, the meridional projections of the blade leading and trailing edges, as well as of the hub and shroud, are superimposed on the CFX flow fields. There is a significant agreement for the mean flow velocity magnitude C, as evidenced in Figure 10. Both the main flow features and the velocity magnitude are well approximated by the present method, compared to the fully-3D Navier-Stokes solution, in spite of the underlying inviscid and axisymmetric mean flow assumptions. The same can be said about the radial and axial mean velocity components, presented in Figures 11  and 12, respectively.
The greatest discrepancy is evidenced in the flow near the shroud, where the mean axial velocity is somewhat overpredicted by the present method. This behavior is most likely a consequence of the inviscid flow assumption. Another divergence occurs for the mean radial velocity prediction near the intersection between the leading edge and the shroud, where the magnitude of the flow acceleration in this high-curvature region is slightly underpredicted by the present method. Overall, a good agreement between both approaches is evidenced.   The comparison of the mean tangential velocity C θ is presented in Figure 13. This velocity component illustrates the manner in which the inflow angular momentum is extracted by the runner. As such, the notable agreement achieved validates the proposed methodology by demonstrating that the designed runner accomplishes an essentially complete extraction of the flow angluar momentum. Moreover, not only is the momentum extraction thorough, it is also performed according to the prescribed blade loading distribution, as revealed by the equivalent gradient of C θ on both solutions.

Quantitative Validation
For a quantitative validation of the inverse design method implementation, a comparison of the pressure distribution on the blade along five representative span locations is performed.
The pressure difference ∆p across the blade is presented in Figure 14. As revealed by Equation (13), the pressure difference is proportional to ∇rC θ , i.e., proportional to the blade loading, and to the meridional projection of the relative flow velocity, which is greatest towards the shroud. Both of these factors are encompassed in the pressure difference distributions obtained.  There is an overall agreement between the present method results and those calculated with the Navier-Stokes solver, although there are also some noticeable discrepancies that are probably due to a combination of factors, including some error in the computed flow velocity, the assumptions inherent to Equation (13), and the simplifications of the present model. Perhaps the assumption of infinitely thin blades has the greatest impact: it is likely responsible for the underprediction of the pressure difference in the first 8% of the chord, i.e., form < 0.08; the method then tends to compensate for this by an overprediction of the pressure difference for 0.20 <m < 0.40, especially near the shroud. Based on the five profiles presented in Figure 14, the average relative error of the pressure difference computed with the present method amounts to 12%, whereas the relative error in the computed torque is equal to 3%. Since the torque is calculated by integration of the pressure difference, a smaller error is expected because, to a certain extent, regions where the pressure difference is underpredicted will cancel out the overcontribution of regions where it is overpredicted. Moreover, at least part of the 3% torque overprediction by the current method is explained by the neglect of all energy losses in the flow.
A comparison of the pressure profiles over the blade is presented in Figure 15. Whereas in the Navier-Stokes solver the pressure is actually computed on each of the blades surfaces, in the current implementation Equation (14) is used to estimate the average meridional pressure p, and then the computed pressure difference across the blade is symmetrically superimposed: p ± blade = p ± 1 2 ∆p. Although there is some agreement in the pressure profiles over the blade between the inverse design method implementation and the Navier-Stokes solution, there are at least two significant discrepancies. The first one has to do with the leading and trailing edges: On the one hand, as already discussed, the assumption of infinitely thin blades hinders the rapid build-up of the pressure at the leading edge, and renders the formation of a stagnation point impossible. On the other hand, there is a pressure jump at the trailing edge near the shroud that is not supposed to be there according to the simplified model. Two possible explanations for that discrepancy are the accumulation of viscous effects or the influence of secondary structures, although no conclusive explanation has been found.
The second significant discrepancy has to do with the repartition of the pressure difference between the suction and pressure sides. It is evident that the simplified approach used to calculate the pressure on each surface, namely a symmetric repartition of the pressure difference, is only approximative. The Navier-Stokes solution suggests that, relative to the mean meridional pressure, the pressure drop on the suction side is greater than the increment on the pressure side.
The energetic efficiency of the runner, computed using the Navier-Stokes solution, is equal to η e = 97.8%. This value does not consider the energy losses in the guide vanes, which were not modeled, and by definition it does not include the disk friction or leakage discharge losses either. Acknowledging this, the calculated efficiency is considerably high, implying that the implemented inverse design method is actually capable of providing a good hydraulic design that follows the blade loading distribution prescribed as input.

Discussion and Outlook
Even though the computational cost of the ANSYS CFX simulation is about 7000 times greater than the present method, the latter is clearly not expected to replace the former: the Navier-Stokes simulations are still necessary to calculate the turbine efficiency, to assess the risk of cavitation, and to study the flow behavior at off-design operating conditions. On the contrary, the two approaches may be complementary: the present inverse design method may be used to quickly identify suitable runner designs that would then be evaluated thoroughly by means of the more expensive and reliable solver.
The blade angle β b distribution provided by the inverse design method, presented in Figure 9, is non-monotonic and varies considerably along the span, suggesting that it would not be straightforward to arrive at this geometry by direct parametric manipulation of the blade angles, i.e., by the direct design method. Furthermore, the manner in which these angles affect the flow is unintuitive, whereas the effect of the blade loading on the pressure difference distribution and flow velocity components is more intelligible. It might therefore be advantageous to integrate the present inverse method into the design methodology of Francis runners.
The main assumptions of the model, namely considering the flow as inviscid and the blades as infinitely thin on the periodic flow equations, are expected to introduce only minor error into the design: The blade camber surface angle computed by the method mostly depends on the mean velocity components, which are well predicted as evidenced in Section 3.2.1; this accuracy is partly due to the blockage factor B f , which effectively corrects these velocity components by accounting for the blade thickness. The most important discrepancies were found on the pressure distribution estimated by the method, which has no effect on the output geometry. The pressure distribution is subject to the error introduced by the approximative nature of Equation (14), and by the lack of a stagnation region on the leading edge, itself a consequence of the assumption of infinitely thin blades.
One possible improvement to the method would be to implement a better approximation of the pressure distribution. In the momentum balance used to compute the meridional projection of the pressure gradient, Equation (14), only the axisymmetric mean velocity components were considered, i.e. the gradients in the tangential direction were neglected. Nevertheless, the flow solution does include the complete 3D velocity field by means of the tangential harmonic decomposition, so it is in fact possible to compute all the terms in the momentum balance to derive a better approximation of the pressure gradient. The problem is, however, that the velocity gradients in the tangential direction are subject to the noise introduced by the Gibbs phenomenon near the blade, so a filtering of the spurious oscillations would prove necessary.
Another possible improvement would be to introduce more flexibility and generality into the design method. In the current implementation, it is assumed that the blade leading edge is well aligned with the incoming flow, that the angular momentum at the trailing edge is null, and that the extracted work is constant along the blade span. Although these assumptions are reasonable, they limit the design space and may prevent a widespread use of the method in real applications.
A third development direction lies in the input and output capabilities of the code. It would be beneficial to allow the user to specify an arbitrarily parametrized meridional channel geometry, and to facilitate several output geometry formats compatible with CAD and CFD software, since the current implementation only allows exporting nodal coordinates. Even if these coordinates are sufficient to reconstruct the surfaces and assemble the geometry to be used in the CFD evaluation, this process is unnecessarily cumbersome and could be significantly shortened.