3.1. Optimization results

The results obtained from the optimization procedure are shown in figure 5, and indicate that the algorithm improves the desired wake characteristics of the porous cylinder arrangements with successive generations. Each symbol on the plots represents a unique arrangement of cylinders within the confined patch diameter $D$ as described in § 2.2, with the axes depicting the corresponding drag and enstrophy values. The drag values shown represent the average drag coefficient $C_d = F_{Drag}/(0.5\rho U_\infty ^2\times 9d)$ for each 9-cylinder array. We also note that the drag axis is inverted, i.e. the individuals are arranged such that the drag decreases along the positive vertical axis on the plots. This is necessary since the optimizer attempts to minimize wake enstrophy, while at the same time trying to maximize the drag acting on the structure, which are two mutually conflicting objectives, as can be observed from figure 5; in general, the wake enstrophy increases with increasing drag, and vice versa.

Figure 5. The evolution of a population of 128 distinct cylinder array arrangements with respect to the chosen fitness values, i.e. drag and average enstrophy, across successive generations: (a) generation 1, (b) generation 11, (c) generation 22, and (d) generation 33. The drag values shown represent the average drag coefficient for the arrays, $C_d = F_{Drag}/(0.5\rho U_\infty ^2 \times 9d)$. Each ‘$\times$’ symbol represents a unique arrangement of 8 cylinders around a central cylinder. The two most extreme Pareto-optimal individuals for each generation are shown as inset.

In figure 5(a), we observe a dense clustering of drag and enstrophy values for the initial 128 individuals, a subset of which is shown in figure 2. Although there is some variation in enstrophy, the drag coefficient is limited to approximately 0.7. Within the next 10 generations, the optimizer is able to discover cylinder arrangements that significantly increase drag, as can be seen in figure 5(b). This indicates that the arrangement of cylinders within a porous patch must be considered, in addition to patch porosity, for obtaining the desired wake characteristics. The influence of projected surface area and velocity flux through the interior of the array is explored at a later point.

Comparing the distribution of individuals from generations 11, 22 and 33 in figure 5, we do not observe a significant change in the overall distribution of the results, although there is an increase in the number of high-drag individuals in later generations. Given this convergence, the optimization procedure was stopped after 33 generations. Individuals that comprise the final Pareto front are identified in figure 5(d). This front consists of various distinct cylinder arrangements, all of which are optimal in their own right, i.e. no other solution is better than these Pareto-optimal individuals in both fitness metrics (i.e. lower enstrophy and higher drag simultaneously). Four of these Pareto-optimal individuals were selected for further analysis, namely, the two individuals at the extreme ends of the Pareto front, as well as two intermediate individuals. These individuals are depicted in figure 6, along with images of the physical models used for investigating the wake flow using the PIV experiments.

Figure 6. (a) Results from generation 33 of the optimization process. The insets show the Pareto-optimal arrangements that were selected for further analysis. (b) Images of the corresponding physical models used in the PIV experiments (5 mm diameter borosilicate glass rods).

3.2. Examining the wake flow

We now consider wake characteristics obtained using 2-D Navier–Stokes simulations for the four different arrangements shown in figure 6. Figure 7(a) depicts the array with the lowest average wake enstrophy (individual A), whereas figure 7(d) shows the array with the highest wake enstrophy (individual D). Qualitatively, we observe that there is an increasing tendency to form distinct high-intensity vortices as we go from figure 7(a) to figure 7(d). This trend is expected, given that higher enstrophy corresponds to higher vorticity magnitude. In addition to enstrophy, the configurations shown in figure 7 are also arranged in order of increasing drag, which was the second objective considered by the optimization procedure.

Figure 7. The wake generated by each of the four Pareto-optimal arrangements shown in figure 6 when a uniform inflow is imposed from left to right. The data were obtained using DNS, and the colours indicate vorticity. Corresponding animations are provided in supplementary movie 1, available at https://doi.org/10.1017/jfm.2023.255.

Considering the internal porous structure of these arrays, we note that the low-enstrophy low-drag case (individual A) has noticeable gaps within the array that resemble channels, which provide a path for the flow to pass through the interior of the array, as observed in the visualization of speed shown in figure 8(a). This has the dual impact of reducing drag, as well as disrupting the vortices shed in the wake, which in turn reduces the average wake enstrophy. Individual B in figure 8(b) also has small gaps that allow flow to pass through the porous structure, but to a lesser extent. Additionally, we observe a low-speed region in front of the array, which indicates that a majority of the freestream flow gets redirected around the array instead of passing through it. We observe that the high-drag high-enstrophy cases (figures 8c,d) have cylinders that display comparatively tighter clustering, with minimal streamwise flow allowed through the arrays. Quantitative details regarding the impact of the arrays’ internal porous structure are discussed at a later point.

Figure 8. Flow speed normalized with respect to the freestream velocity ($\lVert \boldsymbol {u}\rVert /U_\infty$). The data were obtained using DNS and time-averaged over one shedding time period for each of the four selected individuals A–D. Localized high-speed regions are visible as bright spots in the arrays’ interior.

