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

#### **Spyros A. Kinnas**

Ocean Engineering Group, Department of Civil, Architectural, and Environmental Engineering, The University of Texas at Austin, Austin, TX 78712, USA; kinnas@mail.utexas.edu; Tel.: +1-512-475-7969

Received: 23 January 2020; Accepted: 4 March 2020; Published: 11 March 2020

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

**Keywords:** vorticity; viscous flow; turbulent flow; mixture model; cavitating flow

#### **1. 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.

**Figure 1.** Typical domain and grids for flow around hydrofoil, from RANS with 150K cells (**left**) and from VISVE with 28K cells (**right**). Note the image on the right is drawn at a larger scale than that on the left.

#### **2. 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*<sup>3</sup> coordinate system. The Navier-Stokes equations are as follows:

$$
\rho \frac{\partial \vec{q}}{\partial t} + \rho (\vec{q} \cdot \nabla) \vec{q} = -\nabla p + \nabla \cdot \mathbf{T} - \nabla \Phi \tag{1}
$$

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

$$\mathbf{T} = \begin{bmatrix} \tau\_{11} & \tau\_{12} & \tau\_{13} \\ \tau\_{21} & \tau\_{22} & \tau\_{23} \\ \tau\_{31} & \tau\_{32} & \tau\_{33} \end{bmatrix} \quad \text{and} \quad \vec{\nabla} \cdot \mathbf{T} = \begin{bmatrix} \frac{\partial \tau\_{11}}{\partial x\_{1}} + \frac{\partial \tau\_{21}}{\partial x\_{2}} + \frac{\partial \tau\_{31}}{\partial x\_{3}} \\ \frac{\partial \tau\_{12}}{\partial x\_{1}} + \frac{\partial \tau\_{22}}{\partial x\_{2}} + \frac{\partial \tau\_{32}}{\partial x\_{3}} \\ \frac{\partial \tau\_{13}}{\partial x\_{1}} + \frac{\partial \tau\_{22}}{\partial x\_{2}} + \frac{\partial \tau\_{33}}{\partial x\_{3}} \end{bmatrix} = \begin{bmatrix} \frac{\partial \tau\_{11}}{\partial x\_{1}} + \frac{\partial \tau\_{12}}{\partial x\_{2}} + \frac{\partial \tau\_{31}}{\partial x\_{3}} \\ \frac{\partial \tau\_{21}}{\partial x\_{1}} + \frac{\partial \tau\_{22}}{\partial x\_{2}} + \frac{\partial \tau\_{32}}{\partial x\_{3}} \\ \frac{\partial \tau\_{31}}{\partial x\_{1}} + \frac{\partial \tau\_{32}}{\partial x\_{2}} + \frac{\partial \tau\_{33}}{\partial x\_{3}} \end{bmatrix}$$

with

$$\tau\_{i\bar{j}} = \mu \left[ \frac{\partial u\_{\bar{i}}}{\partial \mathbf{x}\_{\bar{j}}} + \frac{\partial u\_{\bar{j}}}{\partial \mathbf{x}\_{\bar{i}}} \right] - \delta\_{i\bar{j}} \frac{2}{N} \left[ \mu \vec{\nabla} \cdot \vec{q} + \rho k \right] \tag{2}$$

where *μ* = *μ<sup>m</sup>* + *μ<sup>t</sup>* is the total dynamic viscosity, with *μ<sup>m</sup>* being the molecular viscosity and *μ<sup>t</sup>* 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 <sup>i</sup>* = −2*ρk*, where *u <sup>i</sup>* is the turbulent velocity in the *i th* direction.

Equation (1) can also be written as follows, as shown in the Appendix A:

$$
\rho \frac{\partial \vec{q}}{\partial t} + \rho (\vec{q} \cdot \vec{\nabla}) \vec{q}
\quad = \quad - \vec{\nabla} p + 2 \vec{\nabla} \cdot (\mu \vec{\nabla} \vec{q}) + \vec{\nabla} \times (\mu \vec{\omega})
$$

$$
$$

where *ω*- = (*ω*1, *<sup>ω</sup>*2, *<sup>ω</sup>*3) = ∇ × - *<sup>q</sup>* is the vorticity of the flow, and the tensor <sup>∇</sup>- *q* is the *gradient* of *q*, defined as follows

$$
\vec{\nabla}\vec{q} = \begin{bmatrix}
\frac{\partial u\_1}{\partial x\_1} & \frac{\partial u\_2}{\partial x\_1} & \frac{\partial u\_3}{\partial x\_1} \\
\frac{\partial u\_1}{\partial x\_2} & \frac{\partial u\_2}{\partial x\_2} & \frac{\partial u\_3}{\partial x\_2} \\
\frac{\partial u\_1}{\partial x\_3} & \frac{\partial u\_2}{\partial x\_3} & \frac{\partial u\_3}{\partial x\_3}
\end{bmatrix}
$$

*J. Mar. Sci. Eng.* **2020**, *8*, 191

