Porous cylinder arrays for optimal wake and drag characteristics

Abstract

The root systems of mangroves, a tree species found in intertidal tropical and subtropical coastal zones, provide a natural barrier that dissipates wave energy effectively and reduces sediment erosion. Here, we use a combination of experiments and numerical simulations to examine the wake and drag characteristics of porous arrays of cylinders, which serve as simplified models of mangrove root networks. Optimal arrangements of the arrays are obtained by coupling Navier–Stokes simulations with a multi-objective optimization algorithm, which seeks configurations that minimize wake enstrophy and maximize drag on the porous structure. These optimal configurations are investigated using particle image velocimetry, and the internal and external flows around the porous arrays are analysed using a combination of proper orthogonal decomposition and Lagrangian particle tracking. Large variations in drag and enstrophy are observed by varying the relative positions of the cylinders, which indicates that the geometrical arrangement of porous arrays plays a prominent role in determining wake and drag characteristics. A sensitivity analysis suggests that enstrophy is more sensitive than drag to specific cylinder placement, and depends on distinctive flow patterns that develop in the interior due to interactions among neighbouring cylinders. Arrays with higher drag involve a combination of larger projected frontal area and minimal flux through the interior, leading to increased wake enstrophy, which is unfavourable for particle deposition and erosion. Based on the analysis of characteristics associated with the optimal arrays, several manually designed arrays are tested, and they display the expected behaviour with regard to drag and enstrophy.

1. Introduction

Coastal vegetation such as mangrove forests is one of the most effective natural forms of coastal protection. Mangroves usually grow in tropical and subtropical regions, and provide many ecosystem services, such as carbon dioxide sequestration (Sanderman et al. 2018), nurturing juvenile fish (Mumby 2006), wave energy attenuation, and reducing coastal erosion (Furukawa, Wolanski & Mueller 1997; Barbier et al. 2011; Mitsch & Gosselink 2015; Narayan et al. 2016). The root systems of these mangrove forests have inspired novel designs for artificial sea walls (Takagi 2017) with the aim of mitigating storm damage along coastlines. Compared to impervious breakwater structures, porous bio-inspired designs require less material, and are less disruptive to natural habitats as they allow the natural ebb and flow of sediments. A recent experimental study demonstrated that mangrove-inspired roots models could also help to reduce sediment erosion (Kazemi, Castillo & Curet 2021). While significant effort in the literature has been directed towards studying the influence of porous vegetation on the flow (Leonard & Reed 2002; Rominger & Nepf 2011; Kazemi, Van De Riet & Curet 2018), the design of optimal porous breakwater structures requires further examination. Moreover, a detailed investigation of flow patterns that develop within such porous structures can help us to understand their influence on sediment deposition as well as the intensity of vorticity and enstrophy generated in the wake.

A number of studies have investigated various aspects of flow patterns generated by isolated finite-sized vegetation clusters. Anderson & Charters (1982) found that while flexible plants reduce the turbulence intensity in the oncoming flow, they also introduce significant unsteadiness due to vortex shedding by individual roots in the system. Zong & Nepf (2012) conducted experiments using cylinders to show that the flow velocity decreases immediately behind the vegetation cluster, due to increased drag exerted by the porous obstacle. Ricardo, Sanches & Ferreira (2016) showed that the incoming flow enters the porous patch and passes through to the wake, which can delay the onset of the von Kármán vortex street up to a certain distance behind the patch (e.g. $L1$ – see figure 1). This length $L1$ is a function of the wake velocity, which depends on the freestream velocity and the extent of flow blockage (Chen et al. 2012). Over this distance, the mean velocity is lower in magnitude compared to the freestream velocity, which leads to enhanced sediment deposition. Beyond this region, the von Kármán vortex street leads to a reduction in sediment deposition (Leonard & Luther 1995).

Figure 1. Schematic showing the horizontal cross-section through a simplified model of a mangrove root system, with a uniform inflow imposed from left to right ($U_{\infty }$). The green cylinders represent auxiliary mangrove roots located around a main trunk indicated by the brown cylinder in the centre. Here, $D$ represents the maximum nominal diameter of the array, whereas $d$ represents the diameter of each individual cylinder in the array. Also, $r_{i}$ and $\theta _{i}$ indicate the radial distance and azimuthal angle with respect to the central cylinder, and their values determine the overall layout of the array. Length $L1$ represents the steady wake region between the end of the patch and the onset of vortex shedding.

