VIScous Vorticity Equation (VISVE) for Turbulent 2-D Flows with Variable Density and Viscosity

: The general vorticity equation for turbulent compressible 2-D ﬂows with variable viscosity is derived, based on the Reynolds-Averaged Navier-Stokes (RANS) equations, and simpliﬁed versions of it are presented in the case of turbulent or cavitating ﬂows around 2-D hydrofoils.


Introduction
The vorticity equation has been utilized by several authors in the past to analyze the viscous flow around bodies. Vortex element (or particle, or blub) and vortex-in-cell methods have been used for several decades for the analysis of 2-D or 3-D flows, as described in Chorin [1], Christiansen [2], Leonard [3], Koumoutsakos et al. [4], Ould-Salihi et al. [5], Ploumhans et al. [6], Cottet and Poncet [7], and Cottet and Poncet [8]. Those methods essentially decouple the vortex dynamics (convection and stretching) from the effects of viscosity.
In recent years, the VIScous Vorticity Equation (VISVE), has been solved by using a finite volume method, without decoupling the vorticity dynamics from the effects of viscosity. This method has been applied to 2-D and 3-D hydrofoils, cylinders, as well as propeller blades as described in Tian [9], Tian and Kinnas [10], Wu et al. [11], Wu and Kinnas [12], Wu et al. [13], Wu and Kinnas [14]. The major advantage of this approach is that it requires a significantly smaller computational domain than RANS over which VISVE must be solved, as shown in Figure 1, due to the fact that the vorticity vector vanishes much closer to the body surface (in the order of 1-3 maximum body thickness), as opposed to the velocity vector which vanishes much farther (in the order of 5-10 body lengths) from the body surface. Due to the significantly smaller domain, VISVE requires a much smaller number of cells (5-50 times smaller) than RANS for the same accuracy of the predictions, and subsequently, a significantly smaller CPU time, as reported in Wu et al. [11].
All applications of VISVE mentioned above have addressed laminar flow of an incompressible fluid. However, in the case of incompressible turbulent flow the vorticity equation must be modified since the viscosity (molecular + turbulent) varies within the flow. In addition, in the case of cavitating flow, treated via a mixture model, both the density and the viscosity vary within the flow. In this work the general equation in terms of the mean vorticity is derived in the case of turbulent flows with variable density and viscosity, and then simplified in the case of turbulent or cavitating flows around 2-D hydrofoils.

Review of Reynolds-Averaged Navier-Stokes Equations
We consider flow of a fluid, where q = (u 1 , u 2 , u 3 ) is the mean velocity in a x 1 , x 2 , x 3 coordinate system. The Navier-Stokes equations are as follows: where ρ is the density of the fluid, p is the mean pressure, Φ is the potential of a conservative body force per unit volume (Φ = ρgz in the case the body force is the gravity, with the gravity acceleration g pointing in the -z direction.), and T is the (symmetric) tensor of viscous stresses where µ = µ m + µ t is the total dynamic viscosity, with µ m being the molecular viscosity and µ t being the turbulent viscosity of the fluid, under the Boussinesq approximation for the turbulent stresses; δ ij is the Kronecker delta; k is the turbulent kinetic energy; N = 2 for two-dimensional flows, and N = 3 for three-dimensional flows. The second term in Equation (2) guarantees that τ ii = −ρu i u i = −2ρk, where u i is the turbulent velocity in the i th direction. Equation (1) can also be written as follows, as shown in the Appendix A: where ω = (ω 1 , ω 2 , ω 3 ) = ∇ × q is the vorticity of the flow, and the tensor ∇ q is the gradient of q, defined as follows With the help of some identities of vector calculus, as shown in the Appendix A, Equation (3) becomes: Dividing Equation (4) by ρ and rearranging will give us: where ν = µ ρ is the kinematic (total) viscosity

Flow with Constant Density
In that case, since ∇ · q = 0 and ρ = constant (In this work we consider incompressible flow when ρ = constant, as opposed to the more general definition of incompressible flow when Dρ Dt = 0), Equation (5) becomes: In the case of Newtonian fluid in laminar flow (µ = µ m =constant and k = 0) and incompressible flow, in which case ∇ · q = 0, Equation (4) reduces to its most common form:

Vorticity Equation in 2D
We consider 2D flow, in which case x 1 = x, x 2 = y, x 3 = z; u 1 = u, u 2 = v, u 3 = w; ∂ ∂y = 0 and v = 0. The velocity has two components: q = (u, w), and the vorticity has only one component in the y direction ω 2 = ω = ∂u ∂z − ∂w ∂x . As shown in the Appendix A the vorticity equation will become as follows: As shown in the Appendix A Equation (8) may also be written as follows:

2D Flow of Fluid with Constant Density
In that case ρ = constant, ∇ρ = 0, ∇ 1 ρ = 0, and ∇ · q = 0. Then Equation (9) reduces to: · ∇u (10) or, by expressing ω = ∂u ∂z − ∂w ∂x in the second term on the RHS of Equation (10), and by making use of the continuity equation, ∂u ∂x + ∂w ∂z = 0, we get: It is worth noting that Equation (10) is valid for turbulent flows (with the vorticity ω being the mean vorticity), but the turbulent kinetic energy k is not involved in the equation.
In the case of laminar flow of Newtonian fluid ν = µ m ρ = constant, Equation (10) becomes the vorticity equation in its most common diffusion equation form: or, since ∇ · q = 0:

