Ducted hydrokinetic turbine

Abstract

A ducted turbine (10) includes a turbine (12) having a plurality of blades (14) extending from a hub (36) that rotate together about a turbine axis (X). The turbine (12) is surrounded by a coaxial duct (16) having an inner surface (22) radially outboard of the turbine. The duct (16) has a diameter (D) that changes along the axial direction such that its inner surface (22) has a curvilinear profile in an axially extending plane. The duct (16) accelerates and conditions fluid flowing past the turbine blades (14) to enhance energy extraction efficiency. The geometry of the duct (16) and the turbine blades (12) can be optimized using computational fluid dynamics (CFD) with a blade-resolved Reynolds-averaged Navier-Stokes (RANS) solver. To handle the high-dimensional design space of the blade and duct geometry, a gradient-based optimization can be performed with gradients computed by the adjoint method. The optimized geometries can provide the ducted turbine (10) with an efficiency significantly higher than that of a typical unducted turbine.

Description

[0001]DUCTED HYDROKINETIC TURBINE This invention was made with government support under DE-AR0001438 awarded by the U.S. Department of Energy. The government has certain rights in the invention. TECHNICAL FIELD This disclosure is related to turbines that convert fluid movement to rotational movement and, in particular, to hydrokinetic turbines designed to extract energy from a flowing liquid such as water. BACKGROUND The increasing demand for renewable energy has motivated extensive research on hydrokinetic energy conversion systems that extract energy from natural riverine and oceanic flows. Various types of conversion systems have been investigated, including oscillating hydrofoils and turbines with horizontal or vertical axes. Horizontal-axis turbines have been studied the most because of the relatively mature technology. A popular benchmark for horizontal-axis hydrokinetic turbines is the Bahaj model, which has been experimentally tested in a cavitation tunnel and a towing tank. The unducted Bahaj model generates power with an efficiency of about 46%, which is the ratio of generated power to the inflow power, at the optimal operating condition. This is the typical efficiency of well-designed hydrokinetic turbines. One theoretical limit on turbine efficiency is known as Betz's limit, which is derived based on the one-dimensional momentum theory in an unbounded flow domain. Betz’s limit on turbine efficiency is 59.3%, but there is no general consensus on whether this should be considered as a hard upper bound on the efficiency of practical energy conversion systems in unbounded flow due to simplifications in the theory. Research on ducted hydrokinetic turbines has been quite limited due to the complexity of the turbine-duct interaction, which involves flow features that cannot be captured by analytical models. For example, the available analytical models cannot account for flow separation, which can significantly affect the mass flow through the duct. Additionally, any research conducted on ducted wind turbines is largely inapplicable to hydrokinetic turbines, the blades of which are shorter and require a larger sectional thickness to sustain the higher loads in water. As a result, most or all of the available information on ducted hydrokinetic turbines is merely the result of empirical data on turbines designed by human intuition or a grid search in a low-dimensional design space. SUMMARY Embodiments of a method of optimizing the geometry of a ducted turbine including a duct and turbine blades within the duct include: (a) determining an efficiency of a low-fidelity model of the ducted turbine; (b) changing the geometry; (c) determining a change in the efficiency of the model after step (b); and (d) repeating steps (b) and (c) until the change in the efficiency is below a threshold value to arrive at an optimized geometry for the ducted turbine. The method may additionally include any technically feasible combination of one or more of the following features: step (b) comprises changing at least one of the following: a duct length, a duct diameter, a blade twist angle, a blade cross-section; the ducted turbine includes a turbine hub within the duct with the turbine blades extending from the turbine hub, step (b) comprising changing at least one of the following: a hub length, a hub diameter, a duct thickness; an outlet diameter of the duct and a tip-gap ratio of the ducted turbine remain unchanged during step (b); the low-fidelity model includes a blade-fixed rotating reference frame enclosing the turbine blades and a stationary inertial reference frame outside the duct; use of a first-order numerical scheme to compute convective terms of Reynolds- averaged Navier-Stokes (RANS) equations; parameterizing the geometry of the ducted turbine before step (a) using free-form deformation (FFD), wherein control points defining FFD boxes are moved relative to one another in step (b) such that a computational mesh defined in and around the model is deformed when the geometry is changed; the FFD boxes include a first FFD box enclosing the duct and excluding the turbine blades, second FFD boxes enclosing each of the turbine blades and excluding the duct, and a third FFD box enclosing the first and second FFD boxes; and determining the efficiency of a high-fidelity model of the ducted turbine having the optimized geometry of step (d). Embodiments of a method include: (a) determining an optimized geometry for a turbine using a low-fidelity model of the turbine, and (b) determining an efficiency of the turbine using a high-fidelity model of the turbine with the optimized geometry. The turbine may be part of a ducted turbine. Step (a) may include solving steady RANS equations using a multiple reference frame (MRF) model, and / or step (b) may include solving unsteady RANS equations using a rotating- sliding mesh (RSM) model. Embodiments of a ducted turbine include a duct and turbine blades within the duct. The duct has a length, a throat diameter, and an exit diameter, and each turbine blade has an effective radius, a blade twist profile, and a blade section profile. The length, the throat diameter, the blade twist profile, and the blade section profile are configured to provide the ducted turbine with an efficiency greater than 46% as determined by solving unsteady Reynolds-averaged Navier–Stokes (RANS) equations using a rotating-sliding mesh (RSM) model. The ducted turbine may additionally include any technically feasible combination of one or more of the following features: the length of the duct is less than three times the effective radius of each turbine blade; the length of the duct is less than two times the effective radius of each turbine blade; the throat diameter is greater than 65% of the exit diameter; the throat diameter is greater than 70% of the exit diameter; each blade twist profile has a first twist angle in a range from 20 degrees to 30 degrees at r / R = 0.4 and a second twist angle in a range from 10 degrees to 20 degrees at r / R = 1, where r is a radial distance from a rotational axis of the turbine blades and R is the effective radius of the turbine blades; the first twist angle is in a range from 25 degrees to 30 degrees and the second twist angle is in a range from 10 degrees to 15 degrees; each blade twist profile has a first twist angle greater than 20 degrees at r / R = 0.4 and a second twist angle less than 20 degrees at r / R = 1, where r is a radial distance from a rotational axis of the turbine blades and R is the effective radius of the turbine blades; and a hub from which the turbine blades extend within the duct extends from within the duct to a nose extending beyond an inlet of the duct. BRIEF DESCRIPTION OF THE DRAWINGS FIG. 1 schematically illustrates an unducted turbine operating in a uniform inflow in an unbounded fluid domain. FIG.2 schematically illustrates a ducted turbine operating in the same uniform inflow in an unbounded fluid domain. FIG.3 schematically depicts an illustrative turbine blade having a baseline geometry with a foil-shaped cross-section and twist angles defined at multiple locations along the blade that are permitted to vary in an optimization method. FIG.4 schematically depicts the turbine blade of FIG.3 and illustrates cross-sectional areas of the foil-shaped cross-section defined at multiple locations along the blade that are permitted to vary in the optimization method. FIG. 5 schematically depicts an illustrative ducted turbine with a duct having a baseline geometry with diameters defined at multiple locations along the duct that are permitted to increase or otherwise vary in the optimization method. FIG.6 schematically depicts the baseline geometry of FIG.5 with the diameters permitted to decrease or otherwise vary in the optimization method. FIG. 7 schematically depicts the baseline geometry of FIGS.5 and 6 with a length that is permitted to increase or otherwise vary in the optimization method. FIG.8 schematically depicts the baseline geometry of FIGS.5-7 with the length permitted to decrease or otherwise vary in the optimization method. FIG.9 is a flowchart depicting an illustrative method including optimization followed by reevaluation. FIG.10 depicts a baseline ducted turbine including a duct and turbine blades for use in the optimization method. FIG.11 depicts twist angle as a function of radial location along the turbine blades of FIG. 10 with two different baseline geometries. FIG. 12 depicts a multiple reference frames (MRF) approach for use in the optimization method. FIG. 13 depicts a rotating-sliding mesh (RSM) approach for use in the optimization method. FIG. 14A is an axial view of an unstructured computational mesh for use in the optimization method from an inlet end of the ducted turbine. FIG.14B is a side view of the unstructured computational mesh of FIG.14A. FIG.14C is an enlarged view of a portion of FIG.14B. FIG.15A illustrates an example of a two-dimensional surface to be deformed in a manner controlled by free-form deformation (FFD) with control points constructed to enclose the surface. FIG.15B illustrates an example of the surface of FIG.15A after FFD deformation. FIG.16A illustrates an example of a three-dimensional surface to be deformed in a manner controlled by FFD with control points constructed to enclose the surface. FIG.16B illustrates an example of the surface of FIG.15A after FFD deformation. FIG.17A is a perspective view of an FFD setup including a parent box enclosing a duct and turbine blades of a ducted turbine and a child box enclosing the duct. FIG.17B is a side view of the FFD setup of FIG.17A. FIG.17C is an end view of the FFD setup of FIG.17A. FIG. 17D is a perspective view of a portion of the FFD setup of FIG. 17A illustrating a child box enclosing one of the turbine blades. FIG.17E is an axial end view of a portion of the FFD setup of FIG.17A illustrating child boxes enclosing each of the turbine blades. FIG.18 illustrates a portion of the undeformed computational mesh of FIG.14C. FIG. 19 illustrates an example of the computational mesh of FIG. 18 after deformation associated with radial expansion and axial shortening of the duct. FIG. 20 includes power coefficient curves for an unducted Bahaj turbine obtained from RANS-MRF and URANS-RSM compared with experimental results. FIG. 21 illustrates power coefficients during optimization of an unducted turbine having two different turbine blade baseline, along with the power coefficients as reevaluated after optimization. FIG.22 includes power coefficient curves for the two geometries optimized in FIG.21, for the Bahaj model of FIG.20, and for an optimized unducted turbine. FIG.23 illustrates two optimized duct geometries relative to the baseline duct geometry. FIG. 24 includes plots of twist angle and chord length for two optimized turbine blade geometries relative to their respective baseline geometries. FIG.25 illustrates variations in the power coefficient, flow rate, duct thrust coefficient, and blade thrust coefficient of a ducted turbine geometry during optimization. FIG. 26A is a flow visualization of streamlines in a flow tube passing the blades of an optimized unducted turbine. FIG.26B is the flow visualization as in FIG.26A for a baseline ducted turbine geometry. FIG.26C is the flow visualization for the baseline geometry of FIG.26B after optimization. FIG.27 illustrates baseline and optimized geometries of a ducted turbine duct and a turbine hub. FIG.28A depicts a non-conformal donor mesh and a target mesh. FIG.28B depicts a connection between the donor mesh of FIG.28A and a supermesh. FIG.28C depicts a connection between the target mesh of FIG.28A and the supermesh of FIG.28B. FIG.29A depicts a portion of a grid search attempting to increase efficiency of a baseline duct of a ducted turbine including five different duct outlet areas. FIG.29B depicts another portion of the grid search including five different duct outlet areas at the duct throat. FIG. 29C depicts another portion of the grid search including five different sets of duct inlet areas, outlet areas, and finite thicknesses. FIG. 29D depicts another portion of the grid search including five different duct outlet areas and a flange at the duct outlet. FIG.29E depicts another portion of the grid search including five different duct inlet areas. FIG.30 depicts calculated efficiencies for the twenty-five duct geometries of FIGS.29A- 29E with multiple different turbine blade geometries and tip-gap ratios. DESCRIPTION OF EMBODIMENTS Described below is a ducted hydrokinetic turbine with improved efficiency over unducted turbines such as the Bahaj model, along with a method of optimizing the geometry of the duct, the turbine blades, and other features such as the hub from which the turbine blades extend. The method employs computational fluid dynamics (CFD) for optimization. In the examples described below, twenty-one individual parameters are used to optimize and describe the shape of the duct and the turbine blades, including duct length and multiple duct radii, along with pitch, spanwise twist, and chord distributions of the turbine blades. The method also employs gradient-based optimization with gradients computed by a discrete adjoint method coupled with steady Reynolds- averaged Navier-Stokes (RANS) blade-resolved simulations. FIG.1 schematically illustrates an unducted turbine 12’ operating in a uniform inflow ^^^in an unbounded fluid domain. The turbine 12’ includes a plurality of turbine blades 14’ that may be interconnected by a central hub (not shown) for rotation together about a turbine axis X. The turbine 12’ converts inflow power (energy) into rotational power. The effectiveness of this conversion is characterized by a power coefficient, ^^^^ ൌభ(1) ఘ^^ where ^^^is the power coefficient or power given by ^^ ൌ ^^Ω, in which Q is torque and Ω is the rotational speed of the turbine. The coefficient ^^^is non-dimensionalized by the fluid density ^^ and a reference area ^^ ൌ ^^ ^^ଶ. Generally, for purposes of calculating turbine efficiency, the reference area A can be chosen as either the blade swept area Ab or the maximum projected area of the device Amax, which are identical in the case of the unducted turbine 12’ of FIG.1. FIG.2 schematically illustrates a ducted turbine 10 operating in the same uniform inflow ^^^in an unbounded fluid domain. The ducted turbine 10 includes a rotatable turbine 12 having a plurality of turbine blades 14 that may be interconnected by a central hub (not shown) for rotation together about the turbine axis X. The ducted turbine 10 further includes a duct 16 having an inlet end 18, an outlet end 20, and an inner duct surface 22 extending between the inlet and outlet ends. The duct 16 and its inner surface 22 are coaxial with the rotatable turbine 12 and extend circumferentially around the turbine 12. The duct surface 22 is radially outboard of the turbine blades 14 and has a diameter D that changes along the axial direction between the inlet and outlet ends 18, 20 such that the duct inner surface has a contoured profile in an axially and radially extending plane. The turbine 14 may be located along a radially extending plane at a throat T of the duct 16, at which the inner surface 22 is at its minimum diameter. The diameter D of the duct 16 as defined at the inner surface 22 decreases from the inlet end 18 to the throat T and increases from the throat to the outlet end 20. The portion of the duct 16 between the inlet end 18 and the throat T may be referred to as the intake 24, and the portion of the duct between the outlet end 20 and the throat may be referred to as the diffuser 26. The power coefficient ^^^of the ducted turbine 10 is also governed by Equation (1). However, in the case of the ducted turbine 10, the blade swept area Abis different from the maximum projection area of the device Amax, which leads to a different power coefficient depending on which value is used. For purposes of this disclosure, the reference area A of the ducted turbine 10 is ^^ ൌ ^^୫ୟ^. ^^୫ୟ^is the appropriate choice for a ducted turbine in order to have a meaningful comparison of its performance with unducted turbines. Using ^^ ൌ ^^୫ୟ^effectively offers a comparison of generated power when the inflow power is the same for unducted and ducted turbines. On the other hand, using ^^^for a ducted turbine results in a larger value of ^^^(sometimes even above 1) that can be misleading when compared to the efficiency of an unducted turbine. The power coefficient ^^^is also generally a function of two other non-dimensional parameters—namely the tip-speed ratio λ and the Reynolds number Re, ^^ ൌஐோand ^^ಮ^^ೌ^^^ ^^ ൌఔ, (2) where R is the radius of the ^^ is the fluid kinematic viscosity. For purposes of the examples presented herein, ^^^ൌ 1.4 m / s, ^^ ൌ 1 ൈ 10ି^mଶ⁄s , ^^ ൌ 1.853 mଶ, and ^^^^௫ൌ ^^4 / ^^^ ^^ ൌ 1.536 m, such that ^^ ^^ ^ 2 x 106for both ducted and unducted number, the fluid flow is considered fully turbulent, and the dependence of ^^^on Re in this range is expected to be relatively weak. Below, ^^^is evaluated for a broad range of tip-speed ratios λ at this Reynolds number for both ducted and unducted turbines. The rotational speed of the turbine 14 is held constant at Ω = 17.5 rad / s for purposes of each optimization. The optimization method presented herein can of course be employed with different fluid flowrates, kinematic viscosities, reference areas, and turbine speeds within practical limits. While the method is demonstrated below as applied to a ducted turbine, the same method can be used to optimize an unducted turbine, but with a simpler setup. FIGS.3 and 4 schematically depict illustrative turbines blades 14 having a base or proximal end 28 and extending to a tip or distal end 30. Each blade 14 has a foil-shaped cross-section defining a leading edge 32 and a trailing edge 34 and a blade axis Z extending through a centroid of the foil-shaped cross-section. The blade axis Z is not necessarily a straight line. In each of FIGS. 3 and 4, a baseline turbine blade 14* is illustrated in phantom as a baseline blade geometry to which modifications are made for purposes of the optimization. FIG.3 illustrates the case where a blade twist angle ^^^is defined at a plurality (n) of radially spaced locations i along the turbine blade 14 relative to the baseline blade 14*, with no change in the size or shape of the blade cross- section relative to the baseline geometry. In the illustrated example, n = 8 with ^^^defined at the base 28 of the blade 14 and ^^଼defined at the blade tip 30. Each radially spaced location i along the turbine blade 14 may be identified by a ratio r / R, where r is the radial distance from the turbine axis X to the respective location i. Each radially spaced location i also has an associated chord length ci—i.e., the distance between the leading and trailing edges 32, 34. FIG.4 illustrates the case where a blade cross-sectional area ^^^is defined at the plurality (n) of radially spaced locations i along the turbine blade 14 relative to the corresponding area ^^^^of the baseline blade 14*, with no change in the twist angle or foil shape relative to the baseline blade. The cross-sectional area ^^^of the blade 14 at each location i is either expanded ^ ^^^ / ^^^^^ 1^ or contracted ^ ^^^⁄ ^^^^^ 1^ relative to that of the baseline blade 14* during optimization. The twist angle and areal changes relative to the baseline blade 14* are illustrated separately in FIGS.3 and 4 for simplicity in explanation and to set-up the optimization problem. It should be understood that one or both of the twist angle θ and cross-sectional area b may vary relative to the baseline blade 14* for each identified radial location i along the turbine blade 14. FIGS. 5-8 schematically depict illustrative ducted turbines 10, including respective turbines 12 and ducts 16 including the features discussed in conjunction with FIG.2. In each of FIGS. 5-8, a baseline duct 16* is illustrated in phantom as a reference duct geometry to which modifications are made for purposes of the optimizations. FIGS.5 and 6 illustrate the case where a duct diameter ^^^is defined at a plurality (m) of axially spaced locations j along the inner surface 22 of the duct 16 relative to that of the baseline duct 16*, with no change in the length l of the duct 16 relative to the baseline duct. In the illustrated example, m = 4, d1 is at the inlet end 18, d3 is at the throat T, and each of d2 and d4 is between the throat and one of the ends 18, 20 of the duct 16. In FIG.5, the duct 16 is radially expanded relative to the baseline duct 16*. In FIG.6, the duct 16 is radially contracted relative to the baseline duct 16*. During radial expansion or contraction of the duct 16, the length R of the turbine blades 14 is respectively increased or decreased from that of the baseline blades 14* to maintain a constant tip-gap ratio 2R / dT at the duct throat T. In FIGS. 5 and 6, dT = d3. FIGS.7 and 8 illustrate the case where the length l of the duct 16 is increased or decreased with respect to a length lBof the baseline duct 16*, where length is defined as the axial distance between the inlet and outlet ends 18, 20 of the respective duct 16. For purposes of lengthening or shortening the duct 16 relative, the diameter of the inner surface 22 is bounded between the diameter dTat the throat T and Dmax= Dexit. In FIG. 7, the length l of the duct 16 is increased relative to the length lBof the baseline duct 16* (l > lB). In FIG. 8, the length l of the duct 16 is decreased relative to the length lBof the baseline duct 16* (0 < l < lB). The ratio lT / l is held constant based on that of the baseline duct 16* (lT / lB), where lTis the axial distance between the throat T and the inlet end 18 of the duct 16. The radial expansion or contraction and the axial lengthening or shortening relative to the baseline duct 16* are illustrated separately in FIGS.5-8 for simplicity in explanation and to set-up the optimization problem. It should be understood that one or both of the duct length l and the diameter dj of each identified location j may vary relative to the baseline duct 16*. A constrained optimization problem (3) for the ducted turbine 10 can be stated as: maximize ^^^(3a) by varying^^^ , ^^^^ ^^^ ^ୀ^ , ^ ^^^^^ୀ^, ^^ (3b) subject to 0.8 ^^^^ಳ^^ 1.2, (3d) Here,^^^^^^ୀ^are the twist angles at the illustrated in FIG.3, where n = 8;^^^^^^^ୀ^are the cross-sectional areas of the turbine blade at the plurality n of radial locations i along the blade as illustrated in FIG.4, where n = 8; ^ ^^^^^^ୀ^are the diameters at the inner surface 22 at the plurality m of axially spaced locations j along the duct 16 as illustrated in FIGS.5 and 6, where m = 4; and l is the length of the duct 16 as illustrated in FIGS. 7 and 8. The inequalities (3c), (3d), and 3(f) provide respective bounds on the blade twist angles ( ^^^^^, ^^^^௫), the blade section areas ( ^^^^^, ^^^^௫), and the length ( ^^^^^, ^^^^௫) of the duct 16. The bounds on the blade section area and duct length are provided as bounds on an areal scaling factor ^^⁄ ^^^and the duct lengthening factor ⁄^ ^^^^ ^^ . The cross-sectional at the plurality n of locationsi along the baseline blade 14* can be expressed as^^^^^^^^ୀ^ , such that ^^^⁄^^^^defines the scalingfactor at each location i. The bounds in (3c), (3d), and (3f) may be set up in a trial-and-error manner to ensure that the optimized variables do not touch the bounds. The full turbine geometry should be considered to morph smoothly with the variation of the design variables. This may be implemented through a Free-Form Deformation method that is discussed later. The inequality (3e) ensures that the diameter dj at every axial location along the duct 16 (excluding the throat T) is greater than or equal to the diameter dT at the throat and less than or equal to the diameter Dexitat the duct outlet 20. For the examples of FIGS.5 and 6, where dT= d3and m = 4, 0 ^ ^^ଷ^ ^ ^^^^^ୀ^,ଶ,ସ^ ^^^^୧^, ensuring that d3is always located at the throat T and that no duct diameter diameter ^^exitൌ ^^^^௫ൌ ^^4 / ^^^ ^^. The diameter Dmax at the outlet end 20 of the duct 16 is used as the maximum duct to the need for a large exit area that reduces flow velocity through the expansion, which in turn increases the flow momentum extraction at the turbine blades. The constraint (3g) keeps the tip-gap ratio 2 ^^⁄^^்constant throughout the optimization based on the baseline blades 14* and duct 16*. In this example, the tip-gap ratio, which is the ratio of the effective diameter (2R) of the turbine 14 to the diameter d3at the throat T, is set at 91%. This leads to the radius R of the turbine blade 14 changing with any variation in the throat diameter due to radial expansion or contraction from the baseline duct 16*. For purposes of the optimization examples, lT is set at 26.4% of the duct length based on the baseline duct 16*. A gradient-based optimization method may be applied to the solve problem (3) set forth above. FIG.9 is a flowchart depicting an illustrative method 100 including an optimization process 110 followed by a high-fidelity reevaluation process 130. These processes may be implemented via open-source software, such as DAFoam for the optimization process 110 and OpenFOAM for the reevaluation process 130, or other suitable software. DAFoam and OpenFOAM were used for the examples presented in this disclosure. The optimization process 110 includes defining a baseline geometry for the turbine (step 112) and parameterizing surfaces of the turbine geometry (step 114), followed by an iterative process 116. Parameterizing includes defining surface features of the turbine geometry in a manner such that a change in one of the design variables causes a corresponding change in the fluid flow field in which the turbine is being modeled. One parameterizing technique used in this disclosure is free-form deformation (FFD), which is subsequently described in further detail. The iterative process 116 generally includes determining the efficiency Cpof a model of the turbine (step 118), changing the turbine geometry (step 120), determining a change in Cp(step 122) after changing the turbine geometry, and then repeating steps 120 and 122 until the change in Cpis below a threshold value—i.e., until the efficiency reaches a plateau within pre-defined limits, such as when the efficiency changes by less than 0.1% or less than 0.05% relative to the previous efficiency calculation. As illustrated in FIG. 9, the iterative process 116 may include a gradient computation 124. In the gradient computation, the efficiency Cp is determined as a function of all of the identified design (geometry) variables in the vicinity of the current value of the design parameters. In other words, the effect of each identified design variable on the efficiency Cp is established. This informs the next iterative change in geometry. In a simple example with only a single design variable, if an increase in the duct throat diameter in the most recent geometry change results in a higher calculated efficiency—a positive gradient—the throat diameter is incrementally increased in the next geometry change. Of course, with 21 design variables, the gradient computation is much more complex and may be performed using the adjoint method as described in greater detail below. Once the efficiency Cp plateaus at step 122, the computational mesh may be reconstructed around the model in a remeshing step 128 to perform further geometry changes and gradient computations, unless the model has already been previously remeshed and no change in Cpresulted with the new computational mesh, as provided for in step 126. The continued and cumulative changes in the turbine geometry during the iterative process can lead to increased computational error as the computational mesh becomes progressively more deformed, sometimes leading to a false plateau in calculated turbine efficiency. Remeshing around the most recent geometry and restarting the iterative process can lead to further increases in efficiency with additional geometry changes. Once a plateau in Cpis reached and remeshing results in no further increases, the turbine geometry is optimized. Once optimized, the identified geometry is used in a higher-fidelity model to more accurately determine the efficiency of the turbine with the optimized geometry. This dual-tiered modeling is unconventional on its own and represents an advancement in the art by addressing one of the problems with modeling the complex geometries of turbines. Modeling unducted wind turbines is known to be very complex and difficult to perform with accuracy. Adding a duct—i.e., an additional surface or surfaces along which fluid flow must be considered and which is not moving with the turbine blades—adds a layer of complexity. While others seek to develop more accurate, and often more complex, models for turbines, the disclosed method takes advantage of the lower computational cost of a low-fidelity model (e.g., steady RANS-MRF) during the optimization process 110 and strategically only employs a higher-fidelity model (e.g., unsteady RANS-RSM) when necessary. In other words, the illustrated and described method intentionally uses a lower accuracy model to its own advantage. It has been determined that a lower accuracy model results in a highly accurate optimization, which those of ordinary skill in the art would not expect. Indeed this technique has resulted in ducted turbine geometries with efficiencies significantly greater than the 46% efficiency of a typical unducted turbine. The efficiencies of the optimized turbine designs in the examples disclosed herein are greater than 50%, greater than 52%, and / or in a range from 53% to 55%. The results also reveal new and unusual component geometries that provide these high efficiencies. For instance, the optimizations conducted in this disclosure effectively eliminate extraneous duct length and employ throat diameters that are a larger proportion of the exit diameter than expected. Duct lengths may be less than 3 times the effective radius R (from rotational axis to blade tip) of the turbine blade, or less than twice the radius R of the turbine blade. For such a low-profile (i.e., axially shallow) ducted turbine with a duct length less than the diameter of the rotating turbine to have such a high efficiency is an unexpected result. Unusual blade twist profiles have also been revealed as efficient, with a relatively high root pitch but a relatively low degree of twist along the high torque portion of the blades. Baseline geometries for the turbine blade 14* and the duct 16* are first identified in accordance with step 112 of FIG.9 as starting points for the optimization. Due to some possibility of the existence of local optima in the design space, there is no guarantee a priori that optimizations starting from different baseline designs will converge to the same ducted turbine geometry. For this reason, two different baseline geometries, referred to herein as baseline geometry A and baseline geometry B, are used in the examples below. Both baseline geometries use a baseline duct 16* having the same geometry, as depicted in FIG. 10. But the baseline turbine 14* for each baseline geometry has a different geometry, the twist-angle profiles of which are provided in FIG. 11. For the baseline duct 16*, the duct length is lB= 2.107 m, and the distance between the inlet end 18* and the throat T is lT = 0.557 m, such that the ratio ^^் / ^^^ൌ 26.4%. The inlet diameter is ^^^^ൌ 1.27 m, the throat diameter is ^^ଷ^ൌ 0.96 m, and the outlet diameter is ^^^௫^௧ൌ 1.536 m. For both baseline turbine blades 14*, ^^ ൌ 0.44 m, such that a radial gap G between the blade tips and the inner surface of the duct 16* at the throat T is ^^ ൌ 0.04 m. The baseline turbine 14* for baseline geometry A is the original Bahaj model (Bahaj, et al., “Power and thrust measurements of marine current turbines under various hydrodynamic flow conditions in a cavitation tunnel and a towing tank,” Renewable energy 32 (2007) 407-426) with a 0.44 m radius and a 20° root pitch at ^^ / ^^ ൌ 0.2 and the twist profile illustrated in FIG. 11. Blade geometry A has a 15° twist from base-to-tip. The turbine 14* for baseline geometry B is a modified Bahaj model with a 0.44 m radius, a 45° root pitch, and a modified twist profile as illustrated in FIG. 10, in which the blade has a 28° twist from base-to-tip. The particular twist profile of baseline geometry B is obtained by matching the local angle of attack of each blade section in the duct to that of the unducted turbine counterpart through an iterative procedure. Evaluated by the unsteady RANS solver, baseline Geometries A and B at λ = 5.5 respectively yield efficiencies of CP = 28% and CP = 45%. The turbine hub is not included in this example since it is difficult to make it compatible with the current geometry parametrization using the Free-Form Deformation method. The optimization setup, starting respectively from baseline Geometries A and B, is summarized in TABLE I. TABLE I Optimization A Optimization B Duct^ ^^^^ସ^ୀ^, ^^ (5 variables) r- Stokes equations, ∇ ⋅ ^^ ൌ 0, (4a) ^ ∂ ^^ 1 ^∇ ⋅ ^ ^^ ^^^ ൌ െ ∇ ^^ ^ ∇ ⋅ ^ ^^∇ ^ ^(4b) ∂^^ ^^^ ,where ^^ is the flow velocity and p is the pressure. Equation Set (4) needs to be solved with boundary conditions of inflow velocity ^^^at the far field and no-slip boundary conditions on the blade and duct surfaces. As in many engineering problems of high Reynolds number, direct numerical simulation of Equation Set (4) is prohibitive for modern computers. To address this problem, one method is to use the Reynolds-averaged Navier-Stokes (RANS) equations with grids only resolving the averaged components of the flow. To Equation Set (4), one can apply the Reynolds decomposition^^ ൌ ^ഥ^ ^ ^^ᇱ(and the same for pressure), with the ^ഥ^ denoting the averaged velocity (in a timewindow or by an ensemble) and ^^ᇱthe zero-mean turbulent fluctuation, leading to the unsteady RANS equation: ^∇ ⋅ ഥൌ 0ப^^^ഥ^^^^̅^ ^ ^(5)whereത^ത^തᇱത^ത^തᇱis the example, we use the ^^ െ ^^ Shear Stress Transport (SST) turbulence model together with the automatic near-wall treatment (see Appendix A for details on both). Rotation of the turbine blades can be handled in simulations by two blade-resolved approaches: a Multiple Reference Frames approach (MRF) and a rotating-sliding mesh approach (RSM). The former may be applied in generating a steady solution of the RANS equations with a different reference frame (i.e., a solution with time-derivative terms set to zero), while the latter may be used directly in the unsteady solution of Equation (5). These two approaches are associated with lower and higher fidelity and are used respectively in the optimization process 110 and the final reevaluation process 130 of FIG.9. The MRF approach employs an efficient model for modeling turbomachinery flow in which the computational mesh remains stationary and the rotational effect is handled through a rotational reference frame. In particular, with reference to FIG.12, the fluid domain is separated into two regions: a rotational region FR surrounding the turbine blades with a blade-fixed x’-y’-z’ reference frame, and the remaining stationary region FS with an inertial reference frame x-y-z. In both regions, fluid flow is considered steady with respect to the corresponding reference frame, so only steady RANS equations need to be solved. To be more specific, in the rotational region FR, the blades 14 are stationary and experience a steady inflow. The flow velocity in the blade-fixed reference frame can be expressed by ^^ோൌ ^^ െ ^^ ൈ ^^, (6) where U is the velocity in the inertial reference frame, ^^ is the rotation vector of the turbine blades, and r is the distance vector from the axis of rotation to the point of interest (position vector). The steady RANS equations in the rotational region also need to be established with the blade-fixed reference frame, which requires further formulations of both Equation (5) and the ^^ െ ^^ SST model equations. Although the OpenFOAM / DAFoam code solves this complete set of equations for this example problem, only the rotational-region formulation regarding Equation (4) is presented here to provide key insights of the method. Combining Equation (4) and Equation (6), the following is obtained (see Appendix B for detailed derivation): ∇ ⋅ ^^ ൌ 0, ^ப^^^^ ^^^^ ^^ ^ ^ ^^^ ^^ (7) The steady = 0. Therefore, in the MRF method, steady versions of Equations (4) and (7) are solved in stationary and rotating regions, with the additional procedure regarding Reynolds decomposition and steady k − ω SST equations. This can be done using the SIMPLE algorithm implemented as simpleFoam in OpenFOAM. Although the RANS-MRF method provides an efficient numerical solution for the turbine problem (i.e., only two steady RANS equations need to be solved), its accuracy can be compromised mainly because of two issues. First, the rotational region FR and stationary region FSare usually chosen in a subjective manner. There is no guarantee that the rotational region covers all of the flow features resulting from the rotating and discrete blades. Any mismatch between the choice of the region and the nature of the flow can lead to errors at the interface and thus in the final results. Second, for many designs, the flow can be unsteady in nature, especially when flow separation occurs from the duct and / or blade surfaces. Assuming a steady state, as in the RANS- MRF method, can lead to significant errors for this type of unsteady flow. Because of the above reasons, the MRF model is considered a low-fidelity model, and the RANS-MRF approach is considered a low-fidelity technique in the context of this disclosure. The rotating-sliding mesh (RSM) approach allows for direct simulation of the unsteady RANS (URANS) equation (5) with mesh domains that exhibit relative motions, which is required for modeling rotating geometries. The underlying idea of RSM model is to allow a region of computational mesh M to rotate with the turbine blades, as illustrated in FIG. 13. The rotational mesh creates a technical problem in that the mesh at the rotating / non-rotating interfaces I becomes non-conformal—i.e., the nodes on opposite sides of each interface I do not match up. The data transfer across the interface, therefore, may be handled by a special interpolation method involving a “supermesh,” as described in more detail in Appendix C. URANS coupled with the RSM model is, in principle, much more accurate than the RANS- MRF method as it resolves the unsteady nature of the fluid flow. This is a significant advantage in simulating flow around the ducted turbine 10 because the possible flow separation from the duct surface 22 can be captured much better. On the other hand, in the URANS-RSM approach, a small time step is needed to resolve the blade rotation, yet a long simulation time is required for the solution to reach a quasi-stationary state. The computational cost of the URANS-RSM approach is therefore much greater than that of the RANS-MRF approach. In this disclosure, the URANS- RSM approach is considered a high-fidelity technique and may be limited to use in the reevaluation process 120 of FIG.11 after the optimization process 110. The unstructured computational mesh can be generated using the OpenFOAM meshing tool snappyHexmesh. An overview of the mesh is illustrated in FIG.14A-14C, which include an axial view from the inlet end of the duct in FIG. 14A, a side view in FIG. 14B, and an enlarged view of a portion of FIG. 14B in FIG. 14C. The total box size of the computational domain is 10.4D × 10.4D × 23.7L, where ^^ ൌ^^4 / ^^^ ^^ ൌ 1.536 m is the maximum diameter Dexitof the duct, and L = 2.107 m is the length duct 16* for a total mesh size of 16 × 16 × 50 meters. In the base 3D-mesh M illustrated in FIG.14, each 3D unit is 0.4D × 0.4D × 0.3L, with y-z cells (unit faces) at 0.4D × 0.4D and with x-y and x-z cells at 0.4D × 0.3L. This mesh size, together with the boundary conditions shown in FIG.14 is sufficient to avoid a blockage effect. To model fluid flow near the turbine and immediately downstream of the outlet end of the turbine with higher accuracy, a refinement region RR with dimensions of 2D × 2D × 2.5L is defined around the ducted turbine. Level-4 refinement is applied within the refinement region RR—i.e., each cell of the base mesh is divided into (2)4cells in each direction. Prism layers are then added around the duct and blade surfaces, with the cells further continuously refined as the mesh approaches those surfaces. Two prism layers are used around the duct surfaces, and three prism layers are used around the blade surfaces with expansion ratios of 1.1. The use of the prism layers provides better resolution for boundary layers and is critical to obtain well-convergent results in the grid sensitivity study shown later. In this example, three mesh resolutions, M0, M1, and M2, are used with an increasing number of cells—i.e., further refinement—from M0 to M2. The coarsest mesh M0 is that just described and illustrated in FIG.14. The M0 mesh is used in the RANS-MRF in the optimization process and has 2-3 million cells (the exact number depends on turbine geometry and re-mesh procedure in optimization). For the M1 and M2 meshes, the full M0 mesh, including the base mesh M and refinement region RR, is uniformly refined in each direction by factors of about 1.3 and 1.6, respectively, leading to 4-5 million cells for M1 and 7-8 million cells for M2. All meshes M0, M1, and M2 are used in the URANS-RSM reevaluation process, including the grid sensitivity study. The constrained optimization problem (3) can be converted to an unconstrained optimization problem, with the constraints (3c) to (3f) implemented as terms with Lagrangian multipliers and the constraint (3g) handled through Free-Form Deformation (FFD) points, as discussed later. Sequential Quadratic Programming (SQP), implemented in SNOPT (Sparse Nonlinear OPTimizer), can be used to solve problem (3), where the most critical (and challenging) task is to obtain the gradient of the objective function ∇CP (or simply ∇P) with respect to all 21 design variables. Below, two components (see FIG.11) involved in the computation of gradient are discussed: (1) geometry optimization so that the turbine geometry morphs smoothly with changes of design variables; and (2) an adjoint method for gradient computation which requires only one CFD run to compute the vector ∇ ^^ ∈ ℝଶ^. The gradient information is then used in the SQP algorithm to obtain the next design points. The procedures are then iterated until either convergence or the prescribed number of iterations is reached. The geometry deformation to the next design point, as well as the gradient computation, requires a method to deform the computational mesh together with the geometry, which will also be discussed below. To control the smooth morphing of the full geometry of the ducted turbine with the variation of design variables in the constrained optimization problem (3), the geometry is parameterized. In this example, Free-Form Deformation (FFD) is used for the geometry parametrization, implemented in PyGeo package of DAFoam. The principle of the FFD method is to enclose the geometry of interest in an FFD box with a specified number of control points (also known as FFD points). The FFD points are analytically connected to the enclosed geometry surface using tri-variate B-spline, with details presented in Appendix D. Through the control of FFD points, both global and local smooth deformation of the enclosed geometry are enabled. FIGS. 15A and 15B illustrate a two-dimensional example of geometry deformation controlled by FFD, where the FFD box FB, and control points CP are constructed to enclose the surface S to be deformed (FIG. 15A). The deformed surface S’ is illustrated in FIG. 15B. FIGS. 16A and 16B illustrate a three-dimensional example of geometry deformation controlled by FFD, where the FFD box FB, and control points CP are constructed to enclose the surface S to be deformed (FIG.16A), and the deformed surface S’ is illustrated in FIG.16B. FIGS. 17A-E illustrate an overview of an FFD setup for a ducted turbine 10. In this example, two levels of FFD boxes are used, with one cube-shaped parent box BP enclosing the entire ducted turbine geometry and two child boxes BD, BB. Child box BDis an octagonal tube shape and encloses the duct 16 of the ducted turbine 10. Child box BTis a truncated wedge shape having a quadrilateral cross-section that decreases in size from one lengthwise end to the other and encloses a single turbine blade 14. The above-mentioned 21 design variables can be represented by 21 degrees of freedom (DoF) associated with the FFD control points. The child FFD box BBfor the blades 14 is associated with 32 FFD control points placed on eight sections (i = 1 to 8) of one blade 14. For each four FFD control points at one section, two DoFs are provided. The first DoF is a linked (or simultaneous) rotation of the four FFD control points of the section, which controls the twist angle θi of the section, where the center of rotation is a quarter chord. The second DoF is a linked scaling (expansion or contraction) of the four FFD control points of the section, controlling overall area or scale biof the section. The FFD control points on all blades 14 are also linked to ensure the same deformation for all blades. The child FFD box BD for the duct 16 contains 112 FFD control points placed on 7 sections. Only 1 DoF is provided to the FFD control points controlling the duct length l. During the process of scaling (lengthening or shortening) the duct 16, all duct FFD control points move in the axial direction of the duct 16 with perturbations proportional to their axial distances from the throat T. This movement is to ensure that the throat T is consistently located with respect to the inlet end of the duct as a percentage of the overall duct length (26.4% in the present example). It should be noted that seven sections are not necessarily needed for the duct FFD box BD, as this choice is mainly for convenience during setup. The parent FFD box BPhandles the duct diameter changes, the constraint (3g) on the tip- gap ratio, and the condition ^^exitൌ ^^4 / ^^^ ^^. Twenty-eight FFD control points are placed on seven sections in the parent box, with the control points for the three sections nearest the outlet end of the duct 16 packed closely together in the axial direction (see FIG. 17B). These control points remain stationary in optimization, such that ^^exitൌ^^4 / ^^^ ^^ is guaranteed. The 16 FFD control points on the four sections of the parent box BPnearest the inlet end of the duct are used to control the duct diameters—i.e., design variables dj. As these FFD control points move radially, the entire geometry—including the child FFD boxes BD and BB—deforms according to the analytical link of boxes BD and BB to the control points of the parent box BP, since the children boxes BDand BBare all enclosed in the parent box BPby definition. As such, the constraint (3g) is automatically satisfied since the blade 14 expands or contracts proportionally to the duct throat T. Such an implementation of FFD for the complex geometry of a ducted turbine is previously unknown in the art. To compute the derivative of an objective function—in this case, the power coefficient CP—with respect to design variables, an adjoint method is used. In order to clearly explain the adjoint method, we first introduce notations as follows: Let ^^ ≡ ^^ ^^^, ^^^^଼^ୀ^, ^ ^^^^ସ, ^^^ ∈ ℝே^with ^^௫ൌ 21 as the design of the RANS-MRF equation. Here ^^ ∼ ^^^10^^ includes three velocities and pressure at each cell in the computational grid. The goal is to compute ^^ ^^^ / ^^ ^^ ∈ ℝଶ^. If a finite-difference method is used to compute the derivative, one needs at least 22 CFD simulations for each derivative computation, even with the lowest order approximation. This is computationally prohibitive for the problem at hand. For the adjoint method to compute ^^ ^^^ / ^^ ^^, the first step is to write the function ^^^^ ^^, ^^^ ^^^^, and express its total derivative with respect to x as ௗ^ು ப^ು ப^ௗ ^^ ดௗ ^^ൌดப ^^^ดುப^^ ௗ^^^, (8) where ^^ ^^^ / ^^ ^^ should be (i.e., geometry) are varied, with flow solution s remaining unchanged. ^^ ^^^ / ^^ ^^ is the change of CP as the flow solution s changes with a fixed turbine geometry. These partial derivatives are relatively easy to compute, with more details presented in Appendix E. The term that is difficult to compute in Equation (8) is ^^ ^^ / ^^ ^^. To compute it, one needs to further involve the RANS-MRF state equations in terms of their discretized residual form ^^^ ^^, ^^^ ^^^^ ൌ 0. Here, ^^^ ^^, ^^^ ^^^^ ∈ ℝெconsidering the same number of equations as the number of unknowns in s. Since ^^^ ^^, ^^^ ^^^^ should remain zero with a change of x (if flow solution is correctly obtained), then ௗ^^ ப ^^ௗ ^^ப ^^ௗ^ ^^ൌ 0 ⇒ப^ ^^ ௗ^^^ൌ െப^ ^^. (9) ெൈே^ெൈெ ெൈே^ெൈே^Direct solution of Equation (9) gives ௗ^^ି^ ൌ െ ^ப ^^ ப ^^ ௗ^^ᇣபᇤ^^^ᇥப^^^. (10) ெ ெൈே^ It is worthwhile discussing the computational cost associated with Equation (10) at this point. The matrix multiplication in Equation (10) leads to a computational complexity of ^^^ ^^ଶ^^௫^, which is very expensive since ^^ ∼ ^^^10^^ and ^^௫are also large. This has to be added by the cost to invert a ^^ ൈ ^^ matrix, which is, in general, more expensive. Even if one uses some iterative solver for linear systems to solve Equation (9), the procedure needs to be repeated for ^^௫times since ^^ ^^ / ^^ ^^ (as well as the RHS) has ^^௫columns. The computation is therefore also very expensive. On the other hand, the computational cost can be significantly reduced by substituting Equation (10) to Equation (8) and considering a re-grouping of the multiplications:ௗ^ดುௗ^^ൌ ப^ดುப^^െ ^ ப^ดುப^^^ப ^^ି^ ᇣபᇤ^^^ப ^^ᇥ^ப^ ^^. (11) Instead of computing is first computed. This computation can be performed by solving the so-called adjoint equation, with the adjoint being equivalent to the transpose of a real matrix in this case: ப^் ் ^^ᇣபᇤ^^^ᇥ^^^ൌ ^ப^ುᇣபᇤ^^ᇥ ^ (12) whose solution transpose provides ப ப ^^ ି^ ^^்ൌ^ು^^^^^^ (13) as the parenthesis term in Equation The solution of Equation (12) involves solving a linear system only once, instead of ^^௫times as needed for Equation (9), and hence is much less expensive. It is also less expensive when compared to the direct computation of Equation (10). The computational cost to solve Equation (12) is generally similar to the RANS-MRF computation. Therefore, in each iteration of the optimization, the computational cost is on the same order as one RANS-MRF solution. The only remaining component is the calculation of partial derivatives ∂R / ∂s in Equation (12), which can be found in Appendix E with other derivatives mentioned above. In both the partial derivative computation in the adjoint method (Appendix E) and the matching to the next design point, turbine geometry is deformed, which leads to the deformation of the entire computational volume mesh. The volume mesh deformation is computed based on the analytic inverse distance weighting method, implemented in the IDWarp package of DAFoam. The principle of the method can be understood as follows. Given a 2D surface (e.g., a blade surface) with a surface mesh of N nodes, the surface geometry deformation results in the movement of each node. Two quantities (Mi, bi) for each node are assigned with i = 1, 2, ...,N, where bi is the translation distance of the node and Mi is the rotation matrix such that ^^new^ ^^old new old^ , with ^^^and ^^^the normal vectors at the node. In particular, both ^^new^ and ^^o^ldare computed by a weighted average of the normal vectors for all surrounding of the node. After (Mi, bi) are obtained for i = 1, 2, ..., N, the deformation of any volume mesh can be computed by summing the contribution from each surface node—i.e., Δ ^^ ൌ ∑ே^ୀ^  ^^^^ ^^^^^ ^ ^^^െ ^^^ with ^^ any volume node and Δ ^^ its movement. The weighting factor ^^^has the empirical form that grows in a polynomial form with the inverse distance between and surface nodes. FIG.18 illustrates the undeformed computational mesh associated with a portion of the ducted turbine with a baseline geometry (see FIG.14C), and FIG.19 illustrates an example of the computational mesh of FIG.18 after deformation associated with radial expansion and axial shortening of the duct. Before illustrating the results of the optimization and reevaluation processes, two additional studies are described. The first is experimental validations of the RANS-MRF and the URANS-RSM solvers, and the second is about a grid-search study with five design parameters using low-fidelity and coarse-grid (even coarser than M0) flow solvers conducted before the gradient-based optimization. Since there are no known systematic measurements of the performance of ducted hydrokinetic turbines, the published experimental results of the unducted Bahaj turbine are used for validation. FIG. 20 illustrates the power coefficient CPfor the unducted Bahaj turbine for a range of λ values, obtained from the RANS-MRF and the URANS-RSM, in comparison with Bahaj’s experimental results. The RANS-MRF solver is run with a 2nd-order numerical scheme for the convection term in the RANS equations, which differs from the optimization process as described later in detail. It is clear from FIG.20 that both solvers predict similar values of CP for the Bahaj unducted turbine, with both predictions relatively close to the actual results in the chosen range of λ. As demonstrated later, RANS-MRF and URANS-RSM do not predict similar values for ducted turbines. The relatively small mismatch between the simulations and experiments at higher values of λ in FIG. 20 may be related to the fact that the simulations are conducted at a slightly higher Reynolds number (2×106vs.1×106) and that no hub is included in the simulations. However, the influences of those two factors are rather insignificant considering the favorable comparison in FIG.20. The grid-search study is detailed in Appendix F and does not provide a successful design of the ducted hydrokinetic turbine with improved efficiency relative to the unducted Bahaj model. This failure may be due in part to the low-dimensional parameter space, which is insufficient to explore effective designs. It also implies that designs with more parameters conducted using a gradient-based method as disclosed herein are necessary for designing complex geometries such as ducted turbines. The constrained optimization problem (3) can be solved using the above-described methods and processes. For this example, the processes were implemented in DAFoam. In the RANS-MRF solver, a first-order numerical scheme is used for computing the convective term— i.e., to construct the flux in cell faces. The first-order scheme is more dissipative than the normally used second-order scheme, but the former is important to obtain convergent flow solutions for many duct designs—especially those associated with flow separation. Specifically, after extensive tests, it was discovered that the second-order scheme shows more fluctuating flow solutions in many cases, which in turn affects the accuracy of the adjoint method, thereby prohibiting an accurate gradient computation. On the other hand, while the first-order scheme may provide even less accurate solutions for cases with flow separation, the obtained CP is usually low for the converged flow solution and the optimization can drive the design to new ones with improved efficiency. As such, the use of a lower order numerical scheme for computing the convective term of the RANS equation is surprisingly beneficial, given its inherently lower accuracy. In each RANS-MRF simulation, the solution is considered to be convergent as the residuals stop dropping, which in general provides residuals of momentum equation to be ^^^10ିସ∼ 10ିହ^. FIG. 21 shows the change in CP for Optimization A, starting with baseline geometry A, and Optimization B, starting with baseline geometry B, as the optimization process progresses through multiple iterations. The two starting points yield CP= 25% and 41%, respectively, as evaluated by the RANS-MRF. Both of these values are lower than other reports (28% and 45%) based on URANS-RSM. From there, the optimization process of both baseline geometries leads to a fast increase in CPat the beginning until both CPvalues effectively plateau near CP= 48%. The small dips in both curves in FIG. 21 between 20 and 30 iterations correspond to a restarting / remeshing procedure, which is necessary because mesh deformation (induced by geometry deformation) can deteriorate the quality of the flow solution and adjoint solver. A manually implemented restarting / remeshing procedure can improve the situation, leading to additional increases in CP. Both Optimizations A and B are halted when the CP value plateaus even with further restarting / remeshing. In practice, it was found that this convergent situation corresponds to the SNOPT’s optimality metric of approximately 10-1.7, which is consistent with cases in an unducted wind turbine optimization. For such optimality conditions, although ∇P does not completely vanish, the benefit of further optimizing the ducted turbine geometry is compromised by the mesh deformation so that some practical optimal points are reached. The two optimized geometries yield similar values for CP, with Optimized Geometry A at CP = 0.4782 and Optimized Geometry B at CP= 0.4822, achieved respectively at λ = 6.18 and λ = 6.39 (since the blade radius R is optimized, which affects λ). Since RANS-MRF is the low-fidelity solver, for reasons discussed above (including the first-order convection scheme), Optimized Geometries A and B undergo the higher-fidelity efficiency reevaluation process using URANS-RSM. The values obtained from URANS-RSM on the M0 grid are added to FIG.21, with Reevaluation A at CP = 54.6% for Optimized Geometry A and Reevaluation B at CP= 52.9% for Optimized Geometry B. Here, the difference between the two solvers as applied to ducted turbines, due to the complexity of the flow, is remarkably larger than that for unducted turbines such as that illustrated in FIG. 20 where the difference between RANS and URANS is less than 1.5% over the full range of λ. A grid sensitivity study was also conducted, which evaluates CPfor the Optimized Geometries A and B using URANS-RSM on meshes M0, M1, and M2 with increasing resolution. In these URANS solutions, an adaptive time step is used, with an average time step of 1 × 10−5seconds—i.e., about 0.01° turbine rotation for 1 time step—with total simulation times of 15 seconds to reach the quasi-stationary state. With all simulations conducted on the NSF Stampede2 cluster, the simulations on the M2 grid take about 240 hours on 576 CPUs. The results are illustrated in TABLE II, which indicates that CPvalues vary by only O(1%) in terms of the absolute value—i.e., not very sensitive to the large range of variation of grid resolutions. Based on the results from M2, the two final designs yield a similar CP ≈ 54% that is much higher than 46% of conventional unducted turbines, such as the Bahaj model. TABLE II Reevaluated Optimized Geometry A Reevaluated Optimized Geometry B 7 5 7 The efficiency CPof the two optimized designs is further evaluated at a range of values of λ using URANS-RSM on the M0 grid. The results are illustrated in FIG. 22 together with the efficiencies CP of the unducted Bahaj model (from FIG. 20) and an unducted turbine optimized with the same setup. The unducted turbine is optimized for fixed Ω = 21 rad / s, corresponding to λ = 6. Notably, the optimized ducted turbine geometries not only work well for the designed value of Ω, but also perform with high efficiency for a large range of λ. Moreover, the maximum CP for each geometry is in fact not achieved at the designed tip-speed ratio λ (marked by stars in the figure), but at some larger value of λ. Considering FIG.22, the maximum CPfor the two optimized geometries A, B are 56% and 54%, achieved at λ = 6.94 and 6.99, respectively. The optimized geometries A and B are illustrated in FIG.23 together with the geometry of the common baseline duct 16*. In general, the optimized geometries have a duct length that is shorter than that of the baseline duct 16* and are radially expanded relative to the baseline duct. In particular, for the optimized geometries A and B, the respective duct lengths lA and lB are 0.857 m and 1.327 m, and the respective radii RAand RBat the throat are 0.559 m and 0.554 m. This represents respective 59.3% and 37.0% reductions in length relative to the baseline duct, and respective 16.5% and 15.4% increases in throat radius relative to the baseline duct 16*. The twist and chord length (i.e., the distance profiles of the optimized and baseline geometries) are illustrated in FIG.24. The chord lengths of the optimized blade geometries do not vary significantly from the baseline geometries, but the optimization led to significant changes in the twist profiles of the turbine blades 14. The blade of optimized geometry A has a root pitch of 27° and a 14° base-to-tip twist, and the blade of optimized geometry B has a root pitch of 40° and a 28° base-to-tip twist. The degree of twist thus did not change significantly from the respective baseline geometries. Optimized geometry B is essentially the same twist profile as baseline geometry B with the root pitch shifted downward by 5°. Optimized geometry A has little or no twist for r / R < 0.3 but converges with optimized geometry B for r / R > 0.4. Notably, although the twist profiles of the two baseline geometries are drastically different, the corresponding optimized geometries converge to very similar twist distributions, particularly for r / R > 0.35, where most torque is expected to be generated. The two different optimized blade geometries A and B may indicate that the optimization finds two different local optima in the design parameter space. Given the comparable performance despite the different duct lengths, it may be inferred that the blade twist profile and duct throat size are the driving design parameters. The flow fields of unducted turbines, the baseline ducted turbines, and the optimized ducted turbines are analyzed in order to understand the major flow mechanism leading to the improvement in efficiency. Optimized geometry A is used as an illustration. FIG. 25 illustrates variations of some relevant performance metrics together with the variation in CPin the first 22 iterations of the optimization process—before the first restarting / remeshing apparent in FIG. 21. The metrics in FIG. 25 include the flowrate J passing the turbine blades and the thrust coefficient ^^்ൌ ^^ / ^0.5 ^^ ^^^ଶ^^^ on the duct and blades. In FIG.25, the 22 iterations of the optimization process are also divided into three stages. Stage I is the fast growth stage of CP, stage II is the slow growth stage of CP, and stage III is the plateau stage of CP. In stage I, while the thrust coefficient CTremains almost unchanged for the duct, the blade CTincreases rapidly. This is the most favorable situation to improve CPsince clearly more and more loadings from the total are distributed on the blades. In stage II, this favorable variation of CTcannot be maintained—i.e., its potential has been exhausted in stage I, and the opposite trend is observed with increased duct CTand decreased blade CT. The further slow increase of CPin stage II therefore must be associated with a different mechanism that is perhaps the more effective transition from blade loading to rotational motion (or torque). Both duct and blade CT remain unchanged in stage III as CPplateaus. The overall increase in CTin the whole process is only 14%, in contrast to the 87% increase in CP. This is another favorable feature of the current optimization since a much higher CT would be much more demanding on supporting structures. Finally, the flowrate J is highly correlated with CPover the whole optimization process and increases constantly until its simultaneous plateau with CP. The above analysis prompted further study of the flowrate metric, which remains consistent with the trend of CPduring in the optimization process. FIGS.26A-C include flow visualizations of streamlines in a flow tube passing the blades of the optimized unducted turbine (FIG.26A), the baseline ducted turbine geometry A (FIG. 26B), and the optimized ducted turbine geometry A (FIG. 26C), obtained in the quasi-steady solution from the URANS solver. To facilitate a meaningful comparison, only the streamlines in the flow tube that pass the turbine blades are shown. Since the inflow velocity is fixed at 1.4 m / s for all three cases, the flowrate in each illustrated flow tube is proportional to the inlet area AI of the respective tube. Here, each inlet area AIcorresponds to the cross-sectional area taken upstream of the turbine that includes all of the streamlines that pass the turbine blades, which is different from the area of the inlet end of the duct. From FIG.26, it is clear that the optimized geometry in FIG.26C corresponds to the case with the largest tube inlet area AI. That is, fluid from a larger area upstream of the turbine makes its way through the turbine of FIG.26 than through the unducted and baseline geometry turbines. Physically, this indicates that a well-designed geometry draws a larger volume of water (relative to unducted and baseline geometries) into the throat, which is accompanied by a higher flowrate across the blades. The volumetric flowrate across the blades of the optimized ducted turbine is 2.25 m3 / s, compared to 1.80 m3 / s for the baseline ducted turbine and 1.55 m3 / s for unducted turbine. The flowrate is a more effective indicator of ducted turbine performance than flow velocity at the throat. As an example, the baseline design A is associated with an accelerated flow speed at the throat, but no improved efficiency. This analysis explains for the first time that improved flow conditioning provided by the duct of a ducted turbine enhances its performance. In summary, gradient-based geometry optimization can be conducted on ducted turbines, such as hydrokinetic turbines, using computational fluid dynamics and the adjoint method. This optimization was conducted on ducted turbine baseline geometries having drastically different turbine blade geometries with different respective performances. The resultant geometries for both cases yield similar performance with CP≈ 54% as evaluated by the high-fidelity URANS solver, and both geometries capture similar critical geometrical features in terms of the duct throat area and blade twist profile. This value of CP is significantly higher than standard unducted turbines such as the Bahaj model. While the optimization may be performed at a single rotational speed, the optimized geometries also yield high performance over a wide range of rotational speeds. By evaluating the flow mechanism associated with the performance improvements, it has been discovered that the CPamong different geometries is correlated to the flowrate of fluid passing the turbine blades. The optimized geometry corresponds to the case with a maximum flowrate due to the suction of a well-designed duct. With reference to FIG.27, the full methodology can be extended to include other turbine features in the optimization, including a turbine hub 36 of sufficient size to house a generator, a larger design space of the duct to exhibit arbitrary thickness profiles (e.g., with an air-foil shape), and realistic constraints such as a minimum thickness of the duct and its maximum curvature at the leading edge. To accommodate these additional features, the above-described FFD parameterization may be replaced by ESP (Engineering Sketch Pad) parameterization. ESP allows parameterization of the geometry with flexible geometric features (e.g., coordinates, polynomial coefficients, twist and chord length of blade sections) and is thus suitable for the more complex geometry associated with the addition of other turbine features. This optimization illustrated in FIG. 27 is performed for a turbine with rated power of 5KW, which requires the cylindrical part of the hub to have a radius of at least 0.2 meters and a length of at least 0.78 meters. Baseline geometries 16*, 36* and optimized geometries of the duct 16 and hub 36 obtained using the above described methods are illustrated in FIG.27. The baseline geometry 16* of the cross-section of the duct is a foil shape, and the baseline geometry 36* of the turbine hub has a hemispherical nose at one end of a cylindrical body and a trailing end that tapers to a generally conical shape at the opposite end of the body. While only the inner surface of the duct 16 is considered in the above examples, the outer surface of the duct is additionally considered in the example of FIG.27 such that the duct has a thickness equal to the difference of the radius at the outer surface and the radius at the inner surface, which can both vary in the optimization problem with a constraint set for the minimum required thickness and for the curvature at the inlet end 18. For the hub 36, the cylindrical portion of the baseline geometry 36* is held constant while the lengths and diameters of the nose and trailing end are permitted to vary with some constraints, such as the respective diameters being equal to that of the cylindrical portion where they meet the cylindrical portion. The optimized geometries include certain unconventional features. First, the nose of the hub 36 protrudes beyond inlet end 18 of the duct 16, resulting in a stagnation point farther upstream than is conventional and thus more acceleration of the fluid flow to the turbine blades. Second, the cross-section of the duct 16 is deformed from the baseline foil shape to a thin-wall shape constrained by the finite leading curvature and minimum thickness for manufacturing considerations. This allows for a more cambered duct surface, given the fixed total device projected area, to generate more acceleration of the fluid flow. The optimized turbine yields an efficiency of about 50%, evaluated by both RANS and URANS CFD simulations. Appendix A. k − ω SST turbulence model and automatic wall treatment The approximation of the Reynolds stress term in Equation (5) can be obtained by the Boussinesq hypothesis, which assumes that the Reynolds stress is related to an eddy-viscosity ^^௧and the mean velocity gradients, namely (for incompressible flow)െത^ത^തᇱത^ത^തᇱൌ ^^ ൬ப^ഥ^^ ப^ഥೕଶ పఫ ௧ ப௫ೕ ப௫^^ െଷ^^ ^^^^(A.1)where ^^ ൌ^ത^ത^തᇱത^ത^തపᇱis the ^^ and ^^௧, two-equation turbulence models based on ^^ െ ^^ model or ^^ െ ^^ model have been developed. Here, ^^ is the turbulence dissipation rate, and ^^ ൌ ^^ / ൫ ^^ ^^ఓ൯ with ^^ఓ= 0.09 is the so- called specific turbulence dissipation rate. These models include transport for ^^ and ^^ (or ^^), and ^^ is obtained as a function of the computed ^^ and ^^ (or ^^) It has been identified that the ^^ െ ^^ model is robust for regions far from the wall, but not accurate when integrating down to the wall; and some wall functions are necessary to simulate the viscous sublayer. On the other hand, the ^^ െ ^^ model offers a more accurate resolution of the viscous sublayer but is less accurate for the far field—e.g., it is over-sensitive to the freestream turbulent conditions. On top of both models, the ^^ െ ^^ based Shear Stress Transport (SST) model is proposed, which connects the near-wall region predicted by the ^^ െ ^^ model and far-field by the ^^ െ ^^ model using some blending functions. In principle, the ^^ െ ^^ SST model captures the advantages of both models. In addition to the ^^ െ ^^ SST model, an automatic wall treatment is applied for near-wall simulation. The automatic wall treatment offers a wall ^^ା-insensitive simulation—i.e., great flexibility in the first-cell size at the wall. The principle is to set up different wall boundary conditions for ^^ and ^^ depending on the ^^ାvalue. It has been verified that the flow prediction in the near-wall region is insensitive and robust for ^^ାfrom O(0.1) to O(100), and the automatic wall treatment has been practically applied for cases with the first cell up to the range in the logarithmic layer— ^^ା^ 300. Appendix B. Derivation of equations in the rotating reference frame Transformations of velocity and acceleration are first stated in inertial and rotating reference frames ^^ ൌ ^^ோ^ ^^ ൈ ^^ (B.1) ^^ ^^ ^^ ^^ ^^ ^^ where variables with the case ^^ ^^ / ^^ ^^ ൌ 0 that is consistent with the instant application. Derivations of Eq. (B.2) can be found in standard textbooks on dynamics. Substituting Eq. (B.1) into Eq. (4a), and making use of the fact that divergence of a curl of a vector is zero, the continuity equation is obtained in rotating reference frame ∇ ⋅ ^^ோൌ 0. In order to derive momentum equations in the rotating reference frame, consider Eq. (4b) with material derivatives replacing the unsteady and convection terms, to which Eq. (B.2) can then substituted to obtain ^ௗ ^^^ ^ 2 ^^ ൈ ^^ ^ ^ ௗ௧ ோ ோ^^ ൈ ^^ ൈ ^^ ൌ െఘ∇ ^^ ^ ∇ ⋅ ^ ^^∇ ^^^ (B.3) The viscous ∇ ⋅ ^ ^^∇^ ^^ ൈ ^^^^ ൌ 0. The material derivative in Eq. (B.3) is then expanded to obtain the momentum equations in the rotating reference frame: ப^^ೃௗ௧ ^ ∇ ⋅^^^ோ ^^ோ^ൌ െ^ ఘ∇ ^^ ^ ∇ ⋅^^^∇ ^^ோ^െ 2 ^^ ൈ ^^ோ െ ^^ ൈ ^^ ൈ ^^. (B.4) be simplified further for easier implementation in the finite volume method. For this purpose, the convection term is rearranged as ∇ ⋅ ^ ^^ோ^^ோ^ ൌ ∇ ⋅ ^ ^^ோ^ ^^ െ ^^ ൈ ^^^^ ∇⋅^^^^∇ ⋅ ^^ ^^^ ⋅ ∇^ ^^ ^^^where the ^^ோ⋅ ∇^ ^^ ൈ ^^^ ൌ ^^போ^ ப௫^^^^^^Ω^^^^The final equations in the ∇⋅ ^^ ൌ 0^ which is Eq. (7). For , are consistent with the outer stationary region, and one unified code is used for the solution with some minor modifications—in terms of forcing ^^ ൈ ^^, and the ^^ோterm leading to a correction in computing the cell face flux—for the rotating region. Appendix C. Data Transfer at the Interface in the RSM The interfaces involved in the turbine problem are illustrated in FIG.13. For cells on each side of the interface I, the flux value needs to be constructed on the cell face toward the interface. This cannot be obtained from a standard interpolation method because the cells on two sides of the interface are non-conformal—i.e., with non-overlapping cell faces from the two opposite sides. Given the donor and target meshes at the interface, the goal is to compute the flux value at the cell face of the target mesh. This requires reconstruction of the cells, as well as cell-centered properties, on the donor side into a “supermesh” whose faces at the interface contain all nodes of the target mesh, so that the cell face flux on the target side can be computed by interpolation between the supermesh and target mesh. In general, the interpolation via supermesh involves a Galerkin projection method. A simple example is depicted in FIGS. 28A-28C to illustrate the gist of the approach. With donor mesh TD and target mesh TT that are non-conformal as in FIG. 28A, the goal is to construct a supermesh TS that serves as a common ground for interpolation. The supermesh TS has to satisfy the following properties: 1. TScontains all the nodes on TDand TTand intersections of their edges; and 2. The intersection of every face cell element of TS with any element of TD or TT must either be zero or the whole element. The connection between TDand TSis depicted in FIG. 28B, with the shaded areas denoting cell faces of TD. Since each cell face of TD is now perfectly split into triangle cell faces of TS, properties can be assigned to the center of the volume supermesh (associated with triangle faces) according to the area ratio of the triangles. In other words, the property at TDis distributed into fractions according to the area fraction of triangles in TS. The properties at TS can then be used for interpolation with that at the target cell TT since each target cell face is also perfectly split into triangle faces of TS, as in FIG.28C. Specifically, the flux value associated with the each triangle area is obtained by interpolating between the volume cell of TSand the corresponding volume cell at the TT side. These fluxes are then added to form the value for a single cell face of TT . Appendix D. Analytical connection between FFD points and enclosed geometry While the FFD method is versatile in handling geometries in different dimensions, the following focuses on the situation consistent with the above disclosure, which involves the FFD control points creating a lattice box in ℝଷthat encloses a three-dimensional geometry within the box. The analytical connection is to establish a mapping ℝଷ→ ℝଷfor each point inside the box. The mapping can be established by making use of the tri-variate B-spline function: ^^^ ^^, ^^, ^^^ ൌ∑ேೠି^ ேೡି^ ே^ି^^ୀ^  ∑^ୀ^  ∑^ୀ^   ^^^,^ೠ^ ^^^ ^^^,^ೡ^ ^^^ ^^^,^^^ ^^^ ^^^,^,^ (D.1)where ^^^ ^^, the geometry of interest is a part of this volume, parameterized by u, v, and w (all ∈ ^0,1^). The indexes i, j, and k loop in three directions of the lattice, and ^^^,^,^∈ ℝଷare coordinates of the lattice FFD control points. ^^௨, ^^௩, and ^^௪are prescribed degrees of the B-spline basis functions in three directions. The B-spline basis function is defined recursively by, taking the i-direction as an example, ^^^,^^ ^^^ ൌ ^1 if ^^^ ^ ^^ ^ ^^^ା^0otherwise, (D.2)where ^^^is the so- (D.1), the sensitivity of the geometry to a particular FFD control point, in terms of derivative, can be analytically expressed as ப^^௨,௩,௪^ൌ ^^^,^ೠ^ ^^^ ^^^,^ೡ^ ^^^ ^^^,^^^ ^^^. (D.4) Since ^^ ∈ ℝଷand ^^ ∈ ℝଷ, It essentially means that the three piecewise derivatives of one ℝଷvector with respect to another ℝଷvector are equal with one another and also equal to the RHS. In our case, linked FFD control points are used as one degree of freedom. The derivative of the geometry coordinates with respect to the particular DoF motion can be computed by summing the contribution of each FFD point using Eq. (D.4). Appendix E. Partial Derivative Computation in the Adjoint Method This section presents the underlying principle of how the partial derivatives in the adjoint method can be computed. All these derivatives are computed by backward automatic differentiation in DAFoam with graph coloring method for acceleration. Four types of partial derivatives are involved: ப^ดು ப^ು ப ^^ ப ^^ப^^,ดப ^^,ப^ ^^,ப^ ^^(E.1) ^ൈெ ^ൈே^ெൈெ ெൈே^For the former two, consider ^^ ൌ Ω^ ^^ ൈ ^^ ^^ ^^, where ^^ is the pressure, ^^ is the distance to the rotating axis, and the integration is over the blade surface with differential element ^^ ^^. It is clear that ^^^depends on both the state variable ^^ (i.e., pressure ^^) and design variable ^^ (i.e., the integration surface). The derivative ∂ ^^^⁄ ∂ ^^ can be simply computed by perturbing the pressure field in the flow solution. For computing ∂ ^^^⁄ ∂ ^^ , we write ^^ ൌ Ω ∑^  ^^^ൈ ^^^Δ ^^^, which is in discrete form. A perturbation in ^^ results in perturbations of the points, which deforms the geometry. This, in turn, leads to the surface mesh deformation the ^^^and Δ ^^^. The partial derivative ∂ ^^^⁄∂ ^^ can therefore be computed using the chain rule to all processes. For the latter two, we consider the discretized governing equation in residual form, i.e. ^^ ൌ 0. Perturbation of ^^ changes the velocity field, pressure field, and the velocity flux at cell faces so that ^^ is perturbed. The derivative ∂ ^^⁄ ∂ ^^ can then be computed correspondingly. Perturbation of ^^ leads to the deformation of all computational mesh in the full fluid domain. The construction of velocity flux or, more precisely, the coefficients in the algebraic governing equation is thus affected, leading to the computation of ∂ ^^⁄∂ ^^ correspondingly. Appendix F. Discrete grid-search in a lower-dimensional design space This section describes an effort in a “brute force” grid search of a lower dimensional design space. Taking the baseline duct geometry as a reference, five sets of duct geometries were first created, each set representing one type of duct with five varying geometries. The five sets are as follows: (1) varying the outlet area from the baseline geometry (FIG.29A); (2) varying the outlet area from the baseline geometry with the throat T as the duct inlet (FIG.29B); (3) varying both the inlet and outlet areas from the baseline geometry with a finite- thickness duct, where the outer surface of the duct is a straight surface between duct ends (FIG.29C); (4) varying the outlet area from the baseline geometry with a flange at the outlet end (FIG 29D); and (5) varying the inlet area from the baseline geometry (FIG.29E). These represent 25 individual duct geometries in total, with the ratio of maximum and minimum duct areas ranging from 1.25 to 2.75. For each of the duct designs, other design variables considered are listed below with discrete values: • λ: 3, 3.5, 4, 4.5, 5, 5.5, 6, 6.5, 7, 7.5 • Blade root pitch: 10°, 20°, 30°, 40°, 45° • Twist profiles: Baseline blade geometries A and B • Tip-gap ratio: 0.09R, 0.2R The RANS-MRF solver is used for each design to evaluate its performance on a coarse grid with 1-2 million cells (even coarser than M0). In total, about 450 cases are run with results of their efficiencies Cpillustrated in FIG. 30. Since it is difficult to visualize results in the design space of more than three dimensions, all results are simply plotted as functions of λ—i.e., for each λ, a large number of Cp values resulting from varying other design variables are stacked. While the goal is not to further distinguish designs at each λ, it is clear that the maximum Cpamong all O(450) cases is only 38%, as evaluated by the RANS-MRF. A few geometries having relatively high Cp were selected for the reevaluation process using the higher-fidelity URANS solver on the M0 grid, with those results also included in FIGS.29A-E. The maximum efficiency computed by the URANS is 45%, which happens to correspond to the baseline geometry B in terms of both the duct and the blade geometry. This is still lower than the 46% efficiency of the unducted Bahaj model. It can be concluded that a “brute force” grid search does not provide any ducted turbine geometries with higher efficiency than that of standard unducted turbines. In other words, arriving at a ducted turbine geometry that optimizes turbine efficiency is demonstrably not a matter of routine experimentation involving the tuning of a result-effective variable. Indeed, no result-effective variables have been identified in the prior art that have a predictable effect on ducted turbine efficiency, which is due in part to the interdependence of multiple variables on each other. The methods disclosed herein address this by providing a geometry optimizing technique that iteratively changes multiple variables simultaneously to approach and identify an optimized geometry. In various embodiments, a system for optimizing the geometry of a ducted or unducted turbine includes at least one processor and memory (implemented as one or more non-transitory computer-readable mediums) storing or having instructions that, when executed by the at least one processor, cause the system to perform one or more of the above-described method steps and / or one or more method steps based on or derivable from the above-discussed turbines, turbine blades, and / or ducts. It is to be understood that the foregoing is a description of one or more embodiments of the invention. The invention is not limited to the particular embodiment(s) disclosed herein, but rather is defined solely by the claims below. Furthermore, the statements contained in the foregoing description relate to particular embodiments and are not to be construed as limitations on the scope of the invention or on the definition of terms used in the claims, except where a term or phrase is expressly defined above. Various other embodiments and various changes and modifications to the disclosed embodiment(s) will become apparent to those skilled in the art. All such other embodiments, changes, and modifications are intended to come within the scope of the appended claims. As used in this specification and claims, the terms "e.g.," “for example,” “for instance,” “such as,” and “like,” and the verbs “comprising,” “having,” “including,” and their other verb forms, when used in conjunction with a listing of one or more components or other items, are each to be construed as open-ended, meaning that the listing is not to be considered as excluding other, additional components or items. Other terms are to be construed using their broadest reasonable meaning unless they are used in a context that requires a different interpretation. In addition, the term “and / or” is to be construed as an inclusive OR. Therefore, for example, the phrase “A, B, and / or C” is to be interpreted as covering all the following: “A”; “B”; “C”; “A and B”; “A and C”; “B and C”; and “A, B, and C.”

