Utilities

zCFD is bundled with a number of utility functions to streamline specific workflows:

Windfarm Modelling

Overview

Windfarms consist of a set of specific turbine types at geographic locations. zCFD allows users to specify the turbine characteristics and locations, and will automatically create turbine models within an existing mesh and update the local flow conditions to match the thrust coefficient and tip speed ratio for each turbine.

The locations and TRBX turbine definition

zCFD provides a utility function create_trbx_zcfd_data that can be called from within the zCFD environment and Python interactive shell. To start the zCFD command line environment, type:

> source /INSTALL_LOCATION/zCFD-version/bin/activate

The command prompt will now appear with a (zCFD) prefix. To run the create_trbx_zcfd_data utility, start the interactive Python shell:

> (zCFD) > python

The utility takes several arguments, depending on the turbine model to be used. For simple and induction models (where a reference velocity point is used to set the power, thrust and torque for the turbine) the function is used to look up data in a TRBX file format:

import zutil.farm as farm
farm.create_trbx_zcfd_input(case_name='windfarm_name',
                            wind_direction=267.0,
                            reference_wind_speed=10.0,
                            terrain_file=None,
                            report_frequency=100,
                            update_frequency=1,
                            reference_point_offset=1.0,
                            turbine_zone_length_factor=2.5,
                            model='simple',
                            turbine_files=[['xy_file_1.txt','turbine_type_1.trbx'],
                                           ['xy_file_2.txt','turbine_type_2.trbx']])

The utility will read each [<xy_file> , <trbx_file>] pair provided. The terrain file is any file format that can be read by ParaView (STL, VTK) to use as the basis for the ground level at any given location. The default ‘None’ sets ground level to 0. The model uses a reference velocity at a point upstream of the turbine. The point is offset upstream by a factor applied to the rotor diameter. The fluid zone associated with the turbine is a cylinder with diameter equal to the rotor diameter and length given by the turbine zone length factor. The model is either an induction model where the induction factor is automatically calculated from the thrust coefficient and 1D momentum theory, or the simple model where the reference velocity is equal to the velocity at the rotor. The <xy_file> is a 3-column data file in ASCII format with <turbine name>, <easting> and <northing> data such as:

T001 251865.0 646486.0
T002 252240.0 646200.0
...
T015 253308.0 647985.0

The TRBX file is a data file with turbine information provided by manufacturers including a description of the turbine geometry and aerodynamic performance characteristics.

The utility extracts key information from the TRBX file to automatically create a turbine definition (<case_name>_zones.py) with an entry for each turbine in the array:

turb_zone = {
    'FZ_1':{
    'type':'disc',
    'def':'./turbine_vtp/T001-267.0.vtp',
    'thrust coefficient':0.821,
    'thrust coefficient curve':[[0.0,0.0],[1.0,0.0],...,[51.0,0.0],],
    'tip speed ratio':8.79645943005,
    'tip speed ratio curve':[[0.0,0.0],[1.0,0.0],...,[51.0,0.0],],
    'centre':[251865.0,646486.0,338.5],
    'up':[0.0,0.0,1.0],
    'normal':[-0.998629534755,-0.0523359562429,-0.0],
    'inner radius':1.95500004292,
    'outer radius':60.0,
    'reference point':[251745.164456,646479.719685,338.5],
    },
    ...
}

The ‘thrust coefficient’ may be defined as a single value, and if a data point curve is present in the TRBX file this data is also supplied as the tuple array ‘thrust coefficient curve’, where the wind speed in metres per second is the index. If this curve is provided, the utility will interpolate the curve at the reference speed to create the single value otherwise a default value (10 m/s) is used.

The same approach is taken for the ‘tip speed ratio’ single value and curve - which is automatically calculated from the rotor speed array (revolutions per minute) in the TRBX file.

The ‘reference point’ defines the location in the flow domain that is used as the reference value of wind velocity for this turbine. This velocity is used in combination with the thrust coefficient and the tip speed ratio for zCFD to calculate the momentum sources associated with the turbine. By default the reference point is automatically located 1.0 turbine diameters upstream of the disc centre, assuming that the reference wind direction is also the local wind direction. This is easy to modify.

