34.2.2 Pressure Recovery

34.2.2.1 The Problem With Incompressible Fluids

At this point, we have all the equations required in a mathematical sense. We have a total of four equations (see Eq. 34.1, Eq. 34.2, Eq. 34.3, and Eq. 34.4) for a total of four dependent variables, v_x, v_y, v_z, and p. However, at this point, we run into a problem that we have worked around nicely so far. Until now, we considered the pressure to be applied from the outside, i.e., to be an imposed value which was left constant. In fact, in Poiseuille flow (see section 15.4) the pressure is the driving force of the flow. It occurs as a source term on the right-hand side of the Navier–Stokes equation (see Eq. 11.40) and was constant. If the pressure is not applied from the outside, e.g., because the flow is not pressure driven but volume driven by a volume-controlled pump such as a syringe or piston pump, there is no means for us to determine the pressure term. In addition to that we must assume the pressure to change across the domain. We know from the incompressible Bernoulli equation (see Eq. 14.11) that the pressure changes with changing velocity. However, we do not know the velocity distribution because this is the variable we are actually seeking.

Obviously, in a mathematical sense, it should be possible to determine four unknowns from four equations. However, if you look at Eq. 34.1, Eq. 34.2, and Eq. 34.3 you see that we actually only have the gradient, i.e., the derivative of the pressure in the equation. We cannot simply solve the equations for the value of the pressure at a given point as we can for the velocities. These occur both as a variable and derivatives on the left-hand sides of Eq. 34.1, Eq. 34.2, and Eq. 34.3. In addition to this Eq. 34.4 does not even have a distinct variable to solve for. It would be really convenient to rewrite Eq. 34.4 such that it depends on p. Then all four equations would have unknowns for which we could solve them.

Interestingly, we would have a significantly easier time if the fluid in question was compressible. In this case, the density would not be constant and we could bring in an additional equation, e.g., the ideal gas equation (see Eq. 6.8) we discussed in section 6.4 and section 12.8. Such an equation would link the pressure to the density and we would have equations which we could solve for them.

For incompressible fluids, we have to find a way around this problem. All of these approaches will try to rewrite the continuity equation, Eq. 34.4, such that it becomes a function of the pressure p. These approaches are therefore generally referred to as pressure term recovery approaches because this is what they do.

34.2.2.2 Artificial Compressibility Method

The most simple means of reintroducing the pressure into the continuity equation is by simply assuming that the fluid is pseudo-compressible. We introduce an artificial compressibility factor β given by

$\beta =\frac{d\tilde{\rho }}{dp}\to d\tilde{\rho }=\beta dp$

which simply states that a change in the pressure will lead to a change in density. This effect is assumed to be linear and β is assumed to be small, i.e., large changes in the pressure only lead to small changes in the density. Pressure changes will result in (very small) changes of the density. We therefore obtain a weak coupling between p and $\tilde{\rho }$ . As you can see, the density has been written as $\tilde{\rho }$ and is often referred to as the pseudo-density. This is an artificial value which has no further meaning. This value is not required for the momentum equations; we only require it for a small trick in the continuity equation.

We now add a time-dependent term on the left-hand side of the continuity equation, Eq. 34.4, which we also find in the incompressible continuity equation (see Eq. 10.16). We find

$\frac{\partial \tilde{\rho }}{\partial t}+\frac{\partial v_x}{\partial x}+\frac{\partial v_y}{\partial y}+\frac{\partial v_z}{\partial z}=0$

This equation simply states that if we allow a certain compressibility, the pseudo-density should change as a consequence of mass in- and outflux into and out of our control volume. We can now insert Eq. 34.5 into Eq. 34.6 to rewrite the pseudo-density as a function of the pressure to find

$\beta \frac{\partial p}{\partial t}+\frac{\partial v_x}{\partial x}+\frac{\partial v_y}{\partial y}+\frac{\partial v_z}{\partial z}=0$

where we now have a function of the pressure. We can rewrite Eq. 34.7 to solve for the pressure, which we can then insert into the momentum equations thereby recovering the pressure.

Choosing the Value β. There are certain restrictions for the value β. If the value is chosen too small, the equation becomes numerically stiff, i.e., there will be only very small contributions from the artificial compressibility and we only obtain very small increments in our iteration scheme. Therefore it takes a high number of iterations to come up with a value of p, which actually solves the momentum equation and the modified continuity equation (Eq. 34.7). For higher values of β we obtain fast convergence but increase the risk of obtaining physically meaningless results. This is due to the fact that at a certain point changes in the pressure (which occur due to the artificial compressibility) propagate faster through the domain than the momentum. The speed at which pressure changes propagate cannot be faster than the speed of sound. In general, β should be chosen such that the Mach number (see section 9.9.1) M < 1 [1]. It is therefore important to choose β based on the medium. In general, choosing smaller values may lead to a longer calculation time, but it ensures the results are physically meaningful.

34.2.2.3 SIMPLE Algorithm

The SIMPLE algorithm was originally described by Patankar and Spalding [2]. The acronym stands for “semiimplicit method for pressure-linked equations.” As the name suggests, this algorithm gives a semi-implicit correlation for the pressure dependency of the momentum equations and the continuity equation. There have been numerous improvements to the algorithm, all of which make use of one concept: Linking up changes in the pressure as the source terms of the momentum equation with changes of the velocity. If we solve Eq. 34.1, Eq. 34.2, and Eq. 34.3 for the values ${\overline{v}}_x$ , ${\overline{v}}_y$ , and ${\overline{v}}_z$ , respectively, we end up with

${\overline{v}}_x=\frac{1}{\frac{\partial v_x}{\partial x}}\left(\frac{1}{\rho }\left(-\frac{\partial p}{\partial x}+\eta \left(\frac{{\partial }^2{v}_x}{\partial x^2}+\frac{{\partial }^2{v}_x}{\partial y^2}+\frac{{\partial }^2{v}_x}{\partial z^2}\right)+k_x\right)-\left(\frac{\partial v_x}{\partial t}+v_y\frac{\partial v_x}{\partial y}+v_z\frac{\partial v_x}{\partial x}\right)\right)$

${\overline{v}}_y=\frac{1}{\frac{\partial v_y}{\partial x}}\left(\frac{1}{\rho }\left(-\frac{\partial p}{\partial y}+\eta \left(\frac{{\partial }^2{v}_y}{\partial x^2}+\frac{{\partial }^2{v}_y}{\partial y^2}+\frac{{\partial }^2{v}_y}{\partial z^2}\right)+k_y\right)-\left(\frac{\partial v_y}{\partial t}+v_x\frac{\partial v_y}{\partial x}+v_y\frac{\partial v_y}{\partial y}\right)\right)$

