{'fraction_non_alphanumeric': 0.021349797341879535, 'fraction_numerical': 0.00909605052314073, 'mean_word_length': 4.747291440953413, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Computational fluid dynamics approach for understanding oscillating and interacting convective flows

A 2D numerical hydrodynamics approach is considered for modelling recent experimental results on the oscillation and collective behavior of convective flows. Our simulations consider the rising dynamics of heated fluid columns in a gravitational field. Simulations are done on two entirely different length-scales, showing also the generality of the investigated phenomena. 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 yield and for a constant flow yield the frequency decreases as a power law with the increasing nozzle diameter. Concerning the collective behavior of two nearby flow we find a counter-phase synchronization of the oscillations and an increasing trend in the common oscillation frequency when the distance between the flows is decreased. These results are again in agreement with the experimental findings.


I. INTRODUCTION
Recent experimental results have proven that rising gas columns can perform oscillations and their interaction leads to fascinating collective behavior 1 . These oscillations have been reported much earlier by Yuan et. al 2 and can be discussed in analogy with the very similar phenomena known for diffusive flames [3][4][5][6][7][8][9][10][11] . The toy-model presented in 1 is clearly inadequate to explain the fine details observed in our recent experiments, therefore a more sophisticated theoretical approach is needed. The present study contributes in this sense, considering a 2D numerical hydrodynamics approach to this phenomenon.
For validating the numerical results, the present work considers as reference the experiments realized using a controlled flow of Helium into air 1 . In these experiments the Schlieren technique 12,13 was used to visualize the flow, allowing also a digital processing of the oscillations.  11 . 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 11 . 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 numerical simulation approach on this interesting phenomenon.
Due to the complexity of the problems related to flows in different spatial configurations, computation approaches are many times the only theoretical possibilities to realistically model such phenomena (see for example 17,18 ). Even with such a modelling methodology, imposing the right boundary conditions and offering a proper discretization of space and time raises many technical challenges 19 . The incredible revolution we experience nowadays in computational resources and methods, helped us to overcome much of these difficulties and computational fluid dynamics became the primary tool to investigate theoretically problems related to fluid dynamics. However, even with the presently available computational power, many times we are forced to investigate a simpler flow topology and reduce the dimensionality of the problem 20 . This is nowadays a standard procedure for those cases where the problem becomes computationally difficult in 3D.
Usually a two-dimensional simplification is considered when the periodicity and symmetry of the considered flow allows this. Assuming in the followings a cylindrical symmetry for the flow, we consider a two-dimensional numerical fluid dynamics approach to the above described phenomenon. First, we discuss the theoretical background on which our approach is build and the details of the applied numerical method. Using simple and straightforward examples we test the simulation environment and gain confidence in the method. After this methodological part, we approach the proposed problem and compare critically the results of the simulations with the experimental data from 1 . Finally, conclusions are drown and universality features of this intriguing phenomenon is discussed.

II. THE NUMERICAL APPROACH
We present here a 2D numerical approach, that is suitable for modelling the oscillations and collective behavior observed in 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 manner we get a rising gas column which is also realizable in experiments.
Using the same Schlieren technique as before, on Figure 1 and on the movies presented in 21 we show that very similar instabilities and oscillations occur. In these kitchen-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 were done for Helium. This is the main The advantage of the proposed setup is that we do not have to apply the numerical fluid dynamics method for two component gases. We pay however for this simplification by the nonhomogeneous temperature field, therefore extra transports and gradients has to be taken into account.
In our approach the fluid is considered to be an ideal one described by the Navier-Stokes equation. For an incompressible fluid in a gravitational field, the Navier-Stokes equation writes in the following form: Here ρ denotes the density, p is the pressure, g the gravitational acceleration, u is the velocity of the fluid, and µ denotes the fluid's viscosity. The quantities u 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 have to take into account the temperature dependence of the density and describe also the time evolution of the temperature inside the fluid.
The evolution of temperature and the temperature dependence of the density, are considered by the following equations: In the above equations, ρ 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 (1)-(2) was done by using the "FEniCS" software package 22 .
FEniCS is an open-source platform developed for solving Partial Differential Equation (PDE) systems. We chose this platform because it has highlevel 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.

A. The simulation code
FEniCS uses finite element methods to solve PDEs. As an example in Appendix A 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 temperature: 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: where the indexes for the T temperature label the discretized time-steps. We then bring each term to the left hand side of the equation, we multiply the equation by a τ test function, and integrate the equation over the entire simulated domain: The equations from above contains a second-order derivative for the coordinates, which is eliminated by partial integration: Here we denoted by n the unit normal vector to the ∂ Ω surface. The derivative in respect to n is defined as: Rewriting equation (5) using the above result and the fact that under the Dirichlet and free boundary conditions the surface integral disappears, we obtain the final form: The incompressible Navier-Stokes equation (1) was solved using the IPCS (Incremental Pressure Correction Scheme) scheme 23 . The IPCS method consists of 3 steps but before specifying the steps we introduce the following functions and notation: The ⊗ product is defining a matrix with the following elements: We denoted by [...] a square matrix, by [...] T the transpose of a matrix and by : the inner product of matrices: Using the ε and σ functions and the specified notation, the steps of the method will be described in the followings. First we reconsider the Navier-Stokes equation in a discretised form on the Ω velocity space: Here v is a test function, for more information on choosing this one should consult 22 The first step of the method is the calculation of an intermediate velocity u * from which the pressure will be determined. Then, the pressure is determined in the i-th step in equation: In the last step, the velocity in the t + dt time step is determined based on the pressure and the intermediate velocity: In the equation from above q is a test function for the pressure. The method described above for solving the incompressible Navier-Stokes equation is implemented in 2D. 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 tested visually our 2D simulation environment in two simple problems. First we intended to reproduce the Karman vortices in the flow of a fluid around an obstacle (Appendix B). Second, we simulated the expansion and rising of a heated sphere, verifying the code for non-homogeneous temperature conditions as well (Appendix C). The test simulations reproduced the expected realistic behavior for these known problems, giving confidence for the correct implementation of the relevant equations discussed above.

B. Simulating the rising hot air column
In the followings we give 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 2 leading to the flow illustrated in Figure   3 (a). On the sidewalls, the value of the velocity is fixed to 0, on the lower boundary, the xdirection component of the velocity is considered as 0 and the y-direction component is given by the following parabolic-like kernel (see Figure 2) with: Here, the T heating temperature governs the extra heat amount introduced in the system.  . and b.). The gravity acts in the negative direction of the y-axis and the parameters of the simulation were chosen as follows: 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.
C. Numerical results for the oscillation frequency The computed oscillation frequency of the heated fluid column as a function of the Φ parameter is plotted in Figure 5a and Figure 4a. One will observe that the oscillation frequency increases as  Figure 3, the value of c 1 for the different nozzle diameters are given in Table I.  Figure 3, the value of c 1 for the different nozzle diameters are given in Table II. the flow rate Φ increases and this increasing trend can be well approximated by a linear fit in good agreement with the experimental results plotted in Figure 6.
The effect of nozzle diameter on the oscillation frequency was investigated at a constant inflow  Table I and Table II Here as it is illustrated in Figure 8. The temperature profile is adjusted accordingly: Again, for the presented results we considered T heating = T 0 . Experimental results from 1 shows that at small separation distance collective behavior in form of counter-phase synchronization appears. A snapshot for a simulated stable collective behavior is visible in Figure 9, reproducing successfully 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 Figures 10 a and b, for the small and large length-scales, respectively.
conditions in 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 11 and 1 . The computationally derived synchronization parameter is plotted in Figure 11 a. It's values in the neighborhood of -1 indicates that we have counter-phase synchronization for the studied distances. In Figure 11b 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 with what has been reported in our experiments for the Helium columns 1 .

