3.1.1. Macroscale momentum budget

The macroscale momentum budget is computed using the macroscale momentum conservation equation (de Lemos Reference de Lemos2012). The macroscale momentum conservation equations are derived from the Navier–Stokes equations (2.11) by applying the VAT. The authors note that other forms of the conservation equation may be adopted to suit the modelling requirements, such as in Lasseux, Valdés-Parada & Bellet (Reference Lasseux, Valdés-Parada and Bellet2019). The time dependence of the conservation equation is retained. Applying this to the transient flow through the periodic porous medium, the macroscale spatial derivatives are eliminated to result in the following:

(3.1)\begin{equation}\underbrace{{\rho \frac{\partial }{{\partial t}}(\varphi {{\langle {u_i}\rangle }^i})}}_{{inertial}} = \underbrace{{\rho \varphi {g_i}}}_{{applied}} + \underbrace{{\frac{\mu }{{\Delta V}}\int_{{A_{interface}}} {{n_j}{\partial _j}{u_i}\,\textrm{d}S} }}_{{viscous\;drag}} - \underbrace{{\frac{1}{{\Delta V}}\int_{{A_{interface}}} {{n_i}p\,\textrm{d}S} }}_{{pressure\;drag}}.\end{equation}

The operator ${\langle \ \rangle ^i}$ denotes an intrinsic volume average in the fluid domain, ΔV is the total volume of the periodic REV, A_interface is the interfacial area and n_i is the normal vector of A_interface. The value of ΔV for all the cases with a 4 × 4 GPM is 3.2 × 10⁻⁵ m³. The value of A_interface can be calculated as $50(1 - \varphi )\Delta V$. The sum of the macroscale pressure drag and viscous drag in (3.1) (referred to hereafter as macroscale pressure and viscous forces) gives the total drag R_i. In this paper, the term drag is used to denote any force acting on the solid obstacle surface. In this porous medium system, an applied pressure gradient is used to sustain mass flow, all the forces that oppose it are called drag regardless of the convention. The authors note that there is a different naming convention followed in the aerodynamics community. At steady state, the Reynolds average of the drag $\langle {R_i}\rangle $ is balanced solely by the contribution of the applied pressure gradient $\langle \rho {g_i}\rangle $. The applied pressure gradient acts as a source of mechanical energy in the REV. The drag force on the solid obstacles will act as a sink of mechanical energy in the REV. The drag force in the x direction needs to be compensated for by the applied pressure to sustain the flow. In the case of symmetric flow, symmetry in the stress distribution on the solid obstacle surface results in zero pressure and viscous drag forces acting in the y direction after Reynolds averaging. In the case of deviatory flow, stress distribution on the solid obstacle surface is not symmetric about the y axis (Appendix C). Since the conservation of mechanical energy is satisfied by the numerical method in all three directions, a net zero y-direction drag force acts on the solid obstacle surface in deviatory flow. Note that there is no applied pressure gradient in the y direction to counteract the flow.

The question arises as to how a net zero drag force is possible in the absence of symmetric stress distribution on the solid obstacle surface. The answer lies in the balance in the components of drag force themselves. Note that we are still considering the Reynolds-averaged flow field where the inertial term in (3.1) vanishes. The integration of the pressure and viscous stress distributions on the solid obstacle surface yields non-zero, residual y-direction pressure and viscous force values that have ~5 % the magnitude of the x-direction pressure drag force. The y-direction pressure and y-direction viscous drag forces have equal magnitude and act in opposite directions resulting in a net zero drag force (see figure 2a). This unique feature of the Reynolds-averaged deviatory flow is elucidated in Appendix C. It is shown in the following discussion that the pressure force is substantially larger than the viscous force in the x direction, consistent with the expectations at the range of Reynolds numbers used in this work (100 < Re_p < 10 000). It is also shown that the viscous force has a small magnitude when compared to the pressure and inertial forces for the instantaneous flow in all three directions.

Figure 2. (a) A two-dimensional sketch of the forces that are acting on the solid obstacle in deviatory flow (the force vectors are not to scale). (b) The macroscale x-direction momentum budget computed for the laminar case A1. The individual components of the budget are calculated in the units of force per unit volume and then normalized by the applied pressure force per unit volume g ₁.

The symmetry-breaking phenomenon is considered a bifurcation because the deviatory flow solution is equally probable in both the positive and negative y directions (in figure 2a, the positive deviation is illustrated). The polarity of the macroscale y-direction viscous and pressure forces will interchange for the two possible symmetry-breaking solutions. Several modes of deviatory flow are possible from the two solutions, which are discussed at a later stage. Therefore, the interplay between the components in the x and y momentum budget is different. In the x direction, the pressure drag and the applied pressure gradient are dominant. Consider the case A1 (table 1) of laminar flow through the porous medium. The balance of x-direction forces for case A1 is illustrated in figure 2(b). The forces in the y and z directions are five orders of magnitude less than those in the x direction. For case A1, the microscale flow is strictly two-dimensional and the macroscale flow can be considered one-dimensional. Even at laminar Reynolds numbers, the macroscale pressure force is four times the magnitude of the macroscale viscous force. It will become apparent later on that the macroscale pressure force and the microscale pressure distribution dominate the flow properties for all of the simulations. In the y direction, both components of the drag – pressure and viscous – are relevant only in the Reynolds-averaged deviatory flow. It is prudent to bear this in mind while analysing the dynamic response of the fluid flow.

For packed beds, the limit of the laminar flow regime is at a pore-scale Reynolds number of 180 (Seguin et al. Reference Seguin, Montillet and Comiti1998a). The onset of the fully turbulent regime in finite porous media is gradual and can occur in the range 180 < Re_p < 900 (Seguin et al. Reference Seguin, Montillet, Comiti and Huet1998b). The transition regime may not exist for periodic porous media since the flow field is continuously perturbed and fed back into the REV. Transition and intermittency may be possible only through the local re-laminarization of the flow. However, transition effects are not encountered in any of the simulations in this paper. Turbulence has been observed in infinite porous media for Reynolds numbers as low as 478 for circular-cylinder solid obstacles (Uth et al. Reference Uth, Jin, Kuznetsov and Herwig2016). In this work, the dynamic nature of the flow emerges at a Reynolds number of Re_p = 225 (case A2, table 1). The macroscale momentum budget for case A2 is plotted against time in all three directions in figure 3. Note that the applied pressure gradient is maintained at a constant value for this case.

Figure 3. The dynamic macroscale momentum budget for case A2 (turbulent, Re_p = 225). The budget is computed for the (a) x, (b) y and (c) z directions. The time is non-dimensionalized using d/u_m, and the force components using the applied force component.

