High-Precision Numerical Scheme for Vortical Flow

In this study, a new high-precision numerical simulation scheme for vortical flows (vortex-based scheme) is proposed. This scheme identifies a vortical flow in each computational cell, and then, reconstructs a vortical velocity distribution based on the Burgers vortex model. In addition, a pressure distribution in the vicinity of the vortex center is also reconstructed. The momentum transfer is calculated with the reconstructed velocity and pressure distributions, and therefore, the vortex-based scheme can simulate vortical flows more accurately than the conventional schemes. In fact, as the simulation result of inviscid vortex attenuation problem, the vortex-based scheme shows lower simulation error compared to the conventional discretization schemes. Moreover, also in the numerical simulation of the quasi-steady vortical flow, the simulation accuracy of the vortex-based scheme is superior to those of the conventional schemes.


Introduction
Vortex cavitation can be observed in various engineering scenes, e.g.pump sump.However, in usual, the occurrence of the vortex cavitation is not favorable because the vibration of structural components and/or the noise can be induced by the vortex cavitation.Moreover, the cavitation bubbles generated at the vortex core may damage a structural surface when they collapse.Also in sodiumcooled fast reactors [1], the sub-surface vortex cavitation might be induced by high-speed suction flow into an outlet pipe in the reactor vessel.However, the onset condition of the vortex cavitation cannot be clarified easily because the vortex cavitation shows highly complicated behaviors associated with phase change.Therefore, the authors have conducted experimental and numerical studies to evaluate the onset condition of the vortex cavitation.As the experimental study, a fundamental experiment is conducted to investigate the influence of the dynamic viscosity on the onset condition [2], and some scale model tests are also conducted to check the dependency on the local structural geometry and the effect of countermeasure obstacles [3].On the other hand, in the numerical study, a high-precision simulation algorithm for sub-surface vortex flows is under development, in which a unstructured mesh scheme is employed to model accurately the local structural geometry near a sub-surface vortex.
The high-precision numerical simulation algorithm is developed originally for the evaluation of gas entrainment phenomena in sodium-cooled fast reactors [4][5][6].Currently, it is confirmed this algorithm can simulate the gas entrainment phenomena caused by free surface vortices, that is, dynamic vortex behaviors, i.e. vortex development, movement and attenuation, can be simulated.Therefore, the algorithm is considered to be applicable to numerical simulations of the vortex cavitation.However, the numerical simulation of the vortex cavitation is not easy because the sub-surface vortex has highly thin vortex core compared to the free surface vortex.In other words, it is very difficult to reproduce the velocity distribution at the vicinity of the core of the sub-surface vortex.In this case, a very fine mesh at the vicinity of the vortex core and/or a high-precision numerical simulation scheme should be employed to simulate the sub-surface vortex accurately.A very fine mesh is easy to be constructed but requires high computational cost which may exceeds the limit of usable computational resource.Therefore, a highprecision numerical simulation scheme is necessary to simulate the sub-surface vortex accurately.Here, it should be noted the implementation of well-known high order (more than 2nd order) schemes, e.g.K-K scheme [7] are difficult on unstructured meshes because of irregular computational cell arrangement.
In this paper, a high-precision simulation scheme for vortical flows on unstructured meshes is presented.In conventional schemes, the velocity distribution between computational cells is interpolated with the constant, linear or higher order functions based on the discrete values defined at the cells, and therefore, the existence of a vortical flow is considered indirectly through the velocity distribution around the vortex.On the other hand, in our new high-precision scheme, a vortical flow is identified in each cell and the vortical velocity distribution at the vicinity of the vortex center is reconstructed locally based on the Burgers vortex theory [8].Namely, the existence of a vortical flow is considered directly in the new high-precision scheme.In this sense, this new scheme is called "the vortex-based scheme" in this paper.In addition to the velocity distribution, the pressure distribution at the vicinity of the vortex center is also considered based on the Burgers vortex theory to simulate accurately the mechanical balance between pressure gradient and centrifugal force.Then, the calculations of momentum transport through cell faces are performed in consideration of the reconstructed velocity and pressure distributions.The superiority of the vortex-based scheme to the conventional schemes is confirmed by the fundamental verifications, such as the numerical simulation of vortex attenuation.

