2.1. Problem set-up

The problem under consideration is as shown in figure 1. A cylinder with height $D$ and span $W$ is mounted to a spring–mass–damper system with total mass $m$ , spring constant $k$ and damping constant $c$ . The cylinder is exposed to a fluid flow with free-stream velocity $U$ , and the density, dynamic viscosity and kinematic viscosity of the fluid are $\rho$ , $\mu$ and $\nu$ , respectively.

Figure 1. Sketch of the problem set-up. A cylinder with height $D$ and span $W$ is immersed in a fluid flow with free-stream velocity $U$ . The cylinder has a single degree of freedom, the displacement $y$ in the transverse direction, and is attached to a spring–damper system with mass $m$ , stiffness $k$ and damping coefficient $c$ . Finally, $\rho$ , $\mu$ and $\nu$ are the density, dynamic viscosity and kinematic viscosity of the fluid, respectively.

Non-dimensionalising by $D$ , $U$ and $\rho$ , the equations of motion for the fluid are the Navier–Stokes equations

(2.1a) \begin{align} \frac {\mathrm d \boldsymbol{u}}{\mathrm d t} + \boldsymbol{u} \boldsymbol{\cdot } \boldsymbol{\nabla } \boldsymbol{u} &= - \boldsymbol{\nabla } p + \frac {1}{\textit{Re}}{\nabla} ^2 \boldsymbol{u}, \end{align}

(2.1b) \begin{align} \boldsymbol{\nabla } \boldsymbol{\cdot } \boldsymbol{u} &= 0, \\[9pt] \nonumber \end{align}

where $\boldsymbol{u}$ and $p$ are the non-dimensional fluid velocity and pressure, respectively, while $Re = U D/\nu$ is the Reynolds number. The boundary conditions on the cylinder’s surface are $\boldsymbol{u} = \dot {y}$ , where $y$ is the dimensionless cylinder displacement, and the overdot denotes the non-dimensional time derivative.

The equation of motion for the cylinder displacement, in dimensionless form, is

(2.2) \begin{equation} \ddot {y} + \frac {4\pi \xi }{ U^*} \dot {y} + \frac {4\pi ^2}{{U^*}^2} y = \frac {C_{\!y}}{m^*}, \end{equation}

where $C_{\!y}$ is the transverse force coefficient, $m^*$ is the dimensionless mass, $U^*$ is the reduced velocity and $\xi$ is the damping ratio. These dimensionless parameters are defined as

(2.3) \begin{align} C_{\!y} &= \frac {F_y}{\dfrac {1}{2}\rho U^2 D W}, \\[-12pt] \nonumber \end{align}

(2.4) \begin{align} U^* &= \frac {2\pi U}{ D \sqrt {k/m}}, \\[-12pt] \nonumber \end{align}

(2.5) \begin{align} m^* &= \frac {m}{\dfrac {1}{2} \rho D^2 W}, \\[-12pt] \nonumber \end{align}

(2.6) \begin{align} \xi &= \frac {c}{2\sqrt {km}}, \\[9pt] \nonumber \end{align}

where $F_y$ is the transverse force applied to the body.

We define the reduced velocity with respect to the undamped natural frequency of the system in a vacuum, $U^* = U/(f_n D)$ , where $f_n = 1/(2\pi ) \sqrt {k/m}$ . Note that a range of other frequencies are often used to define the reduced velocity, such as the natural frequency in still water, or the natural frequency including the added mass effect (Sarpkaya Reference Sarpkaya2004). We define the reduced velocity as (2.4), since this form of the reduced velocity is determined completely by the system parameters alone.

For the dimensionless mass, previous studies typically consider the mass ratio, which is the ratio between the total oscillating mass and the mass of displaced fluid

(2.7) \begin{equation} m_r = \frac {m}{\rho W A_c}, \end{equation}

where $A_c$ is the cross-sectional area of the cylinder. The dimensionless mass and mass ratio are related as $m^*/m_r = 2A_c/D^2$ . Note that for a fixed total mass, $m^*$ is independent of the cylinder’s shape, while $m_r$ depends on the cylinder’s cross-sectional area. In the present study, optimisations are performed with either a fixed $m^*$ , corresponding to a constant total mass, or for a fixed $m_r$ , corresponding to a fixed density of the solid body.

Numerically solving the Navier–Stokes (2.1), coupled with the equation of motion for the cylinder (2.2), yields a predicted vibration response $y(t)$ and both velocity $\boldsymbol{u}$ and pressure $p$ fields, for a given cylinder shape and combination of the dimensionless parameters $m^*$ , $U^*$ , $\xi$ and $Re$ . The amplitude of vibration is characterised by either the maximum amplitude ( $y_{\textit{max}}$ ), or root mean square (r.m.s.) amplitude ( $y_{\textit{rms}}$ ).

For power generation applications, the attached generator must extract mechanical energy from the system. We model the energy extraction as a linear damper, with damping coefficient $c$ . Assuming perfect efficiency of the generator, and no other sources of damping, the mean power coefficient is

(2.8) \begin{equation} \bar {C}_P = \frac {\bar {P}}{\dfrac {1}{2} \rho U^3 D W} = \frac {c}{\dfrac {1}{2} \rho U D W} \,\overline {\dot {y}^2} = \frac {4\pi m^* \xi }{U^*} \,\overline {\dot {y}^2}, \end{equation}

where $P$ is the dimensional power output, and overbars represent time-averaged quantities. In practice, there may be significant energy losses due to other sources of damping (e.g. friction in the bearing system), as well as losses in the transducer used to convert mechanical energy into electrical power. As a result, the power coefficients and efficiencies reported in this paper refer only to the performance of the cylinder (i.e. the transfer of energy from the fluid to the cylinder), not the performance of a realistic energy conversion system. Additionally, real power conversion devices may not impart a linear damping on the system.

Note that $\bar {C}_P$ is non-dimensionalised by the projected area of the cylinder perpendicular to the flow direction. However, it is typical in power generation applications to non-dimensionalise the power extraction by the total swept area of the device. This represents an efficiency by which power is extracted from the total power available in the fluid. The total swept area for a single cylinder is $\textit{WD}(1+2y_{\textit{max}})$ , and therefore, the efficiency is

(2.9) \begin{equation} \eta = \frac {\bar {{P}}}{\dfrac {1}{2}\rho U^3 \textit{WD}(1+2y_{\textit{max}})} = \frac {\bar {C}_{P}}{1+2y_{\textit{max}}}. \end{equation}

As discussed by Lo et al. (Reference Lo, Thompson, Hourigan and Zhao2024), increasing the maximum amplitude of the cylinder will generally increase the power coefficient. However, this may not lead to an increase in the efficiency, since the swept area also increases. The present study focuses on optimising the power coefficient ( $\bar {C}_p$ ) of a single cylinder, representing the maximum power extraction from a single device. In practical applications, one may wish to extract the most energy from a given environmental fluid flow, in which case the efficiency ( $\eta$ ) should be optimised instead. In such cases, one should consider the efficiency of an array containing multiple FIV devices, in either a tandem or parallel configuration, which is beyond the scope of the present study. Also note that the cost of a single VIVACE device is not likely to depend substantially on the vibration amplitude. Therefore, for a sufficiently abundant source of environmental energy, it may be more cost effective to maximise power coefficient, rather than the efficiency.

2.2. Parametrising the geometry

The Bayesian optimisation technique used in this work requires the geometry to be parametrised. In this section we introduce the various parameters that control the shape of the ellipse. These parameters are chosen to allow various elementary shapes, such as circles, ellipses, squares, D-section cylinders, diamonds, triangles to be represented.

We first consider the following generalisation of the ellipse (superellipse):

(2.10) \begin{equation} |x/a|^{2/m} + |y/b|^{2/n} = 1, \end{equation}

which can be defined parametrically as

(2.11) \begin{equation} x = a\,\mathrm{sign}(\sin (\theta ))\,|\sin (\theta )|^{m}\,,\,y = b\,\mathrm{sign}(\cos (\theta ))\,|\cos (\theta ) |^{n}\,,\, 0\leqslant \theta \leqslant 2\pi . \end{equation}

For $n = m = 1$ , (2.10) describes an ellipse with width $2a$ and height $2b$ . For the present study, the cylinder has unit height in dimensionless coordinates ( $b = 0.5$ ), while $a = 1/(2\textit{AR})$ , where $ \textit{AR}$ is the cylinder aspect ratio. Varying the parameters $n$ and $m$ varies the shape of the cylinder, while maintaining a constant width and height. As shown in figure 2, the superellipse can represent a diamond ( $n = m =2$ ), a rectangle ( $n = m = 0$ ), or a shape with parabolic edges ( $m = 2$ and $n = 1$ or $m = 2$ and $n = 1$ ).

