CFD-based design optimization of ducted hydrokinetic turbines

This section describes the methodology in optimization and CFD evaluations. An overall flowchart is shown in Fig. 5. The whole process involves the optimization and the re-evaluation of the optimized design using higher-fidelity simulations. The optimization and high-fidelity re-evaluation use DAFoam⁴³ and OpenFOAM⁴², respectively. In what follows, we will describe each component of the methodology in subsections. To provide a self-contained but easy-to-follow paper for readers, we put additional details in the Supplementary Information and keep the main paper as concise as possible. We start from CFD models involved in both the optimization and re-evaluations and then follow up with other components in the optimization framework.

CFD models

The governing equations for the flow field around the turbine (Fig. 1) are the Navier–Stokes equations

where (varvec{U}) is the flow velocity and p is the pressure.

We consider the Reynolds-averaged Navier–Stokes (RANS) equations with grids only resolving the averaged components of the flow. One can apply the Reynolds decomposition (varvec{U} = overline{varvec{U}}+varvec{u}’) (and the same for pressure) to Eq. (4). (overline{varvec{U}}) denotes the averaged velocity in a time window or by an ensemble and (varvec{u}’) represents the zero-mean turbulent fluctuation. This leads to the unsteady RANS equation:

where (overline{varvec{u}’varvec{u}’}) is the Reynolds stresses that need to be approximated by turbulence models. In this work, the (k-omega) SST turbulence model is used together with the automatic near-wall treatment (see Supplementary Information A for details on both).

The rotating blades are handled in simulations by two blade-resolved approaches: the Multiple Reference Frames (MRF) and the rotating-sliding mesh approach (RS). The former is used for steady RANS solutions (i.e., a solution with time-derivative terms set to zero) with multiple different reference frames, while the latter is used directly in the unsteady solution of Eq. (5). The MRF is used in optimization. The RS is defined as the higher-fidelity approach and used for optimized result re-evaluations (see Fig. 5). A detailed introduction of the two approaches is included in the following sections.

Multiple reference frames (MRF) method

The MRF method is an efficient method for modeling turbomachinery flow. In the MRF method, the computational mesh stays stationary, and the rotational effect is handled through a rotational reference frame. In particular, the fluid domain is separated into two regions: a rotational region surrounding the turbine blades with a blade-fixed reference frame, and the remaining stationary region with an inertial reference frame, as shown in Fig. 6. In both regions, the flow is considered steady with respect to the corresponding reference frame, so only steady RANS equations need to be solved.

MRF method illustration. The computational domain is split into stationary and rotating (shaded blue) regions. The inertial and rotating reference frames are denoted by gray and yellow axes, respectively.

To be more specific, in the rotational region, the blades are stationary and experience a steady inflow. The flow velocity in the blade-fixed reference frame can be expressed by