The ‘centre’ is the centre of the disc, which is automatically determined from the nominal hub height in the TRBX file as an offset to the ground height at the specified location. The local ground height is automatically determined from a supplied VTK files, or from a previous solver run. Note that the VTK ground data can be created with a single cycle of the solver, and does not need to include any turbines.

The vertical orientation is defined by the ‘up’ vector - normally this will be the unit vector in the z-direction. The ‘normal’ defines the vector perpendicular to the disc. The inner and outer radii are based on the TRBX definition of the size of the disc. No account is made of the hub or tower geometry.

For the blade element model the TRBX file is replaced with a dictionary that contains information about the blades and aerofoil sections used on the rotor, plus information that can be specified to inform the controller how to adjust the rotation rate and blade pitch each cycle:

farm.create_trbx_zcfd_input(
    case_name='windfarm_name',
    wind_direction=280.0,
    reference_wind_speed=8.0,
    terrain_file=None,
    report_frequency=100,
    update_frequency=1,
    turbine_zone_length_factor=2.5,
    model='blade element theory',
    turbine_files=[['data/xy_file_1.txt', bet_dict_1],
                   ['data/xy_file_2.txt', bet_dict_2]],
)

The xy_file is as above and the bet_dict is a Python dictionary with details of the rotor and blade, including sectional aerodynamic characteristics and dimensions:

bet_dict = {
    'name':'T001',
    'status': 'on',
    'number of blades': 3,
    'rotation direction': 'clockwise',
    'hub height': 70.0,
    'inner radius': 4.0,
    'outer radius': 40.0,
    'rated power': 2000000.0,
    'tip speed limit': 80.0,
    'verbose': False,
    'tilt': 6.0,
    'yaw': 0.0,
    'auto yaw': True,
    'cut in speed': 4.0,
    'cut out speed': 25.0,
    'number of segments': 12,
    'up': [0.0, 0.0, 1.0],
    'aerofoil profile':
        {'source': 'http://airfoiltools.com/airfoil/details?airfoil=n63415-il',
         'upper surface': [[0.0, 0.0], [0.003, 0.01287], ..., [0.95028, 0.00931], [1.0, 0.0]],
         'lower surface': [[0.0, 0.0], [0.007, -0.01087], ..., [0.94972, 0.00333], [1.0, 0.0]]},
    'aerofoil cl': [[-90.0, -0.0],[-49.869, -0.202], ..., [49.907, 1.823], [90.0, 0.0]],
    'aerofoil cd': [[-90.0, 2.0], [-70.099, 1.796], ..., [60.11, 1.645], [90.087, 1.996]],
    'blade chord': [[0.0, 0.1], [0.1, 0.1], ..., [0.7, 0.034], [1.0, 0.018]],
    'blade twist': [[0.0, 21.924], [0.1, 21.924], ..., [0.7, 2.461], [1.0, 0.023]],
    'blade pitch range': [-10.0, 10.0],
    'blade pitch': 1.0,
    'blade pitch tol': 0.01,
    'blade pitch step': 0.1,
    'mean blade material density': 200.0,
    'blade pitch tol': 0.01,
    'dt': 0.2,
    'inertia': True,
    'damage ti': 0.15,
    'damage speed': 10.0,
    'friction loss': 0.01,
    'thrust factor': 1.10,
    'tip loss correction': 'acos-fit',
    'rpm ramp':[[4.0,6.0],[13.0,16.0]],
}

The ‘name’ key is used to tag results associated with the turbine.

The ‘status’ key indicates whether the turbine should be active or not. This makes it straightforward to turn a turbine off without regenerating the model. The ‘number of blades’ is typically 3 (which is also the default value) but some turbines have 2 blades. The ‘rotation direction’ is typically clockwise (when viewed from the front) for modern turbines, but the keyword values allow ‘clockwise’ or ‘anticlockwise’. This alters the rotation of the wake.

The ‘hub height’, ‘inner radius’ and ‘outer radius’ define the vertical offset from the ground in metres, and the rotor dimensions respectively as above. The ‘rated power’ (Watts) of the turbine defines an operating limit of power production that will not be exceeded.