${\overline{v}}_z=\frac{1}{\frac{\partial v_z}{\partial x}}\left(\frac{1}{\rho }\left(-\frac{\partial p}{\partial z}+\eta \left(\frac{{\partial }^2{v}_z}{\partial x^2}+\frac{{\partial }^2{v}_z}{\partial y^2}+\frac{{\partial }^2{v}_z}{\partial z^2}\right)+k_z\right)-\left(\frac{\partial v_z}{\partial t}+v_x\frac{\partial v_z}{\partial x}+v_y\frac{\partial v_z}{\partial y}\right)\right)$

We will come to the overline notation in a moment.

Solving the Equations One By One. Mathematically speaking, Eq. 34.8, Eq. 34.9, and Eq. 34.10 are obviously not very useful forms of the momentum equations because we solve for the (unknown) values ${\overline{v}}_x$ , ${\overline{v}}_y$ , and ${\overline{v}}_z$ , using equations which contain gradients and derivatives of these unknown values on the right-hand side. However, in numerical solvers, we usually have to iterate the solution anyway. It therefore seems plausible to use the values from the previous iteration i on the right-hand side of the equations to calculate the values of the iteration i + 1, which are the values on the left-hand side of the equation. In order to make the scheme numerically even less expensive, we do not even solve these equations as a linear system of equations. However, mathematically speaking, this is what we should do as, e.g.,Eq. 34.8 requires the value ${\overline{v}}_y$ and ${\overline{v}}_z$ whereas Eq. 34.9 requires the values ${\overline{v}}_x$ and ${\overline{v}}_z$ , etc. However, as we iterate the solution it seems perfectly plausible to work with “old” values instead of solving linear systems of equations, which is definitely slower. This is why, we usually solve Eq. 34.8, Eq. 34.9, and Eq. 34.10 in a one-by-one fashion.

Satisfying the Momentum Equations:${\overline{\mathit{v}}}_{\mathit{x}}$ ,${\overline{\mathit{v}}}_{\mathit{y}}$ , and${\overline{\mathit{v}}}_{\mathit{z}}$ . Let us assume for a moment that we solved Eq. 34.8, Eq. 34.9, and Eq. 34.10 such that we obtain velocity values ${\overline{v}}_x$ , ${\overline{v}}_y$ , and ${\overline{v}}_z$ ,. These velocity values fulfill the momentum equations, i.e., satisfy the Navier–Stokes equation. This is why we label them with the overline. Obviously, as we have just learned, we will solve these equations one by one, which will yield values that are slightly off, but we ignore this for a moment.

Satisfying the Continuity Equation: ṽ_x, ṽ_y, and ṽ_z. In the next step, we need to solve the continuity equation (see Eq. 34.4) using the values ${\overline{v}}_x$ , ${\overline{v}}_y$ , and ${\overline{v}}_z$ . We will most likely find that these values do not satisfy the continuity equation. This is why we need to correct them, adding small terms ṽ_x, ṽ_y, and ṽ_z, which changes our solution to

$\begin{array}{l}{v}_x={\overline{v}}_x+{\tilde{v}}_x\\ v_y={\overline{v}}_y+{\tilde{v}}_y\\ v_z={\overline{v}}_z+{\tilde{v}}_z\end{array}$

where the values ${\overline{v}}_{\left\{x,y,z\right\}}$ satisfy the momentum equations and the terms ṽ_{x,y,z} are added to make sure the continuity equation is fulfilled. Having added these small correction terms ṽ_{x,y,z} will most likely lead to the problem that the velocity values v_x, v_y, and v_z will (most likely) not fulfill the momentum equations anymore. We can correct this by adding (again) small correction terms. As you can see, this becomes an iterative problem. We have to continue the process until we reach a certain problem-dependent precision upon which we assume the values to be sufficiently refined to satisfy both momentum and continuity equations.

Pressure Correction. Now the obvious question is the following: Where does the “pressure-linked” concept come into play? In fact, we cannot simply add velocity correction terms without a source term in the momentum equation. The additional velocity must originate from some physical source. We define this source to be the pressure, thereby achieving a semi-implicit link between the velocity and the pressure. If you consider Eq. 34.8, Eq. 34.9, and Eq. 34.10 you can see that we can rewrite the function of velocities ${\overline{v}}_x$ , ${\overline{v}}_y$ , and ${\overline{v}}_z$ sas a function of the pressure as

${\overline{v}}_x=c_x+c_{x\mathit{,}p}\frac{\partial \overline{p}}{\partial x}$

${\overline{v}}_y=c_y+c_{y\mathit{,}p}\frac{\partial \overline{p}}{\partial y}$

${\overline{v}}_z=c_z+c_{z\mathit{,}p}\frac{\partial \overline{p}}{\partial z}$

where we write the pressure as $\overline{p}$ , indicating that this is the pressure that satisfies the momentum equations. In analogy to the velocities, we now add a small pressure $\tilde{p}$ , which gives rise to the velocities ṽ_x, ṽ_y, and ṽ_z, which we need to satisfy the continuity equation. We therefore find

$\begin{array}{l}{\overline{v}}_x=c_x+c_{x\mathit{,}p}\frac{\partial \overline{p}+\tilde{p}}{\partial x}\\ =c_x+c_{x\mathit{,}p}\left(\frac{\partial \overline{p}}{\partial x}+\frac{\partial \tilde{p}}{\partial x}\right)\\ ={\overline{v}}_x+{\tilde{v}}_x\end{array}$

$\begin{array}{l}{\overline{v}}_y=c_y+c_{y\mathit{,}p}\frac{\partial \overline{p}+\tilde{p}}{\partial y}\\ =c_y+c_{y\mathit{,}p}\left(\frac{\partial \overline{p}}{\partial y}+\frac{\partial \tilde{p}}{\partial y}\right)\\ ={\overline{v}}_y+{\tilde{v}}_y\end{array}$

$\begin{array}{l}{\overline{v}}_z=c_z+c_{z\mathit{,}p}\frac{\partial \overline{p}+\tilde{p}}{\partial z}\\ =c_z+c_{z\mathit{,}p}\left(\frac{\partial \overline{p}}{\partial z}+\frac{\partial \tilde{p}}{\partial z}\right)\\ ={\overline{v}}_z+{\tilde{v}}_z\end{array}$

