A Computational framework to design optimally loaded supercavitating hydrofoils by differential evolution algorithm and a new viscous lifting line method

Advances in Engineering Software 133 (2019) 28–38

Contents lists available at ScienceDirect

Advances in Engineering Software journal homepage: www.elsevier.com/locate/advengsoft

A Computational framework to design optimally loaded supercavitating hydrofoils by differential evolution algorithm and a new viscous lifting line method

T

Giuliano Vernengo ,a, Luca Bonfigliob ⁎

a b

Dept. of Electric, Electronic and Telecomunication Engineering and Naval Architecture (DITEN), University of Genova, 16145, Italy MIT Sea Grant, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

ARTICLE INFO

ABSTRACT

Keywords: Differential evolution Supercavitating hydrofoil Cavitation Viscous lifting line Surface-piercing hydrofoil

The engineering design of a three dimensional submerged hydrofoil operating at very high speeds is obtained leveraging a Differential Evolution (DE) approach. The final goal is to identify the optimal load distribution over the span of a super-cavitating hydrofoil by using a design by optimization approach driven by hydrodynamic analysis of complex, turbulent, multi-phase flows. We achieve this goal by modeling the load distribution over the hydrofoil by means of a B-spline curve, which provides a rigorous parametric description of the hydrofoil operating conditions through the points of the load distribution control polygon. The parametric model includes design variables representing the most relevant hydrofoil shape parameters. We predict hydrodynamic performance by means of a Viscous Lifting Line method specifically conceived for the application targeted in the present study. This computational model accounts for the strong non-linear hydrodynamic characteristics of super-cavitating hydrofoils. We demonstrate the validity of the proposed design by optimization framework for high speed super-cavitating hydrofoils showcasing two design applications, namely a fully submerged hydrofoil operating close to a rigid boundary and a surface-piercing hydrofoil with variable dihedral angle. A statistical analysis of DE algorithm is performed to assess its performance on such an engineering design problem.

1. Introduction The quest for speed at sea is a challenging problem that has interested naval architects for more than 200 years. A significant technological revolution has been introduced before the World War I, when an aeronautical engineer equipped a small craft with lifting devices designed to sustain its weight when sailing at high speed [1,2]. The possibility of using lifting devices, such as submerged or surface-piercing hydrofoils, represents a great opportunity to drastically reduce the drag of a marine vessel by lifting it above the water surface. A renovated interest in hydrofoil supported vessels is demonstrated by the technology race triggered in 2013 among the America’s Cup sailing boats. The recently revealed hydrofoil supported monohull (AC75) confirms the technological direction taken in this past decade. Even though hydrofoil technology is well known since the beginning of the 20th Century, the speed of marine vehicles did not increase significantly since then. At very high speeds the pressure on lifting surfaces falls below the saturation pressure. The pressure drop triggers phase change which leads to water vaporization at ambient temperature. This

⁎

phenomenon is known as cavitation and it represents one of the main limiting factors for achieving higher speeds at sea. Most of the lifting devices commonly used to support the vessels weight are designed in order to avoid or at least delay cavitation, but when reaching speeds in excess of approximately 50 knots, cavitation can not be avoid, hence it becomes a design condition. A class of hydrofoils and propellers called super-cavitating have been specifically conceived to operate in cavitating flow regime. Conventional span-wise sections of super-cavitating hydrofoils are specifically designed to trigger and maintain a stable cavitating regime by inducing a localized pressure drop at a sharp leading edge and by ensuring flow separation and base ventilation by means of a blunt trailing edge. Examples of conventional super-cavitating hydrofoils were presented by Waid and Lindberg [3], Kermeen [4], Parkin [5]. The loss of performance in sub-cavitating conditions is one of the main flaws of blunt trailing-edge hydrofoils. In this respect, unconventional 2-D hydrofoils able to operate both in sub-cavitating and in super-cavitating regimes have been designed and optimized [6–8]. The design of super-cavitating 2-D and 3-D hydrofoils can be approached by using different computational methods, ranging from

Corresponding author. E-mail addresses: [email protected] (G. Vernengo), [email protected] (L. Bonfiglio).

https://doi.org/10.1016/j.advengsoft.2019.04.006 Received 21 February 2019; Received in revised form 26 March 2019; Accepted 15 April 2019 0965-9978/ © 2019 Elsevier Ltd. All rights reserved.

Advances in Engineering Software 133 (2019) 28–38

G. Vernengo and L. Bonfiglio

relatively simple closed solutions for conventional shapes [9,10] to potential flow-based methods [11–15] up to computationally expensive unsteady viscous solvers [16–18]. Comprehensive surveys of the DE state of the art and possible algorithm variants and hybridization have recently been proposed by Neri and Tirronen [19], Das and Suganthan [20], Das et al. [21]. Comparative studies of the DE with respect to other optimization algorithms such as particle swarm optimization and artificial bee colony, have been presented among the other by Vesterstrom and Thomsen [22], Civicioglu and Besdok [23] while numerical study on the DE control parameters can be found in Peñuñuri et al. [24]. The DE algorithm has been successfully used in several engineering fields, such as e.g. in aerodynamics [25–30], acoustics [29,31], power system planing [32,33], turbine optimization [34,35], water resources optimization [36–38], for the optimization of process chemical processes [39] and systems [40]. However, despite its proven effectiveness and the availability of its formulation, the DE algorithm has not been largely used for marine and offshore related researches and applications. Among the few published examples, Brizzolara [41] has been one of the first to use the DE algorithm to optimize the shape of a Small Waterplane Area Twin Hull (SWATH) ship with respect to resistance in calm water. Due to the good performance of the optimal hull shape found in the aforementioned study the DE algorithm has then been used as a reference for further comparative resistance and seakeeping analyses [42]. A design-by-optimization framework for unconventional super-cavitating 2-D profiles based on DE algorithm has recently been proposed in [8]. Hydrodynamic performance have been evaluated during the evolution of the optimization by anad-hoc developed Boundary Element Method (BEM). The final verification of the optimal candidate shapes, carried out by means of a high-fidelity Unsteady Reynolds Averaged Navier Stokes (URANS) solver, has confirmed the validity of the designs found by the DE algorithm. In this paper we aim to discover the optimal load distribution of a super-cavitating 3-D hydrofoil using a design-by-optimization approach that relies on the convergence characteristics of the Differential Evolution (DE) algorithm [43,44]. The optimal load distribution over the span of the hydrofoil maximizes the hydrofoil hydrodynamic efficiency by reducing the overall drag while ensuring a prescribed design lift. To achieve this target, we construct a parametric model that include both design and operating conditions variables. The design variables describe the shape of the super-cavitating hydrofoil through relevant shape parameters such as span, chord, taper ratio, etc. Operating conditions are described in terms of load distribution across the hydrofoil span, In addition to shape parameters, we model the load distribution along the hydrofoil span by using a B-spline curve. The control points of the B-spline polygon represent the free variables of the optimization problem describing the operating conditions in terms of lift force. The optimization algorithm is driven by the hydrodynamic performance of the hydrofoil that are evaluated by using a recently developed Viscous Lifting Line method [15] able to account for the strongly non-linear characteristics of super-cavitating profiles.

