Computational Fluid Dynamics Approach for Oscillating and Interacting Convective Flows

Abstract

The oscillation and collective behavior of convective flows is studied by a computational fluid dynamics approach. More specifically, the rising dynamics of heated fluid columns is simulated in gravitational field using a simplified 2D geometry. The numerical method uses the FEniCS package for solving the coupled Navier–Stokes and heat-diffusion equations. For the flow of a single heated fluid column, the effect of the inflow yield and the nozzle diameter is studied. In agreement with the experiments, for a constant nozzle diameter the oscillation frequency increases approximately linearly as a function of the the flow rate, while for a constant flow rate the frequency decreases as a power law with the increased nozzle diameter. For the collective behavior of two nearby flows, we find a counter-phase synchronization and a decreasing trend of the common oscillation frequency with the distance between the jets. These results are in agreement with the experiments, and our computational study also suggests that the phenomenon is present on largely different length-scales.

1. Introduction

Our recent experimental results reported that rising gas columns can perform oscillations and their interaction leads to fascinating collective behavior [1]. The oscillations and their related instabilities have been previously known (see for example [2,3,4,5,6]), and this problem is still actively in the focus of the scientific community [7,8]. Besides many refined experiments, theoretical studies based on simple hydrodynamics [1,5], theory of dynamical systems [5], impulse response [2], scaling theory [8], linear stability analyses [2,9], and numerical fluid dynamics [6,7] were considered. Although the emerging oscillations are well-studied, to the best of our knowledge there are no theoretical studies on the interaction and collective behavior of nearby jets.

The collective behavior of convective flows can be discussed in analogy with the very similar phenomena known for diffusive flames [10,11,12,13,14,15,16,17,18]. For interactive jets, the toy-model presented in Ref. [1] is inadequate to explain the fine details of the observed phenomena, therefore a more sophisticated theoretical approach is needed. On the other hand, we also believe that this intriguing phenomena is present on larger length-scales as well, being relevant to industrial processes also. The present study contributes in this sense, by considering a numerical hydrodynamics approach to this puzzling phenomenon.

For experimental results, we consider as a reference our previous study realized with a controlled flow of Helium into air [1]. In these experiments, the Schlieren technique [19,20] was used to visualize the flow, also allowing a digital processing of the oscillations. From the images processed by the Otsu method [21,22], the characteristic frequency and the relevant synchronization order parameter was derived. For a better understanding of the phenomena, some sample movies with original recordings and the ones processed with the Otsu method are provided on our YouTube channel [23] and are uploaded also as Supplementary Materials for this article. For a single flow column, the experiments investigated the effect of the nozzle diameter and flow rate on the observed oscillation frequency. For the collective behavior of two nearby flows, our experiments investigated as a function of the separation distance between the flows (i) the phase difference between the oscillations, (ii) their common oscillation frequency, and (iii) a proper synchronization order parameter.

At constant Helium flow, the oscillation frequency of the rising gas column decreases in form of a power law as a function of the nozzle diameter. This finding is similar with the observed oscillation frequency of the flames of candle bundles as a function of the number of candles in the bundle [18]. For a constant nozzle diameter it was found that the oscillation frequency of the flow increases linearly with the flow yield. For the collective behavior of two nearby and clearly separated flow columns with similar flow parameters, only counter-phase synchronization was observed. This is somehow different from the phenomena observed for candle bundle flames, where both in-phase and counter phase synchronization is present depending on the distance between the flames [18]. The experiments concluded that for short distances, the oscillation frequency of the flow column becomes significantly higher than the frequency observed for non-interacting Helium columns with the same parameters (flow rate and nozzle diameter). All the above summarized results should be a test for any model and numerical approach on this intriguing phenomenon.

