Simulations of stably stratified flow past two spheres at Re = 300

Intricate perturbations induced by a configuration of two spheres immersed in a stratified flow are representative of patterns typical of stratified flows past multiple generic obstacles. The study of such systems is relevant to geophysical flows, such as airflows over complex terrain and the associated mesoscale phenomena. Examples include flows past mountain ranges and the weather conditions within a valley, resulting from the interaction between the large-scale flow and flow driven by local conditions.¹ Features such as lee waves, dividing streamline height, low-pressure zones behind obstacles, and the drag coefficient are affected by the valley's width, driving both constructive and destructive interferences.² A sphere is traditionally used as a canonical shape in experimental investigations to gain insights into the physical processes common to atmospheric flows past orography. Flow structures past a single sphere have strong commonalities with flows past isolated hills, as can be observed, for example, by comparing simulations in Refs. 3 or 4 with those reported in Refs. 5–7. This interest can be extended to experiments involving configurations of multiple obstacles relevant, for example, to the design of wind turbine farms or to oceanographic applications, where spheres can represent archipelagos in different configurations, flows past complex bathymetry,⁸ and marine appliances submersed in thermoclines.⁹ 

Flow structures developing in stratified flows past a single sphere have been explored both experimentally^10,11 and computationally, e.g., Refs. 12 and 13. Numerical studies of particular interest here, which explicitly represent a sphere in a computational domain using body fitted meshes, include Refs. 3, 4, and 14 and references therein. Flow patterns past an isolated sphere have been observed and classified for wide ranges of Froude, $Fr=Vo/Nh$⁠, and Reynolds, $Re=VoL/ν$⁠, numbers; V_o, N, h, L, and ν represent free stream velocity, buoyancy frequency, obstacle height, characteristic length scale, and kinematic viscosity, respectively. Stratified flows produce a gravity-induced restorative buoyancy force, as fluid parcels are displaced from their hydrostatically balanced equilibrium position.⁴ At moderate Reynolds numbers, the three main distinct flow patterns, identified in Ref. 10, are a non-axisymmetric attached vortex, lee-wave instability regime, and two-dimensional vortex shedding. These regimes are briefly reviewed below, in preparation for the discussion of results presented in Secs. III and IV.

At Reynolds number Re = 300, a neutrally stratified flow past a sphere generates periodic shedding consisting of hairpin vortices¹⁵ with time-independent planar symmetry defined by the initially induced vortex. For very low stratification (i.e., $Fr≥10$⁠), the characteristic structures of the neutrally stratified flow are retained; however, the plane of symmetry moves to the horizontal plane coinciding with the sphere's center. As stratification increases, buoyancy forcing flattens the wake and inhibits vortex shedding—bringing the flow to a steady state. The resulting pattern is characterized by a double wake structure, deformed in the vertical doughnut-shaped vortex ring in the lee of the sphere, and gravity wave radiation.¹⁰ This flow regime is called a non-axisymmetric attached vortex. Hereafter, we adopt the nomenclature of flow patterns reported in Ref. 10.

The peak of gravity wave amplitude occurs at $Fr≈1$⁠, corresponding to the establishment of the lee-wave instability pattern. At this level of stratification, fluid parcels below the dividing streamline height (cf. Ref. 15) are forced to pass around the sides of the obstacle rather than over it, developing into a horizontal bi-dimensional recirculation behind the sphere. Fluid parcels above the dividing streamline height contribute to the production of gravity waves in the lee of the sphere. Together, the two flow structures constitute the lee-wave instability regime.

For approximately Fr < 0.56, the flow pattern changes into two-dimensional vortex shedding set in the horizontal plane. The characteristic features of this regime are fluid parcels traveling around the sphere, producing a periodic motion similar to that observed for a flow past a cylinder. In numerical simulations, the two-dimensional vortex shedding regime is slow to develop; therefore, a sufficiently long computational time is required for the regime to become correctly established, cf. Refs. 3 and 16.

The present work builds on our earlier studies of linearly stratified flows past a sphere at Re = 200 and Re = 300 (Refs. 3 and 15) and extends them to flows past two spheres for a range of Froude numbers $Fr∈[0.25,∞]$⁠. A moderate Reynolds number (Re = 300) is selected in order to avoid uncertainty related to turbulence modeling and to facilitate comparisons with cases of neutrally stratified flows. Furthermore, historically the experimental data for stratified flows past a sphere also concentrated on low Reynolds numbers due to ease of measurement. However, our study in Ref. 3 and the study in Ref. 4 show that as stratification increases, the key flow patterns depend on viscosity only slightly. It is therefore likely that the presented results are also relevant to flows with higher Reynolds numbers. The study reported here considers configurations of two identical spheres: obtained by varying the distance between the spheres and their orientation relative to the free stream direction.

Studies of neutrally stratified flows for a variety of geometrical arrangements of two spheres have been documented in, e.g., Refs. 17–22, with Ref. 23 extending the study to a flow past in-line spheres. However, to the best of our knowledge, this is the first comprehensive investigation of stratified flows past two spheres. We present a systematic characterization of flow patterns, complemented by a quantitative analysis of drag coefficients, separation angles, and vortex-shedding frequencies. Simulations suggest that vortex characteristics vary significantly with the relative placement of the spheres when in close proximity to each other. This is particularly pronounced for weakly stratified flows, whereas at lower Froude numbers the flow patterns are dominated by stratification. As stratification increases, buoyancy effects progressively overcome viscosity; however, the effects of vortex interactions and upstream objects blocking the incoming flow remain present.