For sea walls, the drag force and extent of sediment deposition in the wake are two important functional characteristics. These two aspects have been studied primarily with regard to their dependence on root density, or the porosity of uniformly distributed arrays of cylinders, as in several of the studies mentioned previously. Root porosity is usually defined as the fraction of volume occupied by water in a given water–root space. Experiments conducted using model forests have demonstrated the existence of a high-drag region associated with the roots, resulting in reduced velocity within that zone (Maza et al. 2017). Furthermore, Chang & Constantinescu (2015) used three-dimensional simulations to demonstrate the dependence of drag and steady wake length on porosity. Shan, Liu & Nepf (2019) found that random physical arrangements of root models resulted in higher drag compared to uniform grid arrangements, and that the random arrangements resulted in higher variations in local velocity and drag on the individual cylinders. Norris et al. (2019) examined the role of turbulence within actual mangrove canopies, and found that high root densities enhance turbulent mixing, which in turn is detrimental to sediment deposition. Several of these studies remark that for a group of cylinders, the drag force and the length of the sediment deposition region are two competing characteristics that depend on porosity and geometrical arrangement within the porous patch. Ideal designs should result in high drag and long sediment deposition regions. However, the drag force is known to decrease with increasing porosity (Kazemi et al. 2021), whereas the sediment deposition length increases at lower porosity (Kazemi, Van de Riet & Curet 2017). The trade-off between drag force and length of sediment deposition underscores the need for an effective method for optimizing the arrangement of the individual cylinders that constitute the porous patch.

A few studies have investigated the compromise between drag force and sediment deposition length. For instance, Ricardo et al. (2016) examined the wakes for cylinders arranged in random arrays using particle imaging velocimetry (PIV), and found that vorticity cancellation caused by neighbouring cylinders within the array was an important factor in suppressing the wake. However, there are few studies that aim to leverage this finding to optimize the drag or turbulence intensity generated by a porous array of cylinders. Most studies that examine the wake characteristics’ dependence on the relative locations of cylinders consist of cylinders arranged in predetermined geometrical configurations. Moreover, most such arrays are either spread uniformly across the area under consideration in a staggered or inline arrangement, or arranged completely randomly (Chen et al. 2012; Gijón Mancheño et al. 2021). Maza, Lara & Losada (2015) determined using three-dimensional simulations of cylinder array distributions that wave-induced forces and wave attenuation were affected notably by the geometrical arrangement, and that using uniform models as a basis for studying natural vegetation could lead to an incorrect estimation of certain characteristics. This can be attributed to differences in local velocity that are dependent on an individual cylinder's placement within the array, resulting in large variations in drag force experienced by individual cylinders (Tinoco, San Juan & Mullarney 2020).

Other studies of porous arrays have adopted simplified mathematical models to focus on rigid cylinder arrays with various arrangements in uniform flow (Nepf 1999; Nicolle & Eames 2011; van Rooijen et al. 2018). However, these models can involve certain limitations and simplifications with regard to important physical phenomena such as viscous effects and vortex shedding, and often do not provide detailed flow field information. Most of the experimental studies mentioned previously also entail certain limitations as described by Tinoco et al. (2020), such as their time-consuming and labour-intensive nature, which in turn necessitates the use of simplifying parameters such as uniform arrangements. This makes it difficult to examine a large and varied assortment of cylinder arrangements when studying optimal characteristics experimentally. Numerical simulations present a useful complement to experiments, at least at low to moderate Reynolds numbers, since they enable the evaluation of a large number of samples, which is a crucial requirement for conducting optimization studies. In the present work, we couple two-dimensional direct numerical simulations (DNS) of flow around porous cylinder arrays with a multi-objective optimization algorithm, in order to discover geometrical arrangements that optimize drag force and wake enstrophy simultaneously. We note that while the present study is limited to two-dimensional flows owing to the high computational cost involved with optimization, three-dimensional flow structure is important when considering sediment transport and deposition, which is modulated primarily by the near-bed flow. This near-bed flow depends heavily on the formation and interaction of horseshoe vortices near the base of individual cylinders.

