September 6, 2016

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

Wind as per ISO 19901

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 screen.
  --height=H  Maximum height (m). [default: 140.0]
  --speed=S   1h mean wind speed at zr ref. height (m/s).

"""
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.savefig('isowind.png')
    pass

if __name__ == '__main__':
    main()

Probability density function

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

Probability density function

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.png')
    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 flavor of things that come with the territory.

PVP, RF, TI, TLS

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.
  --height=H  Height of structure [default: 140.]
  -i          Plot turbulent intensity.
  -l          Plot turbulent length scale.
  -p          Plot peak velocity pressure.
  -r          Plot terrain roughness.
  --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)
        Lz = Lt * (z / zt)**alpha # m, turbulent length scale
        plt.plot(Lz, z, label='Terrain: %d' %(z0.index(i)),linewidth=2)
    misc()
    plt.xlabel('Turbulent length scale, $L(z)$')
    plt.savefig('tls.png')
    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.png')
    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)
        k1 = 1.0 # recommended turbulence factor, eq. 4.7
        c0 = 1.0 # orography factor for below 3deg gradient, \S 4.3.3
        Iv = k1 / (c0 * np.log(z / i)) # turbulence intensity, \S 4.4
        plt.plot(Iv, z, label='Terrain: %d' %(z0.index(i)),linewidth=2)
    misc()
    plt.xlabel('Turbulence intensity, $I_v(z)$')
    plt.savefig('ti.png')
    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)
        k1 = 1.0 # recommended turbulence factor, eq. 4.7
        c0 = 1.0 # orography factor for below 3deg gradient, \S 4.3.3
        Iv = k1 / (c0 * np.log(z / i)) # turbulence intensity, \S 4.4
        qp = (1.0 + 7.0 * Iv) * 0.5 * cr**2 # in terms of vb
        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.png')
    pass

if __name__ == '__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

  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.)