We first look at the properties of macroscale turbulence before symmetry breaking (cases A2 and A3). For the transient flow, the macroscale pressure and inertial components are dominant in the x and y directions when compared to the viscous component. Since the flow is symmetric, the viscous force in the y direction is virtually zero. The macroscale pressure force is zero in the z direction since the solid obstacles are cylindrical in shape with their axes oriented in the z direction. From the macroscale perspective, the random fluctuation of the components of the budget indicates that the flow at Re_p = 225 is turbulent. The order of magnitude of the forces in the z direction is only one order less than that in the y direction, which is indicative of three-dimensionality. We primarily study the momentum budget in the x and y directions since we only observe small perturbations in the z direction. The development of symmetry breaking is not studied starting from Re_p = 225 because the flow after symmetry breaking is vastly different, as shown in the following sections. Thus, the variation in the flow properties cannot be solely attributed to the symmetry-breaking phenomenon, but also to the development of turbulence.

The von Kármán instability is present at Re_p = 225, which introduces dynamic behaviour to the flow (figure 3). The macroscale force is calculated by summing the forces acting on the 16 solid obstacles present in the REV-T. The pressure force acting on the individual solid obstacles and their Fourier transform (absolute values) are plotted in figure 4. There are two features to note in the pressure force plots in the time and frequency domain: (1) the existence of two distinct, dominant time scales and (2) the phase difference between the drag forces acting on individual solid obstacles. The phase difference can be identified in figures 4(a) and 4(c), where the peaks in the pressure force on individual solid obstacles do not coincide in time. The high-frequency peaks at non-dimensional frequency f = 3.8 in figures 4(b) and 4(d) correspond to the vortex shedding frequency. This is corroborated with the peak frequencies in y velocity that are measured at the midpoint between two solid obstacles. The results are also supported by the visualization of the vortex shedding process using the streamlines and the skin-friction lines on the surface of the solid obstacle. We use the non-dimensional frequency f to analyse the time series of force components. The frequency is defined as the reciprocal of time and it is non-dimensionalized using the characteristic length and velocity. The form of the non-dimensional frequency f is similar to that of the Strouhal number. The Strouhal number is used in studies of the von Kármán vortex shedding frequency. In this paper, other flow instabilities are present due to the presence of an array of solid obstacles in the porous medium geometry. The porous medium geometry introduces a spectrum of frequencies of vortex shedding even at low Reynolds numbers. It is not common practice to use the Strouhal number for analysis in these conditions. However, note that the value of f at the vortex shedding frequency in some cases (e.g. case A2) is similar to the Strouhal number. Note that the Strouhal number does not take a constant value in porous medium flows at a moderate Reynolds number, unlike the flow around a single circular cylinder. Therefore, there is no need to use the Strouhal number in this paper.

Figure 4. The pressure force acting on individual solid obstacles inside the REV-T for case A2 (table 1) in the x direction versus non-dimensional (a) time and (b) frequency, and in the y direction versus non-dimensional (c) time and (d) frequency. See pressure drag definition in (3.1). The colours represent the location of the solid obstacle in the 4 × 4 matrix. Example: (2,3) refers to the solid obstacle in second place in the x direction and third place in the y direction.

The low-frequency peak is observed at two different values of frequency in the x and y directions (f = 0.3 and f = 0.6). The discrepancy is attributed to the finite sample size of the signal used for the fast Fourier transform and the peaks are, therefore, associated with the same phenomenon. The low-frequency oscillations arise from a secondary flow instability. The plots of the pressure force in the frequency domain (figures 4b and 4d) are qualitatively similar to the velocity probe measurements that are reported by Agnaou et al. (Reference Agnaou, Lasseux and Ahmadi2016). The low-frequency peak is a derivative of a Hopf bifurcation in periodic porous media. Since the porosity is low in this case (φ = 0.5), there exists a strong interaction between the primary and secondary flows. The macroscale pressure and inertial components of the macroscale momentum budget are in competition (figure 3), which results in the formation of the secondary flow instability. The presence of phase difference between the pressure forces acting on the solid obstacles alleviates the high magnitude of macroscale adverse pressure gradient that is introduced by the flow instability.

At higher Reynolds numbers Re_p ≥ 300 (cases A3 to A7), the power spectrum of the pressure force is spread over a band of frequencies. Sharp peaks like those seen in figure 4(d) for case A2 become more diffuse at high Reynolds numbers, which is a characteristic of turbulence. The pressure forces for Re_p = 300 and Re_p = 489 (cases A3 and A4) are plotted versus f in figure 5. For case A3, there are three dominant frequencies for the pressure force, one more than for case A2. The peak at f = 0.333 corresponds to the secondary flow instability that was present in case A2. The peak at f = 3.331 corresponds to the vortex shedding process. The frequency band 0.333 < f < 3.331 is also excited as a result of turbulence. The oscillation in the dimensional pressure force has a higher amplitude in case A3 than in case A2. The higher amplitude alludes to an increased strength of the vortex shedding process, which is verified using the magnitude of vorticity associated with the vortices. In figure 5(a), however, the non-dimensionalization of the pressure force using the applied pressure gradient shows a decrease in the magnitude due to the nonlinear increase in mean drag force. For case A4, a band of frequencies is excited with no distinct peaks in the power spectrum. The vortex shedding frequency is no longer one of the dominant frequencies for the oscillation in the pressure force. The maximum power is observed in the frequency range 1 < f < 3 in figure 5(b), which is one order of magnitude less than the mean vortex shedding frequency. The low-frequency oscillations of the flow instability that were present in case A3 are not present in case A4. The change in the dynamics of the flow from case A3 to A4 is brought about by the symmetry-breaking process. The flow after symmetry breaking has only one stagnation point on the solid obstacle, as opposed to the two stagnation points observed in symmetric flow. After symmetry breaking, the stagnation point is not formed by the incidence of a vortex with the solid obstacle (see § 3.1.2).

Figure 5. The pressure force acting in the y direction on the individual solid obstacles inside the REV-T versus non-dimensional frequency for (a) Re_p = 300 (case A3) and (b) Re_p = 489 (case A4).