Fig. 1. Cavity shape experienced in super-cavitating regimes. The vapor pocket detaches from the hydrofoil sharp leading edge, closing far behind the blunt trailing edge. The hydrofoil suction surface is entirely enveloped inside the vapor cavity.

High speeds ensures a low pressure field over the entire suction surface of the hydrofoils. In particular, vapor production is initiated when the pressure field drops below the saturation pressure. The cavitation index σ compares the difference between the local and the saturation pressure with the dynamic pressure and it provides a metric for evaluating the cavitation regime:

=

patm + gh 1 2

pV

V2

(1)

In Eq. (1), patm, ρgh, pV and V∞ represent the atmospheric pressure, the hydrostatic pressure corresponding to the hydrofoil submergence h, the water saturation pressure and the hydrofoil advancing speed, respectively. For a given cavitation index, different hydrodynamic regimes are experienced for different angles of attack α. Where α indicates the inclination of the hydrofoil chord with respect to the upstream flow velocity direction, as displayed in the example of Fig. 1. Fig. 2 presents Computational Fluid Dynamic predictions of the lift coefficient of a 2-D super-cavitating hydrofoil operating in super-cavitating conditions ( = 0.3) for different angles of attack. The lift curve shows a strong non-linear trend with four major slope changes, highlighted in Fig. 2. Variations in the slope of the lift coefficient, ∂CL/∂α have a distinctive physical meaning. Slope value of ( CL/ = 2 ) corresponds to fully wet, base-cavitating flow, occurring when the cavity detaches from the blunt trailing edge developing only in wake of the hydrofoil. Slope value of ( CL/ = 3 ) corresponds to a partial cavitation regime, in

2. Supercavitating hydrofoils: physics backgrounds and numerical performance prediction Super-cavitating hydrofoils are lifting devices designed to operate at high speeds in stable cavitating flow regimes. In design conditions, the hydrofoil suction surface is fully enveloped in a vapor pocket detaching from the hydrofoil leading edge and closing several chords aft the blunt trailing edge (see Fig. 1). Unsteady cavitation might induce vibration, noise, material erosion due to the high temperatures occurring during the vapor bubble collapse and structural damage due to fatigue. Vapor cavity lengths exceeding twice the length of the hydrofoil chord safeguard agains flow reattachment that cause flow unsteadiness [45]. In order to achieve flow stability, we trigger cavity detachment by inducing a localized pressure drop, obtained through a sharp leading edge.

Fig. 2. Lift force expressed in terms of lift coefficient cL =

L 1 / 2 V 2 cb

across dif-

ferent angles of attack α for a cavitation index = 0.3. Solid black line represents CFD results obtained using a Reynolds Averaged Navier–Stokes solver. Different slopes corresponding to different flow regimes are represented by means of solid lines of different colors. 29

Advances in Engineering Software 133 (2019) 28–38

G. Vernengo and L. Bonfiglio

which the cavity detaches from the leading edge but does not extend across the entire hydrofoil chord (base cavitation still interests the wake of the hydrofoil). Slope value of ( CL/ = /4) corresponds to transition regimes, characterized by unstable vapor cavities. Slope value of ( CL/ = /3) corresponds to fully cavitating flow, occurring for high angles of attack. The design by optimization presented in this study is performed by evaluating objective function outputs using a Viscous Lifting Line (VLL) approach [15]. The lifting line approach enables significant computational savings however, the potential flow assumptions at the basis of its formulation, leads to inaccurate results in case of marked non-linear behavior such as the one presented in Fig. 2. The VLL approach used in the present study allows to include viscous effects and non-linearity, increasing the accuracy of the fluid dynamic predictions without the computational burden of 3-D viscous simulations. Viscosity effects are recovered by using a 2-D high-fidelity Unsteady Reynolds Averaged Navier-Stokes (URANS) solver based on the open-source libraries of Open-FOAM [46]. Viscous effects are introduced in the 3-D solution by a particular treatment of the control point position XC(y), i.e. the location where the influence coefficients are computed and the boundary conditions are imposed. Their locations are found based on the local slope of the non-linear lift curve according to Eq. (2).

X C (y ) =

1 CL ( + ) 4

circulation is the same as in CP1 but the span-wise location is a free parameter. This ensures the horizontal tangent at the root of the load distribution necessary to impose the rigid wall boundary condition in case where the hydrofoil is clamped at its root. In cases where the hydrofoil pierces the water surface (free surface), the load distribution vanishes at the root and the circulation intensity at CP1 must take zero values at the free surface. In this scenario, CP2,3,4 are allowed to move both along the horizontal and vertical directions. CP5 is placed at y/ b = 1, i.e. at the tip of the hydrofoil, hence only vertical displacements are allowed. The last control point, CP6, according to the physics, is used to set to zero the load distribution at the tip of the winglet. The distribution of circulation is normalized in order to obtain a maximum value of circulation corresponding to = 1.0 . The second set of parameters represents the plan shape geometry of the hydrofoil, represented in Fig. 5. This set includes the winglet span bWL, the hydrofoil span b, the root chord cRoot, the taper ratio = cTip/cRoot and the sweep angle Λ. The overall lift force developed by the hydrofoil is computed by integrating the distribution of circulation Γ(y) over the span of the hydrofoil. The Kutta–Joukowsky theorem in Eq. (3) establish a direct relationship between the area of the circulation distribution and the hydrofoil lift: i

