Turbulent flows over dense filament canopies

TL;DR

This study uses DNS to analyze turbulent flows over dense filament canopies, revealing how element spacing and height influence flow structures, turbulence characteristics, and the formation of Kelvin-Helmholtz-like rollers.

Contribution

It provides new insights into the flow dynamics over dense canopies, especially regarding the role of spacing and height in turbulence and roller formation, supported by DNS data.

Findings

Kelvin--Helmholtz-like rollers dominate within dense canopies.

Spacing influences dispersive velocity fluctuations and roller wavelength.

Shallow and very small spacings inhibit roller formation.

Abstract

Turbulent flows over dense canopies of rigid filaments of small size are investigated for different element heights and spacings using DNS. The flow can be decomposed into the element-coherent, dispersive flow, the Kelvin--Helmholtz-like rollers typically reported over dense canopies, and the background, incoherent turbulence. The canopies studied have spacings $s^+ = 3$--$50$, which essentially preclude the background turbulence from penetrating within. The dispersive velocity fluctuations are also mainly determined by the spacing, and are small deep within the canopy, where the footprint of the Kelvin--Helmholtz-like rollers dominates. Their typical streamwise wavelength is determined by the mixing length, which is essentially the sum of its height above and below the canopy tips. For the present dense canopies, the former remains roughly the same in wall-units, and the latter, which scales with the drag length, depends linearly on the spacing. This is the result of the drag being essentially viscous and governed by the planar layout of the canopy. In shallow canopies, the proximity of the canopy floor inhibits the formation of Kelvin--Helmholtz-like rollers, with essentially no signature for height-to-spacing ratios $h/s \approx 1$, and no further inhibition beyond $h/s \approx 6$. Very small spacings also inhibit the rollers, due to their obstruction by the canopy elements. The obstruction decreases with increasing spacing and the signature of the instability intensifies, even if for canopies sparser than those studied here the instability eventually breaks down. Simple models based on linear stability can capture some of the above effects.

Tables

Table 1. Table 1: Simulation parameters. N x subscript 𝑁 𝑥 N_{x} and N z subscript 𝑁 𝑧 N_{z} are the number of rows of canopy elements in the streamwise and spanwise directions, respectively. The number of points used to resolve each period of the canopy in the streamwise and spanwise directions are n x subscript 𝑛 𝑥 n_{x} and n z subscript 𝑛 𝑧 n_{z} , respectively. u τ subscript 𝑢 𝜏 u_{\tau} is the friction velocity based on the shear at the canopy tips scaled with the channel bulk velocity. \Rey τ subscript \Rey 𝜏 \Rey_{\tau} is the friction Reynolds number based on u τ subscript 𝑢 𝜏 u_{\tau} and δ 𝛿 \delta . The canopy frontal area density, height, spacing and width are λ f subscript 𝜆 𝑓 \lambda_{f} , h ℎ h , s 𝑠 s and w 𝑤 w , respectively.

	Case	$N_x\times N_z$	$n_x\times n_z$	$u_{\tau }$	${\text{Rey}}_{\tau }$	${\lambda }_f$	$h^{+}$	$s^{+}$	$w^{+}$
Smooth	SC	–	–	0.064	186.3	–	–	–	–
	S10	108$\times$54	12$\times$12	0.071	176.7	4.6	95.5	10.3	5.2
Fixed	S16(H96)	72$\times$36	24$\times$24	0.088	187.7	3.1	101.8	16.4	8.2
height	S24	48$\times$24	24$\times$24	0.102	187.4	2.0	101.2	24.5	12.3
($h^{+}\approx 100$)	S32	36$\times$18	24$\times$24	0.112	186.4	1.5	100.7	32.5	16.3
	S48	24$\times$12	24$\times$24	0.124	180.1	1.0	97.0	47.2	23.6
	H16	72$\times$36	24$\times$24	0.071	184.7	0.5	17.2	16.1	8.1
Fixed	H32	72$\times$36	24$\times$24	0.080	188.8	1.0	34.6	16.5	8.3
spacing	H64	72$\times$36	24$\times$24	0.086	186.2	2.0	68.1	16.2	8.1
($s^{+}\approx 16$)	H96(S16)	72$\times$36	24$\times$24	0.088	185.7	3.1	101.8	16.4	8.2
	H128	72$\times$36	24$\times$24	0.086	184.8	4.1	133.2	16.4	8.2
	G10	432$\times$216	9$\times$9	0.064	175.9		10.1	2.6	0.6
Self-similar	G20	216$\times$108	9$\times$9	0.072	190.7		22.2	5.6	1.2
geometry	G40	108$\times$54	18$\times$18	0.106	188.2	0.85	43.4	11.0	2.5
($h/s\approx 4$)	G60	72$\times$36	18$\times$18	0.127	183.3		64.7	16.0	3.6
	G100	48$\times$24	18$\times$18	0.147	185.7		97.9	24.3	5.4
Different ${\text{Rey}}_{\tau }$	H32180	72$\times$36	12$\times$12	0.075	184.6	1.0	33.9	16.1	8.1
	H32400	162$\times$81	12$\times$12	0.066	399.9		32.0	15.5	7.8

Table 2. Table 2: Most amplified instability wavelengths observed in the DNSs and predicted by the stability analysis, scaled in friction units. The column labelled ‘DNS’ lists the approximate streamwise wavelength associated with the instability in the wall-normal spectra portrayed in figure  10 . SA C0 , most amplified wavelengths from stability analysis on DNS mean profiles without drag on fluctuations; SAM C0 , on synthesised velocity profiles without drag on fluctuations; and SA, on DNS mean profiles with drag on fluctuations.

	Case	DNS	SAC0	SAMC0	SA
	S10	140	105	90	320
Fixed	S16(H96)	160	115	105	385
height	S24	200	140	115	420
	S32	230	152	130	465
	S48	250	152	140	165
	H16	130	95	90	140
Fixed	H32	150	105	95	200
spacing	H64	160	115	105	290
	H96(S16)	160	115	105	385
	H128	160	115	105	560
	G10	–	95	70	–
Self-similar	G20	120	95	90	170
geometry	G40	140	140	115	240
	G60	190	170	130	320
	G100	260	220	152	350
Varying ${\text{Rey}}_{\tau }$	H32180	140	105	–	180
	H32400	140	105	–	180

Equations

$$\left[\mathrm{I}-\Delta t\frac{{\beta }_k}{\text{Rey}}\mathrm{L}\right]{\mathbf{u}}_k^n$$

$$\mathrm{DG}({\varphi }_k^n)$$

$${\mathbf{u}}_{k+1}^n$$

$$p_{k+1}^n$$

$$\mathbf{u}=\mathbf{U}+{\mathbf{u}}^{\prime },$$

$${\mathbf{u}}^{\prime }=\tilde{\mathbf{u}}+{\mathbf{u}}^{\prime \prime },$$

$$\frac{\text{p}\mathbf{u}}{\text{p}t}+\mathbf{u}\text{bcdot}\nabla \mathbf{u}$$

$$\nabla \text{bcdot}\mathbf{u}$$

$$\nu C_{x}U=-\frac{\mathrm{d}P}{\mathrm{d}x}.$$

$$\nu C_y⟨v⟩=-\frac{\mathrm{d}P}{\mathrm{d}y},$$

$$\frac{\text{p}u}{\text{p}t}+U\frac{\text{p}u}{\text{p}x}+v{U}^{\prime }$$

$$\frac{\text{p}v}{\text{p}t}+U\frac{\text{p}v}{\text{p}x}$$

$$\frac{\text{p}u}{\text{p}x}+\frac{\text{p}v}{\text{p}y}$$

$$\left(\frac{\text{p}}{\text{p}t}+U\frac{\text{p}}{\text{p}x}+\nu C_y\right){\nabla }^{2}v-\nu {\nabla }^{4}v=U^{\prime \prime }\frac{\text{p}v}{\text{p}x}-\nu (C_x-C_y)\frac{{\text{p}}^{2}v}{\text{p}{y}^2}.$$

$$\begin{gathered} (\alpha U-i\nu C_y)({\mathrm{D}}^2-{\alpha }^2)\tilde{v}-\alpha U^{\prime \prime }\tilde{v}-i\nu (C_x-C_y)D^2\tilde{v} \\ +i\nu ({\mathrm{D}}^4-2{\alpha }^2{\mathrm{D}}^2+{\alpha }^4)\tilde{v}=\omega ({\mathrm{D}}^2-{\alpha }^2)\tilde{v}, \end{gathered}$$

$$\frac{\mathrm{d}}{\mathrm{d}y}\left(\left[\nu +{\nu }_T(y)\right]\frac{\mathrm{d}U}{\mathrm{d}y}\right)-\nu C_x(y)U-\frac{\mathrm{d}P}{\mathrm{d}x}=0,$$

$$\frac{{\mathbf{u}}^{n+1}-{\mathbf{u}}^n}{\Delta t}=\frac{-{\mathbf{u}}^n}{\Delta t}$$

$$\left[\mathrm{I}-\Delta t\frac{{\beta }_k}{\text{Rey}}\mathrm{L}\right]{\mathbf{u}}_k^n=-\Delta t\frac{{\beta }_k}{\text{Rey}}\mathrm{L}{\mathbf{u}}_{k-1}^n.$$