The main aim of the present work is to discover optimal array configurations that maximize drag and minimize wake enstrophy, and to examine the underlying mechanisms that influence these characteristics. Higher drag leads to a higher velocity deficit in the wake (Maza et al. 2017; Gijón Mancheño et al. 2021), which is preferable for coastal defence structures inspired from mangrove vegetation since it is more effective in mitigating the impact of the incoming flow. However, while high drag leads to increased energy dissipation in the wake, the turbulent kinetic energy produced in this process generates stronger vortices (i.e. higher wake enstrophy) and hinders sediment deposition, and promotes resuspension (Chen et al. 2012; Tinoco & Coco 2018). Thus a multi-objective optimization approach is adopted in the present work in an attempt to attain simultaneously the conflicting metrics of high drag and low enstrophy. The optimal cylinder arrangements obtained in this manner are also examined experimentally using PIV to support the numerical findings.

Details of the optimization algorithm, the simulation methodology and the PIV set-up are provided in § 2. The optimal cylinder arrangements discovered by the optimization procedure are discussed in § 3, and the flow field both within and outside the porous arrays is examined in detail. A short discussion is presented in § 4, followed by concluding remarks in § 5.

2. Methods

In the present work, arrays of circular cylinders are used as simplified models for mangrove roots. This approach is consistent with several of the studies discussed in § 1, which have shown that three-dimensional effects are of second order with regard to the wake structure (Zong & Nepf 2012; Verschoren et al. 2016). The individual mangrove roots are represented as two-dimensional (2-D) rigid cylinders with diameter $d$, which are grouped together in various arrangements to form a porous patch with maximum nominal diameter $D$, as shown in figure 1. The Reynolds number for all the configurations examined in this work is $Re_d = U_{\infty }d/\nu = 500$, based on the diameter of an individual cylinder $d$ and the uniform inflow velocity $U_{\infty }$. The porosity of the 9-cylinder arrays, defined as the solid volume fraction within a given reference volume (Nicolle & Eames 2011; Chang & Constantinescu 2015), is $\phi = 9(d/D)^{2} = 0.316$, and is constant for all the 2-D arrays discussed here since they contain the same number of identical cylinders per unit reference area.

2.1. Numerical simulations

The simulations used in this work solve the 2-D incompressible Navier–Stokes equations. Brinkman penalization is used to enforce the no-slip boundary condition at the fluid–solid interfaces with the help of a penalty forcing term:

(2.1)\begin{equation} \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}= \frac{-\boldsymbol{\nabla} P}{\rho} + \nu\,\nabla^2 \boldsymbol{u} + \lambda \chi (\boldsymbol{u}_{\boldsymbol{s}} - \boldsymbol{u}). \end{equation}

Here, $\lambda =1/{\rm d}t$ is the penalization parameter (with ${\rm d}t$ being the time step size), and $\chi$ is the characteristic function that represents the solid in discretized form on a Cartesian grid. Grid cells with $\chi = 0$ are occupied entirely by the fluid, and those with $\chi = 1$ are occupied entirely by the solid; $\chi$ transitions smoothly from $1$ to $0$ within 2 grid points at the interface, using a discrete representation of the Heaviside function. The pointwise velocity $\boldsymbol {u_s}$ of the discretized solid accounts for the motion of the solid object, if any. In the present work, $\boldsymbol {u_s} = 0$ for the rigid stationary cylinders.

To solve the Navier–Stokes equations, we use an open-source solver based on the remeshed vortex method (Rossinelli et al. 2015; Verma et al. 2017). This method uses the vorticity form of the momentum equation, which is obtained by taking the curl of (2.1):

(2.2)\begin{equation} \frac{\partial \omega}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\omega= \nu\,\nabla^2 \omega + \lambda\, \boldsymbol{\nabla} \times \left(\chi (\boldsymbol{u}_{\boldsymbol{s}} - \boldsymbol{u})\right). \end{equation}

In obtaining this equation, we have used $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u} = 0$, and the fact that the vortex-stretching term $\boldsymbol {\omega }\boldsymbol {\cdot }\boldsymbol {\nabla } \boldsymbol {u}$ is absent in two dimensions. Godunov splitting is used to perform time splitting of advection, diffusion and the penalty forcing term. The sequential steps used for solving (2.2) can be found in Verma et al. (2017). We note that the simulations use the free-space boundary condition, which allows the use of Green's function for solving the Poisson's equation in the vortex method. A validation of the numerical simulations for flow around 2-cylinder configurations is shown in Appendix A.