from which we learn that

${\tilde{v}}_x=c_{x\mathit{,}p}\frac{\partial \tilde{p}}{\partial x}$

${\tilde{v}}_y=c_{y\mathit{,}p}\frac{\partial \tilde{p}}{\partial y}$

${\tilde{v}}_z=c_{z\mathit{,}p}\frac{\partial \tilde{p}}{\partial z}$

This means that the velocity corrections are only dependent on the pressure corrections we make. This gives the semi-implicit link between changes in the pressure and the change in the local velocity fields.

Modified Continuity Equations. We can now rewrite the continuity equation (Eq. 34.4) as

$\begin{array}{l}\frac{\partial v_x}{\partial x}+\frac{\partial v_y}{\partial y}+\frac{\partial v_z}{\partial z}=0\\ \frac{\partial \left({\overline{v}}_x+{\tilde{v}}_x\right)}{\partial x}+\frac{\partial \left({\overline{v}}_y+{\tilde{v}}_y\right)}{\partial y}+\frac{\partial \left({\overline{v}}_z+{\tilde{v}}_z\right)}{\partial z}=0\end{array}$

where we use Eq. 34.11, Eq. 34.12, and Eq. 34.13 for rewriting ${\overline{v}}_x$ , ${\overline{v}}_y$ , and ${\overline{v}}_z$ and Eq. 34.17, Eq. 34.18, and Eq. 34.19 for rewriting ṽ_x, ṽ_y, and ṽ_z, in which case we find

$\frac{\partial {\overline{v}}_x}{\partial x}+c_{x\mathit{,}p}\frac{{\partial }^2\tilde{p}}{\partial x^2}+\frac{\partial {\overline{v}}_y}{\partial y}+c_{y\mathit{,}p}\frac{{\partial }^2\tilde{p}}{\partial y^2}+\frac{\partial {\overline{v}}_z}{\partial z}+c_{z\mathit{,}p}\frac{{\partial }^2\tilde{p}}{\partial z^2}=0$

Eq. 34.20 is an equation that contains the second derivative of the pressure. As we know from the most common numerical scheme for the second-order partial derivative (see Eq. 4.10) we will be able to calculate the pressure at the current cell or node in our domain using the second derivative.

Scheme. In general, the SIMPLE algorithm proceeds in these steps:

1. Solve ${\overline{v}}_x$ , ${\overline{v}}_y$ , and ${\overline{v}}_z$ using Eq. 34.14, Eq. 34.15, and Eq. 34.16.

2. Calculate the new pressure using Eq. 34.20.

3. Update the velocities ${\overline{v}}_x$ , ${\overline{v}}_y$ , and ${\overline{v}}_z$ with the new pressure term.

4. Check for convergence, i.e., check if the change in velocities of the current iteration is below a given threshold.

The SIMPLE algorithm is known for bad convergence and the pressure value usually has to be under-relaxed in order to obtain a stable scheme. This means that the updated value p^(i + 1) must be calculated as

$p^{\left(i+1\right)}=p^{\left(i\right)}+\omega \tilde{p}$

with 0 ≤ u < 1.

Improvements to the Algorithm. As you can see, we do not even calculate the corrected velocity terms ṽ_x, ṽ_y, and ṽ_z in the SIMPLE algorithm, which is the main reason why we need to apply under-relaxation for the pressure correction. If we were able to calculate these values, we could update them alongside the pressure, thereby balancing the momentum equations and obtaining faster convergence. There are several optimizations of the SIMPLE algorithm, e.g., the semi-implicit method for pressure-linked equations revised (SIMPLER), the semi-implicit method for pressure-linked equations-consistent (SIMPLEC), and the pressure-implicit with splitting of operator (PISO), which establish a correlation between ṽ_{x,y,z} and $\tilde{p}$ thereby achieving numerical schemes that do not require under-relaxation.

34.2.2.4 Divergent Navier–Stokes Equation

Another common method for pressure term recovery can be done by a small trick. Instead of taking the continuity equation that essentially states that the divergence of the velocity field must be zero, we demand that the divergence of the momentum field must be zero. We can obtain this by simply forming the divergent Navier–Stokes equation. For this we simply take the divergence of the left- and the right-hand side of the Navier–Stokes equation.

We will perform this operation first strictly mathematically and look at the interpretation once we have obtained the result. The divergent Navier–Stokes equation is found by applying the nabla operator to the left-and right-hand side of the Navier–Stokes equation (Eq. 11.40). We find

$\overrightarrow{\nabla }\left(\rho \left(\frac{\partial \overrightarrow{v}}{\partial t}+\left(\overrightarrow{v}\cdot \overrightarrow{\nabla }\right)\overrightarrow{v}\right)\right)=\overrightarrow{\nabla }\left(\overrightarrow{k}-\overrightarrow{\nabla }p\mathit{+}\eta \Delta \overrightarrow{v}\right)$

We will look at two cases of this equation: the stationary and the instationary flows.

Stationary Flows. In stationary flows, the time derivative on the left-hand side of Eq. 34.21 is zero, in which case we simplify Eq. 34.21 to

$\rho \overrightarrow{\nabla }\left(\overrightarrow{v}\cdot \overrightarrow{\nabla }\right)\overrightarrow{v}=\overrightarrow{\nabla }\left(\overrightarrow{k}-\overrightarrow{\nabla }p+\eta \Delta \overrightarrow{v}\right)$

The right-hand side of Eq. 34.22 can be simplified by noting that

$\begin{array}{l}\overrightarrow{\nabla }\left(\overrightarrow{k}-\overrightarrow{\nabla }p+\eta \Delta \overrightarrow{v}\right)=\overrightarrow{\nabla }\overrightarrow{k}-\Delta p\mathit{+}\overrightarrow{\nabla }\left(\eta \Delta \overrightarrow{v}\right)\\ =\overrightarrow{\nabla }\overrightarrow{k}-\Delta p\mathit{+}\eta \overrightarrow{\nabla }\left(\Delta \overrightarrow{v}\right)\\ =\overrightarrow{\nabla }\overrightarrow{k}-\Delta p\mathit{+}\eta \Delta \left(\overrightarrow{\nabla }\overrightarrow{v}\right)\end{array}$