Due to the complexity of the problems related to flows in different spatial configurations and environments, the computation approaches are often the only theoretical possibilities to realistically model such phenomena (see for example Refs. [24,25]). Even with such a modeling methodology, imposing the right boundary conditions and offering a proper discretization of space and time raises many technical challenges [26]. The incredible revolution we experience nowadays in computational resources and methods have helped us overcome much of these difficulties, and computational fluid dynamics have become the primary tool to investigate theoretical problems related to fluid dynamics. However, even with the presently available computational power, we are often forced to investigate a simpler flow topology and reduce the dimensionality of the problem [27]. This is nowadays a standard procedure for cases where the problem becomes computationally difficult in 3D. A two-dimensional simplification is usually considered when the periodicity and symmetry of the considered flow allows for it. Assuming in the following a cylindrical symmetry for the flow, we consider a two-dimensional numerical fluid dynamics approach for the above mentioned phenomenon. First, we discuss the theoretical background on which our approach is built and the details of the applied numerical method. Using simple and straightforward examples, we thoroughly test the simulation environment to gain confidence in the method. After this methodological part, we approach the proposed problem and compare the results of the simulations with the experimental data from Ref. [1]. Finally, the conclusions are drawn and the universality features of this intriguing phenomenon are discussed.

Before presenting our simulation methodology, we have to mention that we use equations and system parameters in a dimensional form, rather than following the accepted methodology with dimensionless variables. The reason for this is that in our approach we need to take into account the spatial and temporal variation of the density that is connected with the temperature field. In such cases we cannot use a constant Reynolds number, and the numerical advantage of the dimensionless formalism is not obvious. On the other hand, by using dimensional variables and parameters, the connections with the experimental conditions, time, and length-scales are more straightforward.

2. The Numerical Approach

We present here a 2D numerical approach, which is suitable for modeling the oscillations and collective behavior observed in the rising gas columns. In order to further simplify the problem, instead of a Helium column injected from the bottom we consider the flow of the same incompressible fluid as the surrounding, heated in a restricted region at the bottom of the simulated area. In such a manner we get a rising gas column that is also realizable in experiments.

Using the same Schlieren technique as previously, in Figure 1 and in the movies presented in Ref. [28], we show that very similar instabilities and oscillations occur. In these experiments the heating is realized by a simple heating coil in which one controls the dissipated electric power. Unfortunately, in such experiments there is no good control over the flow debit, therefore one cannot conduct such carefully monitored experiments as the ones done for Helium. The numerical results are therefore compared with our earlier experiments [1].

The advantage of the proposed setup is that we do not have to apply the numerical fluid dynamics method for two component gases. We do pay however for this simplification by the non-homogeneous temperature field, therefore extra transports and gradients have to be taken into account.

In our approach, the fluid is considered to be ideal, as described by the Navier–Stokes equation. For an incompressible fluid in a gravitational field, the Navier–Stokes equation is written in the following form:

$$\begin{array}{c}\hfill \rho \frac{\partial \mathbf{u}}{\partial t}+\rho \left(\mathbf{u}·\nabla \right)\mathbf{u}=-\nabla p+\mathbf{g}\rho +\mu \Delta \mathbf{u}\end{array}$$

(1)

$$\begin{array}{c}\hfill \nabla ·\mathbf{u}=0\end{array}$$

(2)

Here $\rho$ denotes the density, p is the pressure, g the gravitational acceleration, u is the velocity of the fluid, and $\mu$ denotes the fluid’s viscosity. The quantities u, $\rho$ and p can be time- and position-dependent in the flow-space.

For the considered problem, the convective flow due to the temperature difference plays a key role, therefore in the Navier–Stokes equation, we will take into account the temperature dependence of the density and also describe the time evolution of the temperature inside the fluid. As previously emphasized, this is the main reason as to why the dimensionless form of the equation does not reduce the numerical complexity of the problem.

The evolution of the temperature and the temperature dependence of the density are approximated by the following equations:

$$\begin{array}{c}\hfill \rho =\frac{{\rho }_0}{1+(T-T_0)·\alpha }\\ \hfill \frac{\partial T}{\partial t}=D·\Delta T-(\mathbf{u}·\nabla )·T\end{array}$$

(3)

In the above equations, ${\rho }_0$ is the density at $T_0$, T is the temperature of the fluid at a given spatial position and in a given time-moment, D is the diffusion constant, and $T_0$ is the ambient temperature. The numerical solution of the coupled systems of partial differential Equations (2) and (3) was done by using the "FEniCS" software package [29]. FEniCS is an open-source platform developed for solving Partial Differential Equation (PDE) systems. We chose this platform because it has high-level programming interfaces (C ++, Python), the shape of the equations in the program code is similar to their symbolic form, and the program is optimized for a wide range of hardware from laptops to high-performance clusters.

2.1. The Simulation Code