(2)

L3D = V

Due to the induced angles generated by 3-D effects, each span-wise section of the super-cavitating hydrofoil could experience a different effective angle of attack αe with respect to the flow direction. According to this variation of the relative flow direction, each hydrofoil span-wise section operates at a different angle of attack and a different flow regime, hence it delivers different performance in terms of lift (CL) and drag (CD) coefficients. Vernengo et al. [15] validated this numerical method against experiments performed on a conventional super-cavitating hydrofoil [4]. The super-cavitating hydrofoil used in the experimental campaing by Kermeen [4] has been used in the study. It is characterized by a rectangular plan shape and a symmetric wedge-like span-wise sections with the apex angle at the leading edge (LE) equal to = 6 deg. Compared to both experimental measurements and higher-fidelity viscous solutions obtained using a URANS based method, the VLL provides accurate results as shown in Fig. 3 with significant computational savings. The accuracy of the solution over a wide range of angles of attack and cavitation indexes justify its employment in an automatic design approach driven by an optimization algorithm.

b 2

(y ) dy

b 2

(3)

In Eq. (3), V∞ represents the free stream velocity. According to Eq. (3), each circulation distribution identified using the DE algorithm, needs to be scaled in order to respect the constraints on the required lift CLTarget. The weighting factor is computed as the ratio between the area of the circulation distribution Γ(y/b)i and the required lift, as L = L3D / LTarget . B-spline curves are independent from the basis used for their representation, this allows us to turn a constrained optimization problem (where the constraint is the target lift) into an unconstrained optimization problem. Consequently, the parametric description provides circulation distributions that are feasible candidates for the DE algorithm. In the present study, we employ the VLL for the characterization of 3-D effects and the shape optimization of a super-cavitating hydrofoil having variable plan-form area and 3-D parameters and a constant span-wise section, defined by Kermeen [4]. Different combinations of parameters representing the circulation distribution Γ(y) produce different outputs in terms of induced drag Di and viscous shape drag Dv. Each span-wise section is in fact subjected to different induced angles depending on the particular circulation distribution. The induced drag Di, formulated in Eq. (4), is affected by the distribution of the induced angle αI(y), that in turn depends on the distribution of the collocation points XC(y). The shape drag Dv, Eq. (5), is proportional to the non-linear 2-D drag coefficient CD2D (y ), predicted using the URANS solver. Moreover, different combinations of parameters defining the plan shape of the hydrofoil produce different outputs in terms of both induced and viscous drag. For instance, increasing the sweep angle induces the center of hydrodynamic pressure to move towards the tip of the hydrofoil and increasing the dihedral angle reduces the vertical component of the hydrodynamic force (i.e. lift). The aforementioned complexities in the design by optimization of super-cavitating hydrofoils need to be addressed by using appropriate optimization algorithms driven by accurate hydrodynamic prediction tools, capable of solving complex cross interactions between parameters of different nature affecting hydrofoils performance.

3. Design by optimization of super-cavitating 3D hydrofoils The lifting line method allows to design a hydrofoil for given performance metrics (direct problem) and to analyze the performance given a specific design (indirect problem). In this paper, we solve the direct problem, searching for an optimal super-cavitating hydrofoil able to meet specific requirements in terms of lift force developed in a given operating condition. The plan-form area of the hydrofoil as well as the load distribution on the hydrofoil span are unknowns of the problem and are variables subject to optimization. The free variables of the optimization process are reported in Table 1 and in Table 2 as well as their bounds of variation. There are two sets of variables considered in the study. The first set describes the distribution of circulation, while the second set includes design parameters representing the main shape variables of the hydrofoil. Fig. 4 presents the B-spline curve used to model the circulation distribution. The B-spline polygon is determined by the position of six control points CP1: 6. The free variables correspond to the coordinates of the control points. With reference to Fig. 4, CP1 correspond to the circulation at the root of the hydrofoil and is free to move only in the vertical direction (corresponding to the circulation intensity). CP2 is free to move only horizontally, i.e. the intensity of the

Di = V

Dv =

30

1 2 V 2

b /2 b /2

(y ) i (y ) dy

b /2 b /2

CD2D (y ) c (y ) l (y ) dy

(4) (5)

Advances in Engineering Software 133 (2019) 28–38

G. Vernengo and L. Bonfiglio

Fig. 3. Lift and drag coefficient as predicted by the Viscous Lifting Line model for the reference super-cavitating hydrofoil [4]. Results are obtained for different angles of attack at two cavitation indexes: = 0.1 and = 0.4 . Solid red line represents predictions obtained using the VLL, while blu circles represent results obtained using a 3-D URANS model and green triangles represent experiments performed by Kermeen [4].

3.1. Definition of the objective function and set-up of the design process The goal of the optimization problem summarized in Eq. (6), is to find the parameter vector x* in the input N-dimensional space , maximizing the hydrodynamic efficiency expressed in terms of lift-todrag ratio, while ensuring requirements (LTarget) in terms of delivered lift force L. Since the DE is formulated for minimization of continuous functions, the objective function is computed as the inverse of the hydrodynamic efficiency. D to find: x * D D so that: (x *) (x ) L L subject to: L = LTarget

x (6)

