E Clément

A Görgen

W Korten

A Bürger

A Chatillon

Y Le Coz

Ch Theisen

M Zielinska

Bertram Blank

P J Davies

M Zielińska

S P Fox

J Gerl

G Georgiev

S Grévy

J Iwanicki

D G Jenkins

F Johnston-Theasby

P Joshi

I Matea

P J Napiorkowski

F De

Oliveira Santos

M G Pellegriti

R Wadsworth

A New Device for Combined Coulomb Excitation and Isomeric Conversion Electron Spectroscopy with Fast Fragmentation Beams

Keywords: Coulomb excitation, electromagnetic transition probability, conversion-electron spectroscopy, E0 transitions PACS: 25.70.De, 29.30.Dn, 27.50+e

published or not. The documents may come

A new device for combined Coulomb excitation and isomeric conversion electron spectroscopy with fast fragmentation beams

Introduction

A new device is presented for the combined spectroscopy of prompt γ rays and isomeric conversion electrons after reactions using fast beams from a fragmentation facility. It has been developed for the physics case of nuclei in the mass region A ≃ 70 close to the N = Z line, where low-lying 0 + states related to shape coexistence and shape isomerism are very common. Such low-lying 0 + states have been observed experimentally in the even-even Ge isotopes from A = 68 to 78, in the Se isotopes from mass A = 70 to 82, and in the Kr isotopes from mass A = 72 to 84 [START_REF] Firestone | Table of Isotopes[END_REF].

The N = Z nuclei 68 Se and 72 Kr are also of major importance for the rapid proton capture (rp) nucleosynthesis process on the neutron-deficient side of the nuclear chart [START_REF] Schatz | [END_REF]. Both nuclei are thought to be waiting points in the rp process: the synthesis of heavier elements can only proceed after β decay of rather long-lived isotopes or the sequential capture of two protons. Nuclear structure information of low-lying isomeric states is also an important input parameter for simulations of the rp process [3].

Nuclei in the middle of the f p shell have a large quadrupole deformation in their ground state, and adding or removing a few nucleons can change the nuclear shape dramatically. The self-conjugate nucleus 72 Kr is predicted to have an oblate deformation in its ground state [START_REF] Bengtsson | Nuclear Structure of the Zirconium Region[END_REF][START_REF] Petrovici | [END_REF][6][7], while the heavier Kr isotopes were shown to have prolate ground-state deformation and excited states with oblate shape [8]. The first experimental evidence for this shape coexistence was the observation of low-lying 0 + 2 states [9]. A similar scenario is expected for the light Se isotopes. However, the situation is less clear experimentally, since no low-lying 0 + 2 state has been observed in 68 Se.

If a low-lying 0 + state in 68 Se is located below or just above the first 2 + state, its decay occurs (exclusively or predominantly) via an E0 transition, i.e. through the emission of a conversion electron. The excited 0 + state then becomes meta-stable with a lifetime of the order of a few to several hundred nanoseconds, depending on the overlap of the wave functions of the two 0 + states. If the nuclei are produced in a high-energy fragmentation reaction (E > 40A MeV), they have a high probability of being fully stripped of their atomic electrons. In this case the E0 decay is blocked (if the excitation energy of the excited 0 + state is insufficient for internal pair creation), and the decay can only proceed once the nuclei have regained their atomic (inner-shell) electrons, e.g. after the implantation in a detector.

Intermediate-energy Coulomb excitation after fragmentation is a well-known technique used to measure the reduced transition probability B(E2; 0 + 1 → 2 + 1 ) for exotic nuclei far from stability, which gives a first indication on the collectivity of the nucleus [10]. Specific equipment has been designed at GSI [11], RIKEN [12], MSU [10], and previously at GANIL [13] for such measurements. Taking advantage of the high beam energy, thick targets can be used to increase the interaction probability. Coulomb excitation events can be separated from other reaction mechanisms by selecting only small scattering angles and thus large impact parameters. Contrary to Coulomb excitation at low energy, where multiple-step excitation is possible, single-step excitations dominate at energies well above the Coulomb barrier, usually simplifying the data analysis. The existence of a low-lying (metastable) 0 + 2 state can, however, influence the measured Coulomb excitation cross sections and has to be taken into account explicitly. This effect is particularly important if the fragmentation products, which are to be excited in the electromagnetic field of a heavy target nucleus, have been produced with a certain probability in the meta-stable state. In that case the isomeric ratio of the beam has to be measured simultaneously with the Coulomb excitation probability in order to extract B(E2) values (see e.g. the discussion in Ref. [14]).

Consequently a new device was designed to perform Coulomb excitation experiments of fragmentation products at intermediate energies with the LISE3 spectrometer [15] at GANIL. In addition to the standard detection systems for the tracking and identification of the incident and scattered ions, and for the γ rays emitted in-flight after Coulomb excitation, a dedicated setup was installed to measure conversion electrons emitted from unscattered beam particles in silicon detectors after implantation of the ions in a plastic scintillator. In this way the probability of the fragmentation products to be in an isomeric state could be measured and taken into account for the determination of B(E2) values obtained from the Coulomb excitation process.

Experimental procedure and description of the setup

The exotic nuclei were produced by the fragmentation of an intense primary 78 Kr 33+ beam at 70.1A MeV, delivered by the coupled GANIL cyclotrons, on a natural Ni production target located between the two superconducting solenoids of the SISSI device [16]. The target of 125 ñm thickness was inclined by 25 • and followed by a carbon stripper foil of 10 mg/cm 2 . The reaction products of interest were selected and purified in the LISE3 spectrometer [15], operated in achromatic mode with a Be degrader of 215 ñm thickness. The thickness of both target and degrader was optimized for the production of nuclei on the N = Z line. With a primary beam intensity of 2.4 eñA the intensity of the secondary 68 Se beam on the Coulomb excitation target was ∼100 particles per second. The secondary beam energy was around 40A MeV, depending slightly on the setting of the spectrometer. The composition of the cocktail beam was regularly measured with reduced intensity using a silicon telescope placed before the Coulomb excitation setup. The complete setup is illustrated schematically in Fig. 1 and described in more detail in the following sections. The associated electronics was based on an analog VXI-VME system and used the standard data-acquisition system of the GANIL facility.

Beam tracking

The incoming particles were tracked by a set of two micro-channel plate detectors (MCP) [17] of 65×65 mm 2 size, measuring secondary electrons produced in a thin emissive foil. The time difference between the first detector and the radio frequency of the cyclotrons provided the time of flight (ToF). The distance between the two detectors was 25 cm and the second detector was located 29 cm upstream from the target. Both detectors have a horizontal and vertical position resolution of about 1 mm, allowing the reconstruction of the position and the incident angle of the incoming particles on the target. In this way it was ensured that all nuclei were Coulomb excited under the same conditions. The reconstructed beam spot on the target and the distribution of incident angles are shown in Fig. 2. The beam profile is characterized by a quasi-Gaussian distribution with a horizontal and vertical width of

Detection of prompt γ rays

The nuclei in the cocktail beam were Coulomb excited on a thick 208 Pb target of 220 mg/cm 2 thickness placed in the center of an array of four large segmented clover detectors from the EXOGAM array [18]. Each detector comprises four HPGe crystals and each crystal is segmented into four sectors. The detectors were mounted without BGO escapesuppression shields due to mechanical constraints. The distance between the front face of the detectors and the target was 50 mm. An absolute photo-peak efficiency of 9% for 1.3 MeV γ rays was measured using a 60 Co source. The Doppler correction of the γ rays which were emitted in flight by the projectile was performed on an event-by-event basis using the position information of the EXOGAM clover detectors and the velocity vector of the outgoing particles.

Detection of γ rays from long-lived isomers

The setup can also be used to measure γ rays of isomeric states with lifetimes longer than the time of flight through the spectrometer, which is ∼ 1.2 ñs.

In this configuration, the standard silicon telescope of the LISE beam line is introduced about 1 m upstream from the Coulomb excitation setup. The incident ions are in this way identified by a ∆E-ToF measurement, and they are sufficiently degraded in energy after passing through the Si detectors to be stopped and implanted in the 208 Pb target, which is surrounded by the highly efficient array of four EX-OGAM clover detectors. This technique can be used to search for long-lived γ-decaying isomers in exotic nuclei. Furthermore, known isomeric γ rays can be used to identify nuclides unambiguously in the ∆E-ToF plot. This was for example possible using the 9/2 + isomer in 69 Se, which has a half live of T 1/2 = 960 ns and decays via a 534 keV M 2 transition [19]. The γ-ray spectrum of 69 Se nuclei implanted into the 208 Pb target is shown in Fig. 3.

Detection of scattered particles

At high beam energies a pure electromagnetic excitation process can only be ensured by selecting large impact parameters corresponding to small scattering angles. The scattered particles were detected and identified in a telescope of two highly segmented double-sided annular silicon strip detectors placed 42 cm downstream from the 208 Pb target. The two silicon detectors were segmented into 16 concentric rings on the front side and 16 azimuthal sectors on the back. Both detectors had active surfaces with inner and outer radii of 11 and 35 mm, respectively. The first silicon detector covered scattering angles between 1.5 • and 4.76 • in the laboratory frame, ensuring that the impact parameter was always greater than 17 fm for all nu- reactions in the secondary target such as transfer can thus be excluded. An identification plot of the energy loss ∆E and the time of flight between the cyclotron RF and the first MCP detector is shown in Fig. 4. Nuclides with the same isospin have a similar time of flight through the spectrometer, but can clearly be separated by the ∆E measurement.

Isomer detection system

Non-interacting beam particles which passed through the central hole in the annular silicon detectors were detected in a telescope of two plastic scintillators of 300 ñm thickness, placed 62 cm behind the 208 Pb target. The second scintillator was tilted by 45 • with respect to the beam axis. The thickness of the scintillators was chosen such that the ions were implanted in the center of the second scintillator. Due to the limited energy resolution, a Z identification via a ∆E -E measurement could not be achieved. However, the detectors were used to determine the total number of implanted nuclei. The exact beam composition was determined in regular intervals of a few hours by introducing a standard stack of silicon detectors into the beam upstream from the secondary target. In this way the absolute rate of the different nuclides incident on the target could be monitored.

The second scintillator was placed in the center of a modified version of the BEST array [20], originally designed for the spectroscopy of transfermium isotopes. Two alcohol-cooled silicon detectors of 1 mm thickness and divided into four quadrants each were facing the tilted surfaces of the second scintillator as illustrated in Fig. 1. Each silicon detector had an active surface of 50 × 50 mm 2 . The distance between the two Si detectors was 50 mm. The beam particles arrived on the secondary target at 40A MeV fully stripped of their atomic electrons. Most of the ions regain electrons only when they are implanted in the plastic scintillator. Ions arriving in an excited 0 + state which decays exclusively via an E0 transition are therefore trapped; the decay can only proceed after implantation, no matter how long the lifetime or the flight time through the spectrometer is. Even if there is an E2 branch depopulating the 0 + state, the partial lifetime can become longer than the time of flight. The array of silicon detectors around the implantation scintillator is therefore well suited to search for low-lying excited 0 + states, populated in the fragmentation reaction, by detecting the conversion electrons of the E0 decay.

