Discretisation of the k-epsilon model for neutral athmospheric flows for OpenFOAM
 

Hello everebody out there,

I'm corrently struggelingwith the k-esplion model applied to neutral athmosperic boundary layers.

In order to check the numerical solution I set up a Couette flow in OpenfFOAM (which matches the zero Pressure gradient assumption) with periodic boundary conditions.

At the walls I applied
nutkAtmRoughWallFunction, the no slip for the velocity, zeroGradient for k
and for espilon a fixeed value of

epsilon[cellI] += w*Cmu75*pow(k[cellI], 1.5)/(kappa_*(y[faceI] + z0_[faceI]));
G[cellI] +=
w
*(nutw[faceI] + nuw[faceI])
*magGradUw[faceI]
*Cmu25*sqrt(k[cellI])
/(kappa_*(y[faceI] + z0_[faceI]));

.

SigmaEsp I set equal to kappa^2/(Cm^0.5 * (C2 - C1)). I thought that the Log-profile would result for the velocity. But I realized that the profiel was too steep. The reason is that the inclination turbulent length scale was to low.

The problem probably is due to a non neglectable discretization error. See

Richards, P. J., and S. E. Norris. "Appropriate boundary conditions for computational wind engineering models revisited." Journal of Wind Engineering and Industrial Aerodynamics 99.4 (2011): 257-266.

or

k-epsilon simulations of the neutral ABL: Achieving horizontal homogeneity on practical grid, Sumner and Masson, 2010.

Did anyone solve this issue.

I tried to implement the solution of Richard and Norris (averaging the turbulent shear stress over the cell faces) in OpenFoam. Unfortunatly it didnt work.

Here the relevant piece of the sourcecode:

surfaceVectorField GradU (fvc::snGrad(U_));

dimensionedScalar dims ( "dims", dimensionSet(0,0,-1,0,0,0,0), scalar(1.0));

surfaceVectorField GradUsqr (fvc::snGrad(U_)*dims);

forAll(GradUsqr, i)
{
GradUsqr[i].component(0) = GradU[i].component(0)*GradU[i].component(0);
GradUsqr[i].component(1) = GradU[i].component(1)*GradU[i].component(1);
GradUsqr[i].component(2) = GradU[i].component(2)*GradU[i].component(2);
}
// GradUsqr = sqr (GradU);
volVectorField snGradUSqr("snGradUSqr",fvc::average( GradUsqr ) );
volScalarField G(GName(), nut_*( mag ( snGradUSqr ) ) );

Maybe someone has a suggestion how to solve this issue.