OpenFOAM 14 Part 7: Physical Properties and Models

This part is built as a sequence of short units, each answering a single question. Sections 1–5 fix the material layer (what physicalProperties contains, how viscosity is defined and measured, how v14 handles units, where to find a real example). Sections 6–8 introduce the stress model (laminar, the idea of closure, the Reynolds number as a guide). Sections 9–15 go through the turbulence branches (RAS with kEpsilon and kOmegaSST, the fields they require, inlet estimates, wall functions, coefficients, LES). Sections 16–18 show the neighbouring layers (non-Newtonian rheology, model availability, thermophysical packages). Sections 19–26 are method: dependency audit, two experiments, diagnostics, an independent assessment, solutions and references.

Every worked task is meant to be done in a terminal, with OpenFOAM 14 sourced, on a copy of a case. The answers to the numerical tasks are in section 25; the deep-dive boxes add the theory, the v14 implementation details and the expected outputs that the main text leaves implicit.

A practical suggestion: open two terminals. In the first, keep your copied laminar channel (the one from the Poiseuille case); in the second, a copy of the official pitzDailySteady tutorial. Many sections ask you to compare the two side by side, and the contrast between a laminar case and a RAS case is the main teaching device of this part.

In v14 there is a single executable, foamRun, which loads the solver module named by the solver entry in system/controlDict (for this course solver incompressibleFluid;). It is the module that decides which files of constant/ it needs. incompressibleFluid constructs a viscosity model from physicalProperties and a momentum-transport model from momentumTransport; a compressible module such as fluid instead builds a thermophysical package from physicalProperties and reads quite different entries.

You can see this happen at the beginning of the log: after the mesh is read, the module prints lines announcing the selection of the viscosity model and of the momentum-transport model (for a laminar case the selected laminar stress model, by default the Newtonian one). Reading these few lines is the fastest way to confirm that OpenFOAM interpreted your dictionaries as you intended.

foamDictionary system/controlDict -entry solver -value
foamRun 2>&1 | tee log.foamRun
grep -i "select" log.foamRun

The table above is therefore a map of responsibility: if a run fails or gives an unexpected result, the layer that owns the choice is the first place to look. A wrong viscosity is a constant/ problem; a missing initial field is a 0/ problem; a missing div(phi,k) entry is a system/ problem, even if the cause was a decision taken in constant/momentumTransport.

Instead of opening the file in an editor, query exactly the entries the solver will read. foamDictionary parses the file with the same reader that OpenFOAM uses, so a syntax error shows up here before a run:

foamDictionary constant/physicalProperties -entry viscosityModel -value
foamDictionary constant/physicalProperties -entry nu -value
foamDictionary constant/momentumTransport -entry simulationType -value

Expected answers for the laminar channel: constant, 1e-05 (or whatever your baseline uses) and laminar. The FoamFile/object entry must be physicalProperties: it is the name under which the file is registered, and it should match the file name.

Answer to task 3. If only nu changes, the mesh (constant/polyMesh), the list of fields in 0/ and their boundary conditions, momentumTransport, fvSchemes, fvSolution and controlDict do not change. What changes is the solution: the viscous term of the momentum equation is multiplied by a different coefficient, so velocity profile development, pressure drop and Reynolds number change. If the case is defined by a fixed velocity at the inlet, the flow rate stays the same and the pressure drop responds; if it were driven by a fixed pressure difference, the flow rate would respond instead.

Typical failure. If nu is misspelt or missing, the run stops while constructing the viscosity model with a fatal IO error stating that keyword nu is undefined in dictionary …/constant/physicalProperties. The message names both the keyword and the file: read it before changing anything else.

The v14 release notes state that physical properties should specify units with the value rather than repeating a dimension set. A plain number is read in the base units of the active unit set (SI by default), so nu 1e-05; means 1 × 10⁻⁵ m²/s. Writing the unit explicitly adds a check and allows conversion on reading:

viscosityModel    constant;
nu                1e-05 [m^2/s];

// equivalent, the value is converted to m^2/s on reading:
// nu             10 [cSt];

Both forms give the same internal value: 10 cSt = 10 × 10⁻⁶ m²/s = 1 × 10⁻⁵ m²/s. If you write a unit of the wrong kind (for example [m/s]), the dimension check of the dimensioned property detects the inconsistency. If you write a correct unit with a wrong magnitude, nothing detects it: see section 4.

In older Foundation versions and in many tutorials found online the property was written with an index dimension set, nu [0 2 -1 0 0 0 0] 1e-05;. Read it as “dimensions of m²/s, value 1e-05”; in v14 the preferred style is the bare value or the value followed by a unit.

For a Newtonian fluid the viscous stress is τ = 2μD, where D = ½(∇U + ∇U^T) is the strain-rate tensor (for incompressible flow the volumetric part vanishes). With constant density, the momentum equation can be divided by ρ:

dU/dt + div(U U) = -grad(p/rho) + div( 2 (mu/rho) D )
                 = -grad(p_k)    + div( 2 nu D )

