This paper focuses on a validation study in conjunction with an in-house Lattice Boltzmann Method (LBM) based code for a multiphase flow application. The LBM technique takes into account one partial differential equation for each fluid component of multiphase flows to describe the evolution of a particle distribution function. In the present work, a numerical study has been carried out to investigate the behaviour of a heavy and elliptic bubble immersed in a lighter fluid which oscillates until the equilibrium state is reached. The benchmark problem of an oscillating bubble has been implemented, validated and investigated in conjunction with the Colour Gradient Model (CGM)

Figueroa-González, A.

Jenkins, Karl W.

Könözsy, László Z.

Moulitsas, Irene

Oliveira, J.

Teschner, Tom-Robin

English

Cranfield CERES

6th European Conference on Computational Mechanics (ECCM 6),7th European Conference on Computational Fluid Dynamics (ECFD 7), 11-15 June 2018, Glasgow, UKValidation of an in-house Lattice Boltzmann solver for a multiphase flow applicationAntonio Figueroa-González, José Oliveira, Tom-Robin Teschner, László Könözsy, Irene Moulit-sas and Karl W. Jenkins6th European Conference on Computational Mechanics (ECCM 6) 7th European Confere ce on Computati nal Fluid Dynamics (ECFD 7) 11 – 15 June 2018, Glasgow, UK    VALIDATION OF AN IN-HOUSE LATTICE BOLTZMANN SOLVER FOR A MULTIPHASE FLOW APPLICATION  Antonio Figueroa-González1, José Oliveira1, Tom-Robin Teschner1, László Könözsy1, Irene Moulitsas1 and Karl W. Jenkins1 1Centre for Computational Engineering Sciences, Cranfield University, College Road, Cranfield, Bedfordshire, MK43 0AL, United Kingdom,  antonio.figo0693@gmail.com, jose.p.s.olive@gmail.com, t.teschner@cranfield.ac.uk, laszlo.konozsy@cranfield.ac.uk, i.moulitsas@cranfield.ac.uk, k.w.jenkins@cranfield.ac.uk   Key Words: Validation and Verification, Multiphase Flows, Lattice Boltzmann Method.  Abstract. This paper focuses on a validation study in conjunction with an in-house Lattice Boltzmann Method (LBM) based code for a multiphase flow application. The LBM technique takes into account one partial differential equation for each fluid component of multiphase flows to describe the evolution of a particle distribution function. In the present work, a numerical study has been carried out to investigate the behaviour of a heavy and elliptic bubble immersed in a lighter fluid which oscillates until the equilibrium state is reached. The benchmark problem of an oscillating bubble has been implemented, validated and investigated in conjunction with the Colour Gradient Model (CGM). 1  INTRODUCTION  Multiphase flows is one of the main research areas of the Lattice Boltzmann Method (LBM). The main reason is the flexibility of the LBM formulation in comparison with Navier-Stokes solvers when dealing with complex physical phenomena such as multiphase flows. The traditional Navier-Stokes equations techniques; Volume of Fluid (VOF), Eulerian-Eulerian and Lagrange-Eulerian approaches involve the solution of a set of non-linear partial differential equations for each phase or component. In contrast, the LBM techniques for multiphase flows deal with only one partial differential equation for each fluid component or phase to describe the evolution of a particle function distribution. As a result, it provides an efficient framework which requires less computational resources to treat complex flows.  Since the emergence of the LBM for multiphase flows in the late 1980s, the validation of different methods was carried out on simple and theoretically well-known solutions such as the Laplace’s law for the calculation of the surface tension [1], the layered two-phase flow in a channel [2], or relying on Jurin’s and Wahsburn’s laws [3]. In this work, the evolution towards equilibrum of a two-dimensional bubble is analysed in the context of investigations performed to validate an in-house parallelised LBM solver for a multiphase flow application.  The numerical experiment consists of a heavy and elliptic bubble immersed in a lighter fluid which oscillates until it reaches an equilibrium state. The initial curved interface induces a damped oscillating movement whose frequency can be calculated analytically for inviscid two-dimensional flows [4, 5]. This theoretical prediction requires to know a priori the surface tension value of the interface which is taken from the numerical method. Hence, this numerical experiment allows to validate the performance of the in-house LBM code to reproduce an unsteady process in which case the surface tension is taken into account.  2422A. Figueroa-González, J. Oliveira, T.-R. Teschner, L. Könözsy, I. Moulitsas and K. W. Jenkins  2 2 COMPUTATIONAL METHODOLOGY  The multiphase Lattice Boltzmann Model (LBM) implemented in this work is the so-called Colour Gradient Model (CGM), firstly proposed by Rothman and Keller [4]. In particular, the version employed is based on the method described by Leclaire et al. [5, 6, 7, 8, 9]. Their version is stable for very high density ratios and allows to control the surface tension and thickness of the interface through model parameters.   The Density Distribution Function (DDF)  tf ki ,x  is associated to each fluid to describe the flow on the lattice, determined by the speed velocity of its links ie , where  1,,1,0  mi  . The lattice D2Q9 is used in this study, hence m = 9, and k indicates the colour of the fluid (r for red and b for blue). Besides, another DDF, the so-called colour-blind function, is defined as birii fff  . The algorithm steps from time t to t+1 are as follows: 1. Single-phase collision:       tftf iii ,,1 xx . 2. Two-phase collision:        tftf iii ,,2 xx  . 3. Re-colouring operation:        tftf iiki ,,3 xx  . 4. Streaming:     tftf kiiki ,1, xex . 5. Apply boundary conditions. 6. Update macroscopic variables: ikik f ,  k k ,   i k ikif eu . The first three steps are further explained hereafter. 2.1 Single-phase collision The single-phase collision is based on the standard BGK operator in LBM for one fluid. The collision frequency eff  determines the relaxation of the DDF to the equilibrium state as           eqiieffiii ffff  1 ,        (1) where eqif  is the equilibrium function, given by           uuueueeeu23293: 2iiiiiiiieqi wGf 	 ,     (2)  where   is the density weighted average kinetic viscosity. Given the red and blue kinetic viscosities r  and b ,   can be calculated as     bbrbrbrr111 ,        (3)  and when   is known, the effective collision frequency is calculated by      5.031eff .         (4) The operators ٔ and ǣ are the tensor product and the tensor contraction respectively. The coefficients ߰௜, ߦ௜, ߶௜, ߮௜ and ݓ௜ are weights which depend on the lattice scheme and are shown in Table 1.  Coefficients ሼ݅ǣ ȁࢋ݅ȁଶ ൌ Ͳሽ ሼ݅ǣ ȁࢋ݅ȁଶ ൌ ͳሽ ሼ݅ǣ ȁࢋ݅ȁଶ ൌ ʹሽ ݓ௜ 4/9 1/9 1/36 ߶௜,  0 1/5 1/20 ߮௜ 1 -1/5 -1/20 2423A. Figueroa-González, J. Oliveira, T.-R. Teschner, L. Könözsy, I. Moulitsas and K. W. Jenkins  3 ߰௜ -8/3 -1/6 1/12 ߦ௜ 0 1/2 1/8 ܤ௜ -4/27 2/27 5/108 ߠ௜ 1/3 1/3 1/12 ࢋ௜ ࢋ଴ ൌ ሺͲǡͲሻ ࢋଵ ൌ ሺͳǡͲሻ ࢋଶ ൌ ሺͲǡͳሻ ࢋଷ ൌ ሺെͳǡͲሻ ࢋସ ൌ ሺͲǡെͳሻ ࢋହ ൌ ሺͳǡͳሻ ࢋ଺ ൌ ሺെͳǡͳሻ ࢋ଻ ൌ ሺെͳǡെͳሻ ࢋ଼ ൌ ሺͳǡͳሻ Table 1: D2Q9 weights and lattice velocities. On the other hand, the tensor G  is         TG   uu .        (5)  The last parameter  required for computing the equilibrium function is given by     bbrbrbrr  ,        (6) where k  relates the pressure kp   with its density as      kkkp 1 .          (7) In the D2Q9 speed model,   is 3/5. To compute the equilibrium function at every node, it is required to know the density and velocity, which are given for multi-component flows by        ikik f ,  k k ,   i k ikif eu .        (8)  In the above equations, only one k  is a free parameter since they are coupled by the density ratio   between fluids as      rbbr1100,          (9) where 0r  and 0b  state for the initial values of the densities. The red fluid is often chosen as the heaviest one, hence 10  rb  . 2.2 Two-phase collision The two-phase collision operation introduces the interaction between fluids and the two-phase collision can be described and modelled by         	 iiiiii BwAff 222FeFF ,     (10) where iB  is shown in Table 1. This multiphase method is named after the colour gradient vector F. It measures the spatial “colour difference” in terms of densities by         brbrF .       (11)  2424A. Figueroa-González, J. Oliveira, T.-R. Teschner, L. Könözsy, I. Moulitsas and K. W. Jenkins  4Numerically, it has been computed through the approximation as            iibiribirii eexexexexF80,      (12) where i  is also shown in Table 1. One of the main advantages of this method in comparison with others multiphase LBM is that it is possible to prescribe the surface tension between fluids  . This is done through the user-defined parameter A as       effA 94 .        (13) 2.3 Re-colouring This step guarantees the immiscibility of the fluids. In this work, the operator of Latva-Kokko [10], proposed by Leclaire et al. [3], is employed.  The re-coloured-red function is            0ueqiibrirrii fff  cos3 ,     (14) while the re-coloured-blue function is obtained through            0ueqiibrirrii fff  cos3 ,     (15) where   0ueqif  is given by Eq. (2) imposing 0u  . On the other hand, i  is the angle between the colour gradient and the link ie . Finally, the parameter   controls the thickness interface and must be taken between 0 and 1. The higher the value, the stronger the effect of the re-colouring and therefore the thinner the interface.  3  NUMERICAL SIMULATIONS OF AN OSCILLATING BUBBLE   The numerical experiment consists of a heavy and elliptic bubble immersed in a lighter fluid which oscillates until it reaches an equilibrium state. It is based in the 2D case reported by Leclaire et al. [8]. Other authors such as Yuan [11] and more recently Saito et al. [12], the latter being a 3D study with a low density ratio, have also employed this numerical experiment as benchmark case. This numerical experiment allows to validate the behaviour of the implemented algorithm in a complex flow; an initially non-equilibrium bubble evolves until becoming circular following an unsteady movement governed by the surface tension. A square of 1282 lattices has been employed. The red fluid which is the heavier one, initially occupies the nodes located at the interior of an ellipse as          1246416642222 yx ,       (16) whereas the blue fluid occupies the remaining nodes. All boundaries are periodic, which implies that the density distribution functions are streamed through the exterior of the boundary by copying them into the unknown directions on the opposite side. The initial curved interface induces a damped oscillating movement whose frequency can be calculated analytically for inviscid two-dimensional flows. According to Toutant [13], and Lemonnier 2425A. Figueroa-González, J. Oliveira, T.-R. Teschner, L. Könözsy, I. Moulitsas and K. W. Jenkins  5 and Jamet [14], for an inviscid and two-dimensional fluid, the nth-mode angular frequency of oscillation is given by         32 11eqbrn Rnnn   ,      (17) where eqR  is the radius of the bubble in equilibrium. This theoretical prediction requires to know a priori the surface tension value of the interface which is taken from the numerical method. Additionally, n = 2 since it corresponds to an ellipsoidal initial shape. Although the fluid is not strictly inviscid in the LBM setting, the viscous effects are small enough to be neglected and demonstrate the agreement of the numerical experiment with the theoretical prediction of the oscillation frequency. The numerical oscillation frequency is extracted from the simulations. This is done by tracking the maximum height of the bubble contour. This contour is the line in which 0 br  . Thus, to find the maximum height of the bubble at each time step, at the end of each iteration, the quantity br    is evaluated at the nodes placed at the vertical symmetry line. A linear interpolation between the nodes in which the quantity br    jumps from positive to negative values is computed to obtain the maximum height. By performing this operation at the end of each iteration, the evolution of the height contour is captured, and the numerical period of the oscillation is predicted by            232 ttTnum  ,       (18)  where  3t  and  2t  are the times corresponding to the third and second extrema of the movement, respectively, recorded after 500 time steps to avoid the initial noise in the data. From numT , it is straightforward to compute the numerical frequency by         numnum T2 .       (19) Finally, the error in percentage is calculated as      thnumthE  ,       (20) where th  is the theoretical frequency given by Eq. (17). The numerical set up employed used in this study are summarised in Table 2.   Ȟ ߩ௥଴ ߙ௕ ߚ ߥ௥ ൌ ߥ௕ ߪ 50 1 4/9 0.99 1/30 ͷ ή ͳͲିଷ,ͳͲିଷ,ͷ ή ͳͲିସ Table 2: Numerical set up parameters for oscillating bubble simulations. Table 3 shows the relative error in percentage between the numerically and the theoretically predicted frequency values. The results show that the in-house code is able to reproduce the oscillation correctly considering the very small viscosity of the heavier and lighter fluids.   ߪ ߱௧௛ ߱௡௨௠ ܴ݈݁ܽݐ݅ݒ݁ܧݎݎ݋ݎ[%] ͷ ή ͳͲିଷ ͳǤͺ͸ ή ͳͲିଷ ͳǤͺʹ ή ͳͲିଷ 2.15 ͳͲିଷ ͺǤͳͷ ή ͳͲିସ ͹Ǥ͸͵ ή ͳͲିସ 6.38 ͷ ή ͳͲିସ ͷǤ͹Ͷ ή ͳͲିସ ͷǤͲͻ ή ͳͲିସ 11.32 Table 3: Oscillation frequency results comparing the theoretical frequency (߱௧௛) with the numerically predicted one (߱௡௨௠) for different surface tension values. 2426A. Figueroa-González, J. Oliveira, T.-R. Teschner, L. Könözsy, I. Moulitsas and K. W. Jenkins  6 The accuracy of the model decreases for the smaller values of the surface tension, as it is expected for multiphase applications of the LBM based Colour-Gradient Method [1,3,6].  Figure 1 shows the evolution of the maximum height of the contour with the number of iterationsfor ߪ ൌ ͷ ή ͳͲെ͵. The numerical frequency of oscillation is extracted from the extremes of this curve.  Figure 1: Maximum height of the contour evolution with the number of time stepsfor ɐ ൌ ͷ ή ͳͲିଷ. The contour is located where the heavy fluid density (red:ɏ୰) is equal to the light fluid density (blue:ɏୠ).  Figure 2 depicts the shape of the bubble at the beginning of the simulations and at the times associated to the extrema of the contour’s maximum height for ߪ ൌ ͷ ή ͳͲെ͵. The top left picture shows the initial bubble shape. Due to its non-equilibrium shape, surface tension induces capillary waves and makes the bubble deform and oscillate about the circular shape, i.e. the equilibrium shape. At the interface regions where the curvature is high, the surface tension pushes the red fluid inwards; causing the low curvature interface to move outwards. Since the viscosity is very low and the surface forces are strong, the induced movement causes the ellipse major axis becomes the minor axis and vice versa, resulting in an oscillating movement which is slowly damped until the bubble adopts a circular shape. It is observed that the interface is very thin due to the high value of ߚ, which is the fundamental parameter in the recolouring step that assures the segregation of the fluids in order to remain immiscible.    ݐ ൌ ͳͲ ݐ ൌ ݐሺଵሻ ൌ ͳͺ͵Ͳ 2427A. Figueroa-González, J. Oliveira, T.-R. Teschner, L. Könözsy, I. Moulitsas and K. W. Jenkins  7   ݐ ൌ ݐሺଶሻ ൌ ͵ͷ͹ͷ ݐ ൌ ݐሺଷሻ ൌ ͷ͵Ͳʹ   ݐ ൌ ݐሺସሻ ൌ ͹Ͳ͹ͳ ݐ ൌ ݐሺଵሻ ൌ ͺͺͳͷ Figure 2: Snapshots of the oscillating bubble at the beginning of the simulations and at the time steps associated to the extrema of the maximum height of the contour for ߪ ൌ ͷ ή ͳͲିଷ. 4  CONCLUSIONS In this study, we have shown the validation of an in-house Lattice Boltzmann Method (LBM) code in conjunction with the Colour Gradient Model (CGM) to solve multiphase flows through an oscillating bubble. Specifically, we employ the modified CGM by Leclaire et al. [5, 6, 7, 8, 9] and test it here for unsteady multiphase flows for different values in surface tension. We found that the unsteady dynamics of this classical multiphase benchmark case is well reproduced for all values of surface tension, where the error shows a tendency to increase with decreasing surface tension. The numerical results of this bubble dynamics problem presented in this paper show that the implemented in-house code is capable of capturing the physically correct behaviour of the surface tension for a moderate density ratio. This work has demonstrated the code capability to describe unsteady multiphase flows, paving the way to tackle more complex flows for industrial and academic research purposes.   2428A. Figueroa-González, J. Oliveira, T.-R. Teschner, L. Könözsy, I. Moulitsas and K. W. Jenkins  8 REFERENCES [1] A. K. Gunstensen, D. H. Rothman, S. Zaleski and G. Zanetti, Lattice Boltzmann model of immiscible fluids. Phys. Rev. A. 43, 4320, 1991. [2] H. Huang, L. Wang and X. Lu, Evaluation of three lattice Boltzmann models for multiphase flows in porous media. Comput. Math. Appl. Vol. 61, pp. 36063617, 2011. [3] S. Leclaire, A. Parmigiani, O. Malaspinas, B. Chopard and J. Latt, Generalized three-dimensional lattice Boltzmann color-gradient method for immiscible two-phase pore-scale imbition and drainage in porous media. Phys. Rev. E. 95, 033306, 2017.  [4] D. H. Rothman and J. M. Keller, Immiscible cellular-automaton fluids. J. Stat. Phys. Vol. 52, pp. 11191127, 1988. [5] S. Leclaire, M. Reggio and J.-Y. Trepainer, Progress and investigation on lattice Boltzmann modeling of multiple immiscible fluids or components with variable density and viscosity ratios. J. Comput. Phys. Vol. 246, pp. 318342, 2013. [6] S. Leclaire, N. Pellerin, M. Reggio and J.-Y. Trepainer, A multiphase lattice Boltzmann method for simulating immiscible liquid-liquid interface dynamics. Appl. Math. Model. Vol. 40, pp. 63766394, 2016. [7] S. Leclaire, M. Reggio and J.-Y. Trepainer, Numerical evaluation of two recoloring operators for an immiscible two-phase flow lattice Boltzmann model. Appl. Math. Model. Vol. 36, pp. 22372252, 2012. [8] S. Leclaire, N. Pellerin, M. Reggio and J.-Y. Trepainer, Enhanced equilibrium distribution functions for simulating immiscible multiphase flows with variable density ratios in a class of lattice Boltzmann models. Int. J. Multiph. Flow. Vol. 57, pp. 159168, 2013. [9] S. Leclaire, M. Reggio and J.-Y. Trepainer, Isotropic color gradient for simulating very high-density ratios with a two-phase flow lattice Boltzmann model. Comput. Fluids. Vol. 48, pp. 98112, 2011. [10] M. Latva-Kokko and D. H. Rothman, Diffusion properties of gradient-based lattice Boltzmann model of immiscible fluids. Phys. Rev. E., Vol. 71, 2005. [11] P. Yuan, Thermal Lattice Boltzmann two-phase flow model for fluid dynamics. Ph.D. thesis, University of Pittsburgh, School of Engineering, 2000. [12] S. Saito, Y. Abe and K. Koyama, Lattice Boltzmann modeling and simulation of liquid jet break up. Phys. Rev. E., Vol. 96, 2017. [13] A. Toutant, Modélisation physique des interactions entre interfaces et turbulence. Ph.D. thesis, Institut National Polytechnique de Toulouse, 2006. [14] H. Lemonnier and D. Jamet, Test-cases for interface tracking methods, Test Case 5: Oscillation of an inclusion immersed in a quiescent fluid (PA), pp. 4147, 2004. 2429

Validation of an in-house lattice Boltzmann solver for a multiphase flow application

http://dspace.lib.cranfield.ac.uk/bitstream/1826/14723/1/Validation_of_an_in-house_lattice_Boltzmann_solver-2018.pdf

Validation of an in-house lattice Boltzmann solver for a multiphase flow application

Abstract

Similar works

Full text

Available Versions

Cranfield CERES