Engineering

How to Calculate First Cell Height for a Target y+ in OpenFOAM

Oktober 29, 2026 / Aldy Dwi Prasetyo

Almost every turbulent CFD case I set up starts with the same question: how thin should the first cell next to the wall be? Get it wrong and the wall shear stress, the pressure drop and the heat transfer you report can all be off, even if the solver converges nicely. The quantity that controls this is y+ (y plus), and the good news is that you can estimate the first cell height before you mesh anything, with a few lines of Python. This tutorial walks through the estimate, shows how to put the number into a snappyHexMesh layer setup, and how to check the actual y+ in OpenFOAM after the run.

What y+ actually measures

y+ is the distance from the wall made dimensionless with the friction velocity and the kinematic viscosity:

y+ = u_tau * y / nu, with u_tau = sqrt(tau_w / rho)

Here y is the distance from the wall to the centre of the first cell, tau_w is the wall shear stress and nu is the kinematic viscosity. In practice y+ tells you which part of the boundary layer your first cell sits in:

  • y+ below about 5: viscous sublayer. The velocity profile is linear and can be resolved directly by the mesh.
  • y+ between about 5 and 30: buffer layer. Neither the linear law nor the log law fits well. Try to keep the first cell out of here when you use classic wall functions.
  • y+ between about 30 and a few hundred: log-law region. This is where standard wall functions are designed to work.

So before meshing you have to make a choice: resolve the wall (target y+ around 1) or model it with wall functions (target y+ around 30 to 100). That choice decides both the first cell height and the boundary conditions you will use.

The problem: y+ depends on a result you do not have yet