Claims

CLAIMS 1. A method of optimizing the geometry of a ducted turbine (10) including a duct (16) and turbine blades (14) within the duct, the method comprising: (a) determining an efficiency of a low-fidelity model of the ducted turbine (10); (b) changing the geometry; (c) determining a change in the efficiency of the model after step (b); and (d) repeating steps (b) and (c) until the change in the efficiency is below a threshold value to arrive at an optimized geometry for the ducted turbine (10).

2. The method of claim 1, wherein step (b) comprises changing at least one of the following: a duct length (l), a duct diameter (D), a blade twist angle (θ), a blade cross-section (b).

3. The method of claim 1, wherein the ducted turbine (10) includes a turbine hub (36) within the duct (16) with the turbine blades (14) extending from the turbine hub, step (b) comprising changing at least one of the following: a hub length, a hub diameter, a duct thickness.

4. The method of claim 1, wherein an outlet diameter (Dexit) of the duct (16) and a tip-gap ratio (2R / dT) of the ducted turbine (10) remain unchanged during step (b).

5. The method of claim 1, wherein the low-fidelity model includes a blade-fixed rotating reference frame enclosing the turbine blades (14) and a stationary inertial reference frame outside the duct (16).

6. The method of claim 1, wherein the method includes use of a first-order numerical scheme to compute convective terms of Reynolds-averaged Navier-Stokes (RANS) equations.

