Computational fluid dynamics simulation of buoyant mixing of miscible fluids in a tilted tube

Buoyancy-driven flows and mixing of fluids with different densities occur frequently both in nature and as part of industrial processes within chemical and petroleum engineering. This work investigates the buoyant exchange flow of two miscible fluids in a long tube with closed ends at varying tilt angles using OpenFOAM. The study focuses on the evolution of the concentration field and front velocities of the mixing zone at different inclinations. Numerical results based on a miscible solver agree with previous experiments and direct numerical simulations. Treating the fluids instead as immiscible with no surface tension leads to unrealistically high front velocities at intermediate inclinations.


Introduction
The gravitational mixing of two miscible fluids, where a light fluid and a heavier fluid displace each other, is present in many natural phenomena [1,2] and industrial processes in chemical and petroleum engineering such as cementing operation. In well cementing by reverse circulation method a heavier fluid displaces a light one [3]. In The density difference between the fluids can lead to Rayleigh-Taylor and Kelvin-Helmholtz instabilities [4], and result in buoyancy-induced turbulence. While gravity-driven flows in nature often take place in relatively large geometries, such flows in industrial processes occur mostly in confined geometries such as tubes and narrow channels. Because of that, there are many numerical and experimental studies on unsteady fluid displacement in long tilted pipes [5][6][7][8][9][10][11][12].
In the mainly experimental studies of Séon et al [5][6][7], where buoyancy-driven exchange flows were examined in long, tilted tubes, three different regimes were identified, depending on the inclination of the tube axis from the vertical, defined by the angle θ. It was found that for inclinations between vertical and 60°, Kelvin-Helmholtz instabilities grow along the interface between the two fluids and the front velocity increases with increasing inclination. For inclinations between 60° and 80°, the more significant transverse component of buoyancy leads to increased segregation and stratification of the fluids. At these inclinations, the front velocity plateaus at a maximum value. In the third regime, at nearly horizontal inclinations, no significant mixing was reported.
There are some computational works on gravity currents as well, which are highlighting differences between 2D and 3D simulations [13,14]. It is suggested that 3D effects cannot be ignored in the long/time dynamics of gravity currents due to the structural differences between 2D and 3D vorticity dynamics. In the study done by Hallez and Magnaudet [10], it is illustrated that the three regimes found in 3D flows at small, intermediate, and large inclinations have no counterpart in two dimensions. It is also shown that direct numerical simulation is a reliable tool for investigating the complex dynamics of IOP Publishing doi:10.1088/1757-899X/1201/1/012021 2 buoyancy-driven mixing of liquids [10]. While previous work has shown that computationally intensive DNS can resolve these challenging exchange flows, the capability of OpenFOAM to do the same has been studied in this paper. Three-dimensional numerical simulations of buoyancy-driven flows of a heavy fluid over a lighter one in a long-inclined tube have been performed using OpenFOAM. The impact of modelling choices on the results has been investigated, using the experimental results of Séon et al [5][6][7] and numerical investigation of Hallez and Magnaudet [10] to validate the computed results.

Numerical method
OpenFOAM was used to simulate the flow using twoLiquidMixingFoam and interFoam solvers. In the "Volume of Fluid" (VOF) method, which these solvers utilize, the volume of a fluid in a cell is computed as , with being the phase fraction in the cell and being the volume of the cell. The governing equations are the Navier-Stokes equations and the momentum and continuity equations are: where S = [(∇ ) + (∇ ) ]/2 is the rate of strain tensor, is the gravity, is the fluid viscosity, is the surface tension force and is the mixture density found from: and the subscripts 1 and 2, signify the light and heavy fluids respectively. Note that twoLiquidMixingFoam is an incompressible multiphase solver for miscible fluids while interFoam solver is for incompressible multiphase immiscible fluids. For immiscible fluids there is also an interfacial tension term in (1), which is set to zero in this study in order to simulate equivalent miscible fluids. Additionally, the molecular diffusivity in twoLiquidMixingFoam is assumed to be negligibly small [14], so diffusion present in the solutions is due to numerical diffusion. Hence, the transport equation for the phase fraction becomes: It should be noticed that the main difference between interFoam and twoLiquidMixingFoam is that interFoam solver includes an interface compression scheme to attempt to maintain a sharp interface between the fluids while in twoLiquidMixingFoam they can mix freely.