Density disappears as an independent parameter: only ν = μ/ρ appears, and the pressure solved for is the kinematic pressure p_k = p/ρ, in m²/s². This is why incompressibleFluid reads nu and not rho, and why 0/p in such a case is in m²/s², not Pa. To report a pressure difference in pascal, multiply by the density you assume for the fluid: Δp [Pa] = ρ × Δp_k [m²/s²].

Answer to task 3. Changing the density of the real fluid, at the same dynamic viscosity, changes ν. In the incompressible channel that change enters only through nu; there is no separate density entry to edit, and the kinematic pressure field remains independent of ρ until you convert it for reporting.

Useful reference values (approximate, 20 °C, 1 atm): water μ ≈ 1.0 × 10⁻³ Pa·s, ρ ≈ 998 kg/m³, ν ≈ 1.0 × 10⁻⁶ m²/s; air μ ≈ 1.8 × 10⁻⁵ Pa·s, ρ ≈ 1.2 kg/m³, ν ≈ 1.5 × 10⁻⁵ m²/s. Note that air is about 15 times more “viscous” than water in kinematic terms even though its dynamic viscosity is about 55 times lower.

Units of practice. 1 P (poise) = 0.1 Pa·s, so 1 cP = 1 mPa·s = 10⁻³ Pa·s (dynamic). 1 St (stokes) = 10⁻⁴ m²/s, so 1 cSt = 10⁻⁶ m²/s (kinematic). Water at about 20 °C is close to 1 cP and 1 cSt, which is the origin of many confusions: the numbers coincide only because ρ ≈ 1000 kg/m³.

Internally every dimensioned quantity carries a dimension set of seven exponents in the order mass, length, time, temperature, quantity (moles), current, luminous intensity. Kinematic viscosity is [0 2 -1 0 0 0 0]; in v14 the same thing can be written with a name, [kinematicViscosity], or with an expression such as [area/time] or [sqr(length)/time]. The v14 example cases have replaced numerical dimension sets with named ones.

A unit instead says how the number is expressed. Units are written in square brackets after a value ([cm/s], [cm s^-1], [mm], [bar], [cSt], [deg], [%]) and the value is converted to the base units of the active unit set when it is read. The unit sets (SI by default; CGS and USCS are also preconfigured) are defined in the units section of $FOAM_ETC/configDict, and a case can override them in system/configDict.

foamUnits (an application in v14, User Guide §4.7.7) lists dimensions and the units that correspond to them. For pressure the guide shows output of this form:

Dimension [pressure]
+ Dimensions = [mass length^-1 time^-2]
+ Exponents = [1 -1 -2 0 0 0 0]
+ Standard Units = [kg m^-1 s^-2] [Pa]
+ Other Units = [psi] [atm] [kPa] [bar] [MPa]

Running foamUnits kinematicViscosity gives the analogous report: exponents [0 2 -1 0 0 0 0], standard unit m²/s and any other units defined for that dimension in your configuration.

What is checked and what is not. Dimensioned properties and fields compare the supplied unit with the dimensions hard-coded in the solver; plain scalar parameters (a relaxation factor, a tolerance) are not dimension-checked. When two quantities with different dimensions are added or equated, OpenFOAM stops with a fatal error stating that the dimensions differ for the operation (the message shows the two dimension sets). A wrong magnitude — 1e-5 typed as 1e-6, a value in cSt written without the unit — passes every check: compare with an independent hand estimate, as in section 3.

The complete sequence, working on a copy (never inside $FOAM_TUTORIALS):

mkdir -p "$FOAM_RUN"
cp -r "$FOAM_TUTORIALS/incompressibleFluid/pitzDailySteady" "$FOAM_RUN"/
cd "$FOAM_RUN/pitzDailySteady"
ls 0 constant system
foamDictionary system/controlDict -entry solver -value
foamDictionary constant/physicalProperties
foamDictionary constant/momentumTransport

According to the v14 User Guide (§2.1), the tutorial runs the incompressibleFluid module through foamRun; physicalProperties contains viscosityModel constant; and nu 1e-05;; momentumTransport contains simulationType RAS; with model kEpsilon; and turbulence on; (the guide also shows a viscosityModel Newtonian; entry in the RAS settings, the laminar part of the stress, and a comment listing other models tested on this case). 0/ contains U, p, k, epsilon, nut and also omega, so that kOmegaSST can be selected without creating a new file.

The geometry is the Pitz–Daily backward-facing step, with a step height of 25.4 mm and an inlet velocity of 10 m/s. With ν = 1 × 10⁻⁵ m²/s the Reynolds number based on step height is 10 × 0.0254 / 10⁻⁵ ≈ 25 400: clearly turbulent, which is why the tutorial uses RAS and not laminar.

Answer to task 4 (one possible formulation). Separating the material from the stress model means a viscosity can be changed without touching the closure, and a closure can be changed without retyping the material data. It also documents the case: anyone can see at once which part is a property of the fluid and which part is a modelling choice.

Much material online refers to older Foundation releases or to the ESI fork (openfoam.com), whose files have the same purpose but different names and keywords. The mapping below is a guide for reading such examples; always rewrite them in the v14 form and verify with the installed tools.