Figure 9 shows velocity data obtained using PIV experiments in the wakes of the four Pareto-optimal arrangements. The variation of the normalized streamwise velocity ($u/U_{\infty }$) along a line cut is shown in the figure, plotted against the normalized downstream distance from the rear edge of the clusters (i.e. $x/D$). This line cut (also referred to as the midline here) was selected to be the straight line in the streamwise direction that passes through the array's centre. The general trend in figure 9 indicates a reduction in the streamwise velocity a short distance behind the cluster, after which there is a gradual recovery up to a maximum value that is lower than the freestream value. Furthermore, we observe that the most noticeable drop in $u/U_{\infty }$ occurs for individual D, which experiences the highest flow-induced drag as well as the highest wake enstrophy. At the same time, individual A with the minimum drag and lowest wake enstrophy displays the least severe drop. Individuals B and C are also shown for comparison, but neither is an extreme example with regard to drag or enstrophy. We observe that the recovered streamline velocity is highest for individuals A and C at $0.72U_{\infty }$, and lower for B at $0.62U_\infty$ and D at $0.58U_\infty$. We also observe that the minimum for individual D occurs farther downstream compared to the other arrays at $x=2.5D$, with velocity drop $-0.35U_\infty$. The minima for B and C occur at approximately $x=1.6D$, and display similar velocity drops $-0.25U_\infty$. These negative values are indicative of the formation of strong recirculation regions in the wakes of these arrays. The velocity profile for individual A does not display a distinctive minimum, but instead we observe a comparatively long region with a moderate velocity drop $-0.12U_\infty$. Overall, the most prominent difference among the velocity profiles is between that of individuals A and D, with the profiles for individuals B and C not differing significantly from that of A. These observations indicate the need for utilizing other aspects besides a one-dimensional line cut in order to fully characterize the arrays’ performance characteristics.

Figure 9. Normalized streamwise velocity along the midline for the four selected Pareto-optimal arrangements. The velocity measurements shown here were obtained experimentally using PIV, and were time-averaged over 8 vortex-shedding cycles.

Figure 10 shows the time-averaged streamwise velocity along the midline from DNS data for the four individuals. We observe a decreasing trend in far-wake velocity going from individual A to individual D, with A being the highest at 0.9$U_{\infty }$, followed by individual B at 0.8$U_{\infty }$, C at 0.75$U_{\infty }$, and finally D at 0.5$U_{\infty }$. Individual A has the largest drop in velocity in the region behind the cluster ($-0.2U_{\infty }$); however, velocity recovery occurs rapidly. The recovery in velocity slows down noticeably as we go from configuration A to D, i.e. with increasing drag and enstrophy. Comparing figures 9 and 10, we observe that the largest drop in velocity occurs closer to the arrays in the simulations compared to the experiments, in addition to faster recovery observed for the simulations. Several potential reasons may explain the differences observed between velocity profiles obtained from experiments and simulations, and these are discussed at a later point. However, the near-wake velocity field from the DNS was confirmed to match 2-D simulation results from a separate open-source solver (OpenFOAM), which uses finite volume methods and body-fitted meshes instead of the vortex methods and Brinkman penalization approach adopted in the present work.

Figure 10. Normalized streamwise velocity along the midline from DNS data, time-averaged over 8 vortex-shedding cycles. The line types correspond to those shown in figure 9.

It is expected that lower flow speeds will favour the sedimentation of suspended particles by increasing the time available for gravitational settling. However, examining only the streamwise component along a one-dimensional line cut may not provide a complete picture of the sedimentation process. For instance, in experiments conducted by Chen et al. (Reference Chen, Ortiz, Zong and Nepf2012), porous cylinder arrangements with higher overall midline streamwise velocity were found to be more conducive to soil deposition, which seems contrary to the expected behaviour. As discussed earlier, increased drag leads to higher velocity deficit and energy dissipation in the wake (Maza et al. Reference Maza, Adler, Ramos, Garcia and Nepf2017; Gijón Mancheño et al. Reference Gijón Mancheño, Jansen, Uijttewaal, Reniers, van Rooijen, Suzuki, Etminan and Winterwerp2021), which can promote gravitational settling. However, this process is also accompanied by increased turbulent kinetic energy production (Chen et al. Reference Chen, Ortiz, Zong and Nepf2012; Tinoco & Coco Reference Tinoco and Coco2018), which in turn may hinder sediment deposition and promote resuspension. Thus it is the combined influence of drag and wake enstrophy over the entire wake volume that determines particle sedimentation levels. It is likely that for the cases explored here with higher drag, high wake enstrophy might make the configurations unfavourable for sediment deposition (Chen et al. Reference Chen, Ortiz, Zong and Nepf2012; Norris et al. Reference Norris, Mullarney, Bryan and Henderson2019) even though the streamwise wake velocity is notably lower than $U_\infty$. On the other hand, for the cases with low enstrophy, the overall speed in the wake may be higher due to low drag, which is also unfavourable for gravitational settling (Yamasaki et al. Reference Yamasaki, de Lima, Silva, Preza, Janzen and Nepf2019). This implies that a balance must be struck between drag and wake enstrophy to promote sediment deposition.

In addition to the streamwise line cuts shown in figures 9 and 10, cross-stream line cuts for the time-averaged streamwise velocity $u$ were examined using both PIV and DNS data. The resulting comparison is shown in figure 11. The plots show reasonably good agreement between data from the experiments and 2-D simulations, and the differences observed may be due to several potential reasons, such as sensitivity to array orientation and relative positioning of the cylinders in the experiments, flow blockage within the experimental channel, and three-dimensional effects that are absent from the 2-D simulations, to name a few. The time-averaged cross-stream velocity profiles shown in figure 11 indicate an increase in velocity deficit going from A to D, and comparable deficits for individuals B and C. We observe that individual D displays the largest reduction in velocity, although this deficit is confined to a relatively narrow cross-stream region. In comparison, individuals B and C experience a less severe reduction in streamwise velocity, albeit a larger cross-stream area is affected, as indicated by the broader profiles. In addition to the velocity profiles, the Strouhal number ($St=fD/U_\infty$) for all four arrays was also computed using both PIV and DNS data, and a comparison is shown in figure 12. We observe good agreement between the values of $St$ obtained from the experiments and simulations. We note that it is difficult to compute enstrophy from PIV data since velocity gradients cannot be computed with high accuracy due to spatial averaging involved in the PIV analysis.

Figure 11. Cross-stream velocity profiles for the selected Pareto-optimal individuals A–D (a–d). Profiles obtained using DNS are shown as solid lines, whereas those obtained using PIV are shown as dash-dotted lines. The profiles were computed using time-averaged streamwise velocity over at least 8 shedding cycles, at a distance $5D$ downstream from the arrays’ rear edges.

