The y-plus (y+) value and its importance for mesh generation (2024)

All wall-bounded turbulent flows require knowledge about the y-plus (y+) values in our simulation, in order to accurately predict turbulent boundary layers and thus quantities such as drag coefficients and separation points. In this article, we take a deep dive into the significance of the y-plus (y+) value and how we can use it in our CFD simulations.

Table of Contents

A primer on turbulent boundary layers

When you are starting out with CFD, chances are you will sooner or later have to deal with turbulent flows, unless you are dealing with very specific niche applications in the stokes flow regime (for example, microfluid devices). Turbulence itself is an endless rabbit hole you can easily get yourself lost in but one of the most fundamental concepts that is of practical use to CFD practitioners is that of the y+ value. In this post, we will look at what the y+ value represents and why it is important for CFD simulations.

Let’s start by considering the flow past a flat plate which is a simplified version of what you would see in any wall-bounded flow. A schematic version of this is shown below.

At first, our flow is laminar and forms well-aligned streamlines in our boundary layer, which is forming as a result of shear stresses at the wall. The freestream velocity which hits the flat plate has to slow down to zero at the wall due to the no-slip condition and this causes our boundary layer to grow in the flow direction, here shown as the x-axis. As the Navier-Stokes equations are inherently non-linear in nature, oscillations are bound to arise. For small boundary layers, these can be damped thanks to a relatively speaking thick viscous sublayer that forms near the wall. Once our boundary layer starts to become larger, the relative thickness of that viscous sublayer diminishes and, as a result, it becomes more susceptible to fluctuations and oscillations and starts to transition to a more turbulent-like flow. After some further distance is covered by the boundary layer, it transitions to a fully turbulent flow.

The mechanism by which a flow transitions from a laminar to a turbulent flow can vary, but broadly speaking, there are 4 types of mechanisms that can cause the transition:

• Natural transition
• Bypass transition
• Crossflow transition
• Separation-induced transition

The most common type you’ll encounter are natural and bypass transition. For a natural transition, the flow generally contains very little fluctuations to begin with and so any instabilities within the boundary themselves can grow and then cause a transition from laminar to turbulence. These fluctuations naturally appear in the boundary layer due to the non-linear behaviour of the Navier-Stokes equations, hence the name natural. Bypass transition, on the other hand, is characterised by a freestream flow which contains a sufficient level of fluctuations; these are then powerful enough to penetrate into the boundary layer from the outside and cause the boundary layer to transition.

Actually, the schematic shown above (and in any other textbook, for that matter!) is misleading and does not really accurately reflect how flow transitions from a laminar to a turbulent state. Take a look at the animation shown below, which focuses on what is labelled the transition region above.

Here, we can see the viscous sublayer, to some extent visualised by the darker blue surfaces. On top of that, we see streaks forming, shown by the green to yellow surfaces. These streaks are typically a sign of bypass transition, where at some location in the boundary layer the viscous sublayer can’t keep the flow aligned and small spots of turbulence appear. These are typically not energetic enough by themselves to cause the transition to a fully turbulent flow and they may disappear and the flow becomes locally laminar again. As the boundary layer grows further downstream, more of these spots are generated and combined and together they are then able to transition into a fully turbulent flow.

Empirical research into turbulent boundary layers done by Schlichting in the 1970s suggests that the height of a laminar and turbulent boundary layer can be described by the following equations

\overline{\delta(x)}_{laminar}=5\frac{x}{\sqrt{Re_x}}

\overline{\delta(x)}_{turbulent}=0.37\frac{x}{Re_x^{0.2}}

Thus, we can see, that for different flow types and Reynolds numbers, we would expect to get a different boundary layer height. This puts us in a difficult position if we wanted to model boundary flows using CFD; how many cells should we place inside the boundary layer to sufficiently resolve it? We have looked at turbulent flows past flat plates above, but what about more complicated arrangements? Ideally, we would like to have a uniform way to describe and view boundary layers, and this is what the y-plus (y+) value helps us to achieve.