In this paper we leverage the Differential Evolution (DE) algorithm as proposed in its original formulation [44] to perform the design by optimization of a complex engineering system. The optimization is driven by hydrodynamic performance predictions obtained using a Viscous Lifting Line (VLL), specifically formulated for the design and analysis of super-cavitating hydrofoils. Being N the number of free variables, the DE algorithm is initialized by using 10 · N designs per iteration, as suggested by Storn and Price [44]. However it has to be noted that a recent study by Neri and Tirronen [47] showed that reducing the population size might lead to

Fig. 4. Definition of the control polygon (red points, black dashed lines) and corresponding B-spline curve (red curve) for a generic distribution of circulation Γ(y/b) over the hydrofoil span b. Bounds of the input space are represented through the shaded regions representing the space containing all the possible position of the circulation distribution control points. variations of the control points are represented in gray in case of vertical control point motion allowed, in light-blue in case of vertical and horizontal displacements of the control point. 31

Advances in Engineering Software 133 (2019) 28–38

G. Vernengo and L. Bonfiglio

Table 2 Parametric model for the description of the hydrofoil geometry. Variable bounds are described by Min and Max columns. Each variable is described according to Fig. 5. Variable

Min

Max

Description

X8 X9 X10 X11 X12 X13

0.2 3.8 0.8 0.5 0.0∘ 5.0∘

0.3 4.2 1.0 1.0 10.0∘ 20.0∘

Winglet span, bWL Hydrofoil span, b Root chord, cRoot Taper ratio, τ Sweep angle, Λ Dihedral angle, Ψ

and the remaining set of design variables related to the load distribution (see Table 1) the hydrodynamic performance of the 3D SCH and the effective angle distribution of its profiles over its span are predicted by using the Viscous Lifting Line solver. The DE optimization algorithm pursues the search towards the optimal design by iteratively changing the design variables. 4. Results of the DE based design by optimization approach for supercavitating hydrofoils We leverage the capabilities of the presented design by optimization approach based on DE algorithm and VLL prediction method to optimize the designs of two different engineering applications consisting of a fully submerged and a surface piercing super-cavitating 3-D hydrofoil. The main conceptual difference is the boundary condition imposed at the hydrofoil root. The fully submerged hydrofoil is modeled using half of the span and imposing a rigid wall boundary condition at its root. The surface piercing hydrofoil span is entirely modeled, fulfilling a free surface (or biplane) boundary condition at its root. This difference is presented in Fig. 7, where the two boundary conditions are imposed by using the method of images. According to this technique, the hydrofoil is mirrored with respect to the plane where the boundary is located. The vortex system necessary to formulate the lifting line method [15] follows the same transformation. In order to model a rigid wall the mirrored vortex system is oriented so that the trailing vortexes at the root will vanish while the opposite happens in the case of the free surface. As a result of the lifting line model, the circulation at the root of the hydrofoil operating close to a rigid wall is positive while that of a surfacepiercing hydrofoil is zero. In the following Sections 4.1 and 4.2 we present results of the design by optimization of a fully submerged hydrofoil and a surface-piercing supercavitating hydrofoil, respectively.

Fig. 5. Schematic representation of the geometrical parameters described the supercavitating hydrofoil geometry. Top panel represents front view showing the span of the hydrofoil b and the dihedral angle ψ. bWL indicates the span of the winglet. Bottom panels represents plan view. Angle Λ defines the leading edge (LE) sweep. Trailing edge (TE) distribution and tip chord cTip are determined from the taper ratio distribution, given a certain root chord value (cRoot). Table 1 Parametric model for the description of the circulation distribution. Variable bounds are described by Min and Max columns. Each variable is described according to Fig. 4. Variable

Min

Max

Description

X1 X2 X3 X4 X5 X6 X7

5%b 40%b 82%b 0.2 0.3 0.3 0.3

38%b 80%b 96%b 1.0 1.0 1.0 0.8

XCP2 XCP3 XCP5 YCP1 YCP3 YCP4 YCP5

4.1. Tests n.1: optimal design of a 3D submerged hydrofoil

2

The optimal distribution of circulation over the span corresponds to the distribution that generates the lower induced angles, hence the lower induced drag. Classic finite wing theory [48] demonstrates that the optimal load distribution has an elliptic shape. In case of a hydrofoil clamped to a rigid wall, the circulation Γ(y/b) becomes a semi-ellipse having the minor axis coincident with the hydrofoil root. The first design case gives a unique possibility of testing the capabilities of the DE algorithm in finding the optimal elliptic load distribution for a complex hydrofoil shape operating in challenging flow conditions. With reference to a benchmark super-cavitating hydrofoil having rectangular plan shape and aspect ratio AR= 4 [4], we optimize the hydrofoil efficiency of a super-cavitating hydrofoil operating in stable cavitating regime at angle of attack = 5 and cavitation index = 0.1. Lift and drag coefficients predicted using the VLL method for the benchmark case correspond to CL = 0.198 and CD = 0.031, respectively. These values accurately match the measurements [15]. The hydrodynamic efficiency of the benchmark hydrofoil is L /DB = 6.387 and it is used to normalize the quantity of interest driving the objective

better DE performance in many applications. A maximum of 40 iterations is allowed. Considering N = 13, a total amount of 5,200 designs are generated and evaluated before convergence is reached. The crossover probability and the weighting factor have been set to CR = 0.7 and F = 0.6, respectively. The DE / rand to best /1/ bin strategy has been followed. Fig. 6 displays the whole process exploited by the computational design framework. In the propose release the shape of the 2D profiles is selected a priori, hence not being involved in the design process, and the corresponding hydrodynamic performance are computed by a CFD solver or taken from EFD results off-line. The geometry of the 3D Super Cavitating Hydrofoil (SCH) is then defined based both on the characteristics of the selected 2D profile and on the set of geometric parameters listed in Table 2. Taken as input such a geometric information 32

Advances in Engineering Software 133 (2019) 28–38

G. Vernengo and L. Bonfiglio

Fig. 6. Flowchart of the process exploited by the computational design framework.

Fig. 7. Different modeling of the boundary conditions representing a fully submerged hydrofoil clamped on rigid wall (left panel) and a hydrofoil piercing the water surface (right panel). Green dots represent the training vortex cores where the trailing vortexes are applied.