Figure 12. Comparison of the Strouhal number ($St = fD/U_\infty$, where $f$ is the frequency of vortex shedding, and $U_\infty$ is the freestream velocity) between the PIV experiments ($\times$) and DNS ($\bullet$). The symbols represent the mean over 8 shedding cycles, and the error bars represent the standard deviation.

To estimate the drag acting on the arrays using the PIV data, a control-volume analysis was performed using the velocity profiles shown in figure 11. A simple schematic of the control-volume set-up is shown in figure 13(a). The average drag force on the arrays was determined as follows (assuming unit span):

(3.1a)\begin{align} F_{Drag} &= \left(-\unicode{x222F} (\rho \boldsymbol{u}) \boldsymbol{u} \boldsymbol{\cdot} \hat{\boldsymbol{n}} \,{\rm d}A \right) \boldsymbol{\cdot} \hat{\boldsymbol{i}} \end{align}

(3.1b)\begin{align} &= 4D\rho U_\infty^2 - \int_{{-}2D}^{2D}\rho\,u^2(y)\, {\rm d}y. \end{align}

Figure 13. (a) Schematic of the control volume used for estimating the drag force on the porous arrays using PIV data. The outflow velocity profile $u(y)$ was measured at a distance $5D$ downstream of the array edges, and the corresponding profiles are shown in figure 11. (b) Drag coefficient computed using control-volume analysis, with the drag force determined using (3.1).

The corresponding drag coefficient for each 9-cylinder array was then calculated as $C_d = F_{Drag}/(0.5\rho U_\infty ^2 9d)$, and the resulting values are shown in figure 13(b). We observe an increasing trend in $C_d$ going from configuration A to B to D; however, there is a slight decrease for configuration C. The value of $C_d$ for configuration D is considerably higher than that for A, and slightly higher than (but close to) that for B. We note that the trend for $C_d$ computed from DNS using (3.1) displayed the same trend as that observed in figure 13(b), but the magnitude tends to differ between the 2-D DNS and PIV. This can be explained by the differences observed between the cross-stream velocity profiles in figure 11, since these profiles form the basis of the control-volume based $C_d$ calculations. We remark that the lower value of $C_d$ for array C compared to array B in figure 13(b) is due to the control-volume analysis being done in the steady shedding state, which is necessary due to the use of time-averaged velocity profiles. This is also the reason why the far-wake velocity for configuration C recovers to a higher value than for B in figure 9. The $C_d$ computed during optimization (figure 6a) included the unsteady startup phase (i.e. the initial transient), and the subsequent time-averaging resulted in lower drag for array B than for array C (see Appendix A). During optimization, it is not feasible to predict beforehand when a specific array configuration might reach a steady shedding state, or what the resulting shedding frequency might be, making it necessary to use time averaging over a long duration. We note that the very early stages of the startup can be dependent on numerical ramp-up characteristics. To prevent any non-physical effects that might arise from this ramp-up, the initial non-dimensionalized time from 0 to $4t^*$ was excluded from all DNS runs for calculating the time-averaged $C_d$.

3.3. Flow within the porous arrays

The differences in drag and enstrophy values associated with each of the four selected individuals ultimately arise from differences in how the flow behaves around and within the clusters. As mentioned in some of the studies discussed earlier (e.g. Ricardo et al. Reference Ricardo, Sanches and Ferreira2016), vorticity cancellation due to the interference of neighbouring cylinders plays a prominent role in determining wake characteristics. Thus we take a closer look at flow patterns that develop within the porous structures for each of the four selected Pareto-optimal individuals, and examine the resulting impact on enstrophy and drag. The evolution of the vorticity field over one shedding time period is shown in figure 14 for the four selected arrays. One notable observation for individual A is the presence of cohesive structures (i.e. internal boundary layers) that occupy the array's interior, and remain persistent throughout the shedding period without undergoing notable changes in form and strength. These structures exhibit a tendency to align in the streamwise direction and appear to ‘attach’ to neighbouring cylinders. This is likely related to the relative closeness of adjacent cylinders and their specific arrangement within the array. Such internal cohesive structures are found less frequently as we move from individual A to individual D, where we instead observe structures that are less persistent and undergo rapid changes in shape, direction and intensity over time. The changes are most prominently observed for individual D, i.e. the high-drag high-enstrophy configuration.

Figure 14. Time evolution of vorticity over one shedding time period from DNS data, with the corresponding colour bar provided in figure 7. Panels (a) to (d) represent the Pareto-optimal individuals A to D, respectively, in order of increasing drag and enstrophy.

For further analysis, vorticity data collected from the numerical simulations were decomposed into constituent modes using proper orthogonal decomposition (POD), also referred to as principal component analysis. This allows us to isolate the most energetic flow structures (i.e. POD modes) in the interior and the wake of the cylinder arrays. We note that a POD mode's ‘energy’ is not related to the concept of potential or kinetic energy, but rather indicates how much variation in the data can be accounted for using a particular mode. Figure 15 shows the most energetic modes, which contain the highest variation in vorticity over time. Mathematically, these modes can be described as the first eigenvector (i.e. that associated with the largest eigenvalue) of the covariance matrix constructed using the flow snapshots.

Figure 15. The first POD mode of vorticity for the four selected arrangements obtained from DNS data, using columns of $\boldsymbol {V}$ from (3.3). The first modes from all four cases were observed to account for at least 33 % of variation in the data. Note that the colours do not necessarily correspond to positive or negative values of vorticity, since they represent the eigenvector of the covariance matrix, with Euclidean norm equal to 1.

A total of 100 vorticity snapshots represented by $\omega (x,y,t)$, spanning one shedding time period in the steady state, were used for each of the individuals shown in figure 15. Each such 2-D snapshot was rearranged into a one-dimensional row vector of length $M$ (the number of grid points) containing the vorticity values, and the 100 frames were arranged sequentially in rows. The time-averaged vorticity $\bar {\omega }(x,y)$ was subtracted from each frame in order to form the sample matrix

