What is y+ and why is it important?

In setting up cases for Reynolds averaged Navier-Stokes (RANS) analysis for CFD, we usually incorporate a prism layer along wall surfaces to capture boundary later effects in better detail. The technique in such cases would be to size the first layer of mesh for a target y+ value.

The image above I got from here shows an example. When searching for y+ online, one commonly finds the definition:

$\begin{align*} y^+ = \frac{y u_{\tau}}{\nu} \quad , \quad u_{\tau}^2 = \frac{\tau_w}{\rho} = \nu \frac{dU}{dy}\rvert_{y=0} . \end{align*}$

$\tau_w$ can be found with empirical relations, for example by using the $1/7$ law, $C_f = 0.026/\left(Re_x^{1/7}\right)$ . Then, $\tau_w = (1/2) \rho U^2 C_f$ . The Reynolds number $Re_x$ is usually evaluated with the characteristic length like chord for an airfoil. There are some snazzy websites where you can easily compute the cell length for your desired y+ value like here and here. These are all very well and good, but what exactly is this y+ we are calculating and what dictates it’s value?

Before we proceed, it would be timely to visit the post on RANS and the mixing length model as we will be using knowledge on both of them for our discussion later on.

RANS flow in a channel

Suppose we have a two-dimensional pressure driven flow in a channel. The flow is steady and fully developed. This means that the $\partial / \partial t$ and $\partial / \partial x$ terms are all zero except for $\partial p / \partial x$ . The flow is also assumed to be parallel to the channel, so $V=0$ and $U = U(y)$ . From this, we get for the $x$ direction RANS:

$\begin{align*} \underbrace{ \frac{\partial U}{\partial t} }_{\partial / \partial t = 0} + \underbrace{ U \frac{\partial}{\partial x} U }_{\partial U / \partial x = 0}+ \underbrace{ V \frac{\partial}{\partial y} U }_{V = 0} = -\frac{1}{\rho} \frac{\partial}{\partial x} P + \underbrace{ \nu \frac{\partial^2}{\partial x^2} U }_{\partial U / \partial x = 0} + \nu \frac{\partial^2}{\partial y^2} U - \underbrace{ \frac{\partial}{\partial x} \left< u' u'\right> }_{\partial \left<u'u'\right> / \partial x = 0} - \frac{\partial}{\partial y} \left< u' v'\right> \\ \Rightarrow -\frac{1}{\rho} \frac{\partial}{\partial x} P + \nu \frac{\partial^2}{\partial y^2} U - \frac{\partial}{\partial y} \left< u' v'\right> = 0 . \end{align*}$

In the $y$ direction,

$\begin{align*} \underbrace{ \frac{\partial V}{\partial t} }_{V=0} + \underbrace{ U \frac{\partial}{\partial x} V }_{V=0}+ \underbrace{ V \frac{\partial}{\partial y} V }_{V=0} = -\frac{1}{\rho} \frac{\partial}{\partial y} P + \underbrace{ \nu \frac{\partial^2}{\partial x^2} V }_{V=0} + \underbrace{ \nu \frac{\partial^2}{\partial y^2} V }_{V=0}- \underbrace{ \frac{\partial}{\partial x} \left< u' v'\right> }_{\partial \left<u'v'\right> / \partial x = 0} - \frac{\partial}{\partial y} \left< v' v'\right>\\ \Rightarrow -\frac{1}{\rho} \frac{\partial}{\partial y} P - \frac{\partial}{\partial y} \left< v' v'\right> = 0 . \end{align*}$

Differentiating the $y$ direction equation in $x$ , then integrating in $y$ ,

$\begin{align*} -\frac{1}{\rho} \frac{\partial}{\partial x} \frac{\partial}{\partial y} P - \frac{\partial}{\partial x} \frac{\partial}{\partial y} \left< v' v'\right> &= 0 \\ -\frac{1}{\rho} \frac{\partial}{\partial y} \frac{\partial}{\partial x} P - \frac{\partial}{\partial y} \underbrace{\frac{\partial}{\partial x} \left< v' v'\right>}_{\partial \left<v'v'\right> / \partial x = 0} &= 0 \\ \\ \int \frac{\partial}{\partial y} \frac{\partial}{\partial x} P dy &= const \\ \Rightarrow \frac{\partial}{\partial x} P &= const . \end{align*}$

