Comparison between the Lagrangian and Eulerian Approach in Simulation of Free Surface Air-Core Vortices

: The problematic consequences regarding formation of air-core vortices at the intakes and the drastic necessity of a thorough investigation into the phenomenon has resulted in particular attention being placed on Computational Fluid Dynamics (CFD) as an economically viable method. Two main approaches could be taken using CFD, namely the Eulerian and Lagrangian methods each of which is characterized by speciﬁc advantages and disadvantages. Whereas many researchers have used the Eulerian approach for vortex simulation, the Lagrangian approach has not been found in the literature. The present study dealt with the comparison of the Lagrangian and Eulerian approaches in the simulation of vortex ﬂow. Simulations based on both approaches were carried out by solving the Navier–Stokes equations accompanied by the LES turbulence model. The results of the numerical model were evaluated in accordance with a physical model for steady vortex ﬂow using particle image velocimetry (PIV), revealing that both approaches are sufﬁciently capable of simulating the vortex ﬂow but with the difference that the Lagrangian method has greater computational cost with less accuracy.


Introduction
Formation of air-core vortices in dam reservoirs increases the chance of air entrance into the power plant structures resulting in irreversible damage [1][2][3][4][5]. Therefore, provision of a deep insight into the vortex flow field by means of comprehensive investigations facilitates effective prevention and control of the phenomenon.
Various experimental researches aiming to provide knowledge for profound understanding of the flow field employing different types of flow measurement instruments have been conducted so far [6][7][8][9][10][11][12]. Owing to the recent technological advances and development of laboratory set-ups namely the PIV, numerous laboratory experiments have been performed for precise measurement of velocity field and the relevant characteristics of free surface vortices [5,[13][14][15][16].
Alongside the laboratory research, numerical investigations have also been of interest due to the advantages regarding the low economical and computational costs. Numerical simulations are basically conducted using the Eulerian and Lagrangian methods. In the Eulerian method, the flow is treated as a continuous phase mainly based upon the finite volume scheme while in the Lagrangian method the flow is considered as a discrete phase of particles for which the equations of flow are to be solved. Each of the methods has been identified by certain pros and cons, adoption of which is mainly dependent on the purpose of the simulation and the influential variants affecting the phenomenon [17].
The Eulerian-based methods are widely known and desired for simulation of different kinds of flows. Several researchers have managed to simulate the vortex flow field using the Eulerian approach [18][19][20][21][22]. The acceptable performance of the Eulerian approach in simulation of the vortex flow field using the VOF free surface model together with the LES turbulence model has been reported by several researchers [23][24][25].
Nakayama et al. [26] compared the obtained results from simulation of the vortex in a vertical intake using the LES method with some experiments conducted using the PTV method. Moreover, Sungur [27] simulated the vortex flow in asymmetric horizontal intakes using the Flow-3D where the effect of geometric and hydraulic parameters namely the distance of the side wall from the horizontal intake, the diameter of the intake, and the flow rate were investigated; the critical submerged depth for each test was also compared with laboratory data. The obtained results suggested that this model could be regarded as a reliable numerical model for observation of vortex formation in horizontal intakes. In the recent Zi et al. [28] research on numerical investigation of air-core vortices in horizontal intakes as two-phase flow using the LES turbulence method, the flow patterns and the dynamic of the vortices at different stages of formations were scrutinized.
The Lagrangian methods are advantageous in locations where the moving boundary (free surface flow) or complex geometries exist since no grid generation is required [29]. Precise computation of the water free surface could be regarded as the main advantage of the Lagrangian methods for intricate flows namely breaking waves [30], dam break [31,32], dam overtopping [33], and interactions between waves and coastal structures [34,35]. The smoothed particle hydrodynamic (SPH) method is one the major and widely-used Lagrangian methods including micro to macro scale cases or even astronomical issues [36]. Development of parallel processing techniques on graphic cards (GPU) facilitates application of the time-consuming Lagrangian methods for three dimensional simulation of complex flows in spite of the relevant high computational costs [37][38][39][40].
Review of the numerical investigations demonstrates that application of the SPH method for investigation of the vortex flow field has not received adequate attention so far. To this end, the SPH method was incorporated for investigation of the vortex flow using the Lagrangian approach in the present study. Furthermore, the capability of the open source DualSPHysics code in simulation of complex flows using the SPH method along with parallel processing on a powerful graphic card was exploited herein [41][42][43].
According to that mentioned, the main objective of the present study concerned simulation of steady state vortex flow by incorporation of the Lagrangian model (SPH method) which had not yet been considered. The secondary objective was to reach a comparative conclusion on performance of the Eulerian (Finite Volume) and Lagrangian approaches in the simulation of a complex flow such as vortex flow. Experimental study was also conducted on a steady vortex flow over a vertical intake where velocity field was measured using the PIV. The numerical results were then evaluated with the experimental data. The results obtained from the present study could shed light on the performance of the present Lagrangian and Eulerian methods in the simulation of a complex flow field such as vortex flow.

The Lagrangian Approach
In the Lagrangian approach, the flow is considered as a discrete phase of particles moving in space and carrying specific computational information. The Lagrangian method is based on discretization of the integral approximation on the particles. For example, the value of the "f " function in particle "i" can be calculated using Equation (1) in which m j and ρ j are representatives of the mass and density of particle "j" while N corresponds to the total number of the particles in the field of influence of particle "i".
Moreover, W corresponds to the Kernel function for determination of the weight of the particles while h represents the smoothing length specifying the efficacy range of the Kernel function [44]. As a matter of fact, Equation (1) implies that the value of the "f " function for each particle is a weighted averaging of "f " function values for the existing particles under the influence of that specific particle. The Kernel function needs to be even integration which on the relevant domain equals unity. The Quintic function is one subset of the Kernel function as defined in Equation (2) [45] leading to the best balance between stability of solution and the corresponding computational cost [46].
With regard to Equation (2) q = r/h in which "r" represents the distance of the two particles and "h" corresponds to the smoothing length. In addition, the α d coefficient for two and three dimensional corresponds to 7/4πh 2 and 21/16πh 3 respectively.
The smoothed particle hydrodynamics (SPH) could be classified as one of the principal methods with regard to the Lagrangian approach which is currently becoming more and more practical in science and engineering. With respect to the SPH method and assuming negligible compressibility of the fluid (WCSPH), the system of the principal governing equations of the flow (conservation of mass, momentum, and equation of motion) could be written as presented in Equation (3): In which P, ρ, V, r, m, t, g correspond to pressure, velocity, location, mass, time, gravitational acceleration, and ε is a coefficient valued between zero and one. Moreover, the kinematic viscosity of the fluid is assumed to be 10 −6 (m 2 s) and the stress tensor → τ ij in the SPS method could be defined in accordance with Equation (4). The SPS method is applied aiming to exert the LES turbulence model [44].
With regard to Equation (4), ν t = [(C S ∆l)] 2 |S| represents the eddy viscosity, k corresponds to the turbulent kinetic energy, C S is the Smagorinsky constant (C S = 0.0066), ∆l corresponds to particle spacing and |S| = 2S ij S ij 0.5 in which S ij represents the strain tensor in SPS [34]. Furthermore, the pressure is calculated according to Equation (5) by modification of the state equation using artificial reduction of the sound speed where the specific weight and the reference density are calculated as γ = 7, ρ 0 = 1000 kg m −3 respectively, B 0 = ρ 0 c 2 0 /γ and the speed of sound in the reference density can be calculated as: c 0 = c(ρ 0 ) = (∂P/∂ρ) ρ 0 [47].
The scalar field of the density populated by the stiff density field is described by the state equation, the low amplitude oscillations characterized by high frequencies together with the innate irregularities highly increase the chances of pressure fluctuations. In attempts to mitigate the resultant pressure fluctuations, Marrone et al. [48] proposed a formulation for the delta SPH correction in which a corrective term was added to the continuity equation as follows: where the coefficient δ used to tune, modulates the intensity of the diffusion (inactive when δ = 0). A value of δ = 0.1 was recommended for most applications and also utilized in the present study [44]. The solid boundaries in the SPH method are simulated by employing dynamic particles where the continuity and the state equations, excluding the equations concerning the particle motion and the momentum equation, are solved for particles constituting the solid boundary as well as the fluid particles, such that the particles of the entire solid boundary remain still [49]. It is noteworthy that, the density of the boundary particles is increased when the fluid particles are located at a distance less than twice the smoothing length, resulting in higher pressures. Under such circumstances, some kind of repulsive force is exerted on the fluid particles affected by the pressure term in the momentum equation [44,47].
The present study dealt with Lagrangian simulation of the flow by taking advantage of the open source DualSPHysics codes using the CUDA programming language, allowing execution of the main computations on GPU, capable of providing the possibility to perform intricate simulations with a large number of particles rather faster than CPU. The remarkable capability of the DualSPHysics in simulation of complex geometries is considered as one of the practical advantages of the code [44].
Herein, The DualSPHysics v4.4 promoted by implementation of open boundary conditions for inlet/outlet simulation based on definition of a buffer layer in the inlet/outlet boundary was incorporated. The cardinal task of the particles in the buffer layer created/omitted with respect to the boundary type is the reduction of the Kernel function discretization error in the vicinity of the boundaries [50]. The physical properties of the particles in the buffer layer (density and velocity) could be either incorporated as a constant value or calculated based on interpolation from the environment. The open boundary conditions for simulation of the steady state vortex flow were utilized in the present study.

The Eulerian Approach
The flow field in the Eulerian approach is mainly considered as a continuous phase for which the Navier-Stokes equations for the flow field are solved based on methods namely the finite volume scheme (FVM). The three-dimensional Navier-Stokes equations in cartesian coordinates can be written as: where u i denotes the component of the velocity in the direction, x i corresponds to the coordinate of "i" in the Cartesian system, ρ, P and g represent the density, pressure, and gravitational acceleration respectively. Moreover, τ ij is the descriptor of the Reynolds stresses for which a turbulent model needs to be defined. In the present study, the LES turbulence model characterized by the desired performance in the simulation of vortex flows [23][24][25] was solved using the Smagorinsky model by the length scale L = (δxδyδz) 1 3 based on the grid cell dimensions. The eddy viscosity in the method is as follows: In which c refers to the Smagorinsky constant selected by the program based in the range of 0.1 to 0.2 and e ij defines the strain tensor. The eddy viscosity ν t is combined with the dynamic viscosity using the following relationship µ = ρ(ν + ν t ) and is finally utilized for stress calculation according to Equation (9). Furthermore, according to the VOF method as one of the widely used procedures for water surface simulation regarding the Eulerian method, a unit value is assigned to the cells filled by water, a zero value is assigned to empty cells and for the cells with water content between zero and one the assigned value is determined based on the fluid quantity [51].

The Experimental Setup
The experimental set up shown in Figure 1 was used in order to develop a steady state vortex flow. The physical model constituted a vortex cylindrical tank of internal diameter 29.4 cm and height 40 cm made up of Borosilicate glass. The diameter of the intake at the bottom of the container corresponds to 3.5 cm. The discharging flow from the intake is pumped towards the inlets. The input system of the flow is contrived by two aluminum pipes on both lateral sides of the tank on each of which vertical corrugations are hatched in order to assure uniform inflow into the tank. The corrugations on the two side pipes are hatched in opposite directions such that the flow enters the tank tangentially in a clockwise direction. Each corrugation is of 2.5 mm width and 100 mm height in depth. The entry corrugations must not be located higher than the water surface in the tank in order to prevent bubbling and turbulence on the water surface (The fluid in the tank is of a constant 18 cm depth). In addition, the corrugations are hatched 7 cm above the bottom of the container to neutralize the effect of the bottom of the tank on the inflow. follows: In which c refers to the Smagorinsky constant selected by the program based in the range of 0.1 to 0.2 and ij e defines the strain tensor. The eddy viscosity t ν is combined with the dynamic viscosity using the following relationship and is finally utilized for stress calculation according to Equation (9). Furthermore, according to the VOF method as one of the widely used procedures for water surface simulation regarding the Eulerian method, a unit value is assigned to the cells filled by water, a zero value is assigned to empty cells and for the cells with water content between zero and one the assigned value is determined based on the fluid quantity [51].

The Experimental Setup
The experimental set up shown in Figure 1 was used in order to develop a steady state vortex flow. The physical model constituted a vortex cylindrical tank of internal diameter 29.4 cm and height 40 cm made up of Borosilicate glass. The diameter of the intake at the bottom of the container corresponds to 3.5 cm. The discharging flow from the intake is pumped towards the inlets. The input system of the flow is contrived by two aluminum pipes on both lateral sides of the tank on each of which vertical corrugations are hatched in order to assure uniform inflow into the tank. The corrugations on the two side pipes are hatched in opposite directions such that the flow enters the tank tangentially in a clockwise direction. Each corrugation is of 2.5 mm width and 100 mm height in depth. The entry corrugations must not be located higher than the water surface in the tank in order to prevent bubbling and turbulence on the water surface (The fluid in the tank is of a constant 18 cm depth). In addition, the corrugations are hatched 7 cm above the bottom of the container to neutralize the effect of the bottom of the tank on the inflow.

Velocity Measurement
For the purpose of preventing any perturbation in the flow the non-destructive PIV system was utilized in the horizontal planes for measurement of velocity together with a camera capable of filming 240 frames per second with resolution of 1MP (720*1280) and an f/1.8 lens. Moreover, a continuous solid-state laser with wavelength, opening degree, and thickness of 532 Nm, 45 degrees, and 2 mm respectively with adjustable power from zero to 2.15W was also used. During the tests, the plane of the laser was regulated in horizontal extensions/elevations and the camera was placed perpendicular to the laser underneath the tank. The particle trackers in the system were of the Pliolite type with dimensions of 300 to 600 micrometers and density of 1.3 g per cubic centimeter with fair laser reflection. For the values of velocity the open-source PIVlab code was employed [47]. The PIV systematic error was minimized by temporal averaging the convergence test on the tangential and radial velocity fields in 7200 consecutive frames.

Water Surface Measurement
Measurement of the water surface profile was dealt with by installation of a camera in the horizontal elevation in front of the tank with a light source (halogen lamp) placed behind the tank. At first, a picture of the static fluid without any vortices at a depth of 18 cm was taken as the initial picture. Afterwards, another picture of the tank at an arbitrary time was taken, subsequent to the inception of the test and formation of a steady state vortex flow. Taking advantage of the image processing technique and subtraction of the initial picture from the secondary one and with regard to the fact that the camera and the light source are stationary, the remaining part resulted by subtraction of the former and latter images representing the free surface of the water with great quality and accuracy.

The Procedure of the Experiments
By adjusting the flow rate (discharge) and the depth at 0.2 L per second and 18 cm respectively, a steady state vortex flow was formed over the intake. The velocity field of the vortex flow was measured on three horizontal planes located at elevations of 5 cm, 9 cm, and 13 cm from the bottom of the tank respectively using the PIV. The properties of the experiments together with the Froude number of the intake are presented in Table 1. The intake Froude number which stands as one of the most influential dimensionless parameters in definition of the vortex flow was calculated with respect to the following relationship Fr = U 0 / gd in which U 0 and d represent the velocity of the flow at the intake and the diameter of the intake respectively.

The Numerical Model
With respect to the geometry of the physical model, the numerical model was developed as a cylindrical container with diameter and height of 29.4 and 40 cm respectively. An outlet orifice with diameter of 3.5 cm was placed in the bottom of the tank for discharge of the flow. The results of the sensitivity analysis revealed that the optimal mesh dimension and particle size leading to computationally viable solutions for discretization of the flow field corresponding to the Eulerian and Lagrangian methods were 2.5 mm while consideration of smaller mesh dimensions or particles sizes would increase the computational time of the program considerably. Therefore, 900,000 fluid and boundary particles were chosen in total in the Lagrangian model. The total number of meshes in the Eulerian approach were 1,700,000, assuming identical dimensions in the three directions. Cartesian mesh was implemented in the Eulerian model and the Fraction Area/Volume Obstacle Representation (FAVOR) method was also used to model the boundary curvature [24,52]. The intended fluid in both models was water with density and kinematic viscosity of 1000 kg/m 3 and 10 −6 N.s/m 2 respectively.
Simulation of the boundary conditions in numerical modeling is of great importance. The no-slip boundary condition for the walls was implemented for the Eulerian model. Moreover, the volumetric flow rate condition was used for the inlet and outlet. The outlet was modeled with a pipe at which the flow rate was specified in a section 5 cm below the bottom of the tank. The entries were also of rectangular cross sections and dimensions of 2.5 mm × 100 mm from the elevation of 7 cm to 17 cm from the bottom of the container To define the inflow and outflow in the Lagrangian method the open boundary condition was considered. The discharging flow rate from the outlet was equal to the summation of the input flow rates and constant velocity was assumed in the outlet and inlet. The density distribution in the buffer layer for the inlets was considered constant with regards to the assumption of the hydrostatic pressure while the density of the buffer particles in the intake was calculated by interpolation from the surrounding particles. The buffer layer particles were classified into eight layers to assure full Kernel support for the particles in the vicinity of the buffer layer. An h/dp = 2 ratio between the smoothing length and the initial particle spacing was considered in the SPH model. Moreover, the runtime of the program was regarded as 100 s to reassure formation of a steady vortex flow and the outputs were saved every 0.5 s.
The SPH model was executed using parallel processing on the NVIDIA GTX 1080 graphic card with memory of 8 GB, 2560 cuda cores and performance of 8228 GFlops while the Eulerian method was implemented on the Intel core i7 processor with memory of 8 GB and performance of 230.4 GFlops.

The Air Core Vortex
The developed vortex flow in both of the numerical simulations together with the physical model are presented in Figure 2 indicating the sought performance of both approaches in the simulation of the air-core vortices. According to the results obtained from the numerical models, it was observed that both models reported lower depths for the air core. In comparison with the experimental data the Lagrangian model led to even lower air core depths compared to the Eulerian method. However, more accurate investigation of the obtained results concerning the performance of the two approaches is presented in the following.

The Tangential Velocity Distribution
The tangential velocity distribution corresponding to the physical model and numerical simulations at elevations of 5 cm, 9 cm, and 13 cm from the bottom of the tank are presented in Figure 3. The presented graph was made dimensionless by incorporation of the inductive tangential velocity with the inflow (Vin) and the radius of the intake. It should be pointed out that in the physical model the velocity field close to the core could not be assessed due to the light scattering in that region. According to Figure 3 it can be observed that the Eulerian method performed more efficiently in presentation of the tangential velocity component. The maximum difference between the results of the numerical model and the experimental data corresponds to 7.4% and 23.3% for the Eulerian and Lagrangian approach respectively.

The Tangential Velocity Distribution
The tangential velocity distribution corresponding to the physical model and numerical simulations at elevations of 5 cm, 9 cm, and 13 cm from the bottom of the tank are presented in Figure 3. The presented graph was made dimensionless by incorporation of the inductive tangential velocity with the inflow (V in ) and the radius of the intake. It should be pointed out that in the physical model the velocity field close to the core could not be assessed due to the light scattering in that region. According to Figure 3 it can be observed that the Eulerian method performed more efficiently in presentation of the tangential velocity component. The maximum difference between the results of the numerical model and the experimental data corresponds to 7.4% and 23.3% for the Eulerian and Lagrangian approach respectively. the inductive tangential velocity with the inflow (Vin) and the radius of the intake. It should be pointed out that in the physical model the velocity field close to the core could not be assessed due to the light scattering in that region. According to Figure 3 it can be observed that the Eulerian method performed more efficiently in presentation of the tangential velocity component. The maximum difference between the results of the numerical model and the experimental data corresponds to 7.4% and 23.3% for the Eulerian and Lagrangian approach respectively.

The Radial Velocity Distribution
The radial velocity distribution attained from the Eulerian and Lagrangian simulations and the laboratory experiments are compared in Figure 4. It is evident that the radial velocity component in the reaches far from the vortex core are of negligible values while increasing on moving towards the core with an upward trend. The compatibility of the numerical and experimental results is evaluated as favorable in terms of the trivial radial velocities. Furthermore, the observed oscillations in the radial velocity distribution from the laboratory experiments can be interpreted to be in conjunction with either laboratory errors or the existence of secondary flows.

The Free Surface
The results of the numerical and experimental outcomes regarding the water free surface are presented in Figure 5. In the presented graph, the elevation of the free surface (H r ) has been made dimensionless relative to the water elevation over the intake axis (H 0 ) and the depth of the fluid in the vicinity of the borders (H). Comparison of the numerical simulations with the experimental results demonstrated that although the error with regard to the depth of the air core is high, the trends of both methods are similar to the experimental results with the Lagrangian model following more accurately the free surface curvature.
tions and the laboratory experiments are compared in Figure 4. It is evident that the radial velocity component in the reaches far from the vortex core are of negligible values while increasing on moving towards the core with an upward trend. The compatibility of the numerical and experimental results is evaluated as favorable in terms of the trivial radial velocities. Furthermore, the observed oscillations in the radial velocity distribution from the laboratory experiments can be interpreted to be in conjunction with either laboratory errors or the existence of secondary flows.

The Free Surface
The results of the numerical and experimental outcomes regarding the water free surface are presented in Figure 5. In the presented graph, the elevation of the free surface (Hr) has been made dimensionless relative to the water elevation over the intake axis (H0) and the depth of the fluid in the vicinity of the borders (H). Comparison of the numerical simulations with the experimental results demonstrated that although the error with regard to the depth of the air core is high, the trends of both methods are similar to the experimental results with the Lagrangian model following more accurately the free surface curvature.

Sensitivity Analysis
For both the Lagrangian and Eulerian approaches, sensitivity analysis was carried out for different mesh or particle sizes. The obtained results indicated that for particle size of 3.5 mm in the SPH model, no vortex core or even dimple was formed on the free surface and the tangential velocity distribution was not correct. According to Figure 6 it is evident that decreasing the mesh dimension or particle size to 2 mm results in increment of the peak. For mesh dimension or particle size smaller than 2.5 mm the velocity profiles are close for r/rintake > 2.

Sensitivity Analysis
For both the Lagrangian and Eulerian approaches, sensitivity analysis was carried out for different mesh or particle sizes. The obtained results indicated that for particle size of 3.5 mm in the SPH model, no vortex core or even dimple was formed on the free surface and the tangential velocity distribution was not correct. According to Figure 6 it is evident that decreasing the mesh dimension or particle size to 2 mm results in increment of the peak. For mesh dimension or particle size smaller than 2.5 mm the velocity profiles are close for r/r intake > 2. els.

Sensitivity Analysis
For both the Lagrangian and Eulerian approaches, sensitivity analysis was carried out for different mesh or particle sizes. The obtained results indicated that for particle size of 3.5 mm in the SPH model, no vortex core or even dimple was formed on the free surface and the tangential velocity distribution was not correct. According to Figure 6 it is evident that decreasing the mesh dimension or particle size to 2 mm results in increment of the peak. For mesh dimension or particle size smaller than 2.5 mm the velocity profiles are close for r/rintake > 2. Besides attribution of the correct dimensions to the grid and particle dimensions, the computational cost is also another controversial factor needed to be considered. Therefore, the runtime and the maximum difference between the results of the numerical model and the experimental data (error) corresponding to each of the numerical models for different grid and particle sizes are presented in Table 2. The lower the mesh and particle size, the higher is the difference between the two models. The fact that the Lagrangian model takes up to four times longer to run on a much more powerful processor compared to the Besides attribution of the correct dimensions to the grid and particle dimensions, the computational cost is also another controversial factor needed to be considered. Therefore, the runtime and the maximum difference between the results of the numerical model and the experimental data (error) corresponding to each of the numerical models for different grid and particle sizes are presented in Table 2. The lower the mesh and particle size, the higher is the difference between the two models. The fact that the Lagrangian model takes up to four times longer to run on a much more powerful processor compared to the Eulerian model needs to be considered although the Eulerian model shows more accurate results. Hence, the efficiency of the Lagrangian method needs to be clarified prior to making decisions about which approach to take on board.

Conclusions
Formation of air-core vortices on power plant intakes has always been considered as a high-risk issue. The deeper the cognition of the phenomenon, the better the relevant potential hazards can be handled. One of the most economically viable methods for investigation of the vortex flow field implies application of numerical simulations. Two main approaches namely the Lagrangian and Eulerian approaches are the most widely applied methods. Unlike the Lagrangian approach, numerous investigations based on the Eulerian approach have been conducted for vortex flow while the application of the Lagrangian method has not received attention. The present study was performed aiming to provide a sensible comparison between the two methods in a complicated flow condition such as vortex flow. According to the obtained results it was concluded that both the Lagrangian and the Eulerian approaches are capable of simulation of steady vortex flow. However, the Lagrangian method was found to be of much greater computational cost with less accuracy.