With the help of some identities of vector calculus, as shown in the Appendix A, Equation (3) becomes:

$$
\rho \frac{\partial \vec{q}}{\partial t} + \rho (\vec{q} \cdot \nabla) \vec{q} \quad = \quad - \vec{\nabla} p + \mu \nabla^2 \vec{q} + 2(\vec{\nabla} \mu \cdot \vec{\nabla}) \vec{q} + \mu \vec{\nabla} (\vec{\nabla} \cdot \vec{q}) + \vec{\nabla} \mu \times \vec{\omega}
$$

$$
$$

Dividing Equation (4) by *ρ* and rearranging will give us:

$$\begin{split} \frac{\partial \vec{q}}{\partial t} + (\vec{q} \cdot \nabla) \vec{q} &= \quad -\frac{\not\nabla p}{\rho} + \nu \nabla^2 \vec{q} + \nu \nabla (\nabla \cdot \vec{q}) + \frac{2}{\rho} (\nabla \mu \cdot \nabla) \vec{q} + \frac{1}{\rho} \nabla \mu \times \vec{\omega} \\ \quad - \frac{2}{N\rho} \vec{\nabla} [\mu \vec{\nabla} \cdot \vec{q}] - \frac{2}{N\rho} \vec{\nabla} [\rho k] - \frac{\not\nabla \Phi}{\rho} \end{split} \tag{5}$$

where *ν* = *<sup>μ</sup> <sup>ρ</sup>* is the kinematic (total) viscosity

#### *2.1. 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 *<sup>D</sup><sup>ρ</sup> Dt* = 0), Equation (5) becomes:

$$\frac{\partial \vec{q}}{\partial t} + (\vec{q} \cdot \vec{\nabla})\vec{q} = -\frac{\vec{\nabla}p}{\rho} + \nu \nabla^2 \vec{q} + 2(\vec{\nabla}\nu \cdot \vec{\nabla})\vec{q} + \vec{\nabla}\nu \times \vec{\omega} - \frac{2}{N} \vec{\nabla}k - \frac{\vec{\nabla}\Phi}{\rho} \tag{6}$$

In the case of Newtonian fluid in laminar flow (*μ* = *μ<sup>m</sup>* =constant and *k* = 0) and incompressible flow, in which case ∇ · - *q* = 0, Equation (4) reduces to its most common form:

$$
\rho \frac{\partial \vec{q}}{\partial t} + \rho (\vec{q} \cdot \vec{\nabla}) \vec{q} = -\vec{\nabla}p + \mu\_m \nabla^2 \vec{q} - \vec{\nabla}\Phi \tag{7}
$$

#### **3. Vorticity Equation in 2D**

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

$$\begin{aligned} \frac{\partial\boldsymbol{\omega}}{\partial t} + \boldsymbol{\nabla} \cdot (\boldsymbol{\omega}\boldsymbol{\tilde{q}}) &= \quad -\boldsymbol{\nabla} \left(\frac{1}{\rho}\right) \times \boldsymbol{\nabla}\boldsymbol{p} + \boldsymbol{\nabla}\boldsymbol{\nu} \times \boldsymbol{\nabla}^{2}\boldsymbol{\tilde{q}} + \boldsymbol{\nu}\boldsymbol{\nabla}^{2}\boldsymbol{\omega} + \boldsymbol{\nabla}\boldsymbol{\nu} \times \boldsymbol{\nabla}(\boldsymbol{\nabla}\cdot\boldsymbol{\tilde{q}}) \\ &- \boldsymbol{\omega}\boldsymbol{\nabla}^{2}\boldsymbol{\nu} - \boldsymbol{\omega}\boldsymbol{\nabla} \cdot \left[\boldsymbol{\nu}\frac{\boldsymbol{\nabla}\boldsymbol{\rho}}{\rho}\right] + \boldsymbol{\vec{\nabla}}\boldsymbol{\nu} \cdot \boldsymbol{\nabla}\boldsymbol{\omega} + \frac{\boldsymbol{\nu}}{\rho} \boldsymbol{\vec{\nabla}}\boldsymbol{\rho} \cdot \boldsymbol{\nabla}\boldsymbol{\omega} \\ &+ 2\left[\frac{\partial}{\partial\boldsymbol{z}}\left(\frac{\boldsymbol{\nabla}\boldsymbol{\mu}}{\rho}\right) \cdot \boldsymbol{\nabla}\boldsymbol{u} - \frac{\partial}{\partial\boldsymbol{x}}\left(\frac{\boldsymbol{\nabla}\boldsymbol{\mu}}{\rho}\right) \cdot \boldsymbol{\nabla}\boldsymbol{w}\right] \\ &- \boldsymbol{\nabla}\left(\frac{1}{\rho}\right) \times \boldsymbol{\nabla}\left(\mu\boldsymbol{\nabla}\cdot\boldsymbol{\vec{q}}\right) - \boldsymbol{\nabla}\left(\frac{1}{\rho}\right) \times \boldsymbol{\nabla}\left(\rho\boldsymbol{k}\right) - \boldsymbol{\nabla}\left(\frac{1}{\rho}\right) \times \boldsymbol{\nabla}\boldsymbol{\Phi} \end{aligned} \tag{8}$$