This is helpful because we can now integrate the $x$ equation in $y$ :

(1) $\begin{align*} \int _0^y -\frac{1}{\rho} \frac{\partial}{\partial x} P + \nu \frac{\partial^2}{\partial y^2} U - \frac{\partial}{\partial y} \left< u' v'\right> dy &= 0 \\ -\frac{1}{\rho} \frac{\partial P}{\partial x} y + \nu \frac{\partial U}{\partial y} \rvert_y - \nu \frac{\partial U}{\partial y} \rvert_0 - \left< u' v'\right> \rvert_y + \underbrace{\left< u' v'\right> \rvert_0}_{0} &= 0 . \end{align*}$

Here we make the argument that the flow is symmetrical at $y=H$ . Imagine we have a mirror at that location, and a fluid particle is moving towards it from the bottom. There will be another particle approaching the mirror from the top, and the $y$ component velocities will cancel. This means $V=v'=0$ at $y=H$ . Also, another more obvious observation is that $\partial U / \partial y=0$ at $H$ . Now since we defined our integration limits to an arbitrary $y$ , we can let it be $H$ :

$\begin{align*} -\frac{1}{\rho} \frac{\partial P}{\partial x} H + \nu \underbrace{ \frac{\partial U}{\partial y} \rvert_H}_{0} - \nu \frac{\partial U}{\partial y} \rvert_0 &- \underbrace{ \left< u' v'\right> \rvert_H}_{0} + \underbrace{\left< u' v'\right> \rvert_0}_{0} = 0 \\ \Rightarrow \nu \frac{\partial U}{\partial y} \rvert_0 &= -\frac{1}{\rho} \frac{\partial P}{\partial x} H \\ \frac{u_{\tau}^2}{H} &= -\frac{1}{\rho} \frac{\partial P}{\partial x} . \end{align*}$

This is a very useful result since it helps us get rid of the pressure term in equation 1. Substituting the results into equation 1:

(2) $\begin{align*} \frac{u_{\tau}^2}{H} y + \nu \frac{\partial U}{\partial y} - u_{\tau}^2 - \left< u' v'\right>&= 0 \\ \nu \frac{\partial U}{\partial y} - \left< u' v'\right> &= u_{\tau}^2\left(1-\frac{y}{H}\right) \end{align*}$

Whew! After all these manipulations (and latex typing), we finally get an equation for our turbulent flow in a channel. Let’s zoom in into the case when $y$ is very close to the wall (i.e. $y/H \ll 1$ ).

Laminar sublayer

As you can already guess, the first thing we can do is to neglect the $y/H$ term when we are close to the wall. This still does not give us anything useful. At the wall, our velocities go to zero. So we can also assume $\left<u'v'\right> \approx 0$ . This leaves us with:

$\begin{align*} \nu \frac{\partial U}{\partial y} &= u_{\tau}^2 . \end{align*}$

Dividing both sides by $u_{\tau}^2$ and introducing the nondimensional velocities and wall distances (as defined by the terms in the parenthesis):

$\begin{align*} \frac{\partial \left(U / u_{\tau}\right)}{\partial \left(y u_{\tau}/\nu\right)} &= 1 \\ \frac{\partial u^+}{\partial y^+} &= 1 . \end{align*}$

Integrating to an arbitrary $y$ as we did before:

$\begin{align*} \int _0^y \frac{\partial u^+}{\partial y^+} dy^+ &= y^+\rvert_y \\ u^+\rvert_y - \underbrace{u^+\rvert_0}_0 &= y^+\rvert_y \\ \\ \therefore u^+ &= y^+ . \end{align*}$

Finally, a useful result we get from all the math we did. What does this simple expression $u^+ = y^+$ mean? Despite its simplicity, it is a powerful result – it tells us that no matter what flow conditions you are in or how large the channel is (i.e. it can extend to infinity and we essentially get a flat plate), the velocity profile close to the wall will always fall into this distribution. It is a universal result. Note that since we ignored the Reynolds stresses, this velocity profile holds for laminar flows too. So now the question is how far away from the wall can we get before the results breakdown?