FEniCS uses finite element methods to solve PDEs. As an example in Appendix A.1, we illustrate how to solve the simple 2D Poison equation with FEniCS. For our specific problem we first deal with the term describing the evolution of the temperature:

$$\begin{array}{c}\hfill \frac{\partial T}{\partial t}=D·\Delta T-(\mathbf{u}·\nabla )·T\end{array}$$

(4)

This equation contains a time derivative, so in addition to the coordinates we also have to discretize time. This is done by the Euler method, as follows:

$$\begin{array}{c}\hfill \frac{T(t+dt)-T\left(t\right)}{dt}=D·\Delta T\left(t\right)-(\mathbf{u}·\nabla )·T\left(t\right),\end{array}$$

(5)

We then bring each term to the left hand side of the equation, we multiply the equation by a $\tau$ test function, and integrate the equation over the entire simulated domain:

$${\int }_{\Omega }[T(t+dt)-T\left(t\right)-D·\Delta T\left(t\right)·dt+(\mathbf{u}·\nabla )·T\left(t\right)·dt]·\tau d\Omega =0$$

(6)

The equations above contain a second-order derivative for the coordinates, which is eliminated by partial integration:

$${\int }_{\Omega }({\nabla }^2·T\left(t\right))·\tau d\Omega ={\int }_{\partial \Omega }\left(\frac{\partial T\left(t\right)}{\partial \mathbf{n}}\right)·\tau ds-{\int }_{\Omega }\nabla T\left(t\right)·\nabla \tau d\Omega$$

(7)

Here we denoted by n the unit normal vector to the $\partial \Omega$ surface. The derivative with respect to n is defined as:

$$\frac{\partial T}{\partial \mathbf{n}}=(\nabla T)·\mathbf{n}$$

(8)

Rewriting Equation (6) using the above result and the fact that under the Dirichlet and free boundary conditions the surface integral disappears, we obtain the final form:

$${\int }_{\Omega }([T(t+dt)-T\left(t\right)+(\mathbf{u}·\nabla )·T\left(t\right)·dt]·\tau +D·\nabla T\left(t\right)·\nabla \tau ·dt)d\Omega =0$$

(9)

The incompressible Navier–Stokes equation (2) was solved using the IPCS (Incremental Pressure Correction Scheme) scheme [30]. The IPCS method consists of three steps, but before specifying the steps we introduce the following functions and notation:

$$\begin{array}{c}\hfill \left[\epsilon \left(\mathbf{u}\right)\right]=\frac{1}{2}·([\nabla \otimes \mathbf{u}]+{[\nabla \otimes \mathbf{u}]}^T)\\ \hfill \sigma (\mathbf{u},p)]=2·\mu ·\epsilon (\mathbf{u})-p·I\\ \hfill {\langle \mathbf{f},\mathbf{g}\rangle }_{\Omega }={\int }_{\Omega }\mathbf{f}·\mathbf{g}d\Omega \\ \hfill {\langle \left[A\right],\left[B\right]\rangle }_{\Omega }={\int }_{\Omega }\left[A\right]:\left[B\right]d\Omega \end{array}$$

(10)

The ⊗ product defines a matrix with the following elements:

$$\nabla \otimes \mathbf{u}=\left[\frac{\partial u_j}{\partial x_i}\right]$$

(11)

We denoted by $[...]$ a square matrix, by ${[...]}^T$ the transpose of a matrix, and by : the inner product of matrices:

$$\left[A\right]:\left[B\right]\equiv \sum _{i,j}{A}_{ij}{B}_{ij}$$

(12)

Using the $\epsilon$ and $\sigma$ functions and the specified notation, the steps of the method will be described in the following. First we reconsider the Navier–Stokes equation using a set of test functions:

$$\begin{array}{cc}& \rho ·{〈\frac{{\mathbf{u}}^{*}-\mathbf{u}\left(t\right)}{dt},\mathbf{v}〉}_{\Omega }+\rho ·{〈\mathbf{u}\left(t\right)·\nabla \mathbf{u}\left(t\right),\mathbf{v}〉}_{\Omega }+{〈\left[\sigma \left(\frac{\mathbf{u}\left(t\right)+{\mathbf{u}}^{*}}{2},p\left(t\right)\right)\right],\left[\epsilon \left(\mathbf{v}\right)\right]〉}_{\Omega }+\\ \hfill & +{〈p\left(t\right)·\mathbf{n},\mathbf{v}〉}_{\partial \Omega }-{〈\mathbf{n}·{[\nabla \otimes \left(\frac{\mathbf{u}\left(t\right)+{\mathbf{u}}^{*}}{2}\right)]}^T,\mathbf{v}〉}_{\partial \Omega }=\rho ·{〈\mathbf{g},\mathbf{v}〉}_{\Omega }\end{array}$$