As shown in the Appendix A Equation (8) may also be written as follows:

$$\begin{split} \frac{\partial \boldsymbol{\omega}}{\partial t} + \vec{\nabla} \cdot (\boldsymbol{\omega} \vec{\eta}) &= \quad - \vec{\nabla} \left( \frac{1}{\rho} \right) \times \vec{\nabla} p + \nabla^2 (\nu \boldsymbol{\omega}) - 2 \boldsymbol{\omega} \nabla^2 \boldsymbol{\nu} \\ &+ 2 \boldsymbol{\nabla} \boldsymbol{\nu} \times \boldsymbol{\nabla} (\boldsymbol{\nabla} \cdot \vec{\eta}) - \boldsymbol{\omega} \boldsymbol{\nabla} \cdot \left[ \boldsymbol{\nu} \frac{\boldsymbol{\nabla} \rho}{\rho} \right] + \frac{\nu}{\rho} \boldsymbol{\nabla} \rho \cdot \boldsymbol{\nabla} \boldsymbol{\omega} \\ &+ 2 \left[ \frac{\partial}{\partial \boldsymbol{z}} \left( \frac{\vec{\nabla} \mu}{\rho} \right) \cdot \vec{\nabla} \boldsymbol{u} - \frac{\partial}{\partial \boldsymbol{x}} \left( \frac{\vec{\nabla} \mu}{\rho} \right) \cdot \vec{\nabla} \boldsymbol{w} \right] \\ &- \boldsymbol{\nabla} \left( \frac{1}{\rho} \right) \times \boldsymbol{\nabla} (\mu \boldsymbol{\nabla} \cdot \vec{\eta}) - \boldsymbol{\nabla} \left( \frac{1}{\rho} \right) \times \boldsymbol{\nabla} (\rho \boldsymbol{k}) - \boldsymbol{\nabla} \left( \frac{1}{\rho} \right) \times \boldsymbol{\nabla} \boldsymbol{\Phi} \end{split} \tag{9}$$

#### *3.1. 2D Flow of Fluid with Constant Density*

In that case *<sup>ρ</sup>* <sup>=</sup> *constant*, <sup>∇</sup> *<sup>ρ</sup>* <sup>=</sup> 0, <sup>∇</sup>- <sup>1</sup> *ρ* <sup>=</sup> 0, and ∇ · - *q* = 0. Then Equation (9) reduces to:

$$\frac{\partial \omega}{\partial t} + \vec{\nabla} \cdot (\omega \vec{q}) = \nabla^2 (\nu \omega) - 2\omega \nabla^2 \nu - 2 \frac{\partial \vec{\nabla} \nu}{\partial x} \cdot \vec{\nabla} w + 2 \frac{\partial \vec{\nabla} \nu}{\partial z} \cdot \vec{\nabla} u \tag{10}$$

or, by expressing *ω* = *<sup>∂</sup><sup>u</sup> <sup>∂</sup><sup>z</sup>* <sup>−</sup> *<sup>∂</sup><sup>w</sup> <sup>∂</sup><sup>x</sup>* in the second term on the RHS of Equation (10), and by making use of the continuity equation, *<sup>∂</sup><sup>u</sup> <sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>∂</sup><sup>w</sup> <sup>∂</sup><sup>z</sup>* = 0, we get:

$$\frac{\partial \omega}{\partial t} + \vec{\nabla} \cdot (\omega \vec{q}) = \nabla^2 (\nu \omega) + 2 \frac{\partial^2 \nu}{\partial z^2} \frac{\partial w}{\partial x} - 2 \frac{\partial^2 \nu}{\partial x^2} \frac{\partial u}{\partial z} + 4 \frac{\partial^2 \nu}{\partial x \partial z} \frac{\partial u}{\partial x} \tag{11}$$

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 *ν* = *<sup>μ</sup><sup>m</sup> <sup>ρ</sup>* = *constant*, Equation (10) becomes the vorticity equation in its most common diffusion equation form:

$$\frac{\partial \omega}{\partial t} + \vec{\nabla} \cdot (\omega \vec{q}) = \nu \nabla^2 \omega \tag{12}$$

or, since ∇ · - *q* = 0:

$$\frac{D\omega}{Dt} \equiv \frac{\partial\omega}{\partial t} + \vec{q} \cdot \nabla\omega = \nu\nabla^2\omega\tag{13}$$

#### *3.2. 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 *<sup>∂</sup> <sup>∂</sup><sup>x</sup>* << *<sup>∂</sup> <sup>∂</sup><sup>z</sup>* , 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:

$$\frac{\partial \boldsymbol{\omega}}{\partial t} + \nabla \cdot (\boldsymbol{\omega} \vec{q}) = \nabla^2 (\boldsymbol{\nu} \boldsymbol{\omega}) + \nabla \cdot \left[ \nu \frac{\vec{\nabla} \rho}{\rho} \boldsymbol{\omega} \right] \tag{14}$$