2.2. Optimization algorithm

To determine the best possible configuration of cylinder arrangements, we employ a multi-objective genetic optimization algorithm called the non-dominated sorting genetic algorithm II (NSGA II) (Deb et al. 2002). The optimizer attempts to simultaneously maximize the drag acting on the structure and minimize the intensity of vortices in the wake, which is represented by the enstrophy in this case. In general, genetic algorithms operate on a collection of ‘individuals’ and iteratively improve a population's characteristics by assigning each individual a certain ‘fitness’ value or rank. The fitness value represents the eligibility of a particular individual to pass on its characteristics to future generations. In the present work, each ‘individual’ corresponds to one distinct arrangement of 8 cylinders around a central fixed cylinder, and is characterized by 8 pairs of $r_{i}$ and $\theta _{i}$ values as shown in figure 2. A large number of distinct individuals are created and evaluated throughout the optimization process. Fitness values are assigned to each individual based on the average enstrophy observed in the wake, and the net drag experienced by the cylinder array. These quantities are determined via 2-D DNS, where the particular individual (i.e. array arrangement) being evaluated is placed in a uniform inflow, as depicted in figure 1. The average wake enstrophy is computed within a predefined window $\varOmega \in [3.21D \times 2.58D]$ downstream of the arrays, which is sufficiently large to accommodate the spanwise extent of the vortices generated in the wake. Each individual is evaluated for a total of 120 units of non-dimensionalized time ($t^{*} = tU_{\infty }/d$), during which the time-averaged enstrophy and net drag are computed. The enstrophy $E$ is computed as the double integral of the square of the vorticity magnitude $\omega$ within the predefined region $\varOmega$:

(2.3)\begin{equation} E = \iint_\varOmega| \boldsymbol{\omega}| ^2\,{\rm d}A. \end{equation}

Figure 2. Randomly generated individuals from the first generation used to initialize the optimization algorithm. Only 8 out of the 128 individuals generated for the first generation are shown here. Each of these randomized arrangements consists of one central cylinder, and 8 equally spaced cylinders that are also equidistant from the centre. Subsequent generations created by the optimization procedure are not constrained to using equal radii and azimuthal angles for the 8 surrounding cylinders.

The drag is computed using the streamwise component of the penalty force determined from the 2-D simulations:

(2.4a)$$\begin{gather} \boldsymbol{F} = \iint \rho \lambda \chi \boldsymbol{u} \,{\rm d}A, \end{gather}$$

(2.4b)$$\begin{gather}F_{Drag} = \dfrac{\boldsymbol{F} \boldsymbol{\cdot} \boldsymbol{u}}{\lVert \boldsymbol{u} \rVert}. \end{gather}$$

Once the two fitness values, i.e. drag and enstrophy, are determined for each individual using the numerical simulations, the individuals can be ordered using the concept of Pareto dominance. A member $\alpha _1$ of a population dominates another member $\alpha _2$ if: (1) $\alpha _1$ is no worse than $\alpha _2$ for all the objectives under consideration; and (2) $\alpha _1$ is strictly better than $\alpha _2$ in at least one objective. Each individual of the population is assigned a rank according to its level of non-dominance. Individuals with the same non-dominance ranks form a front, and are ranked amongst themselves using the crowding distance parameter, as described by Deb et al. (2002). The crowding distance parameter is an estimate of the average side-length of the cuboid formed by an individual's two closest neighbours. Individuals with a higher crowding distance are preferred because they increase diversity in the solution. Once the individuals are ranked according to their dominance and crowding distance parameter, the best ones are selected to act as parents for creating the next generation of individuals using tournament selection, crossover and mutation (Deb et al. 2002). In the present work, the probability of crossover was set to 0.9, with crossover distribution index 5, whereas the probability of mutation was set to 0.5 with mutation distribution index 10. Furthermore, an additional constraint was applied so that members of the new generation that resulted in overlapping cylinders were discarded, and a new replacement individual was sampled by the optimizer. The initial condition for the optimization algorithm was a set of 128 configurations that consisted of a central cylinder and 8 surrounding cylinders in a regular circular arrangement at a randomized but equal azimuthal angle and radial distance from the centre (figure 2). The optimization algorithm was run for a total of 33 generations, until notable changes in the Pareto front were no longer observed. A flowchart depicting the optimization process is shown in figure 3.