Formulation of High-Precision Numerical
Scheme for Vortical Flow (Vortex-Based Scheme)

The Identification of a Vortical Flow in Each Computational Cell
In the vortex-based scheme, the existence of a vortical flow is considered directly.Therefore, in the first place of the calculation, vortical flows are identified by applying the discriminant [9,10] to each computational cell.Namely, a vortex center exists in a cell when the following discriminant is satisfied in the cell.
where Q and R are the second and third invariants of the velocity gradient tensor In the cell with a vortex center, the direction of the vortex, i.e. the direction of the rotational axis, is obtained as the eigen vector of the eigenvalue equation of the velocity gradient tensor.Then, a plane is defined, which is normal to the direction of the vortex center and passes through the cell center.On this plane, the point with zero in-plane velocity is determined as the vortex center.

The Supplement of the Vortical Velocity Distribution
A strong vortical flow has very large velocity gradient in the immediate vicinity of the vortex center.Therefore, it is highly difficult to simulate such a velocity distribution accurately with limited computer resources (lack of sufficient number of cells).In the vortex-based scheme, the vortical velocity distribution is supplemented in the vicinity of the vortex center by utilizing the Burgers vortex model which is known as an excellent model of the vortices in nature, e.g.free surface vortices, sub-surface vortices or turbulent vortices.
The circumferential velocity   u  distribution of the Burgers vortex model is written as where r is the radial coordinate with the origin at the vortex center.Equation (4) has two parameters, i.e. the circulation Γ and the specific radius r 0 .In the vortexbased scheme, these parameters are determined to make the circumferential velocity distribution consistent with the local velocity distribution.Namely, upon a given (simulated) velocity field, two parameters are adjusted to minimize the difference in the circumferential velocity distributions of Equation ( 4) and the given field at the cells in the vicinity of a vortex center.The calculation procedure is shown in Figure 1 as follows: 1) a cell with a vortex center (called "vortex cell" in this paper) is selected; 2) vertices on the vortex cell and the cells around the vortex cell is selected; 3) when the velocity is given at the cell centers, the circumferential velocities at the vertices   v u  are inter- polated from the cell center velocities   c u  as follows: where the summations are operated on all cells surrounding each selected vertex, and w c is the weighting factor for the summation; 4) the differences between the Burgers model (Equation (4)) and the circumferential velocities at the selected vertices   v u   are calculated and averaged as the form of the weighted root mean square (shown by an overbar) as follows:

Vortex center
Vortex cell where the summation is operated on all selected vertices, and w v is the weighting factor for the summation; 5) the circulation and the specific radius are adjusted iteratively to minimize the averaged differences   The weighting factor in Equation ( 5) is determined by minimizing the cost-function (Co in Equation ( 7)) under the constraint formulated as Equation ( 8) [11], and the weighting factor in Equation ( 6) is determined in a similar manner.
in the numerical simulations on unstructured meshes.In he cell based FVM is employed on an unstructu where cv is the vector joining a cell center to a vertex.
In Figure 1, the calculation procedure on a two-dimensional structured mesh is shown for simplification.However, the calculation on a three-dimensional unstructured mesh can be conducted in the same procedure.r