The v14 release notes also say that the Coeffs suffix has been removed from the names of model coefficient sub-dictionaries (with backward compatibility). If a copied entry is not recognised, OpenFOAM often does not fail: an unknown keyword in a dictionary is simply ignored, and the model runs with its defaults. That is the most dangerous outcome, because nothing warns you. Use #print (section 14) to see what the model actually read.

simulationType selects one of three branches: laminar, RAS or LES. Within the laminar branch there is again a model selection, made by the optional laminar { model …; } sub-dictionary. When it is absent, the default Newtonian laminar model is used (in the Foundation source code this default class has traditionally been named Stokes, after the linear viscous-stress law; the name printed in the log and by foamToC in your build is authoritative). The other laminar models are the generalised-Newtonian and visco-elastic ones of section 16.

The equations solved by incompressibleFluid in this configuration are continuity and the momentum equation with the stress 2νD and ν taken from physicalProperties. No k, epsilon, omega or nut fields are needed, and if they are present in 0/ they are simply not read.

To list the laminar models available in your build:

foamToC -table laminarincompressibleMomentumTransportModel

Worked example for the “why it matters” box. For fully developed flow between parallel plates spaced H, with bulk velocity U_b, the kinematic pressure gradient is −dp_k/dx = 12 ν U_b / H². With ν = 1 × 10⁻⁵ m²/s, U_b = 0.1 m/s, H = 0.01 m and a length of 0.1 m: Δp_k = 12 × 10⁻⁵ × 0.1 / 10⁻⁴ × 0.1 = 0.012 m²/s² (12 Pa for a fluid with ρ = 1000 kg/m³). The drop is entirely viscous and proportional to ν — but only in the fully developed part of the channel (see section 20).

Split the velocity into a mean and a fluctuation and average the momentum equation. The non-linear convective term produces a new unknown, the Reynolds stress tensor ⟨u′u′⟩: six additional unknown components with no additional equation. This is the closure problem.

U = <U> + u'                         (Reynolds decomposition)

d<U>/dt + div(<U><U>) = -grad(<p>) + div[ 2 nu D(<U>) ] - div( <u'u'> )

-<u'u'> = 2 nut D(<U>) - (2/3) k I   (Boussinesq hypothesis)

nuEff = nu + nut

The Boussinesq hypothesis models the Reynolds stresses as an extra viscous stress with a turbulent (eddy) viscosity ν_t, which is a property of the flow, not of the fluid: it varies in space and time, can be hundreds of times larger than ν in a turbulent core and goes to zero at a wall. The momentum equation then has the same form as the laminar one with ν replaced by ν_eff = ν + ν_t. The RAS model’s job is to compute ν_t, usually from one or two additional transported scalars (k and ε, k and ω, or ν̃).

In OpenFOAM this is visible in the field names: nut is ν_t, stored as a field with its own boundary conditions; the effective viscosity nuEff appears in the names of some terms of fvSchemes. LES uses the same idea at a different scale: the filtered equations contain a subgrid stress modelled with a subgrid viscosity, and the large eddies are resolved in time rather than averaged away.

A Reynolds number is meaningful only together with its definition. For internal flows the usual choice is the bulk (mean) velocity and the hydraulic diameter D_h = 4A/P (A cross-section area, P wetted perimeter). For a circular pipe D_h = D; for two infinite parallel plates spaced H, D_h = 2H; for the 10 mm × 5 mm duct of Part 2, D_h = 4 × 50 / 30 ≈ 6.67 mm. Using the channel height instead of D_h changes Re by a factor of two for plates: the Re = 2000 based on H used here would be Re = 4000 based on D_h.

The classical reference values are therefore not interchangeable: about 2300 for a pipe (bulk velocity, diameter) as the usual lower bound of transition; for plane channel flow, linear stability theory gives about 5772 based on centreline velocity and half-height, while in experiments transition can occur at considerably lower values when disturbances are present. Inlet disturbances, roughness and the shape of the inlet profile all move the actual threshold.

Fully developed is a different question. Re says nothing about the distance needed for the profile to stop changing. For laminar internal flow the entrance length grows with Re; a common estimate for pipes is L_e ≈ 0.05 Re D_h. At Re = 400 and D_h = 8 mm this already gives about 0.16 m: a short channel may never be fully developed. Check it in the solution by comparing velocity profiles at two stations, not by Re.

Answer to task 3. At low Re, with a physically laminar flow, a laminar model is not a simplification but the correct physics: there are no turbulent fluctuations to model, the solution can be compared with an analytical one (Poiseuille), and adding a RAS model would introduce a turbulent viscosity that does not exist, extra fields and extra uncertainty.

The standard k–ε model transports the turbulent kinetic energy k (m²/s²) and its dissipation rate ε (m²/s³) and builds the eddy viscosity from them:

nut = Cmu * k^2 / epsilon

Dk/Dt       = div[ (nu + nut/sigmak)   grad(k) ]       + G - epsilon
Depsilon/Dt = div[ (nu + nut/sigmaEps) grad(epsilon) ] + C1*G*epsilon/k - C2*epsilon^2/k

G = nut * 2 |symm(grad(U))|^2