(13)

Here v is the test function. For more information on making this choice, one should consult [29] The first step of the method is the calculation of an intermediate velocity ${\mathbf{u}}^{*}$ from which the pressure will be determined. Then, the pressure is determined in the $t+dt$ step in equation:

$$\begin{array}{c}\hfill {\langle \nabla p(t+dt),\nabla q\rangle }_{\Omega }={\langle \nabla p\left(t\right),\nabla q\rangle }_{\Omega }-\frac{{\langle \nabla {\mathbf{u}}^{*},q\rangle }_{\Omega }}{dt}·\rho \end{array}$$

(14)

In the equation above, q is a test function for the pressure. In the last step, the velocity in the $t+dt$ time step is determined based on the pressure and the intermediate velocity:

$$\begin{array}{c}\hfill {\langle \mathbf{u}(t+dt),\mathbf{v}\rangle }_{\Omega }={\langle {\mathbf{u}}^{*},\mathbf{v}\rangle }_{\Omega }-\frac{dt·{\langle \nabla \left(p\right(t+dt)-p(t\left)\right),\mathbf{v}\rangle }_{\Omega }}{\rho }\end{array}$$

(15)

The method described above for solving the incompressible Navier–Stokes equation is implemented in 2D. FEniCS uses a triangular adaptive grid to solve the 2D partial differential equation. We carefully verified the grid independence of the results, aspects which will be discussed in the next section. Here, in Figure 2, we illustrate the topology of the grid we used in the simulation space.

In order to solve the equations numerically, we need boundary conditions in addition to discretization. We used Dirichlet and free boundary conditions. The Dirichlet boundary condition means that the value of the quantity at a given point is fixed. In the case of the free boundary condition, the derivative as a function of the coordinates of the quantity at the given point is 0.

We visually tested our 2D simulation environment on two simple problems. First we intended to reproduce the Karman vortices in the flow of a fluid around an obstacle (Appendix A.2). Second, we simulated the expansion and rising of a heated sphere, verifying the code for non-homogeneous temperature conditions as well (Appendix A.3). The test simulations reproduced the expected realistic behavior for these known problems, giving confidence for the correct implementation of the relevant equations discussed above.

2.2. Simulating the Rising Hot Air Column

In the followings we provide the details for implementing the simulations, aiming to reproduce the characteristic oscillations observed in a rising gas column. The boundary conditions introduced for velocity and temperature will be justified, and we explain how the time series of the characteristic oscillations were obtained and how the oscillation frequency was calculated.

We consider the inflow geometry presented in Figure 3, leading to the flow illustrated in Figure 4a. On the sidewalls, the value of the velocity is fixed to 0, on the lower boundary, the x-direction component of the velocity is considered as 0, and the y-direction component is given by the following parabolic-like kernel (see Figure 3)

$$\begin{array}{cc}& v_y(x,0)=c_{1}f\left(x,\frac{d}{2}\right)\left(\frac{d}{2}-x\right)\left(x+\frac{d}{2}\right)+c_{2}f\left(x,\frac{H}{2}(1-\frac{1}{30})\right),\end{array}$$

(16)

with:

$$\begin{array}{cc}& f(a,b)=\frac{1}{e^{-c_3·(b+a)}+1}-\frac{1}{e^{c_3·(b-a)}+1}.\end{array}$$

(17)

In the equation from above, d denotes the nozzle diameter, H denotes the width of the simulated space, the parameters $c_1,c_2$ determine the incoming flow rate of the fluid, and $c_3$ is a tuning parameter governing the cut in each profile.

At the upper boundary (height L), free boundary conditions are applied for the y component of the velocity, and for the x component the Dirichlet condition is applied, i.e., $v_x(x,L)=0$. For pressure, Dirichlet boundary conditions were used in the upper part of the simulated volume, $p(x,L)=g·{\rho }_0·l$, and for the free boundary condition for the other boundaries. The temperature on the walls is fixed to $T_0$. On the upper boundary we consider $T_0$ if the y-direction component of the velocity is negative, otherwise free boundary conditions are used. The temperature at the lower boundary is determined by using the following equation:

$$\begin{array}{cc}& T(x,0)=T_0+T_{heati\mathbf{n}\mathbf{g}}f\left(x,\frac{d}{2}\right)\end{array}$$

(18)

Here, the $T_{heati\mathbf{n}\mathbf{g}}$ temperature governs the form of the temperature profile at $y=0$ height. For simplicity reasons we have used in all the presented results $T_{heati\mathbf{n}\mathbf{g}}=T_0$. In the first attempts at the upper part of the simulated box, the free boundary condition was considered for the velocity. However in such cases, unexpected instabilities occurred and after a certain time the heated fluid column was pushed to one of the sidewalls. We have carefully examined this phenomenon and concluded that a self-amplifying effect is responsible for its development. Due to the convective flow, the amount of fluid leaving the simulation box is larger than the volume of fluid flowing into the simulation box through the lower boundary. Since the fluid is incompressible, the fluid must flow back into the simulation box through the upper boundary. Since there is always an asymmetry in the profile of the fluid inflow, this will slightly deflect the outflowing column. In the direction of the deflection, the inflow area decreases, so the asymmetry in the fluid inflow increases. An increase in asymmetry over time will result in the fluid flowing along one of the walls. This is the simple explanation of the observed instabilty.

Two methods were used to eliminate these instabilities. The first method is to flow a fluid of ambient temperature $T_0$ at a constant rate on both sides of the heated air column. Since the flow is two-dimensional, the fluid flowing on a given side can only leave on the same side and this will always provide a minimum distance from the wall for the rising jet. The second method is to allow only the y-direction component of the velocity at the upper boundary. Combining these two methods will eliminate the tendency of the jet to approach one of the sidewalls.

For the upper boundary, a proper boundary condition has to be applied for the inflowing fluid temperature as well. At the upper boundary, an inflow is also necessary in order to respect the incompressibility of the fluid. Since the temperature of the outflowing fluid varies over a wide range we cannot apply the Dirichlet boundary condition to the whole upper boundary because this would cause unmanageable gradients. Avoiding large gradients due to large temperature differences was solved by applying the boundary condition only to those points where the y-direction component of the velocity became negative.

Before performing our large-scale computations we have checked that the used space-discretization (grid) and chosen time-step does not influence the observed trends. Grid independence was proved by reducing and increasing space and/or time discretization consecutively, and comparing the trends and values for the relevant numerical parameters. In Figure 5, we illustrate the grid independence by plotting the time series of the observed oscillations for different grid sizes. More precisely we plot the number of pixels with an intensity above a given threshold detected at bottom of the simulated area (up to height: $L/3$). The observed oscillation frequency is practically independent of the grid size in case of the refined grids considered in the simulations.

The time series for the relevant hydrodynamical parameters were generated by following the temperature distribution in the simulated space. The characteristic frequency was determined by a Fourier transform and from the Power Spectral Density (PSD) the characteristic frequency was calculated. The used signal is the average of the pixel intensities in the lower part of the simulated area (height smaller than L/3). A characteristic signal and the corresponding PSD is sketched in Figure 6.

For realistically chosen parameters, it was shown that the model is capable of producing an oscillation similar to the one observed in the case of the Helium column. Interestingly, it was found that such oscillations are possible even on largely different length-scales. The observed oscillation is shown in Figure 4a,b. where we illustrate the temperature space at subsequent time moments. For the simulations presented in Figure 4, we used the parameters specified in the figure caption.

For a quantitative evaluation of the simulated dynamics, the Otsu method was applied for the 2D temperature field. To obtain the time series, in uniform time intervals the Otsu processed pixels were summed up to a certain height, after this the obtained time series was divided by its average value. The oscillation frequency was calculated in a similar manner with the experiments, based on the above generated time series. In the first step, a Fourier transform was applied to the time series and then the value of the frequency belonging to the largest peak was determined as the relevant oscillation frequency.

With the implemented simulation code we examined how the inflow rate (yield) of the heated fluid column and the nozzle diameter affects its oscillation frequency. We also investigated the collective behavior for the oscillation of two columns placed nearby each other.