#### **4. 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.

**Figure 2.** Vorticity contour plots from RANS (ANSYS/Fluent, **top**) and VISVE (**bottom**) for turbulent flow around hydrofoil. *Re* <sup>=</sup> <sup>2</sup> <sup>×</sup> <sup>10</sup>6.

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.

**Figure 3.** Velocity profiles of turbulent flow around hydrofoil, from RANS (ANSYS/Fluent) and from VISVE, at different locations along the chord, *<sup>c</sup>*, on the suction side. *Re* <sup>=</sup> <sup>2</sup> <sup>×</sup> 106.

**Figure 4.** Vorticity contour plots and velocity profiles for turbulent flow over cylinder from RANS (ANSYS/Fluent) and from VISVE, at t = 4 s and *Re* = 106.

**Acknowledgments:** Support for this research was provided by the U.S. Office of Naval Research (Grant Number N00014-18-1-2276; Ki-Han Kim) and by Phase VIII of the "Consortium on Cavitation Performance of High Speed Propulsors".

**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):

$$
\rho \frac{\partial u\_i}{\partial t} + \rho \vec{q} \cdot \nabla u\_i = -\frac{\partial p}{\partial \mathbf{x}\_i} + \frac{\partial \mathbf{r}\_{ji}}{\partial \mathbf{x}\_j} - \frac{\partial \Phi}{\partial \mathbf{x}\_i} = -\frac{\partial p}{\partial \mathbf{x}\_i} + \frac{\partial \mathbf{r}\_{ij}}{\partial \mathbf{x}\_j} - \frac{\partial \Phi}{\partial \mathbf{x}\_i} \tag{A1}
$$

with the *∂τji ∂xj* 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):

$$
\mu \left[ \frac{\partial u\_i}{\partial \mathbf{x}\_j} + \frac{\partial u\_j}{\partial \mathbf{x}\_i} \right] = 2\mu \frac{\partial u\_i}{\partial \mathbf{x}\_j} + \mu \left[ \frac{\partial u\_j}{\partial \mathbf{x}\_i} - \frac{\partial u\_i}{\partial \mathbf{x}\_j} \right] = 2\mu \frac{\partial u\_i}{\partial \mathbf{x}\_j} + \epsilon\_{kij}\mu\omega\_k \tag{A2}
$$

where *kij* is the Levi-Civita symbol.

Then for given *i* and with *j* = 1, 2, 3, using the Einstein notation we have:

$$\begin{split} \frac{\partial}{\partial \mathbf{x}\_{j}} \left[ \mu \left( \frac{\partial \mathbf{u}\_{i}}{\partial \mathbf{x}\_{j}} + \frac{\partial \mathbf{u}\_{j}}{\partial \mathbf{x}\_{i}} \right) \right] &= \ 2 \frac{\partial}{\partial \mathbf{x}\_{j}} \left( \mu \frac{\partial \mathbf{u}\_{i}}{\partial \mathbf{x}\_{j}} \right) + \epsilon\_{kij} \frac{\partial (\mu \boldsymbol{\omega}\_{k})}{\partial \mathbf{x}\_{j}} \\ &= \ 2 \vec{\nabla} \cdot (\mu \vec{\nabla} \boldsymbol{u}\_{i}) + \epsilon\_{ijk} \frac{\partial (\mu \boldsymbol{\omega}\_{k})}{\partial \mathbf{x}\_{j}} \\ &= \ 2 \vec{\nabla} \cdot (\mu \vec{\nabla} \boldsymbol{u}\_{i}) + \epsilon\_{ij'k'} \frac{\partial (\mu \boldsymbol{\omega}\_{k'})}{\partial \mathbf{x}\_{j'}} + \epsilon\_{ik'j'} \frac{\partial (\mu \boldsymbol{\omega}\_{j'})}{\partial \mathbf{x}\_{k'}} \quad (\text{i} \neq k' \neq j' \neq i) \\ &= \ 2 \vec{\nabla} \cdot (\mu \vec{\nabla} \boldsymbol{u}\_{i}) + \epsilon\_{ij'k'} \frac{\partial (\mu \boldsymbol{\omega}\_{k'})}{\partial \mathbf{x}\_{j'}} - \epsilon\_{ij'k'} \frac{\partial (\mu \boldsymbol{\omega}\_{j'})}{\partial \mathbf{x}\_{k'}} \quad (\text{i} \neq k' \neq j' \neq i) \\ &= \ 2 \vec{\nabla} \cdot (\mu \vec{\nabla} \boldsymbol{u}\_{i}) + [\vec{\nabla} \times (\mu \vec{\omega})]\_{i} \end{split} \tag{A3}$$

where [∇ × - (*μω*- )]*<sup>i</sup>* means the *i th* component of ∇ × - (*μω*-)

Based on Equations (2) and (A3) we get:

$$\frac{\partial \mathsf{T}\_{lj}}{\partial \mathbf{x}\_{j}} = 2\vec{\nabla} \cdot (\mu \vec{\nabla} u\_{l}) + [\vec{\nabla} \times (\mu \vec{\omega})]\_{i} - \frac{2}{N} \frac{\partial}{\partial \mathbf{x}\_{i}} \left[\mu \vec{\nabla} \cdot \vec{q} + \rho k\right] \tag{A4}$$

which eventually leads to Equation (3).

*Appendix A.2. How to Get from Equation (3) to Equation (4)*

$$\vec{\nabla} \cdot (\mu \vec{\nabla} \vec{q}) = \begin{bmatrix} \vec{\nabla} \cdot (\mu \vec{\nabla} u\_1) \\ \vec{\nabla} \cdot (\mu \vec{\nabla} u\_2) \\ \vec{\nabla} \cdot (\mu \vec{\nabla} u\_3) \end{bmatrix} \tag{A5}$$

Using vector identities we get:

$$\nabla \cdot (\mu \nabla u\_i) = \nabla \mu \cdot \nabla u\_i + \mu \nabla^2 u\_i \tag{A6}$$

or

$$\vec{\nabla} \cdot (\mu \vec{\nabla} \vec{q}) = (\vec{\nabla} \mu \cdot \vec{\nabla}) \vec{q} + \mu \nabla^2 \vec{q} \tag{A7}$$

Then

$$\vec{\nabla} \times (\mu \vec{\omega}) = \mu \vec{\nabla} \times \vec{\omega} + \vec{\nabla} \mu \times \vec{\omega} \tag{A8}$$

$$
\nabla \times \vec{\omega} = \nabla \times (\nabla \times \vec{q}) = \nabla(\nabla \cdot \vec{q}) - \nabla^2 \vec{q} \tag{A9}
$$

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*<sup>1</sup> = *x*, *x*<sup>2</sup> = *y*, *x*<sup>3</sup> = *z*; *u*<sup>1</sup> = *u*, *u*<sup>2</sup> = *v*, *u*<sup>3</sup> = *w*; *<sup>∂</sup> <sup>∂</sup><sup>y</sup>* = 0 and *v* = 0. The velocity has two components: *q* = (*u*, *w*), and the vorticity has only one component in the *y* direction *ω*<sup>2</sup> = *ω* = *<sup>∂</sup><sup>u</sup> <sup>∂</sup><sup>z</sup>* <sup>−</sup> *<sup>∂</sup><sup>w</sup> ∂x* .

First:

$$
\vec{\nabla} \times \left(\frac{\partial \vec{q}}{\partial t}\right) = \frac{\partial (\vec{\nabla} \times \vec{q})}{\partial t} = \frac{\partial \vec{\omega}}{\partial t} \tag{A10}
$$

Then by using the identity:

$$(\vec{q}\cdot\vec{\nabla})\vec{q} = (\vec{\nabla}\times\vec{q})\times\vec{q} + \vec{\nabla}\left(\frac{q^2}{2}\right) = \vec{\omega}\times\vec{q} + \vec{\nabla}\left(\frac{q^2}{2}\right) \tag{A11}$$

we get

$$\vec{\nabla} \times \left[ (\vec{q} \cdot \vec{\nabla}) \vec{q} \right] = \vec{\nabla} \times (\vec{\omega} \times \vec{q}) \tag{A12}$$

due to the identity:∇ × - <sup>∇</sup>*f* = 0

Then using the following identities:

$$\vec{\nabla} \times (\vec{\omega} \times \vec{q}) = (\vec{q} \cdot \vec{\nabla})\vec{\omega} - (\vec{\omega} \cdot \vec{\nabla})\vec{q} + \vec{\omega}(\vec{\nabla} \cdot \vec{q}) - \vec{q}(\vec{\nabla} \cdot \vec{\omega}) \tag{A13}$$

$$
\vec{\nabla} \cdot \vec{\omega} = \vec{\nabla} \cdot (\vec{\nabla} \times \vec{q}) = 0 \tag{A14}
$$

We get:

$$\vec{\nabla} \times \left[ \text{LHS of equation (5)} \right] = \frac{\partial \vec{\omega}}{\partial t} + (\vec{q} \cdot \vec{\nabla})\vec{\omega} - (\vec{\omega} \cdot \vec{\nabla})\vec{q} + \vec{\omega}(\vec{\nabla} \cdot \vec{q}) \tag{A15}$$

Now, in 2-D, since the vortex stretching term (*ω*- · <sup>∇</sup>- )*q* = 0, the above equation becomes:

$$\vec{\nabla} \times \left[ \text{LHS of equation (5)} \right]\_{2D} = \frac{\partial \omega}{\partial t} + (\vec{q} \cdot \vec{\nabla})\omega + \omega(\vec{\nabla} \cdot \vec{q}) = \frac{\partial \omega}{\partial t} + \vec{\nabla} \cdot (\omega \vec{q}) \tag{A16}$$

We now take the ∇×of each term on the RHS of Equation (5):

• The ∇×of the 1st term on the RHS of (5):

$$\nabla \times \left(\frac{\nabla p}{\rho}\right) = \nabla \left(\frac{1}{\rho}\right) \times \nabla p + \left(\frac{1}{\rho}\right) \nabla \times \nabla p = \nabla \left(\frac{1}{\rho}\right) \times \nabla p \tag{A17}$$

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):

$$\nabla \times (\nu \nabla^2 \vec{q}) = \nabla \nu \times \nabla^2 \vec{q} + \nu \nabla^2 (\nabla \times \vec{q}) = \nabla \nu \times \nabla^2 \vec{q} + \nu \nabla^2 \omega \tag{A18}$$

• The ∇×of the 3rd term on the RHS of (5):

$$\vec{\nabla} \times \left[ \nu \vec{\nabla} (\vec{\nabla} \cdot \vec{q}) \right] = \vec{\nabla} \nu \times \vec{\nabla} (\vec{\nabla} \cdot \vec{q}) + \nu \vec{\nabla} \times \left[ \vec{\nabla} (\vec{\nabla} \cdot \vec{q}) \right] = \vec{\nabla} \nu \times \vec{\nabla} (\vec{\nabla} \cdot \vec{q}) \tag{A19}$$

• The ∇×of the 4th term (divided by 2) on the RHS of (5):

We only handle this term in 2D. In that case:

$$\left\{\vec{\nabla} \times \left[\left(\frac{\vec{\nabla}\mu}{\rho} \cdot \vec{\nabla}\right)\vec{q}\right]\right\}\_{2D} = \frac{\partial}{\partial z}\left(\frac{\vec{\nabla}\mu}{\rho} \cdot \vec{\nabla}u\right) - \frac{\partial}{\partial x}\left(\frac{\vec{\nabla}\mu}{\rho} \cdot \vec{\nabla}w\right) = \tag{A20}$$

$$\mathcal{L} = \frac{\partial}{\partial \mathbf{z}} \left( \frac{\nabla \mu}{\rho} \right) \cdot \nabla u + \frac{\nabla \mu}{\rho} \cdot \frac{\partial \nabla u}{\partial z} - \frac{\partial}{\partial \mathbf{x}} \left( \frac{\nabla \mu}{\rho} \right) \cdot \nabla w - \frac{\nabla \mu}{\rho} \cdot \frac{\partial \nabla w}{\partial \mathbf{x}} = \tag{A21}$$

$$\mathcal{I} = \frac{\partial}{\partial z} \left( \frac{\vec{\nabla} \mu}{\rho} \right) \cdot \vec{\nabla} u - \frac{\partial}{\partial x} \left( \frac{\vec{\nabla} \mu}{\rho} \right) \cdot \vec{\nabla} w + \frac{\vec{\nabla} \mu}{\rho} \cdot \vec{\nabla} \left( \frac{\partial u}{\partial z} - \frac{\partial w}{\partial x} \right) = \tag{A22}$$

$$\mathcal{L} = \frac{\partial}{\partial z} \left( \frac{\nabla \mu}{\rho} \right) \cdot \nabla \mu - \frac{\partial}{\partial x} \left( \frac{\nabla \mu}{\rho} \right) \cdot \nabla w + \frac{\nabla \mu}{\rho} \cdot \nabla \omega \tag{A23}$$

• The ∇×of the 5th term on the RHS of (5):

$$\vec{\nabla} \times \left(\frac{1}{\rho} \vec{\nabla} \mu \times \vec{\omega}\right) = \vec{\nabla} \left(\frac{1}{\rho}\right) \times \left(\vec{\nabla} \mu \times \vec{\omega}\right) + \frac{1}{\rho} \vec{\nabla} \times \left(\vec{\nabla} \mu \times \vec{\omega}\right) \tag{A24}$$

then

$$\vec{\nabla} \times (\vec{\nabla}\mu \times \vec{\omega}) = (\vec{\omega} \cdot \vec{\nabla})\vec{\nabla}\mu - (\vec{\nabla}\mu \cdot \vec{\nabla})\vec{\omega} + \vec{\nabla}\mu(\vec{\nabla} \cdot \vec{\omega}) - \vec{\omega}(\vec{\nabla} \cdot \vec{\nabla}\mu) \tag{A25}$$

with ∇ · *ω*- = 0, and since *ω*- · <sup>∇</sup>-= 0 in 2D we get:

$$[\vec{\nabla} \times (\vec{\nabla} \mu \times \vec{\omega})]\_{2D} = -\vec{\nabla} \mu \cdot \vec{\nabla} \omega - \omega \nabla^2 \mu \tag{A26}$$

In addition, in 2D:

$$(\nabla \mu \times \vec{\omega})\_{2D} = \left(\frac{\partial \mu}{\partial \mathbf{x}}, 0, \frac{\partial \mu}{\partial z}\right) \times (0, \omega, 0) = \left(-\omega \frac{\partial \mu}{\partial z}, 0, \omega \frac{\partial \mu}{\partial \mathbf{x}}\right) \tag{A27}$$

and

$$\left[\vec{\nabla}\left(\frac{1}{\rho}\right)\times(\vec{\nabla}\mu\times\vec{\omega})\right]\_{2D}=\left(\frac{\partial(1/\rho)}{\partial x},0,\frac{\partial(1/\rho)}{\partial z}\right)\times\left(-\omega\frac{\partial\mu}{\partial z},0,\omega\frac{\partial\mu}{\partial x}\right)=\tag{A28}$$

$$=-\omega \left[ \frac{\partial (1/\rho)}{\partial x} \frac{\partial \mu}{\partial x} + \frac{\partial (1/\rho)}{\partial z} \frac{\partial \mu}{\partial z} \right] = -\omega \vec{\nabla} \left( \frac{1}{\rho} \right) \cdot \vec{\nabla} \mu \tag{A29}$$

Finally

$$\left[\nabla \times \left(\frac{1}{\rho} \nabla \mu \times \vec{\omega}\right)\right]\_{2D} = -\omega \nabla \left(\frac{1}{\rho}\right) \cdot \nabla \mu - \frac{\vec{\nabla}\mu \cdot \vec{\nabla}\omega}{\rho} - \frac{\omega \nabla^2 \mu}{\rho} = \tag{A30}$$

$$=-\omega \nabla \left(\frac{\vec{\nabla}\mu}{\rho}\right) - \frac{\vec{\nabla}\mu \cdot \vec{\nabla}\omega}{\rho} \tag{A31}$$

The combination of the ∇×of the 4th and the 5th terms will give us:

$$2\frac{\partial}{\partial z}\left(\frac{\vec{\nabla}\mu}{\rho}\right)\cdot\vec{\nabla}\mu - 2\frac{\partial}{\partial x}\left(\frac{\vec{\nabla}\mu}{\rho}\right)\cdot\vec{\nabla}w + \frac{\vec{\nabla}\mu\cdot\vec{\nabla}\omega}{\rho} - \omega\vec{\nabla}\left(\frac{\vec{\nabla}\mu}{\rho}\right) \tag{A32}$$

*J. Mar. Sci. Eng.* **2020**, *8*, 191

The last two terms above can also be rewritten by using *<sup>μ</sup>* <sup>=</sup> *νρ* and <sup>∇</sup>- (*νρ*) = *<sup>ν</sup>*∇ *<sup>ρ</sup>* <sup>+</sup> *<sup>ρ</sup>*∇ *ν*, which renders the RHS of Equation (8)

*Appendix A.4. How to Get from Equation (8) to Equation (9)*

In 2D *v* = 0 and *<sup>∂</sup> <sup>∂</sup><sup>y</sup>* <sup>=</sup> 0, and *<sup>ω</sup>* <sup>=</sup> *<sup>∂</sup><sup>u</sup> <sup>∂</sup><sup>z</sup>* <sup>−</sup> *<sup>∂</sup><sup>w</sup> ∂x* [∇ *<sup>ν</sup>* · <sup>∇</sup> *<sup>ω</sup>*]2*<sup>D</sup>* <sup>=</sup> *∂ν ∂x ∂ω <sup>∂</sup><sup>x</sup>* <sup>+</sup> *∂ν ∂z ∂ω <sup>∂</sup><sup>z</sup>* <sup>=</sup> *∂ν ∂x ∂ ∂x ∂u <sup>∂</sup><sup>z</sup>* <sup>−</sup> *<sup>∂</sup><sup>w</sup> ∂x* + *∂ν ∂z ∂ ∂z ∂u <sup>∂</sup><sup>z</sup>* <sup>−</sup> *<sup>∂</sup><sup>w</sup> ∂x* = *∂ν ∂x ∂*2*u <sup>∂</sup>x∂<sup>z</sup>* <sup>−</sup> *∂ν ∂x ∂*2*w <sup>∂</sup>x*<sup>2</sup> <sup>+</sup> *∂ν ∂z ∂*2*u <sup>∂</sup>z*<sup>2</sup> <sup>−</sup> *∂ν ∂z ∂*2*w <sup>∂</sup>x∂<sup>z</sup>* <sup>=</sup> *∂ν ∂x ∂ ∂z* −*∂<sup>w</sup> <sup>∂</sup><sup>z</sup>* <sup>+</sup> ∇ · - *q* <sup>−</sup> *∂ν ∂x ∂*2*w <sup>∂</sup>x*<sup>2</sup> <sup>+</sup> *∂ν ∂z ∂*2*u <sup>∂</sup>z*<sup>2</sup> <sup>−</sup> *∂ν ∂z ∂ ∂z* −*∂<sup>u</sup> <sup>∂</sup><sup>x</sup>* <sup>+</sup> ∇ · - *q* = −*∂ν ∂x* <sup>∇</sup>2*<sup>w</sup>* <sup>+</sup> *∂ν ∂z* <sup>∇</sup>2*<sup>u</sup>* <sup>+</sup> *∂ν ∂x <sup>∂</sup>*(∇ · - *q*) *<sup>∂</sup><sup>z</sup>* <sup>−</sup> *∂ν ∂z <sup>∂</sup>*(∇ · - *q*) *<sup>∂</sup><sup>x</sup>* <sup>=</sup> + ∇ *<sup>ν</sup>* × ∇<sup>2</sup> *q* , <sup>2</sup>*<sup>D</sup>* <sup>−</sup> [∇ *<sup>ν</sup>* <sup>×</sup> <sup>∇</sup>- (∇ · - *q*)]2*<sup>D</sup>* (A33)

or

$$\left[\vec{\nabla}\boldsymbol{\nu}\times\boldsymbol{\nabla}^{2}\vec{q}\right]\_{2D}=[\vec{\nabla}\boldsymbol{\nu}\cdot\vec{\nabla}\boldsymbol{\omega}]\_{2D}+[\vec{\nabla}\boldsymbol{\nu}\times\vec{\nabla}(\vec{\nabla}\cdot\vec{q})]\_{2D}\tag{A34}$$

where we have made use of *<sup>∂</sup><sup>u</sup> <sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>∂</sup><sup>w</sup> <sup>∂</sup><sup>z</sup>* <sup>=</sup> ∇ · - *q*

We finally get to Equation (9) by considering the identity:

$$
\nabla^2(\nu\omega) = \nu\nabla^2\omega + \omega\nabla^2\nu + 2\nabla\nu \cdot \nabla\omega \tag{A35}
$$

#### *Appendix A.5. How to Get from Equation (9) to Equation (14)*

In that case we assume that the hydrofoil is placed along axis *x*, with the inflow also along *x* axis, and that *<sup>∂</sup> <sup>∂</sup><sup>x</sup>* << *<sup>∂</sup> <sup>∂</sup><sup>z</sup>* , especially within the narrow region close to the hydrofoil and its wake where *ω* = 0.

Then, since we ignore *<sup>∂</sup> <sup>∂</sup><sup>x</sup>* then all <sup>∇</sup>- *<sup>F</sup>* point in the *<sup>z</sup>* direction and thus <sup>∇</sup>- *<sup>F</sup>* <sup>×</sup> <sup>∇</sup>- *G* ≈ 0. In addition the following approximations can be made:

$$-2\omega \nabla^2 \nu \quad \approx \quad -2\omega \frac{\partial^2 \nu}{\partial z^2} \tag{A36}$$

$$-\omega \vec{\nabla} \left( \nu \frac{\partial \vec{\nabla} \rho}{\rho} \right) \quad \approx \quad -\omega \frac{\partial}{\partial z} \left( \frac{\nu}{\rho} \frac{\partial \rho}{\partial z} \right) \tag{A37}$$

$$\frac{\nu}{\rho} \nabla \rho \cdot \nabla \omega \quad \approx \quad \frac{\nu}{\rho} \frac{\partial \rho}{\partial z} \frac{\partial \omega}{\partial z} \tag{A38}$$

$$2\frac{\partial}{\partial z}\left(\frac{\vec{\nabla}\mu}{\rho}\right)\cdot\vec{\nabla}u \quad \approx \quad 2\frac{\partial}{\partial z}\left(\frac{1}{\rho}\frac{\partial\mu}{\partial z}\right)\frac{\partial u}{\partial z} \approx 2\frac{\partial}{\partial z}\left(\frac{1}{\rho}\frac{\partial\mu}{\partial z}\right)\omega = \tag{A39}$$

$$2\frac{\partial}{\partial z}\left(\frac{\partial \nu}{\partial z} + \frac{\nu}{\rho}\frac{\partial \rho}{\partial z}\right)\omega = 2\omega \frac{\partial^2 \nu}{\partial z^2} + 2\omega \frac{\partial}{\partial z}\left(\frac{\nu}{\rho}\frac{\partial \rho}{\partial z}\right) \tag{A40}$$

By adding the approximate expressions for all the terms shown above we get

$$
\frac{\nu}{\rho} \frac{\partial \rho}{\partial z} \frac{\partial \omega}{\partial z} + \omega \frac{\partial}{\partial z} \left( \frac{\nu}{\rho} \frac{\partial \rho}{\partial z} \right) = \frac{\partial}{\partial z} \left( \omega \frac{\nu}{\rho} \frac{\partial \rho}{\partial z} \right) \tag{A41}
$$

Finally we get to Equation (14) by making the following approximation:

$$\frac{\partial}{\partial z}\left(\omega \frac{\nu}{\rho} \frac{\partial \rho}{\partial z}\right) \approx \vec{\nabla} \cdot \left[\nu \frac{\vec{\nabla} \rho}{\rho} \omega\right] \tag{A42}$$

#### **References**


c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Journal of Marine Science and Engineering* Editorial Office E-mail: jmse@mdpi.com www.mdpi.com/journal/jmse

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18