function:

QoI : =

(D / L) (D / L)B

Fig. 8. Submerged super-cavitating hydrofoil optimization history driven by the normalized hydrofoil efficiency as defined in Eq. (7). Benchmark design is represented by means of purple square. Feasible designs are represented by red circles, while the Optimum design by a blue triangle.

(7)

A fixed value corresponding to = 0 deg is imposed to variable X13 representing the hydrofoil dihedral angle. Fig. 8presents the evolution of the objective function during the optimization process. The computational time required to perform the optimization presented in Fig. 8 amounts to approximately 6 h on whole design process has taken about 6 h on a single core, 2.5 GHz Intel Core i5 processor. The trend of the optimization history shows that the generated solutions are confined within a narrow range of the objective function values as 3000 design evaluations are achieved. This demonstrates the successful convergence of the optimization algorithm. Compared to the benchmark (FObj = 1), the optimized design improves the hydrodynamic efficiency of about 18%.

As expected, the optimization algorithms drives the circulation distribution towards an elliptic shape as demonstrated in Fig. 9. Even though the parametric models allows for elliptic circulation distributions, this particular shape can be generated only by hydrofoils having an elliptic plan-shape. The limited number of parameters used for the description of the 3-D hydrofoil shape and specifically the linear chord distribution over the hydrofoil span do not allow to obtain elliptic plan form areas, hence we cannot obtain perfectly elliptic circulation distributions. Load distributions differences, highlighted in Fig. 9, are due 33

Advances in Engineering Software 133 (2019) 28–38

G. Vernengo and L. Bonfiglio

Fig. 9. Distribution of circulation of the optimum hydrofoil close to wall boundary (light blue curve) compared to the elliptic distribution (black dashed curve) and to the reference hydrofoil (purple curve).

Fig. 12. Surface piercing, super-cavitating hydrofoil optimization history driven by the hydrofoil efficiency L/D. Blue circles represent feasible designs while red triangle represents optimum hydrofoil. Green empty circles represents the best-so-far designs identified for each generation.

Fig. 10. Distribution of induced (αi) and effective angles (αe) over the span of the optimum fully submerged, super-cavitating hydrofoil identified using the DE algorithm. The small discrepancy between the effective and the 2-D angle of attack distributions demonstrates the significant reduction of 3-D effects.

Fig. 11. Comparison of the hydrofoil sections at the sweep and dihedral angles.

y b

Fig. 13. Distribution of circulation of the optimum surface-piercing supercavtitaing hydrofoil (red curve) compared to an equivalent elliptic distribution having the same maximum value (black dashed curve). Higher load closer to the hydrofoil root induce more significant 3-D effects as demonstrated in Fig. 14.

= [0.0, 0.5, 1.0], excluding

to the combination of the other hydrofoil shape parameters, namely X8 to X12, i.e. they are related to the fact that the actual hydrofoil plan shape is not elliptic. Nevertheless, the trend of the optimization algorithm in searching for a elliptic distributions, confirms the consistency of the optimization problem and the convergence characteristics of the DE algorithm. The optimal load distribution results in very small induced angles αi, as shown in Fig. 10, demonstrating a significant reduction of 3-D effects. The hydrofoil sections at y/ b = (0.0, 0.5, 1.0) are compared in Fig. 11 excluding both the sweep and the dihedral angles. As a result, the distribution of the effective angle is very close to the distribution of the 2-D angle of each section α2D, being e = 2D + i .

Fig. 14. Distribution of induced (αi) and effective angles (αe) over the span of the optimum surface-piercing super-cavitating hydrofoil identified using the DE algorithm. 3-D effects are particularly pronounced at the hydrofoil root (close to the free surface) as demonstrated by higher discrepancy between the effective and the 2-D angle of attack distributions.

by optimization of a surface-piercing super-cavitating hydrofoil. The optimization framework is based on DE algorithm driven by an objective function predicted using the aforementioned VLL. Here, the dihedral angle Ψ has been included as free variable, and it is exactly evaluated with respect to the free surface level. Variables bounds are the

4.2. Tests n.2: optimal design of a 3D surface-piercing hydrofoil The second application presented in this paper consists of the design 34

Advances in Engineering Software 133 (2019) 28–38

G. Vernengo and L. Bonfiglio

Fig. 15. Comparison of the hydrofoil sections at the sweep and dihedral angles.

y b

same as previously (Test no. 1 in Section 4.1). Fig. 12 presents the optimization history in terms of absolute values of hydrofoil efficiency L/D. As demonstrated in Fig. 12, the DE algorithm is able to discover a design having approximately 33% of improvement in the hydrodynamic efficiency. The mean value of L/D at the end of the initial Design of Experiments is approximately L/D≅4, while at the end of the optimization process we discovered an hydrofoil having L /D = 5.31. As in designing the submerged super-cavitating hydrofoil, the optimization history reveals a stable convergence to the optimum, completing the whole design process with approximately the same computational burden. In Fig. 13 we compare the optimal distribution of circulation to an equivalent elliptic shape having the same maximum value. The load vanishes at the free surface ( y/ b = 0 ) since it is almost equivalent to an open end, such as for instance the tip of the hydrofoil. With respect to the elliptic load distribution, DE discovers an optimal super-cavitating surface-piercing hydrofoil for which the peak of the circulation is shifted towards the root while decreasing the load at the hydrofoil tip, consistently. Having circulation peak closer to the free surface ( y/ b = 0 ) corresponds to higher induced angles in proximity of the hydrofoil root (αi≅2∘). Induced angles decrease moving towards the hydrofoil tip (αi≅0∘), as shown in Fig. 14. The hydrofoil sections at y/ b = (0.0, 0.5, 1.0) are compared in Fig. 15 excluding both the sweep and the dihedral angles. Compared to the submerged hydrofoil, 3-D effects are more relevant in case of hydrofoils piercing the free surface. Table 3 summarizes the sets of parameters as well as the hydrodynamic performance of the two optimized hydrofoils corresponding to the fully submerged and the surface-piercing hydrofoils, shown in Fig. 17. Optimization trends for the parameters describing the shape of the 3-D hydrofoil are displayed in Fig. 16. Both for the fully submerged and the surface-piercing optimal hydrofoils are characterized by similar values of root chord (X10) and taper ratio (X11). These values correspond to the lower bounds of design space indicating the possibility of discovering further improvements for hydrofoils out from the design space. The DE algorithm drives variable X13, corresponding to the