(3.2)\begin{equation} \boldsymbol{X} = \begin{bmatrix} \boldsymbol{x}(t_{1}) \\ \boldsymbol{x}(t_{2}) \\ \vdots \\ \boldsymbol{x}(t_N) \end{bmatrix}_{N\times M}, \end{equation}

where $\boldsymbol {x}(t)$ is the row vector $\boldsymbol {x}(t) = \omega (x,y,t) - \bar {\omega }(x,y)$, and $N=100$ is the number of snapshots. The eigenvectors $\boldsymbol {v}_i$ of the covariance matrix $\boldsymbol {X^{\rm T}X}$ correspond to the POD modes, while the eigenvalues $\lambda _i$ denote the contribution of each mode towards the total variation in the data about the time mean. These were obtained using singular value decomposition of the matrix $\boldsymbol {X}$, such that

(3.3)\begin{equation} \boldsymbol{X} = \boldsymbol{U \boldsymbol{\varSigma} V^{\rm T}} \end{equation}

where $\boldsymbol {U}$ and $\boldsymbol {V}$ are the left- and right-singular matrices of $\boldsymbol {X}$, and $\boldsymbol {\varSigma }$ is a diagonal matrix consisting of the singular values. The columns of $\boldsymbol {V}$ are the eigenvectors $\boldsymbol {v}_i$ of $\boldsymbol {X}^{\rm T} \boldsymbol {X}$. The elements $\sigma _{i}$ of $\boldsymbol {\varSigma }$ are arranged in descending order of magnitude, and are related to the corresponding eigenvalues as $\lambda _{i} = \sigma _{i}^{2}$. Thus the POD modes are obtained in order of their contribution to variation in the data, from highest to lowest.

The first POD modes obtained in this manner for all four selected configurations are shown in figure 15, and they highlight regions that correspond to the most significant variation in the vorticity field over one shedding time period. We note that any persistent vorticity structures that maintain a consistent form and strength in time are not expected to appear in these modes, since they get subtracted along with the time-averaged vorticity upon forming the sample matrix $\boldsymbol {X}$. For the sake of convenience, we refer to the vortices shed by the arrays as a whole as primary vortices, and those shed by individual cylinders as secondary vortices. For individual A, we observe a prominent region of high variation just downstream of the array, which is caused by consistently shedding primary vortices of opposite sign, as observed in figure 14(a). The other regions of high variation correspond to oscillatory internal jets (secondary vorticity) emerging from the rear and top of the array. These jets exist independently of the formation of the primary vortices to some extent, since they are driven by streamwise flow fed in through the array. At the other extreme, the first POD mode for individual D displays a more diffuse distribution, with fewer isolated regions of high variation. This indicates that a large spatial region of the flow for this array undergoes considerable variation in time. Moreover, it is evident from figure 14(d) that the formation of the secondary vortices is driven primarily by the shedding of the two primary vortices along the top and bottom of the array. Minimal internal flow through array D does not allow for the formation of independently evolving secondary vortices (and the resulting jets), as in the case of array A. The POD modes for individuals B and C display behaviour that lies in between these two extremes, with array B displaying high variation streamwise jets forming via internal flow near the top and bottom of the array, and array C forming a combination of internal flow-driven and primary-vortex-shedding-driven secondary vortices.

The observations discussed here suggest that for a given number of cylinders, arrays that produce lower enstrophy have arrangements that allow flow to pass through the interior to some extent in the streamwise direction. Moreover, these cylinders are arranged such that the boundary layers generated by upstream cylinders attach more readily to downstream neighbours. On the other hand, the clusters with higher drag consist of arrangements that are aligned in the cross-stream direction, which increases the projected frontal area. The corresponding arrangements also allow minimal flow through the porous structure in the streamwise direction (figure 8), resulting in shedding patterns that more closely resemble those of rigid impervious objects.

To quantify the observations regarding the dependence of drag on internal flow through the arrays and cylinder alignment in the cross-stream direction, we compute the internal flux and the projected surface area for the four selected arrays. Additionally, the representative spacing for each array is also calculated as the average of the radial distances of the 8 cylinders from the central cylinder. The average internal flux was calculated from DNS data using the time-averaged speed over approximately 8 shedding time periods, and defining the convex hulls $\varOmega$ manually for each array:

(3.4)\begin{equation} \psi = \dfrac{\displaystyle\iint_{\varOmega} \lVert\bar{\boldsymbol{u}}(x,y)\rVert \,{\rm d}A}{U_\infty(A_\varOmega-9{\rm \pi} d^2/4)} , \end{equation}

where $A_\varOmega$ is the area of the convex hull. The resulting plots are shown in figure 16 for the average spacing $r_{avg} = \sum _{i=1}^n r_{i}/n$ normalized by cylinder diameter $d$, projected frontal area $S$ normalized by $d$ (assuming unit span), and the non-dimensional internal flux $\psi$. We observe that there is no clear dependence of drag on average spacing, with $r_{avg}$ for individual D being smaller than values for B and C, although array D experiences the largest amount of drag (figure 13b). On the other hand, the projected area seems to correlate better with drag, with drag increasing for higher values of $S/d$. We note that the projected area for individual D is slightly lower than that for individual B, although its drag is slightly higher than that of B. This will be examined further using the internal flux. The plot of internal flux indicates that configuration A allows the most flow to pass through the porous array, and amounts to approximately 42 % of the freestream flux through a comparable area. In general, internal flux tends to decrease with increasing drag, except for configuration C, which was observed to have the minimum value $\psi =0.18$ despite experiencing lower drag than D. This indicates that drag may be dependent on a combination of internal flux and the projected frontal area, since the value of $S/d$ is smaller for configuration C than for D. Similarly, individual B ($\psi =0.30$) is observed to allow more flow to pass through the interior than individual D ($\psi =0.25$), which may explain why the drag experienced by B is slightly lower than by D, despite having a slightly larger projected area.