As mentioned before, the phase difference in the drag forces between individual solid obstacles arises from asynchronous vortex shedding. The phase difference in the drag forces on individual solid obstacles has an influence on the macroscale flow field in a manner similar to wave interference phenomena. It causes a reduction in the amplitude of oscillation of the total drag force inside the REV. The phase difference in the vortex shedding increases the dispersion in the microscale velocity fluctuations. This is an important consideration in the definition of the macroscale Reynolds stress terms used in turbulence modelling. Although the order in which volume and Reynolds averaging are applied to dependent quantities in the governing equations does not change the result (Pedras & de Lemos Reference Pedras and de Lemos2001), the order in which these two averaging procedures are applied to model parameters, such as the macroscale turbulence kinetic energy, may change the result. The vortex shedding inside the porous medium provides a physical reason for why the two definitions of the macroscale turbulence kinetic energy that are used in porous medium turbulence modelling simulate entirely different physical phenomena (reported in de Lemos (Reference de Lemos2012) and Pinson, Grégoire & Simonin (Reference Pinson, Grégoire and Simonin2006)). The two commonly used definitions are ${\langle \langle {u^{\prime}_i}{u^{\prime}_i}\rangle /2\rangle ^i}$ and $\langle {\langle {u^{\prime}_i}\rangle ^i}{\langle {u^{\prime}_i}\rangle ^i}\rangle /2$. The prime denotes a fluctuation term obtained from Reynolds averaging. It is noted in the work of Kuwata, Suga & Sakurai (Reference Kuwata, Suga and Sakurai2014) that the first definition is the total turbulent kinetic energy and that the second definition is a macroscale turbulent kinetic energy. According to this, the total turbulent kinetic energy can be calculated as the sum of a macroscale and a microscale turbulent kinetic energy. The macroscale turbulent kinetic energy contributes a very small fraction of the total turbulent kinetic energy. This observation is also reported by Kuwata & Suga (Reference Kuwata and Suga2015). The dominance of microscale turbulent kinetic energy is attributed to large-scale oscillations due to vortex shedding. The phase difference in the vortex shedding behind the individual solid obstacles is also a factor, which increases the dispersion of the velocity fluctuations inside the REV. These observations can be extended to the other terms of the Reynolds stress tensor (RST) as well. Since the first definition $({\langle \langle {u^{\prime}_i}{u^{\prime}_j}\rangle /2\rangle ^i})$ is more commonly used for modelling, we use this definition of the Reynolds stress later in § 3.2. It encompasses all the relevant physics observed inside the porous medium. The presence of phase difference also validates the need for large REVs to simulate turbulent porous media flows. The use of a domain with a single solid obstacle will only simulate a special case where all of the vortex systems are acting in phase. The size of the REV-T should be sufficiently large that the influence of phase difference becomes invariant, which we demonstrate in Appendix A.

To summarize the discussion about the transition from symmetric to deviatory flow, changes in both the flow topology and the flow dynamics are expected. From the observed change in the power spectrum before and after symmetry breaking, the underlying physics of the phenomenon is related to the following:

• • The presence of a flow instability before symmetry breaking that causes large-scale pressure oscillations and its disappearance after symmetry breaking.

• • The increase in the amplitude of the oscillation of pressure until symmetry breaking occurs, followed by a decrease.

• • The reduction in the dominance of the vortex shedding process on the pressure force after symmetry breaking.

The results suggest that the symmetry breaking is brought about by the amplification of the flow instability with the increase in Reynolds number, which ultimately leads to the breakdown of symmetry. To confirm this, the transient stages of symmetry breaking are simulated using the following methodology. At the beginning of the simulation, the applied pressure gradient is maintained at a constant value such that the volumetric flow rate corresponds to a Reynolds number of 300. Then, a step function is used to change the applied pressure gradient at tu_m/d = 0 to that of a Reynolds number of 489. The constant value of applied pressure gradient that sustains the Reynolds number of 500 is determined a priori to use in the present simulation. The discrepancy between the prescribed Reynolds number (Re_p = 500) and the actual Reynolds number (Re_p = 489) is a result of the pathway that the bifurcation takes in this simulation. The step function is used to avoid any oscillations that will be introduced by using a control system to maintain the flow rate. This enabled the study of the transient stages towards deviatory flow without having to account for the large-scale noise that would have come from the time-advancement algorithm. However, it is not possible to segregate the momentum contribution of the symmetry-breaking phenomenon and the change in applied pressure gradient, since they occur concurrently. Therefore, the analysis of the macroscale momentum budget is supplemented by an inspection of the microscale flow field (§ 3.1.2) to confirm the observations.

All of the components of the macroscale momentum budget at tu_m/d = 0 amplify in magnitude as a result of the increased supply of momentum (figure 6). The three-dimensionality of the flow increases as a result of the increase in the Reynolds number (see figure 6c). However, the magnitudes of the forces in the z direction are two-orders-of-magnitude less than those in the x and y directions. It is sufficient to take only the x and y directions into consideration for this macroscale analysis, while bearing in mind that the microscale flow is three-dimensional. Both the macro- and microscale turbulence are strongly anisotropic, such that the streamwise direction is dominant. Before symmetry breaking, the streamwise direction is aligned with the x direction. After symmetry breaking, the deviatory flow results in a macroscale flow angle (θ_macro) between the streamwise direction and the x direction in the xy plane. The macroscale flow angle θ_macro is calculated as the deviation of the direction of the macroscale velocity vector from the x direction.

Figure 6. The macroscale momentum budget for the symmetry-breaking process (Re_p = 300 to 489). The budget is computed for the (a) x, (b) y and (c) z directions. The time axis is shifted such that the step change in applied pressure gradient occurs at t = 0.

Among the components of the macroscale momentum budget, the macroscale inertial and pressure forces dominate the dynamics of the macroscale flow. The viscous force is essential for the formation of the vortices and the flow instabilities, but it is not the driving component in symmetry breaking. In the x direction, the change in the applied pressure gradient $(\rho {g_i})$ leads to an increase in the macroscale x-direction pressure force and its amplitude of oscillation (figure 6a). The macroscale x-direction pressure force does not grow after tu_m/d = 5 onwards as long as the applied pressure gradient is unchanged. The amplitude of oscillation of the macroscale x-direction pressure force decreases after tu_m/d = 4, which is followed by the onset of macroscale symmetry breaking at tu_m/d = 7.3 (see macroscale y-direction pressure force in figure 6b). In the y direction, the change in the applied pressure gradient increases the amplitude of oscillation of the macroscale y-direction pressure force. Small-amplitude oscillations in the macroscale y-direction pressure force are first observed in 0 < tu_m/d < 7.3 as a direct result of the change in applied pressure gradient. Large-amplitude oscillations in the macroscale y-direction pressure force are observed for tu_m/d > 7.3 because of the deviatory flow from symmetry breaking. The change in the amplitude of oscillation of the pressure force due to symmetry breaking is more prevalent in the y direction (figure 6b) than in the x direction (figure 6a). In the y direction, amplitude of oscillation of the pressure force increases by one order of magnitude due to the deviatory flow. Therefore, the transformation of the microscale flow field in the xy plane after symmetry breaking is more evident in the y-direction macroscale momentum budget than in the x direction.

The macroscale y-direction pressure force oscillates about zero both before and after symmetry breaking. This is counterintuitive to the observations in figure 2(a), where a non-zero y-direction pressure force is balanced by a non-zero y-direction viscous force in order to sustain the deviatory flow. If the components of the y-direction macroscale momentum budget have a zero mean value after symmetry breaking, how can the deviatory flow solution exist? To answer this question, the pressure forces that are acting on the individual solid obstacles are examined (figure 7). The 4 × 4 matrix of solid obstacles is indexed by dividing the matrix into rows and columns. The rows are aligned with the x axis and the columns are aligned with the y axis.