The ‘tip speed limit’ is an environmental limit on the rate or rotor rotation, and is typically set to 80m/s. Note that the tip speed is based on the absolute rotation rate of the rotor, and the actual tip speed relative to the air may be increased by angular induction.

The model can produce various diagnostics. Setting ‘verbose’ to True will output these.

The rotor ‘tilt’ is an upward rotation of the entire rotor, to allow for tower clearance under blade deflection. The value is typically a few degrees. The ‘yaw’ is a rotation angle in degrees from the prescribed wind direction, with the option to specify ‘auto yaw’ which will automatically over-ride the ‘yaw’ value with a rotation into the local (disc average) wind direction. A warning is issued if ‘auto yaw’ requests a ‘yaw’ value exceeding 10 degrees in either direction.

The manufacturer typically specifies a ‘cut in speed’ and ‘cut out speed’ range for the local reference wind speed (for blade elements this is the average windspeed on the rotor disc). The rotor is automatically stowed outside these limits.

The ‘number of segments’ defines the resolution of the actuator disc discretisation (not the computational domain mesh) with a default value of 12.

The ‘aerofoil profile’ is used to integrate the blade unit sectional area between the ‘upper surface’ and ‘lower surface’ over the range [0,1]. Once integrated, the value is stored in the turbine zone dictionary. The optional ‘source’ can be provided as a reference for the data. The sectional area multiplied by local radius is integrated over the length of each blade and multiplied by the ‘mean blade material density’ (kilograms per cubic metre, typically 200.0) to provide a value for the ‘rotor moment’.

The turbine controller assumes that the physical time between iterations is ‘dt’. This does not have to be equal to the CFD timestep. If ‘inertia’ is True then the available torque is divided by the ‘rotor moment’ to determine the rotor acceleration. Otherwise the rotor is accelerated with a 10% increase in rotational rate capped per iteration. When the rotor is positively accelerated, half of the available torque is used to accelerate the rotor and half is supplied to the generator (unless doing so would exceed the rated power). This is an arbitrary value and it is expected that users will modify this behaviour when modelling more sophisticated control laws.

The aerodynamic characteristics of the aerofoil section are critical to the blade element theory model and are not calculated automatically. These should be specified as points on the ‘aerofoil cl’ (local lift coefficient) and ‘aerofoil cd’ (local drag coefficient) curves over a range of angles of attack in degrees covering at least (-30,30). Values beyond this range will be extrapolated as needed.

The ‘blade chord’ is the local chord length in metres used to scale the aerofoil along the normalised blade length [0,1]. The ‘blade twist’ is the twist angle in degrees of the blade over the normalised blade length [0,1].

The turbine controller will automatically optimise the power production (up to the limit) by pitching the blades over the range ‘blade pitch range’ in degrees using a Golden Section Search to a tolerance of ‘blade pitch tol’ degrees. All blades are pitched equally and the optimum pitch value accounts for the contribution of all blades averaged over the disc. The ‘blade pitch’ can be set initially and a step size limit ‘blade pitch step’ can be set to the range over which the blade pitch is optimised to within ‘blade pitch tol’ on each iteration. Typically the blades are pitched to their maximum setting on startup and then the pitch value is reduced as the rotor accelerates. A negative value of ‘blade pitch’ may be achieved at high rates of rotation as the local angle of attack of the blade will still be positive. A warning is issued if the aerofoil angle of attack is less than zero anywhere on the rotor.

The turbine controller will determine whether (anywhere on the disc) the ‘damage ti’ (turbulence intensity) and ‘damage speed’ are simultaneously exceeded. If so, a warning is issued and the rotor is stowed. It is intended that this function may be modified by users to provide more sophisticated controller rules and wind farm behaviour.

To account for frictional losses in the system, 1% (or other value of ‘friction loss’) is applied to the rate of rotor rotation each iteration.

The ‘blade element theory’ fluid zones follow the ‘simple’ and ‘induction’ model example above, with inclusion of the ‘blade element theory’ parameters.

A ‘thrust factor’ can be applied to the rotor to include an approximation for the effect of the tower and the hub on the obstruction to airflow past the turbine.