Figure 16. (a) Measure of average array spacing, calculated as the mean radial distance of the cylinders from the central cylinder. (b) Projected surface area $S$ for the arrays normalized by cylinder diameter $d$ (assuming unit span). (c) Normalized internal flux through the arrays, calculated from DNS data using (3.4).

3.4. Sensitivity analysis

To gain further insight into the role that the positioning of the cylinders plays in determining drag and enstrophy, we conduct a sensitivity analysis for the four selected Pareto-optimal cases. The position of each cylinder within an array was perturbed in the radial and azimuthal directions by small amounts, hence changing the corresponding $r_{i}$ and $\theta _{i}$ values one cylinder at a time. Cases where this resulted in a collision with neighbouring cylinders, or in exceeding the maximum patch diameter $D$, were not considered. The enstrophy and drag were calculated for all the perturbed cases, and detailed results for one of the scenarios are shown in figures 17 and 18. Additionally, the most- and least-sensitive cylinders with regard to change in enstrophy and drag are shown in figure 19 for each of the four arrays.

Figure 17. Percentage change in enstrophy (purple squares) and drag (green circles) when the position of one cylinder in the high-drag arrangement (individual D) was changed with respect to $r$ ($\Delta r=d/4$) in DNS. (a–d) Data for cylinders 1, 2, 6 and 7, respectively, which are labelled in (e). The remaining cylinders occupy positions where they could not be perturbed in the radial direction without encountering collisions. The most sensitive (dark red) and least sensitive (light red) cylinders with respect to changes in $r$ are highlighted in (e).

Figure 18. Percentage change in enstrophy (purple squares) and drag (green circles) when the position of one cylinder in the high-drag array (individual D) was changed with respect to azimuthal angle $\theta$ ($\Delta \theta = {{\rm \pi} }/{32}$) in DNS. (a)–(h) Data for cylinders 1-8, respectively, as shown in (i). The most sensitive (dark blue) and least sensitive (light blue) cylinders with respect to changes in $\theta$ are highlighted in (i).

Figure 19. The most and least sensitive cylinders for each of the four selected arrangements. Sensitivity with respect to enstrophy is examined in (a–d), and that with respect to drag is examined in (e–h). The dark red cylinders are most sensitive to changes in $r$ (radial direction), and the dark blue cylinders are most sensitive to changes in $\theta$ (azimuthal direction). The lighter coloured cylinders are least sensitive with respect to $r$ (light red) and $\theta$ (light blue). The values next to the cylinders indicate the absolute percentage changes in drag or enstrophy upon perturbation.

Figure 17 shows the percentage change in drag and enstrophy for individual D when each cylinder's position was changed slightly by $\pm \Delta r$ and $\pm 2\Delta r$ (where $\Delta r = d/4$) in the radial direction with respect to the central cylinder. Similarly, figure 18 shows the percentage change in drag and enstrophy when each cylinder's position was changed by $\pm \Delta \theta$ and $\pm 2\Delta \theta$ (where $\Delta \theta = {{\rm \pi} }/{32}$) in the azimuthal direction about the centre of the array. The cylinders showing the highest and lowest sensitivity overall (i.e. for drag as well as enstrophy) have been highlighted in the respective figures. As observed in figure 17, only the cylinders labelled 1, 2, 6 and 7 were able to be perturbed in the radial direction. The largest change in enstrophy was observed for cylinder 2, when its radial distance was increased by $2\Delta r$. The largest change in drag was also observed for cylinder 2 at this perturbation distance. Hence cylinder 2 is labelled as the most sensitive cylinder overall for individual D with respect to radial perturbation (figure 17e). On the other hand, cylinder 6 exhibits the smallest maximum change in drag, but has a maximum change in enstrophy comparable to that of cylinder 1. Hence cylinder 6 is labelled as the overall least sensitive cylinder for individual D with respect to radial perturbations. Similarly, in figure 18, cylinder 2 was found to have the largest maximum change in enstrophy as well as drag when perturbed in the azimuthal direction, whereas cylinders 7 and 8 both showed similar and overall smallest changes in enstrophy and drag (figure 18i). We note that none of the perturbed array configurations for the four selected individuals resulted in higher drag and lower enstrophy simultaneously, compared to their respective base unperturbed configurations. More specifically, none of the perturbed individuals was better than the original Pareto-optimal individuals in both fitness metrics, which indicates that the Pareto-optimal solutions obtained using the multi-objective optimization are robust.

In figure 19, the sensitivity of the four selected Pareto-optimal individuals is considered separately with regard to drag and enstrophy. The cylinders that have the highest and least impact on each of these two aspects upon perturbation are highlighted. Figures 19(a–d) depict sensitivity with regard to enstrophy, and figures 19(e–h) depict sensitivity with respect to drag. One notable observation is that across all the examined cases, the perturbations tended to have a greater influence on wake enstrophy than on drag. This indicates that enstrophy is more sensitive than drag to specific cylinder placement within porous arrays.

For individual A in figure 19(a), we observe that the cylinders most sensitive with respect to enstrophy are the ones whose perturbation disrupts boundary layer attachment among neighbouring cylinders (as observed in figure 14a), leading to an increase in wake enstrophy. The dark blue cylinder in figure 19(a) blocks flow from passing through the interior of the array when perturbed in the azimuthal direction, which also increases enstrophy due to weaker secondary vortices (and hence weaker vorticity cancellation). The phenomenon of vorticity cancellation, especially at the low local Reynolds number in the interior of the arrays, is similar to that observed by Moulinec, Hunt & Nieuwstadt (Reference Moulinec, Hunt and Nieuwstadt2004), who found that oppositely signed vorticity in close proximity led to cancellation via diffusion. On the other hand, the least sensitive cylinder in figure 19(a), when perturbed in the azimuthal direction, was observed to retain its tendency to generate boundary layers that remain attached to neighbouring cylinders. Examining the drag sensitivity for individual A, we observe in figure 19(e) that perturbing the dark red cylinder in the radial direction increases drag by $25\,\%$. This may be related to an overall increase in the effective frontal surface area as the cylinder in question is moved radially outwards. Moreover, increasing its distance from the centre disrupts the continuity of the attached boundary layer observed in figure 14(a), thereby increasing drag. The same cannot be said for the rest of the cylinders for individual A, since even after perturbation they remained sufficiently close to their neighbours so as not to disrupt the continuous boundary layer. Meanwhile, the light blue cylinder caused no significant variation in drag when perturbed in the azimuthal direction, and it did not induce any major changes to the overall flow pattern, presumably because the internal channels that allow flow to pass through the array remain mostly unchanged.