where we have used the fact that the order of the nabla and Laplace operators can be switched. This allows us to use the continuity equation again as we have another term $\overrightarrow{\nabla }\overrightarrow{v}\mathit{=}0$ on the right-hand side of the equation. Also, we can move the viscosity η out of the differentials as we restrict ourselves to Newtonian media for which the viscosity is a constant (see section 9.5.2). This allows us to write the right-hand side as

$\overrightarrow{\nabla }\left(\overrightarrow{k}-\overrightarrow{\nabla }p+\eta \Delta \overrightarrow{v}\right)=\overrightarrow{\nabla }\overrightarrow{k}-\Delta p$

Combining Eq. 34.22 and Eq. 34.23 yields the divergence Navier–Stokes equation:

$\rho \overrightarrow{\nabla }\left(\overrightarrow{v}\cdot \overrightarrow{\nabla }\right)\overrightarrow{v}=\overrightarrow{\nabla }\overrightarrow{k}-\Delta p$

We can rearrange Eq. 34.24 slightly, yielding

$\begin{array}{l}\Delta p=\overrightarrow{\nabla }\overrightarrow{k}-\rho \overrightarrow{\nabla }\left(\overrightarrow{v}\cdot \overrightarrow{\nabla }\right)\overrightarrow{v}\\ =\overrightarrow{\nabla }\left(\overrightarrow{k}-\rho \left(\overrightarrow{v}\cdot \overrightarrow{\nabla }\right)\overrightarrow{v}\right)\end{array}$

Eq. 34.25 is commonly referred to as the pressure-Poisson equation.

Interpretation. We will now look into what exactly the physical implications of this introduction of divergence actually are. If you look at Eq. 34.25 you see that it looks as if we had just applied the nabla operator to an equation that was rearranged such that the pressure gradient was on the left-hand side. Now, the pressure is prepended by the Laplace operator, which usually describes diffusion effects. This is not a bad interpretation. In fact, one can think of the incompressibility as being a constraint for the pressure. However, the pressure is not a dependent variable. By taking the gradient of the Navier–Stokes equation we generate an artificial pressure field. As this field is derived from the Navier–Stokes equation it satisfies this equation implicitly. If you look carefully at Eq. 34.25 you can see that the right-hand side contains the divergence of the volume forces and the convection term. These are the terms that would contribute to the generation of pressure. Volume forces obviously create pressure, as does in- and outflowing mass due to the incompressibility of the fluid. The other terms of the Navier–Stokes equation canceled during the derivation because they do not contribute to the change of pressure. The generated pressure field thus ensures that the Navier–Stokes and the continuity equation are satisfied.

Usage. Using Eq. 34.25 instead of the continuity equation (Eq. 34.4) to recover the pressure term is a very common implementation. In the following we will also make use of this equation to recover the pressure.

Instationary Flows. Having derived the divergent Navier–Stokes equation for stationary flows, we now turn to instationary flows. For these cases, we know that $\frac{\partial \overrightarrow{v}}{\partial t}=0$ , in which case we cannot apply the simplification made in Eq. 34.22. Instead, we need a numerical scheme for replacing the time derivative. We can choose any suitable method discussed in section 27.2 for doing so. For the sake of simplicity, we will use a forward Euler scheme (see section 27.2.2). For this we find

${\overrightarrow{v}}^{\left(t+1\right)}={\overrightarrow{v}}^{\left(t\right)}+h_t\frac{\partial {\overrightarrow{v}}^{\left(t\right)}}{\partial t}$

As you can see, we have implemented an index variable t in time and stepped the velocity forward from t to t + 1 by a forward Euler scheme. Rearranging the Navier–Stokes equation (see Eq. 11.40) we find the first-order time derivative to be

$\frac{\partial \overrightarrow{v}}{\partial t}=-\left(\overrightarrow{v}\cdot \overrightarrow{\nabla }\right)\overrightarrow{v}+\frac{1}{\rho }\overrightarrow{k}-\frac{1}{\rho }\overrightarrow{\nabla }p\mathit{+}\frac{\eta }{\rho }\Delta \overrightarrow{v}$

which we insert into Eq. 34.26 to find

${\overrightarrow{v}}^{\left(t+1\right)}={\overrightarrow{v}}^{\left(t\right)}+h_t\left(-\left({\overrightarrow{v}}^{\left(t\right)}\cdot \overrightarrow{\nabla }\right){\overrightarrow{v}}^{\left(t\right)}+\frac{1}{\rho }{\overrightarrow{k}}^{\left(i\right)}-\frac{1}{\rho }\overrightarrow{\nabla }{p}^{\left(t\right)}\mathit{+}\frac{\eta }{\rho }\Delta {\overrightarrow{v}}^{\left(t\right)}\right)$

As we have done in case of the stationary flows, we now apply the divergence operator to both sides of the equations, in which case we find

$\overrightarrow{\nabla }{\overrightarrow{v}}^{\left(t+1\right)}=\overrightarrow{\nabla }{\overrightarrow{v}}^{\left(t\right)}+h_t\left(-\overrightarrow{\nabla }\left(\left({\overrightarrow{v}}^{\left(t\right)}\cdot \overrightarrow{\nabla }\right){\overrightarrow{v}}^{\left(t\right)}\right)+\frac{1}{\rho }\overrightarrow{\nabla }{\overrightarrow{k}}^{\left(t\right)}-\frac{1}{\rho }\Delta p^{\left(t\right)}\mathit{+}\frac{\eta }{\rho }\overrightarrow{\nabla }\Delta {\overrightarrow{v}}^{\left(t\right)}\right)$

As before, we want to isolate the pressure term on the left-hand side, in which case we find

$\Delta p^{\left(t\right)}=\frac{\rho }{h_t}\left(\overrightarrow{\nabla }{\overrightarrow{v}}^{\left(t\right)}-\overrightarrow{\nabla }{\overrightarrow{v}}^{\left(t+1\right)}\right)-\rho \overrightarrow{\nabla }\left(\left({\overrightarrow{v}}^{\left(t\right)}\cdot \overrightarrow{\nabla }\right){\overrightarrow{v}}^{\left(t\right)}\right)+\overrightarrow{\nabla }{\overrightarrow{k}}^{\left(t\right)}+\eta \Delta \overrightarrow{\nabla }{\overrightarrow{v}}^{\left(t\right)}$

Obviously, Eq. 34.30 looks significantly more complicated than the respective equation we found for the stationary case (see Eq. 34.25). However, there are two small modifications we will make to this equation later, which will make it slightly less complicated to work with. For the moment we leave it as it is.

34.2.3 Schemes for Stationary and Instationary Flows