= [0.0, 0.5, 1.0], excluding

Table 3 Comparison of the free variables and performance between the optimum designs from Test n.1 and Test n.2, respectively. Variable

Opt. Wall

Opt. Free Surface

X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 X11 X12 X13 CL CDi CDv CD L/D

2.56E−01 6.41E−01 8.69E−01 8.69E−01 6.98E−01 4.61E−01 2.45E−01 2.98E−01 4.19E+00 8.04E−01 5.02E−01 9.33E+00 0.00E+00 1.98E−01 6.41E−04 2.56E−02 2.63E−02 7.54E+00

0.16E−01 7.27E−01 8.71E−01 9.29E−01 7.28E−01 4.68E−01 2.58E−01 2.98E−01 4.20E+00 8.00E−01 5.31E−01 1.87E+00 6.85E+00 1.98E−01 3.32E−03 3.40E−02 3.73E−02 5.31E+00

Fig. 16. Scatter-plots representing the trends of the input parameters during the optimization performed using the DE algorithm. Blue dots represent feasible designs, red dots represent for each variable the value corresponding to the optimal hydrofoil.

35

Advances in Engineering Software 133 (2019) 28–38

G. Vernengo and L. Bonfiglio

4.3. Sensitivity analysis on DE parameters Choosing the combination of DE control parameters does not follow a general rule but it is recognized to be strongly problem-dependent [49]. This means that there should be a setting that works on a specific problem that instead might not be optimal for a different one. In this perspective, in fact, lot of effort has been done to create variants of the DE algorithm where control parameters are set by automatic criteria and self-adapted during the optimization process [50–52]. A systematic sensitivity analysis on the two main control parameters of the DE algorithm has been carried out considering the design of surface-piercing supercavitating hydrofoil, i.e. on the problem defined in the Test n.2. Combinations of the weighting factor F and the crossover probability CR have then been considered in this analysis. According to the selected strategy, the first parameter is used to tune the mutation operation carried out by the DE by the following Eq. (8):

vi, G = xi, G + F ·(xbest , G

xi, G ) + F ·(x rn1, G

x rn2, G )

(8)

being vi,G the mutant design vector with respect to each target design vector xi,G at the Gth iteration of the algorithm. xbest,G, xrn1,G and xrn2,G are the best design and two randomly selected designs of the same iteration. The uniform crossover is then applied after the mutation to each pair of target and mutant designs according to Eq. (9):

vi,jG, if randj [0, 1) ui,jG

=

CR or j = jrand

xi,jG, otherwise j = 1, 2, ..,D

(9)

where jrand is a random integer in the range [1, D] and CR = [0, 1) manages the portion of parameters retained from the mutant design. Three values for each of the two parameters have been chosen corresponding to CR = [0.5; 0.7; 0.9] and F = [0.4; 0.6; 0.8]. According to the formulation proposed by Zaharie [53] and considering the selected population size (NP = 130), the three weighting factors are much higher than the minimum critical scaling factor for all the values of CR. Each of the nine combinations of CR and F has been tested on a complete optimization problem as described in the previous Section 4.2. All the processes have been run for 40 iterations resulting in 5200 overall design evaluations. Fig. 18 displays the results of such a sensitivity analysis in terms of the best-so-far design, i.e. the best design found at the end of each iteration. Each sub-plot corresponds to one of the three values of F. Each curve refers to a specific value of CR. The same optimum design is found by using both F = 0.4 or F = 0.6. The parameter CR affects the convergence rate of the process but not the final result over the selected number of iterations. Considering e.g. the cases having the lower values of F, increasing CR seems to fasten the convergence of the optimization while the opposite happens for the intermediate F = 0.6, at least before reaching 30 iterations. Analyzing the highest value of F, it is evident the stagnation of the optimization processes despite the values of CR.

Fig. 17. Comparison of the two optimal SC-hydrofoils resulting from the two applications of the optimization based design framework. Top: perspective view. Bottom: plan view.

surface-piercing dihedral angle towards low values, correctly capturing the decrease in the net vertical force (lift) due to the inclination of the hydrofoil with respect to the water surface. Sweep angle Λ presents an analogous trend. As expected both the hydrofoil and the winglet spans, b and bWL, respectively, reach the maximum allowed values in order to reduce the effects of the finite wing.

Fig. 18. Best-so-far design evolution for different combinations of F and CR. Curves corresponding to CR = [0.5; 0.7; 0.9] are shown on each sub-figure. 36

Advances in Engineering Software 133 (2019) 28–38

G. Vernengo and L. Bonfiglio

prediction tools is used for this time in the context of shape optimization in a multi-dimensional setting. The optimal shape discovered using the DE algorithm for the optimization of a fully submerged hydrofoil. is compared to theoretical results, demonstrating the consistency of the optimization settings as well as the definition of the input space. In particular the load distribution over the span of the optimal hydrofoil fits almost entirely the elliptic distribution of circulation, that is the one minimizing induced angles and drag. Compared to a benchmark hydrofoil, the optimal design discovered using the DE presents an hydrodynamic efficiency improvement of about 18%. When used to optimize the shape of a surface-piercing hydrofoil, the DE algorithm discovered a different circulation distribution and a set of geometrical parameters leading to a significant improvement of the hydrodynamic efficiency. In particular, considering an initial mean value of the lift-to-drag ratio of approximately L/D≅4, the optimized design bring the performance to an improvement of about 33%. The differential Evolution algorithm employed in this study revealed very successful, providing impressive improvements in the hydrodynamic efficiency even in the setting of a complex engineering problem such the one presented in this paper.