The blade element theory model assumed 2D flow over each aerofoil section, and a 3D correction can be applied at the rotor tips. Options for the ‘tip loss correction’ keyword are ‘elliptic’ (which is not recommended for most modern turbines), ‘acos-fit’, ‘acos shift-fit’ and ‘f-fit’. Definitions for these functions are provided in the Reference Guide.

In addition to the tip speed limit imposed by environmental considerations, some manufacturers also apply a calibrated ramp factor to the revolutions per minute that the rotor can achieve for any given freestream hub-height wind speed. This is to avoid overspeeds at low wind speeds. The ‘rpm ramp’ keyword takes a pair of points [[cut in speed, lowest rpm],[full speed, rpm at full speed]] defining the ends of the linear limit curve applied.

For all models, the turbine_vtk/<turbine>.vtp file defining the fluid zone for each turbine is also automatically created.

The utility also automatically creates a set of monitor points for each turbine, all in a single file (<case_name>_probes.py):

turb_probe = {
              'report' : {
                           'frequency' : 100,
                           'monitor' : {
                                         'MR_1' : {
                                         'name' :'probe1@MHH@87',
                                         'point' : [251865.0,646486.0,338.5],
                                         'variables' : ['V', 'ti'],
                                                  },
                                         ...
                                       }
                         }
             }

The ‘frequency’ is the number of solver cycles between outputs, and the ‘monitor’ defines the name of the probe using the WindFarmer standard notation.

Because the zone and probe files are automatically created, the following lines must be added to the end of the standard zCFD parameter definition file <case_name>.py to insert the data:

z = zutil.get_zone_info('<case_name>_zones')
for key,value in z.turb_zone.items():
parameters[key]=value

p = zutil.get_zone_info('<case_name>_probes')
for key,value in p.turb_probe.items():
parameters[key]=value

When run, zCFD will include the probe data in the <case_name>_report.csv file. Note that this utility may take a few seconds to run, especially if there are large numbers of turbines, or points on the mesh boundary.

Writing WindFarmer Data Files

In order to export the data from a zCFD run in a format that can be read by WindFarmer, we provide the utility write_windfarmer_data, with usage:

import zutil.farm as farm
farm.write_windfarmer_data(case_name='windfarm_name',
                               num_processes=48,
                               up = [0,0,1])

The up vector is used to check that the orientation expected by WindFarmer is that same as the orientation used in the simulation. In most cases this will be the z-axis.

The utility will output the probe information plus additional fields, calculated automatically.

Propellor Modelling

Overview

Propellors are similar to wind turbines except that they are designed to produce thrust rather than torque from lift. zCFD includes a blade element fluid zone model based on an actuator disc, for propellors:

propellor_zones = {
    'FZ_1': {
        'name': 'R01',
        'type': 'disc',
        'model': 'blade element propellor',
        'number of blades': 2,
        'status': 'on',
        'outer radius': 0.175,
        'inner radius': 0.039025,
        'def': './propellor_vtp/R01.vtp'
        'rotation direction': 'clockwise',
        'verbose': True,
        'update frequency': 1,
        'centre': [0.0, 0.0, 0.0],
        'reference plane': True,
        'normal': [0.0, 0.0, 1.0],
        'up': [1.0, 0.0, 0.0],
        'blade chord': [[0.2333, 0.1468], [0.2667, 0.1681], ... , [1.0, 0.0581]],
        'blade twist': [[0.2333, 28.6202], [0.2667, 25.5228], ... , [1.0, 0.0]],
        'aerofoil cd': [[-60.0, 1.4552], [-50.0, 1.0992], ... , [60.0, 2.3542]],
        'aerofoil cl': [[-60.0, -0.7978], [-50.0, -0.8195], ... , [60.0, 1.3818]],
        'number of segments': 24,
        'tip loss correction': 'elliptic',
        'tip loss correction radius': 0.8,
        }
    }

The ‘name’ key is used to tag results associated with the propellor.

The ‘type’ defines the fluid zone type for propellors.

The ‘model’ should be set to ‘blade element propellor’.

The ‘number of blades’ can be any positive integer - typically 2, 3 or 4.