Figure 7. The pressure forces acting in the (a) x and (b) y directions on the individual solid obstacles inside the REV-T versus time for the transient simulation (Re_p = 300 to 489).

After symmetry breaking has developed (tu_m/d > 10), the y-direction pressure forces on the individual solid obstacles in a single row have different magnitudes and signs. Two solid obstacles in the row have a positive y-direction pressure force and the other two have a negative y-direction pressure force, which results in a zero sum. The y-direction pressure forces on all of the solid obstacles in a single column have similar magnitudes and sign. We noted earlier that the y-direction pressure force on the solid obstacle can be positive or negative for a corresponding positive or negative macroscale flow angle. Since there are four columns of solid obstacles in this REV-T, the deviatory flow behind each solid obstacle can assume either solution behind each column. In this paper, each unique combination of the y-direction pressure forces on the individual solid obstacles is called a mode, analogous to wave modes. Modes are formed in the x direction alone, along which the pressure gradient is applied. The modes of deviatory flow need not be symmetric in the x direction. The modes are illustrated in § 3.1.2. If the REV-T consists of one solid obstacle, only a single mode of the symmetry-breaking phenomenon with a unidirectional y-direction pressure force can be formed. However, the case of unidirectional deviatory flow is only a subset of the possible solutions in periodic porous media. The number of possible modes of deviatory flow is decided by the number of solid obstacles in the REV in the direction of applied pressure gradient. This confirms the requirement that the size of the REV-T must be greater than one pore size.

The onset of macroscale symmetry breaking is evident in the x-direction pressure forces at tu_m/d = 7.3 (figure 7a). Since the magnitude of the change of the x-direction pressure force at tu_m/d = 0 is large, the amplitude of oscillation of the x-direction pressure force is small in relation. Therefore, the y-direction pressure forces are used to analyse the symmetry-breaking phenomenon behind individual solid obstacles. There exists a phase difference between the pressure forces on the individual solid obstacles (figure 7). The onset of microscale symmetry breaking is marked by the first increase in the y-direction pressure force on the solid obstacles. The flow around the solid obstacle at location (2,3) is the first to deviate from symmetric behaviour (figure 7b). The deviation occurs sooner in the microscale (tu_m/d = 6) than it does in the macroscale (tu_m/d = 7.3), implying that symmetry breaking begins as a localized, microscale phenomenon.

The flows around the individual solid obstacles break symmetry at different times because of the phase difference in the vortex shedding. The solid obstacle at location (2,3) is defined to have the leading phase at the time of microscale symmetry breaking. The increasing order of phase lag for the solid obstacles in the third row of the REV-T is the following: (2,3) < (3,3) < (1,3) < (4,3). This order is identified by determining the solid obstacle locations where the y-direction pressure force will peak next after the solid obstacle at (2,3). The order in which the flow around each solid obstacle breaks symmetry is the same as that of the phase lag. Thus, the phase lag determines the time at which symmetry breaking occurs for each solid obstacle. Therefore, the phase of the vortex shedding cycle at which the flow breaks symmetry is different for each solid obstacle. The phase difference is only partly responsible for the formation of the different modes. It is shown in § 3.1.2 that the flow field around the neighbouring solid obstacle plays a role in the formation of the modes as well.

The formation of the phase difference and the randomness in the vortex motions behind the solid obstacles suggest the independence of the microscale flow behind each solid obstacle. The formation of modes of deviatory flow along the rows of solid obstacles suggests the dependence of the microscale flow behind a solid obstacle on the flow around the neighbouring obstacles. The flow patterns around the solid obstacles in the direction of the applied pressure gradient are not identical, but the flow patterns are identical in the normal direction.

3.1.2. Microscale flow field and turbulent structure visualization

The macroscale momentum budgets for cases A1–A4 in table 1 revealed that the flows before and after symmetry breaking are dominated by the pressure forces. The dynamic behaviour of the flow is determined by the microvortices, formed as a result of the viscosity of the fluid. However, the viscous force has little influence in the macroscale momentum budget. In this section, the microscale flow field is visualized to connect the observations of the macroscale flow to the microscale flow physics. The laminar flow patterns at Re_p = 100 are similar to those in figure 1(a), with distinct primary and secondary flow regions. The secondary flow region consists of a recirculating (attached) vortex system. Upon transition to turbulence, three-dimensional features appear in the microscale flow at Re_p = 225, consistent with the macroscale momentum budget in figure 3. Even though the microvortices at Re_p = 225 and at Re_p = 100 appear similar, they begin to deform in the z direction at Re_p = 225 due to the vortex stretching process in turbulence. The two-dimensional turbulent structures in figure 8(a) are the microvortices, which possess a swirling characteristic. Note that we are using the term $\rho {g_1}\Delta V$ to scale the microscale pressure distribution to be consistent with the definition of pressure drag force in (3.1). The pressure drag force can be calculated from the pressure distribution by integrating the microscale pressure distribution over the area of the solid obstacle.

Figure 8. Coherent turbulent structures visualized using the iso-surfaces of the Q-criterion for (a) Re_p = 225 (case A2) and (b) Re_p = 300 (case A3). The Q-criterion is normalized using the maximum value and the iso-surfaces of Q/Q_max = 0.001 are plotted. A sub-volume of dimensions (3s, 2s, 2s) of the REV-T is shown.

The microvortices at Re_p = 225 evolve in time due to turbulence production and dissipation. The microvortices are stationary in the void space between two solid obstacles. The microvortices interact with the primary flow to form flow stagnation points on the solid obstacle surface at the converging portion of the GPM. This leads to the formation of the secondary flow instability. The converging geometry intrinsically introduces a local streamwise favourable pressure gradient. However, the flow stagnation reduces the favourability of the pressure gradient in the converging section, which is the source of macroscale pressure drag. The ratio of the stagnation pressure to the applied pressure gradient increases with Reynolds number. The presence of a substantial stagnation pressure in the converging region supports the notion that the flow instability should arise from the competition between the pressure and inertial components of the macroscale momentum budget.

