Build an incompressible turbulence model¶

The pybFoam.turbulence submodule wraps OpenFOAM’s incompressible RAS / LES factory. This how-to reads the shipped pitzDaily case (backward-facing step, k-epsilon RAS), constructs the turbulence model from Python, and inspects its bound state — nut, k, epsilon, correct(), divDevReff().

Prerequisite: finished Read, modify, and write OpenFOAM dictionaries.

Clone the case¶

from pybFoam import clone_example

case = clone_example("pitzDaily")

Build mesh and read U¶

from pybFoam import (
    Time,
    argList,
    createPhi,
    dictionary,
    fvMesh,
    volVectorField,
)
from pybFoam.meshing import generate_blockmesh

time = Time(argList([str(case), "-case", str(case)]))
generate_blockmesh(time, dictionary.read(str(case / "system" / "blockMeshDict")))
mesh = fvMesh(time)

U = volVectorField.read_field(mesh, "U")
phi = createPhi(U)

Construct the transport and turbulence models¶

singlePhaseTransportModel(U, phi) reads the Newtonian nu from constant/transportProperties. incompressibleTurbulenceModel.New dispatches on simulationType and RASModel from turbulenceProperties to instantiate kEpsilon, kOmegaSST, etc.

from pybFoam.turbulence import (
    incompressibleTurbulenceModel,
    singlePhaseTransportModel,
)

transport = singlePhaseTransportModel(U, phi)
model = incompressibleTurbulenceModel.New(U, phi, transport)

Inspect the model¶

Each accessor returns a tmp<volScalarField>; call () to get the underlying field, then ["internalField"] for the NumPy view.

import numpy as np

nut = np.asarray(model.nut()()["internalField"])
k = np.asarray(model.k()()["internalField"])
eps = np.asarray(model.epsilon()()["internalField"])

print(f"nut       : shape={nut.shape}   min={nut.min():.3e}   max={nut.max():.3e}")
print(f"k         : shape={k.shape}   min={k.min():.3e}   max={k.max():.3e}")
print(f"epsilon   : shape={eps.shape}   min={eps.min():.3e}   max={eps.max():.3e}")

nut       : shape=(12225,)   min=0.000e+00   max=0.000e+00
k         : shape=(12225,)   min=3.750e-01   max=3.750e-01
epsilon   : shape=(12225,)   min=1.486e+01   max=1.486e+01

Update the eddy viscosity¶

correct() is what a PISO/PIMPLE solver calls at the end of every outer iteration. For k-epsilon it solves the k and epsilon transport equations and recomputes nut = Cmu * k^2 / epsilon.

model.correct()

nut_after = np.asarray(model.nut()()["internalField"])
print(f"nut after correct() : min={nut_after.min():.3e}   max={nut_after.max():.3e}")

nut after correct() : min=1.938e-07   max=2.148e-03

divDevReff in a momentum equation¶

In a turbulent solver the laminar diffusion laplacian(nu, U) is replaced by the model’s stress divergence:

UEqn = fvVectorMatrix(fvm.ddt(U) + fvm.div(phi, U) + model.divDevReff(U))

UEqn_term = model.divDevReff(U)
print(f"divDevReff term : {type(UEqn_term).__name__}")

divDevReff term : tmp_fvVectorMatrix