ENH: New atmospheric boundary layer (ABL) model suite (Part 1)

Please refer to the header file documentation for complete set of details.

  ENH: add new fvOptions for ABL modelling

    - atmAmbientTurbSource
    - atmBuoyancyTurbSource
    - atmCoriolisUSource
    - atmLengthScaleTurbSource
    - atmPlantCanopyTurbSource
    - atmPlantCanopyUSource
    - atmPlantCanopyTSource
    - atmNutSource

  ENH: add new boundary conditions for ABL modelling
       with PatchFunction1 and TimeFunction1 support

    - atmAlphatkWallFunction
    - atmEpsilonWallFunction
    - atmNutkWallFunction
    - atmNutUWallFunction
    - atmNutWallFunction
    - atmOmegaWallFunction
    - atmTurbulentHeatFluxTemperature

  STYLE: change names of nutkAtmRoughWallFunction -> atmNutkWallFunction by
         ensuring the bitwise backward compatibility

  ENH: add new variable-scaling force computation method to actuationDiskSource

  ENH: review actuationDiskSource and radialActuationDiskSource

  ENH: add new function object, ObukhovLength

  ENH: add new ABL tutorials/verifications

    - verificationAndValidation/atmosphericModels/atmFlatTerrain
      - verification with the Leipzig field experiment
      - illustration of precursor/successor field mapping
    - verificationAndValidation/atmosphericModels/atmForestStability
      - verification with the Sweden field experiment
    - update incompressible/simpleFoam/turbineSiting
This commit is contained in:
Kutalmis Bercin
2020-05-05 14:20:45 +01:00
committed by Andrew Heather
parent 70cd6c6176
commit 41e264f27d
386 changed files with 53513 additions and 349 deletions

View File