Turbulence
For most of the simulations, no explicit turbulence model is applied. However, the computational mesh is fine enough to resolve some of the fluctuations caused by the buoyant mixing. The use of secondorder discretization schemes, as used in this work, has been shown to have an effect similar to a subgrid turbulence model, so-called implicit LES [15]. Hence, the simulations can be considered to have an implicit turbulence model applied. The influence of applying an explicit turbulence model to the simulation is also studied by directly using a large eddy simulation (LES) model. A LES model applies a spatial filter to give the filtered Navier-Stokes equations. The filter width, Δ, is typically taken as the grid size, so there is no filter being directly applied in the numerical code, but the filtering process gives rise to new terms in the momentum equations representing the unresolved part of the turbulence. This is denoted the subgrid scale (SGS) stress tensor, where the bar shows the filtered velocity. The simplest and most used model is the Smagorinsky model [16], in which is calculated from: Here, ���� is the strain rate sensor for the filtered velocity, and is the subgrid-scale eddy viscosity. This eddy viscosity is added to the kinematic fluid viscosity to model the mixing effect of the turbulence. To model the eddy viscosity a one-equation model is applied in this work. This model solves one additional equation for the subgrid kinetic energy, k [17], where is the model constant and the eddy viscosity is modelled as: is also a model constant.
As it is shown in equation 9, the turbulence modelling gives rise to an additional term in the phase fraction equation due to the turbulent mixing, where is the turbulence Schmidt number and it is set to one.

Geometry, mesh and boundary condition
The computational domain is a cylindrical tube of diameter equal to 0.02 m and length 96 . It is shown in Figure 1, along with the computational grid. The grid is based on a similar grid structure as used by Etrati and Frigaard [18]. It comprises of 1242 cells along the tube axis, 24 cells over the tube radius with stretching towards the wall, and 64 cells along the azimuthal direction. Some finer grids were also studied but they led to no significant difference in front velocity. Gravitational acceleration was applied in the vertical (Z) direction. Different inclinations were considered by altering the appropriate vector components of the gravitational acceleration. To ensure a stable solution an adjustable time step was used with a maximum Courant number of 0.5. Second-order discretization schemes were used for the spatial terms and a first-order implicit scheme for the time discretization. In all simulations a no-slip condition is applied to all solid boundaries. The pressure boundary condition was specified as zero gradient. We adjust the flow parameters to the experimental conditions done by Séon et al [6][7][8]