standard constants: Cmu = 0.09, C1 = 1.44, C2 = 1.92, sigmak = 1.0, sigmaEps = 1.3

G is the production of k by the mean shear. The constants are those of the standard model; in OpenFOAM they are defaults in the implementation and are printed with #print (section 14).

foamToC -table RASincompressibleMomentumTransportModel prints the names registered in that run-time selection table. According to the v14 User Guide the incompressible RAS list includes LRR, LamBremhorstKE, LaunderSharmaKE, LienCubicKE, LienLeschziner, RNGkEpsilon, SSG, ShihQuadraticKE, SpalartAllmaras, kEpsilon, kEpsilonLopesdaCosta, kOmega, kOmega2006, kOmegaSST, kOmegaSSTLM, kOmegaSSTSAS, kkLOmega, qZeta, realizableKE and v2f. Your installation is the final authority. foamSearch "$FOAM_TUTORIALS" momentumTransport RAS/model instead scans all momentumTransport files of the tutorials and prints the distinct values of the RAS/model entry: it tells you which models are actually used in official examples.

turbulence off. With turbulence off; the model equations are not solved, so k, ε and ν_t keep the values they have. It is a tool for special procedures (for example initialising a flow field), not a way to obtain a laminar simulation: for that, use simulationType laminar.

Each new file must have the right dimensions and a boundary condition for every patch of the mesh:

In pitzDailySteady, as described by the v14 guide: k and epsilon have internal and inlet values 0.375 and 14.855; on the walls k uses kqRWallFunction, epsilon uses epsilonWallFunction and nut uses nutkWallFunction; nut is calculated on inlet and outlet; frontAndBack is empty for all fields because the case is two-dimensional. Read the outlet conditions of k and epsilon in your copy (an outflow condition that switches to a fixed value if backflow occurs is the usual pattern).

The v14 release notes also report that the inlet turbulence boundary conditions have been simplified into turbulentKineticEnergy, turbulentEpsilon and turbulentOmega, which, as their names and the conditions they replace suggest, derive inlet values from turbulence-intensity and length-scale (or mixing-length) assumptions. Check their actual parameters with foamInfo turbulentKineticEnergy before using them.

A k–ε run additionally needs entries such as these (illustrative pattern; compare with your copy of pitzDailySteady):

// system/fvSchemes - entries a k-epsilon run additionally needs (illustrative)
divSchemes
{
    // ... existing entries of the case ...
    div(phi,k)          bounded Gauss limitedLinear 1;
    div(phi,epsilon)    bounded Gauss limitedLinear 1;
}

// system/fvSolution - the turbulence fields need a linear solver
solvers
{
    // ... existing p and U entries ...
    "(k|epsilon)"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-05;
        relTol          0.1;
    }
}

In a steady SIMPLE case the turbulence fields also usually have relaxation factors and residual controls in fvSolution: see Part 6.

k–ω SST transports k and the specific dissipation rate ω (1/s), related to ε by ω = ε/(C_μ k). It blends a k–ω formulation near walls with a k–ε-like behaviour away from them and limits the eddy viscosity in strong adverse pressure gradients (the “shear-stress transport” limiter). Its wall treatment uses omegaWallFunction on omega. Spalart–Allmaras is a one-equation model for a modified viscosity ν̃ (nuTilda), from which ν_t is computed; its boundary conditions are of a different kind (no k, no ε/ω).

Finding the examples. foamSearch lists the distinct values of an entry and how often they occur, which confirms that kOmegaSST is used; to get the actual tutorial paths, search the files, or ask foamInfo, which also lists example cases:

foamToC -help
foamToC -table RAS
foamToC -scalarBCs | grep -i wallfunction
grep -rl --include=momentumTransport kOmegaSST "$FOAM_TUTORIALS"

Answer to task 2. epsilon is replaced by omega; k and nut remain but their boundary conditions must be checked again. Answer to task 3: the fields in 0/ and their boundary conditions (especially walls and inlet), fvSchemes (new div(phi,omega)), fvSolution (solver, relaxation, residual control for omega), the near-wall mesh (y+) and the observables you will compare.

A concrete warning from the official tutorial. The v14 guide reports that on pitzDailySteady the kEpsilon run converges in 285 iterations and realizableKE in 255, while kOmegaSST with the same numerics does not converge and runs to the end time (2000); it converges in 285 iterations only after changing the gradient scheme of U to cellLimited. Same mesh, same case, different model: the numerical settings had to change too.

Where the k = 3/2 (U I)² formula comes from: with isotropic fluctuations u′ = v′ = w′ = U I, k = ½(u′² + v′² + w′²) = 3/2 (U I)². The other quantities follow from k and a length scale l:

epsilon = Cmu^0.75 * k^1.5 / l
omega   = k^0.5 / (Cmu^0.25 * l)  =  epsilon / (Cmu * k)
nut     = Cmu * k^2 / epsilon
(Cmu = 0.09)