The Momentum Transport Calculation
The finite volume method (FVM) is employed frequently that case, each cell is generally selected as the control volumes for discretization.Also in our high-precision numerical simulation algorithm, such simulation method (cell based FVM) is employed, and therefore, the following description is given based on the cell based FVM.However, the concept of the vortex-based scheme is applicable also to other methods, e.g. the finite difference method.When t red mesh, it is difficult to formulate high order momentum transport schemes, i.e. third or higher order schemes, due to the irregular cell arrangement.Therefore, the first and second order upwind schemes are used frequently for the momentum transport calculation.In the first order upwind scheme, the momentum on a cell face is determined to be same as that in the upwind cell of the cell face.On the other hand, in the second order upwind scheme, the momentum on a cell face   where is the vector joining the upwi cf r fac nd cell center to the cell e center.In Equation ( 9), the momentum gradient can be calculated by the Gauss-Green method [11] or Least-square method [12].The first or second order upwind scheme works well in the numerical simulations of smooth flows without discontinuity (shock) and strong vortices.However, as mentioned briefly in Section 2.2, the velocity distribution in the vicinity of a vortex center is very sharp and a very high resolution mesh is necessary to simulate the vortical velocity distribution accurately by the first or second order upwind scheme.Instead of such a very high resolution mesh, which is not applicable to practical simulations, the vortex-based scheme supplements the vortical velocity distribution in the vicinity of a vortex center by the Burgers vortex model (Equation ( 4)) and calculates the momentum transport based on the supplemented velocity distribution.Namely, when the velocity distribution in the vicinity of a vortex center is given by Equation ( 4), the momentum in the cir- where f  is the fluid density on the cell face and f u  is the c mferential velocity at the cell face center c culated by Equation (4).Then, the momentum transport is calculated based on Equation (10).

ircu al-
In the vicinity of a vortex center, not only the velocity gradient but also the pressure gradient is large due to the mechanical balance between the circumferential velocity and the pressure gradient (see Equation (11)).Therefore, the pressure should be supplemented to simulate vortical flows accurately.In the vortex-based scheme, the pres-sure distribution in the vicinity of a vortex center is calculated in consideration of the mechanical balance equation of the circumferential velocity and the pressure gradient, written as where p is the pressure.Substituting the circumferential velocity distribution of the Burges model (Equation ( 4)) into Equation ( 11) leads the pressure distribution. (12)  where 0 R r r  .The pressure calculated by Equation ( 12) is applied to the faces on a vortex cell and the faces on the cells around the vortex cell.Then, the pressure gradient term in the momentum transport equation (Navier-Stokes equation) is calculated at each cell based on the applied pressure values.

The Reproducibility of a Vortical Flow
As the first verification, the reproducibility of a vortex by the vortex-based scheme is addressed.In this calculation, the Burgers vortex is set on two-dimensional structured meshes, and then, the identification and the velocity supplement calculations are conducted to check the reproducibility of the vortex.As shown in Figure 2, the simulation domain is 1.00 × 1.00 and subdivided into 8 × 8 square cells.The velocity distribution of the Burgers vortex is applied to all cells with various vortex center positions, specific radii and/or circulation values.
Table 1 shows the calculation results.It is evident that the vortex center position and the circulation are well re-  r position, circulation and/or speific radius of the original Burgers vortex.As for the cal-produced in all calculation cases regardless of the variations of the vortex cente c culated specific radius, almost the same values are evaluated regardless of the variations of the vortex center position and/or circulation.However, the calculation accuracy (reproducibility) is affected significantly by the specific radius of the original Burgers vortex.In other words, the calculation error is about 25% when the original specific radius is the same as the cell size (0.125), and the error is reduced to about 6.8% when the original specific radius is twice the cell size.In practical simulations, this kind of local (only in the vicinity of a vortex) calculation error may be negligible even if it reaches about 25%, because much larger error sources can exist in such simulations.However, the calculation results show that the mesh resolution should be at least half of the specific radius of a vortex to perform accurate simulations.This constrain may seem very difficult to be satisfied.However, the authors' previous research shows that the cell size has to be about 1 20 of the specific radius to reproduce the vortical velocity distribution when the second order upwind scheme is employed [4].Therefore, the necessary number of cells for accurate simulations can be reduced significantly by employing the vortex-based scheme.

