View source for Bubble rising in a spiralling path in a large tank by Lionel Gamet

[[category:laminar]]
[[category:multiphase]]
* '''contributor''': Lionel Gamet
* '''affiliation''': IFP Energies nouvelles, France
* '''contact''': <mail address='lionel.gamet@ifpen.fr' description='author'>click here for email address</mail>
* '''OpenFOAM version''': v2312
* '''Published under''': CC BY-NC-SA license ([https://creativecommons.org/licenses creative commons licenses])

Go back to [https://wiki.openfoam.com/Multiphase_modeling Multiphase modeling].

=[https://github.com/jnmlujnmlu/OpenFOAMTeaching/tree/master/LionelGamet Bubble rising in a spiralling path in a large tank]=

A video of the case can be downloaded [https://github.com/jnmlujnmlu/OpenFOAMTeaching/blob/master/LionelGamet/Bubble26_CanoLozano_video_isoAdvectorPlicRDF5_level4.avi here].

The starting case case can be downloaded [https://github.com/jnmlujnmlu/OpenFOAMTeaching/blob/master/LionelGamet/Bubble26_CanoLozano.tar.gz here].

Reference OpenFOAM results can be found [https://github.com/jnmlujnmlu/OpenFOAMTeaching/blob/master/LionelGamet/Bubble26_CanoLozano_Results.tgz here].

==Introduction==

This case is a reference test case for VoF simulations. It is the case number 26 of the 3D quantitative benchmark configuration by [1]. The case is also available as an example from the Basilisk website [http://basilisk.fr] and is detailed in the article of Cano-Lozano [1]. It consists in a single rising bubble in a large tank. In the physical conditions of this case, the rising bubble undergoes a spiralling path.   

==Setting up the test case in OpenFOAM==

The bubble rises along +z direction and is initialized as a sphere at z_0/D_0 = 3.5, where D_0 is the bubble initial diameter. The fluid domain is of size 32x32x128 D_0.

The density ratio &rho;_1/&rho;_2 between the fluids is 1000 and the dynamic viscosity ratio &mu;_1/&mu;_2 is 100. Index 1 refers to the continuous liquid phase while index 2 refers to the gas phase. The chosen Bond/Eotvos number Bo = &rho;_1 g D_0^2= 10 and Galilei number Ga = (&rho;_1 g^{1/2} D_0^{3/2})/&mu;_1 = 100.25 classify the current bubble in the oscillatory dynamics regime, with dominant inertial forces [2]. In the simulations, the gravity g and first phase density &rho;_1 are taken as unity, which gives a surface tension &sigma; = 0.1Nm^{-1} and a rise velocity of the order of unity.

The computational grid is obtained by local refinements over a uniform background grid of 40x40x160 cells. The background grid defines the refinement level 0, which thus corresponds to 1.25 cells per bubble diameter. A computational grid with refinement level up to 4 in regions where the bubble can be present is then created with the ''snappyHexMesh'' mesh generator. The level 4 corresponds to a division of cells by a factor 2^4, and so to 20 cells per bubble initial diameter in the refined regions. In order to reduce the number of grid cells, the refinement at the maximum level has been limited to regions in the centerline of the fluid domain, along the bubble rising direction. A refinement cylindrical region of diameter 2D_0 is imposed for 2 &#8804; z/D_0 &#8804; 32. Then a cone of diameter varying between 2 and 4D_0 is used above for 32 &#8804; z/D_0 &#8804; 64. The top of the fluid domain is refined within a cylindrical region of diameter 4D_0 for 64 &#8804; z/D_0 &#8804; 126. The transition between levels is done through buffer layers of two cells (parameter ''nCellsBetweenLevels'' equal to 2 in ''snappyHexMeshDict''). This method conducts to an overall grid size of 9.6 million cells.

A finer grid up to level 5, with thus one more level of refinement than discussed in the former paragraph, is also proposed. The fine grid has almost 73 million cells.

This test case uses incompressible VoF solvers in OpenFOAM, mainly ''interFoam'' or ''interIsoFoam'' solvers. 

The bubble is initialized as a sphere in 3D using the ''setAlphaField'' utility.

'''NB: We use the invert option in ''setAlphaFieldDict'' to reverse the initial field from a droplet to a bubble.'''

In the ''fvSchemes'' file, a Crank-Nicolson second order time scheme with blending coefficient 0.9 is chosen. ''Gauss limitedLinearV 1'' is used to treat the convective term, and ''Gauss vanLeer'' is used for the &alpha; convective term for MULES simulations. The Gauss linear scheme is used by default for all gradient terms.

In the ''fvSolution'' file, the ''GAMG'' implicit solver is used for pressure terms, while the smooth solver is used for the velocity. The PIMPLE algorithm uses only 1 ''nOuterCorrectors'' and 3 PISO correctors (''nCorrectors''=3). It was found that ''momentumPredictor'' needed to be set to true to get a correct solution in terms of rising velocity and bubble sphericity, in particular with isoAdvector.

Both isoAdvector and MULES numerical parameters are present and appear separately in ''fvSolution'', so that bot ''interFoam'' and ''interIsoFoam'' can be run from the same input files.

A constant CFL of 0.05 is used. Computations are run up to time t = 140 s.

Post-processing quantities of interest are described in details in [1]. These are the bubble trajectory, the bubble rise velocity, the bubble sphericity, bubble volume and area. All these quantities are computed through a coded functionObject inlined in the ''controlDict'' file. The results appear in the log file of the solver. We use a ''sampledIsoSurfaceCell'' object, defined as the isosurface &alpha; = 0:5, to compute the bubble area. We also use volume integrals of the gas fraction overall the computational domain to compute the bubble volume, centroid and velocity. A bubble equivalent diameter named D_A is computed from the bubble volume. A bubble equivalent diameter named D_B is computed from the bubble area. The sphericity is defined in 3D as the ratio D_A^2/D_B^2. This number takes the value 1 at t = 0 as the bubble is initialized as a perfect sphere, and then decreases with time as the bubble rises and deforms.

The bubble shape surface is output at ''writeInterval'' frequency through a surfaces sampling functionObject, based upon an isosurface &alpha; = 0:5. Interpolated and non-regularised iso bubble shapes are output.

Finally, the wake of the bubble can be visualized with the &lambda;2 vortex criterion. runTimePostProcessing visualization is also present in this test case, to generate isosurface images of the &lambda;2 vortex criterion. This is currently commented in the controlDict as it requires to be improved and more validated. This will be the object of future releases of this test case.

==Running the case==

[[File:Gamet1.png|540px|right|Bubble 26 trajectory for isoAdvector iso-Alpha. Courtesy from [4].]]

The Allrun_sensitivity script is at the highest level, above Allrun. It is used to run a grid sensitivity for all solvers, ''interFoam'', ''interIsoFoam'', and also the enhanced version of isoAdvector by H. Scheuer [3], now partially integrated in the v2006 with the ''plicRDF'' reconstruction scheme. Allrun_sensitivity should be run first to create running directories for the different VoF solvers. It will generate 6 running directories (3 solvers x 2 levels of grid refinement, as discussed in the former paragraph). The grids are the same for all solvers. Users might thus want to only generate the grid for one of the solvers (e.g. ''interFoam''), and then create links in the others.

The ''Allrun.pre'' script should be run first to create the mesh. The background grid is first constructed by running ''blockMesh'', then s''nappyHexMesh'' is used to refine near the domain centerline. ''renumberMesh'' is used to reduce the matrix bandwith. It is a good practice to always run ''renumberMesh''. Then the ''setAlphaField'' utility is run to initialize the bubble as a sphere.

The Allrun script should then be run to call the VoF solver. These computations are very long and might benefit from AMR in future releases of the test cases. Results are provided for reference in the Wiki. Post-processing quantities
of interest (bubble volume, area, centroid, velocity and sphericity) are extracted from the solver log file through a grep command.

The bubble trajectory is a 3D spiral, drawing almost a perfect circle in the x-y plane for the fine grid case, as can be seen from figure 1 taken from the article in press by Gamet et al. [4], with the standard isoAdvector solver
''interIsoFoam''. A sample wake visualization, using the &lambda;2 vortex criterion, is shown on figure 2 for the level 4 coarse grid.

==References==
[1] J. Cano-Lozano, C. Martínez-Bazán, J. Magnaudet, and J. Tchoufag, Paths and wakes of deformable nearly spheroidal rising bubbles close to the transition to path instability," Physical Review Fluids, vol. 1, no. 5, 2016.

[2] M. K. Tripathi, K. C. Sahu, and R. Govindarajan, Dynamics of an initially spherical bubble rising in quiescent liquid," Nature Communications, vol. 6, 2015.

[3] H. Scheuer and J. Roenby, Accurate and effcient surface reconstruction from volume fraction data on general meshes," J. Comp. Phys., vol. 383, pp.1-23, 2019.

[4] L. Gamet, M. Scala, J. Roenby, H. Scheuer, and J.-L. Pierson, Validation of volume-of-fluid OpenFOAM isoAdvector solvers using single bubble benchmarks, Submitted to Computers and Fluids, 2020.

[[File:gamet2.png|540px|center|Bubble 26 wake visualization. isoAdvector plicRDF solver on coarse grid. Courtesy from [4]. The bubble rises from left to right.]]