Figure 3. Flowchart depicting the multi-objective optimization of the porous cylinder arrays.

2.3. Particle image velocimetry

The behaviour of the optimal cylinder arrangements discovered via numerical optimization was also examined using a series of PIV experiments carried out in a closed-loop water flume with dimensions $2\ {\rm m}\times 0.25\ {\rm m}\times 0.25\ {\rm m}$ (figure 4). The optimized cylinder arrays were mounted at distance $0.7\ {\rm m}$ ($27.96D$) from the inlet. The blockage ratio of the test section, defined as the patch diameter divided by the channel width, was approximately $10\,\%$. The inflow velocity was set to $0.1\ {\rm m}\ {\rm s}^{-1}$, which corresponds to Reynolds number $Re_d = U_\infty d/\nu = 500$ based on individual cylinder diameter $d$. The diameter of each individual cylinder in the experiments was $d=5\ {\rm mm}$, and the maximum allowable nominal patch diameter was set to $D = 26.7\ {\rm mm}$. The PIV images were recorded after reaching a quasi-steady state for the periodic vortex shedding, i.e. approximately three minutes after starting the inflow. The flow was seeded with hollow glass spherical particles of diameter $10\ \mathrm {\mu } {\rm m}$, which were illuminated with a continuous 5 W laser with $532\ {\rm nm}$ wavelength. The particles were recorded using a high-speed camera (Photron Fastcam Mini UX50) with resolution $1280\times 1024$ pixels, at 125 frames per second, using shutter speed $1/250\ {\rm s}$. A total of 1000 images were recorded for each experimental run, which is equivalent to approximately 8 vortex-shedding cycles. The background was subtracted and a fast Fourier transform correlator was used for the PIV analysis. The velocity vector field was computed from particle displacement within a fixed elapsed time (8 ms) from two sequential images. Measurements of the freestream velocity indicated that the turbulence intensity in the freestream was less than $0.5\,\%$.

Figure 4. (a) Particle image velocimetry set-up used for analysing the wakes of the cylinder arrays. (b) Schematic of the basic components used in the set-up. (c) An image capture showing the seeded particles used for velocimetry.