Figure 2. Effect of the parameters $m$ and $n$ on the shape of the generalised superellipse. (a) Effect of varying both $m$ and $n$ , with $m = n$ ; (b) effect of varying $m$ , with $n =1$ ; and (c) effect of varying $n$ , with $m = 1$ .

Previous studies have identified that the forebody and afterbody play different roles in generating FIV motion (Zhao, Hourigan & Thompson Reference Zhao, Hourigan and Thompson2018; Chen & Li Reference Chen and Li2023; Gupta et al. Reference Gupta, Lo, Zhao, Thompson and Hourigan2025). Therefore, we allow the shape parameters to differ for the front and rear surfaces. First, we restrict $b = 0.5$ for both the forebody and afterbody. Then, we define $a_{\!f}$ , $m_{\!f}$ and $n_{\!f}$ as the parameters for the forebody ( $x\lt 0$ ), and $a_a$ , $m_a$ and $n_a$ as the parameters for the afterbody. We note that $a_a$ and $a_{\!f}$ are constrained by the cylinder aspect ratio, $a_{\!f} + a_a = 1/\textit{AR}$ . We also define a parameter $\alpha$ , which controls the relative width of the forebody and afterbody

(2.12) \begin{equation} {(1 - \alpha ) a_{\!f} = (1+\alpha ) a_a\,,\, -1\leqslant \alpha \leqslant 1.} \end{equation}

As shown in figure 3, using either $\alpha = -1$ or $\alpha = 1$ results in a shape with no forebody or no afterbody, respectively, including D-section ( $n = m = 1$ ) and triangular ( $n = m = 2$ ) cylinders.

Figure 3. Effect of the parameter $\alpha$ on the cylinder shape. For (a,c) $\alpha = -1$ , the cylinder has only an afterbody, while for $\alpha = 1$ , the cylinder has only a forebody. This allows the cylinder to be either a D-section (a,b) or triangular cylinder (c,d).

Finally, we introduce a parameter to control the camber of the cylinder. We define $x_{1,f}(y)$ and $x_{1,a}(y)$ as the forebody and afterbody of the uncambered cylinder, which are computed using the six parameters $ \textit{AR}$ , $\alpha$ , $n_{\!f}$ , $m_{\!f}$ , $n_a$ , $m_a$ as discussed above. We then define the following cambered curve:

(2.13) \begin{equation} x_c = \beta (1 - 2 y_c)^2, -0.5\leqslant y_c \leqslant {0.5}, \end{equation}

where $\beta$ represents the amount of camber. We then compute a normalised arc-length parametrisation of this curve, $x_c(s)$ , $y_c(s)$ , with $-0.5\leqslant s \leqslant 0.5$ . Finally, the cylinder shape is given by

(2.14) \begin{equation} x = x_c(s) + n_x(s) x_1(s), y = y_c(s) + n_y(s) x_1(s), \end{equation}

where $n_x$ and $n_y$ are components of a unit vector. For thin cylinders ( $ \textit{AR}\gt 3$ ), this unit vector is chosen to be the unit normal to the camber line, so that the uncambered cylinder represents a thickness distribution about the camber line. To avoid self-intersections, we use $n_x = 1$ and $n_y = 0$ for thick cylinders ( $ \textit{AR}\leqslant 1$ ). For $1\leqslant \textit{AR} \leqslant 3$ , the angle of the unit vector is linearly interpolated between these values.

Figure 4 demonstrates the application of camber to a cylinder. Panel (a) shows the un-cambered cylinder, with parameters $ \textit{AR} = 5$ , $\alpha = 0$ , $r_{m,f} = 0.5$ , $r_{n,f} = 0.5$ , $r_{m,a} = 2$ and $r_{n,a} = 1$ . Panel (b) shows both the camber curve ( $x_c$ , $y_c$ ), and the cambered cylinder ( $x$ , $y$ ).

Figure 4. Addition of camber to the cylinder. (a) shows the uncambered cylinder $(x_1,y_1)$ , while (b) shows the camber curve $(x_c,y_c)$ and the cambered cylinder $(x,y)$ .

2.3. OpenFOAM simulations

Numerical simulations are performed using the open-source software OpenFOAM 10.0 (Weller et al. Reference Weller, Tabor, Jasak and Fureby1998), with the PIMPLE algorithm (a combination of the Pressure-Implicit with Splitting of Operators (PISO) and Semi-Implicit Method for Pressure Linked Equations (SIMPLE) algorithms) used for pressure–velocity coupling. The equations of motion are discretised using the finite-volume scheme. The linear interpolation scheme (equivalent to second-order centred differencing) was used for all spatial derivatives, aside from the advection term, for which a second-order linear upwind scheme was used. The first-order backwards Euler scheme was used for temporal derivatives. An adaptive time-step scheme was used, with a maximum Courant number of $ \textit{Co}_{\textit{max}}=0.9$ .

A dynamic mesh, using the sixDoFRigidBodyMotion utility in openFOAM, was used to account for the motion of the cylinder. First, the viscous and pressure forces applied to the cylinder are computed from the finite-volume discretisations of the velocity and pressure fields. Then, the position of the cylinder is updated using (2.2). Grid points within 30 diameters of the cylinder move rigidly with the cylinder, while the position of the remaining grid points are updated using a Laplacian mesh motion solver.

We consider a two-dimensional square computational domain, with boundaries located at a distance $L_B$ from the cylinder’s rest position. Specifically, there is an inlet located at $x = -{L_B}$ , an outlet at $x = {L_B}$ and upper and lower free-slip walls at $y=\pm {L_B}$ . Boundary conditions are $\boldsymbol{u} = (0,\dot {y},0)$ at the cylinder’s surface, $\boldsymbol{u} = (1,0,0)$ at the inlet, $p = 0$ at the outlet and $\partial u_x/\partial y = 0$ on the free-slip walls.

To achieve reasonable mesh quality for an arbitrary cylinder shape, unstructured meshes are generated programmatically using the open-source software GMSH (Geuzaine, Remacle & Dular Reference Geuzaine, Remacle and Dular2009). An example mesh is shown in figure 5. Increased resolution is used near the body and in the wake region behind the body. The algorithm used by GMSH attempts to produce elements with a pre-defined target cell size at each location in the computational domain. These target cell sizes are illustrated in figure 5: $l_1$ is the element size near the solid body; $l_2$ is the target element size at a distance $0.5$ from the body; $l_3$ is the target element size in the wake (at a distance $6$ from a line segment passing through the points $(0,0)$ and $(0,5)$ ); and $l_4$ is the target element size in the far-field. Linear interpolation is used to determine the target element size at an arbitrary location in the domain.

Figure 5. Sketch illustrating the unstructured mesh used in this study (not to scale). The parameters $l_1$ , $l_2$ , $l_3$ and $l_4$ represent approximate cell sizes (side lengths of the triangular cells) at various distances from the cylinder.

An inflation layer containing high aspect ratio quadrilateral elements is used near the boundary of the body, while triangular elements are used for the remaining region. The thickness of the first layer is $0.1 l_1$ , and the growth rate is $1.1$ . Curvature-based mesh adaptation is used with $n_c$ boundary elements per radian of curvature. A boundary-layer fan is employed at sharp corners (with angle greater than $0.2$ radians). The number of boundary-layer fan elements per radian of angle is set to $n_{\!f}$ .

To ensure accurate determination of the objective function, it is important that the flow reaches a statistically steady state. Figure 6 plots time histories of cylinder displacement for several different cylinder shapes. While the time taken to reach a steady state differs for each cylinder, in all cases a steady state is approximately reached after $t=200$ , and this steady state is close to periodic. Temporally averaged statistics (mean power coefficient, r.m.s. displacements) are obtained over a time interval spanning the first local maximum of cylinder displacement after $t = 200$ , to the last local maximum before $t = 500$ . For the maximum cylinder displacement, the average value of each local maximum between $t = 200$ and $t = 500$ is reported.

Figure 6. Time histories of cylinder displacement for a variety of different cylinders: (a) a circular cylinder at $U^* = 6$ and $\xi =0$ , (b) the OpU10 cylinder (see figure 9) at $U^* = 11$ and $\xi = 0$ , (c) the semi-elliptical cylinder SE6 (see table 7) at $U^* = 6$ and $\xi = 0.0293$ and (d) the triangular cylinder T7 (see table 7) with $U^* = 7$ and $\xi = 0.193$ . For all cases, $Re = 100$ and $m^* = 4.7124$ .

To reduce the startup time, the cylinder is given an initial velocity of $\dot {y} = 2$ at $t = 0$ . In some cases, this may trigger ‘hard’-oscillating behaviour (Novak & Tanaka Reference Novak and Tanaka1974), in which the cylinder only oscillates given a sufficiently large initial displacement. For cases where the oncoming flow is relatively steady, it is reasonable to assume that a power generating device could be manually started in this manner. For devices which experience substantial gust loads, the transient behaviour must be considered (Alhussein et al. Reference Alhussein, Dalaq and Daqaq2024), which could result in a different optimal shape to those obtained in the present study.