Continuing the example above (U = 2.0 m/s, I = 0.05, k = 0.015 m²/s²) with an assumed l = 1 mm: ε ≈ 0.302 m²/s³, ω ≈ 224 1/s, ν_t ≈ 6.7 × 10⁻⁵ m²/s, that is about 67 times the viscosity of water. A frequently quoted rule for fully developed internal flow takes l ≈ 0.07 D_h; it is an empirical rule of thumb, to be justified like I.

Reverse-engineering the official tutorial. pitzDailySteady has U = 10 m/s and k = 0.375 m²/s²: 1.5 (10 I)² = 0.375 gives I = 0.05 (5%). Its ε = 14.855 m²/s³ gives l = 0.09^0.75 × 0.375^1.5 / 14.855 ≈ 0.00254 m, i.e. one tenth of the 25.4 mm step height; and ω = 14.855/(0.09 × 0.375) ≈ 440.2 1/s, exactly the value the guide reports for 0/omega. The tutorial values are therefore consistent with I = 5% and l = 2.54 mm.

Answer to task 2. I describes the fluctuations of the incoming flow: it depends on what is upstream (a long pipe, a fan, a honeycomb, a contraction) and is rarely measured. Guessing it fixes k at the inlet, and through l also ε or ω and so ν_t; the downstream development, transition and separation can depend on it. Answer to task 3: k is a property of the flow (energy of velocity fluctuations, m²/s²), ν is a property of the fluid (m²/s): different dimensions and different physical meaning, so neither can be derived from the other.

Near a smooth wall the mean velocity, scaled by the friction velocity u_τ = √(τ_w/ρ), follows a universal profile in the variable y⁺: a viscous sublayer (u⁺ ≈ y⁺, roughly y⁺ < 5), a buffer layer, and a logarithmic layer (u⁺ = (1/κ) ln(E y⁺), with κ ≈ 0.41 and E ≈ 9.8, roughly from y⁺ ≈ 30 upwards). A wall function uses this law to supply the wall shear stress and the near-wall values of the turbulence quantities when the first cell is too coarse to resolve the gradient.

In OpenFOAM the wall function is a boundary condition: nutkWallFunction computes ν_t at the wall from k; epsilonWallFunction (or omegaWallFunction) sets the near-wall dissipation; kqRWallFunction is essentially a zero-gradient condition for k. The v14 guide lists the ν_t wall functions nutWallFunction, nutRoughWallFunction, nutUSpaldingWallFunction, nutkWallFunction and nutkAtmWallFunction; the constants E, κ and C_μ can be overridden per boundary condition. Check the exact names in your build with foamToC -scalarBCs | grep -i wallfunction.

Worked example. If τ_w/ρ = 0.25 m²/s², then u_τ = 0.5 m/s; with ν = 1 × 10⁻⁵ m²/s and the first cell centre at y = 0.5 mm, y⁺ = 0.5 × 5 × 10⁻⁴ / 10⁻⁵ = 25. To compute y⁺ on a result, the post-processing function must be able to construct the momentum-transport model, so run it through the solver module:

foamPostProcess -solver incompressibleFluid -func yPlus -latestTime
foamPostProcess -solver incompressibleFluid -func wallShearStress -latestTime

Answer to task 3. A 0/nut written for kEpsilon assumes that k exists and that the wall treatment is consistent with the first-cell y⁺ of that mesh; after switching to SpalartAllmaras (no k) or to a low-Re treatment requiring y⁺ ≈ 1, the same entries are wrong or inconsistent even if the case starts. The patch names in 0/nut must also match constant/polyMesh/boundary exactly (task 1), otherwise the run stops reporting that no patchField entry exists for that patch.

#print is a v14 dictionary directive: placed at the end of a dictionary, it makes OpenFOAM print the dictionary as it is used, including the default entries you did not write. According to the v14 release notes it replaces the writeOptionalEntries function object; in older releases the per-model printCoeffs switch served a similar purpose for model coefficients. In the RAS dictionary of a kEpsilon case it reveals the default coefficients (for example Cmu 0.09;); the same directive inside a PIMPLE dictionary reveals defaults such as nCorrectors 1;. Run the case (a single iteration is enough) and read the printed dictionary at the start of the log. The guide also mentions expanding the file with foamDictionary as an alternative.

If an override is justified, write it in the RAS dictionary or in a sub-dictionary named after the model (in v14 without the old Coeffs suffix):

RAS
{
    model           kEpsilon;
    turbulence      on;

    kEpsilon                // optional, only for a documented override
    {
        Cmu         0.09;
    }
}

Then run again with #print and verify that the printed value is yours: an override with a misspelt keyword is silently ignored. Coefficients are not properties of the fluid: changing C_μ recalibrates the model, and every result obtained with it must say so.

LES applies a spatial filter of width Δ to the equations. Motions larger than Δ are computed; the effect of the smaller ones appears as a subgrid stress, modelled for example with a subgrid viscosity. In the Smagorinsky model ν_sgs = (C_s Δ)² |S|, with |S| the magnitude of the resolved strain rate; WALE, kEqn, dynamicKEqn and dynamicLagrangian are other choices listed by the v14 guide, together with the hybrid RAS/LES models SpalartAllmarasDES, SpalartAllmarasDDES, SpalartAllmarasIDDES and kOmegaSSTDES, and DeardorffDiffStress.