The quest for generality- how to describe boundary layers uniformly

As alluded to above, turbulent boundary layers come in all sorts of shapes and forms. Let’s try to think of the end goal here first; if we were to be able to describe all boundary layers in the same way, that is, if we can find a way to collapse all boundary layers onto themselves, we would only need to work out things such as mesh requirements once and then apply them for our particular boundary layer flow (e.g. scale our mesh accordingly). This is where the y-plus (y+) value comes in.

Let’s consider two different boundary layer flows depicted below

We saw that, apart from the downstream location x at which we consider the height of our boundary layer, the only other parameter influencing the height of the boundary layer is the Reynolds number itself, which we can write as

Re_x=\frac{Ux}{\nu}

Let’s look at the parameters involved:

• x: As we travel further downstream of our boundary layer, the Reynolds number Re_x – here defined on the downstream location x – is going to increase monotonically. However, referring back to the two boundary layers shown above, for the same location x, we can have different heights of the boundary layer and therefore x itself is not a good choice for trying to describe general boundary layers.
• \nu: The kinematic viscosity is a material property and we can’t influence it. For most flows it is sufficient to consider the viscosity to be constant and thus, regardless of where we are in the boundary layer, we can’t associate a change in the boundary layer with a change in viscosity. This parameter is also not suitable to try to describe a general boundary layer.
• U: Finally we have the velocity. But which one to take? There are two characteristic velocities we can specify; one is at the wall and the other is in the freestream. But you can probably see that neither is going to work. At the wall, the velocity is zero (due to the no-slip condition) and in the freestream, it is some finite, non-zero value (if it was zero, there wouldn’t be any boundary layer to begin with). Even if this value is non-zero, it is constant and a change in boundary layer height could not be attributed to a change in the freestream velocity. So we seem to have reached a deadlock here and the only way for us to find a suitable value to describe changes in boundary layers has to be constructed in some form, which we will discuss below.

Even if the velocity itself can’t be used, there is one quantity related to the velocity that does offer information about changes in the boundary layer and that is its gradient, i.e. the velocity gradient \partial U/\partial y. We can see that for the same location x and the same assumed kinematic viscosity \nu, changes in the boundary layer height still produce different velocity gradients. Thus, if we were to use the velocity gradient, then we would have the means to describe changes in the boundary layer locally, which is what we are after. There is then, however, the issue that the velocity gradient does not have the same units as velocity and so we need to modify it first. In order to do so, we introduce the wall shear stresses \tau_w as

\tau_w=\rho\nu\frac{\partial U}{\partial y}

We can see that the wall shear stresses already include the gradient we have identified as the metric that we want to track. Thus, let’s take the wall shear stresses and see how we can modify them to have units of velocity which is then a velocity we can use to describe arbitrary boundary layers.

Looking at the units for the wall shear stresses, we have \tau_w=[N/m^2]=[kg/(m s^2)]. We already have meters and seconds in there (not necessarily in the right place) but we also have gained kilograms kg. To get rid of it, we would need to divide by some mass, but if we think about it, there is no specific mass in our flow we can attribute here, however, we do have the mass per volume, or density, in units of kg/m^3, thus dividing the wall shear stresses by density gives us

\frac{\tau_w}{\rho}=\nu\frac{\partial U}{\partial y}=\left[\frac{m^2}{s^2}\right]

This looks promising, we have the right units but squared, taking the square root provides us with

u_*=\sqrt{\frac{\tau_w}{\rho}}=\left[\frac{m}{s}\right]

Thus, we have found a velocity that is capable to describe changes in the boundary layer locally, as it has the velocity gradient information encoded through the wall shear stresses. This velocity has its own name and is referred to as the friction velocity. This name makes somewhat sense, as the wall shear stresses essentially produce friction on the surface (and we use the same wall shear stresses to define the skin friction coefficient c_f as c_f=\tau/w/q, where q is the dynamic pressure).