The Simulation of Attenuation Behaviors of an Inviscid Vortex
, viscosity de ation errors of each scheme.
tex ted with four es, i.e. the first order upw d order upwi ird order M and vortex-b schemes, and th ulation the pe It is well known that the momentum transport schemes e.g. the first order upwind scheme, have the numerical pendent on trunc Therefore, even in the numerical simulations of an inviscid vortex, the vortical velocity distribution decays temporally and the total kinetic energy in the simulation domain is not conserved.In this sense, the conservativeness of the total kinetic energy can be an indicator of the influence of the numerical viscosity, that is, the simulation accuracy.Here, attenuation behaviors of an inviscid vor-compared in terms of the conservativeness of the total kinetic energy.The numerical simulations are conducted on 1.00 × 1.00 two-dimensional domain subdivided into square cells (shown in Figure 3).As for the mesh resolution, four structured meshes, i.e. 8 × 8, 16 × 16, 32 × 32 and 64 × 64 square cells, are employed.At the initial state, the velocity distribution of the Burgers vortex with the circulation of 0.100 and the specific radius of 0.050 is applied to the meshes, and then, the simulations of 100 time-marching with the time increment of 0.01 are conducted.The temporal variations of the velocity distributions and the total kinetic energies are investigated in each simulation for the comparison of four schemes.
In Figure 4, the temporal variations of the velocity distributions on the 32 × 32 mesh are shown.The simulation result with the first order upwind scheme (Figure 4(a)) shows rapid decay of the initial vortical velocity distribution.When the second order upwind schemes is employed, such a decay behavior is suppressed and are simula schem ind, secon ased nd, th e sim USCL results are ak of the circumferential velocity in the vicinity of the vortex center is maintained appreciably even at t = 0.40.The simulation result with the third order MUSCL scheme is almost the same as the result with the second order upwind scheme.The vortex-based scheme shows superior simulation accuracy to these three schemes.The decay of the circumferential velocity distribution in the vicinity of the vortex center is apparently small compared to the simulation result with the second order upwind scheme.The temporal variations of the total kinetic energies on the 32 × 32 mesh are shown in Figure 5.In the simulation with the first order upwind scheme, the total kinetic energy decreases very rapidly along with the rapid decay of the vortical velocity distribution as shown in Figure 4(a).The second order upwind scheme highly improves the total energy conservation compared to the  first order upwind scheme.The loss of the total kinetic energy at t = 1.0 is about 30% of its original value (t = 0.0), which is about 48% in the simulation with the first order upwind scheme.The third order MUSCL scheme shows slightly higher total kinetic energy conservation than the second order upwind scheme, and the loss of the total kinetic energy at t = 1.0 is about 26%.However, the temporal variation of the total kinetic energy is similar to that in the simulation results with the second order upwind scheme.Compared to the other three schemes, the l kinetic energy t t = 1.0 is about 17%.Therefore, the vortex-based mentation (Section 2.4) is important as well as the velocity supplementation (Section 2.3), because the loss of the total kinetic energy at t = 1.0 increases to about 23% when the pressure supplementation is not employed.
In the vortex-based scheme, two parameters, i.e. circulation and specific radius, dominate the vortical velocity distribution.Here, the temporal behaviors of these two parameters are investigated as shown in Figure 6.When a relatively coarse mesh (16 × 16 mesh) is employed, th specific radius at the initial state is evaluated as a consitex and the specific radius becomes larger with time.On the other hand, the circulation is evaluated accurately at vortex-based scheme shows much better total energy conservation.In fact, the loss of the tota derably larger value than that of the original Burgers vora scheme has a superior ability to simulate vortical flows accurately.It should be noted that the pressure supple-the initial state and does not decrease significantly with time.When a finer mesh (32 × 32 mesh) is employed, the ) circulation.
evaluation results of these two parameters changes completely.In the finer mesh case, the specific radius is evaluated accurately throughout the simulation because the enhanced mesh resolution improves the reproducibility of the vortical flow (as shown in the previous section).Therefore, the loss of the kinetic energy is induced mainly by the change in the circulation.In fact, the circulation decreases temporally with the decay of the peak circumferential velocity.These evaluation results of the specific radius and the circulation show the characteristics of the vortex-based scheme can b anged by mesh resolution.In general, the mesh resolution has a significant effect on the simulation accuracy.Figure 7 shows the loss of the total kinetic energy (ratio from its original value as the functions of the cell size Δ) in the simulations results on 8 × 8, 16 × 16, 32 × 32 and 64 × 64 meshes.The simulation accuracy is highly enhanced with the mesh resolution when the first order upwind, second order upwind or third order MUSCL schemes is employed.On the other hand, the vortex-based scheme gives a superior simulation accuracy of vortical flows on both the coarse and refined meshes.In other words, the vortex-based scheme shows roughly the same simulation accuracy both on the coarse and fin eshes.Therefor he vortex-based scheme is especi useful in the simulations s.