34.2.3.1 Numerical Schemes

We are now ready to implement the solver. We require Eq. 34.1, Eq. 34.2, Eq. 34.3, and Eq. 34.25 and suitable numerical schemes for them. We will use a stepwidth h_xyz in space and a stepwidth h_t in time. There are two flow cases for which we could implement our solver: stationary and instationary. We will start with the stationary flows first.

34.2.3.2 Stationary Flows

Momentum Equations. For stationary flows, the time derivatives in Eq. 34.1, Eq. 34.2, and Eq. 34.3 are zero and we can directly solve the equations for the unknown velocities v_x, v_y, and v_z. Similar to the approach taken in section 33.3.3.2 we will choose suitable numerical schemes for the first-order partial derivatives and implement an iterative scheme in order to avoid making the left-hand sides of our solutions nonlinear. We will use the following numerical schemes for the pressure terms:

$\frac{\partial p}{\partial x}=\frac{p^{\left(x+1,y,z\right)}-p^{\left(x-1,y,z\right)}}{2h_{xyz}}$

$\frac{\partial p}{\partial y}=\frac{p^{\left(x,y+1,z\right)}-p^{\left(x,y-1,z\right)}}{2h_{xyz}}$

$\frac{\partial p}{\partial z}=\frac{p^{\left(x,y,z+1\right)}-p^{\left(x,y,z-1\right)}}{2h_{xyz}}$

as well as the following schemes for the first-order differentials:

$\frac{\partial v_x}{\partial x}=\frac{{v_x}^{\left(x+1,y,z\right)}-{v_x}^{\left(x-1,y,z\right)}}{2h_{xyz}}$

$\frac{\partial v_x}{\partial y}=\frac{{v_x}^{\left(x,y+1,z\right)}-{v_x}^{\left(x,y-1,z\right)}}{2h_{xyz}}$

$\frac{\partial v_x}{\partial z}=\frac{{v_x}^{\left(x,y,z+1\right)}-{v_x}^{\left(x,y,z-1\right)}}{2h_{xyz}}$

$\frac{\partial v_y}{\partial x}=\frac{{v_y}^{\left(x+1,y,z\right)}-{v_y}^{\left(x-1,y,z\right)}}{2h_{xyz}}$

$\frac{\partial v_y}{\partial y}=\frac{{v_y}^{\left(x,y+1,z\right)}-{v_y}^{\left(x,y-1,z\right)}}{2h_{xyz}}$

$\frac{\partial v_y}{\partial z}=\frac{{v_y}^{\left(x,y,z+1\right)}-{v_y}^{\left(x,y,z-1\right)}}{2h_{xyz}}$

$\frac{\partial v_z}{\partial x}=\frac{{v_z}^{\left(x+1,y,z\right)}-{v_z}^{\left(x-1,y,z\right)}}{2h_{xyz}}$

$\frac{\partial v_z}{\partial y}=\frac{{v_z}^{\left(x,y+1,z\right)}-{v_z}^{\left(x,y-1,z\right)}}{2h_{xyz}}$

$\frac{\partial v_z}{\partial z}=\frac{{v_z}^{\left(x,y,z+1\right)}-{v_z}^{\left(x,y,z-1\right)}}{2h_{xyz}}$

and the following schemes for the second-order differentials

$\frac{{\partial }^2{v}_x}{\partial x^2}=\frac{{v_x}^{\left(x+1,y,z\right)}+{v_x}^{\left(x-1,y,z\right)}-2{v_x}^{\left(x,y,z\right)}}{h_{xyz}^2}$

$\frac{{\partial }^2{v}_x}{\partial y^2}=\frac{{v_x}^{\left(x,y+1,z\right)}+{v_x}^{\left(x,y-1,z\right)}-2{v_x}^{\left(x,y,z\right)}}{h_{xyz}^2}$

$\frac{{\partial }^2{v}_x}{\partial z^2}=\frac{{v_x}^{\left(x,y,z+1\right)}+{v_x}^{\left(x,y,z-1\right)}-2{v_x}^{\left(x,y,z\right)}}{h_{xyz}^2}$

$\frac{{\partial }^2{v}_y}{\partial x^2}=\frac{{v_y}^{\left(x+1,y,z\right)}+{v_y}^{\left(x-1,y,z\right)}-2{v_y}^{\left(x,y,z\right)}}{h_{xyz}^2}$

$\frac{{\partial }^2{v}_y}{\partial y^2}=\frac{{v_y}^{\left(x,y+1,z\right)}+{v_y}^{\left(x,y-1,z\right)}-2{v_y}^{\left(x,y,z\right)}}{h_{xyz}^2}$

$\frac{{\partial }^2{v}_y}{\partial z^2}=\frac{{v_y}^{\left(x,y,z+1\right)}+{v_y}^{\left(x,y,z-1\right)}-2{v_y}^{\left(x,y,z\right)}}{h_{xyz}^2}$

$\frac{{\partial }^2{v}_z}{\partial x^2}=\frac{{v_z}^{\left(x+1,y,z\right)}+{v_z}^{\left(x-1,y,z\right)}-2{v_z}^{\left(x,y,z\right)}}{h_{xyz}^2}$

$\frac{{\partial }^2{v}_z}{\partial y^2}=\frac{{v_z}^{\left(x,y+1,z\right)}+{v_z}^{\left(x,y-1,z\right)}-2{v_z}^{\left(x,y,z\right)}}{h_{xyz}^2}$

$\frac{{\partial }^2{v}_z}{\partial z^2}=\frac{{v_z}^{\left(x,y,z+1\right)}+{v_z}^{\left(x,y,z-1\right)}-2{v_z}^{\left(x,y,z\right)}}{h_{xyz}^2}$

As you can see, we have used a central finite difference scheme (see section 4.2.4.3) for the first-order partial differentials. Until now, we have preferably used backward or forward finite difference schemes (see section 4.2.4.4). The central finite difference scheme has two main advantages, which makes it the more suitable scheme in this application. First, it has no directional bias. Until now, we have primarily discussed flow scenarios where we started from a given end of the domain and “calculated our way forward” toward the end. For these scenarios forward finite difference schemes are most commonly used. We have done the same in section 33.3.3.1. The second advantage of the central difference scheme is the fact that it is numerically more stable due to the fact that it is a second-order scheme. This is advantageous since it is very likely that the solver we are about to implement will be used in cases where numerical stability may become a problem. Therefore, choosing a numerically more stable scheme makes the solver more robust to use. However, the central finite difference scheme comes at a disadvantage, which we will discuss in a moment.

