Technical University of Denmark
waveTrainTutorial.tar.gz (1.49 MB)

Tutorial for "multiphaseStabilizedTurbulence" in OpenFOAM® - v1912

Download (1.49 MB)
online resource
posted on 2020-04-27, 07:04 authored by Bjarke Eltard LarsenBjarke Eltard Larsen, David R. FuhrmanDavid R. Fuhrman


This data set contains a tutorial demonstrating the usage of the multiphaseStabilizedTurbulence method included in the official OpenFOAM® - v1912 release (OpenFOAM, 2019):

The new method fundamentally solves the problem of over-production (exponential growth) of turbulent kinetic energy (hence eddy viscosity) in the nearly potential flow region beneath non-breaking surface waves using Reynolds-averaged Navier-Stokes (RANS) models. This is accomplished by modifying the eddy viscosity, as described in Larsen & Fuhrman (2018, referred to as LF18 below).
Additionally, the method also includes a buoyancy production term (so-called buoyancy modification, suggested in this context by Devolder et al. 2017) in the k-equation, which causes a sink in turbulence specifically near the air-water interface (this term is also included in the model of LF18). It is emphasized that the problem of turbulence over-production will occur throughout the entirety of the nearly-potential flow region beneath non-breaking surface waves using standard turbulence closure models i.e. this problem is not somehow limited to the near-surface region, as described in LF18 and further underlined in Fuhrman & Larsen (2020). Achieving formal stability hence requires modification of the eddy viscosity, as described above and in LF18

General usage

The method can be used together with the standard two-equation RANS turbulence models using the same setups as would normally be used. It is utilized by adding the following to constant/fvOptions:

type multiphaseStabilizedTurbulence;
active yes;

// Optional coefficients
lambda2 0.1; // A value of 0 sets the nut correction to 0
Cmu 0.09; // from k-epsilon model
C 1.51; // model coefficient from k-omega model
alpha 1.36; // 1/Prt

// Notes:
//lambda2 defines the threshold of the rotation-to-strain rate ratio which will be formally stable (see e.g. LF18, Eq. 2.38).

// C depends on turbulence model closure coefficients. Below C is written in terms of the closure coefficients in OpenFOAM for the different models (Note that the naming of the closure coefficients is different from LF2018):

//kOmega: C = beta/(Cmu*gamma) = 1.54; // See LF18, Eq. 2.35, taking Cmu=beta^* and gamma=alpha

//kOmegaSST: C = beta1/(betaStar*gamma1) = 1.50; // See LF18, Eq. A.8
//kEpsilon: C = C2/C1= 1.33; // See LF18, Eq. A.16
//RNGkEpsilon: C = 1.26*C2/C1 = 1.49; // See LF18, Eq. A.23
// The value given in the default fvOptions file above comes from the closure coefficients of the Wilcox (2006) k-omega model.

// alpha controls the strength of the buoyancy production term, corresponding to variable alpha_b^* in LF18. If alpha=0 is used, then the buoyancy production term will be inactive. (Note for clarity that alpha here is different from the closure coefficient alpha in the Wilcox (2006) k-omega model.)

To demonstrate the usage of the new method a tutorial is provided here. The tutorial is a simple propagating wave case. The domain is exactly one wave length long with cyclic boundaries, such that the wave propagates out one end of the domain and enters the other end. The bottom has a slip condition, making the computational situation as close to potential flow as possible, such that a turbulence closure model should ideally not influence the physics of the wave propagation.
The wave is the same as the one used in the test cases in LF18 (e.g. their Section 3.2) and corresponds to the incoming wave of the Ting and Kirby (1994) spilling breaking wave experiment. The wave is a stream function wave with wave height H=0.125 m, water depth h=0.4 m and period T=2 s. The wave field has been initialized with waves2FOAM (Jacobsen et al., 2012), but the simulation is set up to run with interFoam.
To run the tutorial download and unzip the setup file, and just type the following in the terminal:

interFoam >& log&

The tutorial demonstrates that the initial turbulence levels decay as the wave propagates. In contrast, if lambda2 is set to 0.0 in the fvOptions file, then this will correspond to a “standard” turbulence model, and the turbulent kinetic energy will increase exponentially. The tutorial is set up to utilize the k-omega turbulence model, but will run similarly for any of the turbulence models indicated above.
For animations of the difference between standard and stabilized turbulence models for both propagating waves and breaking waves see Larsen & Fuhrman (2019):

In past simulations of free-surface waves (both breaking and non-breaking) using RANS models, there have been a marked a collective tendency to over-estimate the turbulence levels. In cases involving breaking waves, this has even been most pronounced prior to breaking with turbulence levels pre-breaking of same order of magnitude as in the surf-zone. The underlying cause was originally identified by Mayer & Madsen (2000) who showed that turbulence production exists in potential flow waves. They further showed that if omega falls below a certain threshold, the model becomes unstable (i.e. the models are conditionally unstable), leading to an exponential growth in both the eddy viscosity and turbulent kinetic energy.
In Larsen & Fuhrman (2018) it was proved that (1) standard two-equation closures (k-omega and k-epsilon) are unconditionally unstable (they will always become unstable) and (2) they can be simply and elegantly stabilized, thus solving this long-standing wide-spread problem.

Please cite as:

Larsen, B.E. and Fuhrman, D.R. (2020) Tutorial for "multiphaseStabilizedTurbulence" in OpenFOAM® - v1912. figshare. Media. DOI:

Devolder, B., Rauwoens, P. & Troch, P. (2017) Application of a buoyancy-modified k-omega turbulence model to simulate wave run-up around a monopile subject to regular waves using OpenFOAM. Coast. Eng. 125, 81-94. DOI: 10.1016/j.coastaleng.2017.04.004

Fuhrman D.R. & Larsen, B.E. (2020) A discussion on “Numerical computations of resonant sloshing using the modified isoAdvector method and the buoyancy-modified turbulence closure model” [Appl. Ocean Res. (2019), 93, article no. 101829, doi:10.1016.j.apor.2019.05.014]. Appl. Ocean Res. 102159. DOI: 10.1016/j.apor.2020.102159

Jacobsen, N.G., Fuhrman, D.R. & Fredsøe, J. (2012) A wave generation toolbox for the open-source CFD library: OpenFoam. Int. J. Numer. Meth. Fluids. 70, 1073-1088. DOI: 10.1002/fld.2726

Larsen, B.E. & Fuhrman, D.R. (2018) On the over-production of turbulence beneath surface waves in Reynolds-averaged Navier-Stokes models. J. Fluid Mech. 853, 419-460. DOI: 10.1017/jfm.2018.577

Larsen, B.E. & Fuhrman, D.R. (2019) Animations comparing surface wave simulations with standard and "stabilized" two-equation turbulence models. figshare. Media. DOI: 10.11583/DTU.8180708

Mayer, S. & Madsen, P.A. (2000) Simulation of breaking waves in the surf zone using a Navier-Stokes solver. In: Proc. 27th Int. Conf. Coast. Eng., pp. 928-941, Sydney, Australia. DOI: 10.1061/40549(276)72

OpenFOAM (2019)

Ting, F.C.K. & Kirby, J.T. (1994) Observations of undertow and turbulence in a laboratory surf zone. Coast. Eng. 24, 51-80. DOI: 10.1016/0378-3839(94)90026-4

Wilcox, D.C. (2006). Turbulence modeling for CFD. 3rd ed. DCW Industries. La Canada CA.


SWASH: 8022-00137B