There is a stark contrast between the turbulent structures observed at Re_p = 225 (figure 8a) and Re_p = 300 (figure 8b). At Re_p = 225, semi-infinite two-dimensional turbulent structures with three-dimensional deformations are visible. The turbulent structure has infinite dimension in the z direction because of the periodic boundary condition. The deformations that are present in the turbulent structures arise from turbulent vortex stretching. At Re_p = 300, the turbulent structures are finite and three-dimensional and they assume the shape of hairpin vortices. The hairpin vortices are oriented in a reverse direction when compared to the observations in classic external flows like turbulent boundary layer flow (Eitel-Amor et al. Reference Eitel-Amor, Órlú, Schlatter and Flores2015). Unlike the boundary layer flow, the head-to-tail direction of the hairpin vortex aligns with the streamwise direction of the flow. Reverse hairpin vortices have also been observed in periodic porous media by Jin & Kuznetsov (Reference Jin and Kuznetsov2017) for spherical solid obstacles. To understand the unique vortex shape, the microscale turbulent structures are classified into the following based on the swirl: microvortices and turbulent eddies. The regions of swirl-dominated flow are identified by determining the centre of swirling flow using the algorithm developed by Sujudi & Haimes (Reference Sujudi and Haimes1995). The microvortices are characterized by swirling flow. The turbulent eddies do not possess a swirling vortex core. The turbulent structures that are observed in the Q-criterion plots (figure 8) in the regions without vortex core lines (figure 9a) are turbulent eddies.

Figure 9. (a) Instantaneous vortex core lines before and after symmetry breaking projected on a plane at z = 0 overlaid on contours of z-direction vorticity. (b) A sketch illustrating the process of formation of the reverse hairpin vortices in porous media.

The direction of the hairpin vortex structure is dependent on the microscale spatial inhomogeneity in turbulence dissipation, brought about by distinct primary and secondary flow regions. At Re_p = 300 (case A3), the microvortex core lines are concentrated in the secondary flow region (figure 9a). The swirling motion of the microvortex is localized in the secondary flow region. The primary and secondary flow regions are separated by a region of strong shearing flow (see vorticity contours in figure 9a). In order to understand the reversed direction of the hairpin vortices, consider an unperturbed vortex filament that is present in the secondary flow region (figure 9b). The microvortex is perturbed and stretched in the z direction by turbulence. As it stretches, the perturbed microvortex is subjected to the swirling secondary flow region and the highly dissipative primary flow region simultaneously. The low velocity and the rotational nature of the secondary flow region are favourable for the sustenance of the microvortex. The high rate of strain and the presence of a pressure gradient in the primary flow region are adverse to the sustenance of the microvortex. Therefore, the microvortex breaks up in the primary flow region, while it continues to sustain in the secondary flow region. This results in the formation of the reverse hairpin vortex with the head in the secondary flow region and the tail in the primary flow region.

The head of the reverse hairpin vortex will reduce in size as the vortex is gradually weakened due to energy dissipation in the flow. The phase difference in the microvortex motions behind the different solid obstacles is visualized by comparing the size of the head of the reverse hairpin vortices. In figure 8(b), the three-dimensional turbulent structures for the two adjacent solid obstacles show that the head of the reverse hairpin vortex is larger for the upstream vortices. Therefore, the upstream vortices are lagging behind resulting in a phase difference in the vortex motions. This observation is supported by the fact that the upstream vortices have a lower magnitude of pressure associated with them when compared to their neighbours. The phase difference in the pressure forces was also observed in figure 7(b). The reverse hairpin vortices in figure 8(b) possess some degree of order in their arrangement, forming a ‘knit’ pattern of three-dimensional turbulent structures. The formation of the pattern of similar turbulent structures that arise from the microvortices is consistent with the formation of significant peaks in the frequency distribution of the pressure forces (figure 5a). The microvortex patterns behind all of the solid obstacles have a similar size, shape and turnover time. The summation of the influence of an infinite number of these vortices on the solid obstacles leads to the excitation of the pressure force at the frequency of vortex shedding. The inherent randomness in the microvortex pattern is responsible for the excitation of the other frequencies.

Before symmetry breaking, microvortices are not transported beyond a single pore space (figure 9a). The microvortices transform into turbulent eddies in the shear-flow region. These turbulent eddies are observed to move downstream beyond a single pore space while continuing to diminish in size. The turbulent eddies have a random shape that does not bear resemblance to the reverse hairpin vortices. On this basis, the domain of influence of the microvortex motion behind one solid obstacle is limited to its neighbours. An indirect influence is present from the propagation of turbulent eddies until they are completely dissipated. After symmetry breaking, the distinction between the primary and the secondary flow regions is lost, as is the distinction between regions of microvortices and turbulent eddies. The flow separation point advances to a downstream location resulting in the formation of a small vortex wake and a weak shear layer. The microvortices propagate downstream and transform into turbulent eddies under the influence of the near-wall dissipation and the pressure gradient from the solid-obstacle geometry. The degree of confinement of the turbulent structures in the pore space is reduced by the occurrence of symmetry breaking evidenced by the microvortices travelling a greater distance before completely dissipating.

At Re_p = 489, the deviatory flow results in the formation of fine-scale turbulent structures (figure 10), when compared to Re_p = 300. Fine-scale turbulent structures are formed because of the reduction in the width of the microvortex wake after symmetry breaking. Fine-scale turbulent structures are also formed because of increased shredding of the turbulent structures at higher Reynolds numbers (Wood et al. Reference Wood, He and Apte2020). The ‘knit’ pattern in the microvortices (figure 8b) is no longer present due to increased turbulence intensity. This is reflected in the loss of dominance of the vortex shedding frequency on the pressure force and the distribution of the pressure force across a band of frequencies (figure 5b). The overall distribution of turbulent structures is chaotic. The deviatory flow is not evident from the Q-structures because of the formation of a mode of deviatory flow that increases the tortuosity of the streamlines (figure 10b). For this reason, the visualization of three-dimensional Q-structures is supplemented by the visualization of two-dimensional projections of the flow streamlines. The streamline plot shows that the flow around each solid obstacle consists of a single stagnation point and a detached vortex system. Separation bubbles are formed when the location of the stagnation point changes between adjacent solid obstacles. The low-pressure turbulent structures formed in the detached vortex system (at obstacles 1 and 4 in figure 10c) are microvortices with a vortex core line inside of them. The elongated, high-pressure turbulent structures near the stagnation point (at obstacles 2 and 3 in figure 10c) are microvortices that are impinging from the upstream neighbouring solid obstacle, marking the beginning of their transformation to turbulent eddies.

Figure 10. (a) Coherent turbulent structures visualized using the iso-surfaces of the Q-criterion for Re_p = 489 (case A4). The Q-criterion is normalized using the maximum value and the iso-surfaces of Q/Q_max = 0.02 are plotted. A sub-volume of dimensions (3s, 2s, 2s) of the REV-T is shown. (b) Instantaneous three-dimensional flow streamlines projected on a plane at z = 0 overlaid on contours of static pressure. (c) Coherent turbulent structures visualized using the Q-criterion for a sub-volume of dimensions (2s, 0.5s, 2s).