The upper silicon detector suffered from lowenergetic δ electrons produced by the ions passing through the first scintillator. The lower silicon detector is mostly shielded by the tilted scintillator, resulting in much cleaner spectra. In cases of longlived isomers decaying by low-energy γ rays, it was possible to detect these γ rays in the silicon detectors, for example the 5.7 ñs isomer in 73 As, decaying by a 360 keV γ ray.

Commissioning with stable 78 Kr and 72 Ge beams

The commissioning of the setup was performed using 78 Kr and 72 Ge beams. The B(E2; 0 + 1 → 2 + 1 ) values are known and very large. The first excited state in 72 Ge is a 0 + 2 state at 691 keV, making 72 Ge the ideal case to test the setup for the detection of isomeric conversion electrons. Instead of using the primary 78 Kr beam directly, a secondary beam of 78 Kr was produced under conditions similar to those of the later measurement with radioactive nuclides. The beam was tracked by the MCP detectors and the velocity vector of the incident beam particles reconstructed. The γ rays emitted after Coulomb excitation of the beam were detected in coincidence with a scattered particle identified as 78 Kr in the annular silicon telescope. The Doppler correction was performed using the angular information of both the annular silicon detector and the segmented outer electrodes of the clover detectors. The calculated value of the beam velocity after passing the 208 Pb target was used since the lifetime of the 2 + 1 state in 78 Kr is longer than the time of flight through the target. For nuclei with shorter lifetimes, the velocity at mid-target was applied. Due to the large collectivity of 78 Kr, the 2 + 1 state is strongly populated. The prompt Doppler-corrected γ-ray spectrum gated on 78 Kr is presented in Fig. 5. Since there is only one γ ray emitted in the vast majority of Coulomb excitation events, a large fraction of the background was suppressed by requiring that only one germanium crystal gave signals. Indications for the population of the 2 + 2 state and two-step processes populating the 4 + 1 state are presented in the inset of Fig. 5. The 3 -state in 208 Pb at 2.6 MeV was also populated, but the peak is too spread out in energy to be seen in the spectrum which was Doppler corrected for 78 Kr. A resolution of 30 keV was achieved for the 2 + 1 → 0 + 1 transition at 455 keV, mainly limited by the velocity spread in the target and the opening angle of the clover detector segments.

In the same way, the prompt Doppler-corrected γ-ray spectrum gated on 72 Ge was obtained. As previously discussed, 72 Ge has a first excited state of spin-parity 0 + . As the nuclei are fully stripped, the 72 Ge nuclei arrive with a certain percentage in the isomeric state; this isomeric ratio has to be taken into account in the Coulomb excitation analysis. The unscattered 72 Ge ions were implanted into the scintillator where they regained their atomic electrons allowing the decay of the 0 + 2 state to proceed by conversion electrons. The spectrum obtained with the lower silicon detector is presented in Fig. 6. It shows very clearly the 0 + 2 → 0 + 1 E0 decay of 691 keV in 72 Ge. The spectrum is gated on the time of flight for nuclei with T z = 4. However, the energy resolution of the scintillator was insufficient to resolve the 72 Ge fragments from neighboring 74 As and 70 Ga. Since the beam composition is known, the isomeric ratio can be determined, and it was found that about 3% of the incident 72 Ge ions arrive in the isomeric 0 + 2 state.

The measurements with 78 Kr and 72 Ge have shown that the setup is well suited to perform Coulomb excitation experiments at intermediate energy and measure simultaneously the ratio of ions arriving in an isomeric 0 + state. This allows determining B(E2) values for exotic nuclei in the presence of low-lying 0 + states.

Measurements with Radioactive Beams

As the next step the spectrometer setting was optimized for several radioactive nuclei with mostly known B(E2) values. The prompt Dopplercorrected γ-ray spectra in coincidence with 66 Ge, 62 Zn and 68 Se are presented in Fig. 7. Due to electronics and data acquisition problems in this part of the experiment, only a small fraction of the recorded events were correlated, and the signals from the outer contacts of the clover detectors were not coded at all, leading to a severe loss of statistics and a degraded energy resolution. The 2 + 1 → 0 + 1 transitions are nevertheless clearly visible for 66 Ge and 62 Zn. In the case of 68 Se, no significant peak is visible above background level at the 2 + energy of 853 keV. An upper limit of 40 counts can be extracted from the spectrum gated on 68 Se for the 2 + 1 → 0 + 1 transition. In order to determine unknown B(E2) values and to minimize systematic errors, nuclei with well known B(E2) values are used as a reference to normalize the unknown values. For this procedure, the γ-ray intensities must be corrected for angular distribution effects and for the energy dependence of the γ-ray efficiency. As the B(E2) of 66 Ge and 62,64 Zn nuclei are known, they could be used as a reference to determine an upper limit for the B(E2) value for 68 Se. The data from the 78 Kr and 72 Ge measurements could not be used for normalization as the experimental conditions were not identical due to the electronics and data acquisition problems. The B(E2) values were extracted using an approximation of the exact Coulomb excitation cross sections at relativistic energies via the so-called "recoil correction" method developed by Alder and Winther [21]. It was shown [22] that this correction is the best approximation in the energy range between 20 and 40A MeV used in the present experiment.

The results of the present measurement are summarized in table 1. The systematics of B(E2; 2 + 1 → 0 + 1 ) values for the Zn, Ge, and Se chains is presented in Fig. 8. The values derived from the present measurement are shown as full symbols. The results are fully compatible with the known B(E2) values for 64 Zn and 66 Ge (open symbols). The value for 62 Zn, however, is found to be smaller than the literature value based on the measured lifetime of 4.2(3) ps The data illustrate the strong collectivity of the isotopes at mid-shell, which decreases toward the shell closures at 28 and 50. A pronounced irregularity is observed for the Se isotopes, which show an unexpected staggering at N = 36 and 38 [26]. However, the relevant lifetimes have been re-measured recently and the data indicate that 70 Se is less collective than previously thought [27]. The upper limit for the B(E2) value of 68 Se supports the interpretation of a smooth and symmetric behavior with a maximum of collectivity at N = 40 and 42, which is decreasing with neutron number. The new value for 62 Zn, B(E2; 0 + 1 → 2 + 1 ) = 0.074(15) e 2 b 2 , suggests that also for the Zn isotopes the collectivity decreases stronger with neutron number than previously thought.

The array of plastic scintillators and silicon detectors was used to search for the predicted 0 + 2 shape isomer in 68 Se by detecting conversion electrons emitted from the 68 Se ions implanted into the plastic scintillator. Assuming that a similar fraction of 68 Se ions is produced in an isomeric state as for 72 Ge (∼3%), one would expect to observe about 400 conversion electrons in the lower silicon detector, ignoring the upper one which suffers from a higher background of δ electrons. Since no sign of an E0 transition was observed in the spectrum conditioned on the ToF for T z = 0 nuclei and an energy loss for Z ≈ 34, the presence of a 0 + 2 state as the first excited state in 68 Se can be excluded. Even a 0 + 2 state situated slightly above the 2 + 1 state should have been observed due to the expected long partial lifetime in that case. If a low-lying 0 + 2 state exists in 68 Se, it can be presumed to lie well above the 2 + 1 state, to which it can decay through γ rays. This is the case in neighboring 70 Se, where an excited 0 + 2 state has been observed at 2011 keV excitation energy, decaying to the 2 + 1 state at 945 keV [26]. In order to search for possible isomers decaying by γ rays, a stack of silicon detectors was inserted into the beam line before the 208 Pb target, which was used to identify the ions event by event and to degrade their energy so that they were implanted in the Pb target. With this method the 960 ns isomer in 69 Se was observed (see Fig. 3). Again, the spectrum conditioned on 68 Se showed no sign of an isomer, so that the presence of a γ-decaying isomer with a lifetime longer than ∼ 1 ñs in 68 Se can be excluded.

The experimental results for 68 Se and in particular the non-observation of a low-lying 0 + 2 state stress the importance of further experimental studies for the understanding of the predicted shape coexistence in the light Se isotopes.

Summary

A new setup was developed for intermediateenergy Coulomb excitation combined with spectroscopy of isomeric conversion electrons. Coulomb excitation of stable and radioactive beams in the A ∼ 70 region at 40A MeV has been performed following the fragmentation of an intense 78 Kr primary beam of 70A MeV. The setup was tested using secondary 78 Kr and 72 Ge beams. The latter served also to test the setup for the detection of isomeric conversion electrons from the 0 + 2 state in 72 Ge. It was found that 3% of the secondary beam arrived in the isomeric 0 + 2 state, the decay of which is blocked during the flight through the spectrometer. B(E2) values were measured for several radioactive secondary beams. The B(E2) value for 62 Zn was found to be lower than the literature value. Only an upper limit for the B(E2) value could be derived for the N = Z nucleus 68 Se. The existence of a low-lying, long-lived 0 + 2 state with a lifetime longer than ∼ 1 ñs in 68 Se could be excluded.

Fig. 1 .

 1 Fig. 1. Schematic view of the experimental setup (not to scale).

Fig. 2 .

 2 Fig. 2. Beam profile of the 68 Se projectiles on the 208 Pb target reconstructed from the MCP beam tracking detectors (left) and distribution of the incident angles of the projectiles on the target (right).

  σ h = 6 mm and σ v = 15 mm, respectively, and a mean incident angle of 2 • .

Fig. 3 .

 3 Fig. 3. Spectrum showing the long-lived 9/2 + isomer (T 1/2 = 960 ns) in 69 Se after implantation of the ions into the thick secondary target.

Fig. 4 .

 4 Fig. 4. Identification matrix of the energy loss in the silicon detector against the time of flight between the cyclotron RF and the beam tracking detector. The setting of LISE3 was optimized for 68 Se.

Fig. 5 .

 5 Fig. 5. Prompt γ-ray spectrum in coincidence with scattered 78 Kr projectiles. The inset shows the same spectrum on an enlarged scale. At low energy the spectrum is dominated by Pb x-rays from the target.

Fig. 6 .

 6 Fig.6. Conversion-electron spectrum in coincidence with 72 Ge nuclei showing the E0 decay of the first excited 0 + state. The peak at low energy is due to atomic electrons.

Fig. 7 .

 7 Fig.7. Prompt Doppler corrected γ spectra for 66 Ge (top), 62 Zn (center), and 68 Se (bottom). The energies of the 2 + 1 → 0 + 1 transitions are marked by arrows. Other peaks are due to x-rays from the 208 Pb target, the pile-up of x-rays, and the 511 keV annihilation line, which appears as a doublet due to the Doppler correction.