Let’s step back at this point and see what we have achieved; for any location x in our boundary layer, we have found a velocity u_* that describes locally the shape of the boundary layer. We found this value by examining which values describe the boundary layer height itself (through Schlichting’s empirical correlations for the boundary layer height) which showed us that the Reynolds number is the main scaling factor and through it, we found the velocity gradient and thus the wall shear stresses to be the main driver for changes to the boundary layer height. The friction velocity provides us with a convenient tool to describe then this boundary layer. Now we want to return to our original question, can we find a way to collapse all boundary layers onto themselves? Doing so would require some non-dimensional coordinate system in which all boundary layers would look the same. This would also mean that the wall-normal direction (here y, which is typically chosen as the wall-normal direction in the turbulence literature) would need to be non-dimensionalised in some form. Let us call this new, non-dimensional wall-normal direction y^+. In order to non-dimensionalise y, we take the friction velocity – which we have already defined to be representative of our local boundary shape and size – and see what we need to multiply it with in order to get units of [1/m]. If we have found such a quantity, then we can multiply that with the y coordinate and have found a non-dimensional form of y.

Since u_* has units of [m/s], in order to get to [1/m], we need to multiply it with some properties that provide us with units of [s/m^2], as this would provide us with

\left[\frac{m}{s}\right]\cdot\left[\frac{s}{m^2}\right]=\left[\frac{1}{m}\right]

As it so happens, we already have a property that gives us these units, just in its reciprocal value, which is the kinematic velocity \nu=\left[\frac{m^2}{s}\right]. Thus, by dividing it by the friction velocity and then multiplying that with the y coordinate, we have found to non-dimensionalise our wall-normal direction as

y^+=y\cdot\frac{u_*}{\nu}

This is our non-dimensional wall-normal direction and we have termed it y+, as alluded to above already. But what have we gained with this knowledge? Well, let’s see. Referring back to the image shown above of the turbulent boundary layer, we now have a way to plot our turbulent boundary layer in a non-dimensional form in the wall-normal direction. If we do the same for the velocity, then we can plot non-dimensional velocity profiles and see how they behave for different flows. To non-dimensionalise our velocity, we already have done the heavy lifting and found the friction velocity, thus, we can define – in analogy to our non-dimensional y^+ value – a non-dimensional velocity u^+ as

u^+=\frac{U}{u_*}

Then what we can do is plot u^+ over y^+ and observe the boundary layer velocity profile. This is reproduced below.

This plot is given in a semi-log scale (y^+ is in log scale while u^+ is given on a linear scale). Let’s investigate this figure a bit further. We can see two different regimes; close to the wall (i.e. small y^+ values), we see that we have a linear velocity profile where we have u^+=y^+. Due to the semi-log nature of the plot, this linear relationship appears exponentially in the plot and this linear relationship is experienced approximately until y^+\approx 5. Then, at some further distance away from the wall, typically around y^+\approx 30, the behaviour changes from a linear to a logarithmic one, here given as u^+=1/\kappa ln(y^+)+B. And again, any logarithmic relationship in this semi-log plot will appear linearly. The values for \kappa and B are determined in experiments through measurements (of note here is \kappa\approx 0.41 which is named after Theodore von Karman). Between this linear and logarithmic profile is a transition region where the velocity gradient changes its behaviour gradually.

So what’s the deal with this plot then? Well, take any boundary layer flow, not just flat plates, measure the velocity profile, either in CFD or in experiments, and you will get the same plot as above (well, sometimes you have to play around with the B value to collapse all graphs onto each other but they all experience the same behaviour, i.e. having a linear and logarithmic profile with a blend in-between at the same position). This universal behaviour is of such importance that we have given it its own name and we refer to it in the literature as the law of the wall. Getting to this point required knowledge of the y^+ value (and, by extension, of the u^+ value). But apart from drawing pretty graphs that all collapse on top of each other, what practical use can we derive from this? That is the question we will answer in the next section.