2.3. Numerical Results for the Oscillation Frequency

The effect of flow yield and nozzle diameter was examined on two different length-scales. To study the flow yield we used the parameter sets (a) and (b) introduced above, and the nozzle diameters were $d=0.08$m and $d=8$m, respectively. For constant $c_3$, $c_2$, and d parameters, the yield (flow debit) of the heated fluid only depends on $c_1$:

$$\begin{array}{cc}& \Phi ={\int }_{-\frac{d}{2}}^{\frac{d}{2}}{v}_y(x,0)dx=\\ & ={\int }_{-\frac{d}{2}}^{\frac{d}{2}}[c_{1}f\left(x,\frac{d}{2}\right)\left(\frac{d}{2}-x\right)\left(x+\frac{d}{2}\right)+c_{2}f\left(x,\frac{H}{2}(1-\frac{1}{30})\right)]dx\end{array}$$

(19)

The computed oscillation frequency of the heated fluid column as a function of the $\Phi$ parameter is plotted in Figure 7a and Figure 8a. One can observe that the oscillation frequency increases as the flow rate $\Phi$ increases, and this increasing trend can be well approximated by a linear fit in good agreement with the experimental results plotted in Figure 9.

The effect of nozzle diameter on the oscillation frequency was investigated at a constant inflow yield. Since the yield $\Phi$ depends on d according to Equation (19), for different nozzle diameters we must rescale the parameters $c_1$ so that the flow rate remains constant. For the smaller length-scale simulations, we used ${\Phi }_1=0.076$m${}^2$/s flow yield and for the larger scale simulations, we used ${\Phi }_2=29$m${}^2$/s flow yield. To keep the flow yield for different nozzle diameters constant, we varied the value of the $c_1$ parameter. The $c_1$ values for the different nozzle diameters are shown in

Table 1. Value of the $c_1$ parameter for different nozzle diameters d, in order to keep the flow rate ${\Phi }_1=0.076$m${}^2/$s.

d [m]	0.04	0.06	0.08	0.1	0.12	0.14
$c_1[$m${}^{-1}·$s${}^{-1}]$	6781	1952	800	398	223	136

and

Table 2. Value of the $c_1$ parameter for different nozzle diameters d, in order to keep the flow rate $\Phi =29$m${}^2/$s.

d [m]	7	8	9	10	11	12
$c_1[$m${}^{-1}·$s${}^{-1}]$	0.45	0.3	0.21	0.153	0.11	0.088

.

For both length scales, a decreasing trend of the oscillation frequency as a function of the nozzle diameter was observed. The results in such sense are plotted in Figure 7b and Figure 8b, the trend is in good agreement with the experimental results shown in Figure 10.

2.4. Numerical Results for the Collective Behavior

We now turn our attention to reproducing the experimentally observed collective behavior in form of anti-phase synchronization.

The dimensions of the simulation boxes used to study the collective behavior are as follows: 46 m wide (H = 46 m) and 30 m high (L = 30 m) for the large length-scale and 0.3 m wide (H = 0.3 m) and 0.15 m high (L = 0.15 m) at the smaller length-scale. At the lower boundary, the x component of the inflow fluid velocity is 0, and the y component is given by the following kernel function:

$$\begin{array}{cc}& v_y(x,0)=c_{1}f\left(x-d_0,\frac{d}{2}\right)\left(d_0-x+\frac{d}{2}\right)\left(x-d_0+\frac{d}{2}\right)+\\ & +c_{1}f\left(x+d_0,\frac{d}{2}\right)\left(\frac{d}{2}-d_0-x\right)\left(x+d_0+\frac{d}{2}\right)+c_2\mathbf{f}\left(x,\frac{H}{2}-c_4\right)\end{array}$$

(20)

Here the value of $c_4$ is $0.5$m for the large scale system and $0.005$m for the small scale system.

This leads to an inflow profile with two peaks, where the centers are separated at a distance of $2d_0$, as illustrated in Figure 11. The temperature profile is adjusted accordingly:

$$\begin{array}{cc}& T(x,0)=T_0+T_{heating}f\left(x-d_0,\frac{d}{2}\right)+T_{heating}f\left(x+d_0,\frac{d}{2}\right)\end{array}$$

(21)