@ -0,0 +1,46 @@
/*--------------------------------*- 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 T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 0 1 0 0 0];
internalField uniform 300;
boundaryField
{
bottom
{
type atmTurbulentHeatFluxTemperature;
heatSource flux;
alphaEff alphaEff;
Cp0 1005.0;
q uniform 0.0;
value uniform 300;
}
top
{
type fixedValue;
value uniform 300;
}
"inlet|outlet|left|right"
{
type cyclic;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,41 @@
/*--------------------------------*- 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 binary;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (17.5 0.0 0.0);
boundaryField
{
bottom
{
type fixedValue;
value uniform (0 0 0);
}
top
{
type slip;
}
"inlet|outlet|left|right"
{
type cyclic;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,47 @@
/*--------------------------------*- 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 alphat;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -1 0 0 0 0];
internalField uniform 0.0;
boundaryField
{
bottom
{
type atmAlphatkWallFunction;
Cmu 0.09;
kappa 0.4;
Pr 0.9;
z0 uniform 0.028;
Prt uniform 0.74;
value uniform 0.028;
}
top
{
type calculated;
value uniform 0;
}
"inlet|outlet|left|right"
{
type cyclic;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,44 @@
/*--------------------------------*- 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.1;
boundaryField
{
bottom
{
type atmEpsilonWallFunction;
kappa 0.40;
Cmu 0.03;
z0 uniform 0.028;
value uniform 0.028;
}
top
{
type slip;
}
"inlet|outlet|left|right"
{
type cyclic;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,41 @@
/*--------------------------------*- 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.285;
boundaryField
{
bottom
{
type kqRWallFunction;
value $internalField;
}
top
{
type slip;
}
"inlet|outlet|left|right"
{
type cyclic;
}
}
// ************************************************************************* //

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 volScalarField;
object leafAreaDensity;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 -1 0 0 0 0 0];
internalField uniform 0.14;
boundaryField
{
bottom
{
type fixedValue;
value $internalField;
}
top
{
type fixedValue;
value uniform 0;
}
"inlet|outlet|left|right"
{
type cyclic;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,46 @@
/*--------------------------------*- 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.0;
boundaryField
{
bottom
{
type atmNutkWallFunction;
boundNut false;
kappa 0.40;
Cmu 0.03;
z0 uniform 0.028;
value uniform 0.028;
}
top
{
type calculated;
value uniform 0;
}
"inlet|outlet|left|right"
{
type cyclic;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,44 @@
/*--------------------------------*- 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.077821;
boundaryField
{
bottom
{
type atmOmegaWallFunction;
kappa 0.40;
Cmu 0.03;
z0 uniform 0.028;
value uniform 0.028;
}
top
{
type slip;
}
"inlet|outlet|left|right"
{
type cyclic;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,38 @@
/*--------------------------------*- 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 binary;
class volScalarField;
object p_rgh;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0.0;
boundaryField
{
"bottom|top"
{
type fixedFluxPressure;
rho rhok;
value uniform 0;
}
"inlet|outlet|left|right"
{
type cyclic;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,36 @@
/*--------------------------------*- 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 Cd;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 0 0 0 0 0];
internalField uniform 0.2;
boundaryField
{
"bottom|top"
{
type fixedValue;
value uniform 0;
}
"inlet|outlet|left|right"
{
type cyclic;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,36 @@
/*--------------------------------*- 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 qPlant;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -3 0 0 0 0];
internalField uniform Q_PLANT;
boundaryField
{
"bottom|top"
{
type fixedValue;
value uniform 0;
}
"inlet|outlet|left|right"
{
type cyclic;
}
}
// ************************************************************************* //

View File

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

View File

@ -0,0 +1,75 @@
#!/bin/sh
cd "${0%/*}" || exit # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions # Tutorial run functions
. ${WM_PROJECT_DIR:?}/bin/tools/CleanFunctions # Tutorial clean functions
#------------------------------------------------------------------------------
./Allrun.pre
cp -rf $FOAM_TUTORIALS/resources/dataset/atm-Arnqvist-2015 system/.
RASmodel="kEpsilon" # kOmegaSST
declare -A stabilityStates
declare -A Lmaxs
declare -A qPlants
stabilityStates[0]="veryStable"
stabilityStates[1]="stable"
stabilityStates[2]="slightlyStable"
stabilityStates[3]="neutral"
stabilityStates[4]="slightlyUnstable"
stabilityStates[5]="unstable"
Lmaxs[0]="5.0"
Lmaxs[1]="13.0"
Lmaxs[2]="25.5"
Lmaxs[3]="41.0"
Lmaxs[4]="80.75"
Lmaxs[5]="200.0"
qPlants[0]="-20.0"
qPlants[1]="-9.0"
qPlants[2]="-5.0"
qPlants[3]="0.0"
qPlants[4]="15.0"
qPlants[5]="60.0"
for i in "${!stabilityStates[@]}"
do
state=${stabilityStates[$i]}
Lmax=${Lmaxs[$i]}
qPlant=${qPlants[$i]}
echo " # Computations for the atmopsheric stability = $state:"
echo " ## Lmax = $Lmax [m], qPlant = $qPlant [-]"
echo ""
sed -e "s|RAS_MODEL|$RASmodel|g" constant/turbulenceProperties.temp > \
constant/turbulenceProperties
sed -e "s|L_MAX|$Lmax|g" constant/fvOptions.temp > constant/fvOptions
sed -e "s|Q_PLANT|$qPlant|g" 0.orig/qPlant.temp > 0/qPlant
sed -e "s|Q_PLANT|$qPlant|g" system/setFieldsDict.temp > \
system/setFieldsDict
runApplication decomposePar
runParallel -s parallel renumberMesh -overwrite
runParallel -s $i setFields
runParallel $(getApplication)
runParallel postProcess -funcs \
"(minMaxComponents(U) minMaxMagnitude(U))" -latestTime
runParallel redistributePar -reconstruct -latestTime
# Collect results into $resultDir
resultDir="results/$state"
mkdir -p $resultDir
timeDir=$(foamListTimes -latestTime)
mv -f $timeDir postProcessing log.* $resultDir
cp -rf system/fv* constant/fv* constant/*Properties 0 $resultDir
cleanTimeDirectories
rm -rf processor*
done
#------------------------------------------------------------------------------

View File

@ -0,0 +1,14 @@
#!/bin/sh
cd "${0%/*}" || exit # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions # Tutorial run functions
#------------------------------------------------------------------------------
restore0Dir
rm -f 0/qPlant.temp
runApplication blockMesh
runApplication renumberMesh -overwrite
#------------------------------------------------------------------------------

View File

@ -0,0 +1,106 @@
/*--------------------------------*- 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;
location "constant";
object fvOptions;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
pressureGradient
{
type vectorSemiImplicitSource;
volumeMode specific;
selectionMode all;
injectionRateSuSp
{
U ((0 1.978046e-03 0) 0);
}
}
atmCoriolisUSource1
{
type atmCoriolisUSource;
atmCoriolisUSourceCoeffs
{
selectionMode all;
Omega (0 0 5.65156e-05);
}
}
atmAmbientTurbSource1
{
type atmAmbientTurbSource;
atmAmbientTurbSourceCoeffs
{
selectionMode all;
kAmb 1.0e-04;
epsilonAmb 7.208e-08;
}
}
atmBuoyancyTurbSource1
{
type atmBuoyancyTurbSource;
atmBuoyancyTurbSourceCoeffs
{
selectionMode all;
rho rho;
Lmax L_MAX;
n 3.0;
beta 3.3e-03;
}
}
atmLengthScaleTurbSource1
{
type atmLengthScaleTurbSource;
atmLengthScaleTurbSourceCoeffs
{
selectionMode all;
Lmax L_MAX;
n 3.0;
}
}
atmPlantCanopyUSource1
{
type atmPlantCanopyUSource;
atmPlantCanopyUSourceCoeffs
{
selectionMode all;
}
}
atmPlantCanopyTSource1
{
type atmPlantCanopyTSource;
atmPlantCanopyTSourceCoeffs
{
selectionMode all;
}
}
atmPlantCanopyTurbSource1
{
type atmPlantCanopyTurbSource;
atmPlantCanopyTurbSourceCoeffs
{
selectionMode all;
rho rho;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,23 @@
/*--------------------------------*- 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 uniformDimensionedVectorField;
location "constant";
object g;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -2 0 0 0 0];
value (0 0 -9.81);
// ************************************************************************* //

View File

@ -0,0 +1,34 @@
/*--------------------------------*- 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;
// Thermal expansion coefficient
beta beta [0 0 0 -1 0 0 0] 3e-03;
// Reference temperature
TRef 300;
// Laminar Prandtl number
Pr 0.9;
// Turbulent Prandtl number
Prt Prt [0 0 0 0 0 0 0] 0.7;
// ************************************************************************* //

View File

@ -0,0 +1,38 @@
/*--------------------------------*- 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
{
Cmu 0.09;
betaStar 0.09;
sigmaEps 1.30;
sigmaK 1;
C1 1.44;
C2 1.92;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,396 @@
#!/bin/bash
cd "${0%/*}" || exit # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions # Tutorial run functions
#------------------------------------------------------------------------------
# Benchmark dataset:
# Arnqvist, J., Segalini, A., Dellwik, E., & Bergström, H. (2015).
# Wind statistics from a forested landscape.
# Boundary-Layer Meteorology, 156(1), 53-71.
# DOI:10.1007/s10546-015-0016-x
#------------------------------------------------------------------------------
# Simple linear interpolation on an ascending table
linearInterp() {
inpFile=$1
zRef=$2
yi=($(awk -v zRef="$zRef" '
BEGIN{
x1=0; y1=0;
x2=$1; y2=$2;
}
{
if((x1 < zRef) && (x2 < zRef) && (getline != 0))
{
x1=x2; y1=y2; x2=$1; y2=$2
}
}
END {
yi = (y2-y1)/(x2-x1)*(zRef-x1) + y1;
print yi
}
' $inpFile))
echo "$yi"
}
plotU() {
timeDir=$1
treeH=$2
zRef=$3
smple=$4
shift 4
declare -a inps=("${!1}")
echo " Plotting the ground-normal normalised "
echo " streamwise flow speed profile."
for i in "${!inps[@]}"
do
inp="results/${inps[$i]}/$timeDir/$smple"
echo $inp
# Store the ground-normal height
z=($(awk '{ printf "%.16f\n", $1 }' $inp))
# Compute velocity magnitude by excluding the ground-normal component
magU=($(awk '
{
mag = sqrt($2*$2 + $3*$3);
printf "%.16f\n", mag
}' $inp))
file="z-magU-$i.dat"
[ ! -f $file ] && for ((j = 0; j < "${#magU[@]}"; j++)) \
do printf "%.16f %.16f\n" "${z[$j]}" "${magU[$j]}" \
>> $file; done
# Find velocity magnitude at zRef height by linear interpolation
uRef=$(linearInterp "$file" "$zRef")
# Find z normalised by the reference tree height
zNorm=($(awk -v treeH="$treeH" '
{
zNorm = $1/treeH; printf "%.16f\n", zNorm
}' $file))
# Find U normalised by uRef
uNorm=($(awk -v uRef="$uRef" '
{
uNorm = $2/uRef; printf "%.16f\n", uNorm
}' $file))
fil="zNorm-uNorm-$i.dat"
[ ! -f $fil ] && for ((j = 0; j < "${#zNorm[@]}"; j++)) \
do printf "%.16f %.16f\n" "${zNorm[$j]}" "${uNorm[$j]}" \
>> $fil; done
done
outName="plots/uNorm-zNorm.png"
gnuplot<<PLT_U
set terminal pngcairo font "helvetica,20" size 1000, 800
set xrange [0:3.5]
set yrange [0:7]
set grid
set key bottom right
set key samplen 2
set key spacing 0.75
set xlabel "U/U_{ref} [-]"
set ylabel "z/z_{ref} [-]"
set offset .05, .05
set output "$outName"
veryStable="zNorm-uNorm-0.dat"
stable="zNorm-uNorm-1.dat"
slightlyStable="zNorm-uNorm-2.dat"
neutral="zNorm-uNorm-3.dat"
slightlyUnstable="zNorm-uNorm-4.dat"
unstable="zNorm-uNorm-5.dat"
bench="system/atm-Arnqvist-2015/uNorm-zNorm.dat"
plot \
bench every ::0::5 u ((\$1)/13.8953):2 t "Vs" w p ps 2 pt 6 lc rgb "#000000", \
bench every ::6::11 u ((\$1)/8.13953):2 t "S" w p ps 2 pt 5 lc rgb "#000000", \
bench every ::12::17 u ((\$1)/6.36047):2 t "Ss" w p ps 2 pt 4 lc rgb "#000000", \
bench every ::18::23 u ((\$1)/5.62791):2 t "N" w p ps 2 pt 3 lc rgb "#000000", \
bench every ::24::29 u ((\$1)/5.31395):2 t "Su" w p ps 2 pt 2 lc rgb "#000000", \
bench every ::30::35 u ((\$1)/5.06977):2 t "U" w p ps 2 pt 1 lc rgb "#000000", \
veryStable u 2:1 t "Very stable" w l lw 2 lc rgb "#009E73", \
stable u 2:1 t "Stable" w l lw 2 lc rgb "#F0E440", \
slightlyStable u 2:1 t "Slightly stable" w l lw 2 lc rgb "#0072B2", \
neutral u 2:1 t "Neutral" w l lw 2 lc rgb "#D55E00", \
slightlyUnstable u 2:1 t "Slightly unstable" w l lw 2 lc rgb "#CC79A7", \
unstable u 2:1 t "Unstable" w l lw 2 lc rgb "#440154"
PLT_U
rm -f {zNorm-uNorm-,z-magU-}*.dat
}
plotK() {
timeDir=$1
treeH=$2
zRef=$3
smple=$4
shift 4
declare -a inps=("${!1}")
echo " Plotting the ground-normal normalised "
echo " turbulent kinetic energy profile."
for i in "${!inps[@]}"
do
inp="results/${inps[$i]}/$timeDir/$smple"
echo $inp
# Store the ground-normal height profile
z=($(awk '{ printf "%.16f\n", $1 }' $inp))
# Store the turbulent kinetic energy profile
k=($(awk '{ printf "%.16f\n", $5 }' $inp))
file="z-k-$i.dat"
[ ! -f $file ] && for ((j = 0; j < "${#k[@]}"; j++)) \
do printf "%.16f %.16f\n" "${z[$j]}" "${k[$j]}" \
>> $file; done
# Find k at zRef height by linear interpolation
kRef=$(linearInterp "$file" "$zRef")
# Find z normalised by the reference tree height
zNorm=($(awk -v treeH="$treeH" '
{
zNorm = $1/treeH; printf "%.16f\n", zNorm
}' $file))
# Find k normalised by kRef
kNorm=($(awk -v kRef="$kRef" '
{
kNorm = $2/kRef; printf "%.16f\n", kNorm
}' $file))
fil="zNorm-kNorm-$i.dat"
[ ! -f $fil ] && for ((j = 0; j < "${#zNorm[@]}"; j++)) \
do printf "%.16f %.16f\n" "${zNorm[$j]}" "${kNorm[$j]}" \
>> $fil; done
done
outName="plots/kNorm-zNorm.png"
gnuplot<<PLT_K
set terminal pngcairo font "helvetica,20" size 1000, 800
set xrange [0:2]
set yrange [0:7]
set grid
set key bottom right
set key samplen 2
set key spacing 0.75
set xlabel "k/k_{ref} [-]"
set ylabel "z/z_{ref} [-]"
set offset .05, .05
set output "$outName"
veryStable="zNorm-kNorm-0.dat"
stable="zNorm-kNorm-1.dat"
slightlyStable="zNorm-kNorm-2.dat"
neutral="zNorm-kNorm-3.dat"
slightlyUnstable="zNorm-kNorm-4.dat"
unstable="zNorm-kNorm-5.dat"
bench="system/atm-Arnqvist-2015/kNorm-zNorm.dat"
plot \
bench every ::0::5 u ((\$1)/5.050476682):2 t "Vs" w p ps 2 pt 6 lc rgb "#000000", \
bench every ::6::11 u ((\$1)/4.24970097):2 t "S" w p ps 2 pt 5 lc rgb "#000000", \
bench every ::12::17 u ((\$1)/3.897762917):2 t "Ss" w p ps 2 pt 4 lc rgb "#000000", \
bench every ::18::23 u ((\$1)/3.788680963):2 t "N" w p ps 2 pt 3 lc rgb "#000000", \
bench every ::24::29 u ((\$1)/4.038160328):2 t "Su" w p ps 2 pt 2 lc rgb "#000000", \
bench every ::30::35 u ((\$1)/4.198358216):2 t "U" w p ps 2 pt 1 lc rgb "#000000", \
veryStable u 2:1 t "Very stable" w l lw 2 lc rgb "#009E73", \
stable u 2:1 t "Stable" w l lw 2 lc rgb "#F0E440", \
slightlyStable u 2:1 t "Slightly stable" w l lw 2 lc rgb "#0072B2", \
neutral u 2:1 t "Neutral" w l lw 2 lc rgb "#D55E00", \
slightlyUnstable u 2:1 t "Slightly unstable" w l lw 2 lc rgb "#CC79A7", \
unstable u 2:1 t "Unstable" w l lw 2 lc rgb "#440154"
PLT_K
rm -f {zNorm-kNorm-,z-k-}*.dat
}
printObukhovLength() {
timeDir=$1
treeH=$2
zRef=$3
smple=$4
shift 4
declare -a inps=("${!1}")
echo " Printing the Obukhov length at a given reference height."
for i in "${!inps[@]}"
do
inp="results/${inps[$i]}/$timeDir/$smple"
echo $inp
# Store the ground-normal height profile
z=($(awk '{ printf "%.16f\n", $1 }' $inp))
# Store the Obukhov length profile
OL=($(awk '{ printf "%.16f\n", $2 }' $inp))
file="z-OL-$i.dat"
[ ! -f $file ] && for ((j = 0; j < "${#OL[@]}"; j++)) \
do printf "%.16f %.16f\n" "${z[$j]}" "${OL[$j]}" \
>> $file; done
# Find the Obukhov length at zRef height by linear interpolation
OLRef=$(linearInterp "$file" "$zRef")
echo "${inps[$i]} = $OLRef" >> "plots/ObukhovLength.dat"
done
rm -f z-OL-*.dat
}
plotVeer() {
timeDir=$1
treeH=$2
zRef=$3
smple=$4
shift 4
declare -a inps=("${!1}")
echo " Plotting the ground-normal normalised veer profile."
for i in "${!inps[@]}"
do
inp="results/${inps[$i]}/$timeDir/$smple"
echo $inp
# Store the ground-normal height
z=($(awk '{ printf "%.16f\n", $1 }' $inp))
# Store streamwise and spanwise velocity components
u=($(awk '{ printf "%.16f\n", $2 }' $inp))
v=($(awk '{ printf "%.16f\n", $3 }' $inp))
fileu="z-u-$i.dat"
[ ! -f $fileu ] && for ((j = 0; j < "${#z[@]}"; j++)) \
do printf "%.16f %.16f\n" "${z[$j]}" "${u[$j]}" \
>> $fileu; done
filev="z-v-$i.dat"
[ ! -f $filev ] && for ((j = 0; j < "${#z[@]}"; j++)) \
do printf "%.16f %.16f\n" "${z[$j]}" "${v[$j]}" \
>> $filev; done
# Find u and v at zRef height by linear interpolation
uRef=$(linearInterp "$fileu" "$zRef")
vRef=$(linearInterp "$filev" "$zRef")
# Find z normalised by the reference tree height
zNorm=($(awk -v treeH="$treeH" '
{
zNorm = $1/treeH; printf "%.16f\n", zNorm
}' $fileu))
# Find veer
veer=($(awk -v uRef="$uRef" -v vRef="$vRef" '
{
x = $2/sqrt($2*$2 + $3*$3);
xR = uRef/sqrt(uRef*uRef + vRef*vRef);
veer = -1*(atan2(sqrt(1 - x*x), x) - atan2(sqrt(1 - xR*xR), xR))*180/atan2(0, -1);
printf "%.16f\n", veer
}' $inp))
fil="zNorm-veer-$i.dat"
[ ! -f $fil ] && for ((j = 0; j < "${#zNorm[@]}"; j++)) \
do printf "%.16f %.16f\n" "${zNorm[$j]}" "${veer[$j]}" \
>> $fil; done
done
outName="plots/veer-zNorm.png"
gnuplot<<PLT_VEER
set terminal pngcairo font "helvetica,20" size 1000, 800
set xrange [-5:20]
set yrange [0:7]
set grid
set key bottom right
set key samplen 2
set key spacing 0.75
set xlabel "alpha/alpha_{ref} [-]"
set ylabel "z/z_{ref} [-]"
set offset .05, .05
set output "$outName"
veryStable="zNorm-veer-0.dat"
stable="zNorm-veer-1.dat"
slightlyStable="zNorm-veer-2.dat"
neutral="zNorm-veer-3.dat"
slightlyUnstable="zNorm-veer-4.dat"
unstable="zNorm-veer-5.dat"
bench="system/atm-Arnqvist-2015/veer-zNorm.dat"
plot \
bench every ::0::5 u 1:2 t "Vs" w p ps 2 pt 6 lc rgb "#000000", \
bench every ::6::11 u 1:2 t "S" w p ps 2 pt 5 lc rgb "#000000", \
bench every ::12::17 u 1:2 t "Ss" w p ps 2 pt 4 lc rgb "#000000", \
bench every ::18::23 u 1:2 t "N" w p ps 2 pt 3 lc rgb "#000000", \
bench every ::24::29 u 1:2 t "Su" w p ps 2 pt 2 lc rgb "#000000", \
bench every ::30::35 u 1:2 t "U" w p ps 2 pt 1 lc rgb "#000000", \
veryStable u 2:1 t "Very stable" w l lw 2 lc rgb "#009E73", \
stable u 2:1 t "Stable" w l lw 2 lc rgb "#F0E440", \
slightlyStable u 2:1 t "Slightly stable" w l lw 2 lc rgb "#0072B2", \
neutral u 2:1 t "Neutral" w l lw 2 lc rgb "#D55E00", \
slightlyUnstable u 2:1 t "Slightly unstable" w l lw 2 lc rgb "#CC79A7", \
unstable u 2:1 t "Unstable" w l lw 2 lc rgb "#440154"
PLT_VEER
rm -f {zNorm-veer-,z-u-,z-v-}*.dat
}
#------------------------------------------------------------------------------
# Requires gnuplot
command -v gnuplot >/dev/null || {
echo "gnuplot not found - skipping graph creation" 1>&2
exit 1
}
# Requires awk
command -v awk >/dev/null || {
echo "awk not found - skipping graph creation" 1>&2
exit 1
}
# The latestTime in postProcessing/samples
timeDir=$(foamListTimes -case results/veryStable/postProcessing/samples \
-latestTime 2>/dev/null)
[ -n "$timeDir" ] || {
echo "No results found in postProcessing - skipping graph creation" 1>&2
exit 1
}
timeDir="postProcessing/samples/$timeDir"
# Settings
declare -a stability
stability[0]="veryStable"
stability[1]="stable"
stability[2]="slightlyStable"
stability[3]="neutral"
stability[4]="slightlyUnstable"
stability[5]="unstable"
sample1="lineZ1_U.xy"
sample2="lineZ1_ObukhovLength_T_Ustar_k_p_rgh.xy"
treeHeight=20
zRef=40
# Postprocessing
mkdir -p plots
plotU $timeDir $treeHeight $zRef $sample1 stability[@]
plotK $timeDir $treeHeight $zRef $sample2 stability[@]
printObukhovLength $timeDir $treeHeight $zRef $sample2 stability[@]
plotVeer $timeDir $treeHeight $zRef $sample1 stability[@]
#------------------------------------------------------------------------------

View File

@ -0,0 +1,101 @@
/*--------------------------------*- 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;
vertices
(
(0 0 0)
(200 0 0)
(200 200 0)
(0 200 0)
(0 0 6000)
(200 0 6000)
(200 200 6000)
(0 200 6000)
);
blocks
(
hex (0 1 2 3 4 5 6 7) (3 3 160) simpleGrading (1 1 80)
);
edges
(
);
boundary
(
top
{
type patch;
faces
(
(4 5 6 7)
);
}
bottom
{
type wall;
faces
(
(0 1 2 3)
);
}
inlet
{
type cyclic;
neighbourPatch outlet;
faces
(
(0 4 7 3)
);
}
outlet
{
type cyclic;
neighbourPatch inlet;
faces
(
(1 2 6 5)
);
}
left
{
type cyclic;
neighbourPatch right;
faces
(
(0 4 5 1)
);
}
right
{
type cyclic;
neighbourPatch left;
faces
(
(3 7 6 2)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //

View File

@ -0,0 +1,104 @@
/*--------------------------------*- 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 buoyantBoussinesqSimpleFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 100000;
deltaT 1;
writeControl timeStep;
writeInterval 50000;
purgeWrite 0;
writeFormat ascii;
writePrecision 16;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable false;
functions
{
ObukhovLength1
{
// Mandatory entries
type ObukhovLength;
libs (fieldFunctionObjects);
// Optional entries
U U;
result1 ObukhovLength;
result2 Ustar;
rhoRef 1.0;
kappa 0.4;
beta 3e-3;
// Optional (inherited) entries
writeControl writeTime;
}
fieldAverage1
{
type fieldAverage;
libs (fieldFunctionObjects);
writeControl writeTime;
fields
(
U
{
mean on;
prime2Mean off;
base time;
}
ObukhovLength
{
mean on;
prime2Mean off;
base time;
}
Ustar
{
mean on;
prime2Mean off;
base time;
}
);
}
#includeFunc "turbulenceFields"
#includeFunc "samples"
}
// ************************************************************************* //

View File

@ -0,0 +1,27 @@
/*--------------------------------*- 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 2;
method hierarchical;
coeffs
{
n (1 1 2);
}
// ************************************************************************* //

View File

@ -0,0 +1,68 @@
/*--------------------------------*- 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;
grad(U) cellMDLimited Gauss linear 1;
}
divSchemes
{
default none;
div(phi,T) bounded Gauss upwind;
div(phi,U) bounded Gauss upwind;
div(phi,k) bounded Gauss upwind;
div(phi,epsilon) bounded Gauss upwind;
div(phi,omega) bounded Gauss upwind;
div((nuEff*dev(T(grad(U))))) Gauss linear;
div((nuEff*dev2(T(grad(U))))) Gauss linear;
}
laplacianSchemes
{
default Gauss linear corrected;
}
interpolationSchemes
{
default linear;
}
snGradSchemes
{
default corrected;
}
fluxRequired
{
default no;
p_rgh;
p_rghPotential;
}
wallDist
{
method meshWave;
}
// ************************************************************************* //

View File

@ -0,0 +1,76 @@
/*--------------------------------*- 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_rgh
{
solver PCG;
preconditioner DIC;
tolerance 1e-15;
relTol 0.1;
maxIter 10;
}
"(U|k|epsilon|omega)"
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-12;
relTol 0.1;
}
T
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-15;
relTol 0.1;
maxIter 10;
}
}
SIMPLE
{
nNonOrthogonalCorrectors 0;
pRefCell 0;
pRefValue 0;
}
relaxationFactors
{
fields
{
p_rgh 0.15;
}
equations
{
U 0.1;
k 0.3;
omega 0.3;
epsilon 0.3;
T 0.015;
}
}
cache
{
grad(U);
grad(T);
}
// ************************************************************************* //

View File

@ -0,0 +1,34 @@
/*--------------------------------*- C++ -*----------------------------------*/
type sets;
libs (sampling);
interpolationScheme cellPoint;
setFormat raw;
writeControl writeTime;
fields
(
T
p_rgh
U
k
ObukhovLength
Ustar
turbulenceProperties:R
);
sets
(
lineZ1
{
type midPoint;
axis z;
start (0 0 0);
end (0 0 6001);
nPoints 200;
}
);
// *********************************************************************** //

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 setFieldsDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
defaultFieldValues
(
volScalarFieldValue qPlant 0
volScalarFieldValue plantCd 0
volScalarFieldValue leafAreaDensity 0
);
regions
(
boxToCell
{
box (0 0 8.60799) (200 200 18.2192);
fieldValues
(
volScalarFieldValue qPlant Q_PLANT
);
}
boxToCell
{
box (0 0 0) (200 200 18.2192);
fieldValues
(
volScalarFieldValue plantCd 0.2
volScalarFieldValue leafAreaDensity 0.14
);
}
patchToFace
{
patch bottom;
fieldValues
(
volScalarFieldValue plantCd 0.2
volScalarFieldValue leafAreaDensity 0.14
);
}
);
// ************************************************************************* //

View File

@ -0,0 +1,9 @@
/*--------------------------------*- C++ -*----------------------------------*/
type turbulenceFields;
libs (fieldFunctionObjects);
field R;
writeControl writeTime;
// *********************************************************************** //