This website exists to create a community of like-minded CFD enthusiasts and I’d love to start a discussion with you. If you would like to be part of it, sign up using the link below and you will receive my OpenFOAM quick reference guide, as well as my guide on Tools every CFD developer needs for free.

Join now

Y-plus (y+) and mesh generation

Now that we have understood what the y^+ value is, let us now turn to a practical application of this value; mesh generation. In mesh generation, especially in cases where we want to resolve turbulence and we have solid walls and expect wall-bounded turbulent boundary layers, the y^+ value plays an important role and we have to match this in our mesh generation stage in order to resolve the turbulence itself. Let us now look at different turbulence modelling approaches and see what our y^+ requirements are.

Direct Numerical Simulations (DNS)

Turbulence is a tricky thing in CFD; there are two ways of dealing with it, either to model it (entirely or partially) or to resolve it. Regardless of the approach, it is a complex endeavour. Depending on the approach chosen, though, the complexity shifts from one domain to another. If we choose to resolve all turbulence, then we have essentially zero modelling requirements and resolving all turbulence is as simple as solving the Navier-Stokes equations as they are. The complexity, however, introduces itself in the stringent mesh requirements and the amount of computational resources required to resolve these turbulence scales. Direct numerical simulations, or DNS, as the name suggests, resolve all turbulence and thus while writing a DNS Navier-Stokes code is relatively straightforward, we put rather drastic restrictions on the mesh.

First of all, the y^+ value here should be y^+\lessapprox 1. But this is not our only requirement. Since we want to resolve all turbulence, we need to have a mesh which is sufficiently small to capture all turbulence scales. This gives us a chicken/egg problem; in order to resolve all scales, we need to have a mesh that allows us to do so, but in order to design a mesh that captures all turbulence scales, we need to know exactly what the length scales involved are, which typically we obtain through simulations. While this seems restrictive at first, we have ways to deal with this, i.e. run a simplified turbulence simulation (e.g. RANS) and estimate the smallest scales or look up mesh requirements in the literature. The latter approach is typically taken, as DNS simulations are not routinely done, only in academic applications, but when they are performed, they are typically for simplified problems for which some body of knowledge exists in the literature that can be used.

To summarise, for DNS simulations, we have y^+\lessapprox 1 and then the remaining mesh needs to be sufficiently small enough to capture the smallest length and time scales in the turbulent flow, which are also known as the Kolmogorov (micro) scales. We can calculate these Kolmogorov (micro) scales and visualise them during post-processing, which gives us an indicator of how well we are resolving our turbulent flow.

Scale-resolved Simulations (SRS)

Scale-resolved or scale-resolving simulations are those where we partially resolve the turbulence. Typically candidates here are Large-Eddy Simulations (LES), Detached-Eddy Simulations (DES) (as well as its derivatives such as Delayed Detached-Eddy Simulations (DDES) and Improved Delayed Detached-Eddy Simulations (IDDES)) and Scale-Adaptive Simulations (SAS).

Depending on the turbulence modelling approach, we have different requirements. But before we look at them, we need to conceptually introduce two additional quantities; the x^+ and z^+ value.

Both the x^+ and z^+ value are entirely analogous to the y^+ value and are defined as

x^+=x\cdot \frac{u_*}{\nu}

z^+=z\cdot\frac{u_*}{\nu}

The x^+ is the non-dimensional coordinate which is parallel to the wall and travels along the turbulent boundary layer (i.e. with the flow). The z^+ is the non-dimensional coordinate direction which is also parallel to the wall and normal to the x^+ direction. This coordinate system is shown below for a wall-bounded turbulent boundary layer flow.

With these definitions at hand, we can look at different grid requirements for the different scale-resolving turbulence modelling techniques.

Large-Eddy Simulations (LES)

If we said that DNS requires the least amount of modelling (i.e. none) but has significant restrictions on the size of the mesh as a result, Large-Eddy Simulations (LES) remove some of the meshing restrictions which allow us to have somewhat coarser grids and thus require less computational resources. For wall-bounded flows, we generally say that we have the following restrictions