Using the given equations, Eq. 34.1, Eq. 34.2, and Eq. 34.3 can be rewritten as

$\begin{array}{l}{v}_x^{\left(i+1,x,y,z\right)}\frac{v_x^{\left(i,x+1,y,z\right)}-v_x^{\left(i,x-1,y,z\right)}}{2}+v_y^{\left(i,x,y,z\right)}\frac{v_x^{\left(i,x,y+1,z\right)}-v_x^{\left(i,x,y-1,z\right)}}{2}+v_z^{\left(i,x,y,z\right)}\frac{v_x^{\left(i,x,y,z+1\right)}-v_x^{\left(i,x,y,z-1\right)}}{2}\\ =-\frac{1}{\rho }\frac{p^{\left(x+1,y,z\right)}-p^{\left(x-1,y,z\right)}}{2}+\frac{\eta }{\rho h_{xyz}}\left(\begin{array}{l}{v}_x^{\left(i,x+1,y,z\right)}+v_x^{\left(i,x-1,y,z\right)}+v_x^{\left(i,x,y+1,z\right)}+v_x^{\left(i,x,y-1,z\right)}+\\ v_x^{\left(i,x,y,z+1\right)}+v_x^{\left(i,x,y,z-1\right)}-6v_x^{\left(i+1,x,y,z\right)}\end{array}\right)\\ +\frac{h_{xyz}}{\rho }{k}_x^{\left(x,y,z\right)}\end{array}$

$\begin{array}{l}{v}_x^{\left(i,x,y,z\right)}\frac{v_y^{\left(i,x+1,y,z\right)}-v_y^{\left(i,x-1,y,z\right)}}{2}+v_y^{\left(i+1,x,y,z\right)}\frac{v_y^{\left(i,x,y+1,z\right)}-v_y^{\left(i,x,y-1,z\right)}}{2}+v_z^{\left(i,x,y,z\right)}\frac{v_y^{\left(i,x,y,z+1\right)}-v_y^{\left(i,x,y,z-1\right)}}{2}\\ =-\frac{1}{\rho }\frac{p^{\left(x,y+1,z\right)}-p^{\left(x,y-1,z\right)}}{2}+\frac{\eta }{\rho h_{xyz}}\left(\begin{array}{l}{v}_y^{\left(i,x+1,y,z\right)}+v_y^{\left(i,x-1,y,z\right)}+v_y^{\left(i,x,y+1,z\right)}+v_y^{\left(i,x,y-1,z\right)}+\\ v_y^{\left(i,x,y,z+1\right)}+v_y^{\left(i,x,y,z-1\right)}-6v_y^{\left(i+1,x,y,z\right)}\end{array}\right)\\ +\frac{h_{xyz}}{\rho }{k}_y^{\left(x,y,z\right)}\end{array}$

$\begin{array}{l}{v}_x^{\left(i,x,y,z\right)}\frac{v_z^{\left(i,x+1,y,z\right)}-v_z^{\left(i,x-1,y,z\right)}}{2}+v_y^{\left(i,x,y,z\right)}\frac{v_z^{\left(i,x,y+1,z\right)}-v_z^{\left(i,x,y-1,z\right)}}{2}+v_z^{\left(i+1,x,y,z\right)}\frac{v_z^{\left(i,x,y,z+1\right)}-v_z^{\left(i,x,y,z-1\right)}}{2}\\ =-\frac{1}{\rho }\frac{p^{\left(x,y,z+1\right)}-p^{\left(x,y,z-1\right)}}{2}+\frac{\eta }{\rho h_{xyz}}\left(\begin{array}{l}{v}_z^{\left(i,x+1,y,z\right)}+v_z^{\left(i,x-1,y,z\right)}+v_z^{\left(i,x,y+1,z\right)}+v_z^{\left(i,x,y-1,z\right)}+\\ v_z^{\left(i,x,y,z+1\right)}+v_z^{\left(i,x,y,z-1\right)}-6v_z^{\left(i+1,x,y,z\right)}\end{array}\right)\\ +\frac{h_{xyz}}{\rho }{k}_z^{\left(x,y,z\right)}\end{array}$

which we rewrite to

$\begin{array}{l}{v}_x^{\left(i+1,x,y,z\right)}\left(\frac{v_x^{\left(i,x+1,y,z\right)}-v_x^{\left(i,x-1,y,z\right)}}{2}+6\frac{\eta }{\rho h_{xyz}}\right)+v_y^{\left(i,x,y,z\right)}\frac{v_x^{\left(i,x,y+1,z\right)}-v_x^{\left(i,x,y-1,z\right)}}{2}+\\ v_z^{\left(i,x,y,z\right)}\frac{v_x^{\left(i,x,y,z+1\right)}-v_x^{\left(i,x,y,z-1\right)}}{2}\\ =-\frac{1}{\rho }\frac{p^{\left(x+1,y,z\right)}-p^{\left(x-1,y,z\right)}}{2}+\frac{\eta }{\rho h_{xyz}}\left(\begin{array}{l}{v}_x^{\left(i,x+1,y,z\right)}+v_x^{\left(i,x-1,y,z\right)}+v_x^{\left(i,x,y+1,z\right)}+v_x^{\left(i,x,y-1,z\right)}+\\ v_x^{\left(i,x,y,z+1\right)}+v_x^{\left(i,x,y,z-1\right)}\end{array}\right)\\ +\frac{h_{xyz}}{\rho }{k}_x^{\left(x,y,z\right)}\\ \end{array}$

$\begin{array}{l}{v}_x^{\left(i,x,y,z\right)}\frac{v_y^{\left(i,x+1,y,z\right)}-v_x^{\left(i,x-1,y,z\right)}}{2}+v_y^{\left(i+1,x,y,z\right)}\left(\frac{v_y^{\left(i,x,y+1,z\right)}-v_y^{\left(i,x,y-1,z\right)}}{2}+6\frac{\eta }{\rho h_{xyz}}\right)+\\ v_z^{\left(i,x,y,z\right)}\frac{v_y^{\left(i,x,y,z+1\right)}-v_y^{\left(i,x,y,z-1\right)}}{2}\\ =-\frac{1}{\rho }\frac{p^{\left(x,y+1,z\right)}-p^{\left(x,y-1,z\right)}}{2}+\frac{\eta }{\rho h_{xyz}}\left(\begin{array}{l}{v}_y^{\left(i,x+1,y,z\right)}+v_y^{\left(i,x-1,y,z\right)}+v_y^{\left(i,x,y+1,z\right)}+v_y^{\left(i,x,y-1,z\right)}+\\ v_y^{\left(i,x,y,z+1\right)}+v_y^{\left(i,x,y,z-1\right)}\end{array}\right)\\ +\frac{h_{xyz}}{\rho }{k}_y^{\left(x,y,z\right)}\end{array}$