We used the same simulation parameters as before and fixed $d=7$m, $c_1=0.45$m${}^{-1}·$s${}^{-1}$ values for the large length-scale and $d=0.04$m, $c_1=6400$m${}^{-1}·$s${}^{-1}$ values for the small length-scale system. Again, for the presented results we considered $T_{heating}=T_0$. The experimental results from Ref. [1] show that at a small separation distance, collective behavior in form of counter-phase synchronization appears. A snapshot for a simulated stable collective behavior is visible in Figure 12, successfully reproducing this counter phase synchronization on the smaller length-scale. Similar behavior is observable for the larger length-scales as well.

For the pictures processed with the Otsu method, the collective oscillation of nearby heated fluid columns are shown in Figure 13a,b, for the small and large length-scales, respectively.

For the indicated separation distances, an almost perfect counter-phase synchronization develops. For larger separation distances, the phases of the oscillations will begin to shift relative to each other and no clear phase-difference blocking is observable.

For the simulations performed on the smaller length-scale, corresponding to the experimental conditions in Ref. [1], we computed the synchronization order parameter, which is meant to characterize the collective oscillation. We used the same z synchronization order parameter as the one used in Refs. [1,18]. The computationally derived synchronization parameter is plotted in Figure 14a. Its values in the neighborhood of −1 indicates that we have counter-phase synchronization for the studied distances. In Figure 14b we also show the oscillation frequency of the two synchronized heated fluid columns as a function of their separation distance. This frequency decreases as we increase the separation distance between the columns, similarly to what has been reported in our experiments for the Helium columns [1].

Similarly with the case of the oscillations for a single flow, we have tested the grid independence of the results for synchronization. On Figure 15, we illustrate the observed collective behavior for two different grid sizes, using the same values for all the other simulation parameters.

3. Discussion and Conclusions

In our previous study [1], we used both experimental and theoretical approaches to investigate whether the hydrodynamic instabilities that occur in rising gas columns are also responsible for the oscillations observed for the diffusion flames [18]. It was shown that this is indeed the case: Helium columns ascending in air from a circular nozzle produce similar oscillations with the ones observed in diffusion flames. In addition, the similar collective behavior of these oscillations (counter-phase synchronization) for Helium columns and flickering candle flames suggest that the hydrodynamic processes by their own are enough to explain these phenomenon. For modeling the observed oscillations, a simplified but analytically treatable hydrodynamic approach was used. The model predicted the right trends for the oscillation frequencies as a function of the relevant parameters, but was unsuitable to approach the collective behavior.

In this study, we offered improved modeling by considering a 2D numerical hydrodynamics computer simulation where, for computational simplicity, heated fluid columns were considered instead of ascending Helium columns. This approach proved to be successful for reproducing the experimentally observed features. For a constant nozzle diameter, the numerics led to an oscillation frequency that increased roughly linearly with the flow yield, which is in agreement with the experimental results. For constant flow yield, the numerical results suggested a decreasing trend of the oscillation frequency as a function of the nozzle diameter, confirming the experimental results. The exact shape of the simulated trend was however slightly different from the one observed in the experiments. The main reason for this discrepancy is most likely the reduction of the real 3D problem to a 2D topology. Finally, the presented computer simulations were successful also in reproducing the counter-phase synchronization of the two heated fluid columns placed nearby each other. The computed trends for the synchronization order-parameter and the collective frequency were also in agreement with the experimental results obtained for Helium columns rising in air.

From a more general physical point of view it is important to notice once again the generality of the spontaneous synchronization phenomena in interacting oscillatory systems. Similarly with the analogous candle flame synchronization, the investigated fluid-dynamical system offers yet another fascinating example in this sense. The Navier–Stokes equation for incompressible fluids coupled with the classical heat-diffusion equation, and by considering a temperature dependent density in a 2D approach, is seemingly enough for reproducing the experimentally observed trends. The use of a 2D topology is based on the assumption of the jet’s cylindrical symmetry. Future experiments will decide whether this is a reasonable approximation. However, performing a realistic 3D fluid dynamical simulation with the FEniCS method was not viable with our available computational resources.

It worth mentioning here that the computer simulations were performed both in the laboratory and for a much larger length-scale than the experiments. The qualitative agreement between the results (trends and collective behavior) on these different length-scales suggests that the investigated phenomenon is more general than it was thought to be, and might have further, yet unexplored, connections.