50\le x^+\le 150 \\ y^+\le1\\15\le z^+\le 30

So we see that in terms of y^+, we still are restricted to somewhat fine grids. However, we no longer say that we have to resolve all turbulent scales (in fact, LES only requires us to resolve about 80% of all turbulence kinetic energy). This means we have a reduced grid size and faster computations but what we lose in turbulence resolution, we have to model and thus we shift the complexity from resolving all turbulence to modelling at least some part of that. However, if we have a well-written CFD solver which already provides us with functions/classes to compute things like the strain rate tensor, implementing LES would probably not require more than a few lines of code, at least for the simplest LES (subgrid-scale) models.

LES simulations, similar to DNS simulations, are still more of an academic research tool and it is likely that there is published material out there on meshing requirements for a particular problem you may be interested in solving, So if you are dealing with LEs simulations, it is still worthwhile to check if there is published material available for best practices on meshing for that particular application/geometry.

Detached-Eddy Simulations (DES) and Scale-Adaptive Simulations (SAS)

Both Detached-Eddy Simulations (DES) and its derivatives, as well as Scale-Adaptive Simulations (SAS) are conceptually very similar (and, in fact, are only different in the way their limiter is defined). So it makes sense to discuss their meshing requirements together.

In general, we still want to resolve the wall-normal direction well and thus we still have y^+\approx 1, however, both DES and SAS can be considered to be hybrid LES/RANS approaches, i.e. we resolve our flow in the far-field with LES resolution while we use RANS to resolve the nearfield, i.e. the inner part of the turbulent boundary layer. Thus, since we employ RANS close to the wall, we do not have many constraints in terms of x^+ and z^+, only in y^+ and therefore we relax our grid requirements even further, resulting in less computational cost but ultimately requiring more complexity in terms of modelling.

However, for us, the importance here is in terms of grid requirements that we are only concerned about the y^+ value and that it should be close (or less than) unity.

Reynolds-Averaged Navier-Stokes (RANS) Simulations

RANS, then, requires the most modelling, it is the most complex of all of the previously mentioned approaches in terms of modelling but it also comes with the least amount of restriction in terms of grid resolutions which makes it very popular and fast to compute solutions. Coupled with this is that RANS is an inherently steady-state approach (while DNS, LES, DES and SAS are all unsteady approaches), thus we are not just reducing the complexity in terms of meshing but we are also reducing the dimension of the simulation by one which provides additional (and significant) computational savings.

Let us return to the law of the wall, for convenience reproduced below

We saw that there are two distinct regions in any turbulent boundary layer flow; one in which we exhibit a linear relation in our velocity profile (up to y^+\approx 5) and one in which we have a logarithmic relationship for the velocity profile (starting after y^+\approx 30). In between we also saw that there is a transition region.

When generating a computational grid to resolve turbulent boundary layer flows, we get two choices; either resolve the computational grid down to y^+\approx 1 or make it coarse and aim for y^+>30. If we aim for y^+\approx 1, then we just have to use our standard RANS model and we are done, no further modelling is required. If we want to save some additional computational time, however, and aim for y^+>30 as the height of our first cell, then we need to use so-called wall-functions, which account for the fact that the velocity profile does not go to zero in a straight line in our semi-log plot. This means that our turbulence model needs to be able to handle wall-functions and not all RANS turbulence models do. But we need to be aware of that and then either implement wall function modelling into our solver or make use of the appropriate wall functions in the solver we are using.

Thus, to summarise for RANS, we can either have y^+ values around unity for which no wall functions are required and values which are greater than 30, where wall functions are required but where we gain additional computational savings due to a lower grid resolution. We can still get reasonably good results in some applications if we go up to y^+=300 so there is some considerable freedom in generating our grids and we can finetune our grid based on the amount of time and computational resources available, which makes RANS a very popular modelling approach for turbulent flows.