Numerical study of fractal-tree-generated non-equilibrium turbulence

Abstract

Self-similar fractal tree models are numerically investigated to elucidate the drag coefficient, non-equilibrium dissipation behaviour and various turbulence statistics of fractal trees. For the simulation, a technique based on the lattice Boltzmann method with a cumulant collision term is used. Self-similar fractal tree models for aerodynamic computations are constructed using parametric L-system rules. Computations across a range of tree-height-based Reynolds numbers $Re_H$, from 2500 to 120 000, are performed using multiple tree models. As per the findings, the drag coefficients ($C_D$) of these models match closely with those of the previous literature at high Reynolds numbers ($Re_H \geqslant 60\,000$). A normalization process that collapses the turbulence intensity across various tree models is formulated. For a single tree model, a consistent centreline turbulence intensity trend is maintained in the wake region beyond a Reynolds number of 60 000. The global and local isotropy analysis of the turbulence generated by fractal trees indicates that, at high Reynolds numbers ($Re_H \geqslant 60\,000$), the distant wake can be considered nearly locally isotropic. The numerical results confirm the non-equilibrium dissipation behaviour demonstrated in previous studies involving space-filling fractal square grids. The non-dimensional dissipation rate $C_\epsilon$ does not remain constant; instead, it becomes approximately inversely proportional to the local Taylor-microscale-based Reynolds number, $C_\epsilon \propto 1/Re_\lambda$. We find significant one-point inhomogeneity, production and transverse transport of turbulent kinetic energy within the non-equilibrium dissipation near wake region.

1. Introduction

The fluid mechanics characteristics of trees have gained considerable attention owing to their crucial role in urban micrometeorology simulation (Takahashi et al. 2013; Kamiya et al. 2019), particularly for the open vegetation environments in urban streets (Gromke & Blocken 2015a , b ; Oshio, Kiyono & Asawa 2021). Most studies conducted on trees have focused on calculating the drag coefficient and aerodynamic porosity, both experimentally (Hagen & Skidmore 1971; Grant & Nickling 1998; Gillies, Nickling & King 2002; Guan, Zhang & Zhu 2003; Dong et al. 2008; Bitog et al. 2011; Manickathan et al. 2018) and numerically (Wilson 1985). This is because trees are considered to have a good blocking effect on airflow due to their large size and complex three-dimensional (3-D) geometry. However, few studies have focused on wake turbulence, such as dissipation characteristics. Obtaining a better understanding of the turbulent characteristics of the wake region of trees can help us perform an accurate simulation of the influence of trees on the flow field in the urban environment and in micrometeorological prediction. Therefore, it is essential to comprehensively analyse the turbulent characteristics of the wake region of trees.

Owing to the complex 3-D structure of trees, it can be expected that the wake turbulence is complex and the wake region near the tree (near wake region) is highly inhomogeneous. Turbulence generated by fractal square grids, which are representative of complex-shaped objects, has been reported to exhibit non-equilibrium dissipation properties. This non-equilibrium dissipation anomaly, which contradicts the classical dissipation theory, has been widely studied in the case of planar grid turbulence, which provides a valuable reference for our study of the wake turbulence of trees; therefore, it must be reviewed. In 2007, Hurst & Vassilicos (2007) conducted the first experiments on turbulence generated by planar fractal grids. Seoud & Vassilicos (2007) investigated this turbulence in a wind tunnel and discovered unusual properties, including a non-equilibrium dissipation anomaly that contradicts the classical dissipation theory. This discovery sparked ongoing interest in fractal multiscale grid turbulence, leading to extensive experimental studies (Seoud & Vassilicos 2007; Mazellier & Vassilicos 2010; Suzuki et al. 2010b ; Valente & Vassilicos 2011a ,b; Krogstad & Davidson 2011; Discetti et al. 2011; Valente & Vassilicos 2012a ,b; Krogstad & Davidson 2012; Gomes-Fernandes, Ganapathisubramani & Vassilicos 2012; Discetti et al. 2013; Hearst & Lavoie 2014b ; Valente & Vassilicos 2014; Hearst & Lavoie 2014a ,2015; Gomes-Fernandes, Ganapathisubramani & Vassilicos 2015; Laizet, Nedić & Vassilicos 2015; Nagata et al. 2017) and numerical studies (Nagata et al. 2008a ,b; Laizet & Vassilicos 2009; Laizet, Lamballais & Vassilicos 2010; Suzuki et al. 2010a ; Laizet & Vassilicos 2011; Laizet et al. 2012; Laizet & Vassilicos 2012; Zhou et al. 2014). One of the aims of the studies of fractal-generated turbulence (Seoud & Vassilicos 2007; Nagata et al. 2008a ,b; Mazellier & Vassilicos 2010; Valente & Vassilicos 2011a ,b) was to generate unexplored turbulent flow conditions at high Reynolds numbers that differ from the conventional energy dissipation mode $\epsilon \sim k^{3/2} /l$ . Here, $\epsilon$ represents the energy dissipation rate, $k$ denotes the turbulent kinetic energy (TKE) and $l$ denotes some local correlation length scale. This equation has been deemed ‘a cornerstone assumption of turbulence theory’ (Tennekes & Lumley 1972).

According to Kolmogorov’s similarity hypotheses (universal equilibrium hypotheses) (Kolmogorov 1941), the small-scale turbulent motions are universally similar if they are sufficiently isolated from the influence of boundaries and the energy input scale. These small scales are where the turbulence is isotropic and homogeneous. That is, in all turbulent motions, including boundary layer flows, turbulent wake, pipe flows, etc., the small-scale motions have the same statistical properties at sufficiently high Reynolds numbers. In particular, the relation $\epsilon \sim k^{3/2} /l$ has been widely used in the modelling of various turbulent flows, such as decaying homogeneous turbulence (Batchelor 1953), shear layers, jet flow, wakes (Townsend 1956) and some Reynolds-averaged Navier–Stokes (RANS) models. Specifically, one-equation RANS models which directly use the $\epsilon \sim k^{3/2} /l$ relationship to calculate the energy dissipation, such as Prandtl’s one-equation model, may produce significant errors when calculating non-equilibrium dissipation turbulence. Conversely, some two-equation RANS models can provide better accuracy; for instance, the well-known $k - \epsilon$ model (Launder & Spalding 1974). In a recent study conducted on grid turbulence by Bos (2020), the two-equation RANS simulation yielded good results for kinetic energy and successfully fit the curve by the $x$ position where the energy peaks and the peak value. This study also elucidated the influence of the initial conditions.

Seoud & Vassilicos (2007) first discovered a long wake region of fractal square grids where the turbulence is approximately isotropic but the $C_\epsilon$ = constant does not hold, i.e. the dissipation behaviour is not in a state of equilibrium. Here, $C_\epsilon$ denotes a non-dimensional dissipation rate. According to Taylor (1935), if the turbulence is isotropic, the energy dissipation rate can be expressed by $\epsilon = 15 \nu u ^{\prime2} / \lambda ^2$ . Here, $\nu$ is the kinematic viscosity, $\lambda$ is the Taylor microscale, $u$ is the velocity fluctuation in the $x$ -direction and $^{\prime}$ represents the root-mean-square. In the case of isotropic turbulence, $\epsilon \sim k^{3/2} /l$ becomes $\epsilon = C_\epsilon u^{\prime3} / L$ , because in isotropic turbulence $u^{\prime}=v^{\prime}=w^{\prime}$ (where $v$ and $w$ are the velocity fluctuations in the $y$ - and $z$ -directions, respectively). Here, $L$ denotes the integral length scale. Applying the relation derived by Taylor, $\epsilon = C_\epsilon u^{\prime3} / L$ is equivalent to $15 L/ \lambda = C_\epsilon Re_\lambda$ , where $Re_\lambda$ denotes the Taylor-microscale-based Reynolds number.

Gomes-Fernandes, Ganapathisubramani & Vassilicos (2012) performed flume experiments of fractal square grids using particle image velocimetry measurement, and introduced a method to estimate the location of $x_{peak}$ along the centreline. Here, $x_{peak}$ represents the $x$ -coordinate where the turbulence intensity reaches its maximum, and it marks the boundary between the turbulence-developing and turbulence-decaying regions in the wake of the fractal square grid. In their study, $x_{peak}$ was used to normalize centreline statistics, and they found that the decay of statistical metrics such as turbulence intensity and $Re_\lambda$ occurs beyond $x_{peak}$ , the relationship $L / \lambda \approx$ constant and $C_\epsilon \propto 1 / Re_\lambda$ also occurs beyond $x_{peak}$ .

Valente & Vassilicos (2012 a) investigated turbulence generated by regular grids and identified a transition from non-equilibrium dissipation to $C_\epsilon$ = constant. They discovered a long $C_\epsilon \sim 1 / Re_\lambda$ non-equilibrium dissipation region, extending up to $x \approx 5x_{peak}$ in the wake region. Simultaneously, they observed that at a location that was sufficiently downstream (when $x \geqslant 5x_{peak}$ ), the flow field returned to a state where $C_\epsilon$ remained approximately constant. This indicates that within the turbulence-decaying region, there are also two regions with significantly different properties: a non-equilibrium dissipation region near the upstream and a $C_\epsilon$ = constant region farther downstream.

Hearst & Lavoie (2014a ) confirmed the existence of these two regions in turbulence generated by grids made of many repetitions of a fractal square grid and observed significant inhomogeneity in one-point statistics within the near wake region. Specifically, they reported that the boundary between these two regions was approximately $x /L_0 \approx 22$ (where $L_0$ denotes the size of each square fractal element). Furthermore, they observed that in the non-equilibrium dissipation region, the flow field was highly inhomogeneous, while after entering the $C_\epsilon$ = constant region, the flow field became approximately homogeneous. There is reportedly an approximate correlation between inhomogeneity in one-point statistics and non-equilibrium dissipation in turbulence generated by fractal square grids. It is important to note that this correlation is not universal. Non-equilibrium dissipation scaling has also been found in direct numerical simulation (DNS) of periodic/homogeneous forced turbulence (Goto & Vassilicos 2015) and periodic/homogeneous freely decaying turbulence (Goto & Vassilicos 2016b ). Additionally, DNS of periodic/homogeneous freely decaying turbulence shows a non-equilibrium dissipation region in time followed by a region where $C_\epsilon$ is constant (Goto & Vassilicos 2016b ; Valente, Onishi & da Silva 2014).

An inhomogeneous flow field typically leads to significant production and transverse transport of TKE, as observed in the production and triple-correlation transport terms in the TKE equation. There is an approximate correlation between the non-zero production and triple-correlation transport terms and non-equilibrium dissipation (Valente & Vassilicos 2011b ; Nagata et al. 2013; Hearst & Lavoie 2016). Specifically, Valente & Vassilicos (2011b ) analysed the TKE budget in the non-equilibrium dissipation region behind a fractal square grid and reported significant transverse transport of the TKE. Similarly, Hearst & Lavoie (2016) analysed the TKE budget in the non-equilibrium dissipation region behind grids made of many repetitions of a fractal square grid (the same as used by Hearst & Lavoie (2014a )). They observed that, in the flow field closer to the upstream, the relevant production term and the triple-correlation transport term are non-negligible compared with the dissipation term. It should be noted that these results are based on one-point statistical analysis and cannot be directly used to explain the non-equilibrium and equilibrium dissipation since one-point inhomogeneities do not affect the two-point energy equation on the scales for which the inertial energy cascade occurs (Valente & Vassilicos 2015). Nagata et al. (2017) pointed out that the cascade mechanism of turbulence dissipation is not affected by one-point inhomogeneities.