Numerical Simulation of Quasi-Steady Vortex
To check the applicability of the vortex-based scheme to the numerical simulation on a three-dimensional unstructured mesh, the quasi-steady vortex is simulated in reference to the experiment by Moriya [13].lindrical tank and the outlet nozzle is installed on the bottom of the cylindrical tank.The uniform flow through the inlet slit generates a vortical flow in the cylindrical tank and the strength of this vortex is enhanced by the downward flow towards the outlet nozzle.The working fluid is water at room temperature.Two cases of unsteady numerical simulations are conducted with the second order upwind scheme and with the vortex-based scheme until quasi-steady states are obtained.As shown in Figure 8, the mesh resolution is rather low, that is, the horizontal cell size at the center of the cylindrical tank is about 5.0 mm, which is not small enough to the measured specifi gasquid inte ental aparatus is not considered, and only the comparison of the simulation results is conducted.All structural walls are modeled as non-slip walls.As the other boundary conditions, a uniform inlet velocity with 8.33 × 10 −4 m 3 /s is applied to the inlet a pressure condition is applied to the outlet.
Figure 9 shows the horizontal velocity distribut t the height of 0.15 m from the bottom of the cylindrical tank.It is evident that the vor scheme gives larger velocity in the vicinity of center compared to the secon fference in th simulation can be observ-c radius (about 6.0 mm).In this study rface located in the Moriya's experim

Conclusion
In this paper, the vortex-based scheme is developed to achieve high-precision numerical simulations of strong vortical flows which may cause the vortex cavitation.In flow is identified in each computational cell, and then, the Burgers vortex model is applied to supplement a vortical velocity distribution in the vicinity of the vortex center.The pressure distribution in the vicinity of the vortex center is also supplemented.The momentum transfer is calculated with those supplemented velocity and pressure distributions.The simulation result of the inviscid vortex attenuation problem shows the vortex-based scheme leads lower si mulation err d order schethe vortexbased scheme can simulate vortical flows more accurately than the conventional schemes.The higher simulation accuracy of the vortex-based scheme than the conventional scheme is shown also in the numerical simulation of the quasi-steady vortex experiment.The vortexbased scheme succeeds in simulating the vortical flow on a three-dimensional unstructured mesh more accurately than the second order upwind scheme.This result implies that the vortex-based scheme is applicable to practical simulations.

Figure 1 .
Figure 1.Procedure to determine vortex parameters: (a) selection of a vortex cell; (b) listing up of vertices; (c) interpolation of circumferential velocity at each vertex (red vectors); (d) calculation of differences at each vertex (Burgers vortex velocity in horizontal is shown by blue curve) and (e) adjustments of parameters.
min the upwin cell of the cell d(9) face:

Figure 4 .
Figure 4. Temporal variations of velocity distributions on 32 × 32 mesh: (a) first order upwind scheme; (b) second order upwind scheme and (c) vortex-based scheme.
Figure 5. Simulation mesh and initial velocity distribution on 32 × 32 mesh.

Figure 6 .
Figure 6.Temporal behaviors of vortex parameters: (a) specific radius and (b
cheme.The di results with these two schemes s e ed clearly in the comparison of the circumferential velocity distributions (shown in Figure10).The vortex-based scheme gives larger peak of the circumferential velocity and smaller specific radius.Namely, the peak of the vortical velocity in the vicinity of the vortex center is simulated more sharply by employing the vortex-based scheme instead of the second order upwind scheme.Therefore, it is confirmed that the vortex-based scheme can simulate vortical flows accurately even on three-dimensional unstructured meshes.