$\begin{array}{l}{v}_x^{\left(i,x,y,z\right)}\frac{v_z^{\left(i,x+1,y,z\right)}-v_z^{\left(i,x-1,y,z\right)}}{2}+v_y^{\left(i,x,y,z\right)}\frac{v_z^{\left(i,x,y+1,z\right)}-v_z^{\left(i,x,y-1,z\right)}}{2}+\\ v_z^{\left(i+1,x,y,z\right)}\left(\frac{v_z^{\left(i,x,y,z+1\right)}-v_z^{\left(i,x,y,z-1\right)}}{2}+6\frac{\eta }{\rho h_{xyz}}\right)\\ =-\frac{1}{\rho }\frac{p^{\left(x,y,z+1\right)}-p^{\left(x,y,z-1\right)}}{2}+\frac{\eta }{\rho h_{xyz}}\left(\begin{array}{l}{v}_z^{\left(i,x+1,y,z\right)}+v_z^{\left(i,x-1,y,z\right)}+v_z^{\left(i,x,y+1,z\right)}+v_z^{\left(i,x,y-1,z\right)}+\\ v_z^{\left(i,x,y,z+1\right)}+v_z^{\left(i,x,y,z-1\right)}\end{array}\right)\\ +\frac{h_{xyz}}{\rho }{k}_z^{\left(x,y,z\right)}\end{array}$

As you can see, we implemented an iteration variable i. The respective dependent variables to solve for are $v_x^{\left(i+1,x,y,z\right)}$ , $v_y^{\left(i+1,x,y,z\right)}$ , and $v_z^{\left(i+1,x,y,z\right)}$ . We avoided making the equations nonlinear by choosing the central finite difference scheme for the first-order partial derivatives, which only depends on the neighboring values of the iteration i; we use this to update the value of the iteration i + 1 of the convection terms on the left-hand side. As you can see, we implemented the simplification that we already introduced during the discussion of the SIMPLE algorithm, i.e., we do not treat Eq. 34.31, Eq. 34.32, and Eq. 34.33 as a linear system of equations but simply solve them sequentially. Obviously, the equations are interdependent and the left-hand sides should be written (for all three equations) using $v_x^{\left(i+1,x,y,z\right)}$ , $v_y^{\left(i+1,x,y,z\right)}$ , and $v_z^{\left(i+1,x,y,z\right)}$ instead of only containing the respective variable of iteration i + 1. However, we need to iterate the solution anyways and it therefore is fine to solve the equations one by one simply updating the values. Solving Eq. 34.31, Eq. 34.32, and Eq. 34.33 yields

$v_x^{\left(i+1,x,y,z\right)}=\frac{1}{\frac{v_x^{\left(i,x+1,y,z\right)}-v_x^{\left(i,x-1,y,z\right)}}{2}+6\frac{\eta }{\rho h_{xyz}}}\left(\begin{array}{l}-\left(\begin{array}{l}{v}_y^{\left(i,x,y,z\right)}\frac{v_x^{\left(i,x,y+1,z\right)}-v_x^{\left(i,x,y-1,z\right)}}{2}+\\ v_z^{\left(i,x,y,z\right)}\frac{v_x^{\left(i,x,y,z+1\right)}-v_x^{\left(i,x,y,z-1\right)}}{2}\end{array}\right)-\\ \frac{1}{\rho }\frac{p^{\left(x+1,y,z\right)}-p^{\left(x-1,y,z\right)}}{2}+\frac{\eta }{\rho h_{xyz}}\left(\begin{array}{l}{v}_x^{\left(i,x+1,y,z\right)}+v_x^{\left(i,x-1,y,z\right)}+\\ v_x^{\left(i,x,y+1,z\right)}+v_x^{\left(i,x,y-1,z\right)}+\\ v_x^{\left(i,x,y,z+1\right)}+v_x^{\left(i,x,y,z-1\right)}\end{array}\right)+\frac{h_{xyz}}{\rho }{k}_x^{\left(x,y,z\right)}\end{array}\right)$

$v_y^{\left(i+1,x,y,z\right)}=\frac{1}{\frac{v_y^{\left(i,x,y+1,z\right)}-v_y^{\left(i,x,y-1,z\right)}}{2}+6\frac{\eta }{\rho h_{xyz}}}\left(\begin{array}{l}-\left(\begin{array}{l}{v}_x^{\left(i,x,y,z\right)}\frac{v_y^{\left(i,x+1,y,z\right)}-v_y^{\left(i,x-1,y,z\right)}}{2}+\\ v_z^{\left(i,x,y,z\right)}\frac{v_y^{\left(i,x,y,z+1\right)}-v_y^{\left(i,x,y,z-1\right)}}{2}\end{array}\right)-\\ \frac{1}{\rho }\frac{p^{\left(x,y+1,z\right)}-p^{\left(x,y-1,z\right)}}{2}+\frac{\eta }{\rho h_{xyz}}\left(\begin{array}{l}{v}_y^{\left(i,x+1,y,z\right)}+v_y^{\left(i,x-1,y,z\right)}+\\ v_y^{\left(i,x,y+1,z\right)}+v_y^{\left(i,x,y-1,z\right)}+\\ v_y^{\left(i,x,y,z+1\right)}+v_y^{\left(i,x,y,z-1\right)}\end{array}\right)+\frac{h_{xyz}}{\rho }{k}_y^{\left(x,y,z\right)}\end{array}\right)$