Non-equilibrium dissipation has been observed across various turbulence. For instance, Valente & Vassilicos (2011b ) and Hearst & Lavoie (2016) found significant TKE production and transverse transport in the non-equilibrium near wake region of fractal square grids. The transition from non-equilibrium dissipation to $C_\epsilon$ = constant in the wake region has been observed with regular grids (Valente & Vassilicos 2012a ; Nagata et al. 2017). Additionally, the non-classical scaling has been observed with normal grids (Krogstad & Davidson 2011), bluff bodies (Obligado, Dairay & Vassilicos 2016) and the instantaneous fluctuations in DNS (Goto & Vassilicos 2015). Fractal grids are considered to amplify the perturbations around the equilibrium, since the largest structures of a fractal grid enable the perturbations to be observed better than a classical grid with the same solidity. However, the non-classical effects (non-equilibrium dissipation) are not due to the fractal nature of the turbulence generator (Horiuti & Tamaki 2013; Bos & Rubinstein 2017), and are probably a common phenomenon in turbulence. Based on these insights obtained from previous studies, an important question arises: Do trees with more complex fractal structures induce highly inhomogeneous, non-equilibrium wake turbulence? How extensive is such a non-equilibrium dissipation region?

Another valuable perspective is from Yoshizawa (1994), the concept of a non-equilibrium spectrum. Bos & Rubinstein (2017) applied this concept to present a framework which allows one to interpret the non-equilibrium scaling. For the ideal case of perfectly isotropic turbulence, the scaling can be related to subdominant corrections to Kolmogorov’s equilibrium spectrum, as predicted by Yoshizawa (1994) and assessed in closure and DNS (Horiuti & Tamaki 2013; Fang & Bos 2023). It is worth noting that in another related study, Goto & Vassilicos (2016a ) used Yoshizawa’s spectrum incorporating non-equilibrium effect but derived a $C_\epsilon$ expression different from that of Bos & Rubinstein (2017). This might suggest that Yoshizawa’s concept of a non-equilibrium spectrum in itself is not sufficient to derive the non-equilibrium dissipation scaling. Further data are required to validate these theoretical perspectives, especially in the context of tree-generated turbulence, for which experimental and numerical data remain scarce.

Due to the large size of trees, conducting high-Reynolds-number wind tunnel experiments is challenging. While scale models, such as Bonsai trees, can be used to obtain drag coefficients, investigating non-equilibrium turbulence may not be straightforward because the Kolmogorov scale also decreases correspondingly, requiring very fine spatiotemporal resolution to capture it accurately. Furthermore, the complex physical shape of tree also incurs a high cost of performing high-Reynolds-number simulations of large-eddy simulations or DNS, especially with conventional numerical computation methods based on multi-central-processing-unit (multi-CPU) parallelization. Therefore, in research fields such as environmental engineering, 3-D fractal trees are typically simulated using the two-equation RANS models (Gromke & Blocken 2015a ,b; Oshio et al. 2021). Although some two-equation RANS models present better accuracy than one-equation models, without DNS or quasi-DNS data for reference, it is difficult to determine the accuracy of the tree calculations, such as the energy dissipation $\epsilon$ . Reliable DNS results can serve as a valuable database for RANS simulations of trees, aiding in validating $\epsilon$ and providing useful references.

In this study, we investigated the turbulence generated by fractal trees under high-Reynolds-number conditions using quasi-DNS simulations (the mesh size in the wake region is $\Delta x / \eta \approx 2.0$ , where $\eta$ denotes Kolmogorov scale) based on a cumulant lattice Boltzmann method (LBM) (Geier et al. 2015) on multi-graphics-processing-unit (multi-GPU) parallelization. Since it is an entirely explicit solver, it is suitable for large-scale computations that require supercomputers. The adaptive mesh refinement (AMR) method was also introduced. The highest-resolution mesh was applied near the boundary layer of tree branches, a higher-resolution mesh in the wake region and a lower-resolution mesh in the remaining unimportant flow field areas. To render the analysis more reproducible and reducible, the tree models generated by the parametric L-system, as employed in Segrovets, Onishi & Kolomenskiy (2022), were utilized in this investigation. To our knowledge, there has been only one DNS study on 3-D fractal-tree-generated flow to date. Segrovets et al. (2022) used DNS to study the effect of tree shape on the drag coefficient at $Re_H=2500$ . More DNS or quasi-DNS studies of 3-D fractal tree flows under high Reynolds number conditions are still required.

The rest of this paper is divided into four sections. Section 2 describes the numerical method, the parametric L-system and the mesh configuration of the boundary layer and wake region. Section 3 presents the results of the simulation validation, confirming the mesh convergence of the drag coefficient and detailing the procedure for determining the energy dissipation rate. Section 4 details the results, comparing the calculated drag coefficients $C_D$ with those measured in previous studies, collapsing statistics from different tree models, providing global and local isotropy parameters at various Reynolds numbers and confirming the non-equilibrium dissipation behaviour of fractal trees at high Reynolds numbers. Finally, § 5 concludes the study.

2. Numerical method

2.1. Lattice Boltzmann method

A code based on the LBM with a cumulant collision term was employed for the numerical computations. This code has been successfully utilized for calculating turbulence around objects with complex geometries (Hasegawa et al. 2019; Ohashi et al. 2021; Watanabe, Hu & Aoki 2023). The LBM is a weakly compressible approximate solution method for incompressible fluids and is widely used for large-scale calculations of isothermal, single-phase incompressible fluids. It treats the fluid as a collection of virtual particles and computes the temporal evolution of their velocity distribution function. To accurately simulate flows at high Reynolds numbers, the D3Q27 cumulant collision model (Geier et al. 2015) was utilized. This model has high numerical stability and computational accuracy. It ensures minimal numerical viscosity when the mesh resolution is sufficiently high, enhancing the simulations’ accuracy. The LBM is fully explicit and does not involve iterative calculations such as Poisson’s equation for pressure in the conventional difference method, thus achieving high computational performance in large-scale calculations.

2.2. Adaptive mesh refinement

A high-resolution grid is necessary to accurately resolve a fractal tree’s boundary layer and wake region. However, employing a uniform high-resolution mesh throughout the entire computational domain would lead to numerous grid points, which would require impractical amounts of computational resources. To maintain accuracy and reduce computational costs, an AMR method was utilized. The AMR method utilizes the octree data structure. In other words, with this structure, high-resolution grids can be allocated to arbitrary regions by recursively subdividing the computational grids (Wahib, Maruyama & Aoki 2016). The finest grid was allocated to the region near the tree surface, and a marginally finer grid was allocated to the wake region. A coarse grid was applied to other distant areas. In the simulation program, most computational subroutines were executed on the GPU, while the CPU handled mesh generation and data output. The LBM fluid calculations were executed at the terminal leaf of the octree. To ensure continuity of memory access for GPU-based calculations, an $8\times 8\times 8$ cell was placed at each leaf. In addition, to perform computations with multiple GPUs, the region was partitioned using a space-filling curve to balance the computational load and memory. More details about the AMR–LBM can be found in a previous paper (Watanabe & Aoki 2021).

2.3. Interpolated bounce-back method

The LBM utilizes orthogonal grids for computation. For object shapes that do not conform to grid directions, like the surface of a fractal tree, we avoided the traditional stair-stepped representation of object shapes. Instead, a second-order accurate interpolated bounce-back method (Bouzidi, Firdaouss & Lallemand 2001) was introduced to establish the boundary conditions. By employing the D3Q27 type of LBM, we could precisely represent the object’s surface in 26 directions. This capability substantially enhanced the accuracy of the fractal tree model. The fluid forces acting on the object surface were calculated based on the momentum exchange (Wen et al. 2014) between the velocity distribution function and the object at the object boundary. At each grid point, the momentum change of the velocity distribution function that bounced back at the object boundary was calculated and then integrated over the object surface to obtain the forces acting on the object.

2.4. Fractal tree model

In this study, a unique algorithm called the L-system (Prusinkiewicz et al. 1996) was chosen to produce fractal tree representations. The L-system is a graphically interpreted parametric algorithm for representing fractal forms from a small set of parameters.

2.4.1. Parametric L-system

Using the parametric L-system, some important parameters, such as the number of branch-generating generations, length, diameter and angle of branches, can be determined. The algorithm of the parametric L-system is as follows:

(2.1) \begin{align} \begin{split} \omega &: A(100, w_0),\\ p_1 &: A(s, w) : s \geqslant min \, \rightarrow \,\, !(w)F(s)\\ [&+(\alpha _1)/(\varphi _1)A(s*r_1,w*q^e)]\\ [&+(\alpha _2)/(\varphi _2)A(s*r_2,w*(1-q)^e)]. \end{split} \end{align}

Here, the first line of (2.1) represents the parameters of the trunk of a tree. The letter $A$ indicates the apex, i.e. a state of not-yet-produced branch (cylinder) with a non-dimensional length of 100 and a diameter of $w_0$ . The second and subsequent lines represent branch-generation rules. The letter $p_1$ signifies replacement of apex $A$ with an internode $F$ (branch) and two new apices $A$ . If the length $s$ of an apex is larger than or equal to the limit minimum value, it generates a branch according to the parameters $(s,w)$ contained in that apex. It also generates two new apices with parameters augmented by $s*r_1,s*r_2,w*q^e$ and $w*(1-q)^e$ for the consecutive branch to be generated. Orientation of the new branch’s apices with respect to the previous branch is determined by rotations $+$ and $/$ . The axis conventions are defined in figure 1, where $\boldsymbol {H}$ represents the direction of the last branch; in other words, it is a vector that passes in the longitudinal direction of the cylinder. Here $\boldsymbol {L}$ and $\boldsymbol {U}$ represent the direction to the left and up. These vectors are perpendicular to each other and $\boldsymbol {H} \times \boldsymbol {L} = \boldsymbol {U}$ .

Figure 1. Axis conventions ( $+, -, \&, ^{\wedge }, \backslash , /$ ) for L-system.