The ‘status’ can be set to on or off, as a quick way to disable a propellor without re-generating the model.

The ‘inner radius’ and ‘outer radius’ refer to the rotor disc.

The ‘def’ is a path to a file containing the VTP file defining the disc fluid zone geometry. This can easily be created interactively in ParaView or else using the ParaView Simple library in Python:

line = Line()
line.Resolution = 10
line.Point1 = (0.0, 0.0, -0.5 * 0.175)
line.Point2 = (0.0, 0.0,  0.5 * 0.175)
tube = Tube(Input=line)
tube.NumberofSides = 128
tube.Radius = 0.175
transform = Transform(Input=tube)
transform.Transform = "Transform"
transform.Transform.Rotate = [0.0, 0.0, 0.0]
transform.Transform.Translate = [0.0, 0.0, 0.0]
writer = CreateWriter('./propellor_vtp/R01.vtp')
writer.Input = transform
writer.UpdatePipeline()

where the above code block would produce a fluid zone for a propellor with normal [0.0,0.0,1.0]. Note that general use of the ParaView rotate command should follow the (default) ZXY convention, for which the following is a reasonable estimate:

def rotation_matrix(axis, theta):
    # Rotation matrix for counterclockwise rotation about axis by theta (radians)
    axis = np.asarray(axis)
    axis = axis / math.sqrt(np.dot(axis, axis))
    a = math.cos(theta / 2.0)
    b, c, d = -axis * math.sin(theta / 2.0)
    aa, bb, cc, dd = a * a, b * b, c * c, d * d
    bc, ad, ac, ab, bd, cd = b * c, a * d, a * c, a * b, b * d, c * d
    return np.array([[aa + bb - cc - dd, 2 * (bc + ad), 2 * (bd - ac)],
                     [2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)],
                     [2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]])

def trans_rot(v, vec, ang):
    v2 = np.dot(rotation_matrix(vec, ang), v)
    return [v2[0], v2[1], v2[2]]

# ParaView rotation emulator - note that ParaView (i.e. VTK) uses the order ZXY for rotation.
def pv_rotate(pt,args):
    pt1 = trans_rot(pt,[0.0,0.0,1.0],math.radians(args[2]))
    x1 = trans_rot([1.0,0.0,0.0],[0.0,0.0,1.0],math.radians(args[2]))
    y1 = trans_rot([0.0,1.0,0.0],[0.0,0.0,1.0],math.radians(args[2]))
    pt2 = trans_rot(pt1,x1,math.radians(args[0]))
    y2 = trans_rot(y1,x1,math.radians(args[0]))
    pt3 = trans_rot(pt2,y2,math.radians(args[1]))
    return pt3

The ‘rotation direction’ can be clockwise or anticlockwise. zCFD will apply the necessary orientations to the forces and wake rotation.

The ‘verbose’ keyword (True or False) specifies the level of output from the model.

The ‘update frequency’ for propellors should be 1 - meaning that the model is evaluated on every solver cycle.

The ‘centre’ is the centre of the disc rotor in 3D, and the ‘normal’ defines the orientation.

The ‘up’ vector is used to define $theta_0$ on the rotor: the radial vector from the centre to the rotor tip at $theta_0$ is defined by the cross product of the up vector and disc normal. All forces on the disc are evaluated by incrementing $theta$ over the range $[0,2pi]$.

The ‘blade chord’ and ‘blade twist’ are each defined as a list of points over the domain [0,1]. Note that the twist is measured in degrees and the chord is scaled by the outer radius in zCFD, so the input has to be scaled for a radius of 1.

The ‘aerofoil cl’ and ‘aerofoil cd’ curves are lists of points defining the aerofoil characteristics. zCFD interpolates values between these points and will extrapolate if necessary.

The ‘number of segments’ used in the model should be chosen to match the mesh spacing in the CFD mesh around the propellor. zCFD does not generate a new mesh in order to apply the propellor model.

A ‘tip loss correction’ can be applied to the blade element model, and the only option is ‘elliptic’. The ‘tip loss correction radius’ is the radius from which the tip correction is applied and has a default value of 0.8.