$$\left[\mathrm{I}-\Delta t\frac{{\beta }_k}{\text{Rey}}\mathrm{L}\right]\mathbf{u}\approx \left[\mathrm{I}-\Delta t\frac{{\beta }_k}{\text{Rey}}{\mathrm{L}}_{xz}\right]\left[\mathrm{I}-\Delta t\frac{{\beta }_k}{\text{Rey}}{\mathrm{L}}_y\right]\mathbf{u},$$

$$\left[\mathrm{I}-\Delta t\frac{{\beta }_k}{\text{Rey}}{\mathrm{L}}_y\right]\mathbf{u}={\left[\mathrm{I}-\Delta t\frac{{\beta }_k}{\text{Rey}}{\mathrm{L}}_{xz}\right]}^{-1}\mathrm{RHS}.$$

$$\left[\mathrm{I}-\Delta t\frac{\beta_{k}}{\Rey}\mathrm{L}_{y}\right]\mathbf{u}^{n}_{k}=\left\{\begin{array}[]{ll}-\Delta t\frac{\beta_{k}}{\Rey}\mathrm{L}_{y}\mathbf{u}^{n}_{k-1},&y=y_{tips}\\[2.0pt] 0,&y\neq y_{tips}.\end{array}\right.$$

$$(\alpha U-i{C}_y)$$

Full text

1 Introduction

The present work studies flows over dense canopies of filaments of small size. Canopy flows are mainly studied in the context of natural vegetation canopies, but they also encompass engineering flows where the canopy parameters may be very different from those of natural canopies. Many of the key findings from natural canopy studies have been summarised in the reviews by Finnigan (2000), Belcher et al. (2012) and Nepf (2012). In engineering applications, filament canopies can, for instance, be used to enhance heat transfer (Fazu & Schwerdtfeger, 1989; Bejan & Morega, 1993) and for energy harvesting (McGarry & Knight, 2011; Elahi et al., 2018). Depending on the geometry and spacing of their elements, canopies can be classified as sparse, dense or transitional (Nepf, 2012). Dense canopies typically have small element spacings compared to the lengthscales in the overlying flow, and thus prevent turbulent eddies from penetrating efficiently within the canopy. Sparse canopies, on the other hand, have large element spacings and consequently, turbulent eddies are essentially able to penetrate the full height of the canopy (Poggi et al., 2004; Nepf, 2012; Sharma & García-Mayoral, 2018, 2019). Transitional, or intermediate, canopies would lie between these two regimes. In the present study, we assess how the flow within and above dense canopies is affected by canopy parameters, such as the element height and spacing. The canopies considered have spacings $s^{+}\approx\mathcal{O}(10)$, which should be small enough to limit the penetration of the overlying turbulence within them. The frontal area density $\lambda_{f}$ is also a commonly used measure to categorise canopies (Finnigan, 2000; Poggi et al., 2004; Huang et al., 2009; Nepf, 2012). Canopies with $\lambda\gg 0.1$ are classified as dense, with $\lambda\approx 0.1$ as transitional and with $\lambda\ll 0.1$ as sparse. However, in addition to $\lambda_{f}$, the lengthscales of the overlying turbulence should also be considered when determining the canopy regime. A given canopy geometry with a fixed $\lambda_{f}$ may have element spacings much smaller than any overlying turbulent eddy at a particular Reynolds number, thereby not allowing turbulence to penetrate within the canopy. As the Reynolds number is increased, however, the size of these eddies will eventually become comparable to the element spacing, allowing turbulence to penetrate efficiently within the canopy. To assess this effect, we also study canopies with self-similar geometries, which have a fixed $\lambda_{f}$, but different sizes in friction units.

We also place attention on the effect of canopy parameters on the Kelvin–Helmholtz-like, mixing-layer instability characteristic of dense canopy flows (Raupach et al., 1996; Finnigan, 2000; Nepf, 2012). This instability originates from the inflection point in the mean velocity profile at the canopy-tip plane (Raupach et al., 1996). Kelvin–Helmholtz instabilities manifest as spanwise coherent rollers whose streamwise scale is determined by the shear-layer thickness (Michalke, 1972; Brown & Roshko, 1974). Ghisalberti & Nepf (2004) noted that, while in free-shear flows the shear-layer thickness, and consequently the instability wavelength, continues to grow downstream, in fully developed canopy flows this thickness is constant and is set by the net canopy drag. Therefore, a fixed instability wavelength is generally associated with dense canopy flows. Several studies have shown that some aspects of this instability can be captured using a mean-flow linear stability analysis (Raupach et al., 1996; White & Nepf, 2007; Singh et al., 2016; Zampogna et al., 2016; Luminari et al., 2016). Some studies have also suggested that at the high Reynolds numbers of natural canopy flows these instabilities can be distorted by the ambient turbulence fluctuations and lose their spanwise-coherent nature (Finnigan et al., 2009; Bailey & Stoll, 2016). The importance of this instability decreases as the element spacing is increased, and sparse canopies do not exhibit a notable signature (Poggi et al., 2004; Pietri et al., 2009; Huang et al., 2009; Sharma & García-Mayoral, 2019).

Based on the observations of previous studies, we would expect that the effect of increasing the canopy height for a fixed element spacing on the instability and the surrounding flow would eventually saturate. Ghisalberti (2009) and Nepf (2012) proposed that the effective canopy height perceived by the overlying flow would be a function of the canopy shear-layer thickness as it determined the extent to which the Kelvin–Helmholtz-like instabilities penetrated within the canopy. In flows over arrays of cuboidal posts, Sadique et al. (2017) found that the mean-velocity profiles over them became independent of the element heights at large element aspect ratios. They concluded that the overlying flow only interacted with the region near the element tips, and that the height below this ‘active’ region was dormant, and did not have a significant effect on the overlying flow. For their geometries, Sadique et al. (2017) observed the height of this active region to be related to the element width. A similar observation was also made by MacDonald et al. (2018), who performed direct numerical simulations (DNSs) of flows over spanwise-aligned bars. They found that the gap between the bars was the relevant lengthscale for the overlying flow, and that increasing the height of the bars beyond a certain height-to-gap ratio did not affect the overlying flow, or cause an increase in the drag they produced.

In the present work, we conduct a systematic range of DNSs changing the canopy height and spacing separately in order to study their individual effects on the surrounding turbulence and on the Kelvin–Helmholtz-like instability. We also consider canopy geometries with constant $\lambda_{f}$ for which the height and spacing are changed simultaneously in a fixed proportion. The canopies consist of rigid, prismatic filaments with small element spacings and large height-to-spacing ratios. The element spacings considered, $s^{+}\approx 3$–$50$, are much smaller than those typical of most natural canopy flows and would, for instance, be representative of flows over engineered canopies such as those mentioned previously in this section. We also assess how models based on linear stability analysis capture some of the effects of the canopy parameters on the Kelvin–Helmholtz-like instability.

The paper is organised as follows. The numerical methods used for the simulations and the canopy parameters are discussed in §2. The results from the DNSs, detailing the effect of the canopy parameters on the overlying turbulence and the Kelvin–Helmholtz-like instabilities are discussed in §3. The results from linear stability analysis and a model to capture the instabilities are presented in §4. The conclusions are summarised in §5.

2 Methodology

We conduct direct numerical simulations of symmetric channels with rigid canopy elements on both walls. The streamwise, wall-normal and spanwise coordinates are $x$, $y$ and $z$, with the associated velocities $u$, $v$ and $w$, and $p$ is the kinematic pressure. The wall-normal origin, $y=0$, is defined at the tip plane of the canopies protruding from the bottom wall. The channel height, $2\delta$, is defined as the distance between the tip planes of the canopies on the top and bottom walls. The canopy elements, therefore, extend below $y=0$ and above $y=2\delta$ and have a height $h$. A schematic representation of the channel is portrayed in figure 1. The size of the domain is a standard $2\pi\delta$ in the streamwise direction and $\pi\delta$ in the spanwise direction. We use the channel half-height as the length scale in outer units, which implies that $\delta=1$ in outer scaling. The domain-to-canopy height ratio for most cases considered here is $(\delta+h)/h\approx 3$. We will show in §3 that the height of the roughness sublayer scales with the canopy spacing rather than their height, as in the configurations of Sadique et al. (2017) and MacDonald et al. (2018), and that outer-layer similarity is recovered well below the channel half-height. The channel height to element spacing ratio for most canopies considered is $\delta/s\gtrsim 10$, and for the canopy with the largest element spacing is $\delta/s\approx 4$. The flow is incompressible and the density is always scaled with the fluid density, implying that $\rho=1$. The simulations are run at a constant flow rate, with the viscosity, $\nu$, adjusted to obtain a friction Reynolds number $Re_{\tau}=u_{\tau}\delta/\nu\approx 185$ for most of the cases, where $u_{\tau}$ is the friction velocity calculated at the canopy tips. In order to ascertain the effects of the Reynolds number on the flow, a simulation at $\Rey_{\tau}\approx 405$ was also conducted. The simulation parameters are given in table 1 for reference. Scaling with $u_{\tau}$ and $\nu$ is referred to as in friction or wall units, and scaling with the channel bulk velocity, $U_{b}$, and $\delta$ is referred to as in outer units.

The numerical method used to solve the three-dimensional Navier–Stokes equations is adapted from Fairhall & García-Mayoral (2018). A Fourier spectral discretisation is used in the streamwise and spanwise directions. The wall-normal direction is discretised using a second-order centred difference scheme on a staggered grid. The grid in the wall-normal direction is stretched to give a resolution $\Delta y^{+}_{min}\approx 0.33$ at the canopy-tip plane, stretching to $\Delta y^{+}_{max}\approx 3.3$ at the channel centre. The grid within the canopies preserves the resolution of $\Delta y^{+}_{min}\approx 0.33$ near the canopy-tip plane, and for the tallest canopies considered stretches to $\Delta y^{+}_{max}\approx 4$ at the base of the canopy, where the flow is quiescent. The wall-normal grid distribution for a representative canopy simulation is provided in figure 21 in appendix A for reference. To resolve the element-induced flow while avoiding excessive computational costs, the domain is divided into three blocks in the wall-parallel directions (García-Mayoral & Jiménez, 2011; Fairhall & García-Mayoral, 2018; Abderrahaman-Elena et al., 2019). In the central block, the resolutions in the streamwise and spanwise directions are $\Delta x^{+}\approx 6$ and $\Delta z^{+}\approx 3$, respectively, sufficient to resolve the turbulent eddies. The blocks including the canopy elements and the roughness sublayer have a finer resolution than the central block. In the fine blocks, the limiting resolution is not the one required to resolve the turbulent scales, but that required to resolve the obstacles or the element-induced flow. The resolutions in these blocks are given in table 1. The height of the fine blocks is chosen such that the element-induced flow decays to zero well within the fine-block region, and this is verified a posteriori. The time advancement is carried out using a three-step Runge–Kutta method with a fractional step, pressure correction method to enforce continuity (Le & Moin, 1991)

[TABLE]

where $\mathrm{I}$ is the identity matrix and L, D and G are the Laplacian, divergence and gradient operators respectively. N is the advective term which is dealiased using the $2/3$-rule (Canuto et al., 2012). The Runge-Kutta coefficients, $\alpha_{k}$, $\beta_{k}$, $\gamma_{k}$ and $\zeta_{k}$, for each substep, $k$, are taken from from Le & Moin (1991). The time step is $\Delta t$.

The canopy elements are represented using an immersed-boundary method adapted from García-Mayoral & Jiménez (2011). Further details about the immersed boundary method and validation studies are provided in appendix A. The parameters of the different simulations conducted are summarised in table 1. The simulation denoted by ‘SC’ is of a turbulent channel flow with smooth walls. The canopy-flow simulations are divided into three groups. The canopy elements studied in each group are prismatic, with a square top-view cross section, and their arrangement is illustrated in figure 2. The first group, denoted by the prefix ‘S’, consists of canopies with a fixed height, $h^{+}\approx 96$, and element spacings ranging from $s^{+}\approx 10$ to $48$. The second group, marked by the prefix ‘H’, consists of canopies with a fixed element spacing, $s^{+}\approx 16$, and element heights ranging from $h^{+}\approx 16$ to $128$. The element width-to-spacing ratio for the canopies of S and H is $w/s=1/2$. The final group, denoted by the prefix ‘G’, consists of self-similar elements with a fixed height-to-spacing ratio $h/s\approx 4$, and $w/s=2/9$. The heights for the canopies of G range from $h^{+}\approx 10$ to $100$, with the element spacings varying in proportion to their height. These canopies have a constant $\lambda_{f}=0.85$ and are used to study the effect of changing the canopy size for a fixed geometry. Two additional simulations, H32180 and H32400, are conducted to check the dependence of the results on the friction Reynolds number. The canopy geometries for both these simulations have $s^{+}\approx 16$, $h^{+}\approx 32$ and $w/s=1/2$, with friction Reynolds numbers $\Rey_{\tau}\approx 180$ and $400$. We also conducted several simulations to assess whether the wall-parallel resolutions used in the simulations are sufficient to resolve the element-induced flow. The simulation S24 was run at resolutions of $12$, $24$ and $36$ points per element spacing, and G100 at $9$, $18$ and $27$ points per element spacing. Different resolution sets are used for the geometries of S and G as they have different element width-to-spacing ratios. The rms velocity fluctuations obtained from these simulations are portrayed in appendix A. The simulation results are grid independent at a resolution of $24$ points per element spacing for the geometry of case S24, and $18$ points per spacing for that of case G100. The simulations with $9$ and $12$ points per spacing tend to under-predict the fluctuations within the canopies, with the maximum deviation observed in the wall-normal fluctuations of $20\%$ within the canopy. This discrepancy reduces to $4\%$ outside the canopy. These resolutions are only used for the densest canopy cases, where the fluctuations within the canopy are already very small, and therefore, higher resolution simulations would not change the trends observed. The higher Reynolds number simulation, case H32400, is also simulated using $12$ points per spacing. For this simulation, using a higher resolution would be computationally restrictive. Note that the same resolutions are used for cases H32180 and H32400 to avoid grid related discrepancies in the comparison of their results.

2.1 Reynolds number effect

To analyse the influence of the Reynolds number in our subsequent DNSs, we compare the results of cases H32180 and H32400, which have the same canopy height and spacings in friction units, but different friction Reynolds numbers. The velocity fluctuations and the Reynolds shear stresses within the canopy, and above it up to a height of $y^{+}\approx 10$, of these simulations essentially collapse, as shown in figure 3. This suggests that the flow in the region near the canopy-tip plane scales in friction units, similar to the near-wall region in smooth-wall flows (Moser et al., 1999). Scaling in friction units over conventional rough surfaces has also been noted by Chan et al. (2015). Beyond $y^{+}\gtrsim 10$, we observe that the magnitude of the peaks in the fluctuations and the Reynolds shear stresses are larger for case H32400 compared to case H32180. The increase in magnitude of the near-wall peaks in the velocity fluctuations at friction Reynolds numbers larger than $Re_{\tau}\approx 180$ is consistent with that observed in smooth-wall flows (Moser et al., 1999; Sillero et al., 2013), also included in figure 3 for reference. Further away from the canopy tips, at $y^{+}>50$, the rms velocity fluctuations from the canopy simulations coincide with those from the smooth-wall simulations at their corresponding Reynolds numbers, which indicates the recovery of outer-layer similarity.

In addition to the rms fluctuations being similar for these simulations, the distribution of energy in different scales is also similar. This is illustrated by the pre-multiplied spectral energy densities at $y^{+}\approx 15$, portrayed in figure 4. This height roughly corresponds to the location where the magnitude of the fluctuations peaks in smooth-wall flows (Jiménez & Pinelli, 1999). The results of H32180 and H32400 suggest that the effect of the canopy scales in friction units, and therefore the results presented in the following sections for flows at $\Rey_{\tau}\approx 180$ should also be relevant for higher Reynolds number flows.

3 Effect of canopy parameters on the surrounding turbulence

In this section, we discuss the results obtained from the DNSs, aiming to characterise the effect of the canopy parameters on both the Kelvin–Helmholtz-like instability and the fluctuating flow within and above the canopies.

3.1 Height of the roughness sublayer

Before discussing the effect of the canopy on the overlying flow, we first define the extent of the region affected by the canopy, that is, the height of the roughness sublayer. Over conventional rough surfaces, with heights comparable to or smaller than their spacings, the height of the roughness sublayer is generally observed to be a function of the roughness height (Raupach et al., 1991; Flack et al., 2007; Abderrahaman-Elena et al., 2019). Jiménez (2004) reviewed the effect of various roughness geometries on turbulent flows and noted that, in flows over closely packed spanwise aligned grooves, the flow within each groove would be isolated from the overlying flow due to the ‘sheltering’ effect of the preceding obstacle. The overlying flow in this case would not interact with the full height of the groove. This sheltering effect was also noted by Sadique et al. (2017) for high-aspect-ratio prismatic roughness and by MacDonald et al. (2018) for spanwise aligned grooves with large spacings, and has been used to model cuboidal roughness by Yang et al. (2016). As the element spacings of the canopies studied here are small, this sheltering effect should result in the overlying flow only interacting with the region near the canopy-tip plane. In order to determine the height of this region, we examine the element-coherent flow induced by the canopy elements. The footprint of the element-induced flow can be observed in the instantaneous realisations of the velocity above the canopy-tip plane, portrayed for the canopies of families H and G in figure 5. We isolate the element-induced flow using the standard triple decomposition of Reynolds & Hussain (1972)

[TABLE]

where $\boldsymbol{u}$ is the full velocity, $\boldsymbol{U}$ is the mean velocity obtained by averaging the flow in time and space, and $\boldsymbol{u^{\prime}}$ is the full space- and time-fluctuating signal. The latter is decomposed into the element-induced, dispersive velocity, $\boldsymbol{\widetilde{u}}$, which is obtained by ensemble-averaging the flow in time alone, and the element-incoherent fluctuating velocity $\boldsymbol{u^{\prime\prime}}$, which includes the contributions from the background-turbulence and the Kelvin–Helmholtz-like instability. The rms fluctuations of $\boldsymbol{\widetilde{u}}$, therefore, result from fluctuations in space alone.

We observe that the element-induced fluctuations, for all the canopies studied here, decay exponentially above the canopy-tip plane, and become negligible at a height of one element spacing above regardless of the canopy depth, as shown in figure 6. This suggests that the height of influence of the element-induced flow is determined by the spacing between the elements rather than their height.

Even though the element-induced fluctuations only extend to one element spacing above the canopy-tip plane, their influence on the background-turbulence extends to a height of approximately $2$–$3$ element spacings, as can be observed in figure 7. Abderrahaman-Elena et al. (2019) observed a similar effect over conventional cubical roughness, where the element-induced fluctuations only extended to $y\approx h$, but the effect of the roughness on the overlying flow extended to $y\approx 3h$ above them. At heights of $y/s>2$–$3$ above the canopy-tip plane, the full rms velocity fluctuations collapse with those of smooth-wall turbulence, as shown in figure 7, which is indicative of the recovery of outer-layer similarity. This is verified by a comparison of the pre-multiplied spectral energy densities of the canopy and smooth-wall cases in figure 8, which shows that the energy densities of the canopies of family S collapse with those of the smooth-wall case for $y^{+}\gtrsim 90$. This corresponds to a height of about $2s$ for case S48, the canopy with the largest spacing. Although not shown, the pre-multiplied spectral energy densities of the canopies of families H and G collapse with the smooth-wall spectra for $y/s\gtrsim 3$ as well. Previous canopy studies have proposed that the influence of the canopy elements on the flow above them is set by the wall-normal extent of the Kelvin–Helmholtz-like instability, which is determined by the canopy shear-layer thickness (Ghisalberti, 2009; Ghisalberti & Nepf, 2009; Nepf, 2012). It will be demonstrated in §4 that the shear-layer thickness of the present canopies also depends on the element spacing.

3.2 Effect of element height and spacing

In this section, we discuss the effect of the element height and spacing on the element-induced and on the full velocity fluctuations, both within and above the canopy. We observe that the element-induced fluctuations are largest near the canopy-tip plane and decay below it, as shown in figure 6, because they are obstructed by the canopy elements. This effect is more intense for smaller element spacings, and eventually results in the fluctuations vanishing completely well above the canopy floor for the simulations with the smallest spacings in families S and G. For the canopies of family H, which have a constant element spacing, the change in element height does not have a noticeable effect on the element-induced fluctuations, as shown in figures 6($b,e,h$). For the canopies of families S and G, the intensity of the element-induced velocity fluctuations within the canopy increases with element spacing, when scaled using either the friction velocity or the channel bulk velocity. These results suggest that, for canopy elements with a given width, the magnitude of the element-induced fluctuations is governed mainly by the element spacing.

As discussed in §3.1, in canopies with very small element spacing the height of the roughness sublayer is small, and we would expect such canopies not to disrupt the overlying turbulence significantly, regardless of their depth. In the literature on conventional roughness, small roughness elements that have a negligible effect on the overlying turbulent flow are termed ‘hydraulically-smooth’, as the flow over them remains essentially smooth-wall like (Nikuradse, 1933; Raupach et al., 1991). Roughness elements with a characteristic size of a few wall-units, $h^{+}\lesssim 5$, typically fall into the hydraulically smooth category (Raupach et al., 1991; Jiménez, 2004; Flack et al., 2007). Of the canopies studied here, we observe that the overlying flow for canopy G10, which has an element spacing of $s^{+}\approx 2.6$, is essentially smooth-wall-like above the canopy-tip plane. This is evidenced by the collapse of the rms velocity fluctuations, Reynolds shear stresses, and the mean velocity profile of this case with those of the smooth-wall case, as shown in figures 7($c,f,i$) and 9($c,f$). In addition, the magnitude of the velocity fluctuations below the canopy-tip plane is negligible. This suggests that the overlying turbulent flow essentially perceives the canopy-tip plane as an impermeable wall, and has little or no interaction with the canopy region below this plane.

For canopies with larger element spacings, we begin to observe deviations from smooth-wall-like behaviour in the overlying flow. Above the canopy-tip plane, an increase in the element spacing causes a reduction in the intensity of the streamwise velocity fluctuations and an increase in the intensity of the wall-normal and spanwise ones, as can be observed in figure 7 for the canopies of families S and G. For canopies with large element spacings, such as those of S48 and G100, the peak in $u^{\prime}$ typical of smooth-wall flows, is significantly reduced. These changes in the velocity fluctuations are accompanied by a reduction in the streamwise coherence in the flow with increasing element spacing as can be observed in the instantaneous realisations of the wall-normal velocity for the canopies of family G, portrayed in figure 5. Near-wall turbulence over smooth walls is characterised by streaks and quasi-streamwise vortices, which are predominantly streamwise-coherent (Kline et al., 1967; Jiménez & Pinelli, 1999). The decrease in $u^{\prime}$ and increase in $v^{\prime}$, $w^{\prime}$ above the canopy with increasing element size is also commonly reported over conventional rough surfaces (Ligrani & Moffat, 1986; Orlandi & Leonardi, 2006). Several authors have attributed these changes in the velocity fluctuations and the loss of streamwise coherence to the roughness elements modifying the near-wall cycle and turbulence becoming more ‘isotropic’ (Jiménez, 2004; Flores & Jiménez, 2006; Flack et al., 2007; Abderrahaman-Elena et al., 2019). We also observe an increase in the Reynolds shear stresses above the canopy tip plane with increasing element spacing, shown in figures 9($d,f$), with an associated increase in the drag. The drag increase caused by rough surfaces is generally expressed in terms of the downward shift in the logarithmic region of the mean-velocity profile compared to that for a smooth wall (Hama, 1954). This shift can be observed for the canopies of families S and G in figures 9($a,c$).

Focusing now on the flow within the canopy, increasing the element spacing results in an increase in the magnitude of all the components of the full velocity fluctuations, as shown in figure 7, which is consistent with the observations of Green et al. (1995), Novak et al. (2000), Poggi et al. (2004) and Pietri et al. (2009). The wall-parallel velocity fluctuations, $u^{\prime}$ and $w^{\prime}$, decay rapidly below the canopy-tip plane, and their magnitude reaches a plateau in the core of the canopy, before dropping again near the canopy base to meet the no-slip condition. The abrupt changes in the velocity fluctuations near the element tips are typical of textures with perfectly flat and aligned tips and have also been observed over conventional cuboidal rough surfaces (Leonardi & Castro, 2010; Abderrahaman-Elena et al., 2019) and permeable substrates (Kuwata & Suga, 2017). However, this effect would likely be smeared out over canopies with irregularly aligned tips. The height over which the fluctuations decay within the canopy and the magnitude of the fluctuations in the core of the canopy appear to correlate with the element spacing. Note that this plateau in $u^{\prime}$ and $w^{\prime}$ within the canopy is asymptotic and requires a sufficiently large canopy depth to occur. Thus, this plateau is essentially absent for the canopy of S48, because of its low canopy height-to-spacing ratio, $h/s\approx 2$. The wall-normal fluctuations within the canopy do not exhibit this plateau and decay gradually below the canopy-tip plane to meet the impermeability condition at the canopy base. Let us also note here that the element-induced flow accounts for less than $30\%$ of the magnitude of the streamwise velocity fluctuations and less than $10\%$ of the cross-velocity fluctuations within the canopy, which is consistent with the observations of Poggi & Katul (2008). This implies that the velocity fluctuations deep within the canopy result mainly from the penetration of the overlying, element-incoherent velocity fluctuations. This will be discussed further in §3.3.

Although, as discussed in the preceding paragraphs, the element spacing has a leading-order effect on the fluctuating flow, their height, $h$, also plays a secondary role. In order to assess the effect of height, we consider the canopies of family H, which have fixed element width and spacing, but different canopy heights. As noted previously, the differences between the element-induced fluctuations for the fixed-spacing canopies of family H are negligible. However, we do observe changes in the full rms velocity fluctuations for these cases, implying that the height affects the element-incoherent flow. Above the canopy tips, we observe a decrease in $u^{\prime}$ and an increase in $v^{\prime}$ and $w^{\prime}$ with increasing canopy height, similar to the effect of increasing element spacing, as shown in figures 7($b,e,h$). Within the canopy, $u^{\prime}$ and $w^{\prime}$ for all the cases collapse to the same curves, only departing to meet the no-slip condition at the canopy base. The corresponding magnitude of $v^{\prime}$ within the canopy, however, increases with canopy height up to $h/s\approx 6$, and saturates for $h/s\gtrsim 6$. This saturation is also observed for the effect of the canopy on the flow in general as illustrated in figures 7($b,e,h$) and 9($b,e$), which show that the velocity fluctuations, Reynolds shear stresses and the mean velocity profiles for cases H96 and H128, with $h/s\approx 6$–$8$ are essentially the same. The changes in the element-incoherent flow observed for different element heights likely result from a modulation of the Kelvin–Helmholtz-like instability, discussed in §3.3, which is essentially independent of the element-induced flow.

3.3 Effect of canopy parameters on the shear-layer instability

The variations observed in the velocity fluctuations for the fixed-element-spacing simulations, discussed above, may result from the growth of the Kelvin–Helmholtz-like, shear-layer instability, typically reported in dense canopy flows (Finnigan, 2000; Nepf, 2012). In order to assess the presence of this instability in the flow, we compare the pre-multiplied spectral energy densities of the wall-normal velocity at $y^{+}\approx 15$ in figures 10($f$–$j$). For case H16 we observe that the spectral energy densities of the fluctuations above the canopy are similar to those above smooth walls. As the height of the canopy is increased, we observe a progressive increase in the energy in long spanwise wavelengths, $\lambda_{z}^{+}>100$, for a narrow range of streamwise wavelengths, $\lambda_{x}^{+}\approx 150$–$250$. This range of streamwise wavelengths remains roughly constant for increasing canopy heights. Such a signature in the spectral energy densities has been previously associated with the presence of spanwise-coherent, Kelvin–Helmholtz-like instabilities over riblets (García-Mayoral & Jiménez, 2011), transitional roughness (Abderrahaman-Elena et al., 2019) and permeable substrates (Gómez-de-Segura & García-Mayoral, 2019). This signature in the spectral energy densities is also reflected in the instantaneous realisations of the wall-normal velocity, portrayed in figure 5, which show increased spanwise coherence with increasing element height. The shear-layer instability is known to generate strong wall-normal fluctuations and, hence, its signature is most clear in the wall-normal spectra (García-Mayoral & Jiménez, 2011; Gómez-de-Segura & García-Mayoral, 2019). Ghisalberti (2009) and Nepf (2012) concluded that canopies only exhibit a shear-layer instability if their height is larger than the wall-normal extent of the instability as otherwise the rollers would be constrained by the lack of canopy depth. Above a short canopy, like that of H16, the proximity of the impermeability condition at the base of the canopy would inhibit the instability by blocking the wall-normal fluctuations. Similarly, Huerre (1983) and Healey (2009) showed that the confinement of a free-shear layer also results in stabilisation of the associated Kelvin–Helmholtz instability. Increasing the canopy height weakens this effect, leading to a stronger signature of the instability, observed in figures 10($f$–$j$). This enhanced signature of the instability is likely responsible for the increase in the cross-velocity fluctuations for the canopies family H with increasing height, discussed in §3.2. For $h/s>6$, the instability no longer perceives the canopy base and the effect of the height on the instability and, consequently, the velocity fluctuations, saturates. The Kelvin–Helmholtz-like instability has also been reported to cause an increase in the Reynolds shear stresses, with an associated increase in the friction drag, over surfaces such as riblets and permeable substrates (García-Mayoral & Jiménez, 2011; Gómez-de-Segura & García-Mayoral, 2019). The increase and saturation of the Reynolds shear stresses with increasing canopy height can be observed in figure 9($e$) for the canopies of family H, and is concurrent with the effect of the element height on the instability, discussed above. This trend in the Reynolds shear stress has a corresponding effect on the drag exerted on the overlying flow, which is illustrated by the downward shift in the mean velocity profiles, portrayed in figure 9($b$). The above discussion suggests that the secondary effect that the height has on the full velocity fluctuations within and above the canopy is mainly through its influence on the Kelvin–Helmholtz-like instability.

The increase in intensity of the Kelvin–Helmholtz-like rollers with increasing canopy height also contributes to the increase in the wall-normal velocity fluctuations within the canopy observed in figure 7($e$). This is demonstrated by the wall-normal spectral energy densities of the flow within the canopies, portrayed at $y^{+}\approx-10$ for all the canopies of family H, in figures 11($a$–$e$). Note that in calculating the spectra for a region with solid obstacles, we have implicitly assumed that the obstacles are fluid regions with zero flow velocity. As discussed in the previous paragraph, for case H16 the instability is inhibited by the proximity of the canopy-base wall, and the flow above shows similarities to a smooth-wall flow. The energy density within the canopy at $y^{+}\approx-10$ for this case also shows some regions overlapping with the smooth-wall spectra, with additional energy in the wavelengths associated with the Kelvin–Helmholtz-like instability. The smooth-wall spectra displayed for reference are at $y^{+}\approx 1$, which is as low as possible while yielding a non-negligible energy, since no direct comparison with $y^{+}\approx-10$ is possible. We also observe some energy in the spanwise wavelength corresponding to the canopy spacing and a broad range of streamwise wavelengths. These regions can be attributed to the modulation of the element-induced flow by the larger scale fluctuations induced by the instability or the overlying turbulence (Abderrahaman-Elena et al., 2019). This suggests that the fluctuations within a short canopy result mainly from the penetration of the overlying turbulence, with additional contributions from the Kelvin–Helmholtz-like instability and the element-induced flow. As the canopy height is increased, and the instability becomes stronger, the deviations in the spectral energy densities from smooth-wall flow become more prominent. The fluctuations within the canopy in cases H32 to H128 arise mainly from large spanwise wavelengths, likely originating from the Kelvin–Helmholtz-like instability near the canopy-tip plane, along with a contribution of the modulated element-induced flow discussed above. The increasing signature of the instability within the canopy with increasing element height can also be observed in the instantaneous realisations of the wall-normal velocity at $y^{+}\approx-10$ portrayed in figure 12. The presence of large spanwise wavelengths deep within the canopy can also be noted for the canopies of family S, whose spectral energy densities and realisations of wall-normal velocity at $y^{+}\approx-40$ are portrayed in figures 11($f$–$j$) and 12, respectively. This suggests that, in the present dense canopies, the background turbulence is not able to penetrate far below the canopy tips, and that the velocity fluctuations deep within originate mainly from the footprint of the Kelvin–Helmholtz-like rollers above.

It is also worth noting that even in canopies with small element spacings, such as that of case S10, although the fluctuations of the wall-parallel velocities decay rapidly below the canopy-tip plane, the fluctuations of the wall-normal velocity decay more slowly, as shown in figure 7. This is also the case for the velocity fluctuations of the canopies of family H, for which the wall-normal velocity fluctuations within the canopy require larger canopy heights to saturate compared to the wall-parallel fluctuations. The presence of wall-normal fluctuations deep within the canopy are a reflection of the canopy layout being able to obstruct the wall-normal flow less efficiently than the tangential flow. It will be shown in §4 that, for the present canopies, the effective drag coefficient in the tangential directions can be up to three times larger than in the wall-normal direction. In the core region of a tall canopy, the only mechanism to inhibit the velocity fluctuations is the canopy drag. As the canopy geometries studied here exert more drag on the wall-parallel flow than the wall-normal flow, $u^{\prime}$ and $w^{\prime}$ decay faster than $v^{\prime}$ within the canopy.

The spacing between the canopy elements affects the excitation of the Kelvin–Helmholtz-like instability through its influence on the canopy drag. White & Nepf (2007) and Sharma et al. (2017) have shown that the canopy drag governs the instability through two competing effects, the shear at the canopy tips and the canopy drag. A small element spacing results in a large drag within the canopy, which in turn results in a larger shear at the canopy tips that enhances the instability, but at the same time it also inhibits the fluctuations more strongly, which weakens the instability. To study this effect, we now compare the pre-multiplied spectral energy densities of the wall-normal velocities for the canopies of family S, which have a constant height and different element spacing. For the canopy with the smallest spacing, S10, the signature of the Kelvin–Helmholtz-like instability in the energy densities is weak and the distribution of energy in different wavelengths is similar to that over smooth walls, as shown in figure 10($a$). This suggests that the large drag exerted on the fluctuations by this canopy inhibits the formation of the shear-layer instability. As the element spacing is increased, the drag on the fluctuations reduces, and there is a progressive increase in the energy in wavelengths associated with the Kelvin–Helmholtz-like instability for the canopies of S16–S48, as can be observed in figures 10($b$–$e$). The change in element spacing also has an effect on the instability wavelength, which will be discussed later in this section. In addition to the increase in the energy associated with the instabilities, the increase in element spacing also results in a progressive decrease in the overlapping regions in the energy densities of the canopy and smooth-wall flows, with a reduction in the energy in wavelengths $\lambda_{x}^{+}\gtrsim 700$, $\lambda_{z}^{+}\approx 50$–$100$. If the element spacing was increased further, the Kelvin–Helmholtz-like instability would eventually weaken, and for sparse enough canopies the flow within would begin to resemble smooth-wall flow perturbed by the element-induced flow of the isolated canopy elements. Such sparse canopies are beyond the scope of the present work, but have been discussed in Sharma & García-Mayoral (2019), as well as in the previous studies of Poggi et al. (2004); Pietri et al. (2009) and Huang et al. (2009).

Let us now focus on the self-similar canopy geometries of family G. Although these canopies have the same $\lambda_{f}$, increasing the size of the canopy elements produces similar effects on the pre-multiplied spectral energy densities as the canopies of family S, discussed in the previous paragraph. For the densest canopy, G10, the spectral energy densities at $y^{+}\approx 15$ collapse with those over a smooth wall, as shown in figures 13($a$–$d$). As the size of the canopy is increased, we observe a stronger signature of the Kelvin–Helmholtz-like instability in the energy densities, portrayed in figures 10($k$–$o$). The associated increase in spanwise coherence in the flow can also be observed in the instantaneous realisations of the wall-normal velocity shown in figure 5. Note that here, we observe the combined effects of increasing canopy height, as in family H, and increasing spacing, as in family S, on the instability. So far we have mainly discussed the spectral energy densities of the wall-normal velocity fluctuations, as they have the strongest signature of the Kelvin–Helmholtz-like instability. For completeness, we now use the canopies of family G to illustrate the effect of increasing the canopy size on the spectral energy densities of the streamwise and spanwise fluctuations, and the Reynolds shear stresses, portrayed in figure 13. The distinct region in the wall-normal spectral energy densities associated with the Kelvin–Helmholtz-like instability is not so apparent in the energy densities of the other velocity fluctuations and the Reynolds shear stresses. Nevertheless, as the canopy size increases, we observe an increase in the energy in streamwise wavelengths associated with the instability, along with a general increase in the energy in shorter and wider wavelengths compared to smooth-wall flows. This is consistent with the gradual shortening and widening of the eddies observed in figure 5.

The results discussed in this section suggest that the growth of the Kelvin–Helmholtz-like instability depends on both the canopy height and the element spacing. The streamwise wavelength of the instability, however, seems to depend mainly on the element spacing. We observe that, in the canopies of family H the streamwise wavelength of the instability is roughly constant regardless of the canopy height, $\lambda_{x}^{+}\approx 150$, as can be observed in figures 10($f$–$j$). For the canopies of family S, however, the increase in the element spacing results in an increase in the streamwise wavelength of the instability from $\lambda_{x}^{+}\approx 140$ for case S10 to $\lambda_{x}^{+}\approx 280$ for case S48, as shown in figures 10($a$–$e$). Similarly, for the canopies of family G there is a progressive increase in the streamwise wavelengths associated with the instability with the increase in canopy size. The streamwise wavelength of the Kelvin–Helmholtz-like instability is determined by the shear length, typically defined in the literature as $L_{s}=U/(dU/dy)$ calculated at the canopy-tip plane (Raupach et al., 1996; Finnigan, 2000; Nepf, 2012). Previous studies have shown that the shear length $L_{s}$ in canopy flows can be determined by the effective streamwise canopy drag coefficient (Finnigan, 2000; Nepf et al., 2007; Ghisalberti, 2009; Nepf, 2012). Intuitively, in tall, dense canopies, we would expect this drag coefficient, and by extension the shear-layer thickness, to be a function of the element spacing. A dependence of the canopy shear-layer thickness on the element spacing was also observed by Novak et al. (2000) in their study of natural canopy flows. Therefore, the canopies of family H, which have a constant element spacing, have similar mean drag coefficients and, consequently, similar instability wavelengths, as observed in the spectral energy densities of the fixed-spacing canopies. For the canopies of families S and G, increasing the element spacing decreases the canopy drag coefficient, thereby resulting in the larger wavelengths observed for the shear-layer eddies. The effect of the canopy spacing on the drag coefficients and the instability wavelengths will be discussed further in §4.

Finally, let us discuss the effect of the Reynolds number on the Kelvin–Helmholtz-like instability. It was shown in §2.1, that the turbulent fluctuations over dense canopies scale in friction units, and therefore similar results are obtained when simulating canopies with the same height and spacing in friction units at different Reynolds numbers. It can be observed in figure 4($b$) that the signature of this instability in simulations H32180 and H32400 are essentially the same, and that the associated streamwise wavelength for both cases is roughly $\lambda_{x}^{+}\approx 150$. As discussed above, the wavelength and amplification of the instability are governed by the shear at the canopy tips. As the canopy parameters for cases H32180 and H32400 are kept constant in friction units, we can also expect the shear at the canopy tips to also be similar. Therefore, the instability characteristics for both these canopies are essentially the same when scaled in friction units. We have observed similar behaviours for Kelvin–Helmholtz-like instabilities originating over riblets (García-Mayoral & Jiménez, 2012) and permeable substrates (Gómez-de-Segura & García-Mayoral, 2019).

4 Linear analysis of Kelvin–Helmholtz-like instabilities

The results from DNS discussed in §3 show that the flow in the region near the canopy-tip plane can be dominated by the presence of spanwise-coherent structures originating from a Kelvin–Helmholtz-like instability. This instability can be captured by a two-dimensional, mean-flow linear stability analysis, even in turbulent flows (Jimenez et al., 2001; White & Nepf, 2007; García-Mayoral & Jiménez, 2011; Zampogna et al., 2016; Gómez-de-Segura & García-Mayoral, 2019). In this section, we discuss the methodology and results from such an analysis conducted on the velocity profiles obtained from the DNSs. As the Kelvin–Helmholtz-like instability is an inviscid phenomenon, several of the studies just cited use an inviscid analysis to capture it. The inclusion of viscosity, however, inhibits the growth of smaller wavelengths in the flow, and consequently, results in the most amplified wavelength being slightly larger compared to that of an inviscid analysis (Jimenez et al., 2001; Gómez-de-Segura & García-Mayoral, 2019). In this section, we present the results only from viscous analysis. The results from an inviscid analysis are presented in appendix B for reference. In addition, we show that some of the key features of this instability can be captured by linear analysis performed on velocity profiles modelled a priori, which would not require any information from the DNSs.

For the purpose of the stability analysis, we model the effect of the canopy using a drag force in the Navier–Stokes equations, which results in the following governing equations

[TABLE]

where $C_{i}$ is the effective canopy drag coefficient in each $i^{th}$ direction, has dimensions of inverse length squared –being essentially the inverse of a permeability–, and is assumed to be homogeneous over the entire canopy region, as in Singh et al. (2016). Given the density of the canopies considered, with maximum spacings $s^{+}=\mathcal{O}(10)$, we assume that inertial effects in the flow deep within the canopy are small and can be neglected. In addition, the element width of the canopies is also small, $w^{+}\approx 1$–$25$. For such canopies, the canopy drag can be assumed to depend linearly on the velocity (Tanino et al., 2005; Tanino & Nepf, 2008).

In the core of the canopy, away from the shear effects at the canopy base and top, the mean momentum equation would reduce to a balance between the canopy drag and the mean pressure gradient

[TABLE]

Equation (9) is essentially Darcy’s equation for flow within permeable substrates (Darcy, 1856), and has been used by Zampogna & Bottaro (2016) to model flow deep within densely packed, rigid fibres. The streamwise drag coefficient, $C_{x}$, can be obtained by substituting the values of $U$ and ${\mathrm{d}P}/{\mathrm{d}x}$ obtained from the DNSs into equation (9). From dimensional arguments, equation (9) predicts that the drag coefficient would scale as $C_{x}\sim 1/s^{2}$. This scaling is demonstrated in figure 14($a$), which suggests that equation (9) provides a reasonable approximation for the flow deep within the present canopies, excluding the sparsest canopy S48. Although we can expect the flow within the canopy to be Darcy-like in the wall-normal direction as well, we cannot use the DNS results to obtain $C_{y}$, as there is no mean flow in this direction. In order to obtain $C_{y}$, we consider separately the Stokes flow along infinitely long canopy elements driven by a constant pressure gradient. The equation for such flow is $\nu(\p^{2}_{x}+\p^{2}_{z})v=\mathrm{d}P/\mathrm{d}y$. The wall-normal drag coefficient is then obtained as

[TABLE]

where the angled brackets represent a spatial average. The estimated values of $C_{y}$ are portrayed in figure 14($b$) for reference. It may be noted that the ratio of the streamwise to wall-normal drag coefficients for the present canopies is $C_{x}/C_{y}\approx 2$–$3$, which shows that the streamwise flow is more obstructed than the wall-normal flow for the layouts considered. It is worth noting here that the canopy drag coefficients, $C_{x}$ and $C_{y}$, and the ratio between them also depends on significantly on the plan view arrangement and the resulting porosity of the canopy (Van der Westhuizen & Du Plessis, 1996; Zampogna & Bottaro, 2016). This is evidenced by the different ratios of the drag lengthscale and the element spacing for the canopies of families S and G, portrayed in figure 14, which have width-to-spacing ratios, $w/s=1/2$ and $2/9$, respectively. The dependence of the $C_{x}$ on $w/s$ can also be predicted using two-dimensional Stokes flow simulations, as shown in figure 14($a$). For more complicated canopy arrangements, such as staggered or random, we would expect the drag to be a function of the planar layout of the elements.

In order to conduct the stability analysis, we linearise the equations (7) and (8) around the mean flow, $U(y)$, yielding

[TABLE]

These equations are used to obtain a modified Orr-Sommerfeld equation (Drazin & Reid, 2004; White & Nepf, 2007; Singh et al., 2016; Zampogna et al., 2016),

[TABLE]

Assuming wavelike solutions of the form $v=\widetilde{v}e^{i(\alpha x-\omega t)}$, equation (14) reduces to the eigenvalue problem

[TABLE]

where the prime superscript denotes differentiation with respect to $y$, and $\mathrm{D}$ represents the operator $\mathrm{d}/\mathrm{d}y$. Equation (15) is solved to obtain the complex frequency, $\omega$, for real values of the streamwise wavenumber, $\alpha$, subject to no-slip and impermeability boundary conditions at the top and bottom walls. The instability is then amplified for positive values of the imaginary part of $\omega$.

The growth rates for different perturbation wavelengths are portrayed in figures 15($a$–$c$), and the wavelengths with the highest growth rates are summarised in table 2. The most amplified wavelengths predicted by the stability analysis only match those observed in the DNSs for canopies with high values of $\delta/h$. The wavelengths predicted for cases H16, H32, G10, G20 and G40 show reasonable agreement with those observed in the DNSs. For canopies with larger heights, however, the analysis predicts wavelengths larger than those observed in the DNSs. For the fixed-spacing canopies of family H, the predicted instability wavelength also increases with increasing canopy height, whereas the DNSs show that the instability wavelength for these cases is essentially independent of the height. The contours of the instability stream function for case H96 for the most amplified wavelength, $\lambda_{x}^{+}\approx 385$, portrayed in figure 16($a$), show that it has a large wall-normal span, extending up to $y^{+}\approx 120$. Such an instability was also reported by Singh et al. (2016), who performed stability analyses similar to the one conducted here, except that the canopy was represented by a drag force depending quadratically on the velocity. Singh et al. (2016) noted that their analysis predicted two instability modes, one similar to the Kelvin–Helmholtz instability and another originating from the canopy drag included in the analysis. They only considered canopies with low $\delta/h$, and observed that the second instability mode, similar to the large-wavelength modes obtained from the present stability analysis, was dominant for canopies with high drag and spanned the entire height of the channel.

It is worth noting here, that in the region near the interface between the canopy and the free-flow, the assumption of a constant drag coefficient given by equation (9) would no longer be valid, as shear and advective effects become stronger. We have conducted some exploratory analyses accounting for this variation in the drag coefficient and these do not provide improved estimates for the instability wavelength compared to the results of the constant-drag analysis presented here. In the present analysis, we have also assumed that the drag coefficient experienced by the perturbations is the same as that experienced by the mean flow. However, we have recently reported for sparser canopies that different wavelengths in the flow can perceive drag coefficients different from that for the mean flow (Sharma & García-Mayoral, 2019). In such a case, the drag coefficient would have to be calculated on a mode-by-mode basis for the different wavelengths. A wavelength dependent drag coefficient would lead to the drag for a given wavelength being dependent on a convolution from all other wavelengths. Such an analysis, however, is beyond the scope of the present work. An alternative approach would be to model the canopy as a permeable substrate, which naturally yields wavelength-dependent equations for the flow within(Zampogna et al., 2016; Abderrahaman-Elena & García-Mayoral, 2017; Sharma et al., 2017; Gómez-de-Segura & García-Mayoral, 2019). In order to illustrate that applying constant drag coefficient on the perturbations may be a coarse assumption, we also present results from stability analyses with no drag on the perturbations. This is also a rather coarse assumption, but we observe that excluding the drag on the perturbations in the stability analysis yields better estimates for the instability wavelengths observed in the DNSs, as shown in figure 10. For the canopies of family H, the stability analysis without drag on the fluctuations shows that the most amplified wavelength does not vary significantly with the canopy height. For the canopies of families S and G, this analysis shows an increase in the most amplified wavelength with increasing element spacing, owing to the increase in the shear-layer thickness. The results from this analysis are portrayed in figures 15($d$–$f$), and the most amplified wavelength for each case is listed in table 2. While neglecting the drag acting on the fluctuations yields better estimates for the most amplified wavelengths for canopies with small spacings, the predictions for larger spacings differ by up to a factor of two from the DNS observations. This is likely due to the assumption that the mean flow is homogeneous in the tangential directions, implicit in the stability analysis, which breaks down for such cases. We have not observed any significant signature of the Kelvin–Helmholtz-like instability for the sparser canopies studied in Sharma & García-Mayoral (2019) despite the presence of an inflection in the mean velocity profiles. There may also be some distortion of the instability by the ambient turbulent fluctuations in the DNSs (Rogers & Moser, 1994; Raupach et al., 1996).

We have also performed stability analyses on the cases with different Reynolds numbers, H32180 and H32400. These analyses predict similar instability wavelengths and growth rates in viscous units, as shown in figure 17, which is consistent with the observations in the DNSs, discussed in §3.3. These results emphasize that the strength of inflection in the mean-velocity profile and the shear-layer thickness for both these cases is similar when scaled in friction units and, hence, so is their effect on the instability.

The results obtained from the DNSs and the stability analysis suggest that there is a dependence on the element spacing of the most amplified wavelength, related to the effect of the spacing on the shear-layer thickness. The usual definition of the shear-layer thickness, $L_{s}=U/(dU/dy)$, misses the contribution of the part of the shear layer above the canopy. Regarding the latter, García-Mayoral & Jiménez (2011) studied the formation of Kelvin–Helmholtz-like instabilities over riblets, and noted that the shear-layer thickness above was given by the height at which the vorticity gradient, $\mathrm{d}^{2}U/\mathrm{d}y^{2}$, concentrated. In smooth-wall flows, this height is roughly $y_{c}^{+}\approx 5$–$10$. For the present cases, we observe that the instability wavelengths predicted by the stability analysis correlate well with the full shear length $L_{s}+y_{c}$, if we take $y_{c}^{+}\approx 5$, as shown in figure 18($a$). This suggests that the shear-layer semi-thickness above the canopies, $y_{c}$, is roughly constant for most of the geometries considered here, and remains close to the smooth-wall value, while the semi-thickness below has the standard form $L_{s}=U/(dU/dy)$, measured at the canopy tips plane. The only notable deviation is for the sparsest canopy studied, S48. For canopies with large element spacings, we observe that the peak in $\mathrm{d}^{2}U/\mathrm{d}y^{2}$ moves closer to the canopy-tip plane, so $y_{c}^{+}=5$ may no longer be a reasonable approximation for the shear-layer semi-thickness above. Regarding the height of the shear layer within the canopy, $L_{s}$, we observe that it is set by the mean canopy drag coefficient, $L_{s}\propto\sqrt{1/C_{x}}$, as also noted in the studies of aquatic canopy flows by Nepf et al. (2007) and White & Nepf (2007). The drag coefficient on the mean flow, in turn, depends on the element spacing and the width-to-spacing ratio, as shown in figures 18($b$) and ($c$). The correlation of $L_{s}$ with $s$, therefore, explains the dependence of the most amplified wavelength on the element spacing observed in the DNSs and the stability analysis.

4.1 Analysis on modelled velocity profiles

In this section we introduce a simple model for the mean velocity profile in dense canopy flows and discuss the results from their stability analysis. As we only consider canopies with small element spacings, $s^{+}=\mathcal{O}(10)$, the magnitude of inertial effects within the canopies are also small and are thus neglected in the model. The results discussed in §3 also suggest that, for very dense canopies, turbulence and, consequently, the Reynolds shear stresses do not penetrate within (Nepf et al., 2007), and are smooth-wall-like above the canopy-tip plane. The mean velocity above the canopy could then be modelled using a smooth-wall eddy viscosity, with the canopy-tip plane acting as the location of the smooth-wall (Jimenez et al., 2001; García-Mayoral & Jiménez, 2011; Gómez-de-Segura et al., 2018; Gómez-de-Segura & García-Mayoral, 2019). The equation for the mean velocity can then be written as

[TABLE]

where $C_{x}(y)$ is the average streamwise canopy drag coefficient, which is assumed constant within the canopy and zero outside, and $\nu_{T}(y)$ is the height-dependent eddy viscosity proposed by Cess (1958) to approximate turbulent smooth-channel flow, and is non-zero only outside the canopy. The drag coefficients, $C_{x}$, used to obtain the velocity profiles are those given by equation (9) and portrayed in figure 14($a$). These have been obtained using the data from the direct numerical simulations, but can also be obtained from Stokes-flow simulations as shown in figure 14($a$), which are significantly less computationally intensive.

The most amplified wavelengths predicted by the stability analysis conducted on these modelled velocity profiles, with no drag applied on the fluctuations, are in reasonable agreement with those obtained from the same no-drag analysis on profiles obtained from the DNSs. The growth rates predicted are portrayed in figures 15($g$–$i$). The wavelengths with maximum growth rates are also summarised in table 2. We have also conducted stability analyses on these modelled velocity profiles including the effect of the eddy viscosity. The results are portrayed in appendix B, and they are essentially the same as the ones obtained using molecular viscosity alone, apart from a slight reduction in the instability growth rates, which suggests that although $\nu_{T}$ is important for setting the shape of the mean velocity profile, its effect on the fluctuations is not significant. This is likely because the Kelvin–Helmholtz-like rollers occur near the canopy tips, where $\nu_{T}$ is small and the molecular viscosity, $\nu$, dominates. It is worth noting that even though this model is able to capture the instability wavelength, the velocity profiles obtained using this model do not match those from the DNSs, apart from those of S10 and G10. This is most likely due to our assumption that the turbulent stresses do not penetrate within the canopy and remain smooth-like, which fails as the element spacing is increased. As discussed previously, the wavelength of the instability is set by the shear length. The shear-layer semi thickness within the canopy, $L_{s}$, is set by the canopy drag coefficient, $C_{x}$. As this drag coefficient is the same both from the DNSs and for the modelled velocity profiles, we expect $L_{s}$ to be similar as well. The shear length above the canopy, however, could differ, as the profiles from DNS would include the effect of the turbulent stresses penetrating into the canopy and deviating from their smooth-wall values, while the modelled velocity profiles do not. The similarity in the instability wavelengths between these analyses therefore suggest that, for most of the dense canopies considered in this work, turbulence is essentially precluded from penetrating into the canopy, and that the shear length above does not vary significantly from its smooth-wall value.

5 Conclusions

In the present work, we have examined the effect of the canopy layout on turbulent flows over canopies of densely packed filaments of small size. Three families of simulations have been conducted, the first with the element height in friction units fixed, the second with the element spacing fixed, and the third with the height-to-spacing ratio fixed. The layouts considered had height-to-spacing ratios greater than one, and elements spacings in the range $s^{+}\approx 3$–$50$. The penetration of turbulent fluctuations within such canopies was limited by their small element spacings. Consequently, the height of the roughness sublayer was also determined by the element spacing, rather than their height, extending up to $y\approx 2$–$3s$ above the canopy tips. The canopy drag coefficient was also found to be determined by the element spacing, $C_{x}\sim 1/s^{2}$. Canopies with small spacings were, therefore, found to suppress the velocity fluctuations within them owing to the large drag exerted, and the fluctuations became more intense as the spacing increased. The intensity of the characteristic Kelvin–Helmholtz-like instability over canopies was observed to be governed by two competing effects resulting from the canopy drag, the inflection at the canopy tips and the drag on the fluctuations. Canopies with large drag had a large shear at the canopy tips, and thus a stronger inflection, which enhanced the instability, but also exerted a large drag on the velocity fluctuations, which suppressed the instability. The instability was found to be inhibited in canopies with $s^{+}\lesssim 10$ and, for the range of canopy spacings considered here, a stronger signature of the instability was observed as the spacing was increased. We also showed that the main contribution to the velocity fluctuations deep within the canopy was the footprint of the Kelvin–Helmholtz-like instability, and that the contribution of the element-induced dispersive flow was negligible. Short canopies with $h/s\sim 1$ were also found to inhibit the instability, owing to the blocking effect of the wall at the canopy base. For height to spacing ratios $h/s\gtrsim 6$, the instability was no longer influenced by the bottom wall, and the effect of the canopy height on the flow within and above the canopy saturated. Increasing the canopy height for a fixed spacing did not change the element-induced velocity fluctuations, and instead affected the surrounding flow through the influence of height on the instability.

Linear stability analysis conducted on the mean velocity profiles obtained from the DNSs is able to capture the approximate wavelength of the instability observed in the DNSs for canopies with small element spacings. The analysis fails for larger element spacings, for which the assumption of the flow perceiving the canopy in a homogenised fashion breaks down. We showed that the shear-layer thickness, which determines the instability wavelength, has two components, one within the canopy and the other above. The latter is set by the height above the canopy tips at which the vorticity gradient concentrates, and is essentially constant for the present canopies, $y_{c}^{+}\approx 5$. The shear layer thickness within the canopy follows the conventional definition, $L_{s}=U/(dU/dy)$, and is determined by the canopy drag, thus depending linearly on the canopy spacing. We have also proposed a simplified model to capture the most amplified instability wavelength over dense canopies. The model assumes that the turbulence above the canopy does not penetrate within and remains smooth-wall-like, and uses the mean streamwise drag coefficient of the canopies to synthesise an approximate mean-flow profile. The stability analysis conducted using these synthesised profiles yields similar results to those conducted using the mean profiles from DNS.

A.S. was supported by an award from the Cambridge Commonwealth, European and International Trust. Computational resources were provided by the “Cambridge Service for Data Driven Discovery” operated by the University of Cambridge Research Computing Service and funded by EPSRC Tier-2 grant EP/P020259/1. The authors report no conflicts of interest.

Appendix A Immersed boundary method and validation

For the present work, a modified version of the immersed boundary algorithm proposed by García-Mayoral & Jiménez (2011) is used. The algorithm of García-Mayoral & Jiménez (2011) was based on a direct-forcing approach, which applies a body force within the immersed-boundary points to drive the velocity at these points to zero (Mittal & Iaccarino, 2005). The condition to implement at the points within the canopy elements is

[TABLE]

Following García-Mayoral & Jiménez (2011), this condition can be approximated by modifying the right-hand-side of equation (1)

[TABLE]

The original code of García-Mayoral & Jiménez (2011) used a collocated grid for the wall-normal coordinate and was extended by Fairhall & García-Mayoral (2018) to employ a staggered grid. Fairhall & García-Mayoral (2018) also split the Laplacian operator on the left-hand-side of equation (1) into its wall-parallel and wall-normal components following Kim & Moin (1985)

[TABLE]

where $L_{xz}$ includes the wall-parallel components of $L$, and $L_{y}$ the wall-normal one. Splitting the Laplacian in this manner still retains the second-order temporal accuracy of the code (Kim & Moin, 1985). Equation (1) can then be written as

[TABLE]

In the present work, we implement a modified version of the immersed boundary algorithm used by García-Mayoral & Jiménez (2011) into the above algorithm, which offers an improvement in the accuracy for the velocities within the immersed boundary regions. This implementation is summarised below. The right-hand-side of equation (20) is then transformed to physical space, and modified to satisfy the following conditions within the immersed boundary points

[TABLE]

where $y_{tips}$ denotes the wall-normal plane at the canopy tips. At the interfaces, the condition imposed by equation (21) yields $\mathbf{u}^{n}_{k}=\mathcal{O}(\Delta t^{2})$, which is of the order of the temporal discretisation error of the code. As a staggered grid is used, the wall-normal grid points for the streamwise velocities are offset by half a grid spacing from those of the wall-normal velocity. The element tips are aligned with the grid for the streamwise velocity. For the wall-normal velocity, the interface condition is set at the grid point just below the canopy-tip plane, which enforces near-zero wall-normal velocity at the canopy tips through continuity. Away from the interfaces, within the immersed boundary region, the condition set by equation (21) results in an exponential decay of the velocity in the wall-normal direction from its $\mathcal{O}(\Delta t^{2})$ value at the interface. To illustrate this decay, the results from a simple, one-dimensional implementation of this algorithm are shown in figure 19. The velocity is assumed to vary only in the wall-normal direction. For these test cases, random initial conditions are used to mimic turbulent fluctuations in the flow. The velocity fields obtained for these cases, after one timestep, using the immersed boundary conditions of equation (21) are compared to those obtained using equation (18) in figure 19.

Although both algorithms result in small velocities within the solid regions, the implementation proposed here results in a smoother and much faster decay of the velocity within the solid. In the DNS code, however, the velocity correction step introduces an error of order $\Delta t^{2}$ in all the immersed boundary points. Even so, in experience it was observed that the proposed algorithm is a more stable numerical implementation of the immersed boundaries compared to the one proposed by García-Mayoral & Jiménez (2011). This is likely due to the present method not generating sharp gradients in the velocity field within the solid obstacles at the pressure calculation step. The velocity within the canopy elements, or the permeability error, in the DNSs is observed to be less than $0.1u_{\tau}$, for all the conducted simulations, and is much smaller than velocity in the ‘fluid’ points surrounding the elements. This is illustrated in the instantaneous realisations of the velocity fields from one of the simulations are portrayed in figure 20. For completeness, the wall-normal grid distribution for case S10 is portrayed in figure 21. In order to validate the implementation of the immersed boundaries, we have replicated the DNS for the collocated roughness elements with height $h^{+}\approx 12$ of Abderrahaman-Elena et al. (2019). The mean velocity profiles, rms fluctuations and the Reynolds shear stresses obtained from the present simulations show good agreement with the results of Abderrahaman-Elena et al. (2019), as shown in figure 22. We have also conducted a grid-dependence analysis, for which the results are portrayed in figure 23. Velocity rms fluctuations for cases S48 and G100 are shown for different wall-parallel resolutions.

Appendix B Comparison of inviscid and viscous stability analysis

Here, we provide the governing equations used to perform a stability analysis with a turbulent viscosity varying in the wall-normal direction, $\nu_{T}(y)$. The modified Orr-Sommerfeld equation is then given by

[TABLE]

A similar equation, excluding the canopy drag terms, has also been used by Reynolds & Hussain (1972), Del Alamo & Jiménez (2006), Pujals et al. (2009) and Gómez-de-Segura et al. (2018). In figures 24($a$–$f$), we compare the results obtained from viscous and inviscid analysis conducted using the velocity profiles from the DNSs. In In figures 24($g$–$i$), we show the results from inviscid and viscous stability analyses, using both molecular and turbulent viscosities, performed using the velocity profiles obtained from equation (16).