7. The method of claim 1, further comprising parameterizing the geometry of the ducted turbine (10) before step (a) using free-form deformation (FFD), wherein control points defining FFD boxes are moved relative to one another in step (b) such that a computational mesh defined in and around the model is deformed when the geometry is changed.

8. The method of claim 7, the FFD boxes include a first FFD box (BD) enclosing the duct (16)and excluding the turbine blades (14), second FFD boxes (BB) enclosing each of the turbine blades (14) and excluding the duct (16), and a third FFD box (BP) enclosing the first and second FFD boxes.

9. The method of claim 1, further comprising determining the efficiency (CP) of a high- fidelity model of the ducted turbine (10) having the optimized geometry of step (d).

10. A method, comprising: (a) determining an optimized geometry for a turbine (12) using a low-fidelity model of the turbine; and (b) determining an efficiency of the turbine (12) using a high-fidelity model of the turbine with the optimized geometry.

11. The method of claim 10, wherein step (a) includes solving steady RANS equations using a multiple reference frame (MRF) model, and step (b) includes solving unsteady RANS equations using a rotating-sliding mesh (RSM) model.

12. The method of claim 10, wherein the turbine (12) is part of a ducted turbine (10).

13. A ducted turbine (10), comprising: a duct (16) having a length (l), a throat diameter (DT), and an exit diameter (Dexit); and turbine blades (14) within the duct (10), each turbine blade having an effective radius (R), a blade twist profile, and a blade section profile, wherein the length (l), the throat diameter (DT), the blade twist profile, and the blade section profile are configured to provide the ducted turbine (10) with an efficiency (CP) greater than 46% as determined by solving unsteady Reynolds-averaged Navier–Stokes (RANS) equations using a rotating-sliding mesh (RSM) model.