In OpenFOAM the filter width is itself a run-time-selected delta model, named by the delta entry; the guide gives cubeRootVol (Δ from the cube root of the cell volume) as an example, with its coefficients in a sub-dictionary named after it. The LES dictionary therefore needs model, turbulence, delta and, as needed, model- and delta-named sub-dictionaries. Find working combinations with foamSearch "$FOAM_TUTORIALS" momentumTransport LES/model and foamSearch "$FOAM_TUTORIALS" momentumTransport LES/delta, and list the models with foamToC -table LESincompressibleMomentumTransportModel.

Why it is demanding: LES is always three-dimensional and unsteady; the time step must keep the Courant number low enough to follow the eddies; inflow boundaries must supply realistic fluctuations; results are statistics, accumulated over many flow-through times (for example with a field-averaging function); near walls the energy-carrying eddies become very small, so a wall-resolved LES needs a mesh far finer than a RAS mesh. That is why the course keeps the laminar channel as its learning object.

A generalised-Newtonian model keeps the form τ = 2ν(γ̇)D but lets the viscosity depend on the local strain rate γ̇. The v14 guide (§8.3) describes these laws; the zero-strain-rate viscosity ν₀ is the nu of physicalProperties:

strain rate:     gammaDot = sqrt(2) * |symm(grad(U))|

BirdCarreau:     nu = nuInf + (nu0 - nuInf) * [1 + (k*gammaDot)^a]^((n - 1)/a)     (a = 2 by default)
CrossPowerLaw:   nu = nuInf + (nu0 - nuInf) / [1 + (m*gammaDot)^n]
powerLaw:        nu = k * gammaDot^(n - 1),          limited to nuMin <= nu <= nuMax
HerschelBulkley: nu = min( nu0, (tau0 + k*gammaDot^n) / gammaDot )

Notice that k and n appear in several models with different meanings and dimensions (a time constant in BirdCarreau, a consistency index in powerLaw and HerschelBulkley): this is exactly the trap warned about above. A complete illustrative branch, with arbitrary numbers that do not describe a real material:

// constant/physicalProperties: nu is the zero-strain-rate viscosity nu0
viscosityModel  constant;
nu              1e-03;

// constant/momentumTransport (illustrative values, not a material)
simulationType  laminar;

laminar
{
    model           generalisedNewtonian;
    viscosityModel  BirdCarreau;

    nuInf           1e-05;
    k               1;
    n               0.5;
}

Beyond generalised-Newtonian fluids, the laminar branch also contains visco-elastic models (Maxwell, Giesekus, PTT, with parameters such as nuM and lambda) and the lambdaThixotropic model; the v14 release notes add multi-mode visco-elastic models with several relaxation times. The available viscosity laws are listed by foamToC -table generalisedNewtonianViscosityModel.

The same viscosityModel entry also exists in the RAS and LES dictionaries, where it describes the laminar (molecular) part of the stress: it is Newtonian by default, which is why the v14 guide shows viscosityModel Newtonian; in the pitzDailySteady RAS settings and notes that it could be omitted. Rheology and turbulence closure are therefore two independent choices that can be combined.

OpenFOAM builds models by name at run time: each family has a run-time selection table that maps a word (kEpsilon) to a class. The table name encodes the family — RAS + incompressible + MomentumTransportModel — which is why a compressible model is in a different table (the guide shows, for example, buoyantKEpsilon only among the compressible RAS models). foamToC -help lists the other options of the tool; foamToC -table RAS is the form the guide uses to search for tables whose name contains RAS.

If you write a name that is not in the table, the run stops with a fatal error stating that the type is unknown, followed by the list of valid types (the exact wording depends on the build): the error message is itself a foamToC listing, filtered by the module actually running. It is the definitive answer to “does this module support this model?”.

foamInfo kEpsilon prints the description from the source header, the location of the source files and a list of tutorials that use the model: the quickest answer to task 3. foamSearch works on the files of a directory tree and is the right tool to see which values are used in practice.

In a compressible or thermal module (for example fluid), physicalProperties describes a whole thermophysical package. Each entry of thermoType selects one ingredient; the options below are those listed by the v14 guide (§8.1):

The package then needs the data for the chosen ingredients, grouped in a sub-dictionary: molecular weight, heat capacity and heat of formation, dynamic viscosity and Prandtl number for const transport. An illustrative structure for air:

mixture
{
    specie
    {
        molWeight       28.96;
    }
    thermodynamics
    {
        Cp              1004.5;
        Hf              0;
    }
    transport
    {
        mu              1.8e-05;
        Pr              0.7;
    }
}

Note that here the transport property is the dynamic viscosity mu (Pa·s), not nu: in a compressible module density is a field, so the kinematic shortcut of section 3 is no longer available. With transport sutherland, μ also becomes a function of temperature. A fresh compressible case also needs a temperature field 0/T, a pressure in pascal and thermal boundary conditions.

The fastest way to see what a physics change touches is to compare a case that already has it. Copy a laminar case and a RAS tutorial and diff their structures:

ls 0 constant system
diff <(ls ../laminarCase/0) <(ls 0)
foamDictionary system/fvSchemes -entry divSchemes
foamDictionary system/fvSolution -entry solvers -keywords

-keywords prints only the keywords of a dictionary, which shows at a glance which fields have a linear solver. Answer the four questions of the golden rule in writing for each row of the table. For laminar → kEpsilon, for example: (1) two transport equations for k and ε are added and the momentum equation uses ν + ν_t; (2) the unknowns k, ε and the derived field ν_t are added; (3) initial and inlet values of k and ε, wall functions on every wall, outlet conditions; (4) div(phi,k), div(phi,epsilon), solvers and, for SIMPLE, relaxation and residual controls for the new fields.

For “Newtonian → BirdCarreau”, the extra word “output/validation” in the table is important: with a shear-thinning fluid the analytical Poiseuille solution no longer applies, so the comparison you used to validate the laminar channel must be replaced by an appropriate reference.

cp -a copies the whole directory preserving timestamps and permissions. If channelBaseline already contains result time directories, remove them in each copy (foamListTimes -rm, run inside the copy) so that all three runs start from the same 0/. Then set the values with foamDictionary and verify:

foamDictionary channel_nu_half/constant/physicalProperties   -entry nu -set 5e-06
foamDictionary channel_nu_double/constant/physicalProperties -entry nu -set 2e-05

for c in channel_nu_base channel_nu_half channel_nu_double
do
    printf "%-20s nu = %s\n" "$c" "$(foamDictionary $c/constant/physicalProperties -entry nu -value)"
done

diff -r channel_nu_base/constant channel_nu_half/constant

The diff -r output must show only the nu line. Anything else is an uncontrolled change and invalidates the comparison.

What theory predicts, and when. In the fully developed region of a laminar channel at fixed flow rate, the kinematic pressure gradient is proportional to ν (12 ν U_b/H² for plates, the analogous Hagen–Poiseuille law for a pipe), so HALF and DOUBLE should give about half and twice the BASE gradient there. The total inlet-to-outlet difference also contains the developing region, where the excess pressure loss depends on Re: with ν halved, Re doubles and the entrance length grows, so the total ratio can deviate from 2. If the inlet were a fixed pressure difference instead of a fixed velocity, the flow rate would change instead. These are the reasons for the interpretation warning.

First make the run itself traceable: foamRun | tee captures only standard output, so add 2>&1 to keep error messages in the log too, and check how far the run went:

foamListTimes -rm                # remove old results, keep 0/
foamRun 2>&1 | tee log.foamRun
grep -c "^Time =" log.foamRun    # how many iterations/time steps were run
tail -n 30 log.foamRun

Then measure on patches, not in arbitrary cells. The packaged post-processing functions of v14 include patchAverage (average of fields on a patch), patchDifference (difference of averages between two patches), patchFlowRate and patchFlux (flux through a patch). The flow rate through inlet and outlet must be equal and opposite in sign: that is your flux check. For the pressure, the area-averaged p on inlet minus that on outlet is a well-defined Δp_k (multiply by ρ for pascal).

cd channel_nu_base
foamListTimes                    # which time directories exist?
foamPostProcess -func "patchFlowRate(name=inlet)"  -latestTime
foamPostProcess -func "patchFlowRate(name=outlet)" -latestTime
foamGet patchAverage             # copy the template into system/ and read its arguments

The v14 guide writes the patch argument as name= on the command line and as patch= in an #includeFunc example, so read the argument list in the header of the template that foamGet copies into system/ before running patchAverage or patchDifference. Use the same time (-latestTime of a converged run, or the same end time) for all three cases, and report the commands themselves: this is the “actual measurement procedure” of task 3.

Plausibility (task 4). Re halves when ν doubles (task 1). Compare the ratios Δp_half/Δp_base and Δp_double/Δp_base with 0.5 and 2: a ratio near those values with the flux check satisfied is physically plausible for a laminar channel at fixed flow rate; a large deviation needs an explanation (developing flow, unconverged run, a mesh too coarse for the thinner boundary layers of the lower-ν run).

Copy the tutorial (see the commands in the deep dive of section 5) and fill each item with an actual entry read from the copy, not from memory. A dependency diagram can be as simple as this text sketch; your version must contain the values and names you actually found:

physicalProperties : viscosityModel constant; nu 1e-05
        |
        v
momentumTransport  : simulationType RAS; model kEpsilon
        |
        +--> 0/k, 0/epsilon, 0/nut   (+ 0/U, 0/p)
        |        |
        |        +--> wall patches: kqRWallFunction, epsilonWallFunction, nutkWallFunction
        |        +--> inlet: values of k and epsilon (I, l assumptions)
        |
        +--> fvSchemes : div(phi,k), div(phi,epsilon)
        +--> fvSolution: solvers for k, epsilon; residualControl

Expected findings, according to the v14 guide: viscosityModel constant, nu 1e-05; simulationType RAS, model kEpsilon, turbulence on; fields k, epsilon, nut (and omega prepared for kOmegaSST); kqRWallFunction, epsilonWallFunction, nutkWallFunction on upperWall and lowerWall; inlet k = 0.375 m²/s² and ε = 14.855 m²/s³ (I = 5% and l ≈ 2.54 mm, as derived in section 12); SIMPLE with residualControl for p, U and the turbulence fields. The linear solvers and divergence schemes must be read from your copy.

