An advanced eddy viscosity type LES model in OpenFOAM®

Compared to many other eddy viscosity type Large Eddy Simulation (LES) models in the literature, the sensor-enhanced Smagorinsky model stands out by an improved sensitivity with respect to small-scale fluctuations of the flow structure including laminar-turbulent transition and correct near-wall scaling for wall-bounded flows. The sub-grid activity sensor is based on explicit test filtering and can be interpreted as a dynamic adaptation of the theoretically determined Smagorinsky constant depending on local flow conditions.

A convincing performance in three test cases relevant to technical applications (laminar-turbulent transition in the Taylor-Green vortex, wall-dominated channel flows, free planar jet flow including passive scalar mixing) and the mathematical details are documented in our recent open-access publication:

  • J. Hasslberger, L. Engelmann, A. Kempf, and M. Klein. Robust dynamic adaptation of the Smagorinsky model based on a sub-grid activity sensor. Physics of Fluids, 33:015117, 2021.

The eddy viscosity according to the sensor-enhanced Smagorinsky model can be expressed as


where the notation and rationale are detailed in the paper mentioned above. 

Meanwhile, promising preliminary results have been obtained from additional testing on multi-physics problems (interface-resolved two-phase flows; reactive flows, i.e. combustion; natural convection). Examples:

  • E. Trautner, J. Hasslberger, T. Trummler, P. Cifani, R. Verstappen, and M. Klein. Towards LES of bubble-laden channel flows: sub-grid scale closures for momentum advection. 13th International ERCOFTAC Symposium on Engineering Turbulence Modelling and Measurements, Rhodes, Greece, Sep 2021.

  • L. Engelmann, J. Hasslberger, E. Inanc, M. Klein, and A. Kempf. A-posteriori assessment of Large-Eddy Simulation subgrid-closures for momentum and scalar fluxes in a turbulent premixed burner experiment. Computers and Fluids, 240:105441, 2022.

Fully consistent with the unstructured mesh handling paradigm, implementation in OpenFOAM® can be realized in the following way. Since only OpenFOAM®-native functionality is used, no parallel-specific implementation is needed.

In order to avoid excessive masking of the explicit sub-grid scale model by numerical dissipation, it is generally recommended to run the LES with central differences for the discretization of convective fluxes, i.e. using scheme Gauss linear in OpenFOAM®. It is likely that the Rhie-Chow-like pressure handling in OpenFOAM® and the diffusive nature of eddy viscosity models already ensure the numerical stability of the LES.

The test filter is realized as a simple box filter. Linear interpolation from cell centres to cell faces and surface area weighting generalize the formulation for arbitrarily shaped unstructured meshes:

tmp filteredField = fvc::surfaceSum(mesh().magSf()*fvc::interpolate(unFilteredField))/fvc::surfaceSum(mesh().magSf());
Using OpenFOAM® syntax, the sensor-enhanced Smagorinsky model reads
tau_smago = -2*nu_turbulent_smago*symm(fvc::grad(U));
E = max(E, dimensionedScalar("minE", E.dimensions(), SMALL));
F_filtered = fvc::surfaceSum(mesh.magSf()*fvc::interpolate(F))
sensorFunction = Foam::pow(mag(F_filtered-F),1.5);

This yields the subgrid-scale stress, which is inserted in the momentum equation:

tmp<fvVectorMatrix> tUEqn
fvm::ddt(U) + fvm::div(phi, U)
+ MRF.DDt(U)
+ turbulence->divDevReff(U)
+ fvc::div(tau_sensor)
fvVectorMatrix& UEqn = tUEqn.ref();
if (pimple.momentumPredictor())
solve(UEqn == -fvc::grad(p));

Through the subgrid-scale viscosity, further closure terms can be determined by means of the gradient flux approach. For example, the implementation of the subgrid-scale thermal diffusivity within the energy equation could look as follows: 

alphat = nu_turbulent_sensor/Prt;
alpha_eff = thermo.alpha() + alphat;

In this example, Prt is the turbulent Prandtl number, which describes the ratio of the subgrid-scale viscosity to the subgrid-scale thermal diffusivity.

To check and facilitate the model’s implementation in your custom-built OpenFOAM® solver, we provide the code based on the pimpleFoam solver and an associated test case that represents the computationally inexpensive Taylor-Green vortex. The solver is based on OpenFOAM-v2006 and the source code has to be compiled with wmake before starting the simulation.

When conducting the simulation of Taylor-Green vortex, the flow field in the triply periodic box of side length 2 π  can be initialized by the setExprFields utilityImposing the boundary conditions (type cyclicand meshing (e.g. 32 x 32 x 32 uniform grid defined in system/blockMeshDict) are straight-forward in this case. 

The temporal evolution of domain-averaged kinetic energy and sensor statistics are monitored as follows. First, the corresponding fields are defined in createFields.H:

//Definition of Kinetic Energy
volScalarField KE
//Domain-Averaged Kinetic Energy
scalar avgKE = KE.weightedAverage(mesh.V()).value();
//Domain Maximum of the Coherent Structure Function
scalar maxSensor = 0.0;
//Domain Maximum of the Sensor Function
scalar maxSensorFunction = 0.0;
//createFields for sensormodel
#include "createFields_LES.H"

Then, in the header of the main solver file pimpleFoam_with_sensor_sgs_model.C, the following output streams are defined:

//define output stream for kinetic energy
fileName outputFileKE("avgKE");
OFstream osKE(outputFileKE);
//define output stream for coherent structure function statistics
fileName outputFileSensor("structureFunction");
OFstream osSensor(outputFileSensor);
//define output stream for sensor function statistics
fileName outputFileSensorFunction("sensorFunction");
OFstream osSensorFunction(outputFileSensorFunction);

Finally, the quantities are calculated and written out within the time loop of the solver:

//On-the-fly kinetic energy calculation and output
KE = 0.5*(mag(U)*mag(U));
avgKE = KE.weightedAverage(mesh.V()).value();
Info << "Kinetic energy domain average = " << avgKE << nl << endl;
osKE << runTime.timeName() << "\t" << avgKE << endl;
//On-the-fly global maximum of the coherent structure function
maxSensor = 0.0;
maxSensorFunction = 0.0;
forAll(F, celli){
maxSensor = max(F[celli],maxSensor);
forAll(sensorFunction, celli)
maxSensorFunction =
reduce(maxSensor, maxOp<scalar>() );
Info << "Coherent structure function F domain maximum = " << maxSensor << nl << endl;
osSensor << runTime.timeName() << "\t" << maxSensor << endl;
reduce(maxSensorFunction, maxOp<scalar>() );
Info << "Sensor function domain maximum = "
<< maxSensorFunction << nl << endl;
osSensorFunction << runTime.timeName() << "\t"
<< maxSensorFunction << endl;

After the successful run of the Taylor-Green vortex simulation, your results for the global kinetic energy should be similar to the profiles in the figure below:


The convective and diffusive fluxes were evaluated by second-order central schemes. Temporal advancement was achieved by the second-order Crank-Nicolson scheme. In order to avoid masking of the results by numerical dissipation due to temporal discretization, it has been ensured that the Courant number is always below the value of 0.1. The sensor constant Csensor = 1/0.16 used here is the same as in the original publication mentioned above but it may be worthwhile to recalibrate it for the specific test filter in OpenFOAM®.