14. The ducted turbine (10) of claim 13, wherein the length (l) of the duct (16) is less than three times the effective radius (R) of each turbine blade (14).

15. The ducted turbine (10) of claim 13, wherein the length (l) of the duct (16) is less than twotimes the effective radius (R) of each turbine blade (14).

16. The ducted turbine (10) of claim 13, wherein the throat diameter (DT) is greater than 65% of the exit diameter (Dexit).

17. The ducted turbine of claim 13, wherein the throat diameter is greater than 70% of the exit diameter.

18. The ducted turbine (10) of claim 13, wherein each blade twist profile has a first twist angle (θ) in a range from 20 degrees to 30 degrees at r / R = 0.4 and a second twist angle in a range from 10 degrees to 20 degrees at r / R = 1, wherein r is a radial distance from a rotational axis of the turbine blades and R is the effective radius of the turbine blades.

19. The ducted turbine (10) of claim 18, wherein the first twist angle is in a range from 25 degrees to 30 degrees and the second twist angle is in a range from 10 degrees to 15 degrees.

20. The ducted turbine (10) of claim 19, wherein each blade twist profile has a first twist angle (θ) greater than 20 degrees at r / R = 0.4 and a second twist angle less than 20 degrees at r / R = 1, wherein r is a radial distance from a rotational axis of the turbine blades and R is the effective radius of the turbine blades.

21. The ducted turbine (10) of claim 13, further comprising a hub (36) from which the turbine blades (14) extend within the duct (16), wherein the hub extends from within the duct to a nose extending beyond an inlet (18) of the duct.

Images