Fig. 19. Mean best-so-far design evolution over 20 design processes. Standard deviation is displayed by using the error-bar on each point. The min-max range at each iteration is highlighted by the shaded area. Table 4 Comparison of the optimal designs over 20 repeated design processes. Statistic

Value

Mean Std. Dev Min Q1 Median Q3 Max

8.238E−01 2.668E−03 8.1670E−01 8.215E−01 8.231E−01 8.254E−01 8.280E−01

Acknowledgments This project has been partially funded by DARPA EQUiPS grant HR0011517798. Authors wish to thank Prof. G.E. Karniadakis at Brown University for his valuable scientific guidance. Supplementary material Supplementary material associated with this article can be found, in the online version, at 10.1016/j.advengsoft.2019.04.006

Despite its simplicity, this analysis proves one again that the selection of the DE control parameters might strongly affect the result of the optimization process, surely changing the convergence performance of the algorithm in the search of the optimum. Considering the stochastic evolution of the DE algorithm in combining the agents xi,G, the convergence properties of the design framework have been verified by running 20 identical design processes of the surface-piercing supercavitating hydrofoil. In order to limit the computational burden, the repeated design processes have been carried out for 30 generations by using a reduced population made of 65 designs, i.e. 5 · N. The evolution of the mean best-so-far design is shown in Fig. 19. The mean value, the standard deviation, the minimum value, the first quartile value, the median, the third quartile value and the maximum value found at the end of the eight repeated processes are reported in Table 4. The general trend of the mean best-so-far design highlights the convergence of the process which is further proven by the reduction of both the standard deviation and the min-max range as the iterations increase.

References [1] Crewe PR. The hydrofoil boat; its history and future prospects. Quart Trans Inst Naval Archit 1958;100(4). [2] Clark DJ, Ellsworth WM, Meyer JR. The quest for speed at sea. Technical Digest, Naval Surface Warfare Center, Carderock Division. 2004. [3] Waid R.L., Lindberg Z. Experimental and theoretical investigations of a supercavitating hydrofoil1957. [4] Kermeen R.W. Experimental investigations of three-dimensional effects on cavitating hydrofoils1960. [5] Parkin B.R. Experiments on circular arc and flat plate hydrofoils in noncavitating and full cavity flows1956;(47-6). [6] Brizzolara S. A new family of dual-mode super-cavitating hydrofoils. Fourth international symposium on marine propulsors (SMP), Austin, TX, May. 2015. p. 1–9. [7] Brizzolara S, Bonfiglio L. Comparative CFD investigation on the performance of a new family of super-cavitating hydrofoils. J Phys 2015;656:012147. IOP Publishing [8] Vernengo G, Bonfiglio L, Gaggero S, Brizzolara S. Physics-based design by optimization of unconventional supercavitating hydrofoils. J Ship Res 2016;60(4):187–202. [9] Johnson VE. Theoretical determination of low-drag supercavitating hydrofoils and their two-dimensional characteristics at zero cavitation number. National Advisory Committee for Aeronautics; 1957. [10] Tulin MP. Supercavitating flows - small perturbation theory. J Ship Res 1964;7(3):16–37. [11] Fine NE, Kinnas S. A boundary element method for the analysis of the flow around 3-d cavitating hydrofoils. J Ship Res 1993;37(3):213–24. [12] Mishima S, Kinnas SA. A numerical optimization technique applied to the design of two-dimensional cavitating hydrofoil sections. J Ship Res 1996;40(1):28–38. [13] Young Y, Kinnas S. Analysis of supercavitating and surface-piercing propeller flows via bem. Comput Mech 2003;32(4–6):269–80. [14] Royset J, Bonfiglio L, Vernengo G, Brizzolara S. Risk-adaptive set-based design and applications to shaping a hydrofoil. J Mech Des 2017;139(10):101403. [15] Vernengo G, Bonfiglio L, Brizzolara S. Supercavitating three-dimensional hydrofoil analysis by viscous lifting-line approach. AIAA J 2017:1–15. [16] Bonfiglio L, Brizzolara S. A multiphase RANSE-based computational tool for the analysis of super-Cavitating hydrofoils. Naval Eng J 2016;128(1):47–64. [17] Bonfiglio L, Perdikaris P, Brizzolara S, Karniadakis G. Multi-fidelity optimization of super-cavitating hydrofoils. Comput Methods Appl Mech Eng 2017. [18] Bonfiglio L, Royset JO, Karniadakis G. Multi-disciplinary risk-adaptive design of super-cavitating hydrofoil. 2018 AIAA non-deterministic approaches conference. 2018. p. 1177. [19] Neri F, Tirronen V. Recent advances in differential evolution: a survey and

5. Conclusions In this paper, we have presented a design framework demonstrating the global convergence capabilities of the Differential Evolution algorithms in discovering optimal design solutions for complex engineering systems. We proved the effectiveness of the optimization framework showcasing two different design problems consisting of a fully submerged and a surface-piercing super-cavitating hydrofoil. We consider a multidimensional input design space including parameters representing the shape of the load distribution and relevant geometrical parameters describing the shape of the hydrofoil. We leveraged a novel numerical approach that allows to predict hydrofoil hydrodynamic performance by using a Viscous Lifting Line, specifically designed to capture the nonlinear characteristics of super-cavitating hydrofoils. This hydrodynamic 37

Advances in Engineering Software 133 (2019) 28–38