$v_z^{\left(i+1,x,y,z\right)}=\frac{1}{\frac{v_z^{\left(i,x,y,z+1\right)}-v_z^{\left(i,x,y,z-1\right)}}{2}+6\frac{\eta }{\rho h_{xyz}}}\left(\begin{array}{l}-\left(\begin{array}{l}{v}_x^{\left(i,x,y,z\right)}\frac{v_z^{\left(i,x+1,y,z\right)}-v_z^{\left(i,x-1,y,z\right)}}{2}+\\ v_y^{\left(i,x,y,z\right)}\frac{v_z^{\left(i,x,y+1,z\right)}-v_y^{\left(i,x,y-1,z\right)}}{2}\end{array}\right)-\\ \frac{1}{\rho }\frac{p^{\left(x,y,z+1\right)}-p^{\left(x,y,z-1\right)}}{2}+\frac{\eta }{\rho h_{xyz}}\left(\begin{array}{l}{v}_z^{\left(i,x+1,y,z\right)}+v_z^{\left(i,x-1,y,z\right)}+\\ v_z^{\left(i,x,y+1,z\right)}+v_z^{\left(i,x,y-1,z\right)}+\\ v_z^{\left(i,x,y,z+1\right)}+v_z^{\left(i,x,y,z-1\right)}\end{array}\right)+\frac{h_{xyz}}{\rho }{k}_z^{\left(x,y,z\right)}\end{array}\right)$

Eq. 34.34, Eq. 34.35, and Eq. 34.36 are the required momentum equations.

Pressure-Poisson Equation. The last equation we require is the pressure-Poisson equation (Eq. 34.25), which we rewrite to

$\begin{array}{l}\Delta p=\overrightarrow{\nabla }\overrightarrow{k}-\rho \overrightarrow{\nabla }\left(\overrightarrow{v\cdot }\overrightarrow{\nabla }\right)\overrightarrow{v}\\ =\frac{\partial k_x}{\partial x}+\frac{\partial k_y}{\partial y}+\frac{\partial k_z}{\partial z}-\rho \overrightarrow{\nabla }\left(v_x\frac{\partial }{\partial x}+v_y\frac{\partial }{\partial y}+v_z\frac{\partial }{\partial z}\right)\overrightarrow{v}\\ =\frac{\partial k_x}{\partial x}+\frac{\partial k_y}{\partial y}+\frac{\partial k_z}{\partial z}-\rho \left(\begin{array}{l}{v}_x\frac{{\partial }^2}{\partial x\partial x}+\frac{\partial v_x}{\partial x}\frac{\partial }{\partial x}+v_y\frac{{\partial }^2}{\partial y\partial x}+\frac{\partial v_y}{\partial x}\frac{\partial }{\partial y}+v_z\frac{{\partial }^2}{\partial z\partial x}+\frac{\partial v_z}{\partial x}\frac{\partial }{\partial z}\\ v_x\frac{{\partial }^2}{\partial x\partial y}+\frac{\partial v_x}{\partial y}\frac{\partial }{\partial x}+v_y\frac{{\partial }^2}{\partial y\partial y}+\frac{\partial v_y}{\partial y}\frac{\partial }{\partial y}+v_z\frac{{\partial }^2}{\partial z\partial y}+\frac{\partial v_z}{\partial y}\frac{\partial }{\partial z}\\ v_x\frac{{\partial }^2}{\partial x\partial z}+\frac{\partial v_x}{\partial z}\frac{\partial }{\partial x}+v_y\frac{{\partial }^2}{\partial y\partial y}+\frac{\partial v_y}{\partial z}\frac{\partial }{\partial y}+v_z\frac{{\partial }^2}{\partial z\partial z}+\frac{\partial v_z}{\partial z}\frac{\partial }{\partial z}\end{array}\right)\left(\begin{array}{c}\hfill v_x\hfill \\ \hfill v_y\hfill \\ \hfill v_z\hfill \end{array}\right)\\ =\frac{\partial k_x}{\partial x}+\frac{\partial k_y}{\partial y}+\frac{\partial k_z}{\partial z}-\rho \left(\begin{array}{l}{v}_x\frac{{\partial }^2{v}_x}{\partial x\partial x}+\frac{\partial v_x}{\partial x}\frac{\partial v_x}{\partial x}+v_y\frac{{\partial }^2{v}_x}{\partial y\partial x}+\frac{\partial v_y}{\partial x}\frac{\partial v_x}{\partial y}+v_z\frac{{\partial }^2{v}_x}{\partial z\partial x}+\frac{\partial v_z}{\partial x}\frac{\partial v_x}{\partial z}+v_x\frac{{\partial }^2{v}_y}{\partial x\partial y}+\frac{\partial v_x}{\partial y}\frac{\partial v_y}{\partial x}+\\ v_y\frac{{\partial }^2{v}_y}{\partial y\partial y}+\frac{\partial v_y}{\partial y}\frac{\partial v_y}{\partial y}+v_z\frac{{\partial }^2{v}_y}{\partial z\partial y}+\frac{\partial v_z}{\partial y}\frac{\partial v_y}{\partial z}+v_x\frac{{\partial }^2{v}_z}{\partial x\partial z}+\frac{\partial v_x}{\partial z}\frac{\partial v_z}{\partial x}+\\ v_y\frac{{\partial }^2{v}_z}{\partial y\partial z}+\frac{\partial v_y}{\partial z}\frac{\partial v_z}{\partial y}+v_z\frac{{\partial }^2{v}_z}{\partial z\partial x}+\frac{\partial v_z}{\partial y}\frac{\partial v_z}{\partial z}\end{array}\right)\\ =\frac{\partial k_x}{\partial x}+\frac{\partial k_y}{\partial y}+\frac{\partial k_z}{\partial z}-\rho \left(\begin{array}{l}{v}_x\left(\frac{{\partial }^2{v}_x}{\partial x\partial x}+\frac{{\partial }^2{v}_y}{\partial x\partial x}+\frac{{\partial }^2{v}_z}{\partial x\partial z}\right)+v_y\left(\frac{{\partial }^2{v}_x}{\partial y\partial x}+\frac{{\partial }^2{v}_y}{\partial y\partial y}+\frac{{\partial }^2{v}_z}{\partial y\partial z}\right)+\\ v_z\left(\frac{{\partial }^2{v}_x}{\partial z\partial x}+\frac{{\partial }^2{v}_y}{\partial z\partial y}+\frac{{\partial }^2{v}_z}{\partial z\partial z}\right)+{\left(\frac{\partial v_x}{\partial x}\right)}^2+{\left(\frac{\partial v_y}{\partial y}\right)}^2+{\left(\frac{\partial v_z}{\partial z}\right)}^2+\\ 2\left(\frac{\partial v_y}{\partial x}\frac{\partial v_x}{\partial y}+\frac{\partial v_x}{\partial z}\frac{\partial v_z}{\partial x}+\frac{\partial v_y}{\partial z}\frac{\partial v_z}{\partial y}\right)\end{array}\right)\\ \end{array}$