2D Flow around Hydrofoil
In that case we assume that the hydrofoil is placed along axis x, with the inflow also along x axis, and that ∂ ∂x << ∂ ∂z , especially within the narrow region close to the hydrofoil and its wake where ω = 0. Then, as shown in the Appendix A, Equation (9) reduces to:

Conclusions and Future Work
The vorticity equation (in terms of the mean vorticity) is derived in the case of turbulent flows with variable density and viscosity, and its simplified version in the case of flows around 2-D hydrofoils has been provided. Equation (14) has already been implemented by Yao and Kinnas [15] to address the turbulent non-cavitating flow around 2-D hydrofoils and cylinders, as well as the cavitating laminar flow around 2-D hydrofoils, using the mixture model, by Xing et al. [16].
Representative results from the work of Yao and Kinnas [15] are shown in Figures 2-4 together with results from RANS by using the commercial code ANSYS/Fluent (https://www.ansys.com/ products/fluids/ansys-fluent). In this case Equation (14) was used where the turbulent viscosity was evaluated by synchronous coupling of VISVE with the open source RANS code OpenFOAM (https://www.openfoam.com/). As described in more detail in Yao and Kinnas [15], at every time step the velocities determined by VISVE were passed into the part of OpenFOAM which solved the k, equations and then returned the turbulent kinematic viscosity back to VISVE, in order to solve Equation (14) with the updated values of the kinematic viscosity. In the case of the hydrofoil the velocity profiles predicted by VISVE, shown in Figure 3, are in good agreement to those predicted by RANS. It should be noted that the hydrofoil assumption has also been used in the case of the cylinder, for which results are shown in Figure 4. The results from VISVE compare well with those from RANS, even though with somewhat larger differences downstream of the cylinder, up to the time shown in the figure, but deteriorate for later times as shown in Yao and Kinnas [15], once asymmetry appears between the top and bottom flow (not shown in this paper).
In the future, the author and his students intend to assess numerically the effect of the last three terms in the right-hand side of Equation (10), on the prediction. These terms are currently ignored, but they might affect the accuracy of predictions, especially in the case of a hydrofoil at high angles of attack where the hydrofoil assumptions, made in this paper, would not hold. The ultimate objective of this research is to extend VISVE in the case of turbulent flows around 3-D hydrofoils, and eventually propeller blades.

Conflicts of Interest:
The author declares no conflict of interest.

Appendix A. Proofs
Appendix A.1. How to Get from Equation (1) to Equation (3) We consider the i th component of Equation (1): with the ∂τ ji ∂x j term meant in the Einstein notation with j = 1, 2, 3 We then consider the rate of strain term of the expression for τ ij , as given by Equation (2): where kij is the Levi-Civita symbol. Then for given i and with j = 1, 2, 3, using the Einstein notation we have: Based on Equations (2) and (A3) we get: which eventually leads to Equation (3).
Appendix A.2. How to Get from Equation (3) to Equation (4) Using vector identities we get: Using the equations above we get to Equation (4).
Appendix A.3. How to Get from Equation (5) to Equation (8) We consider 2D flow, in which case x 1 = x, x 2 = y, x 3 = z; u 1 = u, u 2 = v, u 3 = w; ∂ ∂y = 0 and v = 0. The velocity has two components: q = (u, w), and the vorticity has only one component in the Then by using the identity: due to the identity: ∇ × ∇ f = 0 Then using the following identities: We get: Now, in 2-D, since the vortex stretching term ( ω · ∇) q = 0, the above equation becomes: We now take the ∇× of each term on the RHS of Equation (5): • The ∇× of the 1st term on the RHS of (5): The ∇× of the last 3 terms of Equation (5) are treated in a similar manner as that of the first term, and the corresponding resulting terms are the last 3 terms of Equation (8). Please note the first and the last 3 terms in Equation (8) vanish in the case ρ = constant. • The ∇× of the 2nd term on the RHS of (5): ∇ × (ν∇ 2 q) = ∇ν × ∇ 2 q + ν∇ 2 ( ∇ × q) = ∇ν × ∇ 2 q + ν∇ 2 ω (A18) • The ∇× of the 3rd term on the RHS of (5): • The ∇× of the 4th term (divided by 2) on the RHS of (5): We only handle this term in 2D. In that case: • The ∇× of the 5th term on the RHS of (5): with ∇ · ω = 0, and since ω · ∇ = 0 in 2D we get: In addition, in 2D: The combination of the ∇× of the 4th and the 5th terms will give us: The last two terms above can also be rewritten by using µ = νρ and ∇(νρ) = ν ∇ρ + ρ ∇ν, which renders the RHS of Equation (8) Appendix A.4. How to Get from Equation (8) where we have made use of ∂u ∂x + ∂w ∂z = ∇ · q We finally get to Equation (9) by considering the identity: