openfoam/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/plot
2020-06-16 10:51:44 +01:00

340 lines
11 KiB
Bash
Executable File
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

#!/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/atm-HargreavesWright-2007/Ux-HW-RH-Fig6a" \
u 1:2 t "RH" w p ps 2 pt 6 lc rgb "#000000", \
"system/atm-HargreavesWright-2007/Ux-HW-RH-Fig6a" \
u 1:2 t "HW, x=2500m" w p ps 1 pt 5 lc rgb "#E69F00", \
"system/atm-HargreavesWright-2007/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/atm-HargreavesWright-2007/Ux-HW-RH-Fig6a" \
u 1:2 t "RH" w p ps 2 pt 6 lc rgb "#000000", \
"system/atm-HargreavesWright-2007/Ux-HW-RH-Fig6a" \
u 1:2 t "HW, x=2500m" w p ps 1 pt 5 lc rgb "#E69F00", \
"system/atm-HargreavesWright-2007/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/atm-HargreavesWright-2007/Ux-HW-RH-Fig6a" \
u 1:2 t "RH" w p ps 2 pt 6 lc rgb "#000000", \
"system/atm-HargreavesWright-2007/Ux-HW-RH-Fig6a" \
u 1:2 t "HW, x=2500m" w p ps 1 pt 5 lc rgb "#E69F00", \
"system/atm-HargreavesWright-2007/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/atm-HargreavesWright-2007/k-RH-Fig6b" \
u 1:2 t "RH" w p ps 2 pt 6 lc rgb "#000000", \
"system/atm-HargreavesWright-2007/k-HW-Fig6b-2500" \
u 1:2 t "HW, x=2500m" w p ps 1 pt 5 lc rgb "#E69F00", \
"system/atm-HargreavesWright-2007/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/atm-HargreavesWright-2007/epsilon-HW-RH-Fig6c" \
u 1:2 t "RH" w p ps 2 pt 6 lc rgb "#000000", \
"system/atm-HargreavesWright-2007/epsilon-HW-RH-Fig6c" \
u 1:2 t "HW, x=2500m" w p ps 1 pt 5 lc rgb "#E69F00", \
"system/atm-HargreavesWright-2007/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/atm-HargreavesWright-2007/mut-RH-Fig6d" \
u 1:2 t "RH" w p ps 2 pt 6 lc rgb "#000000", \
"system/atm-HargreavesWright-2007/mut-HW-Fig6d-2500" \
u 1:2 t "HW, x=2500m" w p ps 1 pt 5 lc rgb "#E69F00", \
"system/atm-HargreavesWright-2007/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
#------------------------------------------------------------------------------