Fig. 8 .

 8 Fig. 8. Systematics of the B(E2; 2 + 1 → 0 + 1 ) values for Zn, Ge, and Se isotopes. Values taken from Ref. [23] are shown as open, data from the present work as full symbols. The data points are slightly displaced horizontally for clarity of the presentation.

Table 1 Deduced

 1 B(E2; 2 + → 0 + ) values for 62 Zn, 64 Zn, and 66 Ge. The results are compared to literature values. Only an upper limit can be deduced for 68 Se.

	Nucleus	B(E2; 0 + 1 → 2 + 1 ) [e 2 b 2 ]
		Nuclear Data	this work
	62 Zn	0.124(9) [23]	0.074(15)
	64 Zn	0.149(5) [24]	0.165(21)
	66 Ge	0.099(19) [23]	0.093(12)
	68 Se		≤0.157
	[25]. Consequently, only the values for 64 Zn and
	66 Ge were used for the normalization. The resulting
	B(E2; 2 + 1 → 0 + 1 ) value for 62 Zn is 10.2(21) Weis-skopf units. For 68 Se, an upper limit of B(E2; 2 + 1 → 0 + 1 ) < 19 W.u. is extracted.

Acknowledgments

We would like to thank the technical staff at GANIL for excellent support during the preparation and running of the experiment. This work has been supported in part by the European Community FP6 -Structuring the ERA -Integrated Infrastructure Initiative -contract EURONS No. RII3-CT-2004-506065, by the Polish State Committee for Scientific Research under contract No. 2 P03B 007 25, by the German BMBF under contract 06BN-109, and the UK EPSRC.


Nitin Garg

Gaetan K W Kenway

Joaquim R R A Martins

Yin L Young
email: ylyoung@umich.edu

Y Mpo

O Tatin

High-fidelity Coupled Hydrostructural Optimization of a 3-D Hydrofoil

Keywords: hydrostructural design optimization, gradient-based optimization, cavitation, hydrofoil, propulsor