G. Vernengo and L. Bonfiglio experimental analysis. Artif Intell Rev 2010;33(1–2):61–106. [20] Das S, Suganthan PN. Differential evolution: a survey of the state-of-the-art. IEEE Trans Evol Comput 2011;15(1):4–31. [21] Das S, Mullick SS, Suganthan PN. Recent advances in differential evolution–an updated survey. Swarm Evol Comput 2016;27:1–30. [22] Vesterstrom J, Thomsen R. A comparative study of differential evolution, particle swarm optimization, and evolutionary algorithms on numerical benchmark problems. Evolutionary computation, 2004. CEC2004. congress on. 2. IEEE; 2004. p. 1980–7. [23] Civicioglu P, Besdok E. A conceptual comparison of the cuckoo-search, particle swarm optimization, differential evolution and artificial bee colony algorithms. Artif Intell Rev 2013;39(4):315–46. [24] Peñuñuri F, Cab C, Carvente O, Zambrano-Arjona MA, Tapia J. A study of the classical differential evolution control parameters. Swarm Evol Comput 2016;26:86–96. [25] Rogalsky T, Kocabiyik S, Derksen R. Differential evolution in aerodynamic optimization. Can Aeronaut Space J 2000;46(4):183–90. [26] Madavan N. Aerodynamic shape optimization using hybridized differential evolution. 21st AIAA applied aerodynamics conference. 2003. p. 3792. [27] Madavan N. On improving efficiency of differential evolution for aerodynamic shape optimization applications. 10th AIAA/ISSMO multidisciplinary analysis and optimization conference. 2004. p. 4622. [28] Asouti VG, Giannakoglou KC. Aerodynamic optimization using a parallel asynchronous evolutionary algorithm controlled by strongly interacting demes. Eng Optim 2009;41(3):241–57. [29] Clifton-Smith M. Aerodynamic noise reduction for small wind turbine rotors. Wind Eng 2010;34(4):403–20. [30] Carrigan TJ, Dennis BH, Han ZX, Wang BP. Aerodynamic shape optimization of a vertical-axis wind turbine using differential evolution. ISRN Renew Energy 2012;2012. [31] Snellen M, Simons DG, Van Moll C. Application of differential evolution as an optimisation method for geo-acoustic inversion. Proceedings of the 7th European conference on underwater acoust.(Delft, The Netherlands, 2004). 2004. p. 721–6. [32] Wang Y, Li B, Weise T. Estimation of distribution and differential evolution cooperation for large scale economic load dispatch optimization of power systems. Inf Sci 2010;180(12):2405–20. [33] Wang Y, Li B, Zhang K. Estimation of distribution and differential evolution cooperation for real-world numerical optimization problems. Evolutionary computation (CEC), 2011 IEEE congress on. IEEE; 2011. p. 1315–21. [34] Song L, Feng Z, Li J. Shape optimization of turbine stage using adaptive range differential evolution and three-dimensional navier-stokes solver. ASME Turbo expo 2005: power for land, sea, and air. American Society of Mechanical Engineers; 2005. p. 1033–40. [35] Joly MM, Verstraete T, Paniagua G. Differential evolution based soft optimization to attenuate vane–rotor shock interaction in high-pressure turbines. Appl Soft Comput

2013;13(4):1882–91. [36] Suribabu C. Differential evolution algorithm for optimal design of water distribution networks. J Hydroinf 2010;12(1):66–82. [37] Sedki A, Ouazar D. Hybrid particle swarm optimization and differential evolution for optimal design of water distribution systems. Adv Eng Inf 2012;26(3):582–91. [38] Dong X-l, Liu S-q, Tao T, Li S-p, Xin K-l. A comparative study of differential evolution and genetic algorithms for optimizing the design of water distribution systems. J Zhejiang Univ Sci A 2012;13(9):674–86. [39] Babu B, Sastry K. Estimation of heat transfer parameters in a trickle-bed reactor using differential evolution and orthogonal collocation. Comput Chem Eng 1999;23(3):327–39. [40] Angira R, Babu B. Optimization of process synthesis and design problems: a modified differential evolution approach. Chem Eng Sci 2006;61(14):4707–21. [41] Brizzolara S. Parametric optimization of swat-hull forms by a viscous-inviscid free surface method driven by a differential evolution algorithm. Proceedings 25th symposium on naval hydrodynamics, St. Johns, Newfoundland and Labrador. 5. 2004. p. 47–64. [42] Brizzolara S, Vernengo G, Bonfiglio L, Bruzzone D. Comparative performance of optimum high speed swath and semi-swath in calm water and in waves. Transactions 2015;123:273–86. [43] Storn R, Price K. Minimizing the real functions of the icec’96 contest by differential evolution. Evolutionary computation, 1996., proceedings of IEEE International Conference on. IEEE; 1996. p. 842–4. [44] Storn R, Price K. Differential evolution–a simple and efficient heuristic for global optimization over continuous spaces. J Global Optim 1997;11(4):341–59. [45] Acosta AJ. Hydrofoils and hydrofoil craft. Ann Rev Fluid Mech 1973;5:161–84. [46] OpenFOAM Foundation. OpenFOAM user guide. OpenFOAM Foundation; 2014. [47] Neri F, Tirronen V. On memetic differential evolution frameworks: a study of advantages and limitations in hybridization. Evolutionary computation, 2008. CEC 2008.(IEEE world congress on computational intelligence). IEEE congress on. IEEE; 2008. p. 2135–42. [48] Anderson Jr JD. Fundamentals of aerodynamics. Tata McGraw-Hill Education; 1985. [49] Gämperle R, Müller SD, Koumoutsakos P. A parameter study for differential evolution. Adv Intell Syst Fuzzy Syst Evol Comput 2002;10(10):293–8. [50] Tirronen V, Neri F, Rossi T. Enhancing differential evolution frameworks by scale factor local search-part I. Evolutionary computation, 2009. CEC’09. IEEE congress on. IEEE; 2009. p. 94–101. [51] Qin AK, Huang VL, Suganthan PN. Differential evolution algorithm with strategy adaptation for global numerical optimization. IEEE Trans Evol Comput 2009;13(2):398–417. [52] Zhang J, Sanderson AC. Jade: adaptive differential evolution with optional external archive. IEEE Trans Evol Comput 2009;13(5):945–58. [53] Zaharie D. Critical values for the control parameters of differential evolution algorithms. Proceedings of MENDEL. 2. 2002. p. 6267.

38