3. Results

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. (2012), 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. 2017; Gijón Mancheño et al. 2021), which can promote gravitational settling. However, this process is also accompanied by increased turbulent kinetic energy production (Chen et al. 2012; Tinoco & Coco 2018), 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. 2012; Norris et al. 2019) 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. 2019). 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. 2016), 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 (2004), 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. 2004; Ricardo et al. 2016) 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. (2019), 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 1998; Kothyari, Hashimoto & Hayashi 2009). 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 2000; Shadden, Lekien & Marsden 2005; Lekien, Shadden & Marsden 2007). Ku & Hwang (2018) 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. (2021) 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 2011; Giudici et al. 2021). 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 (2000), 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. (2000) 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. (2000) ($\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.

4. Discussion

The present study extends earlier work by Kazemi et al. (2017, 2018), where the porosity and spacing between individual cylinders were varied for an array of nine cylinders. In other experiments involving porous cylinder arrays, Chen et al. (2012) considered cylinders of various shapes, diameters and patch densities to study the influence of flow blockage on the velocity, turbulent kinetic energy and sediment deposition in the wake. They found that flow blockage is an effective predictor of the extent of the diminished velocity region and streamwise exit velocity directly behind the cluster. However, we find that this metric is not the sole predictor of the diminished velocity region, since individuals B and D are similar to each other with regard to flow blockage, but they give rise to diminished velocity regions of different lengths, as is evident in figures 9, 10 and 11. Another experimental study that considered the drag acting on individual cylinders in a group was conducted by Shan et al. (2019), where they placed randomly as well as uniformly distributed tree models in a uniform inflow. They demonstrated that in uniform arrangements, the trees on the front line experienced the highest drag, while the rest were increasingly sheltered, i.e. they did not contribute significantly to the total drag. However, in the case of random distributions, they found significant variation in the locations of cylinders experiencing the largest drag. A similar conclusion may be drawn from the sensitivity analysis presented in figure 19, where the cylinders that most influence the drag are not always the ones in the front of the array. A prior study by Cheer & Koehl (1987) also draws a correlation between drag and spacing in cylinder arrays, albeit at very low Reynolds numbers. They found that bristled appendages of small organisms, which were treated as porous cylinder arrays, functioned to provide locomotion in certain conditions, and were used to catch particulate food in other conditions. They concluded that it was a combination of ‘leakiness’ and frontal area that determined the hydrodynamic forces acting on the appendages (i.e. the simplified cylinder arrays), which is similar to observations presented here. Overall, the results presented here demonstrate that even minor changes to internal flow patterns within a porous array can lead to substantial changes in the wake enstrophy and drag. Thus studying the influence of the relative geometrical arrangement of cylinders can lead to an improved understanding, and consequently to better models of flow in porous arrays.

We note that there is no strict reason for specifically selecting nine cylinders for the arrays in the present study, except as a balance between sufficient generality and too many degrees of freedom. As discussed earlier, mangrove trees usually consist of a central trunk surrounded by a number of auxiliary roots, and the arrays used here are meant to be simplified model representations of this structure. Choosing more cylinders for the arrays would lead to an increased degree of freedom in the array configurations, in turn allowing more diverse configurations to arise. However, the increased dimensionality would pose difficulties with regard to optimization and computational cost. Thus nine cylinders were selected as a reasonable compromise. Furthermore, we note that the stopping criterion for the NSGA II optimization algorithm is generally user-defined. In the present work, there were no significant changes to the Pareto front for several generations leading up to generation 33, which led to the decision to stop further iterations. The convergence of the optimization is further supported by the fact that none of the manually designed arrangements presented in § 3.6 performs significantly better than the Pareto-optimal individuals.

In the sensitivity analysis presented in § 3.4, we observed that enstrophy shows a stronger dependence on $\theta _i$, whereas drag shows a stronger dependence on $r_i$. This can be explained by the fact that changing $\theta _i$ does not impact the projected frontal area of the array significantly, whereas perturbing $r_i$ can result in notable changes to the area, depending on a cylinder's location. This makes drag more sensitive to radial perturbations. Enstrophy, on the other hand, was observed to be sensitive to disruption in boundary layer attachment among neighbouring cylinders. Such disruptions occurred more readily upon perturbing the cylinders in the azimuthal direction rather than in the radial direction, which explains the heightened sensitivity of enstrophy to $\theta _i$.

One aspect that differs between the present study and isolated mangrove trees is the tendency of natural peripheral roots to be distributed in concentric arrangements around the central trunk. This is likely related to the function of peripheral roots as mechanical support for the main trunk against wind and waves, which is not considered as an objective in the present optimization. We also note that the results and observations presented in the current work were obtained at a relatively modest Reynolds number, and may need to be revisited when considering significantly higher Reynolds numbers. Nonetheless, the optimization procedure and analyses presented here, which leverage a complementary combination of numerical and experimental approaches, can be extended to other related scenarios in an effort to minimize the labour-intensive and time-consuming nature of discovering optimal configurations.

5. Conclusion

In the present work, 2-D DNS of flow around porous cylinder arrays have been coupled with a multi-objective optimization algorithm to study the dependence of drag and wake enstrophy on the relative positioning of individual cylinders. In a group of nine cylinders placed in uniform flow, the radial and azimuthal positions of eight surrounding cylinders with respect to a central cylinder are modified by the optimization algorithm in an attempt to simultaneously minimize enstrophy and maximize drag. These characteristics are desirable for coastal defence structures, since low enstrophy is more conducive to promoting sediment deposition and suppressing resuspension, whereas high drag increases velocity deficit in the wake, mitigating the impact of the incoming flow. Drag and enstrophy show a general monotonic relationship, that makes multi-objective optimization a suitable approach for attaining the desired conflicting metrics simultaneously. The optimization process converges to a set of Pareto-optimal individuals after iterating through several generations. Four such individuals were selected for further analysis: a low-enstrophy low-drag individual (A), a high-enstrophy high-drag individual (D), and two individuals with characteristics in between the two extremes (individuals B and C).

Using particle image velocimetry, the streamwise velocity along streamwise and cross-stream line cuts was obtained for the four selected arrays. Individual A was observed to display the least amount of velocity deficit, whereas individual D displayed the largest velocity deficit in the wake. The differences in drag and enstrophy among the four individuals are related to how the flow behaves within and around the arrays, which was analysed further using data from the numerical simulations. Individual A was observed to have noticeable gaps within the array that created a path for the upstream flow to pass through the interior, which had the dual impact of reducing both drag and wake enstrophy. The lower enstrophy was determined to result from vorticity cancellation among boundary layers of oppositely signed vorticity, which formed due to flow being directed through closely spaced neighbouring cylinders in the array's interior. The high-drag configurations were observed to allow minimal flux through the porous interior, with the majority of the freestream being redirected around the arrays instead of passing through them. This resulted in shedding patterns that closely resembled those of rigid impervious objects. The high-drag arrays were also observed to have larger projected frontal areas.

A sensitivity analysis was conducted to determine the influence of an individual cylinder's position on the resultant drag and enstrophy. None of the perturbed individuals were 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. The perturbations tended to have a greater influence on wake enstrophy than on drag, which indicates that enstrophy is more sensitive than drag to specific cylinder placement within porous arrays. The most sensitive cylinders were found to be those that could significantly alter the internal flow through the array when perturbed. The least sensitive cylinders tended to be isolated from other cylinders in the array.

To better understand the impact on particle sedimentation and erosion, tracer particles and Lagrangian coherent structures were used to examine particle transport through and around the porous arrays. The arrangement with lowest enstrophy (individual A) was found to retain the largest number of particles in the steady shedding state, while the arrangement with a large recirculation region and a delayed start for vortex shedding (individual B) was found to retain a comparable number of particles in the initial transient state. High array drag was observed to reduce particle retention in the quasi-steady shedding state, and these observations were supported by the behaviour of the Lagrangian coherent structures, where the ridges did not vary with time for individual A, but displayed significant time variation for individuals B, C and D. Thus although increased drag leads to lower streamwise wake velocity along the midline, the flow speed experienced by particles in the near-wake region can be relatively high on account of strong vorticity (and hence high enstrophy) induced by the primary shedding vortices, making it more difficult for particles to sediment.

Based on the characteristics observed for the four selected Pareto-optimal arrays, several manually designed array configurations were examined, and their performance was compared to the designs obtained via optimization. The relationship between wake enstrophy and boundary layer interactions among neighbouring cylinders, and its dependence on the relative positioning of cylinders within the arrays, was confirmed using simple designs resembling square array configurations, and even using an 18-cylinder configuration that was shown to generate low wake enstrophy. One of the 9-cylinder designs was developed with the aim of achieving both high drag and low wake enstrophy simultaneously, and relied on a large projected frontal area to maximize drag, and a large posterior cavity where the primary vortices being shed from the array generated strong secondary vortices resulting in reduced wake enstrophy. However, none of the manually designed configurations was able to surpass the Pareto front, indicating that it is difficult to exceed the performance of the designs obtained via optimization.

We note that the selection of the Reynolds number ($Re_d = 500$) and porosity ($\phi = 0.316$) in the current work is a limitation not of the optimization method itself, but rather of using 2-D simulations and for keeping computational cost manageable. This is because the overall flow patterns are influenced not by the individual cylinder diameter $d$, but by the array size, which can include up to 5 cylinders lined up in the cross-stream direction, which corresponds to $Re=2500$. Results from 2-D simulations may not be a valid representation at much higher $Re$ values due to the absence of the vortex-stretching term. The porosity value, on the other hand, is determined by the use of 9 cylinders within the given reference area, a number that was selected as a balance between sufficient generality and too many degrees of freedom in the array configurations.

In conclusion, the observations presented here indicate that the internal geometric arrangement of cylinders in porous arrays plays a critical role in regulating the internal flow patterns, which in turn influence the overall wake and drag characteristics. These aspects may not always be accounted for when using randomized array configurations or using pre-specified uniform arrangements, and relative cylinder positions should be factored in alongside other common parameters such as cylinder spacing, number and diameter when studying such porous arrays. The findings from this study can prove to be useful for designing improved coastal protection infrastructure that can dissipate wave energy effectively, promote sediment deposition, and reduce erosion.