The present study considers two-dimensional laminar flow at low Reynolds numbers ( $Re \approx 100$ ), to keep the computation cost per simulation low. The assumption of two-dimensionality restricts the results to wide cylinders ( $W/D \gg 1$ ). For future research, it would be interesting to perform the optimisation at larger Reynolds numbers (i.e. $Re \approx 10^{4}$ ), using an appropriate turbulence model, and validate the results against experimental measurements. Additionally, other effects relevant to practical applications, such as the effects of gust loads or the angle of oncoming flow, should be considered in the future.

2.4. Validation of numerical scheme

The mesh resolution and domain size are validated using an optimised cylinder obtained in the preliminary stages of this project. Table 2 provides the target element sizes for three meshes with increasing resolution, while table 3 reports the amplitude obtained using each cylinder at $U^* = 6$ and $U^*=10$ , as well as the r.m.s. velocity of the cylinder $\overline {\dot {y}^2}$ and Strouhal number based on cylinder displacement $St_y$ . The differences between meshes 2 and 3 are below $1\,\%$ , and therefore mesh 2 is satisfactory for the present study.

Table 2. Parameters for the meshes used in the grid resolution study. The parameters $l_1$ , $l_2$ $l_3$ and $l_4$ are the target cells sizes shown in figure 5, while $n_c$ and $n_{\!f}$ are the minimum number of cells per radian on the curved boundary and at sharp corners, respectively. The approximate number of cells, $N_c$ , is also shown.

Table 3. Mesh resolution study comparing the oscillation amplitude ( $y_{\textit{max}}$ ), the r.m.s. velocity of the cylinder ( $\overline {\dot {y}^2}$ ) and the Strouhal number based on the cylinder displacement $St_{y}$ , obtained on the three different meshes listed in table 2, at various reduced velocities ( $U^*$ ). Simulations are performed for a preliminary optimised cylinder at $Re = 100$ and $m^* = 4.7124$ . Parentheses indicate the relative error compared with mesh 3.

Simulations were performed using the Setonix supercomputer located at the Pawsey Supercomputing Research Centre in Perth, Western Australia. Simulations using mesh 2 require approximately 33 hours of wall clock time, using 8 CPU cores in parallel. Due to the small size of this problem, further increasing the number of CPU cores does not substantially improve the time. Simulations using mesh 1 typically require 6 hours of wall time on a single CPU core. To speed up the optimisation process, an initial optimisation is performed using mesh 1 to efficiently explore the entire design space. Then, a second optimisation is performed over a reduced design space, using mesh 2, to obtain a more accurate solution. Further details are provided in § 4.

Table 4 reports the variation of oscillation amplitude ( $y_{\textit{max}}$ ), r.m.s. velocity of the cylinder $\overline {\dot {y}^2}$ and Strouhal number based on cylinder displacement $St_y$ obtained using various domain sizes ( $L_B$ ). The difference between ${L_B} = 200$ and ${L_B} = 400$ are below $2 \,\%$ , and therefore ${L_B}=200$ is considered satisfactory for this project. Note that ${L_B} =50$ was found to be satisfactory for circular cylinders, which oscillate at a much lower amplitude. The large domain size needed to achieve a domain-independent solution suggests that the effective blockage ratio is related to the amplitude of oscillations, rather than the cylinder diameter.

Table 4. Domain size study comparing the oscillation amplitude ( $y_{\textit{max}}$ ), the r.m.s. velocity of the cylinder ( $\overline {\dot {y}^2}$ ) and the Strouhal number based on the cylinder displacement $St_{y}$ obtained for different domain sizes $L_B$ using the mesh 2 resolution, reduced velocities ( $U^* = 6$ and $U^* = 10$ ). Simulations are performed for a preliminary optimised cylinder at $Re = 100$ and $m^* = 4.7124$ . Parentheses indicate the relative error compared with ${L_B} = 400$ .

Figure 7 presents a validation study, comparing the present numerical scheme with previous studies of Leontini, Thompson & Hourigan (Reference Leontini, Thompson and Hourigan2006) and Rajamuni, Thompson & Hourigan (Reference Rajamuni, Thompson and Hourigan2020), for the VIV of a circular cylinder with $m_r = 10$ , $\xi = 0.01$ and $Re = 200$ . Reasonable agreement between the present simulations and the prior studies is observed, validating the numerical scheme.

Figure 7. Comparison between the present numerical simulations and previous studies of Leontini et al. (Reference Leontini, Thompson and Hourigan2006) and Rajamuni et al. (Reference Rajamuni, Thompson and Hourigan2020), showing the oscillation amplitude $y_{\textit{max}}$ against $U^*$ , for a circular cylinder with $m_r = 10$ , $\xi = 0.01$ and $Re = 200$ .

3.1. Overview of Bayesian optimisation with Gaussian process regression

We now give a brief overview of the BO-GPR method, based on the presentations given in Morita et al. (Reference Morita, Rezaeiravesh, Tabatabaei, Vinuesa, Fukagata and Schlatter2022) and Diessner et al. (Reference Diessner, O’Connor, Wynn, Laizet, Guan, Wilson and Whalley2022). We let $\boldsymbol{q}$ be a vector of design parameters, and $r(\boldsymbol{q})$ be the objective function to be minimised. The first step in Bayesian optimisation (BO) is, given a set of $n$ measurement data $D_n = {(\boldsymbol{q}_i,r_i)}_{i=1}^n$ , to estimate the posterior probability distribution for the measurement at the next sample point $\boldsymbol{q}_{n+1}$ . The value of $\boldsymbol{q}_{n+1}$ that maximises some acquisition function is then selected and the response value $r_{n+1} = r(\boldsymbol{q}_{n+1})$ is obtained. This process is then repeated until an optimum solution is found.

