# Wind

The following offers a way to visualize wind profile for a given Uo and t for a range of heights up to a maximum of h.

## Based on ISO 19901-1:2005

Plot code is as follows.

#!/usr/bin/env python
# -*- coding: UTF-8 -*-

"""Wind speed plots based on ISO 19901-1:2005.
2016 ckunte.

Usage: isowind.py --speed=S [--height=H]
isowind.py --help
isowind.py --version

Options:
-h, --help  Show this help.
--speed=S   1h mean wind speed at zr ref. height (m/s).
--height=H  Maximum height (m). [default: 140.0]

"""
import numpy as np
import matplotlib.pyplot as plt
from docopt import docopt

def main():
args = docopt(__doc__, version='ISO 19901-1 wind speed plots, version: 0.1')
Uo = float(args['--speed']) # m/s 1h mean wind speed at zr
h = float(args['--height']) # m, max. height (e.g., a flare tower)

t = [3., 5., 60., 600., 3600.] # in seconds
t0 = 3600.0 # s (=> 1h)

zr = 10.0 # m, ref. height
z = np.linspace(0.1, h)

Iuz = 0.06 * (1 + 0.043 * Uo) * (z / zr)**(-0.22)
C = 0.0573 * (1 + 0.15 * Uo)**(0.5) # for zr = 10m
Uz = Uo * (1 + C * np.log(z / zr))

for i in t:
u = Uz * (1 - 0.41 * Iuz * np.log(i / t0))
lbl = 'u(z,t) -- %0.0fs' %i
plt.plot(u, z, label=lbl, linewidth=2)

# ref: https://goo.gl/CJgoyu
plt.title('ISO Wind profile (Uo=%.1fm/s at %.1fm reference elevation)' %(Uo,zr))
plt.ylabel('Height, z (m)')
plt.xlabel('Wind speed (m/s)')
plt.legend(loc=0)
plt.grid(True)
plt.savefig('isowind.svg')
pass

if __name__ == '__main__':
main()


## Probability density function

The following code lets me regenerate wind pdf function plots shown below.

Plot code:

#!/usr/bin/env python
# -*- coding: UTF-8 -*-
"""
pdf.py: 2016 ckunte
"""
import numpy as np
import matplotlib.pyplot as plt

def main():
# Weibull distribution: probability distribution function
# (for extra-tropical storms)
## -- begin inputs --
ls = [1.0, 1.12] # Mean / nominal speed
lf = [0.78, 0.94] # Mean / nominal force
## -- end inputs ----

k = 5. # Safety index

x = np.linspace(0.,2.5)

for i, j in zip(ls, lf):
fs = (k / i) * ((x / i)**(k - 1.)) * np.exp(-(x / i)**k)
plt.plot(x, fs, linewidth=2, label='Mean/nominal speed: %0.2f' %i)
ff = (k / j) * ((x / j)**(k - 1.)) * np.exp(-(x / j)**k)
plt.plot(x, ff, linewidth=2, label='Mean/nominal force: %0.2f' %j)

plt.ylabel('Probability density function (pdf) of wind speed')
plt.xlabel('Normalised wind speed')
plt.grid(True)
plt.legend(loc=0)
plt.savefig('pdf.svg')
plt.show()
pass

if __name__ == '__main__':
main()


## Based on EN 1991-1-4:2005

The generic nature of EN’s applicability makes for a complicated recipe. The number of things one needs to determine before one even gets to wind velocity is staggering. Here’s a flavour of things that come with the territory.

Here’s the code for these plots.1

#!/usr/bin/env python
# -*- coding: UTF-8 -*-

"""Wind action plots based on EN 1991-1-4:2005.
2016 ckunte.

Usage: enwind.py ( -i | -l | -p | -r ) [--height=H]
enwind.py --help
enwind.py --version

Options:
-h, --help  Show this help screen.
-i          Plot turbulent intensity.
-l          Plot turbulent length scale.
-p          Plot peak velocity pressure.
-r          Plot terrain roughness.
--height=H  Height of structure [default: 140.]
--version   Show version.

"""
import numpy as np
import matplotlib.pyplot as plt
from docopt import docopt

z0 = [0.003, 0.01, 0.05, 0.3, 1.0]  # m, Roughness length
zmin = [1.0, 1.0, 2.0, 5.0, 10.] # m, reference length scale
zt = 200. # m, reference height, Annex B.1
Lt = 300. # m, reference length scale, Annex B.1
rho = 1.25 # kg/m^3, air density, \S 4.5

args = docopt(__doc__, version='EN wind action plots, version: 0.1')
h = float(args['--height'])

def misc():
plt.legend(loc=0)
plt.grid(True)
plt.ylabel('Height, z (m)')
pass

def tls():
# Wind turbulence, Annex B, EN 1991-1-4:2005
for i, j in zip(z0, zmin):
alpha = 0.67 + 0.05 * np.log(i)
z = np.linspace(j, h)
# turbulent length scale (m)
Lz = Lt * (z / zt)**alpha
plt.plot(Lz, z, label='Terrain: %d' %(z0.index(i)),linewidth=2)
misc()
plt.xlabel('Turbulent length scale, $L(z)$')
plt.savefig('tls.svg')
pass

def tr():
# Terrain roughness, Clause 4.3.2, EN 1991-1-4:2005
for i, j in zip(z0, zmin):
z = np.linspace(j, h)
kr = 0.19 * (i / z0[2])**0.07
cr = kr * np.log(z / i)
plt.plot(cr, z, label='Terrain: %d' %(z0.index(i)),linewidth=2)
misc()
plt.xlabel('Roughness factor, $c_r(z) = v_m(z) / v_b$')
plt.savefig('rf.svg')
pass

def ti():
# Turbulence intensity, \S 4.4, EN 1991-1-4:2005
for i, j in zip(z0, zmin):
z = np.linspace(j, h)
# recommended turbulence factor, eq. 4.7
k1 = 1.0
# orography factor for below 3deg gradient, \S 4.3.3
c0 = 1.0
# turbulence intensity, \S 4.4
Iv = k1 / (c0 * np.log(z / i))
plt.plot(Iv, z, label='Terrain: %d' %(z0.index(i)),linewidth=2)
misc()
plt.xlabel('Turbulence intensity, $I_v(z)$')
plt.savefig('ti.svg')
pass

def pvp():
# Peak velocity pressure, \S 4.5, EN 1991-1-4:2005
for i, j in zip(z0, zmin):
z = np.linspace(j, h)
kr = 0.19 * (i / z0[2])**0.07
cr = kr * np.log(z / i)
# recommended turbulence factor, eq. 4.7
k1 = 1.0
# orography factor for below 3deg gradient, \S 4.3.3
c0 = 1.0
# turbulence intensity, \S 4.4
Iv = k1 / (c0 * np.log(z / i))
# in terms of vb
qp = (1.0 + 7.0 * Iv) * 0.5 * cr**2
plt.plot(qp, z, label='Terrain: %d' %(z0.index(i)),linewidth=2)
misc()
plt.xlabel('Peak velocity pressure $q_p(z)$ in terms of $v_b$')
plt.savefig('pvp.svg')
pass

def main():
if h <= zt:
if args['-p']:
pvp()
elif args['-l']:
tls()
elif args['-r']:
tr()
elif args['-i']:
ti()
else:
print("No option was selected. For help, try: python enwind.py -h")
else:
print("height >", zt, "m for which these plots are not appropriate.")
pass

if __name__ == '__main__':
main()


1. The arguments in square brackets are optional, while at least one in those parenthesis is mandatory. (Now using the smart docopt, which makes coding for arguments super simple.)