Why “laminar → RAS” is not one word. The laminar case has no k, ε, ν_t fields, no wall functions, no inlet turbulence data, no schemes or solvers for these fields; and its mesh was designed for a laminar profile, not for a wall-function y⁺. Changing only simulationType would stop at the first missing file — or, worse, run with inconsistent data.

The first lines after --> FOAM FATAL ERROR or --> FOAM FATAL IO ERROR almost always name the file and the keyword. Typical patterns, one per row of the table (the exact wording may vary slightly between builds):

--> FOAM FATAL ERROR: cannot find file ".../0/k"

--> FOAM FATAL IO ERROR: keyword div(phi,k) is undefined in dictionary ".../system/fvSchemes/divSchemes"

Unknown ... type kEpsilonn
Valid ... types are: ...   (list of the names in the selection table)

Cannot find patchField entry for walls

Read the error bottom-up in context: the lines just above it in the log tell you what OpenFOAM was constructing (the momentum-transport model, a field, a matrix) when it failed. Note also the file path and, for IO errors, the line number. Only then open the indicated file. For “poor convergence” there is no single message: plot the residuals (foamMonitor on the residuals file, if you have one, or grep on the log), look at where in the domain the solution misbehaves, and change one thing at a time — the boundary conditions and the mesh first, the numerical controls last.

Keep a short log of each attempt (what you changed, what happened). With the method of this section, the record is the evidence; with random fixes, it is lost.

To list candidate cases: ls "$FOAM_TUTORIALS/incompressibleFluid", then check each one with foamDictionary <case>/constant/momentumTransport -entry simulationType -value. A skeleton you can fill in, one line per deliverable item:

Case identification : $FOAM_TUTORIALS/incompressibleFluid/<case>, OpenFOAM 14, solver incompressibleFluid
Material            : constant/physicalProperties -> viscosityModel ..., nu ... [m^2/s]
Momentum model      : simulationType ...; model ...
Unknowns            : 0/U 0/p 0/... (reason for each)
Near-wall           : patches ... -> nut ..., k ..., epsilon/omega ...
Numerics            : fvSchemes div(phi,...) ; fvSolution solvers "(...)"
Evidence            : User Guide v14 section ..., tutorial path ...

For the “Evidence” line, an official reference page is a section of the v14 User Guide (for example §8.2 for the turbulence model, §2.1 if you chose pitzDailySteady), and an installed example is a tutorial path with the entry you cite. For the stretch task, the natural pair is kEpsilon vs kOmegaSST on pitzDailySteady: epsilon vs omega, epsilonWallFunction vs omegaWallFunction, the extra scheme and solver entries, and the documented difference in convergence behaviour (section 11) — which is a numerical observation, not a statement about which model is more accurate.

ν. 1.8 × 10⁻⁵ / 1.2 = 1.5 × 10⁻⁵ m²/s — the kinematic viscosity of air at room temperature, the same value used by many OpenFOAM tutorials. Re. 0.05 × 0.008 / 10⁻⁶ = 4 × 10⁻⁴ / 10⁻⁶ = 400. k. U I = 1.5 × 0.03 = 0.045 m/s; 0.045² = 0.002025; × 1.5 = 0.0030375 m²/s².

Section 3, task 3: see the deep dive “the Newtonian law and why the incompressible solver needs no density”. Section 5: foamVersion must report 14; the tutorial value is nu 1e-05; simulationType RAS. Section 11: epsilon is replaced by omega. Section 13: the patch names must be identical in 0/nut and constant/polyMesh/boundary; wall functions are in 0/nut, 0/k and 0/epsilon (or 0/omega). Section 17: the outputs of the foamToC tables and of foamInfo depend on your build; at least kEpsilon and kOmegaSST are in the incompressible RAS table, and pitzDailySteady demonstrates kEpsilon.

Experiment A, a numerical cross-check. If the BASE channel had H = 10 mm, U_b = 0.1 m/s and ν = 1 × 10⁻⁵ m²/s, then Re_H = 100; HALF gives Re_H = 200 and DOUBLE Re_H = 50. The fully developed pressure gradients would be 0.12, 0.06 and 0.24 m²/s² per metre. Use your own case’s values: these numbers only show the expected pattern.

For any unfamiliar incompressible case, these commands answer the five questions above in less than a minute, and are a good habit to carry into Part 8:

foamDictionary system/controlDict -entry solver -value          # which module
foamDictionary constant/physicalProperties                        # (1) material
foamDictionary constant/momentumTransport                         # (2) model
ls 0                                                              # (3) fields
foamDictionary 0/nut -entry boundaryField                         # (4) wall/inlet treatment
foamInfo kEpsilon                                                 # (5) documentation and examples

Replace 0/nut and kEpsilon with the fields and model you actually found. The point is not the commands themselves but the order of the questions: material, model, unknowns, boundaries, evidence.