Computations adopt a non-oscillatory forward-in-time finite volume (NFT-FV) approach which implements the Message Passing Interface (MPI) standard for parallelisation.^3,15 The numerical scheme integrates the Navier–Stokes equations under gravity in the incompressible Boussinesq limit. Finite volume discretization is implemented on flexible unstructured/prismatic meshes that are particularly well suited to represent various configurations of two spheres and to capture details of flow patterns. The validity of the NFT-FV scheme operating on unstructured meshes has been successfully demonstrated for a range of stratified flows past orography^5,6,24 and a single sphere,^3,15 in addition to global simulations of the Earth's atmosphere employing meshes with a combination of a quasi-uniform octahedral distribution in the horizontal and prismatic spacing in the vertical.²⁵ The present work provides additional comparisons with numerical and experimental results available in the literature^17,20 for two spheres immersed in a neutrally stratified flow, where we illustrate the scheme's efficacy together with its ability to reproduce flow intricacies generated by the presence of multiple objects. Direct numerical simulations (DNSs) of stably stratified flows for Froude numbers $Fr∈[0.25,∞]$ past two spheres are also contrasted with the previously documented computations for a single sphere.³ 

The reminder of the paper is organized as follows: Sec. II outlines the governing equations and the numerical approach employed. Section III illustrates the performance of the NFT-FV model for neutrally stratified flows past two spheres. Section IV documents a systematic numerical study of stably stratified flows past two spheres in representative geometrical configurations at Re = 300. Remarks in Sec. V conclude the paper.

The nonhydrostatic NFT-FV model used in this study solves the Lipps–Hemler^26,27 anelastic system, which describes a broad class of mesoscale atmospheric flows.^5,28,29 For a dry atmosphere, the system can be represented by the following set of conservation equations for mass, momentum, and entropy fluctuations:

where $xI(I=1,2,3)$ are the Cartesian coordinate directions, V_I are the corresponding components of velocity $V$⁠, g is the gravitational acceleration, ρ represents density, δ_IJ is the Kronecker delta, and $τ$ signifies the deviatoric stress tensor. Potential temperature is symbolized as θ and is linked to the ideal gas specific entropy, s, by $ds=cp d ln θ$ (c_p is the specific heat at constant pressure). The normalized pressure perturbation, $φ′$⁠, in Eq. (2) is defined as $φ′=(p−pe)/ρ¯$⁠. The density and potential temperature of the static reference state are $ρ¯$ and $θ¯$⁠, respectively, and subscripts “e” and “o” indicate the stably stratified ambient state and constant reference values, respectively. Similarly, the potential temperature perturbation is written as $θ′=θ−θe$⁠, where the ambient profile, $θe(x3)=θo exp (Sx3)$⁠, assumes constant stratification, S. Components of deviatoric stress tensor, $τ$⁠, are defined as $τIJ=μ(∂VI∂xJ+∂VJ∂xI)$⁠, which incorporates the dynamic viscosity, $μ=ρ¯ν$⁠. The last terms in Eqs. (2) and (3) activate wave-absorbers in the proximity of domain boundaries, with α and $α′$ being attenuation factors for the velocity components and the potential temperature perturbation toward the ambient state.³⁰ 

Simulations undertaken herein focus on gravity wave dynamics for low speed atmospheric and hydrodynamic applications, for which displacements of fluid parcels are small in comparison to the scale height, $S−1$⁠, of the stratified environment. For such fluids, the incompressible Boussinnesq limit of the anelastic system (1)–(3) is assumed, such that $ρ¯=ρo, θ¯=θo$⁠, and $θe(x3)=θo+Sx3$⁠.

Solutions are obtained using the high resolution NFT-FV model, detailed in Ref. 15. The model exploits non-oscillatory forward-in-time integration, derived from a truncation-error analysis of uncentred two-time level approximations of an archetype inhomogeneous partial differential equation (PDE) for fluids.³⁰ This semi-implicit scheme is second-order accurate and builds on the two-time-level nonlinear advection algorithm multidimensional positive definite advection transport algorithm (MPDATA) and the corresponding solution to the elliptic Poisson equation for the pressure perturbation. Derivation of MPDATA for unstructured meshes was documented in Refs. 31 and 32. By design, MPDATA relies on an iterative application of the upwind scheme and is genuinely multidimensional and free of splitting errors. It is also conservative and sign preserving, while its monotonicity is ensured in the spirit of Flux Corrected Transport methods. A non-symmetric Krylov-subspace solver (detailed in Ref. 28) is employed to solve the Poisson equation implied by the soundproof mass-continuity constraint.

The finite volume spatial integration of the NFT-FV model adopts edge-based connectivity and operates on arbitrary shaped polyhedral median-dual cells,⁶ constructed from a given primary unstructured mesh. Such meshes are particularly suitable for simulations involving complex and multi-object geometries, and for ensuring adequate mesh resolution needed to capture key flow features, such as evolution of boundary layers and wake structure. They also facilitate MPDATA driven options for mesh adaptivity.^6,33

The implementation of the NFT-FV scheme for unstructured meshes embraces the distributed memory paradigm,³ where efficient partitioning of the computational domain is accomplished by the MeTis library.³⁴ An in-house pre-processing code configures layers of single halo points to perform MPI point-to-point communications between neighboring partitions. These exchanges are the base of NFT-FV code parallelization, with limited implementation of MPI global communications within the elliptic solver. The model uses a collocated arrangement for all depended variables, which simplifies the discretization of differential operators, and has low memory and communication requirements.

We consider flows past two identical spheres with diameter D = 1 located within a cuboidal computational domain. The size of the domain is dependent on the arrangement of spheres. It is chosen to maintain a minimum distance (15D) between any sphere and the external boundaries,¹⁷ to be large enough to adequately simulate the far-wake field (cf. Ref. 3 and references therein). Geometrical configurations of the spheres adopted here follow those reported in Ref. 17 for neutrally stratified flows. The first sphere is always placed at the origin of coordinate system (0, 0, 0) and the second sphere's position is varied between simulations. For neutrally stratified flows, each geometrical configuration, illustrated in Fig. 1, is described by two parameters: the normalized gap between spheres, l/D, measured along the line connecting spheres' centers and the angle, $θ∈[0°,90°]$⁠, this line forms with the direction of the free stream flow. In this study, with the exception of results in Sec. IV D where we introduce a configuration tilted in the vertical plane, we focus on arrangements with mutual locations of spheres changed only in the horizontal plane that are most relevant to applications. The domain is therefore defined as $[−15D,25D+(D+l) cos θ]×[−15D,15D+(D+l) sin θ]×[−15D,15D]$⁠.