Results
In this study to have a quantitative view of fluid mixing, more focus has been on the front velocity, . To calculate this parameter the spatiotemporal diagrams of the front positions shown in Figure 2 are used. These diagrams show the mean concentration level along the tube length over the time. The mean concentration in a section ( = .) is the average value of concentration over the cross section of the tube.
is computed from the local slope of the boundary of the mixing zone, i.e., from the dashed line separating pure fluid region from the mixing zone in Figure 2.C. The front position is defined as the flow region where < 0.99. As it is clear from the diagrams shown in Figure 2, the front velocity specifically for small inclinations (the slope of the boundary of the mixing zone) decreases after a short period of time. In these cases, the front velocity is estimated based on the late time growth of the mixing zone after the slope has stabilized. On the other hand, the spatiotemporal diagrams of experimental studies done by Séon et al [5] show that the front velocity slowly decreases in time during a fairly long, transitional stage to reach a constant value. Owing to the computational cost, it is too difficult to consider the same length as in the experiment which is 4 m equal to 200 . Hallez and Magnaudet [10] reported that the long-time front velocities obtained in a computational domain shorter than the actual tube in the experiments were slightly larger IOP Publishing doi:10.1088/1757-899X/1201/1/012021 5 than their experimental counterparts, which were obtained after a much longer observation time. In fact, if the front of the current is far enough from an end wall the tube length has no effect on the instantaneous dynamics of the flow. This length directly determines the total time over which the flow evolution can be observed before it is influenced by the end walls so that by a longer domain or longer observation time, better agreement between front velocities of simulations and experiments could be obtained [10]. Here, a long enough tube has been chosen for the simulations where the front velocities obtained from the simulations are in a good agreement with the mentioned references. The computed front velocities normalized by the characteristic velocity, , are plotted against the tilt angle in Figure 3.  There are no significant differences between results of two solvers at the tilt angles between 60° and 80°, where the segregation effect is dominant, and the front velocity reaches a plateau close to 0.7 . ~� is the characteristic velocity of the front resulting from a balance between buoyancy and inertia forces. For the miscible solver the maximum front velocity, , is obtained at the tilt angle of 70°, the same as the experimental results and the value reported by Hallez and Magnaudet. The values of for the miscible solver is equal to 0.696 , where the maximum velocity of front obtained from the experiments was 0.703 . The difference is about 1%, which is in the order of magnitude of the experimental accuracy. This difference in Hallez's simulation is 2%. Figure 4 also compares the mixed fluids simulated by both solvers at the tilt angle of 30°. As it is shown there is a clearer border between two fluids in the results obtained by the immiscible solver, while two fluids have smeared into each other when the miscible one is applied. As a result, there is a higher local density difference when the immiscible solver is used leading to the higher front velocity. It can be seen by comparing the thinner plots in Figure 4, that after some time the front position simulated by the immiscible solver has moved faster.
It is also illustrated in Figure 3 for small tilt angles, increases strongly with . It is due to a symmetry breaking effect when the tube is inclined. When the tube is perfectly vertical, there is not "high side" and no "low side", that is preferred. As soon as there is an inclination, one fluid will tend to move "toward the left" and the other "toward the right" side of the pipe, leading to a more effective exchange of the two fluid. This effect suggests analogies with the Boycott effect, namely, with the enhancement of the sedimentation velocity of particles in a tilted tube with respect to a vertical one due IOP Publishing doi:10.1088/1757-899X/1201/1/012021 6 to a similar transverse segregation [7,19]. In Figure 5   However, Figure 6 shows contour fields for turbulence properties and the resulting concentration field for simulation with LES at tilt angle of 30°. It is clear that turbulent energy is produced in the shear layers of the flow, typically at the interface and near walls. This results in higher subgrid eddy viscosity in these regions, and thereby higher mixing. This can be seen from the more smeared out concentration field compared to the results without a turbulence model. The resulting front velocity becomes higher giving a front velocity 0.458 which is 10.6% higher than the result without turbulence model. This is in less agreement with the experimental result and suggests that the current grid is too coarse for a proper LES, since additional viscosity is added implicitly from the numerical discretization.

Summary
Some three-dimensional CFD simulations for buoyant mixing of two miscible fluids in a tilted tube have been carried out using OpenFOAM. Two different solvers, twoLiquidMixingFoam and interFoam, have It is shown that the front velocity increases with tilt angle up to 60°. Afterward there is a little change for the front velocity. Our studies demonstrate that there is a better agreement between the results of miscible solvers and the references values for the smaller tilt angles, and for higher tilt angles the results of both solvers are closer to each other. By applying LES the influence of a turbulence model is investigated and it is shown that there is a less agreement between the obtained results and the experimental ones, however further work on turbulence modelling should investigate the use of finer grids and a more detailed investigation of resolved versus filtered turbulence scales. sustainable use