On the left, from a plot I got here, we see the dashed lines showing the different derived velocity profiles, while the solid line shows a blended velocity profile across all regions. In the diagram on the right, which I took from very a useful set of notes by Brennen, we see plots of actual results and pretty good match in results. Typically, the $u^+ = y^+$ profile is valid for $y^+<5$ . Now we start to appreciate why we need to go through the trouble of sizing our prism layers properly.

Generally there are 2 approaches to RANS simulation: a “low y+” simulation or a “high y+” one. This corresponds to the nondimensional height of the first prism layer closest to the wall. For the “low y+” approach, we are calculating the entire laminar sublayer. To ensure that we capture the velocity profile properly, the guide is normally for $y^+ \leq 1$ for the first layer. However, this can result in a large cell count, which may be too expensive for some scenarios. Hence, turbulence modelers create wall functions such that the entire laminar sublayer up to the log law region is entirely modeled in the first prism layer, reducing the required mesh size. These typically require $30 \leq y^+ \leq 100$ , which can be a bit tricky to achieve in simulations and might require some trial and error.

Log-law layer

Alright, let’s move further up the boundary layer into the turbulent zone. to see where the restriction of $30 \leq y^+ \leq 100$ for the “high y+” simulation come from. In this region above the laminar sublayer, Prandtl assumed that there is constant shear stress across the boundary layer. Recall that total shear stress (laminar plus Reynolds shear stress) is:

$\begin{align*} \tau &= \nu \frac{\partial U}{\partial y} - \left<u'v'\right> \end{align*}$

Constant shear stress means that the $y$ derivative of $\tau$ is zero, so,

(3) $\begin{align*} \frac{\partial }{\partial y} \left[ \nu \frac{\partial U}{\partial y} - \left<u'v'\right> \right] = 0. \end{align*}$

Now, applying the mixing length model,

$\begin{align*} \left<u'v'\right> &= \left< \ell \frac{\partial U}{\partial y} \times \ell \frac{\partial V}{\partial y} \right> \\ \left<u'v'\right> &= \ell^2 \left( \frac{\partial U}{\partial y} \right)^2 \end{align*}$

We have assumed that the variation of $v'$ follows that of $u'$ . Prandtl made one further assumption that $\ell \sim y$ , so we let $\ell = \kappa y$ where $\kappa$ is a constant. Now substituting our results into equation 3 and integrating to an arbitrary $y$ ,

$\begin{align*} \underbrace{\nu \frac{\partial U}{\partial y} \rvert_y}_{\ll 1} - \nu \frac{\partial U}{\partial y} \rvert_0 &- \kappa^2 y^2 \left( \frac{\partial U}{\partial y} \right)^2 \rvert_y + \underbrace{\kappa^2 y^2 \left( \frac{\partial U}{\partial y} \right)^2 \rvert_0}_{0} = 0 \\ \Rightarrow \kappa^2 y^2 \left( \frac{\partial U}{\partial y} \right)^2 &= -u_{\tau}^2 \\ \kappa y \frac{\partial U}{\partial y} &= -u_{\tau} \end{align*}$

Nondimensionalizing the variables as usual, and absorbing the negative sign into $\kappa$ ,

$\begin{align*} \frac{u_{\tau}^2}{\nu} \frac{\partial \left( U / u_{\tau}\right)}{\partial \left( y u_{\tau}/\nu\right)} &= \frac{1}{\kappa y} u_{\tau} \\ \frac{\partial \left( U / u_{\tau}\right)}{\partial \left( y u_{\tau}/\nu\right)} &= \frac{1}{\kappa y u_{\tau} / \nu} \\ \\ \therefore \frac{\partial u^+}{\partial y^+} &= \frac{1}{\kappa y^+} . \end{align*}$

After integration, we finally obtain:

$\begin{align*} u^+ = \frac{1}{\kappa}\log y^+ + B . \end{align*}$

Investigations with experimental results over the years have shown that $\kappa \approx 0.41$ and $B\approx 5.1$ . They generally hold for $30 \leq y^+ \leq 100$ , hence the restriction we have for “high y+” simulations. As a final note, I must mention that advances in CFD software is relaxing the constraint on these requirements, with options like enhanced wall functions in Ansys. However, they are mostly present in propriety software and might not be found in all codes.

What is y+ and why is it important?

Leave a reply Cancel reply