The spheres are immersed in a flow with uniform free stream velocity, $Ve=Vo=(1,0,0)$ ms⁻¹, prescribed at the streamwise boundaries of the domain. No-slip boundary conditions are imposed on the solid surfaces and free-slip conditions on spanwise faces of the domain. Absorbers are applied up to 2D from the streamwise boundaries to attenuate the solution toward ambient flow conditions, with the inverse time scales α and $α′$ in Eqs. (2) and (3) increasing linearly from zero at 2D from the domain boundary to $1/150$ s⁻¹ at the boundary. All simulations are initialized by the solution of potential flow and then run for the non-dimensional time of $T=tVo/D=400$⁠. The Reynolds and Froude numbers for a sphere with constant diameter become $Re=VoD/ν$ and $Fr=2Vo/ND$⁠, respectively. The buoyancy frequency is equal to $N=gS$⁠, where S is the constant stratification introduced in Sec. II.

The primary computational meshes, illustrated in Fig. 2, consist of tetrahedral elements everywhere except for 10 layers of prismatic elements with a triangular base built within $0.2D$ thick regions from the surfaces of the spheres. The prismatic layers vary in thickness, the thinnest (⁠$0.005D$⁠) of which are in the layer closest to the spheres. Nodes are concentrated in the near sphere regions, inter-sphere region, and in the wakes. Mesh resolution increases away from the sphere surfaces, reaching a spacing of $∼1D$ near to the external boundaries. Each geometric configuration has an individually generated mesh, where the number of nodes ranges between 1 192 456 and 1 656 642, with an average of $∼1.5$ × 10⁶ nodes. Simulations are run on 5 nodes of the Lovelace HPC cluster (http://hpc-support.lboro.ac.uk/), each of them equipped with 40 Intel Xeon Gold 6248 CPUs. The computational domain is divided into 200 partitions and the average run-time of simulations is 18.5 h.

The number of prismatic layers adopted results from a mesh sensitivity study performed for a stratified flow past two spheres placed in configuration $θ=0°$ and $l/D=1$ at Fr = 1.667. Two additional meshes have been defined: a coarser one with 7 prismatic layers, the thinnest of which has a thickness of $0.007D$⁠, and a finer mesh with 13 prismatic layers with the minimum thickness of $0.004D$⁠. The two spheres' drag coefficients (C_d) were evaluated after T = 200 and compared to the reference solution obtained on the computational mesh described above with 10 prismatic layers. Drag coefficients resulting from simulations employing the mesh with finer prismatic layers are similar to those computed from the current mesh, having errors relative to the reference solution $<0.01%$ and $<0.1%$ for the first and second sphere, respectively. For the mesh with coarser prismatic layers, the relative errors in the drag coefficients amount to 0.3% for both spheres.

To test the robustness of the NFT-FV scheme, white noise of amplitude $0.2Vo$ was introduced while initializing the simulation of a stratified flow past two spheres placed in the previous tandem configuration at Fr = 1.667. The test was then repeated with the noise amplitude increased to $2Vo$⁠. Simulations were conducted on the computational mesh with 10 prismatic layers and resulting flows at T = 200 were compared with those from the non-perturbed simulation. The difference between perturbed and non-perturbed results for velocity components and potential temperature perturbation is evaluated over the entire domain and their statistics are illustrated in Table I. The low average values and standard deviations demonstrate the robustness of NFT-FV simulations which are generally run for a non-dimensional time of T = 400.

Performance of the NFT-FV model is illustrated first for neutrally stratified flows. For computations presented in this section, the anelastic system (1)–(3) is reduced to the Navier–Stokes equations for incompressible isothermal flows, by omitting the entropy equation (3) and the buoyancy term in the momentum equation (2). We documented simulations of tandem configurations (i.e., $θ=0°$⁠) earlier in Ref. 3, whereas here parallel configurations of spheres (i.e., $θ=90°$⁠) with $l/D∈[0.41,2.5]$ are selected to facilitate comparisons with the existing numerical and experimental data. Interaction between the two wakes occurs only when the spheres are in sufficient proximity and depends on the gap between them.²⁰ For the selected cases, the near-wake fields (i.e., less than 1D from the spheres' surfaces) are individually recognizable, as shown in Fig. 3. To aid vortex identification, Fig. 4 displays the corresponding Q-method,³⁵ i.e., an instantaneous level set of the positive second invariant of the deformation tensor; $Q=0.5(ΩIJΩIJ−dIJdIJ)$⁠, where $dIJ$ is half of the tensorial part of $τIJ$ and $ΩIJ=0.5(∂VI∂xJ−∂VJ∂xI)$ are entries of the rotation tensor.

When the spheres are closely spaced (e.g., $l/D=0.41$⁠, Figs. 3 and 4 upper panels), the wakes interact strongly, starting in the near wake region, merging further downstream into a non-symmetric unified wake in the horizontal plane as the vortex loops shed by each sphere become entangled. Figure 4 top-right panel demonstrates hairpin vortex shedding in the vertical plane, strongly resembling the wake induced by a single bluff body, with symmetry about the z = 0 plane. The resulting vortical structures are consistent with the experimental visualization presented in Figs. 4(a)–4(c) in Ref. 20.

For $l/D=1$⁠, Figs. 3 and 4 middle panels, each sphere produces an identifiable individual wake in the lee, encompassing a single recirculation close to the solid surface followed by a subsequent secondary recirculation formed behind it (see Fig. 3 middle panel). The periodic shedding of vortex loops occurs simultaneously from each sphere. The loops combine downstream to form a wake almost symmetric both in the vertical and horizontal planes (cf. Fig. 4 middle panels, right and left, respectively). The simulated pattern closely matches experimental observations provided by Fig. 4(d) in Ref. 20.

For $l/D=2$⁠, Figs. 3 and 4 lower panels, similar symmetry patterns and vortical structures emerge, however with different wake size and shedding frequency. Additionally, a single recirculation is visible in the horizontal central plane, Fig. 3, behind each sphere. In contrast, two symmetrical recirculations form in this region behind a single sphere.

For all presented cases, a history of the drag coefficient (not shown) generated by each sphere demonstrates periodic oscillations, attaining the same amplitude for both spheres. For the closest arrangements of spheres (⁠$l/D<1$⁠), they are about half a period out of phase, while when the spheres are further apart (⁠$l/D≥1.5$⁠), the histories of the oscillations for both spheres become identical.

Figure 5 compares each sphere's normalized frequency of vortex-shedding, obtained using the NFT-FV model, with those reported in both experimental²⁰ and numerical^17,20 studies. The vortex-shedding frequencies, f, are derived from the frequency spectra of the drag coefficient and are normalized by the frequency obtained for a single sphere fs = 0.138; frequencies are computed using TSA v1.2 package of R software, https://www.rdocumentation.org/packages/TSA. The NFT-FV simulations agree well with previous studies within the full range of selected distances, l/D.

For completeness, Table II documents the computed values of mean drag coefficient (⁠$C¯d$⁠) evaluated within the non-dimensional time interval $T=300–400$⁠. The corresponding amplitude of oscillations (⁠$C′d$⁠), Strouhal numbers computed from the fundamental frequencies (⁠$St=fD/Vo$⁠), and separation angles (⁠$Φs$⁠) for selected simulations with $l/D∈[0.41,2.5]$ are also provided. When the two spheres are in the parallel configuration, both spheres have matching values.

Configurations of two spheres with surface-to-surface distances $l/D∈[0.5,1.5]$ and angles $θ∈[0°,90°]$ are selected for simulations at three stratification levels, $Fr=1.667, 0.625, 0.25$⁠, representative of the patterns identified in the characterization of stratified flows past a single sphere discussed in Sec. I. For very low stratification, the patterns are reminiscent of those obtained for neutrally stratified flows; cf. Sec. III and Ref. 3. For higher levels of stratification, the computational results presented illustrate the most typical flow behavior for each given flow regime. Additionally, Subsection IV D provides examples of flows past tilted configurations, where the azimuthal-like angle (labelled $θ′$⁠) between the spheres' centers is in the vertical plane.

For weak stratification, features of the non-axisymmetric attached vortex regime distinguishing flow past a single sphere can also be observed when a second sphere is placed in the vicinity of the first one. However, simulations demonstrate that the proximity and relative orientation of the spheres determine the extent of wake interactions and deformations of vortices. Figure 6 illustrates a typical wake pattern obtained for the staggered configuration $θ=60°$ and $l/D=0.5$⁠. Buoyancy forces, resulting from the increase in stratification, suppress the hairpin-shape vortex shedding of the neutrally stratified flow until the flow becomes stationary, and both spheres generate double wakes in the lee while the formation of gravity waves begins, as can be observed from Fig. 6.

Figure 7 shows streamlines and contours of the streamwise velocity component in the horizontal plane, z = 0, for parallel (upper panel), staggered (middle panel), and tandem (lower panel) configurations of the spheres. The upper panel shows simulation results with the spheres in close proximity (⁠$l/D=0.5$⁠), for which significant interaction is observed. The interaction between wakes most noticeably results in their mutual flattening in the region of contact between the flows. Indeed, the interacting wake branches are most affected: their thickness is reduced, especially in the near-wake field. As a result, symmetry between the two branches of the individual wakes of each sphere is broken; however, symmetry of the entire flow about the vertical plane passing through the midpoint of the gap between spheres is present. Behind the spheres, in y–z plane flow travels toward the central plane z = 0 and forms two doughnut-shaped vortex rings deformed in the vertical due to stratification. These vortices, along with velocity streamwise component contours, are displayed in Fig. 8. These flows feed two recirculations that form behind the spheres in the x–z plane, not shown. For parallel configurations with distances between the spheres increased to $l/D=1$ and $l/D=1.5$⁠, interactions in the far-wake field are still observed while the near-wake field of each sphere more closely resembles that of a single sphere.

For the staggered configurations (e.g., Fig. 7 middle panel), similarities between the flow behavior behind the first sphere and that of a single sphere progressively reduce with decreasing θ. For angles $θ≥45°$⁠, the majority of interactions occur in the far-wake field for $l/D=1$ and $l/D=1.5$⁠, while for $l/D=0.5$ interactions also occur in the near-wake field affecting the shape of the internal branch of first sphere's wake, as shown in the left panel of Fig. 6. This branch is almost completely (⁠$l/D=0.5$⁠) or partially (⁠$l/D=1$⁠) suppressed for $θ=30°$⁠. The middle panel of Fig. 7 shows the latter case, in which the first sphere's wake is deflected allowing for a single recirculation bubble to form, while the flow pattern behind the second sphere is similar to the patterns observed in the $l/D=0.5$ parallel configuration case discussed earlier. Finally, when the spheres are in a tandem configuration (Fig. 7 lower panel), the flow is symmetric about both the horizontal and vertical planes. Wake branches from the first sphere encompass the second sphere. The flow behind the first sphere is suppressed, forcing a flow in the y–z plane and the subsequent formation of a vortex in the region between the spheres. Its cross section (not shown) just behind the first sphere reveals a doubly reconnected doughnut-shaped vortex, resembling that identified for a single sphere in Ref. 10 for a slightly lower Fr, with two attached vortices located horizontally indicating that the presence of the second sphere adds to the deformation induced by the stratification. However, cross sections taken close to the second sphere show the influence of the first sphere on the stagnation region, which deforms the doughnut-shaped vortex such that the two attached vortices are now placed vertically. The second sphere is shielded from the incoming free stream flow. This decreases the velocity in the gap, and consequently the local Froude and Reynolds numbers. Therefore, the flow in the gap is further suppressed by the stratification and a dividing streamline is formed, below which flow parcels are forced to travel around the second sphere. Two symmetrical recirculations can be observed in its lee in the horizontal plane, reminiscent of those typically formed in the lee-wave instability regime.

Table III reports the drag coefficients from the first and second spheres and the corresponding separation angles. Two values of separation angle are defined on the horizontal plane to address the non-symmetric geometry of the test cases. $Φs1$ is measured clockwise from the upstream centerline of each sphere to the nearest separation point; $Φs2$ is measured anticlockwise from the same direction to the nearest separation point placed on the opposite side of the sphere. When needed, the additional subscripts $,1$ and $,2$ are introduced to denote the first and the second spheres, respectively. Figure 9 illustrates $Φs1$ and $Φs2$ for a single sphere placed in the center of the domain. At the same level of stratification, the drag coefficient for a single sphere amounts to $Cd=0.688$³ and the horizontal separation angle is $Φs1=Φs2=108°$⁠. The data are shown for all simulated geometrical configurations. For the entire range of angles and distances, the drag coefficient computed from the first sphere has the tendency to increase with θ. For low angles (⁠$θ<45°$⁠) the drag coefficient for the first sphere grows as the gap between spheres increases. When the angle is increased to $θ=45°$ results for all selected l/D become very similar, while for $θ>45°$ the drag coefficient increases as the gap reduces. The second sphere's drag coefficient does not show a clear pattern, apart from the case of $l/D=1$ for which the drag decreases as the angle θ increases, until it reaches an almost constant value for $θ≥45°$⁠. This indicates that angle θ affects the flow solution behind second sphere more than the gap between spheres.

For all tandem configurations, the separation angle from the first sphere is always lower than that from the second one with a negligible difference between $Φs1$ and $Φs2$ for each sphere confirming the flow's symmetry. Notably, the separation angles computed for the first sphere nearly match that of a single sphere. For staggered and parallel configurations, the previously described flattening of internal wake branches results in general in $Φs,11>Φs,12$ and $Φs,22>Φs,21$⁠. The flattening of wake branches facilitates the movement of outer separation points toward the aft of the spheres leading the separation angles $Φs,11$ and $Φs,22$ to be generally greater than or equal to the single sphere's value; exceptions are $Φs,11$ for configuration $θ=30°$ and $l/D=1.5$ and $Φs,22$ for configuration $θ=60°$ and $l/D=0.5$⁠. For some configurations of close proximity spheres (⁠$l/D=0.5$⁠), the spread between values of $Φs1$ and $Φs2$ (for either sphere) is larger. In such cases, a recirculation bubble is formed behind the first sphere, similar to that seen in Fig. 7 middle panel, but with its center shifted toward the sphere's centerline, leading to an increased separation angle, $Φs,12$⁠. In the parallel configuration (Fig. 7 upper panel), the flattening of wake branches is associated with $Φs,11∼Φs,22$ and $Φs,12∼Φs,21$⁠.

At medium stratification, flow impinging on a single sphere results in a so-called steady-state lee-wave instability regime. The resulting flow for two staggered spheres is illustrated using the Q-method in Fig. 10. The flows behind each sphere reach steady states that have similar features to those reported in Ref. 3 for a single sphere. Well developed gravity waves can be seen in the lee, strongly interacting with the wakes from both spheres.

The streamwise velocity component plots with associated streamlines in Fig. 11 illustrate the effects of increased stratification for configurations similar to those examined in Fig. 7. In all cases fluid parcels contained below the height of dividing streamline, as discussed in Ref. 15, flow around the spheres forming leeward recirculation bubbles. For the parallel configuration (Fig. 11 upper panel), the pattern is almost symmetric about the vertical plane of mesh symmetry between the two spheres (i.e., y = 0.75) and only one recirculation is visible in the horizontal (z = 0) plane behind each sphere, with the second one being suppressed. Aided by the recirculations, the flow accelerates through the gap between spheres to create a region of locally increased velocity which soon decreases further downstream.

The symmetry noted for the parallel configuration is not present for staggered configurations (Fig. 11 middle panel). The two attached recirculations formed behind each sphere are distorted by the deflection and interaction between the wakes. In the tandem configuration (Fig. 11 lower panel), planar symmetry is present and the two recirculations behind the second sphere are similar to those obtained for a single sphere at Fr = 0.625; however, the recirculations behind the first sphere are suppressed and, therefore, much smaller.

Lee-wave crest instability and overturning are induced by a strong vertical shear. Tiny recirculations, representing overturning motions of fluid parcels, are illustrated for the tandem configuration in the vertical plane (y = 0) in Fig. 12. These recirculations result from the combined effects of buoyancy and shear,¹⁰ leading to the name “lee-wave instability regime.” Only flow parcels above the dividing streamline have sufficient kinetic energy to flow over the top of the spheres. The associated buoyancy-induced gravity waves radiating in the lee of both spheres are presented in the vorticity plot in Fig. 13 for a tandem configuration, along with the results achieved for a single sphere. Gravity waves are shown to propagate in proximity to the first sphere; however, the wave amplitudes start to differ downstream of the second sphere. In Fig. 13 upper panel, regions of high and low amplitudes of the absolute value of vorticity display positions of constructive and destructive interferences, respectively. This interference pattern develops conically immediately behind the spheres, with a cone angle of $∼55°$ from the flow direction. In this region the first destructive interference is noticeable, until eventually the waves dissipate downstream. Similar interference patterns are also observed for staggered configurations (not shown). The superposition pattern is induced by the presence of the second sphere, while lee-wave behavior is more consistent with the pattern of a single sphere. The geometry of the configuration drives the interference influencing frequency, maximum vorticity amplitude, and the shape of the resulting lee waves.

Finally, we record that the flow solution for the tandem configuration and distance $l/D=1$ is a special case, which displays periodic horizontal oscillations in the far-wake field, as shown in Fig. 14. The same behavior is not noted for other values of l/D, where steady state is maintained (e.g., configuration $θ=0°$ and $l/D=0.5$ in the lower panel of Fig. 11). Notably, for neutrally stratified flows at Re = 250, varying the distance between spheres placed in a tandem arrangement also displays a range of flow regimes: from planar steady state to periodic flow.³⁶ 

Lee waves trapped between ridges have been studied in Ref. 2 for elongated obstacles positioned perpendicular to the flow direction. Despite the differences between simulations of stratified flow past spheres and atmospheric flows over mountain ridges, the number of lee-wave crests across the gap between obstacles displays a resemblance. Present results agree with the observation that trapped lee waves increase in number as the gap between obstacles increases.

Table IV documents drag coefficients from the first and second spheres, along with the corresponding separation angles. For reference, at the same stratification level the single sphere value of drag coefficient is $Cd=1.230$ (Ref. 3) and separation angle is $Φs1=Φs2=144°$⁠. Similarly as for a low stratification described above, the presence of a second sphere can affect these values. For $θ=0°$⁠, the second sphere is partially shielded from the incoming flow by the first one; therefore, the drag coefficient of the second sphere is substantially smaller, but tends to grow with the increase in gap between spheres. As the angle increases to $θ=30°$⁠, the drag coefficient slightly increases for the first sphere and dramatically for the second which is now more exposed to the incoming flow. This growth is more pronounced as the gap between spheres increases. In fact, for $l/D=1.5$ the drag coefficient computed for the second sphere is even higher than for the first one. An increase in the second sphere's drag coefficient continues within the interval $30°<θ≤60°$⁠. Although the rate of increase is lower than rates observed for the lower values of θ, the drag from the second sphere is larger than that of the first one, even for the distance $l/D=1$⁠. When $θ=60°$⁠, values of the second sphere's drag coefficient are (almost) the same for all considered gaps between spheres. Finally, when spheres are placed in the parallel configuration, a difference between drag coefficients computed for both spheres is less than 1%, suggesting that their mutual interference is weak.

The separation angle $Φs,11$ remains close to $Φs,22$⁠, and $Φs,12$ has similar values to $Φs,21$ for parallel configuration of spheres, alike that observed for low stratification. Moreover, the pattern of separation angles for which $Φs,11$ and $Φs,22$⁠, respectively, are higher than $Φs,12$ and $Φs,21$ is still apparent when $θ=90°$⁠. These observations can only be partially extended to the staggered configurations. In fact, the same trend in separation angles is observed for the first sphere in the majority of cases, i.e., $Φs,11>Φs,12$⁠, while for the second sphere, there is no clear tendency. For the tandem configurations, values of $Φs1$ and $Φs2$ differ by less than $3°$⁠. This suggests that equal portions of flow circumnavigate the spheres in the two possible directions. Furthermore, for $θ=0°$ the separation angles for the first sphere are greater than those for the second one and all of them are lower than the separation angle for a single sphere. For staggered configurations in which the spheres are close to each other, e.g., $θ=30°$ and $l/D=0.5$ (reported in the middle panel of Fig. 11), the second sphere can force the closest portion of the flow which circumnavigates the first sphere to delay the flow separation leading to an increase in the corresponding separation angle; in this particular case its value is $Φs,11=175°$⁠. The flow deflection is associated with the low separation angle measured on the second sphere, i.e., $Φs,22=87°$⁠.

A flow past a single sphere at Fr = 0.25 is characterized by the well-documented two-dimensional vortex shedding regime. Figure 15 shows a three-dimensional representation of the flow past two spheres placed in a staggered configuration with $θ=60°$ and $l/D=0.5$⁠. Shapes of vortex structures in the near-wake field of each sphere are comparable with those for a single sphere reported in Ref. 3. The far-wake field however, although still periodic, differs greatly due to the presence of the second sphere. The interaction between the wakes of each sphere also results in different frequencies of lee waves and vortex shedding (Fig. 15 right panel) in comparison to flows past a single sphere. The corresponding Strouhal numbers are listed in Table V and will be discussed later. The definition of Strouhal number (cf. Sec. IV) is applied adopting the first two frequencies resulting from the spectrum of drag coefficient] are listed in Table V and will be discussed later.

The generation and development of vortices that are shed behind the two spheres depends strongly on their configuration. Periodic behavior for the tandem configuration with $l/D=0.5$ is illustrated in Fig. 16, displaying instantaneous streamlines and streamwise velocity component contours. The figure also shows how extended shear layers, generated by the first sphere, encapsulate the near-wake field inhibiting vortex development. This is also reflected in the almost constant values of drag coefficients, computed for both spheres, with the maximum amplitude of drag coefficients over their periodic histories not exceeding 1% of $C¯d$⁠. The described flow pattern occurs for all simulated tandem configurations at Fr = 0.25.

Figure 17 presents streamlines and contours of streamwise velocity component in the horizontal plane (z = 0) for a further three representative cases, organized according to angle θ. For staggered configurations (Fig. 17 top and middle panels), some periodicity of the flow in the near-wake field is restored. Figure 17 upper panel (configuration with $θ=30°$ and $l/D=0.5$⁠) demonstrates that the merging of wakes in the near-wake field determines the same Strouhal number for both spheres. Altering the staggered configuration such that the second sphere is no longer directly behind the first (e.g., $θ=60°$ and $l/D=1$⁠, Fig. 17 middle panel), the presence of shear layers in the gap and separate recirculations attached to the spheres trigger wake unification in the far-wake field. The pattern in this case results in multiple vortex shedding frequencies. Table V displays the Strouhal numbers computed from the two highest frequency peaks for the first sphere, while those for the second sphere are not readable from the spectral analysis. Figure 17 lower panel shows the flow for a parallel configuration with $l/D=1.5$⁠. The wakes remain separate in the near-wake field, evoking the same values of drag coefficient and Strouhal numbers for both spheres. They merge gently further downstream.

Figure 18, displaying the y-component of instantaneous vorticity in central vertical plane (y = 0), shows characteristics of the lee-wave pattern demonstrating constructive and destructive interferences. Strong stratification facilitates a pronounced lee-wave propagation upstream and a fast wave dissipation downstream. Furthermore, both a higher frequency of lee waves and reduced area of the interference pattern are observed in comparison to the medium stratification case in the upper panel of Fig. 13. Moreover, the displayed vorticity component shows periodic features in the far-wake field which are symmetric about the centerline of wake. These represent the downstream vortices, also visible in the three-dimensional Q-plots in Fig. 18, for a staggered configuration. Such wake features are also observed for a flow past a single sphere at the same stratification.³ 

Table V provides quantitative data obtained from the simulations. The mean drag coefficient $C¯d$ was evaluated within the non-dimensional time interval $T=300–400$ with the corresponding amplitude, $C′d$⁠, while the Strouhal numbers, St, are computed based on the two dominant frequencies of each drag coefficient spectrum. Additionally, the table provides the values of separation angles $Φs$ (previously defined). For reference, the values for a flow past a single sphere at Fr = 0.25 are $C¯d=1.617$ with amplitude $C′d=0.01$⁠, St = 0.375, and $Φs1≈Φs2≈99°$⁠. The overall patterns of $C¯d$ for both spheres are similar to those obtained for medium stratification (Fr = 0.625); however, there are some notable differences. In particular, the growth of drag coefficient computed for the first sphere for increasing θ is now not as smooth. Values of $C¯d$ remain constant for all l/D when $θ≤30°$ and then increase for $45°$⁠. Further increases of $C¯d$ are moderate with the exception of $l/D=1$ case, for which $C¯d$ increases noticeably for $θ≥60°$⁠. Moreover, for the tandem and parallel configurations, the first sphere's drag coefficients, normalized by $C¯d$ from the single sphere, have similar values to the normalized drag coefficients computed for medium stratification (Fr = 0.625). This demonstrates their relative independence from stratification. Drag coefficients computed from the second sphere have negative values for all tandem configurations, showing that the established flow induces a forcing which would tow the object toward the first sphere. After increasing angle θ, the second sphere's drag coefficient approaches a value similar to that computed for a single sphere and exceeds it at approximately $θ=60°$⁠, for all considered distances l/D, with a tendency to increase as the gap increases.

Interpretation of resulting patterns is aided by the Strouhal number, computed from the history of $C¯d$⁠; therefore, we evaluate St tendencies displayed in Table V further. For tandem configurations, the flow reaches steady state; hence, the Strouhal number is not computed. When the angle between spheres is increased, for $0<θ≤45°$⁠, a single dominant frequency is observed for the first sphere and multiples of this single dominant frequency are also detectable for the second sphere. The values of St₁ are the same for both spheres, outlining the unification of flow patterns (as previously discussed). The configuration $θ=45°$ and $l/D=1.5$ shows features which are characteristic for higher angles between spheres, middle panel of Fig. 17, consisting of multiple vortex shedding frequencies behind at least one sphere; however, the same values of St₁ are found for each of the two spheres. This indicates that the vortex shedding follows an overall fundamental period. Some exceptions are noted, especially for the parallel configuration, in which clear predominant frequencies cannot be detected, since the $C¯d$ harmonics are highly superimposed. However, when the two spheres are widely separated (e.g., $l/D=1.5$⁠, lower panel of Fig. 17), the independence of flows behind each sphere results in the same Strouhal numbers, as expected. Moreover, in general, for the same θ the Strouhal number increases with l/D; this pattern is more pronounced for lower values of θ.

Table V also shows that separation angles for strong stratification are lower in comparison with those for medium stratification, where their values approach the maximum, cf. Ref. 10 for a single sphere. Separation angles for the parallel configuration $Φs,12$ and $Φs,21$ are close to the $Φs1$ and $Φs2$ values for a single sphere. For lower angles θ, computations from the first sphere display the tendency of $Φs,11$ to decrease with increased values of θ and l/D; in particular, as $θ↘0°$ values of $Φs,11$ become closer to that of a single sphere. Separation angles $Φs,12$ are more established than $Φs,11$ and their values are always close to that of a single sphere (although mostly slightly lower). On the other hand, computations of separation angles for the second sphere show greater variability than for the first one; such variability is observed especially for $Φs,22$⁠, since a trend of $Φs,21$ to increase with l/D is noticeable for some angles θ.

A further class of configurations can be established by introducing an azimuthal-like angle $θ′$⁠, defined in the x–z plane counterclockwise from the x-axis to the line connecting spheres' centers. This configuration is analogous to the picture in Fig. 1, replacing y and θ with z and $θ′$⁠, and the spheres' centers now remaining in the x–z plane. Accordingly, the domain is now defined as $[−15D,25D+(D+ l) cos θ′]×[−15D,15D]×[−15D,15D+(D+l) sin θ′]$ and contains 1 301 424 mesh nodes. The spheres appear directly behind each other in the horizontal, but are staggered in the vertical such that the second is either above or below the central horizontal plane (z = 0). To illustrate the flow resulting from these tilted configurations simulations of stratified flows past two spheres at Fr = 1.667, Fr = 0.625, and Fr = 0.25 are documented for configurations described by $θ′=30°$ and $l/D=0.5$⁠.

Figure 19 shows the Q-plots for weak stratification, i.e., Fr = 1.667. The characteristics of the non-axisymmetric attached vortex regime are present, with the double wake structure visible in the horizontal plane (left panel) and recirculations behind the spheres detected in the x–z plane (not shown). The doughnut-shaped vortices are still visible in the lee of both spheres; however, they are distorted in the vertical not only by stratification but also due to compression generated as the flow interacts with the spheres. The absence of the dividing streamline for weak stratification enables the entire flow impacting the first sphere to pass either above or underneath it, depending on the hemisphere encountered. The fluid parcels which rise over the top of the first sphere and then pass below the second one forming the shear layer between the spheres which merges with a shear layer inherently generated by the second sphere. The tilted configuration with $θ′=30°$ facilitates the constructive interference of lee waves generated from the two spheres, as the lee waves in Fig. 19 are more evident than those in Fig. 6 especially in the far-wake field. Drag coefficients of first and second spheres placed in this arrangement are $Cd=0.679$ and $Cd=0.665$⁠, respectively.

For medium stratification (Fr = 0.625), the resulting flow for tilted configurations becomes unsteady and shows all the characteristics of the two-dimensional vortex shedding regime identified for flows past a single sphere at higher stratification levels; evidenced by Q-plots in Fig. 20 where the key features from Fig. 15 can be found. The wake behind the spheres is primarily defined by the interaction of fluid parcels which circumnavigate both spheres. The far-wake field shows a merged wake, characterized by vortices shed by the two spheres. The near-wake field in the two horizontal planes, each coinciding with the center of a sphere, has two non-symmetric recirculations behind the first sphere and one in the lee of the second (not shown). For $θ′=30°$⁠, the fluid parcels below the lower dividing streamline of the first sphere can pass underneath it; however, the parcels above the upper dividing streamline of the first sphere, after passing over its top, are forced to flow around the second sphere, due to insufficient kinetic energy. The parcels above the upper dividing line of the second sphere can rise over its top, while the majority of remaining parcels have lower kinetic energy and travel around the sphere; although a small number of parcels pass through the gap and under the second sphere, their motion is inhibited by parcels with opposing velocities coming from recirculations formed behind both spheres. The overturning of first lee-wave crest is still visible and the lee-wave pattern is shown in Fig. 21, where the y-component of vorticity is displayed. Downstream of both spheres, amplitudes of lee waves from the spheres predominantly sum (e.g., interfere constructively) in the half-plane described by $y=0,z>0$⁠, while, where $y=0,z<0$⁠, destructive interferences of lee waves also contribute, forming a pattern similar to that in the lower panel of Fig. 13 which dissipates quickly downstream. Vortices in the near-wake field are distorted owing to wake interactions. These patterns indicate that the downstream flow is reminiscent of a flow past a single bluff body elongated in the direction perpendicular to the free stream flow. This effectively decreases the Froude number for this configuration, accounting for similarities of flow structures with those of strongly stratified flows past an isolated sphere in the two-dimensional vortex shedding regime.

For the tilted configuration displayed and medium stratification, drag coefficients and related amplitudes averaged over the interval $T=300−400$ are $C¯d=1.312$ and $C′d=0.015$ for the first sphere and $C¯d=1.189$ and $C′d=0.036$ for the second one.

Figure 22 left panel displays the top view of the three-dimensional representation for a stratified flow with strong stratification (Fr = 0.25). The near-wake field shows the same main characteristics as those observed in Fig. 15, while the far-wake field clearly demonstrates the two-dimensional periodic vortex shedding of the flow which resembles that of a single sphere, cf. Ref. 3. At this stratification, the dividing streamline height increases toward higher (and lower) z-coordinates facilitating the planar development of the flow and confining the lee-wave spreading, and this confirms the findings from previously studied configurations. Figure 23 displays the y-component of vorticity. Despite differences in the near-wake field and the absence of destructive interferences downstream, the main flow features displayed in Fig. 18 are also present in this simulation, demonstrating the predominance of stratification in defining the flow patterns also for tilted configurations of spheres. As before, this configuration resembles a single vertically elongated bluff body—effectively decreasing the Froude number.

For the tilted configuration with $θ′=30°$ and strong stratification, the drag coefficient and corresponding amplitude averaged over the interval $T=300–400$ are $C¯d=1.601$ and $C′d=0.013$ for the first sphere and $C¯d=1.189$ and $C′d=0.015$ for the second one.

This numerical study investigates stratified flows past two identical spheres at Re = 300 over a selection of representative configurations. For flows past a single sphere, the experimental and numerical characterizations identified three dominant flow regimes named as follows: a non-axisymmetric attached vortex for weak-, lee-wave instability for medium-, and two-dimensional vortex shedding for strong-stratification. For flows past two spheres these three dominant regimes can still be used for the overall characterization of the main tendencies in flow behavior; however, within these regimes the interference between the flows and the spheres results in numerous variations of structures.

For weak stratification, Fr = 1.667, the wakes behind two spheres are suppressed by the buoyancy forcing and the flow reaches steady state, similarly as for a single sphere. A combination of the buoyancy driven vertical deformation and the flow interactions due to the presence of the spheres results in the creation of doughnut-shaped attached vortex structures which also appear in flows behind a single sphere. The gravity waves radiation begins from both spheres in an analogous way to the single sphere case.

As the stratification increases, gravity waves become more pronounced and their amplitude tends to reach the maximum at $Fr≈1$⁠, similarly as for a single sphere. For the cases illustrated here at Fr = 0.625, the structures and mechanisms of the single sphere lee-wave instability flow regime are also present; however, some of the two recirculations that would appear behind a single sphere in the horizontal central plane are suppressed by the interactions resulting from the presence of two spheres. Tiny recirculations due to lee-wave instability can be observed in the vertical central plane (y = 0) behind the two spheres. The amplitudes of the gravity waves are substantially different when two spheres are present since the waves generated by each sphere contribute to the construction or destruction of the gravity waves induced by the other sphere.

For the stronger stratification, the flow pattern demonstrates two-dimensional vortex shedding set in the horizontal plane, matching behavior observed for a single sphere. The characteristic feature of this regime is fluid parcels traveling around the spheres and producing a periodic motion analogous to that observed for flows past cylinders. The gravity wave interference patterns are still present albeit in a reduced form and generally the waves dissipate downstream sooner than for medium stratification.

Simulations performed with the second sphere behind but tilted in the vertical with respect to the first sphere show how specific configurations of spheres can trigger a flow similar to the flow past a single vertically elongated bluff body, and consequently these simulations result in flows which would typically be expected at higher stratification and bare similarities to cases of hypothetical hills with different heights, cf. Ref. 2. Formation of such patterns is highly dependent on the azimuthal-like angle $θ′$ and the distance between spheres.

In studies of flows past a single sphere,⁴ qualitative agreement was found in DNS and Detached-Eddy Simulations (DES) results at the same Froude but different Reynolds numbers. It is likely that the same is true for flows past configurations of two spheres, making the presented study also useful for larger Reynolds numbers. Additionally, we note that the NFT-FV computations, for a neutrally stratified flow past two spheres in the parallel configuration, agree very well with experimental tests reported in Ref. 20 and improve on the numerical results presented there. This validation, together with validations of the NFT-FV model for the stratified flows past a single sphere in Ref. 3, boosts our confidence in the accuracy of the results obtained in this study.

This work was supported in part by the funding received from the EPSRC Studentship Grant No. 1965773 and Horizon 2020 Research and Innovation Programme (ESCAPE2 Grant Agreement No. 800897).

The data that support the findings of this study are available from the corresponding author upon reasonable request.