The replacement rule $p_1$ was iterated $n$ times for the specified number of branches generated, and then the 3-D geometric data of the fractal tree can be generated. Specific definitions of each sign used in (2.1) are: $A$ , apex, state of not-yet-produced branch, possesses the parameters ( $s, w$ ) and implicit vector parameter $\boldsymbol {P}$ ; $\boldsymbol {P_j^i}$ , vector $\in$ $A$ , specifies origin and orientation via unit vectors $\boldsymbol {U}, \boldsymbol {L}, \boldsymbol {H}$ ; $*$ , multiplication; index $i$ distinguishes between different branch generations; index $j$ distinguishes between branches of the same branch generation; $F$ , branch, if the condition of $p_1$ is met branch and two new apices $A$ become produced; $min$ , threshold limit of branch length; $\omega$ , axiom, used to seed the string rewriting algorithm; $s$ , branch length; $w$ , branch diameter; $\alpha ,\varphi$ , rotation angles; $+$ , clock wise rotation about $\boldsymbol {U}$ ; $-$ , anticlockwise rotation about $\boldsymbol {U}$ ; $\backslash$ , clockwise rotation about $\boldsymbol {H}$ ; $/$ , anticlockwise rotation about $\boldsymbol {H}$ ; &, clockwise rotation about $\boldsymbol {L}$ ; $^{\wedge }$ , anticlockwise rotation about $\boldsymbol {L}$ ; $!(w)$ , set the width of the branch to $w$ ; $[$ , each $[$ produces a new $\boldsymbol {P_j^{i+1}}$ , all actions after $[$ are performed with respect to the $\boldsymbol {P^i}$ state; $]$ , signifies completion of operations on $\boldsymbol {P_j^{i+1}}$ .

Table 1. Parameters for fractal trees shown in figures 2 and 3 using (2.1) (Prusinkiewicz et al. 1996). A minimum value is set to avoid unresolvable branches.

2.4.2. Fractal tree geometry and fractal dimension

In this study, using the parameters stated in table 1, three configurations of the same type of fractal tree geometry were generated using the parametric L-system. Then, all three configurations of trees heights ( $H$ ) were normalized to a length of 1 m while keeping the horizontal and vertical scales the same (figures 2 and 3). The fractal iteration numbers are $n\,=\,4$ , 6 and 8. The values of $L_0$ and $D_0$ , which are the length and diameter of the thickest branch, respectively, change after normalization. Table 2 summarizes the various tree geometries used in this study. As $r_1 = r_2 = 0.8$ , the lengths $L_1, L_2$ and $L_3$ , etc. shorter than $L_0$ are $L_k = r_1L_{k-1} = 0.8L_{k-1}$ . Moreover, as $q = 1-q = 0.5$ , diameters $D_1, D_2$ and $D_3$ etc. are $D_k = q^eD_{k-1} = 0.5^{0.5}D_{k-1}$ . Here, $k$ refers to the current branching generation and $k-1$ refers to the branching generation one branch earlier.

Figure 2. Side view of the fractal tree geometry for $(a)$ Basic $n\,=\,\rm 4$ , $(b)$ Basic $n\,=\,\rm 6$ and $(c)$ Basic $n\,=\,\rm 8$ . The fractal iteration parameter $n$ is the number of branch shapes repeated with different scales. Axis conventions are the same as those used in the simulation.

Figure 3. Upside view of the fractal tree geometry for $(a)$ Basic $n\,=\,\rm 4$ , $(b)$ Basic $n\,=\,\rm 6$ and $(c)$ Basic $n\,=\,\rm8 $ .

Table 2. Fractal trees’ geometry details.

A numerical box-counting method (Da Silva et al. 2006) was utilized to estimate the fractal dimension of the fractal tree models. The tree geometry was first discretized from polygons to points along the tree surface; for the whole tree, there are points of the order of $2.5 \times 10^6$ . The domain was subdivided into boxes with a side length $\delta$ and the number of boxes which include at least one point was counted. The relationship between $H / \delta$ and count number $N$ is illustrated in figure 4. In the range before the boxes became proportional to the smaller branch dimension, all three configurations of tree models exhibited a slope of approximately 1.77. Therefore, the fractal dimension of the tree models was considered to be 1.77. When the side length of boxes became of the order of the smallest branches, the method might be used to estimate the fractal dimension of the surface (leaf, branch, etc.) but not the whole tree (Da Silva et al. 2006). Therefore, approximately when $H / \delta \geqslant 1.0 \times 10^3$ , the counts increase at a lower rate and exhibit a gentler slope. In this study, we name the present type of fractal tree, with a fractal dimension of 1.77, the ‘Basic’ tree. We employ this Basic tree with different iteration numbers for all the numerical simulations.

Figure 4. Numerical box counting dimension for Basic $n\,=\,\rm 4$ , 6 and 8.

2.5. Simulation set-up

The computational domain is $32 H \times 16 H \times 16 H$ . We defined the coordinate system for the numerical simulation as $X,Y$ and $Z$ and the coordinate system for the wake analysis as $x,y$ and $z$ , as shown in figure 5. The origin of the coordinate system, $x,y,z$ , was set as $(X,Y,Z)=(8H,8H,7.5H)$ . The tree’s centre of mass was calculated and placed at $(X,Y,Z)=(8H,8H,8H)$ . We defined the centreline as a line through the tree’s centre of mass ( $y/H\,=\,0$ , $z/H\,=\,0.5$ ) and parallel to the $x$ -axis. Uniform flow ( $U_\infty$ ) in the $x$ -direction was set as the inflow condition. Trees were treated as rigid bodies, and deformation due to fluids is not considered. We use the physical properties of air at room temperature, density $\rho =1.205 \, \mathrm {kg\,m^-{^3}}$ and kinematic viscosity $\nu =1.512 \times 10^{-5} \, \mathrm {m^2\,s^-{^1}}$ .

Figure 5. Computational domain and arrangement of trees. Here $U_\infty$ represents the uniform flow velocity in the $x$ -direction of the inflow condition. The $X,Y,Z$ coordinate system is used for numerical simulation and the $x,y,z$ coordinate system is used for wake analysis.

Figure 6. Computational mesh subdivided hierarchically around Basic $n\,=\,\rm 8$ tree.

After evaluating the resolution dependency of the drag coefficient ( $C_D$ , also see § 3.1), for Basic $n$ = 4, 6 and 8, the finest grid $\Delta x / H = 1/1024,\, 1/1024$ and $\Delta x / H = 1/2048$ is placed near the surface of the tree (approximately 40 meshes from the surface), respectively. For all three tree models, the $\Delta x / H = 1/512$ grid was set in the wake region with a length of approximately $8.0 H$ . Figure 6 represents an example of AMR grids in the case of Basic $n\,=\,\rm 8$ . The total grid numbers were 1 682 158 592, 1 711 257 088 and 2 363 039 744 in the cases of Basic $n\,=\,\rm 4$ , $n\,=\,\rm 6$ and $n\,=\,\rm 8$ , respectively. An outflow boundary condition was imposed at the computational boundary behind the tree, and inflow boundary conditions were imposed at other computational boundaries. The uniform flow condition was imposed at the inflow boundary. The computational boundaries in the $y$ - and $z$ -directions are sufficiently distant from the trees by approximately 8 $H$ , so they hardly affect the flow field around the trees. An interpolated bounce-back condition (non-slip condition) was imposed on the tree surface. The calculations were performed on the TSUBAME 4.0 supercomputer at Tokyo Institute of Technology. For Basic $n\,=\,\rm 4$ and $n\,=\,\rm 6$ , 12 NVIDIA H100 GPUs were utilized. Calculation of a non-dimensional time of 60 for $n\,=\,\rm 4$ and $n\,=\,\rm 6$ took approximately 48 and 54 hr, respectively. In the case of Basic $n$ = 8, 16 NVIDIA H100 GPUs were utilized, which completed the computations within 72 hr to calculate a non-dimensional time of 60. Here, the non-dimensional time was defined as $U_\infty t /H$ where $t$ represents the physical time. As indicated in table 3, calculations were carried out for four cases with different tree-height-based Reynolds numbers $Re_H$ .

Table 3. Uniform flow velocity and $Re_H$ for the four simulation cases. Here $Re_H$ is the Reynolds number based on the tree height, i.e. $Re_H = U_\infty H /\nu$ .

3. Verification calculation

To verify the computation results, the drag coefficient $(C_D)$ was calculated to validate the resolution dependency around the boundary layer. The Taylor microscale $\lambda$ and Kolmogorov scale $\eta$ were calculated to confirm that the simulation was executed on quasi-DNS.

3.1. Resolution dependency of drag coefficient

For Basic $n\,=\,\rm 4$ , 6 and 8, drag coefficients were calculated using several finest mesh sizes (the mesh size near the boundary layer) at $Re_H = 120\,000$ to assess the resolution dependency of the drag coefficient. The drag coefficient is calculated as follows:

(3.1) \begin{align} \begin{split} C_D = \frac {2F}{\rho U_\infty ^2A}. \end{split} \end{align}

Here, $F$ is the drag force in the $x$ -direction acting on the tree and $A$ is the projected area of the tree in the $x$ -direction. The fluid force acting on the tree’s surface is calculated based on the momentum exchange (Wen et al. 2014) between the velocity distribution function and the object at the object boundary, as described in § 2.3. At each grid point, the momentum change of the velocity distribution function that bounced back at the object boundary was calculated and then integrated over the object surface to obtain the forces acting on the object. The $x$ -component of this force is extracted to obtain the drag force, $F$ , in the $x$ -direction. The temporal variation of the drag coefficient is illustrated in figure 7. In the resolution dependency study, only the finest mesh size differed when investigating resolution dependence. Note also that since changing the wake mesh size had little effect on the drag coefficient results, the wake mesh size of $\Delta x/H = 1 /256$ was used instead of $\Delta x/H = 1 /512$ for this section only. The error between these two types of wake mesh sizes in the $C_D$ values was found to be less than 0.5 $\%$ . We did not investigate the resolution dependence of meshes in the wake region because the mesh size of $\Delta x/H = 1 /512$ was the highest resolution we could achieve and improving the wake region mesh size to $\Delta x/H = 1 /1024$ , i.e. almost comparable to that of the Kolmogorov scale, would have resulted in the total number of meshes approaching $10 \times 10^{9}$ , which would incur unrealistic computational costs.

As observed in figure 7, in the flow-developing interval (non-dimensional time $U_\infty t /H \lt 7.5$ ), the flow was affected by the initial conditions and the drag coefficients were not stable. When $U_\infty t /H \geqslant 7.5$ , the drag coefficients remained approximately constant, indicating that the flow around the tree was fully developed. Therefore, the time average of the drag coefficient was calculated from $U_\infty t /H = 7.5$ to the end of the simulation. Notably, the calculation utilized the finest mesh size of $\Delta x /H = 1/4096$ for $n\,=\,8$ , which was only calculated up to approximately $U_\infty t /H = 9.0$ due to the considerable computational cost.

Figure 7. Temporal variation of the drag coefficient of trees at various finest mesh sizes for flows with $Re_H = 120\,000$ for $(a)$ Basic $n\,=\,\rm 4$ , $(b)$ Basic $n\,=\,\rm 6$ and $(c)$ Basic $n\,=\,\rm 8$ .

Table 4 summarized the resolution dependence of the drag coefficients for Basic $n\,=\,\rm 4$ , 6 and 8. As observed in the case of $n\,=\,\rm 4$ , when the finest mesh size was increased from $\Delta x /H = 1/1024$ to $\Delta x /H = 1/2048$ , the drag coefficient $C_D$ changes from 0.682 to 0.704, a change of 3.23 $\%$ . Similarly, for $n\,=\,\rm 6$ , the finest mesh size increases from $\Delta x /H = 1/1024$ to $\Delta x /H = 1/2048$ , which is a change of 1.06 $\%$ and for $n\,=\,\rm 8$ , the finest mesh size increases from $\Delta x /H = 1/2048$ to $\Delta x /H = 1/4096$ , which is a change of 1.06 $\%$ . The variations were all approximately 1 $\%$ –3 $\%$ , indicating that the $C_D$ values converge at the finest mesh sizes of $\Delta x /H = 1/1024$ , $\Delta x /H = 1/1024$ and $\Delta x /H = 1/2048$ for $n\,=\,\rm 4$ , 6 and 8, respectively. Therefore, all simulations for $n\,=\,4$ , 6 and 8 were executed with the mesh sizes described above. In the above cases, the ratios of the finest branch diameters to the finest mesh size, $D_{min} /\Delta x$ , were 38.81, 15.69 and 13.97 for $n\,=\,\rm 4$ , 6 and 8, respectively.

Table 4. Time average drag coefficients ( $C_D$ ) at various finest mesh sizes for Basic $n\,=\,4$ , 6 and 8.

Additionally, to confirm the boundary layer of the tree can be accurately resolved, the non-dimensional distance $y^+$ at a distance of 1 mesh from the tree surface was calculated using

(3.2) \begin{align}y^+ &= \frac {u^*y}{\nu }\nonumber\\&= \sqrt { \frac { \mu \frac {\partial u}{\partial y} \left |_{y=0}\right. }{\rho }} \frac {y}{\nu }. \end{align}

Here, $u^*$ is the friction velocity, $y$ represents the distance from the tree surface (in this paragraph only), $\mu$ is the viscosity coefficient and $\partial u/\partial y$ denotes the velocity gradient perpendicular to the tree surface (in this paragraph only). In this study, $u^*$ was calculated using the tangential velocity and velocity gradient at a distance of 1 mesh from the tree surface. Near the central trunk part of a tree, before the boundary layer separation occurred, the $y^+$ is approximately 2.0–6.0, 4.5–7.2 and 1.8–4.9, respectively, for Basic $n\,=\,\rm 4$ , 6 and 8. Near the thinnest tip branches, before the boundary layer separation occurred, the $y^+$ is approximately 1.5–6.3, 4.2–8.2 and 1.6–6.3, respectively, for Basic $n\,=\,\rm 4$ , 6 and 8 at the finest mesh sizes of $\Delta x / H = 1/1024$ , $\Delta x / H = 1/1024$ and $\Delta x / H = 1/2048$ .

3.2. Energy dissipation estimation

Calculation of the gradient of velocity fluctuations and their mean-squared values is necessary to estimate important parameters such as energy dissipation, Taylor microscale, Kolmogorov scale and integral length scale.

For energy dissipation, in previous experimental studies of grid turbulence, most of the flow field was measured with particle image velocimetr or hot-wire anemometers, which provided only one- or two-dimensional velocity information. In this study, energy dissipation is calculated using the nine components of the gradient of velocity fluctuations using the following equation:

(3.3) \begin{align} \begin{split} \epsilon = 2\nu \left( \overline {s_{11}^2} + \overline {s_{22}^2} + \overline {s_{33}^2} + 2\overline {s_{12}^2} + 2\overline {s_{13}^2} + 2\overline {s_{23}^2} \right), \end{split} \end{align}

where, $s_{ij} = ( (\partial u_i / \partial x_j) + (\partial u_j / \partial x_i) ) / 2$ is the strain rate of velocity fluctuations and $u_i$ is the velocity fluctuation. The overbar represents the time average. We confirmed that the time-averaged values of TKE (TKE $\equiv \frac {1}{2} \overline {q^2}$ , where $\overline {q^2} \equiv \overline {u^2} + \overline {v^2} + \overline {w^2}$ ), averaged over a non-dimensional time length of 1, converge to a certain value around the non-dimensional time $U_\infty t / H = 10$ at several centreline wake locations ( $x /H\,=\,\rm 2.6$ , 4.6 and 6.6). To be on the safe side, the energy dissipation was calculated from non-dimensional time 20–110, 747 instantaneous velocity fields.

Figure 8. Ratio of mesh size in wake region to Kolmogorov scale in the case of Basic $n\,=\,\rm 8$ at $Re_H = 120\,000$ .

Knowing the energy dissipation, the Kolmogorov scale $\eta = (\nu ^3/\epsilon )^{1/4}$ can be computed. Figure 8 illustrates the ratio of mesh size in the wake region to the Kolmogorov scale along the centreline in the case of Basic $n\,=\,\rm 8$ at $Re_H = 120\,000$ . The mesh size in the wake region was $\Delta x / H\,=\,1/512$ , and the average value of $\eta / \Delta x$ from $x /H\,=\,\rm 0.34\;{\rm to}\;8.0$ is approximately 0.498. Where the wake was well developed, the mesh size was approximately twice that of the Kolmogorov scale. The Taylor microscale $\lambda = \sqrt { \overline {u^2} / \overline {(\partial u / \partial x)^2} }$ was also calculated. Here $u$ is the velocity fluctuation in the streamwise direction ( $x$ -direction). The average value of $\lambda / \Delta x$ from $x /H\,=\,\rm 0.34{-}8.0$ was approximately 10.317. The wake mesh size was significantly smaller than the Taylor microscale and approximately twice that of the Kolmogorov scale when $x /H \geqslant 3.0$ ; the simulation can be considered a quasi-DNS. This was also the highest resolution we could achieve, and in any case, it is very difficult to achieve calculations with a wake mesh size fully comparable to the Kolmogorov scale.

The characteristic time such as the Kolmogorov time scale $\tau _\eta = \sqrt {\nu /\epsilon }$ and Taylor time scale $\tau _\lambda = \lambda / u^{\prime }$ are also considered. Using the average value of $\eta$ and $\lambda$ mentioned above, $\tau _\lambda /T_0$ and $\tau _\eta /T_0$ were calculated to be approximately 0.477 and 0.122 at $Re_H\,=\,\rm 120\,000$ , respectively. Here, the representative time $T_0=H/U_\infty$ . In our simulation, at $Re_H\,=\,\rm 120\,000$ , $\Delta t_{wake}/T_0$ was approximately 0.00187 for Basic $n\,=\,\rm 4$ and 6 and 0.000938 for Basic $n\,=\,\rm 8$ , where $\Delta t_{wake}$ represents the time step of the mesh in the wake region. They are considerably smaller than $\tau _\lambda /T_0$ and $\tau _\eta /T_0$ and reasonably tracked the velocity fluctuations.

4. Results

4.1. Comparison of drag coefficient with the literature

Figure 9 illustrates the relationship between $Re_H$ and $C_D$ . Grant & Nickling (1998) designed a direct field measurement and demonstrated a decrease in the $C_D$ values with increasing the Reynolds number at $Re_H = 15\,000$ and 25 000. These are medium Reynolds numbers, and the simulations in this study qualitatively reproduced the observed decrease in $C_D$ values with increasing $Re_H$ . Manickathan et al. (2018) demonstrated that the $C_D$ values remain almost constant when $Re_H \geqslant 60\,000$ for both model and natural trees using wind tunnel measurements. The large spread of points in the graph obtained from previous studies can be attributed to the fidelity of the tree models. For example, the older field measurement study by Grant & Nickling (1998) used artificial Scots pine Christmas trees; Manickathan et al. (2018) used both model trees and natural hardwood and coniferous trees; and Gillies et al. (2002) used real trees with branch trimming. In the case of Gillies et al. (2002), the overall small $C_D$ values they obtained can be considered due to their use of more flexible specimens and the contribution of wall turbulence from the wind tunnel experiments, which contributed to the lower drag estimates in their study. Specifically, wall turbulence may mix with the tree wake, enhancing wake recovery and reducing drag. Gillies et al. (2002) placed the tree on the ground, we speculate that their $C_D$ results were smaller than the present ones due to the influence of the ground. For trees with branches, such as Burning Bush and Colorado Spruce, the $C_D$ values remain approximately constant when $Re_H \geqslant 200\,000$ . Fountain Grass does not exhibit this tendency, likely due to a strong sway in the wind, as it lacks rigid bodies such as trunks and branches.

In this study, all three tree models showed a tendency for $C_D$ to decrease with increasing Reynolds number at low and medium Reynolds numbers ( $Re_H\,=\,2500{-}10\,000$ ). However, after $Re_H$ exceeded 60 000, $C_D$ remained almost constant. These are in qualitative agreement with the study by Grant & Nickling (1998) and Manickathan et al. (2018), despite our tree model not including leaves. Meanwhile, the comparison of $C_D$ between different tree models at the same $Re_H$ revealed that, at low and medium Reynolds numbers ( $Re_H\,=\,\rm 2500{-}10\,000$ ), $C_D$ increases with fractal iteration number ( $n$ ), with a more pronounced difference observed. At high Reynolds numbers ( $Re_H \geqslant\,\rm 60\,000$ ), although $n\,=\,\rm 4$ and $n\,=\,\rm 6$ still exhibited some $C_D$ differences, the difference between $n\,=\,\rm 6$ and $n\,=\,\rm 8$ is minimal. We speculate that, under high Reynolds number conditions (e.g. $Re_H\,=\,\rm 120\,000$ ), if the fractal iteration number ( $n$ ) reaches a certain threshold (e.g. $n\,=\,\rm 6$ ), then $C_D$ will be independent of $Re_H$ and the tree shape and it can be considered a constant value.

Figure 9. Relationship between $C_D$ and $Re_H$ obtained in this study, and comparison with previous studies.

Figure 10 depicts the isosurfaces of the second invariant of the velocity gradient tensor ( $Q$ value) of the flow around Basic $n\,=\,\rm 8$ at $Re_H = 2500$ and $120\,000$ . The $Q$ value was calculated using the following equation with Einstein’s summation notation:

(4.1) \begin{align} \begin{split} Q = \frac {1}{2} (W_{ij} W_{ij} - S_{ij}S_{ij}). \end{split} \end{align}

Here, $W_{ij}$ represents an asymmetric component of the velocity gradient tensor and $S_{ij}$ represents a symmetric component of the velocity gradient tensor. To clearly visualize the vortex structures and fairly compare the $Q$ isosurfaces for $n\,=\,\rm 8$ between different $Re_H$ , we set the threshold ( $Q_{threshold}$ ) to 0.01 $\%$ of the maximum $\overline {Q}$ in the region $0 \lt x/H \lt 1.5$ , $-0.5 \lt y/H \lt 0.5$ , $-0.062 \lt z/H \lt 1.156$ . This region was selected because it exhibits significantly higher centreline turbulence intensity than that in others, as shown in $\S$ 4.2.1. It can be speculated that $\overline {Q}$ in this region is also significantly larger than in others.

Figure 10. Isosurfaces of the second invariant of the velocity gradient tensor of the flow around Basic $n\,=\,\rm 8$ at $Re_H = 2500$ and $Re_H = 120\,000$ . The wake region’s opacity is adjusted to improve the visibility of the tree’s branching part.

Figure 11. Comparison of streamwise evolution of centreline turbulence intensity normalized by its maximum value is given as a function of $x$ scaled by $x_{peak}$ at $(a)$ $Re_H=2500$ , $(b)$ $Re_H=10\,000$ , $(c)$ $Re_H=60\,000$ and $(d)$ $Re_H=120\,000$ .

Specifically, $Q_{threshold} T_0^2$ was set to 37.4 for $Re_H\,=\,\rm 2500$ and 111.8 for $Re_H\,=\,\rm 120\,000$ . Note that the representative time $T_0$ varies with $Re_H$ . Only the grid points with $Q$ values equal to $Q_{threshold}$ were used to generate the isosurface. As shown in figure 10, for $Re_H\,=\,\rm 2500$ , elongated streak structures were concentrated in the wake region behind the centre and lower parts of the tree, with almost no vortices appearing behind the canopy. It can be inferred that the local Reynolds number near the canopy was too low for the eddies generated behind it to exceed the threshold. Additionally, the wake length for $Re_H\,=\,\rm 2500$ is shorter compared with that of $Re_H\,=\,\rm 120\,000$ , primarily because $Re_H\,=\,\rm 2500$ is too small for its wake to be considered turbulent (discussed in § 4.3.1). Meanwhile, for $Re_H\,=\,\rm 120\,000$ , numerous small eddies appeared in the wake region, indicating a highly turbulent field. This confirms the validity of the analysis using the fractal tree model generated by the L-system in this study.

Although we have only computed up to $Re_H=120\,000$ with a maximum iteration number of $n=8$ , we can give a conjecture about what might occur as $Re_H$ and $n$ increase further. If $n$ remains constant, the turbulent scales in the wake region, such as $\eta$ , may decrease as $Re_H$ further increases. However, if both $n$ and $Re_H$ increase, each scenario must be analysed individually as they are not interchangeable. This is because an increase in $Re_H$ increases the local Reynolds numbers, while an increase in $n$ decreases them, causing the opposite effect.

4.2. Turbulence intensity

4.2.1. Centreline turbulence intensity

The centreline turbulence intensities for different tree models at the same Reynolds number as well as different Reynolds numbers for the same tree model were compared. The turbulence intensity was defined as follows:

(4.2) \begin{align} \begin{split} \frac { u_{avg}^{\prime} } { U_\infty } = \frac { \sqrt {\frac {1}{3} \overline {u_i^{2}} } } { U_\infty }, \end{split} \end{align}

where $u_{avg}^{\prime}$ is calculated based on $\overline {u_i^{2}}$ and $u_i$ denotes the velocity fluctuation in the $i$ direction.

In a study on fractal-grid-generated turbulence, Gomes-Fernandes, Ganapathisubramani & Vassilicos (2012) used a method to normalize the centreline turbulence intensity. This method has also been used in the present study, where normalization is performed by dividing the turbulence intensity by its peak value and then dividing $x$ by the $x$ -coordinate where the magnitude of turbulence intensity is maximum ( $x_{peak}$ ). Figure 11 illustrates the results using this normalization. As observed in the figure, the turbulence intensity collapses more effectively when the Reynolds number is high. At low Reynolds numbers (figure 11 a), it does not collapse if $x /x_{peak} \gt\,\rm 1.0$ . The role of viscosity was considered strong, and the flow field was not yet sufficiently developed for turbulence. Hence, this normalization method is unsuitable for cases with low Reynolds numbers. At intermediate Reynolds numbers (figure 11 b), the collapse improves. A more complete collapse appeared at higher Reynolds numbers (figure 11 c,d). As observed in the figure, except for the case of Basic $n\,=\,4$ , the turbulence intensities for $n\,=\,6$ and $n\,=\,8$ collapsed effectively. However, this collapse has only been confirmed for the centreline turbulence intensity of $n\,=\,6$ and 8. Further research must be conducted to determine whether it applies to $n\,=\,9$ , 10 or higher.

Figure 12 illustrates the dependence of turbulence intensity on the Reynolds number for the same tree model. In all three models, $x_{peak}$ moves upstream as the Reynolds number increases. This trend converges at $Re_H =10\,000$ when $n\,=\,\rm 4$ and 8, as there is virtually no difference in the turbulence intensity beyond $Re_H = 10\,000$ . For $n\,=\,\rm 6$ , it converges at $Re_H =60\,000$ . For $n\,=\,\rm 4$ and 8, the peak turbulence intensity increases with an increasing Reynolds number; however, it stagnates at $Re_H =10\,000$ . This suggests that the centreline turbulence intensity hardly changes after a critical Reynolds number (e.g. $Re_H =10\,000$ ). Similarly, for $n\,=\,\rm 6$ , the peak turbulence intensity also grows with increasing Reynolds number and stagnates at $Re_H =60\,000$ . At $Re_H\,=\,\rm 60\,000$ and 120 000, the peak turbulence intensity is significantly greater at $n\,=\,\rm 6$ than at $n\,=\,\rm 4$ and 8. This discrepancy can be attributed to the fact that the peak turbulence intensity does not always occur at $y /H\,=\,\rm 0$ and $z /H\,=\,\rm 0.5$ (the centreline) when viewed along the $y$ - and $z$ -axes. For instance, as shown in figure 13, the peak value of turbulence intensity does not always occur at $y/H\,=\,\rm 0$ , and similarly, the peak value of spatially averaged turbulence intensity does not always occur at $z/H\,=\,\rm 0.5$ , as shown in figure 14 (discussed later).

In the study by Gomes-Fernandes, Ganapathisubramani & Vassilicos (2012), they defined an equation that predicts the peak centreline turbulence intensity and $x_{peak}$ position in the case of grid turbulence. For fractal tree turbulence, an equation for predicting the peak turbulence intensity and location of $x_{peak}$ may exist, which needs to be further investigated.

Figure 12. Comparison of streamwise evolution of centreline turbulence intensity at various $Re_H$ in the case of $(a)$ Basic $n\,=\,4$ , $(b)$ Basic $n\,=\,6$ and $(c)$ Basic $n\,=\,8$ .

4.2.2. Height dependency of turbulence intensity

The height dependency of turbulence intensity for different tree models is presented. As discussed in § 4.2.1, in the case of the same tree model, the turbulence intensity exhibits almost the same behaviour beyond $Re_H = 60\,000$ . Hence, only the case of $Re_H = 120\,000$ is explored in this subsection.

Figure 13 illustrates an example of the turbulence intensity distribution in the spanwise direction for Basic $n\,=\,8$ for $z /H = 0.5$ at various downstream distances. As observed in the figure, the peak value of turbulence intensity does not appear on the centreline ( $y=0$ ). The asymmetric profiles can be attributed to the asymmetric tree shape, as shown in figure 3. The tree branches are located closer to the inflow direction (smaller $x$ coordinate) on the $-y$ side than on the $+y$ side. Simultaneously, the centreline turbulence intensity first grows rapidly to a peak value and then decays gradually along the $x$ -direction, as shown in § 4.2.1. It is reasonable to assume that this is also applicable for the $y \neq 0$ cases. Thus, it can be surmised that the turbulence intensity on the $-y$ side reaches a peak value and starts to decay at a more forward position than that on the $+y$ side, resulting in asymmetric profiles depicted in figure 13. We also considered that the evaluation of turbulence intensity in the canopy part is insufficient only for the centreline, as the thin branches also spread in the spanwise direction. For a given height (e.g. $z/H\,=\,0.5$ , as shown in figure 13), we plotted the transverse profiles of turbulence intensity at various $x$ coordinates. At each $x$ coordinate, we identified the peak turbulence intensity on the line along the $y$ direction. We defined the half-width of turbulence intensity as the range where the turbulence intensity exceeds half of this peak on that line. At each $x$ coordinate, turbulence statistics, including turbulence intensity, the Taylor-microscale-based Reynolds number and isotropy parameters, were averaged over $y$ within the half-width of turbulence intensity. We repeated these steps for different heights to obtain the height dependency of various turbulence statistics.

Figure 13. Turbulence intensity distribution in the spanwise direction for Basic $n\,=\,8$ at $z /H=0.5, Re_H = 120\,000$ at various downstream distances.

Figure 14 illustrates the height dependency of turbulence intensity. Similar to the centreline’s results, it is shown to decrease gradually in the downstream direction in all three models. Near the wake locations behind the fractal tree ( $x /H = 0.6$ ), turbulence intensity is markedly more potent at $z /H\,=\,0.4-0.65$ for $n\,=\,4$ , $z /H\,=\,0.3-0.55$ for $n\,=\,6$ and $z /H\,=\,0.3-0.5$ for $n\,=\,8$ , compared with other heights. This suggests that the turbulence intensity is stronger near the tree trunk and second branch height than at other altitudes. We considered that the large-scale turbulence generated by the trunk is dominant in near-wake locations. This trend is maintained after $x /H = 0.6$ , although the overall turbulence intensity gradually decreases. Generally, at sufficiently distant downstream locations ( $x /H \geqslant 2.2$ ), the height at which turbulence intensity is the strongest is approximately $z /H\,=\,0.35{-}0.5$ . The difference in turbulence intensity between the central part of the trunk and the canopy part also becomes smaller. This may be due to the fractal nature of the tree geometry, which promoted vertical mixing within the wake, as reported in a previous study (Schröttle & Dörnbrack 2013).

Figure 14. Turbulence intensity as a function of height ( $z /H$ ) at $Re_H = 120\,000$ : $(a)$ Basic $n\,=\,4$ ; $(b)$ Basic $n\,=\,6$ ; $(c)$ Basic $n\,=\,8$ . The angle bracket, $\langle \cdot \rangle _{\text {half-width}}$ , represents the spatial average over $y$ within the half-width of turbulence intensity, which defined in § 4.2.2.

4.2.3. Transverse plane average of turbulence intensity

Although we have analysed the centreline turbulence intensity and its height dependency, the complex geometry of the tree presents a highly 3-D flow field. Similar to the method employed by Laizet & Vassilicos (2015), we considered the average turbulence intensity within the $y$ – $z$ plane for further analysis. The $y$ – $z$ plane was selected as a rectangular area that closely follows the tree’s outer boundary, capturing its $x$ -direction projection while minimizing the rectangle’s area, as shown in figure 15 $(a)$ . This selection ensures that the plane captures the tree’s influence on the flow field with the smallest cross-sectional area. The turbulence intensity is averaged over $y$ and $z$ within this plane. Figure 15 $(b)$ shows the streamwise evolution of spatially averaged turbulence intensity.

Figure 15. $(a)$ The $y$ – $z$ plane used to obtain the spatial average. $(b)$ Streamwise evolution of turbulence intensity (with angle brackets $\langle \cdot \rangle _{yz}$ denoting averaging over $y$ and $z$ ) for Basic $n\,=\,4$ , 6 and 8 at $Re_H = 120\,000$ .

The average turbulence intensity for $n\,=\,4$ , 6 and 8 exhibits both a region of rapid growth and a region of decay, with the peak marking the boundary between these regions, as shown in figure 15 $(b)$ . The peaks are located at approximately $x /H\,=\,0.4$ , similar to the results from the centreline analysis (figure 12). Regarding the peak values, $n\,=\,4$ exhibits a higher peak value of 10.2 $\%$ when compared with 7.0 $\%$ and 6.9 $\%$ for $n\,=\,6$ and $n\,=\,8$ , respectively. This difference might be attributed to the dominance of the large-scale turbulence in the wake region of $n\,=\,4$ , which tends to exhibit higher turbulence intensity than the small-scale turbulence. For $n\,=\,6$ and 8, the higher fractal iteration number create thinner and denser branches than $n\,=\,4$ , which might lead to the dominance of small-scale turbulence in their wake regions. Furthermore, figure 14 shows that the turbulence intensity in the canopy region of $n\,=\,4$ ( $z /H\,=\,0.75-0.95$ ) is significantly higher, approximately 10 $\%$ –13 $\%$ , when compared with 8 $\%$ –11 $\%$ for $n\,=\,6$ ( $z/ H\,=\,0.65-0.95$ ) and 6 $\%$ –10 $\%$ for $n\,=\,8$ ( $z /H\,=\,0.6-0.95$ ). Additionally, since all the tree heights were scaled to 1 m, the trunks of $n\,=\,6$ and $n\,=\,8$ are thinner than those of $n$ = 4. The turbulence intensity near the trunks of $n$ = 4 ( $z /H\,=\,0.35-0.65$ ) is approximately 18 $\%$ –20 $\%$ , while it is approximately 16 $\%$ –17 $\%$ for $n\,=\,6$ ( $z /H\,=\,0.3-0.55$ ) and 13 $\%$ –16 $\%$ for $n\,=\,8$ ( $z /H\,=\,0.3{-}0.5$ ). These factors might contribute to the higher peak value observed in $n\,=\,4$ .

In terms of the decay rate, $n\,=\,8$ exhibits a noticeably faster decay in the near wake region (approximately $x /H \leqslant 2.0$ ). For example, at $x /H\,=\,1.0$ , the turbulence intensity for $n\,=\,8$ drops to 3.8 $\%$ , which is 54.8 $\%$ of its peak value, while for $n\,=\,4$ and $n\,=\,6$ , the values are 7.3 $\%$ and 4.7 $\%$ , corresponding to 71.5 $\%$ and 66.4 $\%$ of their respective peak values. As mentioned earlier, in the wake region of $n\,=\,8$ , small-scale turbulence might be dominate. The small-scale turbulence tends to dissipate faster than the large-scale one, which might contribute to the rapid decay of the average turbulence intensity in the near-wake region (approximately $x /H \leqslant 2.0$ ) of $n\,=\,8$ .

Figure 16. Comparison of streamwise evolution of centreline Taylor-microscale-based Reynolds number as a function of $x /x_{peak}$ at $(a)$ $Re_H=2500$ , $(b)$ $Re_H=10\,000$ , $(c)$ $Re_H=60\,000$ and $(d)$ $Re_H=120\,000$ .

4.3. Turbulent Reynolds number

4.3.1. Centreline turbulent Reynolds number

Figure 11(c,d) illustrate an effective collapse in the turbulence intensity between Basic $n\,=\,6$ and $n\,=\,8$ . Therefore, we employed $x_{peak}$ to normalize streamwise distances in all subsequent centreline comparisons between tree models. Figure 16 depicts the centreline Taylor-microscale-based Reynolds number, $Re_\lambda = u^{\prime} \lambda / \nu$ . At low Reynolds number (figure 16 a), in all three models, the $Re_\lambda$ are all small, with peak values below 30 and cannot be considered as even turbulent flow. At medium Reynolds number (figure 16 b), $Re_\lambda$ becomes larger overall.

At high Reynolds number (figure 16 c,d), $Re_\lambda$ increases further, with $n$ = 6 showing the largest $Re_\lambda$ in most decaying regions. The data do not collapse properly, so it can be considered that $x_{peak}$ is not suitable for collapsing centreline $Re_\lambda$ . At $Re_H\,=\,60\,000$ and 120 000, the peak $Re_\lambda$ is significantly greater at $n\,=\,6$ than at $n\,=\,4$ and 8. This discrepancy likely arises because the peak $Re_\lambda$ does not always occur at $y /H\,=\,0$ and $z /H\,=\,0.5$ (the centreline) when viewed along the $y$ - and $z$ -axes.

As shown in figure 16, $Re_\lambda$ remains approximately constant following a certain downstream location, except for the case of $n\,=\,6$ at $Re_H\,=\,2500$ . Specifically, in figure 16(a,b), $Re_\lambda$ remains approximately constant at approximately $x \geqslant\,3$ $x_{peak}$ and 6 $x_{peak}$ for $n\,=\,4$ and 8 at $Re_H=2500$ and at approximately $x \geqslant\,6$ $x_{peak}$ for $n\,=\,4$ , 6 and 8 at $Re_H=10\,000$ . Similarly, in figure 16(c,d), $Re_\lambda$ remains approximately constant at approximately $x \geqslant$ 6 $x_{peak}$ , 12 $x_{peak}$ and 10 $x_{peak}$ for $n\,=\,4$ , 6 and 8 at $Re_H\,=\,60\,000$ and at approximately $x \geqslant$ 7 $x_{peak}$ , 12 $x_{peak}$ and 12 $x_{peak}$ for $n\,=\,4$ , 6 and 8 at $Re_H\,=\,120\,000$ . This phenomenon may be attributed to the gradual decrease in $u^{\prime}$ and the gradual increase in $\lambda$ in the downstream direction, causing $Re_\lambda$ to remain approximately constant.

Figure 17 depicts the dependence of $Re_\lambda$ on the Reynolds number for the same tree model. In all three models, $Re_\lambda$ increases overall with increasing Reynolds number, similar with the turbulence intensity trend illustrated in figure 12. Notably, in figure 12, $u_{avg}^{\prime}$ is normalized by uniform flow velocity $U_\infty$ , so $u_{avg}^{\prime}$ or $u^{\prime}$ will increase with increasing Reynolds number. The position at which $Re_\lambda$ is maximum always appears at approximately $x /H\,=\,0.4{-}0.5$ , independent of the Reynolds number.

Figure 17. Comparison of streamwise evolution of centreline Taylor-microscale-based Reynolds number as a function of $x /H$ at various $Re_H$ in the case of $(a)$ Basic $n\,=\,4$ , $(b)$ Basic $n\,=\,6$ and $(c)$ Basic $n\,=\,8$ .

4.3.2. Height dependency of turbulent Reynolds number

Similar to the height dependency of turbulence intensity depicted in § 4.2.2, the height dependency of $Re_\lambda$ is also given in figure 18. The spatial average, $\langle Re_\lambda \rangle _{\text {half-width}}$ , was taken over $y$ within the half-width of turbulence intensity, which is defined in § 4.2.2.

Figure 18. Taylor-microscale-based Reynolds number as a function of height ( $z /H$ ) at $Re_H = 120\,000$ : $(a)$ Basic $n\,=\,4$ ; $(b)$ Basic $n\,=\,6$ ; $(c)$ Basic $n\,=\,8$ . The angle bracket, $\langle \cdot \rangle _{\text {half-width}}$ , represents the spatial average over $y$ within the half-width of turbulence intensity, which defined in § 4.2.2.

Here $Re_\lambda$ decreased gradually in the downstream direction in all three models. In near-wake locations behind the fractal tree ( $x /H =0.6$ ), $Re_\lambda$ is considerably more potent at $z /H\,=\,0.4{-}0.65$ for $n\,=\,4$ , $z /H\,=\,0.3{-}0.55$ for $n\,=\,6$ and $z /H\,=\,0.3{-}0.5$ for $n\,=\,8$ , compared with other altitudes. This is because the turbulence intensity is stronger near the tree trunk than at the tree crown, as mentioned before. Generally, at sufficiently distant downstream locations ( $z /H \geqslant 2.2$ ), similar with the turbulence intensity, the height at which $Re_\lambda$ is most vigorous is approximately $z /H\,=\,0.35{-}0.5$ in all three models.

4.4. Global and local isotropy

4.4.1. Centreline global and local isotropy parameters

In this section, the isotropy of the wake region of fractal trees is assessed. First, the global isotropy associated with the large scales is assessed, followed by the local isotropy related to the small scales. For clarity, we only present Basic $n\,=\,8$ data. Notably, all the other models also exhibited very similar values and the same trend.

Figure 19 compares the centreline global isotropy parameters $u^{\prime} /v^{\prime}$ , $u^{\prime} /w^{\prime}$ and $v^{\prime} /w^{\prime}$ at various Reynolds numbers for Basic $n\,=\,8$ . At low Reynolds numbers (figure 19 $a$ ), except for $x /x_{peak} =\,4.0{-}6.0$ , $u^{\prime} /v^{\prime}$ , $u^{\prime} /w^{\prime}$ and $v^{\prime} /w^{\prime}$ are significantly away from 1.0, indicating that the wake region is not globally isotropic overall. At intermediate Reynolds numbers (figure 19 b), an improvement in global isotropy is observed, with most parameters lying between 0.9 and 1.2. At high Reynolds numbers (figure 19 c,d), except for the region $x /x_{peak} \lt 1.0$ , all three global isotropy parameters generally hover around 1.0 and are marginally less than 1.2. This result is close to that of grid turbulence (Gomes-Fernandes, Ganapathisubramani & Vassilicos 2012) with $u^{\prime} /v^{\prime}$ around 1.2 when $x /x_{peak} \gt 1.0$ . For Reynolds numbers that are medium and high, the flow is globally anisotropic if $x /x_{peak} \lt 1.0$ , while it is almost globally isotropic if $x/x_{peak} \gt 1.0$ . This concurs with the grid turbulence case.

Figure 19. Global isotropy parameters of centreline for Basic $n\,=\,8$ at $(a)$ $Re_H=2500$ , $(b)$ $Re_H=10\,000$ , $(c)$ $Re_H=60\,000$ and $(d)$ $Re_H=120\,000$ .

In the study of grid turbulence, local isotropy was assessed using two relations ( $K_1, K_3$ ) derived by Taylor (1935). Namely, small-scale turbulence is locally isotropic if the Reynolds number is sufficiently high and $K_1, K_3$ should equal to 1. We also evaluated $K_1$ and $K_3$ , and additionally, we defined $K_2$ and $K_4$ , as follows, to assess the local isotropy of the fractal tree wake region:

(4.3) \begin{align} \begin{split} K_1 = 2 \, \overline {\left (\frac {\partial u}{\partial x}\right )^2} / \overline {\left (\frac {\partial v}{\partial x}\right )^2},\,\,\,\, K_3 = 2 \, \overline {\left (\frac {\partial u}{\partial x}\right )^2} / \overline {\left (\frac {\partial u}{\partial y}\right )^2}, \end{split} \end{align}

(4.4) \begin{align} \begin{split} K_2 = 2 \, \overline {\left (\frac {\partial u}{\partial x}\right )^2} / \overline {\left (\frac {\partial w}{\partial x}\right )^2},\,\,\,\, K_4 = 2 \, \overline {\left (\frac {\partial u}{\partial x}\right )^2} / \overline {\left (\frac {\partial u}{\partial z}\right )^2}. \end{split} \end{align}

Figures 20 and 21 compare the centreline local isotropy parameters $K_1, K_3$ and $K_2, K_4$ at various Reynolds numbers for Basic $n\,=\,8$ .

Similar with global isotropy parameters, at low Reynolds numbers (figures 20 a and 21 $a$ ), $K_1, K_2, K_3$ and $K_4$ are significantly away from 1.0, indicating that the wake region is not locally isotropic. At medium Reynolds numbers (figures 20 band 21 $b$ ), an improvement in the local isotropy is observed, but $K_1$ , $K_2$ and $K_4$ are still considerably away from 1.0, which does not fit the definition of a locally isotropic fluid. At high Reynolds numbers (figures 20 c,d and 21 $c, d$ ), all four parameters generally hover around 1.0, which is close to the value in a previous grid turbulence study (Gomes-Fernandes, Ganapathisubramani & Vassilicos 2012). Generally, the wake region for $x/x_{peak} \gt 1.0$ can be regarded as essentially isotropic at high Reynolds numbers ( $Re_H \geqslant 60\,000$ ).

Figure 20. Local isotropy parameters of centreline for Basic $n\,=\,8$ at $(a)$ $Re_H=2500$ , $(b)$ $Re_H=10\,000$ , $(c)$ $Re_H=60\,000$ and $(d)$ $Re_H=120\,000$ .

Figure 21. Local isotropy parameters of centreline for Basic $n\,=\,8$ at $(a)$ $Re_H=2500$ , $(b)$ $Re_H=10\,000$ , $(c)$ $Re_H=60\,000$ and $(d)$ $Re_H=120\,000$ .

Figure 22. Global isotropy parameters as a function of height ( $z /H$ ) at $Re_H = 120\,000$ for Basic $n\,=\,8$ : $(a)$ $u' /v'$ ; $(b)$ $u' /w'$ ; $(c)$ $v' /w'$ . The angle bracket, $\langle \cdot \rangle _{\text {half-width}}$ , represents the spatial average over $y$ within the half-width of turbulence intensity, which defined in § 4.2.2.

4.4.2. Height dependency of isotropy

The global isotropy parameters at $Re_H = 120\,000$ for Basic $n\,=\,8$ as a function of height ( $z /H$ ) are illustrated in figure 22. The spatial average was taken over $y$ within the half-width of turbulence intensity, which defined in § 4.2.2. As observed in the figure, the distributions of $u^{\prime} /v^{\prime}$ and $v^{\prime} /w^{\prime}$ ranged from approximately 0.8–1.2 and 1.0–1.2, respectively. After $x /H = 2.2$ , no obvious change with altitude could be observed. Here $u^{\prime} /w^{\prime}$ does not exhibit a clear tendency to change with height in the near wake region. After $x /H = 2.2$ , it showed a tendency to decrease from large to small from the crown to the trunk. We confirmed that this trend is also broadly consistent, but not as pronounced, for $n\,=\,6$ and almost absent for $n\,=\,4$ . Overall, the distribution of all three parameters ranged from approximately 0.8–1.2, suggesting that the wake region was essentially global isotropic at all different heights. We also confirmed this trend for Basic $n\,=\,4$ and 6.

The local isotropy parameters at $Re_H = 120\,000$ for Basic $n\,=\,8$ as a function of height ( $z /H$ ) are shown in figures 23 and 24. In these figures, $K_1$ and $K_2$ exhibit approximately the same values at different heights. However, $K_3$ and $K_4$ tended to increase from the crown to the trunk. This trend is also confirmed for $n\,=\,6$ but almost absent for $n\,=\,4$ . This suggests that the shear deformation of streamwise velocity fluctuations ( $\partial u / \partial y$ , $\partial u / \partial z$ ) is stronger near the canopy than the spatial gradient of streamwise velocity fluctuations ( $\partial u / \partial x$ ). Overall, the four parameters ranged from around 0.8–1.2, suggesting that the wake region was essentially local isotropic at all different heights. We also confirmed this for Basic $n\,=\,4$ and 6.

Figure 23. Local isotropy parameters as a function of height ( $z /H$ ) at $Re_H = 120\,000$ for Basic $n\,=\,8$ : $(a)$ $K_1$ , $(b)$ $K_3$ . The angle bracket, $\langle \cdot \rangle _{\text {half-width}}$ , represents the spatial average over $y$ within the half-width of turbulence intensity, which defined in § 4.2.2.

Figure 24. Local isotropy parameters as a function of height ( $z /H$ ) at $Re_H = 120\,000$ for Basic $n\,=\,8$ : $(a)$ $K_2$ , $(b)$ $K_4$ . The angle bracket, $\langle \cdot \rangle _{\text {half\hbox{-}width}}$ , represents the spatial average over $y$ within the half-width of turbulence intensity, which defined in § 4.2.2.

4.5. Dissipation and non-equilibrium nature

This section examines dissipation and non-equilibrium characteristics in the decay region. Specifically, the integral length scale $L$ and non-dimensional dissipation rate $C_\epsilon$ are investigated. The ‘cornerstone dissipation scaling of turbulence theory’, $\epsilon = C_\epsilon u^{\prime3}/{L}$ can be derived to $L / \lambda \propto Re_\lambda$ if the flow field is essentially isotropic and $C_\epsilon$ is a constant, as explained in the introduction. The integral length scale must be estimated to verify if this study follows the relation $L / \lambda \propto Re_\lambda$ .

4.5.1. Integral length scale

The integral length scale in the $x$ -direction, $L_u$ , was calculated by integrating the longitudinal correlation function $f(r, x)$ from $r$ = 0 to the first zero-crossing point of $r$ . Here, $f(r,x) = \overline {u(x)u(x+r)} / \overline {u(x)^2}$ , where $u(x)$ represents the velocity fluctuation in streamwise direction at $x$ .

As mentioned in § 4.4, when the tree-height-based Reynolds number $Re_H \leqslant 10\,000$ , the wake region was not isotropic, and $L / \lambda \propto Re_\lambda$ was not applicable. Hence, this discussion does not include cases with low Reynolds numbers. On the other hand, the wake region can be considered essentially isotropic if $Re_H \geqslant 60\,000$ . In this case, $L / \lambda \propto Re_\lambda$ is applicable and can be verified by plotting the dependency of $L_u /\lambda$ on $Re_\lambda$ . If $L_u / \lambda \propto Re_\lambda$ holds true, then $\epsilon = C_\epsilon u^{\prime3}/{L}$ also holds true with a constant $C_\epsilon$ . Conversely, $C_\epsilon$ is not a constant value but varies with position in the wake region.

Figure 25 shows the variation of $L_u /\lambda$ in the streamwise distance as a function of $x /x_{peak}$ . As observed in the figure, in the decaying region of $x \geqslant x_{peak}$ where $Re_\lambda$ decreases, $L_u /\lambda$ does not significantly decrease for any of the three models. At $Re_H=60\,000$ , it is approximately 3.5 for all three models. This trend is sustained at $Re_H = 120\,000$ . For $n\,=\,4$ , $L_u / \lambda$ remains approximately 4.0 until $x / x_{peak} = 5.0$ . Then, for an unknown reason, it marginally increases, and $L_u / \lambda$ is approximately 6.0. For $n\,=\,6$ , it increases gradually within the 2.5–4.5 range and then becomes approximately 4.5. For $n\,=\,8$ , it approximately fluctuates in the 3.0–4.5 range.

Figure 25. Integral length scale to Taylor microscale ratio $L_u /\lambda$ on the centreline for different tree models in relation to $x /x_{peak}$ at $(a)$ $Re_H=60\,000$ and $(b)$ $Re_H=120\,000$ .

If $Re_\lambda$ and $C_\epsilon$ are constant, it would not be important that $L_u /\lambda$ remains approximately the same value. However, as illustrated in figure 16, at high Reynolds numbers ( $Re_H \geqslant 60\,000$ ), $Re_\lambda$ continues to decrease over a long region beyond $x_{peak}$ . In particular, $Re_\lambda$ keeps decreasing approximately until $x \approx$ 6 $x_{peak}$ , 12 $x_{peak}$ and 10 $x_{peak}$ for $n\,=\,4$ , 6 and 8 at $Re_H\,=\,60\,000$ and until $x \approx$ 7 $x_{peak}$ , 12 $x_{peak}$ , 12 $x_{peak}$ for $n\,=\,4$ , 6 and 8 at $Re_H\,=\,120\,000$ . In the above decay regions of $Re_\lambda$ , we plot the variation of $L_u /\lambda$ as a function of $Re_\lambda$ , as shown in figure 26. Here $L_u /\lambda$ remains approximately constant and does not increase as $Re_\lambda$ increases. This result qualitatively concurs with the findings of fractal-grid-generated turbulence (Seoud & Vassilicos 2007; Mazellier & Vassilicos 2010; Valente & Vassilicos 2011b ; Gomes-Fernandes, Ganapathisubramani & Vassilicos 2012). Therefore, $\epsilon = C_\epsilon u^{\prime3}/{L}$ does not hold, primarily because $C_\epsilon$ is not a constant in the decaying region of the fractal-tree-generated turbulent flows studied here.

Figure 26. Integral length scale to Taylor microscale ratio $L_u /\lambda$ on the centreline for different tree models in relation to $Re_\lambda$ in the decay region of $Re_\lambda$ at $(a)$ $Re_H=60\,000$ and $(b)$ $Re_H=120\,000$ .

4.5.2. Non-dimensional energy dissipation parameter

Figure 27 shows the variation of the non-dimensional energy dissipation parameter $C_\epsilon$ as a function of $Re_\lambda$ in the $Re_\lambda$ decay region. It is observed that $C_\epsilon$ decreases with the increase in $Re_\lambda$ . For $n\,=\,4$ , 6 and 8 at $Re_H=60\,000$ , part of the decaying region occurs at $Re_\lambda \leqslant 100$ , where the low-Reynolds number effect (viscosity effect) may also influence $C_\epsilon$ (Sreenivasan 1984, 1998; Bos, Shao & Bertoglio 2007; Kitamura et al. 2014; Vassilicos 2015; Kitamura et al. 2014). For $Re_\lambda \gt 100$ , we applied the least-squares fit method to evaluate the relationship, $C_\epsilon \propto Re_\lambda ^{-1}$ . We confirmed that for $n\,=\,4$ , 6 and 8 at $Re_H\,=\,60\,000$ , the exponents are approximately $-1.28$ , $-1.33$ and $-1.01$ , respectively. Similarly, for $n\,=\,4$ , 6 and 8 at $Re_H\,=\,120\,000$ , the exponents are approximately $-0.96$ , $-1.17$ and $-0.96$ , respectively. Notably, for $Re_H\,=\,60\,000$ , limited data points used for the fit may affect the representatives of the results. For $n\,=\,4$ at $Re_H\,=\,120\,000$ , we excluded four data points where $C_\epsilon \gt 0.4$ . Overall, for $Re_H\,=\,120\,000$ , $Re_\lambda \gt 100$ , $C_\epsilon$ is approximately proportional to $Re_\lambda ^{-1}$ , which qualitatively agrees with the previous studies conducted on fractal-grid-generated turbulence (Seoud & Vassilicos 2007; Mazellier & Vassilicos 2010; Valente & Vassilicos 2011b ; Gomes-Fernandes et al. 2012).

Figure 27. Non-dimensional energy dissipation parameter $C_\epsilon$ on the centreline in relation to $Re_\lambda$ in the decay region of $Re_\lambda$ for the fractal trees studied at $(a) Re_H=60\,000$ and $(b) Re_H=120\,000$ .

Figure 28. Non-dimensional energy dissipation parameter $C_\epsilon$ on the centreline in relation to $x /H$ for the fractal trees studied at $(a) Re_H=60\,000$ and $(b) Re_H=120\,000$ .

Figure 28 shows the variation of $C_\epsilon$ as a function of $x /H$ . Each fractal tree exhibits a near wake region where $C_\epsilon$ increases rapidly, followed by a downstream region where $C_\epsilon$ remains approximately constant. The boundary between these two regions occurs at approximately $x\,=\,4.0$ $H$ for $n\,=\,4$ , 6 and 8 at $Re_H\,=\,60\,000$ and approximately $x\,=\,4.0-5.0$ $H$ for $n\,=\,4$ , 6 and 8 at $Re_H\,=\,120\,000$ .

4.5.3. The TKE production and transverse transport

This section presents additional analysis of the transverse scans behind all three fractal trees to obtain a full picture of the evolution of wake turbulence. It is worth investigating whether the non-equilibrium dissipation scaling observed here is accompanied by non-zero turbulence production and/or triple correlation transport as in Valente & Vassilicos (2011b ), Nagata et al. (2013) and Hearst & Lavoie (2016). These two parameters are evaluated in the context of the following TKE equation:

(4.5) \begin{align} \begin{split} \frac {U_k}{2} \frac {\partial \overline {q^2}}{\partial x_k} = - \overline {u_i u_j} \frac {\partial U_i}{\partial x_j} - \frac {\partial }{\partial x_k} \left ( \frac {\overline {u_k q^2} }{2} \right ) - \frac {\partial }{\partial x_k} \left ( \frac {\overline {u_k p} }{\rho } \right ) + \frac {\nu }{2} \frac {\partial ^2 \overline {q^2} }{\partial x_k \partial x_k} - \epsilon. \end{split} \end{align}

The left-hand side term represents the advection term (TKE decay). The right-hand side terms represent production, triple-correlation transport, pressure transport, viscous diffusion and dissipation. Here, $U_1$ = $U$ , $U_2$ = $V$ and $U_3$ = $W$ denote the mean velocities, while $u_1$ = $u$ , $u_2$ = $v$ and $u_3$ = $w$ and $p$ represent the fluctuating velocities and pressure.

Referring to the analysis methods of Valente & Vassilicos (2011b ) and Hearst & Lavoie (2016), figure 29 shows the centreline TKE production $\mathcal {P}$ and triple-correlation transport $\mathcal {T}$ , both of which are normalized by the energy dissipation $\epsilon$ . As seen in figure 29 $(a)$ , near the upstream region (around $x /H = 1.0$ ), the TKE production is strong, reaching approximately $-150\%$ of the dissipation for $n\,=\,4$ , 6 and 8, before gradually decreasing. In the range of approximately $x/H\,=\,2.0-4.5$ , the production term is at most approximately $-70\%$ , $-50\%$ and $-15\%$ of the dissipation term for $n\,=\,4$ , 6 and 8, respectively. For $x /H \gt 4.5$ , the production-to-dissipation ratio gradually increases again in absolute value. As shown in figure 29 $(b)$ , the TKE transverse transport remains strong in both the upstream and downstream regions.

Figure 29. $(a)$ Centreline TKE production $\mathcal {P}$ and $(b)$ centreline triple-correlation transport $\mathcal {T}$ in relation to $x /H$ for Basic $n\,=\,4$ , 6 and 8 at $Re_H\,=\,120\,000$ .

Figure 30 $(a, b)$ shows the transverse profiles of $\mathcal {P}$ and $\mathcal {T}$ for $n\,=\,6$ at $Re_H\,=\,120\,000$ , normalized by the spanwise mean dissipation $\langle \epsilon \rangle _{y}$ . For clarity, in figure 30, we only show the results for $n\,=\,6$ . Note that we also confirmed similar behaviour for $n\,=\,4$ and 8. From $x /H\,=\,3.0-7.0$ , the profiles of $\mathcal {P}$ along the $y$ -direction remain inhomogeneous, indicating that the flow field is also inhomogeneous, with the transverse profiles of $\mathcal {P} / \langle \epsilon \rangle _{y}$ showing an overall increase in absolute value. The $\mathcal {T}$ -to-dissipation ratio remains at a non-negligible level. The absolute values of these two ratios can reach approximately 150 $\%$ and 200 $\%$ of $\langle \epsilon \rangle _{y}$ , respectively, which is significantly higher than those observed in the near wake region of the fractal grids (Hearst & Lavoie 2016). Figure 30 $(c, d)$ show the transverse profiles of $\overline {uv} (\partial U / \partial y)$ and $\overline {uw} (\partial U / \partial z)$ , which contribute significantly to $\mathcal {P}$ , also normalized by $\langle \epsilon \rangle _{y}$ . As the turbulence evolves downstream, the profiles of $\overline {uv} (\partial U / \partial y)$ gradually approach zero and become approximately homogeneous, while the increase in $\overline {uw} (\partial U / \partial z)$ exceeds the decrease in $\overline {uv} (\partial U / \partial y)$ . This indicates that the flow field becomes approximately homogeneous in the $y$ -direction farther downstream but remains inhomogeneous in the $z$ -direction, which explains the significant increase in the TKE production-to-dissipation ratio downstream in figure 30 $(a)$ .

Figure 30. Transverse profiles of $(a)$ TKE production $\mathcal {P}$ , $(b)$ triple-correlation transport $\mathcal {T}$ , $(c)$ $\overline {uv} (\partial U / \partial y)$ and $(d)$ $\overline {uw} (\partial U / \partial z)$ at various downstream positions for Basic $n\,=\,6$ at $z /H=0.5, Re_H = 120\,000$ .

Based on the above analysis, the non-equilibrium near wake region (approximately within $x /H \lt 5.0$ at $Re_H = 120\,000$ ) is highly inhomogeneous for the fractal trees, with significant TKE production and transverse transport. This aligns qualitatively with previous studies on fractal-grid-generated turbulence (Hearst & Lavoie 2016), which compared fractal square element grids with regular grids. They found that the fractal square element grids’ non-equilibrium near wake region was significantly more inhomogeneous, with higher TKE production and transverse transport when compared with the regular grids. In this study, the fractal tree exhibits even more significant 3-D inhomogeneities (in both the $y$ - and $z$ -directions) than the fractal grids. Thus, it can be considered that the inhomogeneous shape of the tree leads to greater inhomogeneity and TKE production and transport in the non-equilibrium near wake region.

However, it should be noted that non-equilibrium dissipation is not necessarily due to one-point inhomogeneities (Nagata et al. 2017). While these one-point inhomogeneities cause TKE production and transport, as observed in this study and previous research (Nagata et al. 2013; Valente & Vassilicos 2014; Hearst & Lavoie 2016), they do not affect the two-point energy equation at the scales where the inertial energy cascade occurs (Valente & Vassilicos 2015). As Nagata et al. (2017) pointed out, the cascade mechanism of the turbulence dissipation is not affected by one-point inhomogeneities. The inhomogeneity of a two-point energy balance at the correct length scales is more important than the inhomogeneity of the one-point energy balance. However, analysing two-point statistics is beyond this paper’s scope and requires further investigation.

5. Conclusions

This study successfully employed large-scale numerical simulations for the fluid around fractal trees using AMR–LBM. The main contribution of this study is the elucidation of the relationship of the drag coefficient of trees with the tree shape and Reynolds number, as well as the non-equilibrium dissipation behaviour $C_\epsilon \propto 1 /Re_\lambda$ in the wake region of fractal trees. This behaviour is qualitatively consistent with the non-equilibrium decay turbulence observed in fractal grids (Seoud & Vassilicos 2007; Gomes-Fernandes, Ganapathisubramani & Vassilicos 2012). Moreover, various turbulence statistics and scales were also explored, including turbulence intensity, global and local isotropy parameters and Taylor-microscale-based Reynolds numbers.

Three configurations of the same type of fractal tree geometries with fractal iteration numbers ( $n\,=\,4$ , 6 and 8) were generated through the parametric L-system, with a fractal dimension of 1.77. At low and medium Reynolds numbers, the drag coefficient increases with the fractal iteration numbers ( $n$ ) and decreases with the Reynolds number. However, at high Reynolds numbers ( $Re_H\,=\,120\,000$ ), if the fractal iteration number is greater than a certain value ( $n \geqslant\,6$ ), the drag coefficient becomes independent of the fractal iteration number ( $n$ ) and Reynolds number and tends to a constant value.

The energy dissipation scaling was also investigated. We confirm that, at high Reynolds numbers ( $Re_H \geqslant 60\,000$ ), the ratio of the integral length scale to the Taylor microscale $L_u / \lambda$ remains approximately constant while Taylor-microscale-based Reynolds number $Re_\lambda$ decreases in the wake region of $x / x_{peak} \geqslant 1$ . Consequently, the energy dissipation in the fractal tree’s near wake region does not follow the conventional scaling $\epsilon = C_\epsilon u^{\prime3} / L$ , as $C_\epsilon$ is not a constant, but approximately $C_\epsilon \propto 1 /Re_\lambda$ , which is qualitatively consistent with the non-equilibrium dissipation behaviour reported in the wake region of planar fractal grids. Here $C_\epsilon$ becomes approximately constant farther downstream. These findings can be crucial for improving the accuracy of micrometeorology predictions and human body-temperature simulations in urban environments, as current simulations do not account for the non-equilibrium dissipation properties of tree models. We believe that our results provide a valuable database and will contribute to future research. We find high one-point inhomogeneity within the non-equilibrium near wake region, with significant production and transverse transport of TKE, which qualitatively agrees with previous studies on the near field of fractal grid turbulence (Valente & Vassilicos 2011 b; Nagata et al. 2013; Hearst & Lavoie 2016).

In terms of other turbulence statistics, we found that the normalization method used by Gomes-Fernandes, Ganapathisubramani & Vassilicos (2012) for normalizing the centreline turbulence intensity can also be used to normalize the centreline turbulence intensity in the wake region of fractal trees. At high Reynolds numbers (e.g. $Re_H \geqslant 60\,000$ ), this normalization method can collapse the turbulence intensity of fractal trees with different shapes if the fractal iteration number reaches a certain threshold (e.g. $n\,=\,6$ ). However, this method cannot be used to collapse $Re_\lambda$ . Furthermore, the isotropy characteristics were investigated. At high Reynolds numbers ( $Re_H \geqslant 60\,000$ ), with $x / x_{peak} \geqslant\,1.0$ , the wake region of fractal trees can be considered essentially isotropic at all heights. In addition, the height dependence of turbulence intensity and $Re_\lambda$ was examined in the wake region, the maximum value of both turbulence intensity and $Re_\lambda$ appeared approximately at the height near the trunk and the first branch ( $z /H\,=\,0.35-0.5$ ).