For individual B in figures 19(b) and 19( f), we note that a majority of the cylinders, especially the outermost ones, are restricted in perturbations with respect to $r$; positive perturbations would cause cylinders to exceed the outer patch diameter $D$, whereas negative perturbations would lead to collisions. With respect to drag (figure 19f), one cylinder in particular affects drag significantly when perturbed radially ($29\,\%$), and minimally when perturbed in the azimuthal direction ($4\,\%$). This can be explained by the influence the respective perturbations have on internal flow through the array. Upon decreasing the radial distance of the cylinder in question, it completely blocks off internal flow through the lower half of the array. As a result, the freestream gets redirected around the lower half, behaving more closely to flow around a cylinder. This reduces the interaction of the oncoming flow with the internal cylinders, most of which are placed in the cross-stream direction, resulting in lower drag. We note that azimuthal perturbations for this particular cylinder do not have as significant an impact on the internal flow, resulting in only a $4\,\%$ change. The remaining light blue cylinders in figure 19( f) obstruct the incoming freestream, and the extent of their blockage is not affected by perturbations in $\theta$. We also note that one of the cylinders in the array is minimally relevant for enstrophy, but highly relevant for drag ($1.8\,\%$ versus $31\,\%$), and another cylinder for which the reverse is true ($28\,\%$ versus $5\,\%$). If we consider the vorticity shown in figure 14(b), then the cylinder coloured light red in figure 19(b) is isolated from the other cylinders in the array, and generates a highly unsteady shedding pattern of its own. Hence, even though it changes the flow structure slightly upon perturbation, the structures it affects have a small impact on the overall enstrophy ($1.8\,\%$).

For individual C in figures 19(c) and 19(g), we observe one particular cylinder that has a large impact on enstrophy ($43\,\%$) and a minimal impact on drag ($4\,\%$). This is because the cylinder is in the vicinity of the most consistent internal flow structures that form for this particular array. Increasing the azimuthal angle will hinder the formation of the opposing-vorticity boundary layers found at an angle of approximately ${\rm \pi} /4$ near the top of the array in figure 14(c). However, since the dominant direction of local flow is tangential to the cylinder, its overall contribution to drag does not change significantly upon azimuthal perturbation. In figure 19(g), we observe that the light blue cylinders do not affect drag significantly upon perturbation in the $\theta$ direction, for the same reason as that for individual B, i.e. perturbations in the azimuthal direction do not result in significant changes to the overall flow. The dark red cylinder in 19(g) causes a decrease in drag as its distance from the centre increases. This may be related to the fact that it gets better aligned with the wakes of the cylinders upstream of it, which is also accompanied by a significant reduction in enstrophy (30 % reduction, not annotated on the figure). This reinforces an earlier observation that boundary layers of neighbouring cylinders staying attached is a significant contributor to reducing enstrophy.

Individual D (figures 19d,h) appears overall to be less sensitive to perturbations than the other three individuals, since it has a comparatively large number of lighter coloured cylinders, and the absolute percentage change values are smaller. This may be because the localized flow structures are highly unsteady (figure 14d), so that perturbations to any one specific cylinder will play a smaller role in modifying the overall flow structure. We observe that the same cylinder is most sensitive towards drag (azimuthal direction, 19 %) and towards enstrophy (both radial and azimuthal directions, 17 % and 18.5 %). Decreasing the azimuthal angle for this cylinder would connect it to the lower neighbouring cylinder such that their boundary layers remain attached, thereby reducing both drag and enstrophy. In figure 19(h), the single light red cylinder has a minimal impact on drag as well as on enstrophy (7 %) upon radial perturbation, since it is isolated and does not impact the resultant flow significantly.

We note that while we aimed to perform sensitivity analysis for all of the cylinders in each array, perturbing certain cylinders would have resulted in an overlap or collision with neighbouring cylinders, thus these cylinders were excluded. Furthermore, figure 19 highlights only those cylinders that displayed the highest and lowest sensitivity for a given configuration, although several of the unmarked cylinders were also examined. Overall, the observations presented here indicate that the internal arrangement of cylinders in an array is an important consideration for controlling drag and enstrophy. The sensitivity analysis reveals that enstrophy is minimized by directing incoming flow to pass within the array, which leads to vorticity cancellation (Moulinec et al. Reference Moulinec, Hunt and Nieuwstadt2004; Ricardo et al. Reference Ricardo, Sanches and Ferreira2016) and the suppression of individual vortex shedding due to the attachment of boundary layers among neighbouring cylinders. We also observe that drag appears to be dependent on the projected frontal area, as well as on boundary layer attachment among neighbouring cylinders to a smaller extent. Variation in drag based on cylinder placement was also observed in the study by Shan et al. (Reference Shan, Liu and Nepf2019), who found that random physical arrangements resulted in higher drag compared to uniform grid arrangements. These observations indicate that a high-level metric such as porosity may not be sufficient on its own to predict wake characteristics for porous arrays.

3.5. Particle tracking and Lagrangian coherent structures