III. 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 is responsible also for the oscillations observed for the diffusion flames 11 . It was shown that indeed this is the case: Helium columns ascending in air from a circular nozzle produce similar oscillations with the ones observed in diffusion flames. Also, 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.
In our previous work 1 , for modeling the observed oscillations , a simplified but analytically treatable hydrodynamical 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 the present follow-up, we offered an improved modeling by considering a It worth mentioning here that the computer simulations were performed both for 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.

A. Solving the 2D Poison equation in FEniCS
Here we illsutrate through the 2D Poison equation how to solve a PDE in FEniCS. The Poisson equation can be given in the following form: If we have a simple rectangular space then the above equation can be easily given in the finite element form however, with this simple and intuitive approach, we soon run into problems because even for a circle it is impossible to map the boundary with an acceptably small number of squares. The FEniCS program 22 uses a triangular grid instead of a square grid to cover the simulated space, in which case we can always select the grid so that the grid points are on the boundary surfaces.
Discretization alone does not solve the equation, the next question is how to determine the solution at each lattice point. As a first step, we write the ϕ(x 1 , x 2 ) function in the following form: In the above equation, φ i (x 1 , x 2 ) is a given k-th order polynomial, c i are the coefficients that determine ϕ(x 1 , x 2 ) and N+1 is the number of the grid points. The c i coefficients are determined by multiplying the Poison equation by N+1 different v(x 1 , x 2 ) so-called test functions and integrating the product over the whole domain to obtain N + 1 linearly independent equations from which the c i coefficients can be calculated. All this can be formally given in the following form: We denoted by dΩ = dx 1 dx 2 .The above form of PDE is called the weak formulation of the equation and this is what is calculated by the FEniCS program. The second-order coordinate derivative in the above equation means that the polynomials used need to be twice differentiable. Because the use of polynomials with large degrees requires more memory and computation, we always strive to keep the degree of polynomials to a minimum. In the above equation, the reduction of the order of derivatives can be done by Gauss-Green integration as follows: Since we use Dirichlet boundary conditions for the Poisson problem, the value of v at the boundary is 0, so the equation (24) can be written in the following form: We have seen above how to rewrite the 2D Poisson problem in a form that can be solved with the In these simulations, the length of the simulated volume was considered as 2.2m, the width of the simulated volume is 0.41m and in the middle of the simulated coordinate space a circular obstacle with a radius of 0.05m is placed. The coordinates of the centre of this obstacle was taken at (0.2m, 0.2m). The density of the fluid was taken as unity (1kg/m 2 ), the viscosity is 0.001kg/s, and no gravitational field is considered. For input (the left region of the space in Figure 12) we considered that the velocity in the y direction is 0, the velocity in the x direction has a parabolic profile with a maximum value of 1.5m/s. On the horizontal walls and on the boundary of the obstacle we consider no-slip conditions, thus the velocity is fixed to 0. For the output (rightsight region) we have also imposed for the y direction velocity to be zero, and the pressure at the output is fixed also to 0. Elsewhere there are all free boundary conditions for the pressure. The temperature is fixed at T 0 = 300K at input and on the horizontal walls.
The velocity vector spaces obtained from the simulations are shown in Figure 12 for 4 time moments as it is indicated on the left side of the images. One can observe that the simulation reproduces successfully the expected Karman vortices.
In this second test we aimed to implement the density and temperature evolution of a heated gas sphere in gravitational field. It is expected that the shape of the heated gas will follow the known dynamics of a mushroom cloud in a nuclear explosion.
In the performed simulations, the density at T 0 was chosen as unity (1kg/m 2 ), the ambient temperature T 0 is 300K, the α parameter in equation (2) is α = 0.001 K −1 , the thermal diffusion constant D is 0.3m 2 /s, the initially heated sphere temperature is 600K, the fluid viscosity is taken as µ = 0.05kg/s, the gravitational acceleration is g y = −9.81m/s 2 and the size of the simulated volume is 30m 2 in both the x and y directions.
At the bottom and walls of the simulated space, the velocity is fixed to 0, at the upper boundary, we fix the x component of the velocity to v x = 0. For the pressure, the value of g · ρ 0 · l is fixed for the upper boundary, and free boundary conditions are applied to all the other boundaries. For temperature, free boundary condition is applied at the upper boundary if the y direction component of the velocity is positive otherwise the temperature is fixed to T 0 = 300K units on the upper boundary and on the side-walls. The temperature is fixed to a higher value of 450K units at the bottom-wall of the simulated area. This is necessary in order to make the resulting flow visible in the temperature space. Initially, the velocity in the whole simulated volume is 0 and the volume contains a sphere (disk) with a radius of 5m in which the temperature is 600K. The centre of the sphere is at the coordinates (0m, 7m), the temperature around the sphere is fixed to 300K, and the temperature between the centre and the surface of the disk is given by an interpolation with a sigmoid function.
The time-evolution of the temperature map derived from the simulation is shown in Figure   13. The effect of thermal diffusion can be observed in the first two frames. As a result of this diffusion the initially sharp boundary line between the high and low temperature regions gets blurred. Subsequent frames show the displacement due to convective flow and as a result of this the characteristic mushroom cloud shape is formed.