The three-dimensional turbulent structures observed in figures 10(a) and 10(c) after symmetry breaking bear no resemblance to the reverse hairpin structures and the ‘knit’ pattern observed before symmetry breaking (figure 8b). The reverse hairpin structures are formed by the vortex stretching of a two-dimensional vortex structure that is semi-infinite in the z direction at the time of formation. A virtually one-dimensional flow separation line is formed in the rear of the solid obstacle (figure 11a). After symmetry breaking (Re_p = 489), the size of the vortex core decreases and fine-scale structures are formed. The flow separation line after symmetry breaking is more complex (figure 11b). The flow separation lines are two-dimensional and discontinuous with local flow re-reattachment, which together result in the formation of finite three-dimensional vortex structures. The microvortices are the primary source of turbulent structures in periodic porous media. Therefore, the difference in the large-scale turbulent structures before and after symmetry breaking results in the change in the turbulence intensity and degree of anisotropy.

The origin of microscale symmetry breaking can be traced from the following differences in the microscale flow before and after symmetry breaking:

1. 1. The nonlinear increase in the microscale flow stagnation pressure with increase in the Reynolds number.

2. 2. The change in the location and the number of flow stagnation points. The change in the location of the flow separation points and the microvortex core size. The formation of fine-scale turbulent structures as a result of the change in microvortex size.

3. 3. The disappearance of the order in the flow patterns after symmetry breaking and the increased turbulence intensity and dissipation in the microscale flow.