We now use Lagrangian particle tracking to examine the trajectories of tracer particles seeded into the flow to determine the behaviour of each of the selected Pareto-optimal arrays with regard to particle sedimentation and erosion. We note that sediment transport is a complex process that depends on several flow and sediment properties, particularly the bed shear stress (López & García Reference López and García1998; Kothyari, Hashimoto & Hayashi Reference Kothyari, Hashimoto and Hayashi2009). However, tracers can provide a useful approximation of particle transport around and within the porous arrays. We consider the initial transient state separately from the steady vortex shedding state, to explore how particle transport may differ in systems involving intermittent flow (e.g. wave-dominated coastal regions) from those involving continuous flow (e.g. rivers). In both scenarios, particles were initialized separately upstream of the porous arrays and in the cylinders’ immediate vicinity and wake (see figure 28). Apart from particle tracking, Lagrangian coherent structures (LCS) associated with steady shedding state were computed using the finite time Lyapunov exponent (FTLE) (Haller & Yuan Reference Haller and Yuan2000; Shadden, Lekien & Marsden Reference Shadden, Lekien and Marsden2005; Lekien, Shadden & Marsden Reference Lekien, Shadden and Marsden2007). Ku & Hwang (Reference Ku and Hwang2018) used FTLE analysis to study the transport and dispersion of small particles in stratified flows, demonstrating that the particles could not pass through isosurfaces of the FTLE field. Moreover, Giudici et al. (Reference Giudici, Suara, Soomere and Brown2021) used the concept of LCS to identify successfully areas in the Gulf of Finland where floating items tended to accumulate in large quantities.

3.5.1. Particle tracking during the startup transient state

In order to study particle entrainment, two groups of particles were seeded in the flow field. The first group was initialized in the incoming flow upstream of the array in a region of size ${\sim }0.15D \times 1.1D$ (upstream-seeded – see figure 28). The second group was initialized in the areas surrounding the cylinders and in the near-wake region with dimensions ${\sim }1.6D \times 1.1D$ (vicinity-seeded). The cross-stream width of these regions was varied according to the width of the arrays, which can vary slightly between the selected individuals. The particles’ motions can be observed in the animations provided in supplementary movies 2 and 3, and corresponding snapshots are provided in figure 28 at non-dimensional times $t^{*}=0$ and $t^{*}=40$ (where $t^{*} = tU_{\infty }/d$). The plots in figure 20 show the percentages of particles found in the vicinity of the cylinders (i.e. in a prescribed observation window described by the region where the green-coloured particles are initialized in figure 28), separately for the upstream-seeded particles in figure 20(a) and the vicinity-seeded particles in figure 20(b). The initial results are identical across all four cases until $t^*=2$, since distinct flow features have not yet had time to develop. In figure 20(a), the sudden dip after $t^{*} = 2$ corresponds to the upstream-seeded particles exiting and re-entering the spatial observation window. Individual D ends up losing a significant percentage of the particles after this stage, which is related to most of the flow being directed around the array instead of through the interior. Individuals B and C retain the highest proportions of particles initially ($10 \leq t^{*} \leq 20$), after which individual A outperforms individual C. Similarly, in figure 20(b), individual B retains almost twice as many vicinity-seeded particles in the spatial observation window as the other three arrays. This may be because in the initial transient state, i.e. before vortex shedding begins, the recirculation region is a prominent factor in determining the extent of particle retention. In figure 14, we observe that individuals B and C generate large recirculation regions in their wakes. Overall, individual B is most effective in retaining both the upstream-seeded and the vicinity-seeded particles in the initial transient state, with individual A being a close second.

Figure 20. Percentage of seeded particles that are present in the arrays’ vicinity, i.e. in a prescribed observation window, over the first $40t^*$ of transient flow in DNS. (a) Plots representing the time evolution of the upstream-seeded particles. (b) Plots representing the time evolution of the vicinity-seeded particles.

3.5.2. Particle tracking during the steady vortex-shedding state

Figure 29 shows particle tracers for the four selected individuals, through three vortex shedding cycles once the wake flow reaches a steady shedding state. The plots in figure 21 show the percentage of the initialized particles that remain in the vicinity of the cylinders, i.e. in the observation window described earlier, separately for the upstream-seeded (figure 21a) and the vicinity-seeded (figure 21b) particles. In figure 21(a), we observe that individual C retains the largest proportion of the upstream-seeded particles until about 1.5 time periods, after which individual A takes over. Similarly, in figure 21(b), individuals B and C have the highest retention until one time period, and by the end of three time periods, individual A is the only arrangement to retain approximately 5 % of the vicinity-seeded particles within the observation window, while the other arrangements retain a negligible amount of particles. We will examine the relationship between particle retention and the flow field further with the help of Lagrangian coherent structures in § 3.5.3.

Figure 21. The percentage of initialized particles found in the prescribed observation window over three shedding time periods in the quasi-steady state, for each of the four arrangements using DNS data. (a) Retention of upstream-seeded particles. (b) Retention of vicinity-seeded particles.

3.5.3. Lagrangian coherent structures

We now examine LCS computed at various times during a complete shedding period for each of the four selected arrangements. The FTLE field was computed for each of the four arrays in the steady shedding state, and the resulting images are shown in figure 22. An integration time window of approximately one-tenth of the shedding time period was used at four different stages in a shedding cycle to capture the variation of the LCS with time. The dark lines visible in the FTLE fields correspond to the LCS, where the threshold for distinguishing ridges as LCS is taken to be approximately 70 % of the maximum (Shadden Reference Shadden2011; Giudici et al. Reference Giudici, Suara, Soomere and Brown2021). The most consistent LCS are found in figure 22(a), where the coherent structures maintain their shape throughout the shedding period. This is especially evident in the interior of the array, where the ridges visible in the FTLE field do not vary with time. In figures 22(b), 22(c) and 22(d), the extent of time variation in the LCS keeps increasing as we go from B to D. Since there is minimal flux across LCS, the scenarios with steady LCS are more capable of trapping particles. This implies that the arrangement in figure 22(a) should be highly effective in retaining particles in the steady shedding state, which agrees well with the observations from § 3.5.2.