. The coupled hydrostructural tool has been integrated with the adjoint method to enable efficient high-fidelity, optimization with large number of shape design variables. The coupled hydrostructural optimization study for an unswept, cantilevered, tapered 3-D NACA 0009 hydrofoil made out of solid aluminum was carried out. A total of 210 shape design variables were used in the optimization study, with constraints on the lift coefficient (C L ), cavitation number, and maximum stress (based on the fatigue strength of aluminum). The hydrostructurally optimized foil results in an increase in efficiency (C L {C D ) of 12.4%, a reduction in foil mass of 7.3%, and an increase in cavitation inception speed by 45%, over the original NACA 0009 hydrofoil, while satisfying the cavitation constraint and the maximum stress constraints.

INTRODUCTION

The shape of a vehicle that moves within a fluid (air or water) has a large impact on the vehicle performance. Typically, a low drag and mimimal weight are desired, and in many cases, a certain amount of lift is required. The development of computational fluid dynamics (CFD) based aerodynamic shape optimization over the last few decades has been motivated by these goals, and has had a large influence on how aircraft are designed [START_REF] Hicks | Wing design by numerical optimization[END_REF][START_REF] Jameson | Aerodynamic design via control theory[END_REF]. However, these developments have not been utilized to their full extent in maritime applications, such as marine propulsors, turbines, hydrofoils and control surfaces. As 90% of the world trade still happens by maritime transportation, the energy-efficiency of maritime transport can be critical to the global sustainability as well as economic growth worldwide.

Over the last decade, significant advances have been made in the design and analysis of marine propulsors. However, most of the tools in the literature either only considered hydrodynamic effects, use low-fidelity methods, or use high-fidelity methods with few design variables (typically less than 15). Researchers have performed hydrodynamic shape optimization of hydrofoils using potential flow methods without cavitation consideration [START_REF] Hsin | Application of the panel method to the design of two-dimensional foil sections[END_REF] and with cavitation consideration [START_REF] Brockett | Minimum pressure envelopes for modified NACA-66 sections with NACA A= 0.8 camber and buships type 1 and type 2 sections[END_REF][START_REF] Eppler | Wing sections for hydrofoils-Part 1: Symmetrical profiles[END_REF][START_REF] Sa Kinnas | Nonlinear analysis of viscous flow around cavitating hydrofoils[END_REF][START_REF] Mishima | A numerical optimization technique applied to the design of two-dimensional cavitating hydrofoil sections[END_REF]. Few researchers also carried out high-fidelity hydrodynamic shape optimization for hydrofoils, naval vehicles, and catamarans [START_REF] Wilson | Hydrodynamic shape optimization for naval vehicles[END_REF][START_REF] Campana | Simulation-based design of fast multihull ships[END_REF][START_REF] Mateusz M Pluciński | Optimization of a self-twisting composite marine propeller using genetic algorithms[END_REF][START_REF] Nelson | Simultaneous optimization of propeller-hull systems to minimize lifetime fuel consumption[END_REF]. They used gradient-free algorithms for optimization, which limited the number of design variables to less than 15 due to the large number of function evaluations compounded with the high computational cost of high-fidelity solvers. The challenge of performing high-fidelity hydrody-namic shape optimization with respect to large number of design variables (¡ 200), with consideration for cavitation, has been tackled in Garg et al. [START_REF] Garg | High-fidelity hydrodynamic shape optimization of a 3-D hydrofoil[END_REF][START_REF] Garg | High-fidelity hydrodynamic shape optimization of a 3-d morphing hydrofoil[END_REF] through the use of gradient-based algorithms together with the adjoint method to efficiently compute the required gradients. In Garg et al. [START_REF] Garg | High-fidelity hydrodynamic shape optimization of a 3-D hydrofoil[END_REF], the authors concluded that an efficient high-fidelity hydrodynamic shape optimization method capable of handling large number of design variables is necessary to get a fully converged optimal design for highly loaded conditions, where viscous effects may not be negligible.

In aircraft wing design, Kenway et al. [START_REF] Kenway | Multipoint high-fidelity aerostructural optimization of a transport aircraft configuration[END_REF] carried out coupled aerostructural optimization of the Common Research Model (CRM) wing, with 472 geometric and structural variables and 57 constraints. They minimized the fuel burn subject to lift, pitch moment, geometric and structural constraints. The optimization minimized the fuel burn by up to 11.2%. While the aerostructural optimization has been successfully applied in aircraft wing design, maritime applications bring additional challenges such as higher loading, stronger fluidstructure interaction, as well as the potential susceptibility to free-surface, cavitation, and hydroelastic instabilities.

Previously, some research has been carried out to optimize the design of composite marine propellers, using potential flow methods to compute hydrodynamic performance and the Finite Element Analysis (FEA) to compute the structural performance [START_REF] Motley | Integrated probabilistic design of marine propulsors to minimize lifetime fuel consumption[END_REF][START_REF] Liu | Utilization of bend-twist coupling for performance enhancement of composite marine propellers[END_REF]. Plucinski et al. [START_REF] Mateusz M Pluciński | Optimization of a self-twisting composite marine propeller using genetic algorithms[END_REF] carried out optimization of self-twisting marine propeller for many inflow conditions, using a coupled Boundary Element Method (BEM)-Finite Element Method (FEM). However, as they used gradient-free genetic algorithm as the optimizer, the number of design variables were restricted to the ply angles and their stacking sequence for the composite propellers. They also did not include any shape design variables. In [START_REF] Motley | Performance-based design and analysis of flexible composite propulsors[END_REF], Motley et al. carried out global performance-based design and analysis for rigid and flexible marine propellers, also using 3-D coupled boundary element method-finite element method (BEM-FEM). Thus, there is a need for an efficient, highfidelity 3-D coupled hydrostructural design optimization tool that can handle the large number of design variables and can enforce constraints on the maximum stress, as well as avoid or delay cavitation.

OBJECTIVES

The objective of this work is to present an efficient highfidelity hydrostructural optimization tool for 3-D lifting surfaces operating in viscous and nearly incompressible fluids, with consideration of cavitation and maximum allowable stress. An unswept, cantilevered, tapered NACA 0009 hydrofoil is presented as a canonical representation of more complex lifting surfaces like propellers, turbines, and other control surfaces.

METHODOLOGY

The tool used for optimization is modified from the Multidisciplinary Design Optimization (MDO) of Aircraft Configurations with High-fidelity (MACH) [START_REF] Gaetan | Scalable parallel approach for highfidelity steady-state aeroelastic analysis and derivative computations[END_REF][START_REF] Kenway | Multipoint high-fidelity aerostructural optimization of a transport aircraft configuration[END_REF]. The MACH framework has the capability of performing static aeroelastic (aerostructural) optimization that consists of aerodynamic shape optimization and structural optimization. In [START_REF] Garg | High-fidelity hydrodynamic shape optimization of a 3-D hydrofoil[END_REF], the MACH framework was extended for hydrodynamic shape optimization of lifting surfaces in viscous and nearly incompressible flows, with consideration for cavitation. Fluids and structures are tightly coupled disciplines in marine propulsor design, as slight changes in the geometry can lead to significant changes in the hydrodynamic performance and internal stress distributions. The focus of this work is to present stateof-art high-fidelity hydrostructural optimization, with large number of design variables. The hydrostructural optimization tool can be divided into following components: CFD solver, structural solver, geometric parametrization, mesh perturbation, load and displacement transfer algorithm, and optimization algorithm.

CFD Solver

The CFD solver used in this paper is Stanford University University of Michigan multi-block Automatic Differentiation solver (SUMad) [START_REF] Van Der Weide | Unsteady turbomachinery computations using massively parallel platforms[END_REF]. SUMad is a finite-volume, cell centered multi-block solver for the compressible flow equations, and is already coupled with an adjoint solver (explained in Section 2.6) for optimization studies [START_REF] Lyu | RANS-based aerodynamic shape optimization of a blended-wing-body aircraft[END_REF]. The Jameson-Schmidt-Turkel [START_REF] Jameson | Numerical solutions of the euler equations by finite volume methods using Runge-Kutta time-stepping schemes[END_REF] scheme (JST) augmented with artificial dissipation is used for spatial discretization. An explicit multistage Runge-Kutta method is used for the temporal discretization. The one-equation Spalart-Allmaras (SA) [START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF] turbulence model is used in conjunction with the adjoint optimizer. In Garg et al. [START_REF] Garg | High-fidelity hydrodynamic shape optimization of a 3-D hydrofoil[END_REF][START_REF] Garg | High-fidelity hydrodynamic shape optimization of a 3-d morphing hydrofoil[END_REF], a low-speed (LS) preconditioner was implemented in SUMad, to solve problems involving nearly incompressible flows for Mach numbers as low as 0.01. The LS SUMad RANS solver was validated against the experimental measurements [START_REF] Zarruk | Experimental study of the steady fluid-structure interaction of flexible hydrofoils[END_REF] and verified against commercial CFD software (ANSYS CFX).

Structural Solver

The structural solver used in the MACH framework is the Toolkit for the Analysis of Composite Structures (TACS) [START_REF] Kennedy | A parallel finiteelement framework for large-scale gradient-based design optimization of high-performance structures[END_REF]. TACS is a parallel, 3-D finite-element solver designed specifically for the analysis of stiffened, thin-walled, composite structures using either static linear or geometrically nonlinear strain relationship. In the current work, only linear structural analysis is used. For linear analysis, the governing equation for structural discipline residual can be written as,

Spu o q Ku o ¡ F (1)
where S are the structural residuals, u o is the vector of structural displacements, K is the linear stiffness matrix, and F is the load vector. Similar to SUMad, the gradients in TACS are computed using an adjoint method, explained in Section 2.6, making the cost of the gradient calculation nearly independent of the number of structural design variables. The derivative of the maximum von-Mises stress is analytically evaluated using the Kreisselmeir-Steinhauser (KS) aggregation functions [START_REF] Wrenn | An indirect method for numerical optimization using the Kreisselmeier-Steinhauser function[END_REF][START_REF] Nicholas | Adaptive constraint aggregation for structural optimization using adjoint sensitivities[END_REF]. In the present research, as the foil is assumed to be made of solid homogeneous isotropic material with linear elastic response. There is a constraint on the maximum allowable stress, but there are no explicit structural design variables. However, the structure will undergo changes in shape as per the shape variables in the hydrodynamic optimization.

Geometric Parametrization

The free-form deformation (FFD) volume approach was used to parametrize the geometry [START_REF] Gaetan | Scalable parallel approach for highfidelity steady-state aeroelastic analysis and derivative computations[END_REF]. To get a more efficient and compact set of geometric design variables, the FFD volume approach parametrizes the geometric changes rather than the geometry itself. All the geometric changes are performed on the outer boundary of the FFD volume. Any modification of this outer boundary can be used to indirectly modify the embedded objects. Figure 4 displays an example of the FFD control points used for optimization of the cantilevered, tapered NACA 0009 hydrofoil, which will be explained in detail in Section 4.

Mesh Perturbation

As the geometry is modified during the optimization using the FFD volume approach, the mesh needs to be perturbed to carry out the CFD analysis and FEA analysis for the modified geometry. The mesh perturbation scheme is a hybridization of algebraic and linear-elasticity methods [START_REF] Gaetan | Scalable parallel approach for highfidelity steady-state aeroelastic analysis and derivative computations[END_REF]. In the hybrid warping scheme, a linear-elasticity-based warping scheme is used for a coarse approximation of the mesh to account for large, low-frequency perturbations; the algebraic warping approach is used to attenuate small, high-frequency perturbations. For the results shown in this paper, the hybrid scheme is not required, and the algebraic scheme is used because only small mesh perturbations were needed to optimize the geometry.

Load and Displacement Transfer Algorithm

The CFD and the structural solvers explained in Section 2.1 and 2.2, respectively, are coupled to form a single hydrostructural solver capable of predicting the coupled hydrodynamic and structural performance. In other words, hydrodynamic loads (pressure and shear stresses) computed by the SUMad (hydrodynamic solver) must be transferred to TACS (structural solver), and displacements from TACS must be used to deform the geometry. Kenway et al. [START_REF] Gaetan | Scalable parallel approach for highfidelity steady-state aeroelastic analysis and derivative computations[END_REF] used the load and transfer scheme from the work of Brown [START_REF] Brown | Displacement extrapolation for CFD + CSM aeroelastic analysis[END_REF]. They constructed rigid links between the CFD surface mesh points, and the points on the structural model lying closest to this set of points. The rigid links are used to extrapolate the displacements from the structural surface to the CFD surface.

The consistent force vector is determined by employing the method of virtual work, ensuring that the force transfer is conservative. The integration of force is performed on the CFD mesh and is transferred to the structure through the rigid links. The two major advantage of using this scheme are that (1) it is consistent and conservative by construction, and (2) it may be used to transfer loads and displacements between CFD and structural meshes that are not coincident, such as the cases with different fluid and structural meshes.

To solve the coupled hydrostructural equations, the hydrodynamic analysis is first partially converged, and the forces are evaluated. These forces are then transferred to the structural analysis, and the corresponding displacements are computed. Thereafter, the displacements are transferred back to the hydrodynamic analysis, the geometry and the corresponding mesh is deformed, and a new CFD solution is found. This iterative loop continues until the coupled convergence criterion (25 iterations or when the relative decrease in both the aerodynamic and structural residuals compared to the free-stream values is less than 1.0E-05) is met.

Optimization Algorithm

The evaluation of the CFD solutions and FEA are the most expensive component of hydrostructural shape optimization algorithms, which can take up to several hours, days or even months. Thus, for large-scale optimization problems, the challenge is to solve the problem to an acceptable level of accuracy with as few CFD or FEA evaluations as possible. There are two broad categories of optimization, namely, gradient-free methods and gradient-based methods. Gradient-free methods, such as genetic algorithms (GAs) and particle swarm optimization (PSO), have a higher probability of getting close to the global minima for problems with multiple local minima. However, gradient-free methods can lead to slower convergence and require larger number of function calls, especially with large number of design variables (on the order of hundreds) [START_REF] Lyu | Benchmarking optimization algorithms for wing aerodynamic design optimization[END_REF]. To reduce the number of function evaluations for cases with large number of design variables, gradient-based optimization algorithm should be used. Efficient gradientbased optimization requires accurate and efficient gradient calculations. There are some straight forward algorithms like finite difference; they are neither accurate nor efficient [START_REF] Martins | The complex-step derivative approximation[END_REF].

Finite difference methods are sensitive to the step-size. If the step-size is higher or lower than the optimal step-size, the errors in FD will increase, as noted by [START_REF] Martins | The complex-step derivative approximation[END_REF]. Complex-step methods yields accurate gradients, but are not efficient for large-scale optimization [START_REF] Martins | The complex-step derivative approximation[END_REF]. Thus, for gradient calculations, the adjoint method is used in this paper. The adjoint method is efficient as well as accurate, but is relatively more challenging to implement [START_REF] Lyu | RANS-based aerodynamic shape optimization of a blended-wing-body aircraft[END_REF].

The optimization algorithm used in this paper is called SNOPT (sparse nonlinear optimizer) [START_REF] Gill | SNOPT: An SQP algorithm for large-scale constrained optimization[END_REF]. SNOPT is a gradientbased optimizer that utilizes a sequential quadratic programming method. It is capable of solving large-scale nonlinear optimization problems with thousands of constraints and design variables.

MODEL SETUP AND VALIDATION 3.1 Fluid Model

For all the results presented in this paper, a cantilevered, unswept, tapered NACA 0009 hydrofoil was studied with an aspect ratio of 3.33 (a span length of 0.3 m and a mean chord length of 0.09 m, with maximum chord length of 0.12 m at the root and a minimum chord length of 0.06 m at the tip), Reynolds Number, Re 1.0 ¢ 10 6 , and V 12.4 m/s. In the present work, Reynolds number (Re V c{ν) is defined based on the mid-span chord length, i.e. 0.09 m. The CFD mesh used for all the optimization and validation results is shown in Figure 1. The CFD mesh used is a structured O-grid with 515,520 cells and a y of 1.1. The domain size is 30 chord lengths in all the directions. The tip of the hydrofoil is modeled rounded, which helps in making the mesh continuous and smooth. The applied boundary conditions (B.C.) are shown in blue color in Figure 1.

Structural Model

To carry out the validation study of the hydrostructural solver and the coupled hydrostructural optimization, the cantilevered, tapered 3-D NACA 0009 hydrofoil (model shown in Section 3.1) was studied. Figure 2 shows the structural mesh used for the computation. The mesh has 44,800 second order isoparametric 8-noded elements, with 7 elements along the thickness direction, 40 elements in the chordwise direction, and 160 elements in the spanwise direction. A perfectly clamped boundary condition is set at the root, i.e., the displacement and rotation is zero at the root.

Eigenvalue Validation and Verification

To validate the structural model, a 3-D cantilevered, tapered NACA 0009 hydrofoil (as shown in Section 3.1, 3.2), with Re 1.0 ¢ 10 6 and V 12.4 m/s was studied. The predictions were compared with experimental measurements conducted at the Cavitation Research Laboratory (CRL) variable pressure water tunnel at the University of Tasmania [START_REF] Doolan | Hydroacoustic characterisation of the amc cavitation tunnel[END_REF]. The operating velocity and pressure range in the tunnel was of 2 -12 m/s and 4 -400 kPa, respectively. The tunnel test section is 0.6 m square by 2.6 m long. They tested four foils of similar geometry but with different materials, namely, SS (stainless steel-316L), Al (Aluminum-6061T6), CFRP-00, and CFRP-30 (CFRP are Carbon Fiber Reinforced Polymer composites with the number accompanying them denoting the alignment of unidirectional fibers). The geometry dimensions were selected such that confinement effects are negligible. They reported estimated measurement uncertainty of 0.5% in the force measurement, of 0.001 o in the angle of attack, and the maximum uncertainty of 13% for the tip bending deflections.

Table 1 shows the parameters used in the experiment. All the parameters (geometry and its dimensions, Reynolds number, fluid density, velocity, material of the hydrofoil) from the experiments were matched in the numerical model, except for the Mach number (defined as M V {a, where V is the inflow velocity and a is the speed of sound). However, as the compressibility effects are almost negligible for Mach number less than 0.1, this discrepancy in Mach number (0.008 in the experiments and 0.05 in the numerical model) should not affect the solution. Poisson's ratio 0.33 σ s,f (MPa) Fatigue strength 96.5

In Table 3, comparison of the predicted and measured first mode frequency of the Al foil in air are presented. The predictions from structural solver, TACS, were validated against the experimental measurements [1] and verified against predictions from the commercial Finite Element Analysis (FEA) software (NASTRAN). TACS was found to under-predict the first mode in-air frequency by less than 4.0%, when compared to the experimental measurements, and under-predict by less than 0.5%, when compared to the commercial FEA software.

Table 3.

Comparison of first mode in-air frequency from the structural solver, TACS, for an unswept, cantilevered, tapered NACA 0009 hydrofoil with the experimental results from [START_REF] Zarruk | Experimental study of the steady fluid-structure interaction of flexible hydrofoils[END_REF] and the commercial FEA software (NASTRAN). The structural mesh is shown in Figure 2.

Material Experiment [START_REF] Zarruk | Experimental study of the steady fluid-structure interaction of flexible hydrofoils[END_REF] Nastran TACS Al (Hz) 100.0 98.5 98.4

Hydrostructural Validation

To validate the coupled hydrostructural solver, a 3-D cantilevered, tapered NACA 0009 hydrofoil (as shown in Section 3.1, 3.2), with Re 1.0 ¢ 10 6 and V 12.4 m/s was studied. In Figure 3, the predicted lift force (in Newtons) and the predicted tip bending deformation (in mm) at various angles of attack from modified MACH framework was compared with the experimental measurements from [START_REF] Zarruk | Experimental study of the steady fluid-structure interaction of flexible hydrofoils[END_REF]. The predictions were found to be within the experimental error (0.5% for the lift force and 13% for the bending deformation).

For Re 1.0 ¢ 10 6 , experiments were only conducted until a maximum angle of attack of 6 o , to avoid excessive forces on the foil. the predicted tip bending deformation (in mm) at various angle of attack (α) obtained using the modified MACH framework with the experimental results from [START_REF] Zarruk | Experimental study of the steady fluid-structure interaction of flexible hydrofoils[END_REF] for a cantilevered, tapered Al NACA 0009 hydrofoil. The CFD mesh with 515,520 cells (y 1.1) for MACH framework is shown in Figure 1 and the structural mesh with 44,800 8-noded elements is shown in Figure 2. The conditions selected corresponds to Re 1.0 ¢ 10 6 and V 12.4 m/s, as reported in Table 1.

Grid Convergence Study

To ensure that the structural results are independent of the mesh size, the grid convergence was studied with four different mesh sizes: 25,200 cells, 44,800 cells, 100,800 cells, and 179,200 cells for 3-D cantilevered, tapered NACA 0009 hydrofoil (as shown in Section 3.1), with Re 1.0 ¢ 10 6 , V 12.4 m/s, and α 6 o . To be consistent, a medium-sized CFD mesh of 515,520 cells was used for all the analyses. Grid convergence study of the CFD mesh has already been presented in [START_REF] Garg | High-fidelity hydrodynamic shape optimization of a 3-D hydrofoil[END_REF]. Table 4 shows the comparison of lift coefficient (defined as

C L Lift{0.5ρ f V 2 A), lift force (N),
first mode in-air frequency (Hz), and tip bending deformation (mm) values from coupled hydrostructural simulation. Based on the results and the computational time, the medium-sized structural mesh with 44,800 elements was selected to carry out the hydrostructural analysis and optimization.

Optimization Problem Formulation

To demonstrate the advantages of hydrostructural shape optimization, results are shown for the 3-D cantilevered, tapered NACA 0009 hydrofoil (as shown in section 3.1, 3.2), with Re 1.0 ¢ 10 6 and V 12.4 m/s. The optimization problem setup is described in Table 5. The drag coefficient, defined as C D Drag force{0.5ρ f V 2 A, is minimized for a given C L , a given cavitation number, defined as σ pP ref ¡ P vap q{0.5ρ f V 2 , and the constraint on maximum von-Mises stress (σ v ) of the optimized shape against the fatigue strength (σ f ) of Al with a factor of safety of 1.2. Constraint on the minimum volume and minimum thickness are also detailed in Table 5. C ¦ L is the target C L , d is the volume of the optimized foil, d base is the volume of the original foil, t base is the thickness of the original foil at a given section. The leading edge of the hydrofoil is also constrained to deform to prevent any random behavior at the leading edge.

The influence of the number of design variables was studied in [START_REF] Garg | High-fidelity hydrodynamic shape optimization of a 3-D hydrofoil[END_REF], where the authors concluded at least 200 FFD control points are necessary to achieve a converged optimal hydrodynamic design for the same unswept cantilevered, tapered hydrofoil. Thus, for the results shown in Section 4, a total of 210 design variables were used with 200 FFD control points (10 spanwise ¢ 10 chordwise ¢ 2 thickness) and 10 spanwise twist design variables.

Figure 4 shows a pictorial description of the design variables with 200 FFD control points and 10 twist design variables on the right side. On the left side in the Figure 4, the original NACA 0009 hydrofoil with the deflected shape at C L 0.65 is displayed.

RESULTS AND DISCUSSION

In this section, the single-point RANS-based coupled hydrostructural optimization results are presented for a cantilevered, unswept, tapered aluminum (6061 T6) NACA 0009 hydrofoil at Re 1.0 ¢ 10 6 and V 12.4 m/s. The foil was optimized to achieve the minimum drag coefficient (C D ) for a target lift coefficient (C L ), a cavitation number (σ) of 1.6, and a constraint on the material stress against the fatigue strength of Al with a factor of safety of 1.2. The optimization was carried out with 210 design variables. 4. Comparison of CL, lift force (N), first mode in-air frequency (Hz), and tip bending deformation (mm) values from the coupled hydrostructural simulation for the cantilevered, tapered Al NACA 0009 hydrofoil at Re 1.0 ¢ 10 6 , V 12.4 m/s, and α 6 o using different structural mesh with the modified MACH framework. The same CFD mesh of 515,520 cells was used for all the analyses.

Mesh

Mesh

¡C p σ Cavitation number constraint 1 σ v ¤ σ f {1.2
Stress constraint 1 the design point. In Figure 6, the predicted non-dimensional bending deformation, i. To meet the stress constraint (σ v ¤ σ f {1.2), the optimal foil needed to be thicker near the root, where the stress is the highest, and is balanced by the thinner sections at the tip, as shown in the sectional geometry in Figure 8d. In Figure 8b, the difference between the stress constraint contours, i.e. 1.2 ¢ σ v {σ f of Al, of the original NACA 0009 hydrofoil and the coupled hydrostructural optimization is presented. The red region, indicating where the stress constraint is violated, can be noticed near the root of the original NACA 0009 hydrofoil in the Figure 8b. While meeting the maximum stress constraint, the mass of the optimized foil is also reduced by 7.3%. The mass reduction while satisfying the structural constraints is particularly impressive considering the NACA 0009 hydrofoil is a very thin and efficient foil to start with. The optimized design has a significant camber resulting in higher moment of inertia, which makes the structure more rigid. Thus, the optimized foil has much higher stiffness and lower stress, compared to the original NACA 0009 hydrofoil. Figure 8a indicate that partial leading edge cavitation (as indicated by gray contour region with ¡C p ¥ σ) will develop around the original NACA 0009 hydrofoil at C L 0.65 and σ 1.6, but no cavitation is observed for the optimized foil. The maximum negative pressure coefficient, ¡C p , reduces from 2.1 for the NACA 0009 hydrofoil to 1.0 for the optimized foil, as shown in Figure 8d, which will help to significantly delay cavitation inception. In other words, the cavitation inception speed for the coupled hydrostructural optimized foil will increase from 10.2 m/s to 14.8 m/s, for an assumed submergence depth of 1 m. Note that the performance of the optimized foil exceeds that of the NACA 0009 not only at the design point (CL 0.65), but also for all lift coefficients above 0.25. All the results were obtained using coupled hydrostructural solver for solid Aluminum (Al) hydrofoil at Re 1.0 ¢ 10 6 and V 12.4 m/s.

CONCLUSION

We presented an efficient high-fidelity coupled hydrostructural optimization capable of handling a large number of design variables, with constraints on cavitation and maximum stress. The hydrostructural solver was validated against experimental measurements [START_REF] Zarruk | Experimental study of the steady fluid-structure interaction of flexible hydrofoils[END_REF] for the hydrodynamic load coefficients, tip bending deformations, and in-air natural frequencies for the case of a tapered, cantilevered, solid aluminum NACA 0009 hydrofoil. Good agreements were observed between the numerical predictions and experimental measurements. The maximum difference in the force predictions is 3.5% and in the tip bending deflections is 7.0%.

To provide a canonical representation of a complex hydrodynamic lifting surface,we performed a RANS-based hydrostructural optimization using the adjoint method for the same NACA 0009 hydrofoil at Re 1.0 ¢10 6 and V 12.4 m/s.

The single-point high-fidelity hydrostructural optimization was performed at C L 0.65 with 210 design variables. The optimized hydrofoil resulted in increase in C L {C D of 12.4%, a mass reduction of 7.3%, and an increase in the cavitation inception speed by 45%, when compared to the original NACA 0009 hydrofoil. While the original NACA 0009 hydrofoil fails the maximum stress constraint (i.e., σ v ¡ σ f {1.2) for lift coefficients greater than 0.35, the optimized foil was able to meet the stress constraint up to the design lift coefficient of 0.65. A systematic study of the design space of marine propulsors using the present state-of-art high-fidelity coupled hydrostructural optimization has the potential to drastically improve propulsors's performance. Improved propulsor performance can include improved fuel efficiency and hence a reduction in the operational cost and reduced negative environmental impact, enhanced agility, and reduced structural weight, while ensuring the structural integrity and delaying cavitation inception.

FUTURE WORK

The paper presents a state-of-art high-fidelity coupled hydrostructural design optimization framework, capable of handling large number of design variables with consideration for cavitation and maximum stress. While the optimization tool is demonstrated for hydrofoils, it can be easily extended to other lifting bodies, such as marine propellers and turbines, and to optimize the material properties and stacking sequence of composite propulsors. As shown previously by the authors [START_REF] Garg | High-fidelity hydrodynamic shape optimization of a 3-D hydrofoil[END_REF], multipoint optimization results in a much more efficient design over the wide range of operating conditions, when compared to the single-point optimized design. Thus, coupled hydrostructural multipoint optimization will be carried out to take full advantage of this framework.

Figure 1 .

 1 Figure 1. CFD Mesh (515,520 cells) used for the RANS optimization of a cantilevered, tapered NACA 0009 hydrofoil at Re 1 ¢ 10 6 . a) Front view with the boxed in portion showing the location of the foil. b) Side view with boxed in portion showing the location of the foil. c) The zoomed-in view of the foil, showing the mesh near the foil leading edge (LE) and trailing edge (TE). d) Geometry with the dimensions of the tapered NACA 0009 hydrofoil.

Figure 2 .

 2 Figure 2. Mesh (44,800 elements) used for structural validation as well as hydrostructural optimization. The mesh has 44,800 second order isoparametric 8-noded elements, with 7 elements along the thickness direction, 40 elements in the chordwise direction, and 160 elements in the spanwise direction. A perfectly clamped boundary condition is set at the root, i.e., the displacement and rotation is zero at the root. The red curve depicts the edges of the fluid model.

Figure 3 .

 3 Figure 3. Comparison between the predicted lift force (in N) and

Figure 4 .

 4 Figure 4. Pictorial description of the problem setup used for the optimization study. Original cantilevered Al NACA 0009 hydrofoil with the initial and the deflected shape is displayed on the left side and the shape design variables; 200 FFD control points (10 spanwise ¢ 10 chordwise ¢ 2 thickness, as indicated by the red circles), and twist design variables on the right side.

Figure 5

 5 depicts the difference in 3-D geometry between the cantilevered, tapered NACA 0009 hydrofoil and the singlepoint hydrostructural optimized foil at C L 0.65. In Figure 6, the stress constraint, i.e., σ v ¤ σ f {1.2, of the original NACA 0009 hydrofoil and the hydrostructural optimized foil at C L 0.65, at various C L values is presented. The black dotted line represents the stress constraint, i.e. σ v σ f {1.2. While the original tapered NACA 0009 hydrofoil violates the stress constraint for C L ¡ 0.35, the hydrostructural optimized foil satisfies the stress constraint for C L values up to 0.65, Table

  e., the fraction of the maximum tip bending deformation to the tip chord length (0.06 m), of the NACA 0009 hydrofoil and the hydrostructural optimized foil at C L 0.65, at various C L values is also presented. While the original tapered NACA 0009 hydrofoil has the maximum non-dimensional bending deformation of 0.20, the optimized foil at C L 0.65 has the maximum non-dimensional bending deformation of 0.16.

Figure 5 .

 5 Figure 5.Figure showing the comparison of 3-D geometry

  Figure 5.Figure showing the comparison of 3-D geometry between the original tapered NACA 0009 hydrofoil and the hydrostructural optimized foil at CL 0.65. The computations were carried out for solid Aluminum (Al) hydrofoil at Re 1.0 ¢ 10 6 and V 12.4 m/s.

Figure 7 Figure 6 .

 76 Figure 7 shows the comparison of efficiency (i.e., C L {C D ) at various C L values for the NACA 0009 hydrofoil and the single-point optimized foil at C L 0.65. As can be observed from the figure, the hydrostructural optimized foil performs better than the original NACA 0009 foil for C L ¥ 0.3. Below C L value of 0.3, the optimized foil performance becomes

Figure 8

 8 Figure 8 shows a comparison between the NACA 0009 hydrofoil and the hydrostructural optimized hydrofoil at C L

Figure 7 .

 7 Figure 7.Figure showing the comparison of CL{CD versus CL for the original tapered NACA 0009 hydrofoil and the hydrostructural optimized foil at CL 0.65, at various CL values.

Figure 8 .

 8 Figure 8. Figure showing comparison of the performance of the original NACA 0009 hydrofoil and the single-point hydrostructural optimized hydrofoil at CL 0.65, Re 1.0 ¢ 10 6 , and V 12.4 m/s. A reduction in CD of 11.3% is noted for the optimized foil. Going from top to bottom and from left to right. a) Cp (pressure coefficient) contours plot on the suction side are displayed for NACA 0009 hydrofoil and the optimized foil. Gray contour region along the leading edge of the tapered NACA 0009 foil displays the area with ¡Cp ¥ σ. b) Comparative study of the stress coefficient contours (1.2 ¢ σv{σ f ) of the NACA 0009 hydrofoil and hydrostructural optimization is presented. Excessive stress (σv ¡ σ f {1.2) region can be noticed near the root for the original NACA 0009 hydrofoil. Reduction in mass of 7.3% for the optimized foil can be noted. c) Displays the skin friction (C f ) contour plot for the original tapered NACA 0009 hydrofoil and the optimized foil at CL 0.65. In general, the optimized foil has a lower skin friction, and hence lower drag coefficient. d) Figures show the sectional ¡Cp plots and the geometry profile of the foil at 3 sections along the span of the hydrofoil. Blue solid line represents the NACA

Table 2

 2 displays the material properties of the Aluminum (Al) foil, which is used for validation and optimization studies in this work.

Table 1 .

 1 Parameters used for an unswept, cantilevered, tapered NACA 0009 hydrofoil in the experimental results from[START_REF] Zarruk | Experimental study of the steady fluid-structure interaction of flexible hydrofoils[END_REF].

	Parameter	Experiment [1]
	Geometry	NACA 0009
	Aspect ratio Reynolds number ρ f (fluid density) (Kg/m 3 ) 999 3.33 1 ¢ 10 6
	Mach number	0.008
	V f (m/s)	12.4
	Material	Aluminum (6061 T6)
	Fluid	Water

Table 2 .

 2 Structural

		Description	Al
	ρ s (Kg/m 3 ) Solid density	2700
	E s (GPa)	Elastic modulus 69
	ν s

properties of Aluminum(Al) used in the present research for validation and the optimization study, which corresponds to the values reported for the Al foil in

[START_REF] Zarruk | Experimental study of the steady fluid-structure interaction of flexible hydrofoils[END_REF]

.

Table 5 .

 5 Optimization problem for the cantilevered, tapered Al NACA 0009 hydrofoil (shown in Section 3.1, 3.2). ¢ t base ¤ t i ¤ 1.5 ¢ t base Minimum thickness constraint 400 d ¥ d base

		Size C L	Lift Force (N) Frequency (Hz) Tip bending deformation (mm)
	Coarse	25,200	0.4980 1034.8	99.53	12.46
	Medium 44,800	0.4987 1035.2	98.43	13.02
	Fine	100,800	0.4990 1035.9	97.58	13.06
	Finest	179,200	0.4991 1036.1	97.26	13.36
			Function variables		Description	Qty.
	minimize	C D			Drag coefficient	1
	Design variables f			FFD control points	200
						Twist design variables	10
	Constraint	C ¦ L			Lift coefficient constraint	1
			0.5 Minimum volume constraint	1
						Fixed leading edge constraint 10

ACKNOWLEDGMENTS

Support for this research was provided by the U.S. Office of Naval Research (Contract N00014-13-1-0763), managed by Ms. Kelly Cooper. The authors would like to thank Prof. Paul Brandner of the Australian Maritime College for providing the experimental data and answering questions about the experimental setup. 0009 hydrofoil and the red solid line corresponds to the hydrostructural optimized foil. Grey horizontal line represents the constraint on cavitation number (when ¡Cp σ 1.6). As noted, there is significant decrease in maximum ¡Cp from the NACA 0009 to the optimized hydrofoil. Difference in the sectional shape as fraction of the span (z{s) between the original and optimized foil are also shown in the bottom of each subplot.


Laveena Sharma

Dhiman Chatterjee
email: dhiman@iitm.ac.in

Numerical prediction and optimization of the performance of axial-flow hydrokinetic turbine in an array

Keywords:

Array of hydrokinetic turbines can be used to produce power from flowing rivers akin to the wind turbine farms. However, arriving at optimum configuration (location and number) of turbines through experimental or numerical routes is not easy. Hence in the present work, attempts have been made to utilize computationally less intensive model-based multivariate optimization techniques coupled with some numerical simulations to arrive at suitable turbine configurations in an array. In this work, we have modeled the complex functions, such as turbine power, non-linear interaction of turbines with wakes of preceding turbines and the area utilized by the turbines, by means of simplified mathematical treatment of response surface methodology. The present work has identified three turbine array geometries as potential configurations and based on numerical simulations and optimization techniques has arrived at pareto-domain and has ranked these pareto-solutions. Further simulations were carried out for the best configurations and the numerical results were compared with the optimization results as a part of the validation of the latter. Flow physics was analyzed and final recommendation of turbine array configuration is made.

INTRODUCTION

Hydrokinetic turbine utilizes kinetic energy of flowing water from rivers and seas and converts into this energy into useful mechanical power. The efficiency of this type of "net-zero head" turbine is lower than that of a conventional high head hydraulic turbine; but it is more suited from the perspective of minimum adverse effect on environment and population displacement. This has led to a growth in activities in this research topic in recent times. [START_REF] Khan | Hydrokinetic energy conversion systems and assessment of horizontal and vertical axis turbines for river and tidal applications: a technology status review[END_REF]. Though this technology is an emerging alternative, yet it suffers from low energy conversion coefficient. Thus, akin to wind farms, it may be more desirable to install an array of these turbines rather than a single turbine. However, this nascent technology faces uncertainties and challenges not only for extracting more energy per unit of rotor swept area, but also to extract the maximum energy when such turbines are coupled together at a particular site. One of the reasons that the total power produced is simply not the product of power output from one isolated turbine and the number of turbines is because the turbines which are in the wakes of upstream turbines suffer from loss of power. Thus, efforts have been made to arrange turbines in arrays (rows and columns) or in discrete, optimum locations [START_REF] Chowdhury | Optimizing the arrangement and the selection of turbines for wind farms subject to varying wind conditions[END_REF]. However, it is not easy to arrive at an optimum configuration through experimental or direct numerical simulations of turbine arrays [START_REF] Funke | Tidal turbine array optimsation using the adjoint approach[END_REF]. This paper attempts to use computational fluid dynamics techniques coupled with rigorous optimization algorithm to arrive at optimum turbine array configurations of unidirectional, axial flow hydrokinetic turbine.

NUMERICAL METHODOLOGY

Computational fluid dynamics simulations were carried out for isolated (2.9 million elements) as well as unit (described in the next section) of array of turbines (12.7 million elements) using commercial software, ANSYS. A typical computational domain with boundary conditions is shown in Fig. 1. Performance of isolated turbine was validated against experimental results and the same is reported in one of our earlier works and it was found that there is a deviation of less than 10% between experimental maximum power output and that of numerical, both at a tip speed ratio of 3.5 [START_REF] Chandras | Numerical prediction of the performance of axial flow hydrokinetic turbine[END_REF][START_REF] Chatterjee | Design and development of hydrokinetic turbine suited for Indian condition[END_REF]. It may be mentioned here that the turbine designed and numerically simulated had a hub diameter of 300 mm, a tip diameter of 1600 mm and design speed of 40 rpm. These results as well as numerical methodology are not reproduced here for brevity.

STEPS ADOPTED FOR OPTIMIZATION OF TURBINE ARRAY

Figure 2 gives a summary of the steps adopted in arriving at an optimum turbine array for a given choice of configuration (Fig. 3). The first step involves setting up the objective criteria and is defined as:

𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀 𝑃𝑃𝑀𝑀𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑀𝑀𝑀𝑀𝑃𝑃𝑃𝑃𝑀𝑀 = ∅ � 𝑀𝑀𝑀𝑀𝑀𝑀. 𝑃𝑃 𝑑𝑑 (𝑀𝑀 1 , 𝑀𝑀 2 ,
𝑃𝑃 2 𝑃𝑃 1 = � 𝑉𝑉 2 𝑉𝑉 1 � 3 � 𝑃𝑃 1 , 𝑉𝑉 1 -𝑅𝑅𝑃𝑃𝑇𝑇𝑃𝑃𝑃𝑃 1 𝑃𝑃 2 , 𝑉𝑉 2 -𝑅𝑅𝑃𝑃𝑇𝑇𝑃𝑃𝑃𝑃 2
ub denotes the upper limit of spacing between two turbines such that full power (same as that of the preceding turbine) is obtained from a succeeding turbine. lb is the lower limit of spacing such the succeeding turbine gets completely immersed in the wake of the preceding turbine [START_REF] Frandsen | The necessary distance between large wind farms offshore[END_REF]. V 2 is obtained numerically at the location of the second rotor. Power density (P d ) is defined as the power output per unit area. Power loss (P l ) for the second turbine due to the wake created by the first one is defined as the difference in the ideal power (P 2 ) predicted on the basis of the local velocity at the second turbine location and the actual power (P 2 ') obtained from the turbine. When the incoming flow with a velocity of 1m/s encounters first rotor it loses its kinetic energy which is regained at the downstream of the first rotor. The velocity follows the trend as shown in Fig. 9. Hence, at a distance of 2.5 m from the rotor, the flow should regain the velocity of 0.75 m/s. If a second rotor is placed at this location, it should deliver a power of 0.42 P 1 , but in actual practice the power delivered by the second rotor is less than 0.42 P 1 . This loss is accounted as the power loss due to the wake. Area utilization ratio (A u ) is defined as the ratio of the total area actually occupied by the turbine units in a pattern to the sum of given array area (50 m X 50 m) and virtual area created to avoid degeneracy (area reducing to line at the inline edge when breadth becomes zero). While arranging the rotor units in the above geometric patterns of varying dimensions, and varying spacing between them, it has been found out that it is not possible to accommodate the complete pattern in few cases. For these cases, in order to utilize the space as much as possible, a virtual row of turbines is added and a virtual area is created (Fig. 4). The role of these virtual turbines is to avoid degeneracy during the calculation of the area occupied by the rotor patterns. These virtual units do not contribute to the power density.

Fig. 4. Schematic showing concept of virtual area and spacing between pattern for turbine arrays

A meta-model of the significant functions is developed based on the data obtained through numerical simulations. Random samples (space filling designs, equidistant grids for arrays 1 and 2 and nested designs for Array 3) are taken to capture the feasible domain effectively [START_REF] Jurecka | Robust design optimization based on metamodeling techniques[END_REF]. The model is then validated using Predicted Residual Square Sum of Errors (PRESS) that accounts for the predictive capability of the model. Each observation is removed and the model is refitted using the remaining observations. That refitted model is then used to predict the left-out-sample. This strategy is repeated for each observation and the square sum of errors is then computed, known as PRESS. The model with a least PRESS is then adopted.

The problem is unconstrained and contains a feasible set which is closed and bounded, with the functions being continuous within those bounds. Genetic algorithm is used to circumscribe the pareto-domain which contains solutions (Sol) that are better for at least one of the objectives than all the other solutions that are the also part of the domain. For ranking the obtained pareto-domain, Net Flow Method (NFM) [START_REF] Rangaiah | Multi-Objective Optimization Techniques and Applications in Chemical Engineering[END_REF] is used. The relative weight, and the three parameters, an indifference threshold (Q), a preference threshold (P) and a veto threshold (V) for each function is set according to the understanding of the physical problem. Since the prime objective of a turbine array is to maximize the power density, function f 1 is given a preference over the second and the third functions defined in the objective criteria. NFM suggests that instead of arriving at one best solution corresponding to the highest rank, it is advisable to categorize the solutions into different categories [START_REF] Rangaiah | Multi-Objective Optimization Techniques and Applications in Chemical Engineering[END_REF]. Here, looking at the solution, a decision was made to categorize the solutions into two categories. The solutions with a negative score are eliminated and thus categorized under non-preferred set (NP), whereas the positive score solutions are incorporated under preferred set (P). Overall recommendation is given based on the results obtained for different simulations.

RESULTS AND DISCUSSION

In this section, we shall discuss the performance of turbine arrays for each of the three patterns shown in Fig. 3.

Results of Array 1

Array 1 (Fig. 3a) is essentially an in-line arrangement of turbines. Thus it is expected that the position of the second turbine is very crucial. Figure 5 shows the variation of velocity deficit calculated in the wake of a turbine where U is taken as velocity of the flow after it is disturbed by rotor 1, and U ∞ is the incoming water velocity of 1m/s. This wake has a significant effect on power loss in the second turbine [START_REF] Milborrow | The performance of arrays of wind turbines[END_REF]. This power loss, calculated using the relationship mentioned earlier.

The profile of velocity approaching the second rotor is not uniform. In this work, instead of choosing an average velocity, the velocity is computed at the midspan in the flow passage between two rotor blades. This choice of velocity could have resulted in the peak loss observed in Fig. 6. Inset figure brings out the extent of the wake with the distance between two in-line rotors. With a distance of 2.5 m rotor 2 completely lies in the wake regime of rotor 1, resulting in poor power ratio, whereas as the distance increases rotor 2 has no interaction with the wake regime of the previous rotor, thereby resulting in an efficient power ratio as can be seen at the distance of 8.5 m. As observed from Fig. 7, power ratio i.e. P 2 /P 1 is not sensitive to lateral spacing (x) as compared to in-line spacing (z), where P 1 and P 2 are the power extracted by rotor 1 and rotor 2 respectively. Hence, for Array 1, lateral variation is not considered as a significant factor.

Results of Array 2

Array 2 has a pattern of equilateral triangle of side a and arranged in such a way that each of the 3 turbines face the incoming flow with minimum disturbance from others (Fig. 3b) . This fact is borne by the fact that incoming velocities faced by each turbine is same (Fig. 9) and hence the power produced by these turbines is similar. Thus, for this configuration, the power loss term is almost negligible.

Figure 10 shows the variation of normalized power density and normalized area utilization ratio for Array 2 as functions of the side of equilateral triangle.

𝑃𝑃𝑃𝑃𝑃𝑃𝑀𝑀. 𝑃𝑃𝑢𝑢𝑃𝑃𝑃𝑃𝑇𝑇𝑀𝑀𝑃𝑃𝑃𝑃 𝑃𝑃 𝑃𝑃 = 𝑃𝑃 -𝑃𝑃 min 𝑃𝑃 max -𝑃𝑃 min Points indicate velocities at rotor locations.

It has been observed that power density reduces with increase in length, yet area utilization increases. Hence a pareto-domain gives a solution an acceptable set of feasible solutions. In order to achieve this, a meta-modeling is done using three out of four models -linear, quadratic, power, and exponential. It can be observed here that the models interpolate the values in a close proximity (Fig. 11).

On the basis of the coefficient of determination (R 2 ) value for different models, the ones with highest value are chosen and final selection is made on the basis of PRESS, a statistical measure used for cross validation of model. Figure 12 shows the pareto-domain for the above two functions generated using genetic algorithm in Matlab®. The solutions are found to get converged after 133 generations. Table 2 gives details of parameters used to rank pareto-domains based on NFM. The ranking using NFM is shown in Table 3. This ranking is indicative of the scores obtained as shown in Table 3. It may be added here that these rankings given in Table 3 were arrived at after studying the sensitivity of different weights and thresholds as stated earlier.

Results of Array 3

A unit cell of Array 3 (Fig. 3c) comprises of turbines arranged in the form of an isosceles triangle such that there is a alternating row of turbines in line with the previous one. Thus, it is clear that both in-line and lateral spacing play important roles for such an array. This was verified in our work. Thus for this array, x-and z-separations and all three functions (P d , P l and A u ) become important.

Figure 13 shows the variation of power as functions of inline and lateral spacing. From this figure it is seen that though turbine 2 does not differ significantly with x-or zseparations, turbine 3 does depend on the separation distances. Figure 14 shows the predictive errors for the three functions (power density, area utilization ratio and power loss terms) when calculated using different polynomial functions. The solid line in each of these plots (Fig. 14 a-c) indicates the one function which has been chosen for each of these cases. It may further be mentioned here that though area utilization ratio is a convex function, the other two functions are concave and so we have used genetic algorithm [START_REF] Arora | Introduction to Optimum Design[END_REF]. The fitted curves for each of the functions are given as:

Fig. 14. Predictive capability of four polynomial metamodels used for each of the three functions

Figure 15 shows response surface and contour plots of different metamodels used. Figure 16 shows the pareto-front for Array 3 obtained using genetic algorithm and this result confirms that non-convex part is also captured using the genetic algorithm approach in 169 generations. It shows the conflicting relationship between the three functions. The 3d points represent the relationship between all the three functions, whereas 2d points projected on the planes gives the relationship between the functions taken two at a time. It can be observed that when area utilization increases power density decreases and power loss increases. The units of the functions are mentioned in the earlier section.

The rankings of pareto-solutions were arrived at based on parameters shown in Table 4. Table 5 gives the preferred or non-preferred rankings for different pareto-optimal solutions for Array 3 similar to what was done for Array 2. Based on the optimizations carried out for the three arrays, ranking 1 (in case of Arrays 2 and 3) and upper bound in case of Array 1 has been selected for validation studies. In order to compare the performance of these optimized turbines with that of any configuration, an initial configuration has been chosen. This initial configuration has been chosen such that there is no spacing between the unit patterns as shown in Fig. 3(a-c). These comparisons are reported in the next section on validation.

Validation and comparison of optimized solutions

Detailed CFD simulations were carried out for both optimized solutions and general arrays and Table 6 gives a validation of optimization results against numerical simulations. Thus it is seen that optimization results are in reasonably good agreement with that obtained from numerical simulations. These flow interactions are crucial in achieving the improved performance of optimized configurations as shown in Table 7. In particular attention is drawn to the fact that for each of the cases depicted in Table 7, different numbers of turbines are used and hence from practical point of view of deployment of turbines, it is ideal if fewer number of turbines can produce similar or higher power outputs. It is noted that of all the three configurations chosen, Array 2 produces the maximum output. Array 3, deploying fewer number of turbines, is capable of producing similar power output as Array 2 with similar number of turbines but would require larger space and hence may be problematic for other users of the water passage. Thus, it can be concluded that the optimization approach to arrive at suitable turbine arrays is validated and such an approach can improve the performance of the hydrokinetic power production.

CONCLUSIONS

Hydrokinetic turbines are less efficient and hence in order to produce significant power arrays are essential. The power output from these turbines as well as the area occupied by the turbines determine the successful installation and operation of these turbines. However, carrying out full-fledged computational and/or experimental work is difficult, time consuming and expensive. Hence, in this work, we have used multivariate optimization aided by selective numerical simulations. Significant conclusions are listed below: 1. It is observed here that all the geometries do not respond to the functions in a similar way. For instance for Array 1 geometry, variation of power loss and area utilization with in-line spacing were of similar nature while for Array 3 they were conflicting. Thus, i.

𝐴𝐴𝑃𝑃𝑃𝑃𝑀𝑀𝐴𝐴 1 𝑃𝑃𝑀𝑀𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑀𝑀𝑀𝑀𝑃𝑃𝑃𝑃𝑀𝑀 = ∅{𝑀𝑀𝑀𝑀𝑃𝑃. 𝑃𝑃 𝑙𝑙 (𝑀𝑀), 𝑀𝑀𝑀𝑀𝑀𝑀. 𝐴𝐴 𝑢𝑢 (𝑀𝑀)}𝐴𝐴 ii.

𝐴𝐴𝑃𝑃𝑃𝑃𝑀𝑀𝐴𝐴 2 𝑃𝑃𝑀𝑀𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑀𝑀𝑀𝑀𝑃𝑃𝑃𝑃𝑀𝑀 = ∅{𝑀𝑀𝑀𝑀𝑀𝑀. 𝑃𝑃 𝑑𝑑 (𝑀𝑀), 𝑀𝑀𝑀𝑀𝑀𝑀. 𝐴𝐴 𝑢𝑢 (𝑀𝑀)} iii.

𝐴𝐴𝑃𝑃𝑃𝑃𝑀𝑀𝐴𝐴 3 𝑃𝑃𝑀𝑀𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑀𝑀𝑀𝑀𝑃𝑃𝑃𝑃𝑀𝑀 = ∅{𝑀𝑀𝑀𝑀𝑀𝑀. 𝑃𝑃 𝑑𝑑 (𝑀𝑀, 𝑀𝑀), 𝑀𝑀𝑀𝑀𝑃𝑃. 𝑃𝑃 𝑙𝑙 (𝑀𝑀, 𝑀𝑀), 𝑀𝑀𝑀𝑀𝑀𝑀. 𝐴𝐴 𝑢𝑢 (𝑀𝑀, 𝑀𝑀)} In all these arrays, inter-pattern spacing is held constant. 2. The performance of Array 1 is found to be dependent on two out of three functions  power loss and area utilization  dependent on single variable, in-line spacing. Array 2 performance was found to be univariate parameterized in terms of two significant functions -power density and power loss. Array 3 performance was a multi-variable relation involving all the three functions studied. 3. The mathematical models were found to be in close proximity to the simulated values as shown in the validation, but, the model is strictly valid for the feasible domain set for the work. 4. It is found that for maximum power requirement Array 2 is advantageous while in terms of power per rotor unit, Array 3 is beneficial. Array I is not beneficial in either way. 5. It has been observed that when compared to the arrays without any spacing between patterns, arrays with inter-spacing performs better in terms of power per rotor units.

Fig. 1 .

 1 Fig. 1. Typical computational domain showing turbines, velocity inlet, pressure outlet as well as rotorstator interface.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Steps adopted to arrive at optimum turbine array

Table 1 :

 1 Definition of different parameters given in objective criteria. Spacing is shown in Fig.3.

Fig. 5 .

 5 Fig. 5. Variation of velocity deficit with increasing distance from a turbine Figure 8 (a and b) shows the variation of A u and P d with z. A u , as defined earlier, incorporates three variables -spacing, number of patterns and virtual area. As the spacing increases, number of patterns and virtual area values change. With increase in the spacing between the in-line turbines, the length of the rectangular pattern increases and hence this function increases. But after a certain value (4.5 m) the spacing reduces the number of patterns that can be accommodated in the array, this accounts for a jump as observed in Fig. 8(a). Fig. 8(b) shows the variation of power density for array 1 with in-line spacing, which is found to be negligible, hence this function is ignored safely.

Fig. 6 .

 6 Fig. 6.Variation of power loss with in-line spacing. Inset shows the position of the second turbine in the wake of the first for different separating distance.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Comparative study of the effect of in-line (zdirection) and lateral (x-direction) separation distances on the power produced by a second turbine

Fig. 9 .

 9 Fig. 9. Flow velocities of three rotors for Array 2.Points indicate velocities at rotor locations.

Fig. 10 .

 10 Fig. 10. Variation of normalized power density and area utilization ratio with side of the triangular pattern.Figure12shows the pareto-domain for the above two functions generated using genetic algorithm in Matlab®. The solutions are found to get converged after 133 generations. Table2gives details of parameters used to rank pareto-domains based on NFM. The ranking using NFM is shown in Table3. This ranking is indicative of the scores obtained as shown in Table3. It may be added here that these rankings given in Table3were arrived at after studying the sensitivity of different weights and thresholds as stated earlier.

Fig. 11 .

 11 Fig. 11. Different curve fitting functions used to fit numerical data obtained for Array 2 for two functions  power density (a-c) and area utilization ratio (d-f).

Fig. 12 .

 12 Fig. 12. Pareto optimal domain for two functions of power density (f 1 ) and area utilization ratio (f 2 ).

Fig. 13 .

 13 Fig. 13. Variation of power of rotor 2 and rotor 3 with respect to in-line and transverse spacings

Fig. 15 .

 15 Fig. 15. Response surfaces and contour plots of metamodels adopted.

Fig. 16 .

 16 Fig. 16. Pareto-domain for Array 3 obtained using genetic algorithm.

Fig. 17 .

 17 Fig. 17. Comparison of pressure contours for general configuration for Array 1 (left figures) and optimized configuration (right figures).

Fig. 18 .

 18 Fig. 18. Comparison of pressure contours for general configuration for Array 2 (left figures) and optimized configuration (right figures).

Fig. 19 .

 19 Fig. 19. Comparison of pressure contours for general configuration for Array 3 (left figures) and optimized configuration (right figures).

Table 1 .

 1 Area utilization ratio P l -Power loss due to wake [W] x 1 , x 2 -geometrical dimensions of a pattern corresponding in-line spacing and transverse spacing s -spacing between the patterns as defined in

		𝑃𝑃 𝑑𝑑 =	𝑇𝑇𝑃𝑃𝑇𝑇𝑀𝑀𝑙𝑙 𝐴𝐴𝑃𝑃𝑃𝑃𝑀𝑀𝐴𝐴 𝑃𝑃𝑃𝑃𝑃𝑃𝑀𝑀𝑃𝑃 𝐴𝐴𝑃𝑃𝑃𝑃𝑀𝑀𝐴𝐴 𝑀𝑀𝑃𝑃𝑀𝑀𝑀𝑀
	𝐴𝐴 𝑢𝑢 =	𝑃𝑃𝑀𝑀𝑇𝑇𝑇𝑇𝑀𝑀𝑃𝑃𝑃𝑃 𝐴𝐴𝑃𝑃𝑀𝑀𝑀𝑀 * 𝑁𝑁𝑃𝑃. 𝑃𝑃𝑃𝑃 𝑝𝑝𝑀𝑀𝑇𝑇𝑇𝑇𝑀𝑀𝑃𝑃𝑃𝑃𝑠𝑠 𝐴𝐴𝑃𝑃𝑃𝑃𝑀𝑀𝐴𝐴 𝐴𝐴𝑃𝑃𝑀𝑀𝑀𝑀 + 𝑉𝑉𝑀𝑀𝑃𝑃𝑇𝑇𝑢𝑢𝑀𝑀𝑙𝑙 𝐴𝐴𝑃𝑃𝑀𝑀𝑀𝑀
	𝑃𝑃 𝑙𝑙 = 𝑃𝑃 2 -𝑃𝑃 2 ′ P 2 is obtained using the power velocity relation,

𝑙𝑙𝑙𝑙 ≤ 𝑀𝑀

1 , 𝑀𝑀 2 , 𝑠𝑠 ≤ 𝑢𝑢𝑙𝑙 𝑀𝑀 1 , 𝑀𝑀 2 , 𝑠𝑠 ≥ 0, where, P d -Power density [W/m 2 ] A u -

Table 2 :

 2 Parameters used to rank pareto-solutions using Net Flow Method for Array 2

	Array 2	Weight	Q	P	V
	Power density	0.8	5	12	20
	Area utilization ratio	0.20	0.05	0.15	0.20

Table 3 :

 3 List of pareto-optimal solutions for Array 2 and their ranking using NFM method

	Sol	Objectives	Scoring/Ranking
	a	F1	F2	Score	P/NP	Rank
	8.50	13.22	0.44	3.51	P	4
	2.99	35.50	0.15	7.16	P	2
	3.87	28.13	0.20	1.81	P	7
	6.25	17.96	0.32	-3.17	NP	15
	4.18	26.23	0.21	0.54	P	8
	5.94	18.86	0.31	-3.62	NP	17
	7.36	15.29	0.38	0.42	P	10
	6.91	16.29	0.36	-0.47	NP	11
	5.25	21.2	0.27	-2.13	NP	14
	8.02	14.03	0.42	2.87	P	5
	5.68	19.69	0.29	-3.26	NP	16
	5.54	24.32	0.23	-0.97	NP	12
	4.19	26.17	0.22	0.49	P	9
	3.44	31.29	0.18	3.79	P	3
	2.83	37.29	0.14	8.33	P	1
	4.99	22.27	0.26	-1.67	NP	13
	3.87	28.16	0.20	1.83	P	6

Table 4 :

 4 Parameters used to rank pareto solutions using NFM for Array 3

	Array 3	Weight	Q	P	V
	Power density	0.4	3	4	5
	Power loss	0.4	30	40	50
	Area utilization ratio	0.2	0.03	0.05	0.07

Table 5 :

 5 List of pareto-optimal solutions for Array 3 and their ranking using NFM

	Sol		Objectives		Scoring/Ranking
	x	z	F1	F2	F3	Score P/NP Rank
	2.5	4	11.04	247.37	0.06	4.4	P	5
	2.63	4	11.06	249.4	0.06	4.65	P	3
	5.5	8.49	4.85	420.27	0.16	-8.31	NP	18
	2.99	4.42	10.25	287.82	0.08	-1	NP	13
	3.52	5.87	9.67	320.74	0.11	-0.02	NP	12
	5.41	8.47	4.86	420.83	0.16	-8.39	NP	20
	5.02	7.37	7.9	331.02	0.15	5.88	P	1
	2.68	4.06	10.93	255.89	0.07	4.41	P	4
	5.18	7.91	6.87	361.44	0.15	2.3	P	9
	5.5	8.18	6.11	379.81	0.16	-3.73	NP	15
	2.98	4.61	10.08	297.04	0.08	-4.11	NP	16
	2.92	4.23	10.54	274.46	0.08	2.87	P	6
	5.5	8.49	4.85	420.27	0.16	-8.31	NP	19
	3.88	6.76	9.43	323.36	0.12	1.9	P	10
	5.5	8.36	5.43	401.85	0.16	-6.88	NP	17
	3.36	5.2	9.6	317.57	0.1	-2.45	NP	14
	5.43	7.72	7.27	341.51	0.15	5.3	P	2
	4.03	7.23	8.94	333.68	0.13	2.62	P	7
	4.15	6.23	9	322.33	0.12	1.86	P	11
	5.5	7.99	6.69	360.31	0.16	2.62	P	8

Table 6 :

 6 Validation of models adopted from optimized results with rank 1 of respective arrays In order to appreciate the superior performance of the optimized solutions, flow structures in different turbine arrays are compared and shown in Figs. 17-19. Figure17shows the comparative flow pictures for an initial configuration and an optimized array of turbines. It is seen that downstream turbines, in case of general array (left figures), are immersed in the wake of the upstream turbines. This not only reduces the approaching velocities for these subsequent turbines but, because of non-uniform velocities, produces a further reduction in the output power. The optimized turbine array produces wake which nominally influences the subsequent turbines. Similar results are presented for other two configurations in Figs. 18 and 19.

	Array	functions	model actual variation
	2	Power density	37.22	32.41	14.84%
		Area utilization	0.142	0.147	3.4%
		ratio
	3	functions	model actual
		Power density	7.38	7.93	6.93%
		Power loss	363	361.29	0.47%
		Area utilization	0.154	0.157	1.91%
		ratio

Table 7 :

 7 Comparison of the performance of initial and optimized solutions for turbine array geometry

	Array	Parameters	Initial	Optimized
	1	Rotor units	10	4
		Total Power (W)	257.81	360.18
	2	Rotor units	9	6
		Power	1903.70	1986.65
	3	Rotor units	5	3
		Power	1125.99	1026.33