In BO-GPR, the objective function is approximated as a Gaussian process (GP). A GP is a stochastic process in which every finite collection of random points is represented by a multivariate normal distribution. The GP is defined by a mean function $m(\boldsymbol{q})$ , and a covariance function (kernel) $c(\boldsymbol{q}, \boldsymbol{q}')$ . For the set of sample points $\boldsymbol{q}_{1,2,\ldots ,n}$ (where $\boldsymbol{q}_i$ represents the $i$ th measurement), we have a vector of measurement points $\boldsymbol{r}$ , a vector of mean values $\boldsymbol{m}$ and a covariance tensor $\unicode{x1D63E}$ , with components defined as $r_i = f(\boldsymbol{q}_i)$ , $m_i = m(\boldsymbol{q}_i)$ and $C_{i,j} = c(\boldsymbol{q}_i,\boldsymbol{q}_{\kern-1.5pt j})$ . The prior probability distribution for the measurement points (i.e. the assumed probability distribution prior to obtaining the measurement data) is a multivariate normal distribution

(3.1) \begin{equation} \boldsymbol{r} \sim N(\boldsymbol{m}, \unicode{x1D63E}\kern1pt). \end{equation}

The posterior predictive distribution for the next sample point $r_{n+1} = f(\boldsymbol{q}_{n+1})$ is also a normal distribution

(3.2) \begin{equation} r_{n+1} | D_n = N(\mu _n(\boldsymbol{q}_{n+1}),\sigma _n(\boldsymbol{q}_{n+1})), \end{equation}

with the posterior mean and variance given by

(3.3) \begin{align} \mu _n(\boldsymbol{q}_{n+1}) &= m(\boldsymbol{q}_{n+1}) + \boldsymbol{c}(\boldsymbol{q}_{n+1})^T(\unicode{x1D63E} + \sigma ^2 \unicode{x1D644}\kern1pt)^{-1} (\boldsymbol{r} - \boldsymbol{m}), \\[-12pt] \nonumber \end{align}

(3.4) \begin{align} \sigma _n^2(\boldsymbol{q}_{n+1}) &= c(\boldsymbol{q}_{n+1},\boldsymbol{q}_{n+1}) - \boldsymbol{c}(\boldsymbol{q}_{n+1})^T (\unicode{x1D63E} + \sigma ^2 \unicode{x1D644}\kern1pt)^{-1} \boldsymbol{c}(\boldsymbol{q}_{n+1}), \\[9pt] \nonumber \end{align}

where $\boldsymbol{c}(\boldsymbol{q}_{n+1})$ is a vector of covariances between the candidate point and all design points ( $c_i = c(\boldsymbol{q}_{n+1},\boldsymbol{q}_i)$ ). The result is a probability distribution for the objective function, sampled at the candidate point $\boldsymbol{q}_{n+1}$ .

In this study, we use the zero mean function $m(\boldsymbol{q}) = 0$ for the prior distribution. The sample data $r_i$ are also normalised to have zero mean, and unit standard deviation. For the covariance function $c$ , the Matérn 5/2 kernel is used, following Diessner et al. (Reference Diessner, O’Connor, Wynn, Laizet, Guan, Wilson and Whalley2022),

(3.5) \begin{equation} c(\boldsymbol{q},\boldsymbol{q}') = \sigma _{\!f}^2 \left ( 1 + \frac {\sqrt {5}}{\rho } + \frac {5 d^2}{4\rho ^2} \right ) \exp \left ( -\frac {\sqrt {5}d}{\rho }\right )\! . \end{equation}

Here, $d = \sqrt {(\boldsymbol{q} - \boldsymbol{q}')\boldsymbol{\cdot }(\boldsymbol{q} - \boldsymbol{q}')}$ is the distance between the two sample points, $\sigma _{\!f}^2$ is a scaling factor and $\rho$ is the length scale. Both $\sigma _{\!f}^2$ and $\rho$ are hyperparameters that are estimated from the sample data using maximum likelihood estimation. Note that (3.5) has a single length scale for all parameters. To account for this, we normalise each of our design parameters to lie in the range $0\leqslant q_i \leqslant 1$ .

Given the posterior probability function for the surrogate objective function, an acquisition function is then applied to select the next sample point(s). In this study, we use the expected improvement (EI) acquisition function, which is one of the most commonly used acquisition functions (Morita et al. Reference Morita, Rezaeiravesh, Tabatabaei, Vinuesa, Fukagata and Schlatter2022). Given the best point found so far, $r^* = f(\boldsymbol{q}^*)$ , the EI is defined as

(3.6) \begin{equation} \alpha _{\textit{EI}} = E(\max (f(\boldsymbol{q}^*) - f(\boldsymbol{q}))), \end{equation}

where $E()$ represents the expected value. The value of $\boldsymbol{q}$ that gives the maximum value of $\alpha _{\textit{EI}}$ – i.e. the design point that gives the greatest EI – is selected as the next design point.

In this study, batch optimisation with local penalisation (Gonzalez et al. Reference Gonzalez, Dai, Hennig and Lawrence2016) is used. The acquisition function is initially defined as in (3.6). The first batch member, $\boldsymbol{q}_{n+1}$ , is selected based on the maximum of this acquisition. Then, the acquisition function is penalised in the neighbourhood of this point, and the next batch element selected based on the penalised acquisition function. This is repeated until all points in the batch are selected. Numerical simulations for each batch are then performed simultaneously.

The BO-GPR method was implemented using the Python libraries GPy (http://github.com/SheffieldML/GPy) and GPyOpt (http://github.com/SheffieldML/GPyOpt). The BO-GPR method was coupled non-intrusively to the OpenFOAM simulations using Bash scripts. The overall process is as follows: first, the next batch of design points are obtained using the BO-GPR method implemented in GPyOPT. Then, meshes are procedurally generated using the Python GMSH package (Geuzaine et al. Reference Geuzaine, Remacle and Dular2009). Next, simulations are performed using OpenFOAM. Finally, the results of the OpenFOAM simulations are post-processed in Python to compute the objective function. These results are then supplied to GPyOpt to select the next batch of design points.

4.1. Optimising with a fixed mass ratio

The optimisations presented in the previous section are performed for a fixed dimensionless mass $m^*$ . This assumes that the total oscillating mass is independent of the cross-sectional area of the cylinder. In this section, optimisations are performed with a fixed mass ratio $m_r = 15$ , which assumes that the total oscillating mass is proportional to the cylinder’s cross-sectional area. Three reduced velocities are considered ( $U^* = 6$ , $U^* = 8$ and $U^* = 10$ ), while the Reynolds number is held constant ( $Re = 100$ ). Cylinders optimised for fixed $m_r$ and $U^* = \mathrm{X}$ are denoted OpUX_MR. The shape parameters for these cylinders are included in table 5.

Figure 13 presents a comparison between the OpU8 and OpU8_MR cylinders. For both cylinders, the aspect ratio is approximately $ \textit{AR} = 20$ , which is the maximum $ \textit{AR}$ considered in this study. While both cylinders have the same nominal $ \textit{AR}$ , the OpU8_MR cylinder has a lower area than the OpU8 cylinder. As defined in this study, $ \textit{AR}$ is the ratio between cylinder height and maximum cylinder width. The OpU8_MR cylinder achieves a lower area by having a sharp kink at the cylinder midpoint, so that the cylinder thickness away from the midpoint is reduced. The OpU8 cylinder has a smooth surface, with a greater overall area. For a fixed mass ( $m^*$ ), the OpU8 cylinder outperforms the OpU8_MR, indicating that the smooth surface leads to improved hydrodynamic performance. However, when the mass ratio ( $m_r$ ) is fixed, the thinner OpU8_MR cylinder has a lower mass, and achieves a greater oscillation amplitude, indicating that the reduction in mass due to a thinner profile outweighs the hydrodynamic penalty due to the sharp kink. Overall, however, the difference in performance between the OpU8 and OpU8_MR cylinders is relatively small.

4.2. Effect of Reynolds number

In this section, optimisations are performed for two different Reynolds numbers, $Re = 60$ and $Re = 80$ , to determine the influence of $Re$ on the optimal cylinder shape within the two-dimensional laminar flow regime. The optimised cylinders are labelled OpR60 and OpR80, respectively. These optimisations are performed for $\xi = 0$ , $m^* = 4.7124$ and $U^* = 10$ . Figure 14 plots the optimised cylinders alongside the OpU10 cylinder, and the shape parameters are presented in table 5. Increasing Reynolds number has only a small effect on the optimal cylinder shape, with a slight decrease in the amount of camber as $Re$ increases from $x_{\textit{cam}} = 0.15$ at $Re = 60$ , to $x_{\textit{cam}} = 0.11$ at $Re = 100$ .

Figure 14. Optimised cylinders obtained for three different Reynolds numbers, at $U^* = 10$ , $m^* = 4.7124$ and $\xi = 0$ .

5.1. Vorticity contours

Figure 15 presents visualisations of the vorticity field in the wake of the OpU9 cylinder, at four different values of $U^*$ . Visualisations are provided at various phase angles $\phi = \omega t$ , where $\omega$ is the oscillation frequency and $\phi = 0$ corresponds to the maximum displacement of the cylinder. Similar vorticity visualisations are obtained for all other optimised cylinders, and are not shown. Animations are provided in the online supplementary material (movies 1–4).

Figure 15. Visualisations of spanwise vorticity in the wake of the OpU9 cylinder, at (a) $U^* = 1$ , (b) $U^* = 3$ , (c) $U^* = 5$ and (d) $U^* = 10$ . Here, $\phi$ is the phase angle, with $\phi = 0$ corresponding to the maximum displacement $y$ of the cylinder. Red indicates positive (counter-clockwise) vorticity, while blue indicates negative (clockwise) vorticity.

For $U^*\leqslant 2$ (figure 15 a), the oscillation amplitude is negligible. In this regime, vortices of alternating orientation are periodically shed from the cylinder, similar to the von Kármán vortex street behind a flat plate (Najjar & Vanka Reference Najjar and Vanka1995).

At $U^* = 3$ (figure 15 b), substantial FIVs begin to occur. As shown in figure 12(b), the frequency of vibration matches the natural frequency of the system. Vortex shedding is also synchronised to the cylinder’s vibrations. Referring to figure 15(b), a single vortex is shed from the cylinder every half-period of the motion, after the cylinder begins to travel away from the point of maximum displacement (between $\phi = 0$ and $\phi = \pi /2$ ). However, there is an interesting interaction with a period of three full vortex shedding cycles. First, a vortex $p1$ is shed from the cylinder, at $\phi = \pi /2$ . A second vortex $p2$ is then shed at the next cycle of motion ( $\phi = 5\pi /2$ ). These two vortices merge to form a single vortex $p12$ ( $\phi = 4\pi$ ). Finally, a third vortex is shed at the next cycle of motion ( $\phi = 9\pi /2$ ). The two vortices $p12$ and $p3$ are then advected downstream of the cylinder.

For $4 \leqslant U^*\leqslant 8$ (figure 15 c) the wake structure is reminiscent of the hyperbranch observed by Lo et al. (Reference Lo, Hourigan, Thompson and Zhao2023) for the thin elliptical cylinder. Specifically, a single vortex is shed from the cylinder every half-cycle of the motion (between $\phi = 0$ and $\phi = \pi /4$ ), after the cylinder changes direction at the point of maximum amplitude. At $\phi = \pi /2$ , where the cylinder’s vertical velocity is large, flow appears to be attached to the cylinder, with behaviour similar to an airfoil. Figure 16 shows velocity streamlines near the cylinder at $\phi = \pi /2$ , computed in a frame of reference attached to the cylinder. The streamlines are qualitatively similar to flow over an airfoil. Specifically, flow attaches to the cylinder at the leading edge, and separates at the trailing edge. However, there is a small recirculation region, indicating some flow separation from the rear surface. The similarity to a lifting airfoil suggests that there will be a relatively low drag force and a significant lift force with a component in line with the cylinder’s velocity, which could enhance FIV. This proposal will be investigated in § 5.2. Moreover, the shedding of a single vortex at the peak cylinder displacement resembles the shedding of a starting/stopping vortex from an airfoil. Attached airfoil-like flow has been reported in previous FIV studies of thin elliptical cylinders (Lo et al. Reference Lo, Hourigan, Thompson and Zhao2023) and curved blades (Tucker Harvey et al. Reference Tucker Harvey, Khovanov and Denissenko2019; Reference Tucker Harvey, Khovanov, Murai and Denissenko2020). Lo et al. (Reference Lo, Hourigan, Thompson and Zhao2023) also report a similar pattern of vortex shedding to that observed in the present study.

Figure 16. Velocity streamlines computed in a frame of reference moving with the cylinder, for (a) $U^* = 5$ and (b) $U^* = 9$ , at the point of maximum vertical velocity ( $\phi = \pi /2$ ). A visualisation of the vorticity field is also provided, with red indicating positive (counter-clockwise) vorticity and blue indicating negative (clockwise) vorticity.

For $U^*\geqslant 9$ , (figure 15 d) the shear-layer attached to the cylinder becomes unstable, resulting in the alternate shedding of many small vortices throughout the motion. These small-scale vortex structures (SVS) were also observed by Lo et al. (Reference Lo, Hourigan, Thompson and Zhao2023) for the elliptical cylinder. In addition to the SVS, a strong primary vortex is still shed from the cylinder each half-cycle of motion, as the cylinder begins to travel away from the point of maximum displacement.

Figure 16 plots streamlines near the cylinder for $U^* = 9$ , at $\phi = \pi /2$ . Once again, the flow is qualitatively similar to the flow over an airfoil. Flow attaches at the leading edge, and there is a separation point at the trailing edge. However, flow separation at the rear surface is more severe at $U^* = 9$ than at $U^* = 5$ , resulting in a larger recirculation bubble.

5.2. Quasi-steady analysis

Galloping is a movement-induced instability driven by changes in the transverse force due to variations in the inclination angle between the cylinder and the relative flow (Parkinson Reference Parkinson1989). Galloping is often analysed using a quasi-steady analysis (Parkinson Reference Parkinson1989; Barrero-Gil, Alonso & Sanz-Andres Reference Barrero-Gil, Alonso and Sanz-Andres2010; Javed, Abdelkefi & Akhtar Reference Javed, Abdelkefi and Akhtar2016; Tucker Harvey et al. Reference Tucker Harvey, Khovanov, Murai and Denissenko2020), which assumes the fluid forces are a function only of the relative angle of attack between the cylinder and the fluid. In this section, we present a quasi-steady analysis of the optimised cylinder, showing that a hard-galloping behaviour is possible under the quasi-steady model. This is motivated by two observations. First, the variation of amplitude with $U^*$ (figure 12) resembles galloping behaviour. Second, the proposed airfoil-like behaviour suggests a strong influence of relative angle of attack on the force coefficients, which may be predicted using quasi-steady analysis.

Note that the quasi-steady assumption is only valid at high $U^*$ , when the oscillation frequency is much lower than the vortex shedding frequency. As shown in § 5.1, the frequency of vortex shedding for the optimised cylinder is equal to the oscillation frequency over the range of $U^*$ considered in this study. Therefore, while the quasi-steady model indicates how the relative angle of attack affects the hydrodynamic forces, it does not capture how unsteady effects such as vortex shedding impact the motion. The influence of vortex shedding on the motion is discussed later, in § 5.3.

As shown in figure 17, the relative velocity and inclination angle between the cylinder and the oncoming flow are

(5.1) \begin{align} U_{\textit{rel}} = U\sqrt {1+\dot {y}^2} &= U\sqrt {1+\tan ^2\alpha }, \end{align}

(5.2) \begin{align} \tan \alpha &= \dot {y}. \end{align}

Figure 17. Sketch of the relative velocity experienced by the cylinder. Here, $\dot {y}$ is the velocity of the cylinder, while $U_{rel}$ and $\alpha$ are the magnitude and incidence angle of the relative velocity. Finally, $C_{\!L}$ and $C_{\!D}$ are the lift and drag coefficients, while $C_{\!y}$ and $C_x$ are the transverse and streamwise force coefficients.

Under the quasi-steady assumption, the lift and drag forces depend only on the inclination angle $\alpha$ . The quasi-steady transverse force can be computed from the lift and drag coefficients as

(5.3) \begin{equation} C_{y,q}(\alpha ) = (1+\tan ^2 \alpha ) (C_{L,q} (\alpha ) \cos \alpha - C_{D,q}(\alpha ) \sin \alpha ), \end{equation}

where $C_{L,q}$ and $C_{D,q}$ are the quasi-steady lift and drag coefficients. The instantaneous quasi-steady power coefficient can also be computed, as

(5.4) \begin{equation} C_{p,q}(\alpha ) = C_{y,q} \dot {y} = \tan \alpha C_{y,q}(\alpha ). \end{equation}

The quasi-steady power coefficient represents the transfer of energy from the fluid to the cylinder, due to the transverse force. When this term is positive, the quasi-steady forces enhance the vibration of the cylinder. Note that $C_{\!L}$ and $C_{\!D}$ are non-dimensionalised by the relative velocity $U_{\textit{rel}}$ , while $C_{\!y}$ and $C_{p,q}$ are non-dimensionalised by the free-stream velocity $U$ . Moreover, the effective Reynolds number for computing $C_{\!L}$ and $C_{\!D}$ is

(5.5) \begin{equation} {\textit{Re}}_{\textit{rel}} = \textit{Re}\sqrt {1+\tan ^2 \alpha }, \end{equation}

where $Re$ is the Reynolds number based on free-stream velocity.

The quasi-steady force coefficients are estimated from numerical simulations performed for a stationary cylinder, with fixed inclination angle and Reynolds number determined by (5.5). The inclination angle is varied between $\alpha = 0$ and $\alpha = ( {7}/{16})\pi$ , resulting in Reynolds numbers between ${\textit{Re}}_{rel} = 100$ and ${\textit{Re}}_{rel} = 512$ . Simulations are not performed at $\alpha = \pi /2$ , as this results in an infinite relative Reynolds number.

Figure 18 presents the variation of $C_{L,q}$ , $C_{D,q}$ , $C_{y,q}$ and $C_{p,q}$ against $\alpha$ , for the OpU9 cylinder. The drag coefficient $C_{D,q}$ is largest at $\alpha = 0$ , when the cylinder is perpendicular to the oncoming flow. As $\alpha$ is increased, the profile is more streamlined, resulting in a substantial decrease in drag coefficient. The lift coefficient $C_{L,q}$ is zero at $\alpha = 0$ , due to symmetry. At large $\alpha$ , the cylinder profile is similar to a cambered airfoil, generating a substantial positive lift coefficient.

Figure 18. Variation of the (a) lift $C_{L,q}$ , (b) drag $C_{D,q}$ , (c) transverse $C_{y,q}$ and (d) power $C_{p,q}$ coefficients against inclination angle $\alpha$ under the quasi-steady assumption, obtained from simulations of a stationary OpU9 cylinder. The inclination angle was varied from $0\leqslant \alpha \leqslant ({7}/{16})\pi$ , and symmetry used to compute the force coefficients for negative $\alpha$ . Simulations were not performed at $\alpha = \pi /2$ , since the Reynolds number becomes infinite.

From (5.3), the lift force has a contribution in line with the cylinder’s direction of motion, while the drag force has a component opposed to the cylinder’s motion. For small $\alpha$ , the lift coefficient is small, and therefore the drag term dominates. This produces a negative power coefficient, opposing the cylinder motion. For moderate inclination angles ( $0.84\lt \alpha \lt 1.27$ ) the lift force is large, and dominates the drag term in (5.3), leading to a positive power coefficient, increasing the cylinder motion. Although the lift force remains large for $\alpha \gt 1.27$ , $\cos \alpha$ becomes small (i.e. the lift force is not sufficiently aligned in the transverse direction), so that the contribution of the lift force to the power coefficient decreases, and the power coefficient becomes negative.

Therefore, according to the quasi-steady analysis, there is a range of relative inclination angles ( $0.84\lt \alpha \lt 1.27$ ) where the component of the lift force in the transverse direction exceeds the drag force, leading to a positive instantaneous power coefficient, increasing the energy of the cylinder. This is due to the large lift force and small drag force occurring at large inclination angles, due to the airfoil-like shape.

To determine if galloping is possible, the power coefficient is integrated over a single oscillation cycle (Parkinson & Smith Reference Parkinson and Smith1964; Tucker Harvey et al. Reference Tucker Harvey, Khovanov, Murai and Denissenko2020). We assume the cylinder displacement is sinusoidal

(5.6) \begin{equation} y = A \sin \omega t\,,\,\dot {y} = A \omega \cos \omega t, \end{equation}

so that the instantaneous power coefficient is expressed as

(5.7) \begin{equation} C_{p,q} = C_{y,q}\big (\tan ^{-1}(A\omega \cos ( \omega t) )\big ) A \omega \cos (\omega t). \end{equation}

The cycle-averaged power coefficient is a function only of the combined parameter $A\omega$

(5.8) \begin{equation} \bar {C}_{p,q}(A\omega ) = \frac {A\omega }{2\pi } \int _{0}^{2\pi } C_{y,q}(\tan ^{-1}(A\omega \cos u)) \cos u \, \mathrm{d} u .\end{equation}

Figure 19(a) presents the variation of $\bar {C}_{p,q}(A\omega )$ against $A\omega$ for the OpU9 cylinder. For small-amplitude oscillations ( $A\omega \lt 1.29$ ), the net power added to the cylinder during one cycle is negative, and therefore the amplitude of oscillation will decrease to zero. For $1.29 \lt A\omega \lt 3.59$ , the cycle-averaged power is positive, so that the amplitude will increase until it reaches a value of $A\omega = 3.59$ , For $A\omega \gt 3.59$ , the cycle-averaged power is again negative, and the amplitude of oscillation will decrease.

Figure 19. (a) Variation of the cycle-averaged power coefficient for the quasi-steady model ( $\bar {C}_{p,q}$ ), against $A\omega$ , for the OpU9 cylinder. Dashed lines show the power extracted by mechanical damping ( $\bar {C}_{p,\xi }$ ), for either $m^*\xi /U^* = 0$ (undamped), or $m^*\xi /U^* = 4.5\times 10^{-3}$ (optimally damped). Red circles indicate stable limit cycles. (b) Comparison between the predicted maximum amplitude obtained from the quasi-steady model and numerical simulations.

For an undamped system, this indicates that there are two stable equilibria: a stationary cylinder ( $A\omega = 0$ ), and a stable periodic limit cycle with $A\omega = 3.59$ . There is also an unstable equilibrium at $A\omega = 1.29$ . This indicates hard-galloping behaviour, in which small-amplitude oscillations are stable, and will decay to zero. However, given a sufficient initial displacement ( $A\omega \gt 1.29$ ), the body will gallop with a large amplitude $A\omega = 3.59$ .

Assuming the system vibrates at the natural frequency $\omega = \sqrt {k/m} = 2\pi /U^*$ , the amplitude of oscillation predicted by the quasi-steady model is $A = 0.5714 U^*$ . Figure 19(b) compares this prediction with the numerical simulations, showing reasonable agreement, given the limitations of the quasi-steady model. However, the amplitude obtained from numerical simulations is lower than the predictions obtained from the quasi-static approximation.

5.3. Instantaneous force coefficients

The quasi-steady analysis presented in § 5.2 assumes that the hydrodynamic force applied to the cylinder depends only on the relative angle of attack $\alpha$ . However, the shedding of a primary vortex each half-cycle of motion (figure 15) will also affect the forces, which will not be captured by the quasi-steady model. In this section, we compare the force coefficients obtained from numerical simulations with the predictions of the quasi-static model, to deduce the relative significance of the quasi-steady forces and the vortex-induced forces to the motion.

Bouclin (Reference Bouclin1979) models the VIV–galloping interaction by including an unsteady fluctuating fluid force, representing the vortex-induced forces, in addition to the quasi-steady forces. Therefore, we interpret the difference between the instantaneous and quasi-steady force coefficients as representing the effects of vortex shedding on the instantaneous force and moment coefficients.

Figure 20 presents a comparison between the instantaneous lift ( $C_{\!L}$ ), drag ( $C_{\!D}$ ) transverse force ( $C_{\!y}$ ) and power ( $C_{\!p}$ ) coefficients for $U^* = 4$ , $6$ , $8$ and $10$ , to the quasi-static prediction. Here, the transverse and streamwise force coefficients, $C_{\!y}$ and $C_x$ , were obtained directly from the numerical simulations. Then, the added mass force was subtracted from the transverse force, to obtain the vortex force ( $C^{\prime}_y = C_{\!y} + 2 m^*_{a.m.} \ddot {y}$ , where $m^*_{a.m.}$ is the dimensionless added mass). To estimate the added mass, a numerically simulated pluck test in still fluid is performed using the same numerical set-up described in § 2.3, but with $U=0$ at the inlet. The cylinder is given an initial displacement, and the frequency of the subsequent oscillations is measured. The added mass is then estimated from the relationship $\omega = \sqrt {k/(m+m_{a.m.}}$ , where $\omega$ is the angular frequency of oscillations.

Figure 20. Comparison between the instantaneous (a) lift, (b) drag, (c) transverse force and (d) power coefficients obtained from numerical simulations at various $U^*$ and the quasi-steady model. Solid lines indicate $y\gt 0$ , while dashed lines indicate $y\lt 0$ .

The lift and drag coefficients were then computed according to

(5.9) \begin{align} C_{L} &= (C^{\prime}_y \cos \alpha + C_x \sin \alpha )/(1+\tan ^2 \alpha ), \end{align}

(5.10) \begin{align} C_{D} &= (C^{\prime}_y \sin \alpha - C_x \cos \alpha )/(1+\tan ^2 \alpha ). \\[9pt] \nonumber \end{align}

While the quasi-steady model captures the general behaviour of $C_{\!L}$ and $C_{\!D}$ , particularly for large inclination angles, substantial deviations are observed. This indicates that, while the quasi-steady model provides useful information about the effect of $\alpha$ on the lift and drag coefficients, other factors, such as vortex shedding, are also significant.

In figure 20, solid lines indicate $y\gt 0$ , while dashed lines indicate $y\lt 0$ . For $y\gt 0$ and $\alpha \gt 0$ (or $y\lt 0$ and $\alpha \lt 0$ ), the cylinder is in decelerating motion as it approaches the point of maximum displacement. In this region, the lift and drag coefficients are both slightly lower than the quasi-static prediction. However, the power coefficient is negative, indicating that the drag term dominates over the lift. For $y\gt 0$ and $\alpha \lt 0$ (or $y\lt 0$ and $\alpha \gt 0$ ), the cylinder is accelerating as it travels away from the point of maximum displacement. In this region, both $C_{\!L}$ and $C_{\!D}$ are substantially larger than the quasi-static prediction. The power coefficient is positive, indicating that the lift term dominates over the drag.

Note that $C_{\!y}$ , and therefore also $C_{\!p}$ , contains opposite-signed contributions from $C_{\!L}$ and $C_{\!D}$ (5.3). For the quasi-steady model, these contributions are close to equal for small to moderate $\alpha$ , so that the quasi-steady transverse force contribution $C_{y,q}$ is relatively small. Although the instantaneous lift and drag coefficients approximately follow the quasi-steady prediction, the instantaneous transverse force and power coefficient do not obviously follow the quasi-steady prediction, and are instead dominated by other unsteady contributions, such as vortex shedding.

Figure 20(b) shows a sharp decrease in the power coefficient when the inclination angle exceeds $\alpha = 1.27$ , for both quasi-steady and instantaneous power coefficients. This is due to the poor alignment between the lift force and the transverse direction at high inclination angles. Therefore, the quasi-steady force will limit the maximum inclination angle, which limits the maximum value of $A\omega$ .

Figure 21 presents the time history of the instantaneous force and power coefficients ( $C_{\!y}$ and $C_{\!p}$ ), as well as the quasi-steady predictions ( $C_{y,q}$ and $C_{p,q}$ ), against the phase angle $\phi = \omega t$ , where $\omega$ is the natural frequency. In general, there is little correlation between the instantaneous force coefficients and the quasi-steady predictions, indicating that the transverse force is dominated by vortex shedding or other unsteady effects.

Figure 21. Comparison between the instantaneous transverse force coefficient (left) and power coefficient (right) and the predictions of the quasi-steady model, for (a) $U^* = 4$ , (b) $U^* = 6$ , (c) $U^* = 8$ and (d) $U^* = 10$ .

For moderate and high $U^*$ ( $U^*\geqslant 6$ ), the motion appears to be driven by a large peak in both transverse force and power coefficient, which occurs immediately after the point of maximum displacement (between $0\leqslant \phi \leqslant \pi /4$ and $\pi \leqslant \phi \leqslant 5\pi /4$ ). As shown in figure 15, this corresponds to the shedding of a primary vortex from the cylinder. Therefore, there is a substantial contribution from VIV to the observed motion, which enhances FIVs.

Thus, we have the following interpretation of the motion at high $U^*$ . A primary vortex is shed each half-cycle, as the cylinder accelerates away from the point of maximum displacement, similar to the shedding of a starting vortex. This shed vortex results in a large transverse force in line with the cylinder velocity, adding energy to the cylinder which enhances FIV motion. As the cylinder travels downwards, there is a large inclination angle relative to the flow. This produces both a drag force, which opposes the motion, and a lift force, which enhances FIV motion. For $\alpha \lessapprox 1.27$ , the quasi-steady lift and drag terms are approximately equal, so the quasi-steady effect provides a negligible contribution to the transverse force. However, for $\alpha \gtrapprox 1.27$ , the component of lift in the transverse direction becomes small, so that the quasi-steady drag dominates. This effectively limits the maximum amplitude for a given frequency, so that the amplitude increases approximately linearly with $U^*$ .

This section has demonstrated that, while the quasi-steady model captures some effect of the relative angle of attack, unsteady forces are also substantial. Given the similarities between the crescent cylinders and airfoil/turbine blades, modelling approaches developed for dynamic stall (Bangga, Lutz & Arnold Reference Bangga, Lutz and Arnold2020; Mohamed et al. Reference Mohamed, Melani, Bangga, Aryan, Greco and Bianchini2024) may prove useful in providing further insight into the observed motion. These models include contributions from attached flow, separated flow and vortex shedding all of which are observed in the present flow.

6.1. Power generation according to the quasi-steady model

Assuming a sinusoidal displacement profile $y = A \sin \omega t$ , the power extracted from the cylinder via mechanical damping (2.8) is

(6.1) \begin{equation} \bar {C}_{p,\xi } = \frac {2\pi m^* \xi }{U^*} (A\omega )^2, \end{equation}

which depends on two parameters: the amplitude parameter $A\omega$ , and the mass-damping parameter $m^*\xi /U^*$ . Figure 19 plots $\bar {C}_{p,\xi }$ against $A\omega$ , for $m^* \xi /U^* = 4.5 \times 10^{-3}$ . When $\bar {C}_{p,\xi }\gt \bar {C}_{p,q}$ , the energy added to the cylinder by quasi-steady forces over one cycle is less than the energy extracted by mechanical damping, and so the amplitude will decrease. Likewise, the amplitude will increase for $\bar {C}_{p,\xi }\lt \bar {C}_{p,q}$ . Therefore, there is a stable equilibrium when $C_{p,\xi } = \bar {C}_{p,q}$ and $\partial \bar {C}_{p,\xi }/\partial (A\omega ) \gt \partial \bar {C}_{P,q}/\partial (A\omega )$ , indicated by the red circle in figure 19. The value of $m^* \xi /U^* = 4.5\times 10^{-3}$ corresponds to the maximum value of $\bar {C}_{p,\xi }$ at this stable equilibrium, and therefore the maximum power generation under the quasi-static assumption.

Figure 22 plots a comparison between the predicted power coefficient and vibration amplitude according to the quasi-steady model, and numerical simulations of the OpU10 cylinder, performed for various values of $U^*$ . The quasi-steady model predicts a maximum power coefficient of $\bar {C}_{P,\xi } =0.219$ , occurring at $m^* \xi /U^* = 4.5 \times 10^{-3}$ , and that galloping cannot be sustained for $m^* \xi /U^* \gt 5.2\times 10^{-3}$ .

Figure 22. Comparison of the (a) power coefficient $\bar {C}_p$ and (b) amplitude parameter $A\omega$ predicted from the quasi-steady analysis, and numerical simulations performed at various $U^*$ , for the OpU9 cylinder at $Re = 100$ and $m^* = 4.7124$ .

However, the numerical simulations demonstrate that significant FIVs are sustained for damping coefficients much larger than this, and yield significantly larger maximum power coefficients (up to $\bar {C}_{q,\xi } = 0.364$ ). This shows that the energy added to the cylinder by vortex-induced forces provide a significant contribution to power generation, allowing more power to be extracted via mechanical damping per cycle of motion.

Figure 23 plots the variation of power coefficient $\bar {C}_p$ , efficiency $\eta$ , amplitude $y_{\textit{max}}$ and angular frequency $\omega$ against $\xi$ , for the numerical simulations. Increasing $\xi$ initially increases power extraction, as more power is extracted by mechanical damping. However, damping also decreases the amplitude of vibration, leading to a reduction in the power coefficient for large $\xi$ . For $U^* \geqslant 7$ , the galloping-like behaviour cannot be sustained for $\xi \gtrapprox 0.025$ , leading to a sudden decrease in amplitude and power coefficient. Figure 23(d) shows that the oscillation frequency $\omega$ (solid lines) shifts away from the natural frequency (dashed lines) when this occurs, indicating a desynchronised behaviour. This suggests that the synchronised behaviour that occurs at high $U^*$ is due to the large oscillation amplitude, which significantly alters the vortex shedding. When the damping is sufficiently large to prevent the galloping-like behaviour, vortex shedding instead occurs at a different frequency (presumably the natural frequency of vortex shedding).

Figure 23. Variation of (a) power coefficient $\bar {C}_p$ , (b) efficiency $\eta$ , (c) oscillation amplitude $y_{\textit{max}}$ and (d) angular frequency $\omega$ , against damping coefficient $\xi$ , for numerical simulations of the OpU9 cylinder performed at $Re = 100$ and $m^* = 4.7124$ . Dashed lines in (d) indicate the undamped natural frequency $\omega _n = 2\pi /U^*$ .

At high $U^*$ , the power coefficient profiles display evidence of multiple different flow regimes (sudden jumps to different trend lines), depending on the precise value of $\xi$ . This results in sharp peaks in the power coefficient, meaning that $\xi$ must be precisely tuned to maximise the power coefficient. The overall maximum power coefficient, $\bar {C}_p = 0.364$ is obtained for $U^* = 7$ and $\xi = 0.0220$ . At $U^* = 6$ , there is a much broader peak in $\bar {C}_p$ with a maximum value of $\bar {C}_p = 0.314$ and $\xi = 0.0259$ . The broad peak means that $\bar {C}_p$ is less sensitive to $\xi$ , which will produce more consistent power generation.

The efficiency $\eta$ (2.9) represents the ratio of power extracted from power available in the swept area of the device. This is plotted in figure 23(b), against $\xi$ , for the OpU9 cylinder. The maximum efficiency of $8.25\,\%$ is obtained for $U^* = 6$ and $\xi = 0.0350$ , and the power coefficient at this condition is $\bar {C}_p = 0.289$ . The efficiency decreases as $U^*$ is increased. At $U^* = 7$ , when the power coefficient is maximised, the efficiency is $\eta = 6.71\,\%$ . The reduction in efficiency is due to an increased swept area with increasing $U^*$ .

Table 6. Shape parameters for various cylinders optimised for maximum power coefficient $\bar {C}_p$ , with fixed $U^* = U^*_{{opt.}}$ , $Re = 100$ and $m^* = 4.7124$ . The OpUM_P cylinder was optimised to maximise the average power coefficient across three different conditions ( $U^* = 6,8$ and $10$ ).

6.2. Optimising the power coefficient

We now consider shape optimisation to maximise the power coefficient, using the BO-GPR method. For this case, the objective function to be maximised is the mean power coefficient ( $r(\boldsymbol{q}) = \bar {C}_{p}$ ), and the design parameter vector $\boldsymbol{q}$ contains all shape parameters and the damping coefficient $\xi$ : $\boldsymbol{q} = (\textit{AR},\alpha ,\beta ,m_{\!f},n_{\!f},m_a,n_a,\xi )$ . The remaining details of the implementation are the same as in § 4. For each optimisation, $U^*$ is held constant, and the Reynolds number and dimensionless mass are $Re = 100$ and $m^* = 4.7124$ , respectively. Finally, the damping coefficient is restricted to the range $0\leqslant \xi \leqslant 0.2$ .

Table 6 lists the parameters of the optimised shapes for $U^* = 6$ , $U^* = 8$ and $U^* = 10$ , as well as the mean power coefficient. We use the label OpUX_P to indicate shapes optimised to maximise power coefficient at $U^* = \mathrm{X}$ . The optimised shapes are plotted in figure 24, as well as the OpU6, OpU8 and OpU10 cylinders. The power-optimised cylinder have the same general features as the amplitude-optimised cylinders; namely a thin crescent blade with sharp vertices. However, compared with the amplitude-optimised cylinders, the power-optimised cylinders have a greater curvature. The OpU6_P and OpU8_P cylinders have a mean camber of $x_{\textit{cam}} \approx 0.29$ , approximately thee times greater than the camber of the OpU6 and OpU8 cylinders. The OpU10_P cylinder has a mean camber of $x_{\textit{cam}} \approx 0.1897$ , which is approximately twice the camber of the OpU10 cylinder.

Figure 24. Comparison between cylinders optimised for maximum power coefficient (OpU6_P, OpU8_P, OpU10_P, OpUM_P) and cylinders optimised for maximum amplitude (OpU6, OpU8, OpU10).

Optimising at a single $U^*$ does not guarantee that the cylinder will perform well over a range of different operating conditions. An additional optimisation was performed using the average power coefficient across three different reduced velocities ( $U^*= 6$ , $8$ and $10$ ) as the objective function. This cylinder, labelled OpUM_P, is plotted in figure 24, and its parameters are provided in table 6. This cylinder has a similar shape to the OpU6_P and OpU8_P cylinders, but with a slightly higher camber.

To investigate the influence of the increased camber on the power generation, figure 25 compares the instantaneous and quasi-steady lift, drag and power coefficients for the OpUX_P to the OpU10 cylinder, at $U^* = \mathrm X$ . Increasing camber results in an increased maximum power coefficient according to the quasi-steady approximation, due to an increased lift coefficient. However, the increased drag coefficient reduces the maximum angle of attack for which the quasi-steady power coefficient is positive. Therefore, increasing camber reduces the maximum amplitude, but increases the power generation at smaller amplitudes.

Figure 25. Variation of the instantaneous and quasi-steady lift (a,d,g), drag (b,e,h) and power (c,f,i) coefficients at $Re = 100$ and $m^* = 4.7124$ : (a,b,c) OpU6_P at $U^* = 6$ and $\xi = 0.0681$ ; (d,e,f) OpU8_P at $U^* = 8$ and $\xi = 0.0802$ ; (g,h,i) OpU10_P at $U^* = 10$ and $\xi = 0.0495$ . Plots of the OpU9 cylinder with $\xi = 0.0220$ are also provided for comparison.

Increasing camber may also affect the vortex-induced forces. Figure 26 plots vorticity visualisations for the power-optimised cylinders OpU6_P, OpU8_P and OpU10_P. Vorticity visualisations for the OpU9 cylinder with $\xi = 0.0220$ are also provided for comparison. Animations are provided in the online supplementary material (movies 5–7). For the OpU6_P and OpU8_P cylinders, the increased camber alters the structure of vortices as they are shed from the cylinder, compared with the OpU9 cylinder. The resulting power coefficients (figure 25 c,f) are positive for all inclination angles, while the power coefficients for the OpU9 cylinder are sometimes negative. Therefore, changes to the vortex-induced forces as a result of cylinder camber also appear to enhance power generation.

Figure 26. Vorticity visualisations for the optimised cylinders (a) OpU6_P at $U^* = 6$ and $\xi = 0.0681$ ; (b) OpU8_P at $U^* = 8$ and $\xi = 0.0802$ ; (c) OpU10_P at $U^* = 10$ and $\xi = 0.0495$ . Vorticity visualisations for the OpU9 cylinder with $\xi = 0.0220$ at (d) $U^* = 6$ , (e) $U^* = 8$ and (f) $U^* = 10$ . Red indicates positive (counter-clockwise) vorticity, while blue indicates negative (clockwise) vorticity.

For the OpU10_P cylinder, the vortex structure is similar to the OpU9 cylinder (figure 26 c, f). However, figure 25(g–f) shows that the peak in power coefficient associated with vortex shedding occurs at a larger angle of attack (i.e. delayed later in the vortex shedding cycle), resulting in a much larger peak power coefficient, which enhances the power extracted from FIV.

Figure 27 plots the variation of power coefficient and efficiency against $U^*$ for the optimised cylinders. For the OpU6_P, OpU8_P and OpU_10 cylinders, there is a broad region of high power generation for reduced velocities $U^*\geqslant 6$ , with $\bar {C}_P\approx 0.4$ . Within this region, each cylinder achieves maximum power coefficient at its design $U^*$ , with a decrease in performance away from the design point. The maximum overall power coefficient of $\bar {C}_P =0.513$ is achieved with the OpU_10 cylinder, however, this cylinder performs worse than the OpU6_P and OpU8_P cylinders between $U^*=6$ and $U^* = 9$ . Maximum efficiency of $\eta \approx 12.2\,\%$ is achieved for $U^* = 6$ , and the efficiency decreases as $U^*$ increases beyond this value.

Figure 27. Variation of (a) power coefficient $\bar {C}_p$ and (b) efficiency $\eta$ against reduced velocity $U^*$ , for various optimised cylinders, at $Re = 100$ and $m^* = 4.7124$ .

In practice, devices must perform well over a range of operating conditions. The OpUM_P cylinder was optimised using the average power coefficient obtained at three different reduced velocities ( $U^* = 6$ , $8$ and $10$ ) as the objective function. As expected, the OpUM_P cylinder achieves a more consistent power coefficient for different reduced velocities, but has a lower maximum power coefficient than the single $U^*$ cylinders at their design condition.

As discussed in the introduction (table 1), previous studies have considered a range of different cylinder geometries. Since these studies are at much higher Reynolds numbers than the present study, we cannot compare with their findings directly. Instead, a series of optimisations were performed for several geometries listed in table 7. For these optimisations, the free parameters to be optimised are either the damping coefficient $\xi$ (circle, square, rectangle-triangle) or both $ \textit{AR}$ and $\xi$ (triangle, rectangle, ellipse, semi-ellipse). Of these shapes, the triangular and semi-elliptical cylinders achieve the best performance. This is consistent with prior studies (Sun et al. Reference Sun, Jo and Seok2019b ; Tamimi et al. Reference Tamimi, Armin, Shahvaghar-Asl, Naeeni and Zeinoddini2019; Huang et al. Reference Huang, Luo, Zhou, Wang, Xu and Qin2024) which find that triangular and semicircular afterbodies generate the most power.

Table 7. Optimised parameters $\xi$ and $ \textit{AR}$ for various elementary cylinder geometries at $Re = 100$ and $m^* = 4.7124$ . A sketch of each shape is provided, with the flow direction from left to right.

Figure 28. Variation of (a) power coefficient $\bar {C}_p$ and (b) efficiency $\eta$ against reduced velocity $U^*$ , for various cylinders at $Re = 100$ and $m^* = 4.7124$ .

Figure 28 compares the power coefficient and efficiency obtained using the circular, triangular and semicircular cylinders with the OpU10_P and OpUM_P cylinders. The optimised cylinders (OpU10_P and OpUM_P) achieve a substantially greater power coefficient over a much larger range of $U^*$ compared with the circular, triangular and semi-elliptical cylinders. The maximum power coefficient $\bar {C}_p = 0.513$ is obtained using the OpU10_P cylinder at $U^* = 10$ . This is 6.23 times larger than the maximum power coefficient obtained using the circular cylinder ( $\bar {C}_p= 0.0824$ ), and 1.94 times greater than that obtained using the semi-elliptical cylinder ( $\bar {C}_P = 0.264$ ). The OpUM_P cylinder has an average power coefficient of $\bar {C}_p = 0.404$ over the range $6\leqslant U^* \leqslant 11$ , which is 1.95 times the average power coefficient of the triangular cylinder over this same range ( $\bar {C}_P = 0.208$ ). The efficiency is maximised at $U^* = 6$ , with a maximum efficiency of $\eta = 12.2\,\%$ obtained using the OpUM_P cylinder, which is 2.52 times the efficiency of the circular cylinder ( $\eta = 4.84\,\%$ ), but only slightly greater than the efficiency of the triangular cylinder ( $\eta = 11.6\,\%$ ).

Note that we have not optimised to maximise the efficiency. Despite this, the OpU6_P, OpU8_P and OpUM_P cylinders achieve a better efficiency than other cylinders. For future studies, it may be worth performing a combined optimisation to balance maximising both power coefficient and efficiency.

6.3. Effect of maximum aspect ratio

Figure 29. (a) Variation of oscillation amplitude ( $y_{\textit{max}}$ ) against $U^*$ for modified OpU7 cylinders with various aspect ratios. (b) Variation of mean power coefficient $\bar {C}_P$ against $U^*$ for modified OpUM_P cylinders with various aspect ratios.

Given that the optimised cylinders have aspect ratio close to the upper bound of $ \textit{AR} = 20$ , this suggests that cylinders with a higher $ \textit{AR}$ may achieve even greater performance. Figure 29(a) presents the amplitude response for cylinders with parameters equal to those of the OpU7 cylinder, but with various aspect ratios. This confirms that the oscillation amplitude increases as $ \textit{AR}$ increases, with a zero-thickness blade being the ideal cylinder to maximise the amplitude response. In practice the maximum aspect ratio will be limited by engineering constraints, such as the strength and stiffness of the cylinder.

Figure 29(b) also plots the power coefficient for cylinders of different aspect ratios, with remaining parameters equal to the OpUM_P cylinder. While the maximum power coefficient increases with increasing $ \textit{AR}$ , for some reduced velocities the power coefficient decreases with increasing $ \textit{AR}$ . Given the optimised cylinders all have an $ \textit{AR}$ close to the upper bound of 20, it is likely that a zero-thickness blade will achieve the best overall power generation, however, this needs further investigation.