Figure 22. Forward FTLE field computed from DNS for individuals A–D at $t = t_{o}$, $t_{o} +0.25T$, $t_{o} + 0.5T$ and $t_{o} + 0.75T$, where $T$ is the shedding time period for each of the arrays. The darkest regions represent the ridges of the Lagrangian coherent structures.

3.6. Manually designed arrangements

Based on the characteristics observed for the four selected Pareto-optimal arrays, we now examine several manually-designed arrays and assess their performance in relation to the designs obtained via optimization. Since drag and enstrophy display a general monotonic relationship with each other (figure 6a), many of the array characteristics that maximize drag conflict with those that minimize enstrophy. This makes it difficult to design arrays that attain high drag and low enstrophy simultaneously, which are criteria that can benefit coastal protection designs as discussed earlier. One of the main observations in the previous subsections has been the existence of internal boundary layers within the porous arrays, which remain attached to neighbouring downstream cylinders. Based on experimental observations by Sumner, Price & Païdoussis (Reference Sumner, Price and Païdoussis2000), this phenomenon occurs when the individual cylinders are separated by a relatively small distance (i.e. $P/d = 1 \unicode{x2013} 1.25$), and when the angle that they make with each other relative to the direction of the freestream velocity is small ($\alpha = 0 \unicode{x2013} 20^\circ$).

The manually designed arrays examined here are shown in figure 23, and the resulting drag and enstrophy values are shown in figure 24 along with a few regular arrangements from the initial optimization generation for comparison. The array in figure 23(a) includes a regularly-spaced grid arrangement of individual cylinders in a square shape. The spacing is kept within the $P/d < 1.25$ limit described by Sumner et al. (Reference Sumner, Price and Païdoussis2000) in both the streamwise and cross-stream directions. Based on our prior observations, we expect that this base configuration will lead to low enstrophy, given the existence of internal channels that allow for the formation of internal boundary layers that remain attached to neighbouring cylinders. Furthermore, the relatively small projected frontal area, as well as considerable internal flux through the array, can be expected to lead to low drag. Both these expectations are confirmed in figure 24, which indicates that the performance of the array shown in figure 23(a) resembles closely that of Pareto-optimal individual A, which is located at the left extreme of the Pareto front.

Figure 23. Vorticity from DNS for several manually designed arrays. The corresponding colour bar is provided in figure 7.

Figure 24. The drag and enstrophy values obtained from DNS for the manual arrangements shown in figure 23 ($\times$). The values for two of the configurations from the initial generation are also shown ($+$). The Pareto front from figure 5 is included for comparison.

The second array, shown in figure 23(b), is a slight modification of the base square array, where the cylinders are spaced farther apart in the cross-stream direction. This has two effects: (1) increased flux through the interior results in decreased drag; and (2) the internal boundary layers of oppositely signed vorticity are farther apart, which reduces the extent of vorticity cancellation via diffusion, resulting in increased wake enstrophy. Both these expectations are confirmed from figure 24, which indicates that this design entails lower drag and higher enstrophy than the base square array. The third design, shown in figure 23(c), was obtained by rotating the base square array by $45^\circ$. The increase in projected frontal area can be expected to lead to increased drag. Additionally, the internal boundary layers no longer remain attached to neighbouring cylinders, and exhibit unsteady behaviour, since $\alpha =45^\circ$ now exceeds the range prescribed by Sumner et al. (Reference Sumner, Price and Païdoussis2000) ($\alpha < 20^\circ$). Thus we expect a notable increase in both drag and enstrophy compared to the base square array, which is confirmed from figure 24.

While the first design is relatively simple and attains the expected low-drag low-enstrophy performance, it does not meet the requirement of high-drag low-enstrophy that a coastal protection design might benefit from. The two modifications to the base square array also do not exhibit the required characteristics, since either the drag remains too low or the enstrophy increases noticeably. The fourth design, shown in figure 23(d), was based on a symmetric rearrangement of Pareto-optimal individual D. We observe that while it leads to the expected increase in drag on account of the larger projected frontal area, it also leads to a considerable increase in enstrophy compared to the base square array.

To counter this rise in enstrophy, the design shown in figure 23(e) was developed with the aim of achieving both high drag and low wake enstrophy simultaneously. The projected frontal area for this configuration is similar to that of the D-shaped array described previously, which is beneficial with regard to high drag. However, this configuration no longer allows flow to pass through the array, which removes the possibility of enstrophy reduction via interactions among internal boundary layers. Thus an alternative means of reducing enstrophy was devised in the form of slightly curved boundaries along the top and bottom edges of the array, and the large posterior cavity visible in figure 23(e). The curved boundaries allow the outermost boundary layers to remain attached and give rise to strong primary vortices shed from the top and bottom of the array, which get redirected towards the posterior cavity. As these primary vortices interact with the posterior cavity, they generate strong secondary vortices, as can be observed in figure 23(e). This interaction between the primary and secondary vortices contributes to vorticity diffusion in the same manner as achieved by internal boundary layers of opposing vorticity, and results in lower wake enstrophy. This particular design was observed to perform well with regard to both metrics, with drag being comparable to that of the rotated square array from figure 23(c), and enstrophy being only slightly higher than that for the base square array, indicating that it may be a suitable candidate for coastal structure design. We note that the design lies on the Pareto front in figure 24, indicating that it is difficult to exceed the performance of the designs obtained via optimization.

In addition to the 9-cylinder arrays, an additional 18-cylinder design was devised, as shown in figure 23( f) (with $C_d = F_{Drag}/(0.5\rho U_\infty ^2 \times 18d)$), to determine whether the observations presented here may extend to larger arrays. The design of this porous array aims specifically to reduce wake enstrophy by promoting interactions among boundary layers of opposite vorticity that form in the interior. We note from figure 24 that the enstrophy for this design is comparable to, and lower than, that of several 9-cylinder arrays, including the two regular arrangements from the initial generation.