The transient simulation that is shown in figure 6 is revisited to visualize the stages of symmetry breaking. The flow field is three-dimensional for the entire duration of the simulation. The microscale flow field is visualized using a two-dimensional projection of the instantaneous flow streamlines (figure 12). Microscale flow structures are visualized using three-dimensional coherent structures with the Q-criterion and two-dimensional contour plots of vorticity magnitude (figure 13; supplementary movies 1 and 2 available at https://doi.org/10.1017/jfm.2021.813). The stages of symmetry breakdown are summarized in table 2.

Figure 11. Skin-friction lines (black) plotted on the surface of the solid obstacles for the flow (a) before symmetry breaking and (b) after symmetry breaking. The red and green lines are iso-lines of zero shear stress indicating the flow separation line and the vortex regions, respectively.

Figure 12. The transient stages (a) A, (b) B, (c) C, (d) D and (e) E of the symmetry-breaking process from symmetric to deviatory flow. Instantaneous flow streamlines are seeded on the plane z = 0 and projected on the contours of pressure on the plane at z = −s.

Figure 13. Transition in the three-dimensional flow features before symmetry breakdown due to increase in g_i. Three-dimensional turbulent structures coloured by pressure are visualized in the sub-volume (3s, 2s, 2s) of the REV-T using iso-surfaces of $Q{(d/{\langle {u_1}\rangle ^i})^2} = 50$. Turbulent structures are identified by the vorticity magnitude distribution shown on the boundary of the REV-T. Vorticity contours below a magnitude of 5 are not displayed.

Table 2. The stages of symmetry breaking inferred from the transient simulation illustrated in figures 12 and 13.

The macroscale x-direction velocity increases when the applied pressure gradient is increased. The use of a constant value of the Q-criterion normalized by the maximum value for all of the time steps will reduce the visualization quality. Therefore, the Q-criterion is normalized using the solid-obstacle diameter d and the instantaneous macroscale x-direction velocity ${\langle {u_1}\rangle ^i}$. The microscale flow field at stage A is taken as the reference point before symmetry breaking (figure 12a). The flow streamlines show a large recirculation zone in the secondary flow region with a size similar to the radius of the circular-cylinder solid obstacles. The three-dimensional hairpin structures that were introduced in figure 8(b) are visible as hairpin-shaped streaks in the vorticity magnitude plot (figure 13c). The tail length of the reverse hairpin vortices increases from stage A to B, signifying an increase in vortex deformation due to higher strain rate in stage B compared to A. The three-dimensional nature of the flow is not apparent in the flow streamline plots, which is the case in subsequent time steps as well. However, the flow streamlines in the xy plane clearly indicate the onset and development of symmetry breaking.

Flow recirculation in the secondary flow region before symmetry breaking forms a strong pressure ‘plug’ in the flow. The word ‘plug’ is in quotation marks because it is not a solid plug and therefore allows mass to flow through it. Before symmetry breaking, two stagnation points are formed where the microvortex meets the solid obstacle on its top and bottom surfaces (θ ˜ ±${\rm \pi}$/3 in figure 11a). The stagnation points are located in the converging portion of the GPM (−${\rm \pi}$/2 < θ < 0 and 0 < θ < ${\rm \pi}$/2). The stagnation point creates a region of locally adverse pressure gradient, which is followed by the intrinsically adverse diverging portion of the GPM (${\rm \pi}$/2 < θ < ${\rm \pi}$ and −${\rm \pi}$ < θ < −${\rm \pi}$/2). The magnitude of stagnation pressure increases by a factor of 100 during symmetry breaking, when the applied pressure gradient is only increased by a factor of six (see figure 11b versus figure 11a). A pressure ‘plug’ in the flow is formed by the substantial increase in the stagnation pressure that increases the adversity of the pressure loss in the GPM. An adverse pressure gradient in the turbulent flow increases turbulence dissipation and reduces the size of the vortex structures (Lee & Sung Reference Lee and Sung2008; Tanarro, Vinuesa & Schlatter Reference Tanarro, Vinuesa and Schlatter2020).

Fine-scale turbulence structures are observed at tu_m/d = 4.71 as a result of the increased adverse pressure gradient and the increase in Reynolds number (compare figures 13a and 13c with 13b and 13d). The flow streamlines do not change significantly from stage A to C, implying that the velocity distribution of the large-scale microscale flow is unchanged. Recall that the flow transition from tu_m/d = 3.14 to 4.71 is characterized by an increase followed by a decrease in the amplitude of oscillation of the y-direction pressure force (figure 7). In stage C, the three-dimensionality of the flow field increases and the reverse hairpin vortex structures are no longer visible. The large-scale microvortex breaks down into smaller vortices that continue to possess the large-scale swirling motion. The fine-scale microvortices rotate about the same vortex core, seen at stage A. However, the large-scale swirling motion is weakened by the vortex breakdown and subjected to the continuous interaction of the constituent small turbulent structures.

When the pressure ‘plug’ is formed in the streamwise direction (x direction), the pressure gradient in the lateral direction (y direction) becomes favourable. The flow instabilities (von Kármán and secondary flow instabilities) cause asymmetry in the instantaneous pressure force that is acting on the solid obstacles. The asymmetry in the y-direction pressure, compounded by the increase in the stagnation pressure, will introduce a sufficient lateral favourable pressure gradient to sustain lateral flow. The y-direction pressure force sustains the deviatory flow after symmetry breaking (figure 2a). Therefore, only a trigger is required to induce symmetry breaking, after which the deviatory flow is sustained by the lateral favourable pressure gradient. At stage D (figure 12d), the flow behind the solid obstacles at column 2 of the REV-T begins to deviate from the symmetric flow configuration. Symmetry breaking is triggered by the breakup of the swirling vortex motion at the top. The vortex breakup is caused by the strong turbulent shear and the interaction between the weaker turbulent structures in the secondary flow region. The vortex breakup results in an efflux of fluid from the secondary flow region to the primary flow region. The efflux is compensated by the influx of fluid from the bottom recirculating vortex motion aided by the lateral favourable pressure gradient. The induced fluid flow in the y direction triggers the symmetry-breaking process.

The information about the symmetry breaking of the flow around solid obstacles in column 2 of the REV-T propagates to the neighbouring solid obstacles. The flow around the solid obstacles in column 3, downstream of column 2, is the next to break symmetry. The flow around the solid obstacles in column 1, upstream of column 2, also changes to accommodate the deviatory flow, followed by column 4. Symmetry breaking to deviatory flow propagates away from its location of incidence, similar to the propagation of perturbations in waves. The final mode of deviatory flow field after symmetry breaking is reflected in the flow streamlines at stage E (figure 12e). The mode is decided by the phase of the vortex motion at the time that the symmetry-breaking perturbation reaches the column of solid obstacles. The direction of the lateral favourable pressure gradient at the time chooses the direction of symmetry breaking. This behaviour is more discernible from the y-direction pressure force plots in figure 7(b). The mode of the deviatory flow after symmetry breaking is sustained forever, as long as the Reynolds number (Re_p) is maintained above the critical value for symmetry breaking. An investigation of the modes of deviatory flow due to symmetry breaking in porous media is beyond the scope of this paper.

3.2.1. Porosity

The symmetry-breaking phenomenon is a result of the influence of geometric confinement on the pressure distribution on the solid-obstacle surface. The porosity is the geometric parameter that controls the degree of flow confinement in the porous medium. The Reynolds number of the flow is fixed at Re_p = 1000 and the porosity is decreased in the range 0.80 ≥ φ ≥ 0.43 (table 3). The degree of flow symmetry breaking can be measured by using the macroscale flow angle θ_macro, introduced in § 3.1.1. The macroscale flow angle varies continuously with the porosity between zero before symmetry breaking and a finite value after symmetry breaking (figure 14a). The flow transitions from symmetric to deviatory flow at the critical value of porosity, which lies between 0.67 and 0.72 for the case of the circular-cylinder solid obstacles. Symmetry breaking is dependent on two conditions that are related to the porosity: (1) asymmetry in the stagnation pressures before symmetry breaking and (2) the proximity of the solid-obstacle surfaces. Together, these conditions create the lateral favourable pressure in the secondary flow region. This is why deviatory flow is non-existent at higher porosities even though the von Kármán instability is present.

Figure 14. (a) The macroscale flow angle θ_macro and the principal axis of the macroscale RST projected on the x axis θ_RST and (b) the components of the macroscale RST (note log scale on vertical axis) versus porosity φ for Re_p = 1000. (c) The Cartesian axes (geometric axes), the principal axes of the macroscale RST and the macroscale velocity vector in a three-dimensional unit sphere.

Table 3. The LES cases simulated to analyse the dependence of symmetry breaking on the porosity φ. The solid obstacles are circular cylinders and the Reynolds number (Re_p) is 1000 for all of these cases.

When the porosity is less than the critical value for symmetry breaking (φ < 0.72), two flow configurations are observed within the deviatory flow regime. The first is observed in the range 0.61 < φ < 0.72, where there is flow separation on either side of the plane of geometric symmetry. There exists an asymmetry in the vortex core diameter of the microvortices formed on either side of the plane of geometric symmetry. The microvortices enter the primary flow region and they diminish in size in the converging–diverging section of the solid-obstacle geometry. The second configuration is observed in the range 0.43 < φ < 0.61, where the flow separation occurs on the same side of the plane of geometric symmetry. The degree of asymmetry in the vortex core diameter is less in the second configuration than in the first configuration. In the second configuration, the microvortices impinge on the solid obstacles in the lateral direction.

Deviatory flow is a source of macroscale turbulence anisotropy in homogeneous porous media with symmetric solid obstacles. The non-diagonal terms in the macroscale RST are typically zero in such porous media because the microscale RST is an even function about the plane of geometric symmetry (Chu et al. Reference Chu, Weigand and Vaikuntanathan2018). The microscale and macroscale RSTs are defined in (3.2) and (3.3), respectively. The variable $u^{\prime}$ refers to the microscale velocity fluctuation. The naming convention followed in de Lemos (Reference de Lemos2012) is used here.

(3.2)\begin{gather}{\tau _{micro,ij}} = \langle {u^{\prime}_i}{u^{\prime}_j}\rangle ,\end{gather}

(3.3)\begin{gather}{\tau _{macro,ij}} = \varphi {\langle \langle {u^{\prime}_i}{u^{\prime}_j}\rangle \rangle ^i}.\end{gather}

The double-averaged governing equations contain other terms than ${\tau _{macro,ij}}$ that would also be influenced by the deviatory flow. Terms such as the turbulent dispersion stress ${\langle \langle {}^i{u^{\prime}_i}{}^i{u^{\prime}_j}\rangle \rangle ^i}$ and turbulent RST $\langle {\langle {u^{\prime}_i}\rangle ^i}{\langle {u^{\prime}_j}\rangle ^i}\rangle$ from macroscale velocity fluctuations also pertain to turbulence. However, the turbulent dispersion stress dominates the turbulent Reynolds stress and becomes virtually equivalent to the macroscale RST. This is shown in § 3.1.1 and also in the work of Kuwata & Suga (Reference Kuwata and Suga2015). The qualitative observations made below about the principal axis of the macroscale RST can be expected for the turbulence dispersion stress tensor as well. The non-turbulent stress terms such as dispersion stress and macroscale convection are not analysed here. A comprehensive budget of the double-averaged equations is beyond the scope of this paper but could be a valuable study in future work. In this paper, the analysis of the macroscale RST is adequate to indicate the implications of symmetry breaking on macroscale turbulence anisotropy. When the microscale flow symmetry is broken, the microscale RST ceases to be an even function about the plane of geometric symmetry. After volume averaging, the macroscale RST has non-zero non-diagonal terms (figure 14b). The magnitude of τ_{macro ,12} increases when the porosity decreases in the case of the first flow configuration (0.61 < φ < 0.72). The magnitude of τ_{macro ,12} decreases afterwards when the porosity decreases in the case of the second flow configuration (0.43 < φ < 0.61). The magnitude of the xy microscale RST τ_{micro ,12} is higher in the microvortex region. In the first flow configuration, the magnitude and the degree of asymmetry in the vortex core diameter are high, which results in the higher magnitude of τ_{macro ,12} in the first flow configuration. The components of the macroscale RST τ_{macro ,13} and τ_{macro ,23} are virtually zero, since they are three orders of magnitude less than the diagonal components. They are not exactly equal to zero since the simulation has not been averaged for infinite time.

Symmetry breaking in the microscale flow field causes symmetry breaking between the macroscale velocity vector and the macroscale RST. The principal axes of the macroscale RST are rotated by an angle θ_RST from the Cartesian axes. The orientation of the principal axis of the macroscale RST with respect to the Cartesian axes is determined by computing its eigenvectors. The direction vector of the macroscale velocity is not identical to the orientation of the principal macroscale Reynolds stresses θ_RST (figure 14a). In other words, the macroscale RST is not oriented in the direction of the macroscale flow. This suggests that the macroscale turbulence anisotropy cannot be taken into account by axis transformation. The non-diagonal terms in the macroscale RST will need to be explicitly modelled.

To summarize, the microscale flow breaks symmetry at a critical value of porosity due to a lateral favourable pressure gradient that sustains macroscale deviatory flow. After symmetry breaking, the macroscale velocity vector is oriented at an angle θ_macro with respect to the x axis. The influence of symmetry breaking extends beyond the macroscale velocity vector. If the coordinate axes are rotated by an angle of θ_macro such that the axes are aligned with the macroscale velocity vector, the macroscale RST is not oriented along the transformed axes. The principal axes of the macroscale RST are not aligned with either the macroscale velocity vector or the Cartesian axes (figure 14c). This indicates a complete breakdown of symmetry in both the microscale and macroscale flow without the presence of a single plane of symmetry in the flow.

3.2.2. Reynolds number

The Reynolds number is another critical parameter that determines the possibility of symmetry breaking, since it controls the magnitude of the stagnation pressure in the microscale flow. The porosity is set equal to 0.5 and the Reynolds number is increased from Re_p = 100 to 10 000 (table 1). The critical value of Reynolds number for symmetry breaking lies between Re_p = 300 and 500, and the details of its progression are discussed in §§ 3.1.1 and 3.1.2. A jump discontinuity is observed in the macroscale flow angle as the Reynolds number crosses the critical value. The critical Reynolds number is a necessary condition for symmetry breaking to occur. Unlike the porosity, symmetry breaking cannot be parameterized using the Reynolds number to control the deviatory flow. When the Reynolds number is changed, deviatory flow can either emerge or disappear resulting in a binary change with zero or a constant value of θ_macro depending on the porosity. Whereas, an increase in the porosity can result in a non-binary, incremental change in θ_macro.

The anisotropy in the macroscale RST of turbulence after symmetry breaking is evident for all values of Reynolds numbers above the critical value (figure 15). Microscale turbulence is expected to become more isotropic as the Reynolds number increases (Wood et al. Reference Wood, He and Apte2020). However, the deviatory flow after symmetry breaking leads to an increase in τ_{macro ,12} with the Reynolds number, which results in an increase in the macroscale anisotropy of turbulence. The consistent increase in τ_{macro ,12} with the Reynolds number implies that the macroscale anisotropy is not formed from the microscale turbulence structures. In § 3.1.2, the microscale turbulence structures after symmetry breaking were observed to be three-dimensional fine-scale structures as opposed to the two-dimensional large-scale structures observed before symmetry breaking. Therefore, the origin of macroscale anisotropy is purely the asymmetry in the microscale flow topology arising from symmetry breaking. Here, microscale flow topology refers to the microscale velocity and pressure distribution, the location of the primary and secondary flow regions and the location and the number of flow separation and stagnation points.

3.2.3. Solid-obstacle shape

The shape of the solid obstacles influences the vortex shedding process and the location of flow stagnation points, which are the source of the symmetry-breaking phenomenon. Therefore, deviatory flow cannot occur for all solid-obstacle shapes. Take for instance a square-cylinder solid obstacle forming a porous medium with a porosity of 0.5. This value of porosity is within the critical porosity for symmetry breaking in a porous matrix with circular-cylinder obstacles. Deviatory flow is not observed for the square-cylinder solid obstacles for any of the simulated cases (table 4). The Hopf bifurcation has been reported to occur for square solid obstacles (Agnaou et al. Reference Agnaou, Lasseux and Ahmadi2016). An asymmetric vortex pair is formed behind each solid obstacle (figure 16). Microscale flow symmetry is broken as a result of these instabilities, but it does not get translated into deviatory flow in the Reynolds-averaged macroscale flow field.

Table 4. The LES cases simulated to analyse the dependence of symmetry breaking on the solid-obstacle shape. The solid-obstacle shape is square cylinder for all of these cases.

In the case of a square-cylinder porous medium with low porosity (case C1), the flow separation and stagnation occur at the vertices of the square shape (figure 16b). Thus, the locations of the flow separation and stagnation are predetermined by the geometry. A low value of porosity is required for the formation of a lateral favourable pressure gradient. The microvortices formed in this case are characterized by slow turnover and resemble the flow in an open cavity flow (figure 16b). The primary and secondary flow regions are distinct. The primary and secondary flow regions are separated by a strong shear layer (figure 16a), which ‘channels’ the primary flow with little intervention from the secondary flow. The secondary flow region is confined between the solid obstacles by the primary flow. A sufficient magnitude of lateral favourable pressure gradient is not observed when the porosity is increased (case C2) or when the Reynolds number is increased (case C3). The solid-obstacle shape must possess a greater degree of circularity for the symmetry-breaking phenomenon to occur. The occurrence of deviatory flow for more solid-obstacle geometries is shown in the supplementary material.

Figure 15. The components of the macroscale RST versus Reynolds number for φ = 0.5.

Figure 16. (a) Turbulent structures visualized using the two-dimensional contours of vorticity magnitude for a porous medium with square-cylinder obstacles and porosity φ = 0.5 (case C1). A sub-volume of dimensions (3s, 2s, 2s) of the REV-T is shown. (b) Instantaneous flow streamlines projected on the xy plane for a sub-volume of dimensions (s, s, 2s). The streamlines are plotted on top of contours of static pressure.

Figure 17. The convergence of (a) the applied momentum source and (b) the macroscale turbulence kinetic energy (TKE) with increase in the REV size.

Figure 18. Turbulence energy spectrum for LES cases (a) φ = 0.50, (b) φ = 0.61, (c) φ = 0.72 and (d) φ = 0.80 at Re_p = 1000. The red square data points correspond to Δx_max/s = 0.03, green triangles to Δx_max/s = 0.02, blue circles to Δx_max/s = 0.01 and the dashed line to the −5/3 slope on the log plot.

Figure 19. Turbulence energy spectrum for DNS cases with Incompact3d (a) φ = 0.50, (b) φ = 0.61, (c) φ = 0.72 and (d) φ = 0.80 at Re_p = 1000. The red square data points correspond to Δx_max/s = 0.03, green triangles to Δx_max/s = 0.02, blue circles to Δx_max/s = 0.01 and the dashed line to the −5/3 slope on the log plot.