where (varvec{U}) is the velocity in the inertial reference frame, (varvec{Omega }) is the rotation vector of the turbine blades, and (varvec{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 need to be established with the blade-fixed reference frame, which requires further formulations of both Eq. (5) and the (k-omega) SST model equations. Although the implementation in OpenFOAM/DAFoam solves this complete set of equations, only the rotational-region formulation regarding Eq. (4) is presented here to provide the key insights of the method.

Combining Eq. (4) and Eq. (6), we obtain (see Supplementary Information B for a detailed derivation)

The steady equations solved in the rotational region are Eq. (7), with (partial varvec{U}_R/partial t=0). Therefore, in the MRF method, steady versions of Eq. (4) and Eq. (7) are solved in stationary and rotating regions. Solving these steady equations 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 because of two issues. First, the rotational and stationary regions are usually chosen in a subjective manner. There is no guarantee that the rotational region covers all 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. Secondly, 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 solution, as in the RANS-MRF method, can lead to significant errors for this type of unsteady flow. As a result, the RANS-MRF method is considered a lower-fidelity model in the context of this paper.

Rotating-sliding mesh approach

The rotating-sliding mesh (RS) allows for the direct simulation of the unsteady RANS (URANS) equations (Eq. (5)) with mesh domains that exhibit relative motion. This is needed for modeling rotating geometries. The underlying idea of this method is to allow a region of the computational mesh to rotate with the turbine blades, as illustrated in Fig. 7. The rotating mesh also creates a technical problem that the mesh at the rotating/non-rotating interface becomes non-conformal, i.e., the nodes at two sides of the interface do not match up. The data transfer across the interface, therefore, needs to be handled by a special interpolation method involving a “supermesh”⁴⁵, as described in Supplementary Information C.

Rotating (blue) and stationary (gray) mesh, with non-conformal mesh interfaces.

Coupling URANS with the RS is, in principle, much more accurate than the RANS-MRF method as it captures the unsteady nature of the flow. This can be critical in simulating flow around a ducted turbine because the possible flow separation from the duct surface can be captured better. However, in the URANS-RS, a small time step is needed to resolve the blade rotation, so a long simulation time is required for the solution to reach a quasi-steady state. The computational cost of the URANS-RS is hence much higher than that of the RANS-MRF method. In this paper, the URANS-RS is considered a higher-fidelity method and is only used in re-evaluating the performance of optimized designs.

Mesh configuration

The unstructured computational mesh is generated using the OpenFOAM meshing tool snappyHexMesh. A mesh overview is shown in Fig. 8. The size of the computational domain is (10.4Dtimes 10.4Dtimes 23.7L), where (D=sqrt{(4/pi )A}=1.536) m is the maximum diameter of the duct, and (L=2.107) m is the length of the duct that is taken from the baseline design. This domain size is sufficient to avoid a blockage effect upon tests. For the boundaries of the domain, mixed boundary boundary conditions are applied. To be more specific, fixed velocity values and zero pressure gradient are imposed for flows coming inside the domain, whereas a zero velocity gradient and fixed pressure value are imposed for flows coming out of the domain. On the turbine blades and duct, no-slip boundary conditions are applied.

Mesh configuration for a ducted turbine, with (a) front view, (b) side view, (c) close-up side view to show the refinement region, (d) close-up view at the 0.7R blade section, and (e) close-up view near the duct. The levels of refinements are also shown in the figures.

To model the flow near the turbine and immediately downstream with higher accuracy, a region of (2Dtimes 2Dtimes 2.5L) around the turbine is refined (see Fig. 8c). Within the refinement region, we apply three steps of grid refinement. First, a level-4 refinement is implemented for the full refinement region, i.e., each cell of the original mesh is divided into ((2)^4) cells. Then, further refinements, up to level-6 and level-9 toward respectively the duct and blade surfaces are applied in the sub-regions close to the surfaces. Finally, prism layers are added close to the surfaces, which contain further continuously refined cells toward the surfaces (2 and 3 layers are used for the former and latter, both with the expansion ratio of 1.1). The prism layer provides better resolution for boundary layers and is critical for obtaining well-convergent results in the grid sensitivity study shown later in the Results and Discussion section.

In this work, three grid resolutions, M0, M1, and M2, are used with an increasing number of cells, i.e., further refinement from M0 to M2. The coarsest grid M0 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, but it roughly contains 1.5 million cells in the refinement region and a similar number of cells in the remaining region). In M1 and M2, the full mesh region (including background mesh and refinement region) is uniformly refined in each direction by the 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 grids M0, M1, and M2 are used in the URANS-RS re-evaluation, including the grid sensitivity study.

Optimization

Sequential Quadratic Programming (SQP)⁴⁶, implemented in SNOPT⁴⁷, is used to solve the optimization problem. One of the challenging tasks is to obtain the gradient of the objective function (nabla C_P) with respect to all design variables. The following subsections discuss the components (see Fig. 5) involved in the gradient computation: (1) geometry parametrization and mesh deformation; (2) adjoint method. The gradient information is then used in the SQP algorithm to obtain the next design points. The procedure iterates until satisfying convergence criteria or further design improvements are not achievable.

Geometry parametrization via FFD method

It is necessary to parametrize the geometry to deform the surface mesh. In this work, the Free-Form Deformation (FFD) method⁴⁸ is used for the geometry parametrization, implemented in the package pyGeo^49,50. The principle of the FFD method is to enclose the surface mesh nodes 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 surface nodes using tri-variate B-splines. More details are presented in Supplementary Information D. Controlling the FFD points enables smooth deformation of the enclosed geometry. Fig. 9 shows two examples of geometry deformation controlled by the FFD method in two and three dimensions.

Examples of geometry deformations with FFD points. (a) A 2D NACA0012 airfoil (from the DAFoam tutorial) and (b) a 3D Stanford bunny (figure taken from Kenway et al.⁴⁹).

Figure 10 shows an overview of the FFD setup for our ducted turbine. Two levels of FFD boxes are used, with one parent box (black) enclosing all duct and blade geometries and two children boxes (red and blue) enclosing the duct and blades. The 21 design variables in the optimization problem can now be represented by 21 degrees of freedom (DoF) associated with the FFD points. The child FFD box for an individual blade has 32 FFD points placed on 8 sections. Twist variables rotate the four FFD points about the reference axis located at the quarter chord line. Scale variables scale the cross-section by moving the four control points to expand or contract simultaneously. The FFD points across different blades are linked to ensure the same deformation for all blades.

Geometry parametrization of a ducted turbine through the FFD method. Top: different views of the overall FFD setup, where the parental FFD box (black) controls the radial scales of both duct and blades, the child duct FFD box (red) controls the length of the duct, and the child blade (blue) FFD box controls the pitch/twist angles and sectional areas of the blades. Bottom: Closer views of the FFD setup for the blades.

The child FFD box for the duct contains 112 FFD points placed on seven sections, but overall only one variable is defined to control all FFD points to change the duct length l. When changing the duct length, all duct FFD points move in the axial direction with perturbations proportional to their distances from the throat. This movement is to ensure that the throat is consistently located at 26.4% of the overall duct length. Note that seven sections are not necessary. This choice is mainly for convenience during setup.

The parent FFD box handles the constraint (3f) on the tip-gap ratio and the condition (D_text {exit}=sqrt{(4/pi )A}=1.536) m. Twenty-eight FFD points are placed on seven sections in the parent box, in which the FFD points for the last three sections are closely packed horizontally and remain stationary throughout the optimization, such that (D_text {exit}=sqrt{(4/pi )A}) is guaranteed. The 16 FFD points on the first 4 sections from the duct inlet are used to control the duct diameters, i.e., design variables (d_i). As these FFD points move radially, the child FFD boxes (enclosed in the parent FFD box) deform and move the embedded surface geometry accordingly. Therefore, the constraint (3f) is automatically satisfied since the blade expands/contracts proportionally to the duct throat. This implementation of the complex FFD setup for ducted turbines, we believe, is a novel application.

Adjoint method for derivative computation

To compute the derivative of an objective function (in this case, the power coefficient (C_P)) with respect to design variables, the adjoint method⁴⁶ is used. In order to clearly explain the adjoint method, notations are introduced first as follows: Let (varvec{x}equiv { { theta _i, b_i }_{i=1}^{8}, {d_j}_{j=1}^{4}, l} in {mathbb {R}}^{N_x}) with (N_x=21) as the design variables. Let (varvec{s}in {mathbb {R}}^{M}) be the state variables in the solution of the RANS-MRF equation. Here (M sim {mathscr{O}}(10^7)) includes three velocities and pressure at each cell in the computational grid. Our goal is to compute (dC_P/dvarvec{x} in {mathbb {R}}^{21}). If a finite-difference method is used to compute the derivative, one needs at least 22 CFD simulations for derivative computation, even with the lowest-order approximation. This is computationally prohibitive for our application.

For the adjoint method to compute (dC_P/dvarvec{x}), the first step is to write the function (C_P(varvec{x},varvec{s}(varvec{x}))), and express its total derivative with respect to (varvec{x}) as

where (partial C_P/partial varvec{x}) should be considered as the change of power coefficient (C_P) as the design variables (i.e., geometry) are varied, with flow solution (varvec{s}) remaining unchanged. (partial C_P/partial varvec{s}) is the change of (C_P) as the flow solution (varvec{s}) changes with a fixed turbine geometry. These partial derivatives are relatively easy to compute, with more details presented in Supplementary Information E.

The term that is difficult to compute in Eq. (8) is (dvarvec{s}/dvarvec{x}). To compute it, one needs to further involve the RANS-MRF state equations in terms of their discretized residual form (varvec{R}(varvec{x},varvec{s}(varvec{x}))=0). Here (varvec{R}(varvec{x},varvec{s}(varvec{x})) in {mathbb {R}}^M) considering the same number of equations as the number of unknowns in (varvec{s}). Since (varvec{R}(varvec{x},varvec{s}(varvec{x}))) should remain zero with a change of (varvec{x}) (if the flow solution is correctly obtained), we have

Direct solution of Eq. (9) gives

It is worthwhile to discuss the computational cost associated with Eq. (10) at this point. The matrix multiplication in Eq. (10) leads to a computational complexity of ({mathscr{O}}(M^2 N_x)) that is very expensive since (Msim {mathscr{O}}(10^7)) and (N_x) is also large. This has to be added by the cost to invert a (Mtimes M) matrix, which is, in general, more expensive. Even if one uses some iterative solver for linear systems to solve Eq. (9), the procedure needs to be repeated for (N_x) times since (dvarvec{s}/dvarvec{x}) (as well as the RHS) has (N_x) columns. The computation is, therefore, also very expensive.

On the other hand, the computational cost can be significantly reduced by simply substituting Eq. (10) to Eq. (8) and considering a re-grouping of the multiplications:

Instead of computing Eq. (10), the multiplication grouped in the parenthesis in Eq. (11) is computed first. This computation can be done by solving the so-called adjoint equation (the adjoint is equivalent to the transpose of a real matrix in our case)

whose solution transpose provides

as the parenthesis term in Eq. (11).

The solution of Eq. (12) involves solving a linear system only once, instead of (N_x) times as needed for Eq. (9), and hence is much less expensive (also compared to the direct computation of Eq. (10)). The computational cost to solve Eq. (12) is generally similar to the RANS-MRF computation. Therefore, in each iteration of the optimization, the computational cost is in the same order as one RANS-MRF solution. The only remaining component is the calculation of partial derivatives (partial varvec{R}/partial varvec{s}) in Eq. (12), which can be found in Supplementary Information E with other derivatives mentioned above.

Volume mesh deformation

When marching to the next design point, the turbine geometry is deformed. The entire computational volume mesh is deformed accordingly.

The volume mesh deformation is computed based on the analytic inverse-distance weighting method⁵¹, implemented in the IDWarp package⁵². Given a 2D surface, for example, a blade surface, with N surface mesh nodes, the geometry deformation leads to the movement of each node. Two quantities ((M_i, b_i)) are assigned for each node with (i=1,2,…,N), where (b_i) is the translation distance of the node and (M_i) is the rotation matrix such that (n_i^{new}=M_i n_i^{old}) with (n_i^{new}) and (n_i^{old}) the normal vectors at the node. In particular, both (n_i)’s are computed by a weighted average of the normal vectors for all surrounding cell faces of the node. After ((M_i, b_i)) 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., (Delta varvec{r}=sum _{i=1}^N w_i (M_i varvec{r} + b_i – varvec{r})) with (varvec{r}) being the coordinates of a volume node and (Delta varvec{r}) its movement. The weighting factor (w_i) has the empirical form^51,52 that grows in a polynomial form with the inverse of the distance between the volume and surface nodes. Figure 11 shows the deformation of the computational mesh during a ducted turbine optimization as an example.

Computational mesh deformation during the geometry deformation of a ducted turbine.