ENH: improve/verify atmBoundaryLayerInlet conditions

ENH: add generalised log-law type ground-normal inflow boundary conditions for
  wind velocity and turbulence quantities for homogeneous, two-dimensional,
  dry-air, equilibrium and neutral atmospheric boundary layer (ABL) modelling

  ENH: remove `zGround` entry, which is now automatically computed

  ENH: add `displacement height` entry, `d`

  ENH: add generalised atmBoundaryLayerInletOmega boundary condition

  ENH: add a verification case for atmBoundaryLayerInlet BCs

  DOC: improve atmBoundaryLayerInlet header documentation

  BUG: fix value-entry behaviour in atmBoundaryLayerInlet (fixes #1578)
  Without this change:
  - for serial-parallel computations, if `value` entry is available in
    an `atmBoundaryLayerInlet` BC, the theoretical ABL profile expressions
    are not computed, and the `value` entry content is used as a profile data
  - for parallel computations, if `value` entry is not available, `decomposePar`
    could not be executed.
  With this change:
  - assuming `value` entry is always be present, the use of `value` entry for
    the ABL profile specification is determined by a flag `initABL`
  - the default value of the optional flag `initABL` is `true`, but whenever
    `initABL=true` is executed, `initABL` is overwritten as `false` for the
    subsequent runs, so that `value` entry can be safely used.
  Thanks Per Jørgensen for the bug report.

  BUG: ensure atmBoundaryInlet conditions are Galilean-invariant (fixes #1692)

  Related references:

      The ground-normal profile expressions (tag:RH):
        Richards, P. J., & Hoxey, R. P. (1993).
        Appropriate boundary conditions for computational wind
        engineering models using the k-ε turbulence model.
        In Computational Wind Engineering 1 (pp. 145-153).
        DOI:10.1016/B978-0-444-81688-7.50018-8

    Modifications to preserve the profiles downstream (tag:HW):
        Hargreaves, D. M., & Wright, N. G. (2007).
        On the use of the k–ε model in commercial CFD software
        to model the neutral atmospheric boundary layer.
        Journal of wind engineering and
        industrial aerodynamics, 95(5), 355-369.
        DOI:10.1016/j.jweia.2006.08.002

    Expression generalisations to allow height
    variation for turbulence quantities (tag:YGCJ):
        Yang, Y., Gu, M., Chen, S., & Jin, X. (2009).
        New inflow boundary conditions for modelling the neutral equilibrium
        atmospheric boundary layer in computational wind engineering.
        J. of Wind Engineering and Industrial Aerodynamics, 97(2), 88-95.
        DOI:10.1016/j.jweia.2008.12.001

    The generalised ground-normal profile expression for omega (tag:YGJ):
        Yang, Y., Gu, M., & Jin, X., (2009).
        New inflow boundary conditions for modelling the
        neutral equilibrium atmospheric boundary layer in SST k-ω model.
        In: The Seventh Asia-Pacific Conference on Wind Engineering,
        November 8-12, Taipei, Taiwan.

  Reproduced benchmark:
      Rectangular prism shown in FIG 1 of
        Hargreaves, D. M., & Wright, N. G. (2007).
        On the use of the k–ε model in commercial CFD software
        to model the neutral atmospheric boundary layer.
        Journal of wind engineering and
        industrial aerodynamics, 95(5), 355-369.
        DOI:10.1016/j.jweia.2006.08.002
  Benchmark data:
      HW, 2007 FIG 6

  TUT: update simpleFoam/turbineSiting tutorial accordingly
This commit is contained in:
Kutalmis Bercin
2020-04-07 15:41:37 +01:00
committed by Andrew Heather
parent 5863c94be0
commit 336fb3bddf
44 changed files with 2631 additions and 166 deletions

View File

@ -0,0 +1,66 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v1912 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
inlet
{
type atmBoundaryLayerInletVelocity;
#include "include/ABLConditions"
value uniform (0 0 0);
}
ground
{
type noSlip;
}
top
{
// (HW:p. 365):
// "In addition, as suggested by RH and often ignored by others, a"
// "constant shear stress of rho*(u^*)^2 was applied at the top boundary."
//
// u^* ~ Uref*kappa/ln((Zref+z0)/z0)
// (HW:Table 1):
// Uref = 10 m/s
// Zref = 6 m
// z0 = 0.01 m
// tau = rho*(u^*)^2 = 0.390796574
type fixedShearStress;
tau (0.390796574 0 0);
value uniform (0 0 0);
}
sides
{
type symmetry;
}
outlet
{
type inletOutlet;
inletValue uniform (0 0 0);
value $internalField;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,59 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v1912 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object epsilon;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -3 0 0 0 0];
internalField uniform 0.01;
boundaryField
{
#include "include/ABLConditions"
inlet
{
type atmBoundaryLayerInletEpsilon;
#include "include/ABLConditions"
value uniform 0;
}
ground
{
type epsilonWallFunction;
Cmu $Cmu;
kappa $kappa;
value $internalField;
}
top
{
type zeroGradient;
}
sides
{
type symmetry;
}
outlet
{
type inletOutlet;
inletValue $internalField;
value $internalField;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,19 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v1912 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
kappa 0.40; // (HW:p. 358)
Cmu 0.09; // (HW:p. 358)
flowDir (1 0 0); // (HW:Fig. 1)
zDir (0 0 1); // (HW:Fig. 1)
Uref 10.0; // (HW:Table 1)
Zref 6.0; // (HW:Table 1)
z0 uniform 0.01; // (HW:Table 1)
d uniform 0.0;
// ************************************************************************* //

View File

@ -0,0 +1,55 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v1912 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object k;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 1.3; // (HW:Eq. 6)
boundaryField
{
inlet
{
type atmBoundaryLayerInletK;
#include "include/ABLConditions"
value uniform 0;
}
ground
{
type kqRWallFunction;
value uniform 0.0;
}
top
{
type zeroGradient;
}
sides
{
type symmetry;
}
outlet
{
type inletOutlet;
inletValue $internalField;
value $internalField;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,59 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v1912 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object nut;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -1 0 0 0 0];
internalField uniform 0;
boundaryField
{
#include "include/ABLConditions"
inlet
{
type calculated;
value uniform 0;
}
ground
{
type nutkAtmRoughWallFunction;
kappa $kappa;
Cmu $Cmu;
z0 $z0;
value uniform 0.0;
}
top
{
type calculated;
value uniform 0;
}
sides
{
type symmetry;
}
outlet
{
type calculated;
value uniform 0;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,59 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v1912 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object omega;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 -1 0 0 0 0];
internalField uniform 0.0007;
boundaryField
{
#include "include/ABLConditions"
inlet
{
type atmBoundaryLayerInletOmega;
#include "include/ABLConditions"
value uniform 0;
}
ground
{
type omegaWallFunction;
Cmu $Cmu;
kappa $kappa;
value $internalField;
}
top
{
type zeroGradient;
}
sides
{
type symmetry;
}
outlet
{
type inletOutlet;
inletValue $internalField;
value $internalField;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,51 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v1912 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
inlet
{
type zeroGradient;
}
ground
{
type zeroGradient;
}
top
{
type zeroGradient;
}
sides
{
type symmetry;
}
outlet
{
type uniformFixedValue;
uniformValue constant 0;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,12 @@
#!/bin/sh
cd "${0%/*}" || exit # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/CleanFunctions # Tutorial clean functions
#------------------------------------------------------------------------------
cleanCase0
rm -rf plots
rm -f constant/turbulenceProperties
rm -rf system/atm-HargreavesWright-2007
#------------------------------------------------------------------------------

View File

@ -0,0 +1,21 @@
#!/bin/sh
cd "${0%/*}" || exit # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions # Tutorial run functions
#------------------------------------------------------------------------------
restore0Dir
cp -rf $FOAM_TUTORIALS/resources/dataset/atm-HargreavesWright-2007 system/.
rasModel="kEpsilon" # Tested options="kOmegaSST","kEpsilon"
sed "s|RAS_MODEL|$rasModel|g" constant/turbulenceProperties.template > \
constant/turbulenceProperties
runApplication blockMesh
runApplication renumberMesh -overwrite
runApplication $(getApplication)
#------------------------------------------------------------------------------

View File

@ -0,0 +1,21 @@
#!/bin/sh
cd "${0%/*}" || exit # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions # Tutorial run functions
#------------------------------------------------------------------------------
restore0Dir
rasModel="kEpsilon" # Tested options="kOmegaSST","kEpsilon"
sed "s|RAS_MODEL|$rasModel|g" constant/turbulenceProperties.template > \
constant/turbulenceProperties
runApplication blockMesh
runApplication decomposePar
runParallel renumberMesh -overwrite
runParallel $(getApplication)
#------------------------------------------------------------------------------

View File

@ -0,0 +1,117 @@
#------------------------------------------------------------------------------
Overview
"By setting appropriate profiles for wind velocity and the turbulence
quantities at the inlet, it is often assumed that the boundary layer will
be maintained up to the buildings or obstructions in the flow." (HW:p. 355).
However, it was quantified by (HW:p. 355) that "even in the absence of
obstructions, ..., the velocity and turbulence profiles decay along the
fetch" (HW:p. 355). It was shown by (HW:p. 355) that a set of modifications
were required to maintain a neutral atmospheric boundary layer throughout
an empty and long computational domain of a RANS computation.
Aim:
Verification of the following boundary conditions in terms of the
maintenance of inlet quantities downstream within a RANS computation:
- atmBoundaryLayerInletVelocity
- atmBoundaryLayerInletK
- atmBoundaryLayerInletEpsilon
- atmBoundaryLayerInletOmega
Benchmark (Physical phenomenon):
The benchmark is an empty fetch computational domain, steady-state
RANS simulation involving the following traits:
- External flow
- The surface layer portion of the neutral-stratified equilibrium
atmospheric boundary layer (no Ekman layer)
- Dry air
- Homogeneous, smooth terrain
- Spatiotemporal-invariant aerodynamic roughness length
- No displacement height
- Newtonian, single-phase, incompressible, non-reacting
Benchmark scenario:
- Computational domain: (HW:Fig. 1)
- Benchmark dataset: (HW:Fig. 6) (Obtained by the WebPlotDigitizer-4.2
(Rohatgi, 2019))
Resources:
Computational study (tag:HW):
Hargreaves, D. M., & Wright, N. G. (2007).
On the use of the kε model in commercial CFD software
to model the neutral atmospheric boundary layer.
Journal of wind engineering and
industrial aerodynamics, 95(5), 355-369.
DOI:10.1016/j.jweia.2006.08.002
Wind profile (tag:RQP):
Richards, P. J., Quinn, A. D., & Parker, S. (2002).
A 6 m cube in an atmospheric boundary layer flow-Part 2.
Computational solutions. Wind and structures, 5(2_3_4), 177-192.
DOI:10.12989/was.2002.5.2_3_4.177
Physical modelling:
- The governing equations for:
- Steady-state, Newtonian, single-phase, incompressible fluid flows,
excluding any thermal chemical, electromagnetic and scalar
interactions
- Mathematical approach for the turbulence modelling:
- Reynolds-averaged Navier-Stokes simulation (RANS)
- Turbulence closure model:
- kEpsilon and kOmegaSST linear eddy viscosity closure models
- The sets of input (HW:Table 1):
- Reference height, Zref = 6 [m]
- Aerodynamic roughness height, z0 = 0.01 [m]
- Displacement height, d = 0 [m]
- Reference mean wind speed, Uref = 10 [m/s]
Computational domain modelling:
- Rectangular prism
- (x1, x2, x3)
= (5000, 100, 500) [m]
= (streamwise, spanwise, ground-normal) directions
Computational domain discretisation:
- Spatial resolution:
- (x1, x2, x3) = (500, 5, 50) [cells]
- Refer to the `system/blockMeshDict` for the grading details
- Temporal resolution: Steady state
Equation discretisation:
- Spatial derivatives and variables:
- Convection: Second order
- Others: Second order with various limiters
- Temporal derivatives and variables: First order
Numerical boundary/initial conditions:
- Refer to `0.orig`
Pressure-velocity coupling algorithm:
- SIMPLEC
Linear solvers:
- Refer to `system/fvSolution`
Initialisation and sampling:
- No initialisation/averaging
- Sampling at the end of the simulation via `system/sampleDict`
- Refer to `system/controlDict` for further details
#------------------------------------------------------------------------------

View File

@ -0,0 +1,22 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v1912 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
transportModel Newtonian;
nu 1.5e-05;
// ************************************************************************* //

View File

@ -0,0 +1,42 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v1912 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object turbulenceProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
simulationType RAS;
RAS
{
RASModel RAS_MODEL;
turbulence on;
printCoeffs on;
RAS_MODELCoeffs // Valid only for epsilon-based models
{
Cmu 0.09;
C1 1.44;
C2 1.92;
sigmaEps 1.11; //Original value:1.44
// See p. 358:
// Hargreaves, D. M., & Wright, N. G. (2007).
// On the use of the kε model in commercial CFD software
// to model the neutral atmospheric boundary layer.
// Journal of wind engineering and
// industrial aerodynamics, 95(5), 355-369.
// DOI:10.1016/j.jweia.2006.08.002
}
}
// ************************************************************************* //

View File

@ -0,0 +1,339 @@
#!/bin/sh
cd "${0%/*}" || exit # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions # Tutorial run functions
#------------------------------------------------------------------------------
# Postprocess according to the existence of "epsilon" or "omega"
baseEpsilonOrOmega="epsilon" # options: "epsilon", "omega".
# Note: Benchmark data is available for the standard k-epsilon model from:
# Hargreaves, D. M., & Wright, N. G. (2007).
# On the use of the kε model in commercial CFD software
# to model the neutral atmospheric boundary layer.
# Journal of wind engineering and
# industrial aerodynamics, 95(5), 355-369.
# DOI:10.1016/j.jweia.2006.08.002
# Figure 6.
#------------------------------------------------------------------------------
plotUxUpstream() {
timeDir=$1
zMin=$2
echo " Plotting the ground-normal flow speed profile (upstream)."
outName="plots/Ux_upstream.png"
gnuplot<<PLT_UX_UPSTREAM
set terminal pngcairo font "helvetica,20" size 1000, 800
set xrange [4:18]
set yrange [0:50]
set grid
set key top left
set xlabel "Ux [m s^{-1}]"
set ylabel "Non-dimensionalised height, z/z_{ref}"
set offset .05, .05
set output "$outName"
zRef = 6
inp0="$timeDir/x_0mCell_U.xy"
inp1="$timeDir/x_0mPatch_U.xy"
plot \
"system/benchmark-data-HargreavesWright2007/Ux-HW-RH-Fig6a" \
u 1:2 t "RH" w p ps 2 pt 6 lc rgb "#000000", \
"system/benchmark-data-HargreavesWright2007/Ux-HW-RH-Fig6a" \
u 1:2 t "HW, x=2500m" w p ps 1 pt 5 lc rgb "#E69F00", \
"system/benchmark-data-HargreavesWright2007/Ux-HW-RH-Fig6a" \
u 1:2 t "HW, x=4000m" w p ps 0.5 pt 4 lc rgb "#56B4E9", \
inp0 u 2:((\$1-$zMin)/zRef) t "OF, x=0m (Patch)" w l lw 2 lc rgb "#009E73", \
inp1 u 2:((\$1-$zMin)/zRef) t "OF, x=0m (Cell)" w l lw 2 lc rgb "#F0E440"
PLT_UX_UPSTREAM
}
plotUxMid() {
timeDir=$1
zMin=$2
echo " Plotting the ground-normal flow speed profile (mid-range)."
outName="plots/Ux_mid.png"
gnuplot<<PLT_UX_MID
set terminal pngcairo font "helvetica,20" size 1000, 800
set xrange [4:18]
set yrange [0:50]
set grid
set key top left
set xlabel "Ux [m s^{-1}]"
set ylabel "Non-dimensionalised height, z/z_{ref}"
set offset .05, .05
set output "$outName"
zRef = 6
inp2="$timeDir/x_2500m_U.xy"
inp3="$timeDir/x_4000m_U.xy"
plot \
"system/benchmark-data-HargreavesWright2007/Ux-HW-RH-Fig6a" \
u 1:2 t "RH" w p ps 2 pt 6 lc rgb "#000000", \
"system/benchmark-data-HargreavesWright2007/Ux-HW-RH-Fig6a" \
u 1:2 t "HW, x=2500m" w p ps 1 pt 5 lc rgb "#E69F00", \
"system/benchmark-data-HargreavesWright2007/Ux-HW-RH-Fig6a" \
u 1:2 t "HW, x=4000m" w p ps 0.5 pt 4 lc rgb "#56B4E9", \
inp2 u 2:((\$1-$zMin)/zRef) t "OF, x=2500m" w l lw 2 lc rgb "#0072B2", \
inp3 u 2:((\$1-$zMin)/zRef) t "OF, x=4000m" w l lw 2 lc rgb "#D55E00"
PLT_UX_MID
}
plotUxDownstream() {
timeDir=$1
zMin=$2
echo " Plotting the ground-normal flow speed profile (downstream)."
outName="plots/Ux_downstream.png"
gnuplot<<PLT_UX_DOWNSTREAM
set terminal pngcairo font "helvetica,20" size 1000, 800
set xrange [4:18]
set yrange [0:50]
set grid
set key top left
set xlabel "Ux [m s^{-1}]"
set ylabel "Non-dimensionalised height, z/z_{ref}"
set offset .05, .05
set output "$outName"
zRef = 6
inp4="$timeDir/x_5000mCell_U.xy"
inp5="$timeDir/x_5000mPatch_U.xy"
plot \
"system/benchmark-data-HargreavesWright2007/Ux-HW-RH-Fig6a" \
u 1:2 t "RH" w p ps 2 pt 6 lc rgb "#000000", \
"system/benchmark-data-HargreavesWright2007/Ux-HW-RH-Fig6a" \
u 1:2 t "HW, x=2500m" w p ps 1 pt 5 lc rgb "#E69F00", \
"system/benchmark-data-HargreavesWright2007/Ux-HW-RH-Fig6a" \
u 1:2 t "HW, x=4000m" w p ps 0.5 pt 4 lc rgb "#56B4E9", \
inp4 u 2:((\$1-$zMin)/zRef) t "OF, x=5000m (Cell)" w l lw 2 lc rgb "#CC79A7", \
inp5 u 2:((\$1-$zMin)/zRef) t "OF, x=5000m (Patch)" w l lw 2 lc rgb "#440154"
PLT_UX_DOWNSTREAM
}
plotK() {
timeDir=$1
items=$2
seq=$3
zMin=$4
echo " Plotting the ground-normal turbulent kinetic energy profile."
outName="plots/k.png"
gnuplot<<PLT_K
set terminal pngcairo font "helvetica,20" size 1000, 800
set xrange [1:2]
set yrange [0:50]
set grid
set key top right
set xlabel "k [m^2 s^{-2}]"
set ylabel "Non-dimensionalised height, z/z_{ref}"
set offset .05, .05
set output "$outName"
zRef = 6
inp0="$timeDir/x_0mCell_$items.xy"
inp1="$timeDir/x_0mPatch_$items.xy"
inp2="$timeDir/x_2500m_$items.xy"
inp3="$timeDir/x_4000m_$items.xy"
inp4="$timeDir/x_5000mCell_$items.xy"
inp5="$timeDir/x_5000mPatch_$items.xy"
plot \
"system/benchmark-data-HargreavesWright2007/k-RH-Fig6b" \
u 1:2 t "RH" w p ps 2 pt 6 lc rgb "#000000", \
"system/benchmark-data-HargreavesWright2007/k-HW-Fig6b-2500" \
u 1:2 t "HW, x=2500m" w p ps 1 pt 5 lc rgb "#E69F00", \
"system/benchmark-data-HargreavesWright2007/k-HW-Fig6b-4000" \
u 1:2 t "HW, x=4000m" w p ps 0.5 pt 4 lc rgb "#56B4E9", \
inp0 u $seq:((\$1-$zMin)/zRef) t "OF, x=0m (Patch)" w l lw 2 lc rgb "#009E73", \
inp1 u $seq:((\$1-$zMin)/zRef) t "OF, x=0m (Cell)" w l lw 2 lc rgb "#F0E440", \
inp2 u $seq:((\$1-$zMin)/zRef) t "OF, x=2500m" w l lw 2 lc rgb "#0072B2", \
inp3 u $seq:((\$1-$zMin)/zRef) t "OF, x=4000m" w l lw 2 lc rgb "#D55E00", \
inp4 u $seq:((\$1-$zMin)/zRef) t "OF, x=5000m (Cell)" w l lw 2 lc rgb "#CC79A7", \
inp5 u $seq:((\$1-$zMin)/zRef) t "OF, x=5000m (Patch)" w l lw 2 lc rgb "#440154"
PLT_K
}
plotEpsilon() {
timeDir=$1
items=$2
zMin=$3
echo " Plotting the ground-normal turbulent kinetic"\
"energy dissipation rate profile."
outName="plots/epsilon.png"
gnuplot<<PLT_EPSILON
set terminal pngcairo font "helvetica,20" size 1000, 800
set xrange [0.001:10]
set yrange [0:50]
set grid
set key top right
set xlabel "epsilon [m^2 s^{-3}]"
set ylabel "Non-dimensionalised height, z/z_{ref}"
set offset .05, .05
set logscale x
set output "$outName"
zRef = 6
inp0="$timeDir/x_0mCell_$items.xy"
inp1="$timeDir/x_0mPatch_$items.xy"
inp2="$timeDir/x_2500m_$items.xy"
inp3="$timeDir/x_4000m_$items.xy"
inp4="$timeDir/x_5000mCell_$items.xy"
inp5="$timeDir/x_5000mPatch_$items.xy"
plot \
"system/benchmark-data-HargreavesWright2007/epsilon-HW-RH-Fig6c" \
u 1:2 t "RH" w p ps 2 pt 6 lc rgb "#000000", \
"system/benchmark-data-HargreavesWright2007/epsilon-HW-RH-Fig6c" \
u 1:2 t "HW, x=2500m" w p ps 1 pt 5 lc rgb "#E69F00", \
"system/benchmark-data-HargreavesWright2007/epsilon-HW-RH-Fig6c" \
u 1:2 t "HW, x=4000m" w p ps 0.5 pt 4 lc rgb "#56B4E9", \
inp0 u 2:((\$1-$zMin)/zRef) t "OF, x=0m (Patch)" w l lw 2 lc rgb "#009E73", \
inp1 u 2:((\$1-$zMin)/zRef) t "OF, x=0m (Cell)" w l lw 2 lc rgb "#F0E440", \
inp2 u 2:((\$1-$zMin)/zRef) t "OF, x=2500m" w l lw 2 lc rgb "#0072B2", \
inp3 u 2:((\$1-$zMin)/zRef) t "OF, x=4000m" w l lw 2 lc rgb "#D55E00", \
inp4 u 2:((\$1-$zMin)/zRef) t "OF, x=5000m (Cell)" w l lw 2 lc rgb "#CC79A7", \
inp5 u 2:((\$1-$zMin)/zRef) t "OF, x=5000m (Patch)" w l lw 2 lc rgb "#440154"
PLT_EPSILON
}
plotOmega() {
timeDir=$1
items=$2
zMin=$3
echo " Plotting the ground-normal specific dissipation rate profile."
outName="plots/omega.png"
gnuplot<<PLT_OMEGA
set terminal pngcairo font "helvetica,20" size 1000, 800
set xrange [0.001:10]
set yrange [0:50]
set grid
set key top right
set xlabel "omega [s^{-1}]"
set ylabel "Non-dimensionalised height, z/z_{ref}"
set offset .05, .05
set logscale x
set output "$outName"
zRef = 6
inp0="$timeDir/x_0mCell_$items.xy"
inp1="$timeDir/x_0mPatch_$items.xy"
inp2="$timeDir/x_2500m_$items.xy"
inp3="$timeDir/x_4000m_$items.xy"
inp4="$timeDir/x_5000mCell_$items.xy"
inp5="$timeDir/x_5000mPatch_$items.xy"
plot \
inp0 u 4:((\$1-$zMin)/zRef) t "OF, x=0m (Patch)" w l lw 2 lc rgb "#009E73", \
inp1 u 4:((\$1-$zMin)/zRef) t "OF, x=0m (Cell)" w l lw 2 lc rgb "#F0E440", \
inp2 u 4:((\$1-$zMin)/zRef) t "OF, x=2500m" w l lw 2 lc rgb "#0072B2", \
inp3 u 4:((\$1-$zMin)/zRef) t "OF, x=4000m" w l lw 2 lc rgb "#D55E00", \
inp4 u 4:((\$1-$zMin)/zRef) t "OF, x=5000m (Cell)" w l lw 2 lc rgb "#CC79A7", \
inp5 u 4:((\$1-$zMin)/zRef) t "OF, x=5000m (Patch)" w l lw 2 lc rgb "#440154"
PLT_OMEGA
}
plotMut() {
timeDir=$1
items=$2
seq=$3
zMin=$4
echo " Plotting the ground-normal turbulent viscosity profile."
outName="plots/mut.png"
gnuplot<<PLT_MUT
set terminal pngcairo font "helvetica,20" size 1000, 800
set xrange [0:120]
set yrange [0:50]
set grid
set key bottom right
set xlabel "mu_t [Pa.s]"
set ylabel "Non-dimensionalised height, z/z_{ref}"
set offset .05, .05
set output "$outName"
zRef = 6
inp0="$timeDir/x_0mCell_$items.xy"
inp1="$timeDir/x_0mPatch_$items.xy"
inp2="$timeDir/x_2500m_$items.xy"
inp3="$timeDir/x_4000m_$items.xy"
inp4="$timeDir/x_5000mCell_$items.xy"
inp5="$timeDir/x_5000mPatch_$items.xy"
plot \
"system/benchmark-data-HargreavesWright2007/mut-RH-Fig6d" \
u 1:2 t "RH" w p ps 2 pt 6 lc rgb "#000000", \
"system/benchmark-data-HargreavesWright2007/mut-HW-Fig6d-2500" \
u 1:2 t "HW, x=2500m" w p ps 1 pt 5 lc rgb "#E69F00", \
"system/benchmark-data-HargreavesWright2007/mut-HW-Fig6d-4000" \
u 1:2 t "HW, x=4000m" w p ps 0.5 pt 4 lc rgb "#56B4E9", \
inp0 u $seq:((\$1-$zMin)/zRef) t "OF, x=0m (Patch)" w l lw 2 lc rgb "#009E73", \
inp1 u $seq:((\$1-$zMin)/zRef) t "OF, x=0m (Cell)" w l lw 2 lc rgb "#F0E440", \
inp2 u $seq:((\$1-$zMin)/zRef) t "OF, x=2500m" w l lw 2 lc rgb "#0072B2", \
inp3 u $seq:((\$1-$zMin)/zRef) t "OF, x=4000m" w l lw 2 lc rgb "#D55E00", \
inp4 u $seq:((\$1-$zMin)/zRef) t "OF, x=5000m (Cell)" w l lw 2 lc rgb "#CC79A7", \
inp5 u $seq:((\$1-$zMin)/zRef) t "OF, x=5000m (Patch)" w l lw 2 lc rgb "#440154"
PLT_MUT
}
#------------------------------------------------------------------------------
# Require gnuplot
command -v gnuplot >/dev/null || {
echo "gnuplot not found - skipping graph creation" 1>&2
exit 1
}
# The latestTime in postProcessing/samples
timeDir=$(foamListTimes -case postProcessing/samples -latestTime 2>/dev/null)
[ -n "$timeDir" ] || {
echo "No postProcessing/samples found - skipping graph creation" 1>&2
exit 2
}
timeDir="postProcessing/samples/$timeDir"
zMin=0
mkdir -p plots
# Postprocess flow speed
plotUxUpstream $timeDir $zMin
plotUxMid $timeDir $zMin
plotUxDownstream $timeDir $zMin
# Postprocess turbulence quantities
if [ $baseEpsilonOrOmega == "epsilon" ]
then
items="epsilon_k_nut"
plotEpsilon $timeDir $items $zMin
plotK $timeDir $items 3 $zMin
plotMut $timeDir $items 4 $zMin
elif [ $baseEpsilonOrOmega == "omega" ]
then
items="k_nut_omega"
plotK $timeDir $items 2 $zMin
plotMut $timeDir $items 3 $zMin
plotOmega $timeDir $items $zMin
else
echo "Chosen turbulence model is neither epsilon nor omega based." 1>&2
exit 2
fi
#------------------------------------------------------------------------------

View File

@ -0,0 +1,111 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v1912 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
scale 1;
// x = streamwise
// y = spanwise
// z = wall-normal
nx 500;
ny 5;
nz 50;
xMin 0;
xMax 5000.0;
yMin 0;
yMax 100.0;
zMin 0.0;
zMax 500.0;
// blockMesh calculator input:
// width of start cell = 1.0 (HW:p. 359)
// number of cells = 50
// total length = 500
// blockMesh calculator output:
// cell-to-cell expansion ratio = 1.076030437 (consistent with 1.076 (HW:Fig.1))
zTotalExpansion 36.25795062;
vertices
(
($xMin $yMin $zMin)
($xMax $yMin $zMin)
($xMax $yMax $zMin)
($xMin $yMax $zMin)
($xMin $yMin $zMax)
($xMax $yMin $zMax)
($xMax $yMax $zMax)
($xMin $yMax $zMax)
);
blocks
(
hex (0 1 2 3 4 5 6 7) ($nx $ny $nz) simpleGrading (1 1 $zTotalExpansion)
);
edges
(
);
boundary
(
inlet
{
type patch;
faces
(
(0 4 7 3)
);
}
ground
{
type wall;
faces
(
(0 3 2 1)
);
}
top
{
type patch;
faces
(
(4 5 6 7)
);
}
sides
{
type symmetry;
faces
(
(1 5 4 0)
(3 7 6 2)
);
}
outlet
{
type patch;
faces
(
(2 6 5 1)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //

View File

@ -0,0 +1,62 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v1912 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application simpleFoam;
startFrom latestTime;
startTime 0;
stopAt endTime;
endTime 5000;
deltaT 1;
writeControl timeStep;
writeInterval 500;
purgeWrite 0;
writeFormat ascii;
writePrecision 12;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable yes;
functions
{
#include "sampleDict"
minMax
{
type fieldMinMax;
libs (fieldFunctionObjects);
writeControl writeTime;
fields (U);
}
}
// ************************************************************************* //

View File

@ -0,0 +1,26 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v1912 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object decomposeParDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
numberOfSubdomains 8;
method hierarchical;
coeffs
{
n (8 1 1);
}
// ************************************************************************* //

View File

@ -0,0 +1,58 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v1912 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
ddtSchemes
{
default steadyState;
}
gradSchemes
{
default Gauss linear;
}
divSchemes
{
default none;
div(phi,U) bounded Gauss linear;
div(phi,epsilon) bounded Gauss limitedLinear 1;
div(phi,omega) bounded Gauss limitedLinear 1;
div(phi,k) bounded Gauss limitedLinear 1;
div((nuEff*dev2(T(grad(U))))) Gauss linear;
}
laplacianSchemes
{
default Gauss linear uncorrected;
}
interpolationSchemes
{
default linear;
}
snGradSchemes
{
default uncorrected;
}
wallDist
{
method meshWave;
}
// ************************************************************************* //

View File

@ -0,0 +1,59 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v1912 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
solvers
{
p
{
solver GAMG;
smoother GaussSeidel;
tolerance 1e-8;
relTol 0.01;
}
"(U|k|epsilon|omega)"
{
solver smoothSolver;
smoother GaussSeidel;
tolerance 1e-8;
relTol 0.01;
}
}
SIMPLE
{
nNonOrthogonalCorrectors 0;
consistent true;
pRefCell 0;
pRefValue 0;
}
relaxationFactors
{
equations
{
U 0.9;
"(k|epsilon|omega)" 0.7;
}
}
cache
{
grad(U);
}
// ************************************************************************* //

View File

@ -0,0 +1,79 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v1912 |
| \\ / A nd | Web: www.OpenFOAM.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location system;
object sampleDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// The locations of the sample profiles correspond to:
// Hargreaves-Wright (2007), Fig.6
// DOI:10.1016/j.jweia.2006.08.002
samples
{
type sets;
libs (sampling);
setFormat raw;
interpolationScheme cell;
fields (U k epsilon nut omega);
writeControl writeTime;
sets
(
x_0mPatch // inlet patch face centres
{
type face;
axis z;
start (0 50 0);
end (0 50 500);
}
x_0mCell // inlet-first cell centres
{
type midPoint;
axis z;
start (5.0 50 0);
end (5.0 50 500);
}
x_2500m
{
type face;
axis z;
start (2500 50 0);
end (2500 50 500);
}
x_4000m
{
type face;
axis z;
start (4000 50 0);
end (4000 50 500);
}
x_5000mCell // outlet patch face centres
{
type face;
axis z;
start (4995 50 0);
end (4995 50 500);
}
x_5000mPatch // outlet-first cell centres
{
type face;
axis z;
start (5000 50 0);
end (5000 50 500);
}
);
}
// *********************************************************************** //