The catch is that tau_w comes out of the simulation. To estimate it in advance we use an empirical skin friction correlation. For a turbulent flat plate a common one (from White's Fluid Mechanics) is:

Cf = 0.026 / Re_x^(1/7)

where Re_x = rho * U * L / mu uses a characteristic length L. For internal flows you would use the hydraulic diameter and a pipe friction correlation instead, but the procedure is identical. The steps are:

  1. Compute the Reynolds number from your reference velocity and length.
  2. Estimate the skin friction coefficient Cf.
  3. Compute the wall shear stress: tau_w = 0.5 * Cf * rho * U^2.
  4. Compute the friction velocity: u_tau = sqrt(tau_w / rho).
  5. Solve for the wall distance: y = y+ * mu / (rho * u_tau).
  6. Remember that OpenFOAM evaluates y+ at the cell centre, so the first cell height is roughly 2 * y.

This is an estimate, not an exact answer. Local wall shear on a real geometry can be several times higher near a leading edge or in an accelerating region, and lower near separation. Treat the result as a starting point that you verify after the first run.

A small Python script for the first cell height

Here is the script I use. It only needs Python with NumPy and Matplotlib for the plot.

import math

def first_cell_height(U, L, rho, mu, yplus):
    Re = rho*U*L/mu
    Cf = 0.026/Re**(1/7)          # turbulent flat plate (White)
    tau_w = 0.5*Cf*rho*U**2
    u_tau = math.sqrt(tau_w/rho)
    y = yplus*mu/(rho*u_tau)      # wall-to-centroid distance
    return Re, Cf, tau_w, u_tau, y

def layers(y1, ratio, total):
    # number of geometric layers needed to reach a total thickness
    n = math.log(1 + total*(ratio - 1)/y1)/math.log(ratio)
    return math.ceil(n)

rho, mu = 1.204, 1.825e-5        # air at 20 C, 1 atm
U, L = 10.0, 1.0
for yp in (1, 30, 50):
    Re, Cf, tw, ut, y = first_cell_height(U, L, rho, mu, yp)
    print(f"y+={yp}: Re={Re:.3e} Cf={Cf:.5f} tau_w={tw:.4f} Pa "
          f"u_tau={ut:.4f} m/s first cell={2*y*1e3:.4f} mm")

Example inputs: air at 20 C (density 1.204 kg/m3, dynamic viscosity 1.825e-5 Pa s), free-stream velocity 10 m/s, plate length 1 m. Running the script gives:

y+=1:  Re=6.597e+05 Cf=0.00383 tau_w=0.2308 Pa u_tau=0.4378 m/s first cell=0.0692 mm
y+=30: Re=6.597e+05 Cf=0.00383 tau_w=0.2308 Pa u_tau=0.4378 m/s first cell=2.0772 mm
y+=50: Re=6.597e+05 Cf=0.00383 tau_w=0.2308 Pa u_tau=0.4378 m/s first cell=3.4621 mm

That is a factor of 30 between a wall-resolved mesh and a wall-function mesh, and it shows up directly in your cell count. Note that a Reynolds number of 6.6e5 is close to the transition range on a smooth plate, so assuming a fully turbulent boundary layer here is a deliberate simplification. If transition matters for your case, a plain RANS model with an estimate like this is not the whole story.

The plot below shows how the required first cell height changes with velocity for three y+ targets, using the same air properties and length. Notice the log scale: going from 5 m/s to 50 m/s shrinks the first cell by almost an order of magnitude.

First cell height versus free-stream velocity for y+ targets of 1, 30 and 100, air at 20 C, plate length 1 m

First cell height needed for y+ = 1, 30 and 100 as a function of velocity (air at 20 C, L = 1 m, flat-plate correlation). Plot generated by the author with Python and Matplotlib.

How many prism layers do you need?

The first cell is only the start. The layers should also grow smoothly until they cover a good part of the boundary layer, so the transition to the core mesh is not abrupt. A rough turbulent boundary layer thickness for a flat plate is delta = 0.37 * L / Re^0.2, which for our example gives about 25.4 mm.

With a geometric expansion ratio r and first layer thickness h1, the total thickness of N layers is h1 * (r^N - 1)/(r - 1). Solving for N with r = 1.2 gives:

  • y+ = 1 (first cell 0.069 mm): 24 layers to reach the boundary layer thickness.
  • y+ = 30 (first cell 2.08 mm): 7 layers.

Keep the expansion ratio between about 1.1 and 1.3. Higher values save cells but hurt accuracy in the gradient region, and they make layer addition in snappyHexMesh fail more often.

Putting the numbers into snappyHexMesh

In system/snappyHexMeshDict the layer settings live in addLayersControls. By default snappyHexMesh specifies layer thickness relative to the local surface cell size. For a y+ target you usually want absolute values, so switch relativeSizes off and give the first layer thickness in metres. A wall-function setup for the example above could look like this:

addLayers       true;

addLayersControls
{
    relativeSizes       false;
    layers
    {
        "plate.*"
        {
            nSurfaceLayers  7;
        }
    }
    firstLayerThickness 0.002;
    expansionRatio      1.2;
    minThickness        0.0005;
    nGrow               0;
    featureAngle        130;
    nRelaxIter          5;
    nSmoothSurfaceNormals 1;
    nSmoothNormals      3;
    nSmoothThickness    10;
    maxFaceThicknessRatio 0.5;
    maxThicknessToMedialRatio 0.3;
    minMedialAxisAngle  90;
    nBufferCellsNoExtrude 0;
    nLayerIter          50;
}

Replace "plate.*" with your own wall patch names. After meshing, run checkMesh and open the result in ParaView to confirm the layers were actually added everywhere. snappyHexMesh prints a layer coverage summary at the end of the layer stage, and it is worth reading: sharp edges and tight corners are where layers usually collapse. If you build a structured mesh with blockMesh instead, the same first cell height goes into the simpleGrading or multi-grading of the wall-normal direction.

Matching the wall boundary conditions

The y+ target and the wall treatment must agree. For a kOmegaSST case the usual pairs in the 0/ directory are:

  • Wall functions (y+ around 30 or more): nut as nutkWallFunction, k as kqRWallFunction, omega as omegaWallFunction.
  • Wall resolved (y+ around 1): nut as nutLowReWallFunction (or fixedValue 0), k as fixedValue with a very small value such as 1e-10, and omega still as omegaWallFunction, which switches to its viscous-sublayer form at low y+.

If you cannot keep y+ in one range across the whole surface, which is common on real parts, nutUSpaldingWallFunction is a practical choice because Spalding's law covers the viscous sublayer, buffer and log regions in one expression.

Checking the real y+ after the run

Once the case has converged, compute the actual y+ on the walls. The command depends on which OpenFOAM line you use:

  • OpenFOAM from openfoam.com (for example v2406 or v2412): simpleFoam -postProcess -func yPlus -latestTime
  • OpenFOAM from openfoam.org (versions 11 and 12, modular solvers): foamPostProcess -solver incompressibleFluid -func yPlus -latestTime

Both write a yPlus field you can colour in ParaView, and print the minimum, maximum and average per wall patch. You can also add the function object to system/controlDict so it is written during the run:

functions
{
    yPlus
    {
        type            yPlus;
        libs            ("libfieldFunctionObjects.so");
        writeControl    writeTime;
    }
}

Do not rely on the average alone. Look at the contour plot. A mean of 40 can hide a leading edge at 150 and a separated region at 3. If large areas fall outside your target band, scale the first layer thickness by the ratio of target to actual y+ and remesh. One correction is usually enough.

A few practical notes

  • Use the velocity that the wall actually sees. For a part inside a duct, the local velocity near the wall may be much higher than the inlet value.
  • For heat transfer, wall-resolved meshes (y+ around 1) are usually the safer choice, because thermal wall functions carry additional assumptions.
  • For a first design study, wall functions and a coarser mesh are often good enough. Run a quick case, check the trends, and refine only where it matters.
  • Keep the script with the case files, together with the inputs you used. When someone asks why the first cell is 2 mm, you have the answer in one place.

If you like doing engineering estimates in Python, the same approach works for thermodynamics too: see my refrigeration cycle analysis tutorial with CoolProp.

Wrapping up

Estimating the first cell height takes five minutes and saves a lot of remeshing. Pick the wall treatment first, estimate the wall shear with a correlation, convert the target y+ to a distance, double it for the cell height, and build the layers so they cover the boundary layer smoothly. Then always check the real y+ after the run, because the correlation is only a guess at what the flow will do.

If you want to go deeper into why the near-wall region behaves this way, along with the turbulence models and their equations, my book (see the book overview) Memahami Turbulensi: Konsep, Persamaan, dan Praktik CFD dengan Perangkat Lunak Gratis (in Indonesian) covers the theory and practical CFD with free software. And if this tutorial was useful, follow the blog for more practical CFD and FEA write-ups.

Related posts