Load Models
1. Overview of Load Effects
The following physical effects contribute to system loads:

weight and inertia, governed by line mass, including any pipe contents and external wrapping

hydrostatic forces, dependent on pressure gradients

hydrodynamic forces, dependent on wave, current and structure motions

forced motion of line, dependent on vessel motions

aerodynamic forces, dependent on wind and structure motions
Since the analysis is based on a beam model with pressureindependent cross section properties, the hydrostatic pressure level does not affect the static or dynamic equilibrium condition, hence the pressure gradients are sufficient.
All loads are applied either as distributed line loads (force/length) or as point loads.
The following sections give an outline of essential features of the implemented load models.
With regard to surfacepiercing structures, two classes of hydrodynamic load models are used:

Hydrostatic and hydrodynamic loads are applied on a length equal to the submerged length of the component axis, i.e. the volume forces are concentrated in the pipe axis. Referring to the user manual this applies to the load model:
Load based on Morisons generalized equation (MORI)
. Hydrodynamic loads model are described in section. Section 5. 
Hydrostatic and hydrodynamic forces are distributed according to the actual diameter, or shape of the buoyancy volume. Thus, the cross section can be partly submerged, and structures such as floating hoses and fish cages can be modelled. Referring to the user manual this applies to the load model:
As MORI, but improved by taking into account partially submerged crosssection (MORP)
. Hydrodynamic load models are described in section Section 8.
For submerged structures, and also for surfacepiercing normal pipes and cables, the first class is used, and generally this cross section is referred to in the subsequent discussion. The second class is referred to as a "partly submerged" cross section.
2. Hydrostatic Pressure Effects
Hydrostatic pressure effects on marine risers have been discussed by several authors in the open literature, e.g.Morgan and Peret (1974), Sparks (1984) and Mathisen and Bergan (1986). This discussion normally involves the following notions that need to be defined.

effective tension: force in the pipe that affects stability. This force is relevant for governing the shape of cables and pipes, including buckling analysis, and is used for calculation of geometric stiffness in the finite element method.

axial stress resultant: force found by integrating normal stresses over the cross section. In the presence of external hydrostatic pressure, this force is different from the effective tension.
In principle, two different approaches are possible when analysing a pipe or other slender structure subjected to external and internal hydrostatic pressure.

Volume force model

Use of vertical, conservative forces to represent hydrostatic effects. These forces will be in equilibrium with the effective tension, which means that axial stresses need not necessarily be calculated during an iteration for equilibrium. The theoretical foundation for this way of modelling hydrostatic forces is given by Sparks (1984).


Pressure force model

Calculation of nonconservative hydrostatic forces by considering pressure on the deformed pipe geometry. These forces are in equilibrium with the axial stress resultants. Effective tension needed for calculation of the stiffness matrix must, however, be found by introducing artificial pressure forces in axial direction. This method is described by Mathisen and Bergan (1986).

Both methods are in principle correct, and all nonlinear geometric and material effects can be modelled within both concepts. The choice of method should therefore be made after considering use of computer time and complexity in formulation of the alternatives. For a structure where hydrostatic pressure does not affect the volume or cross section properties the first approach is correct, and it is computationally much simpler. The conclusion is obvious: in the present program hydrostatic effects will be modelled by conservative, vertical forces. No hydrostatic force variation for caused by deformations will be accounted, except for possible increase or decrease of submerged volume and wetted pipe surface at the sea surface.
Hydrostatic pressure is therefore treated in terms of effective weight and effective tension, defined as:
where

\(\mathrm {T\,\qquad}\) : effective tension

\(\mathrm {T_{p}\!\qquad}\) : tension in pipe wall, i.e. resulting force from normal stresses

\(\mathrm {A_{e},A_{i}}\) : external/internal cross sectional

\(\mathrm {P_{e},P_{i}\:}\) : external/internal hydrostatic pressure

\(\mathrm {w\:\qquad}\) : effective weight per unit length, i.e. submerged weight of pipe including content

\(\mathrm {m_{p}\!\!\qquad}\) : mass of pipe per unit length

\(\mathrm {\rho _{i}\,\qquad}\) : density of internal fluid

\(\mathrm {\rho \;\,\qquad}\) : water density

\(\mathrm {g\;\,\qquad}\) : acceleration of gravity

\(\mathrm {v_{i}\;\!\qquad}\) : velocity of internal fluid flow
The basic axial force output from static and dynamic analysis is the effective tension. The axial pipe wall force can, if required, be calculated by post processing according to Equation (2). The last term in Equation (2), which results from a constant velocity internal fluid flow, is further discussed in Section 7.
3. Wave Potential
3.1. Regular Waves
Regular waves are modelled by two alternative methods

Airy linear wave theory

Stokes’ 5th order wave theory
Stokes’ theory is outlined by Skjelbreia and Hendrickson (1960). Both theories describe regular, long crested waves propagating in an arbitrary direction. The main difference between the theories is that the boundary condition pressure (\(\mathrm {p=0.0}\) ) is valid at mean water level for the linear theory, but is approximately correct at the wave surface for Stokes’ theory. This implies that the linear wave theory is strictly valid for infinitesimal wave amplitudes only, while Stokes’ theory describes a wave of finite height. Waveinduced velocities and acceleration in a wave crest can therefore consistently be described by Stokes’ theory only. The wave profile defined by the two theories for waves of identical lengths and heights are shown in Figure 1.
The wave potential \(\mathrm {\displaystyle \phi _{0}}\) for a regular wave according to Airy’s theory can be expressed as follows:
where

\(\mathrm {\zeta_{a}}\) : wave amplitude

\(\mathrm {g\;\,}\) : acceleration of gravity

\(\mathrm {k\;\,}\) : wave number

\(\mathrm {\beta \;\,}\) : direction of wave propagation. (\(\mathrm {\beta =0}\) corresponds to wave propagation along the positive Xaxis)

\(\mathrm {\psi_{\zeta}}\) : a phase angle lag
\(\mathrm {C_\text{1}}\) is given by:
where \(\mathrm {D}\) is the water depth.
In deep water \(\mathrm {C_{1}}\) can be approximated by
We then obtain the following relations for the particle velocities and acceleration in the undisturbed wave field:
where
Using the deep water approximation we have:
Taking into account finite water depth we have:
The surface elevation is given by:
where \(\mathrm {\phi =\psi}\) . The forward phase shift,
\(\mathrm {\phi }\) , is denoted phase angle
in the
program. Similarly the linearized dynamic pressure is given by:
The wave potential described here is, as previously stated, valid to mean water level only. Prediction of velocities and acceleration in the wave crest is, however, important. The present program has four different methods for handling this problem for linear waves. These are illustrated in Figure 2.
Method 1 is the consistent use of linear wave theory assuming that the water level always remains at mean water level, referring to the figure Figure 2.
This method is included in order to have a reference point, particularly for comparing results with frequency domain solutions.
Method 2 is based on the assumption that the potential is correct at some distance from the surface, and that calculated surface values are correct at the instantaneous surface. Consequently, the potential profile is deformed by stretching and compression according to the following formula (Wheeler 1970):
where

\(\mathrm {Z'\quad }\)  modified Zcoordinate to be introduced in place of Z in Equation (3) and Equation (11)

\(\mathrm {\zeta^\textrm{}\quad \,\,}\) actual surface elevation

\(\mathrm {D^\textrm{}\quad }\)  water depth
Method 3 is a parallel move of potential, assuming that the potential itself is correctly described by linear theory, but valid from the instantaneous sea surface. Note that this strategy does not include a redefinition of the potential from variation of effective water depth.
Method 4 allows the user to apply the wave kinematics up to the instantaneous free surface without stretching: for Z levels above mean water level, results from Z = 0 are used.
Using regular wave models, waveinduced velocities and acceleration are calculated at every node at every time step during the time integration. No storage of such data is carried out. During this process, updating of node position can easily be considered simply by calculating the wave potential at the instantaneous structural position determined by the static position plus dynamic displacement.
The program offers this possibility in addition to the traditional method using the static node positions as reference for calculation of wave kinematics. These alternatives are illustrated in Figure 3.
If a major part of the structural motion is due to linear, waveinduced support vessel motion, it is recommended to use wave kinematics at the static position. Otherwise, the forced terminal point motions and the wave kinematics will be inconsistent.
3.2. Irregular Waves
The present section describes the modelling of an undisturbed Airy wave field which is normally used as the basis for wave loading on slender structures. For modelling of disturbed wave kinematics near the support vessel, see section [Wave_Kinematics_Including_Diffraction_Effects].
An irregular sea state is described as a sum of two wave spectra: a wind sea contribution and a swell contribution:
in which \(\mathrm {S_{\zeta,{1}}}\) and \(\mathrm {S_{\zeta,{2}}}\) describe the frequency distribution of wind sea and swell, respectively. Several standard wave spectra are included, such as PiersonMoscowitz, JONSWAP, etc., as well as a general numeric description which can be used for description of measured data or regular waves (see Input to INPMOD in the User Manual).
\(\mathrm {\phi _{1}}\) and \(\mathrm {\phi _{2}}\) describe the directionality of the waves. Unidirectional waves and cosine spreading functions are included.
\(\mathrm {\beta }\) is the direction of wave propagation.
The spectrum directionality parameters satisfy the relations:
In which \(\mathrm {\sigma _{\zeta}^2}\) is the variance of the surface elevation.
In order to generate time series of surface elevation, water particle velocities and acceleration, the shortcrested irregular sea is discretized into a set of harmonic components. Thus, in complex notation the surface elevation is expressed by:
The random phase angles,\(\mathrm {\phi _\textrm{jk}}\) , are sampled from a uniform distribution over \(\mathrm {[\pi ,\pi ]}\) . The positiondependent phase angle is:
The surface elevation can be expressed as:
The velocity and acceleration components are derived from the surface elevation components.
Adding the harmonic components to obtain time series is performed by a Fourier Transformation algorithm. This approach implies a set of relations between time increment,\(\mathrm {\Delta t}\) , frequency increment, \(\mathrm {\Delta \omega }\) , number of time steps, \(\mathrm {N_t}\) and number of components, \(\mathrm {N_\omega }\).
Thus the duration of the time series is limited to:
The number of nonzero components,\(\mathrm {n_\omega }\) , is, however, usually much less than \(\mathrm {N_\omega }\) due to the limited range of wave spectra which is of interest. Typically only wave periods larger than \(\mathrm {t_\textrm{lim}\cong23}\)s. are of interest. With a sampling interval, \(\mathrm {\Delta t=0.5}\)s, the number of nonzero components is:
The frequency increment and the number of components are chosen so that the time increment and the number of time steps coincide with the wave frequency vessel motion simulation parameters, \(\mathrm {\Delta t_{t\textrm{HF}}}\) and \(\mathrm {N_{t\textrm{HF}}}\) , see section [HF_Vessel_Motion_Model].
In order to increase the duration of the simulation period, several sequences, each of duration \(\mathrm {T}\) , can be generated. For each new sequence a new set of random phase angles,\(\mathrm {\psi_{\zeta\,\textrm{jk}}}\) , is sampled.
Velocities, acceleration and surface elevation are calculated for predefined positions using an FFT algorithm. Two different strategies are possible if changed node positions should be considered in a dynamic analysis involving irregular seas.

Calculation of velocities and acceleration by FFT in a four dimensional grid and interpolation during dynamic analysis to obtain values at wanted positions and time steps.

Calculation of velocities and acceleration by using Fourier series and direct calculation of all harmonic components with traditional sine and cosine functions for every wave component at wanted positions and time increments.
The first strategy has been implemented.
In the present program, velocities and acceleration are calculated for predefined points that coincide with nodes in the finite element model at static equilibrium position, but not necessarily covering all nodes. This is feasible since potential profiles are smooth, which makes interpolation of intermediate values an attractive alternative to calculation and storage for every node. In addition, this interpolation is necessary in any case if stretching of the potential profiles is performed in order to describe water particle motions close to the sea surface.
Three options are available with respect to wave particle motion modelling, and, hence, wave force modelling.

Calculate forces to still water level, applied on the part of the riser that is submerged in static equilibrium position.

Calculate forces to still water level, update line submergence according to vertical motion, see method 1 in figure Figure 2.

Calculate forces to actual, instantaneous water level by stretching/shrinking the zcoordinate, see method 2 in figure Figure 2 (Wheeler’s method).

Calculate forces to actual, instantaneous water level by shifting the zaxis upwards or downwards, following the actual wave surface, see method 3 in figure Figure 2.
Since the waveinduced loads are most important for the upper part, and even there may be of minor importance for slender lines, such as anchor lines, options for calculating wave kinematics only down to a specified zcoordinate, and even omitting the wave kinematics are included.
3.3. Wave Kinematics Including Diffraction Effects
Diffraction effects are described in terms of wave kinematics transfer functions relating wave particle velocities in x, y and zdirections to the undisturbed surface elevation.
The transfer functions for calculation of distributed wave velocities and accelerations \(\mathrm {u_\textrm{j}^d}\) and \(\mathrm {\dot {u}_\textrm{j}^d}\) , respectively, are defined as:
where \(\mathrm {\textrm{j}=1,2,3}\) signifies \(\mathrm {x}\) ,\(\mathrm {y}\) and \(\mathrm {z}\) directions, respectively.
Consistent vessel motion transfer functions and wave kinematics transfer functions for given locations in the vessel coordinate system are available from sink/source computer programs. The sink/source computations are normally carried out in the following major operations:

Solution of global system of equations to determine unknown source intensities for all frequencies and wave directions considered.

Computations of vessel motion transfer functions based on results from 1).

Calculation of wave kinematics transfer functions based on results from 1).
In slender structure applications it will normally be sufficient to consider wave kinematics transfer functions for a few locations along the structure near the vessel (35 locations will tentatively be sufficient). For a limited number of wave kinematics transfer functions, the computation time for step 3) is small when compared to the computational efforts required in step 1).
In order to take advantage of the possibility of splitting the transfer function computations as described above, the following strategy is applied:

Vessel motion transfer functions are given as input to INPMOD (i.e. results from step 1) and 2) in transfer function calculations).

The static configuration is computed by running STAMOD. Wave kinematics transfer functions for specified locations along the static configuration can then be efficiently computed (i.e. step 3 in transfer function calculation) prior to dynamic analysis.

Wave kinematics transfer functions are given directly as input to DYNMOD on a separate input file. The locations of transfer functions are specified by zcoordinates in vessel coordinate system. Transfer functions at intermediate nodes in the diffraction area are found by interpolation along the structure, see figure Figure 5.

Wave kinematics time series are generated based on velocity and acceleration spectra computed from interpolated transfer functions in the diffraction area using the FFT technique. Undisturbed wave kinematics are applied outside the diffraction area. Consistent phase angles are applied in the FFT procedure to obtain consistent vessel motions and undisturbed/disturbed wave kinematics.
By this strategy, only step 3) in the calculation of transfer functions must be carried out for each new static configuration which will give a computationally efficient procedure.
The described approach can be extended to a complete 3D description by generation of wave kinematics transfer functions in a grid covering the dynamic motion range of the structure as indicated in figure Figure 6. It should, however, be noted that significantly more input is required and that a more complex data structure and interpolation procedure is necessary (not implemented in present version).
3.4. Second Order Stokes’ Irregular Waves
3.4.1. Stokes’ theory
In this section, Stokes’ theory is presented. Following the governing set of differential equations, it is shown how a Taylor expansion of the velocity potential around the mean surface level leads to a series expansion in terms of wave steepness. Analytical solutions of the velocity potential and wave elevation are given up to and including second order. This is followed by discussions of issues of importance to the application of Stokes’ theory for irregular waves with a target spectrum.
The wave potential
Stokes’ theory of ocean waves assumes irrotational flow of an incompressible fluid, described by a scalar velocity potential \(\mathrm {\Phi(\vec{x},t)}\) which satisfies
with a velocity field \(\mathrm {\vec{u}}\) given by
Boundary conditions
Assuming an impermeable seafloor of constant mean depth \(\mathrm {h}\), we have \(\mathrm {u_z(z=h)=0}\), which we write in terms of the wave potential as
The free surface of the wave is given by \(\mathrm {z=\eta (x,y,t)}\). For sea states with longcrested waves the \(\mathrm {y}\) coordinate is irrelevant, while for shortcrested waves it is sometimes useful to transform between a cartesian and a polar coordinate system, such that \(\mathrm {z=\eta (r,\theta ,t)}\). Now, \(\mathrm {r}\) plays the role which \(\mathrm {x}\) plays in longcrested waves, and \(\mathrm {\eta }\) results from the superposition of several \(\eta _{\mathrm i}(r,t)\equiv\eta (r,\theta _{\mathrm i},t)\).
Two boundary conditions are imposed at the wave surface. The kinematic boundary condition
which states that a particle at the surface moves together with the surface, and the dynamic boundary condition
which follows from assuming constant pressure along the surface. Here, \(\mathrm {g}\) is the acceleration of gravity.
Together with Equation (26), which applies within the bulk of the ocean, Equation (28), specify the set of differential equations from which \(\mathrm {\Phi}\) is obtained.
Series expansion of the surface elevation
In order to obtain analytic results for \(\mathrm {\Phi}\), a perturbative approach is taken.
The velocity potential at the wave surface is extrapolated onto the mean surface level \(\mathrm {z=0}\) using a Taylor expansion, with terms of higher order in the expansion expected to contribute at higher order in the ordering parameter \(\mathrm {\epsilon}\).
The expansion is, to order \(\mathrm {\mathcal{O}(\epsilon^2)}\),
Likewise, the surface elevation is written as a series expansion,
where it is usually assumed that \(\mathrm {\eta ^{(1)}}\) is a linear wave. This is a construction which allows us to study the nonlinear effects which appear at second order and beyond.
The surface elevation is a function of \(\mathrm {x}\), \(\mathrm {y}\) and \(\mathrm {t}\), or equivalently and for many purposes more conveniently, a function of \(\mathrm {r}\), \(\mathrm {\theta }\) and \(\mathrm {t}\). Staying with the former set of coordinates, the surface elevation in Fourier space is a function of wavenumber \(\mathrm {\vec{k}=(k_x,k_y)}\) and angular frequency \(\mathrm {\omega }\). They are related through the dispersion relation for surface waves at finite waterdepth,
which is complete to second order. At third order, the dispersion relation includes an additional term.
The first order contribution \(\mathrm {\eta ^{(1)}}\) is
provided to WaveKin
by its Fourier components
\(a_{\mathrm i}\),
where the sums run over the number of positive frequencies \(\mathrm {N}\) and the number of directions \(\mathrm {M}\) in the discretized directional spectrum, \(\psi_{\mathrm {ik}}=\omega _{\mathrm i}t\vec{k_{\mathrm {ik}}}\cdot \vec{x}\varepsilon _{\mathrm {ik}}\) is the phase function, with angular frequency $_{i}={i}$. The phase is \(\varepsilon _{\mathrm {ik}}\) and the wavenumber \(\vec{k_{\mathrm {ik}}}\) is obtained from \(\omega _{\mathrm i}\) through use of the dispersion relation and wave direction \(\theta _{\mathrm {ik}}\). The two subscripts denote a discretization of the Fourier components of a directional spectrum. In the following, the second subscript \(\mathrm {k}\) is dropped in favor of a compact notation. It should be implicitly understood that for directional spectra, all sums over Fourier components also imply a sum over directions.
For a simulated wave from a given spectrum, the random phase will typically be drawn from a uniform probability distribution. For a directional spectrum this involves \(\mathrm {MN}\) independent and identically distributed random numbers. See Section 3.4.2 for a further discussion on this.
The sum of Equation (3) can be rewritten to run over positive frequencies only. Because \(\mathrm {\eta (x,y,t)}\) is realvalued, \(a_{\mathrm i}=\bar{a}_{\mathrm i}\). This gives
where the sum over directions is now implicit.
Velocities and accelerations are also realvalued; WaveKin
therefore
works with onesided spectra only. Timeseries may have a nonzero mean
value, therefore the onesided spectra include the zerofrequency
components. The zerofrequency component for surface elevation
\(\mathrm {a_0}\) is also included in the onesided spectrum,
though its value should be zero or near zero. See
Section 3.4.1.6 for how symmetries allow summations
of first and second order terms to run over positive frequencies only
and Section 3.4.1.7 for a further discussion on
the zerofrequency components.
The expansion is in terms of an ordering parameter \(\mathrm {\epsilon}\). It can be seen that higher orders in the expansions of the boundary conditions at the wave surface results in terms with higher powers of spatial derivatives, surface elevation, and the velocity potential. In Fourier space, spatial derivatives result in additional factors of \(\mathrm {k}\), and we have constructed the first order surface elevation to be a sum of harmonics with amplitude \(\mathrm {a}\).
The combination of factors of \(\mathrm {k}\) and \(\mathrm {a}\) provides a physical meaning to the ordering parameter. The wave steepness is proportional to \(\mathrm {ak}\), for a harmonic wave elevation \(\mathrm {\eta _k=a\sin kx}\). In order for the series expansion to converge, the linear wave should have a spectral content which satisfies \(a_{\mathrm i}k_{\mathrm i}<\epsilon\) for all \(\mathrm {i}\).
In addition to the Taylor series expansion of the wave surface boundary conditions and the expansion of \(\mathrm {\eta }\), the velocity potential is written as a series expansion
A solution is obtained for each term of the expansion by ordering the terms in the wave surface boundary conditions, setting \(\mathrm {z=\eta }\) in Equation (3). Terms at equal order in \(\mathrm {\epsilon}\) are solved for independently.
First order solution
Solving for the first order velocity potential gives
Second order solution
Following Marthinsen and Winterstein (1992) for longcrested waves and Marthinsen and Winterstein (1992) for shortcrested waves, the velocity potential and second order contribution to surface elevation is given to second order by a set of rather complex expressions, organized in terms of sumfrequency and differencefrequency contributions. The expressions are analytic, and therefore support differentiation, allowing the analytic calculation of velocities and accelerations.
Solving for the second order surface elevation gives
Solving for the second order velocity potential gives
where the latter term is independent of spatial coordinates and therefore irrelevant for kinematics.
The second order kernels \(E_{\mathrm {i,j}}\) and \(P_{\mathrm {i,j}}\) are given by
and
The solutions for both \(\mathrm {\eta ^{(2)}}\) and \(\mathrm {\Phi^{(2)}}\) depend on products, sums and differences of Fourier components. Because the twosided spectrum can be expressed as a onesided spectrum, the second order contributions are naturally separated into sumfrequency and differencefrequency contributions.
Because the sum and differencefrequency components only depend on the sum and difference of frequencies (wavenumbers) and the Fourier components satisfy \(a_{\mathrm i}=\bar{a}_{\mathrm i}\), the calculations of second order contributions to surface elevation, velocities and accelerations can be considerably simplified.
The sumfrequency and differencefrequency transfer functions
The kernels from Equation (41) and Equation (42) can be written in terms of sums and differences of frequencies and corresponding wavenumbers,
and
where \(\mathrm {i,j}>0\).
The symmetry properties of the second order kernels are
where \(\mathrm {i,j}>0\).
For the second order wave elevation of Equation (3) we now have \(\mathrm {\eta ^{(2)}=\eta ^{}\eta ^{}}\) with
The second order velocity potential of Equation (40) can now be written as \(\mathrm {\Phi^{(2)}=\Phi^{}\Phi^{}}\) with
where spaceindependent terms have been omitted.
The diagonal of the differencefrequency kernels
The differencefrequency kernels \(E_{\mathrm {ij}}^\) and \(P_{\mathrm {ij}}^\) are singular for \(\mathrm {i}=\mathrm {j}\). If one chooses to take the analytic limit, one gets a nonzero contribution to the zerofrequency contribution of the differencefrequency transfer functions. For \(E_{\mathrm {ii}}^\) this gives
where
is the group velocity.
With nonzero contributions from the differencefrequency kernels when \(\mathrm {i}=\mathrm {j}\), the resulting transfer functions give finite contributions to the zerofrequency component of the second order elevation and velocity potential. This corresponds to finite meanvalues for \(\mathrm {\eta }\), \(\mathrm {u_x}\) and \(\mathrm {u_y}\). This can optionally be avoided by setting all zerofrequency components to zero, (equivalently, setting the diagonal entries of the differencefrequency kernels to zero) however this procedure lacks physical or theoretical justification.
In particular, the socalled mean sea state setdown leads to an
inconsistency with the boundary condition applied at the sea floor,
i.e. at \(\mathrm {z=h}\), as the sea floor depth
\(\mathrm {h}\) should properly be measured from the mean sea
level, \(\mathrm {<\eta >}\), not from the first order
contribution to the same, \(\mathrm {<\eta ^{(1)}>}\). The
inconsistency might be small for most practical uses, but deserves
mentioning. This is discussed occasionally in the literature, e.g.,
Marthinsen and Winterstein (1992) and
The default behaviour of WaveKin
in RIFLEX
is to set the diagonal
terms equal to zero.
High frequency contributions from the sumfrequency transfer functions
Depending on how the target spectrum is filtered to obtain the first order wave elevation defined by \(\{a_{\mathrm {i}}\}\), the sumfrequency transfer functions may result in significant high frequency contributions to the second order wave elevation and kinematics. The acceleration terms are particularly sensitive to the high frequency part of the spectrum. Care must be taken when defining the first order wave elevation to avoid unphysical high frequency terms arising at second order. Otherwise, unphysically large accelerations could result.
The issue of large accelerations, i.e., rapidly varying velocities, is most acute near the mean sea surface and above (\(\mathrm {z\ge0}\)).
3.4.2. Spectrum
The input to WaveKin
is the first order wave elevation. From this, the
second order contribution to wave elevation is calculated and added to
the first order part. Also, kinematics are calculated complete to second
order.
The input wave elevation may be provided as a directional spectrum, in
which case the 2D array (frequency and direction) of Fourier components
must be accompanied by an array of directions. The generation of the
first order wave elevation is done outside of WaveKin
.
The Fourier components may be taken deterministically from the power spectrum. For a time series of finite length there will be some variability of the power spectrum with respect to the target spectrum, and thus individual Fourier components should be \(\mathrm {\chi^2}\) distributed. However, it should be noted that when \(\mathrm {\eta }\) is given as a sum of directional components this in itself results in a natural source of variability for the Fourier components at any given frequency.
The linear wave and the filtering procedure
By target spectrum is meant the empirically established power spectrum \(\mathrm {S(\omega )}\) for longcrested waves or the directional spectrum \(\mathrm {S(\omega ,\theta )}\) for directional sea states.
In order to generate an irregular sea state, phases are typically drawn from a uniform distribution, with the amplitude of Fourier components fixed by a filtered power spectrum. The resulting wave elevation is then referred to as a linear wave and identified as the first order contribution in the perturbative expansion of \(\mathrm {\eta }\).
Various heuristics are presented in the literature for how the filtering
should be done, see e.g. Marthinsen and Winterstein (1992) and
Marthinsen and Winterstein (1992), but no solid theoretical foundation underlies
these methods. For many purposes, results may be insensitive to the
details of the filtering, however results from WaveKin
indicate that
filtering has a significant effect on wave kinematics near the sea
surface at water depths of \(2040{\mathrm {~m}}\) and
significant wave heights in the range of
\(56{\mathrm {~m}}\).
Equating the linear wave with the first order wave elevation is a matter of convenience, but could underestimate the nonlinear nature of Stokes’ waves. The assumption of uncorrelated phases is justified by very weak wavewave interactions, which in turn applies to low wave steepness. While the first order wave elevation will always have a more narrow banded spectrum than the full wave, it will nevertheless contain steep waves.
Separability of the directional spectrum
A directional spreading function \(\mathrm {D(\theta ,\omega )}\) is defined by
together with the normalization condition
The power spectrum \(\mathrm {S(\omega )}\) for directional sea states is obtained by integrating over the directional coordinate
The simplest approach to describing a directional sea state is to assume a separable directional spectrum,
where the power spectrum \(\mathrm {S(\omega )}\) (e.g., a JONSWAP spectrum) and the directional spreading function \(\mathrm {D(\theta )}\) are given separately and together define the sea state.
A common choice for \(\mathrm {D}\) is a frequencyindependent cosine2s distribution
for \(\mathrm {\theta \theta _0\leq\frac{\pi }{2}}\) and zero otherwise. \(\mathrm {C_s}\) is a normalization factor. The width of the distribution is parametrized by \(\mathrm {s}\) and the mean direction is given by \(\mathrm {\theta _0}\).
3.4.3. Kinematics
Velocities and accelerations are obtained from applying partial derivatives to the wave potential, which is known analytically up to and including second order.
The velocity and acceleration to second order are given by
WaveKin
by default performs differentiation by applying partial
derivatives in Fourier space. This is efficient and accurate, and has
been validated against differentiating along temporal and spatial
coordinates. An option exists to perform differentiation numerically
along the temporal coordinate, which is turned off by default.
Convection terms
Frequently, convective terms are neglected in the calculation of wave
kinematics. In the current implementation of WaveKin
the inclusion of
these terms is optional. If they are omitted, only the contribution from
partial timederivatives are included
Convection terms can optionally be included in the accelerations
returned by WaveKin
. These are calculated by taking spatial
derivatives of the velocity field (in the frequency domain) and
multiplying with the velocity field (in the time domain).
Convection terms contribute to the accelerations of momentumcarrying fluid particles, and thus contribute to the inertia term of Morison’s load model. They emerge when we take the total derivative of velocity with respect to time,
where \(\mathrm {\vec{v}}\) is the velocity vector and \(\mathrm {\vec{u}=\vec{u}(\vec{x},t)}\) is the velocity field.
The inclusion of convection terms is turned off by default in WaveKin
.
This follows an assumption that they are small and can be neglected for
practical purposes. This assumption should be questioned for severe sea
states. The last term in Equation (61) constitutes the convection
terms. It should be seen that its magnitude depends on the gradient of
the velocity field, meaning that rapidly varying velocities, as found
near the sea surface in severe sea states, may lead to significant
contributions from convection terms.
Kinematics near to and above the mean surface level
The velocity potential \(\mathrm {\Phi(x,y,z,t)}\) is extrapolated above \(\mathrm {z=0}\) according to the Taylor expansion of Equation (3). Kinematics are obtained from the extrapolated potential by spatial and temporal differentiation in the same manner as for the velocity potential at and below \(\mathrm {z=0}\). However, it should be noted that the second order contribution at \(\mathrm {z>0}\) only involves the first order velocity potential \(\mathrm {\Phi^{(1)}}\). There are no quadratic transfer functions involved above mean surface level until third order contributions.
High frequencies resulting from providing an improperly filtered spectrum as the input first order wave may cause unphysical accelerations near to and above the mean surface level. To compensate for this, an option is provided to filter out large frequencies prior to extrapolating the kinematics above the mean surface level.
Together with, or instead of this option, there is an option to set a threshold surface level \(z_{\mathrm {lim}}<0\), above which kinematics are extrapolated from the potential at \(z=z_{\mathrm {lim}}\). As the large frequencies are rapidly dampened even at very small values of \(\mathrm {z<0}\), setting \(z_{\mathrm {lim}}\) to only a few centimeters below \(\mathrm {z=0}\) may be sufficient to avoid the excessive spikes in acceleration.
The default behaviour of WaveKin
in RIFLEX
is to use the entire
provided spectrum when extrapolating kinematics above the mean surface
level, and to extrapolate from \(\mathrm {z=0}\).
Kinematics above the sea surface
In order to extrapolate kinematics between two nodes, one of which is
dry, it can be impractical to set kinematics to zero above the sea
surface \(\mathrm {z>\eta }\) (sometimes referred to as a
dry node
). Instead, the option exists to set kinematics at
\(\mathrm {z>\eta }\) equal to the the kinematics at the
surface: \(\mathrm {z=\eta }\). This is the default behaviour
of WaveKin
in RIFLEX
.
Kinematics on the surface and below the mean surface level
For kinematics on the surface and at or above the mean surface level, \(\mathrm {z=\eta \ge0}\), the kinematics are obtained at second order by extrapolation of the first order kinematics at \(\mathrm {z=0}\). For kinematics on the surface and below the mean surface level, \(\mathrm {z=\eta <0}\), there is a choice between using the second order potential, which exists at all \(\mathrm {z\le0}\) or using the first order potential at \(\mathrm {z=\eta <0}\) together with the first order potential extrapolated from \(\mathrm {z=0}\). The latter method provides a consistent calculation of kinematics on the surface regardless of whether it is above or below \(\mathrm {z=0}\). However, it causes a discontinuity in the kinematics on and directly beneath the surface.
The magnitude of the discrepancy between the two methods should be of third order, if the series expansion behaves consistently. Thus, a comparison of the two methods is a test of the series expansion.
4. Current Description
The current velocity is normally assumed to be constant with time. A quasistatic analysis can, however, include the effect from a quasistatic current variation, and it is also possible to describe timedependent current variation in dynamic analysis.
The current velocity at a given position is described by speed and direction. This can be done by input of discrete values and interpolation to actual node positions, or by definition of standard profiles.
The velocity components in global directions are denoted
The vertical component \(\mathrm {U^{c}_z}\) is always assumed to be zero.
5. Hydrodynamic Load Models for Submerged Elements
This section will give an introduction to hydrodynamic load models for submerged or surfacepiercing normal slender structures (i.e. pipes, cables, risers, etc.).
5.1. Generalized Morison’s Equation
The velocity potential \(\mathrm {\phi ^I}\) of the undisturbed wave motion will be written as a superposition of linear waves of different circular frequencies \(\mathrm {\omega _\textrm{i}}\) and wave propagation angles \(\mathrm {\omega _\textrm{j}}\) ,
Here \(\mathrm {g}\) is the acceleration due to gravity, \(\mathrm {t}\) is the time variable and \(\mathrm {(X,Y,Z)}\) is a Cartesian coordinate system with \(\mathrm {Z=0}\) in the mean free surface level.
Positive Zaxis is upwards. The mean water depth is \(\mathrm {D}\) . The wave numbers \(\mathrm {k_\textrm{k}}\) are related to the wave frequency \(\mathrm {\omega _\textrm{k}}\) by the dispersion relationship \(\mathrm {\omega _\textrm{k}^2/g=k_\textrm{k}\textrm{tanh}(k_\textrm{k}h)}\). \(\mathrm {\phi _\textrm{jk}}\) are random phase angles and \(\mathrm {A_\textrm{jk}}\) are determined by a wave spectrum for short crested sea.
When calculating the free surface elevation \(\mathrm {\eta }\) we will use that
This value will be used to position the calculated velocity and acceleration profiles, which enables us in an easy way to calculate the loads on structural surfacepiercing parts.
A schematic picture of a possible flexible riser system is shown in Figure 7. Crosssections of typical flexible pipe configurations are shown in Figure 8.
The hydrodynamic forces are calculated according to a generalized Morison’s equation formulation, with the optional addition of a linear drag force term. A local coordinate system \(\mathrm {(x,y,z)}\) is defined. See figure Figure 7 where \(\mathrm {e(y,z)}\) is in the crosssectional plane. Note that small letters indicate the local coordinate system while large letters \(\mathrm {(X,Y,Z)}\) indicate the global coordinate system.
5.1.1. Bisymmetric cross section
The motion and the water velocity are decomposed into 3 components in the local element coordinate system, with \(\mathrm {x}\) in the axial direction, \(\mathrm {y}\) and \(\mathrm {z}\) normal to the axis, cfr. Figure 8.
It is assumed that both the \(\mathrm {y}\)  and \(\mathrm {z}\) axis are symmetryaxes for the crosssectional surface geometry. The longitudinal, \(\mathrm {F_x}\) , force components and transverse force components \(\mathrm {F_y}\) and \(\mathrm {F_z}\) , on an element of length \(\mathrm {dx}\) subjected to waves and current can then be written:
here dots denote time derivatives and

\(\mathrm {\rho \qquad\qquad\qquad\,}\) : mass density of water

\(\mathrm {A\qquad\qquad\qquad}\) : external crosssectional area

\(\mathrm {m_a^x,m_a^y,m_a^z\quad \;\;\,\,}\) : 2D added mass in local x, yand zdirections

\(\mathrm {(r_x,r_y,r_z)\qquad\:\:}\) : translational structural displacement components in the local coordinate system

\(\mathrm {C_D^x,C_D^y,C_D^z\quad \;\;\:}\) : dimensional quadratic drag coefficients in local x, y and zdirections

\(\mathrm {C_{DL}^x,C_{DL}^y,C_{DL}^z}\) : dimensional linear drag coefficients in local x, y and zdirections

\(\mathrm {(U_x^c,U_y^c,U_z^c)\quad \;\;}\) : current velocity components in the local coordinate system

\(\mathrm {(u_x^I,u_y^I,u_z^I)\qquad\;}\) : velocity component of undisturbed wave field in the local coordinate system
See also the User Manual for the relation between these force coefficients and commonly used nondimensional coefficients.
It should be noted that oscillatory forces due to vortex shedding are neglected. The drag coefficients and added mass coefficients will in general depend on structural form, Reynolds number, roughness ratio, KeuleganCarpenter number, current velocity relative to wave induced velocity, etc. The nature of the orbital motion of the undisturbed local wave field will have an influence. For circular cylinders there is some information available in the literature.
5.1.2. Rotation symmetric cross section
The linear part of the drag force and the inertia forces, which are also linear, are calculated as for the bisymmetric cross section. The longitudinal quadratic drag force is also calculated as for the bisymmetric cross section.
The transverse quadratic drag force is calculated differently, using the total relative velocity normal to the line axis:
where
The hydrodynamic rollmoment (i.e. moment about the xaxis) on an element of length \(\mathrm {dx}\) is given by:
Here \(\mathrm {\ddot \eta _4}\) is the roll angle. Further \(\mathrm {A_{44}^{(2D)}}\) is the twodimensional added moment in roll with respect to the local xaxis. For steady current we can write the Munk moment as
\(\mathrm {M_x}\) is not implemented in the present version.
At every time step, forces and moments are calculated along the whole instantaneously submerged length of the riser. It is earlier discussed how the free surface position \(\mathrm {\eta }\) is calculated, which is important information to calculate the submerged length of the riser.
In the case of additional buoys, the hydrodynamic forces is handled in a similar way as described for the pipe cross sections.
5.2. MacCamyFuchs with Quadratic Drag
For structures with relatively large diameter compared to the wave length, nearfield diffraction effects become more important. An analytical solution for the potential around a vertical, surfacepiercing circular cylinder was developed by MacCamy and Fuchs.
For elements with MacCamyFuchs type loading, the formulation for the inertia and diffraction forces follows (MacCamy and Fuchs, 1954 and Dean and Dalrymple, 1991), with modifications for the irregular wave history. The force per unit length (\(\mathrm {dF}\)) at \(\mathrm {x=0}\) for a given elevation \(\mathrm {z}\) is given by:
where
and
and the wave components are defined as for the wave potential.
The resulting \(\mathrm {C_m}\) and phase angle for a section with diameter 10 m are shown in Figure 9.
The model is only applicable for circular crosssections and approximately vertical sections. The forces are computed based on the zcoordinate in the center of the element and are applied in a lumped manner.
In addition to the inertia force, a quadratic drag force based on linear wave kinematics may also be included when using MacCamyFuchs. The drag force is according to Morison’s Equation, with the same implementation as in section Section 5.1.
5.3. Netloads
The method for calculation of current forces on fish farms is described in detail in Løland, G. (1991). The main equations are shown in this section. For details see Løland (1991) and also Aaarsnes et al (1998), Ormberg (1991).
The method is based on that the mean drag and lift force on a net panel can be written as
where

\(\mathrm {C_D}\) = drag coefficient as a function of the angle between net normal and flow direction

\(\mathrm {C_L}\) = lift coefficient as a function of the angle between net normal and flow direction

\(\mathrm {A}\) = area of net panel

\(\mathrm {U}\) = relative velocity

\(\mathrm {\rho}\) = density of water

\(\mathrm {\alpha}\) = angle between the flow direction and the net normal vector in the direction of the flow
The relative velocity is the difference between the wave and current velocity and the structural velocity. The current velocity is corrected using the reduced current velocity factor (the ratio between reduced current speed and ambient current speed due to upstream net shadowing effect). 
The drag force is defined as the force in the direction of the flow. The lift force is defined as the force normal to the flow direction. In the global axis system, the force components can be written as
\(\mathrm {\boldsymbol{e}_T}\) is the unit tangential vector defining the direction from the upper end to the lower end point of the net panel.
The drag and lift coefficient are function of the solidity ratio. The solidity ratio is degined as the ratio between the area covered by the threads in the screen and the total area of the screen.
\(\mathrm {\lambda}\) is the mesh size and \(\mathrm {D}\) the twine diameter. The relationship between the drag coefficient and the solidity ratio is given as
The factor 0.04 is introduced to take into account the drag on a net panel parallel to the flow. The factor is assumed to be independent of the solidity ratio. This is an approximation. The model test indicated that the factor is almost independent of the solidity ratio, \(\mathrm {S_n}\) , for the solidity ratio range used in the model test. This range was 0.13 to 0.317 as well as a angle \(\mathrm {\alpha}\) from 090 degrees.
5.3.1. References
Milne, P.H. (1970): Fish farming: A Guide to the Design and Construction of NetEnclosures, Department of agriculture and Fisheries for Scotland, Marine ResearchRep. No. 1, 1970
Aarsnes, J.V., Lien E. and Oltedal, G. (1988): SIFICA  Theory manual, MARINTEK report: 513004.05.02, Trondheim
Løland, G. (1991): Current Forces on and Flow Through Fish Farms, Doctoral thesis, NTH, Trondheim, Norway
Ormberg, H. (1991): Nonlinear Response Analysis of Floating Fish Farm Systems. Dr.ing dissertation,The Norwegian Institute of Technology. ISBN 8271193309
Løland (1993): Current forces on, and water flow through and around, floating fish farms ,Aquaculture International 1, Pages 7289
Kristiansen, T. et al (2012): Modelling of current loads on aqua culture net cages Journal of Fluids and Structures Volume 34, October 2012, Pages 218235
6. Forced Vessel Motion
In many cases, the slender structure will be connected to a vessel at one end. Motions of this vessel must therefore be known during a dynamic riser analysis. Generation of motion time series will be consistent with generated time series for waveinduced water particle velocities and accelerations. The rigid body motion responses consist of 6 degrees of freedom: surge, sway, heave, roll, pitch and yaw referred to the global \(\mathrm {(X,Y,Z)}\) coordinate system. The motions are in this chapter denoted by \(\mathrm {\boldsymbol{x}}\) .
The motion model consists of a set of highfrequency (wave frequency) motions in all 6 degrees of freedom and a set of lowfrequency motions in the 3 horizontal degrees of freedom: surge, sway and yaw. These two sets of motions are referred to as HFmotions and LFmotions, respectively. For most dynamic line problems it is sufficient to include only the HFmotions. The effects of typical LFmotions, with periods of 60180 s, can often be covered by suitable selection of static (mean) position.
6.1. HF Vessel Motion Model
The wave frequency, or HF motions are treated as linear responses to
the waves. Thus, the HF motions are described by a set of complex
transfer functions \(\mathrm {j=1,2.....6}\)
where:
6.2. LF Vessel Motion Model
The lowfrequency, or LFmotions, if included, are represented by directly input response spectra for surge, sway and yaw. The LF motions modelled this way are uncorrelated with the wave, and with each other.
Both the transfer functions of HF motions and spectra for LF motions refer to a specified vessel origin point. They are transformed to the structural attachment points on the vessel, (i.e. upper end terminal point for standard systems). This information together with the wave spectrum description, is sufficient to generate spectra and the sample time series of forced vessel motions.
6.3. Frequencydomain Vessel Motion Analysis
6.3.1. HF Responses
In the linearized, frequencydomain analysis the wave frequency motions are completely described by their autospectra, \(\mathrm {S_x}\) .
Response spectra are generated for all 6 rigidbody motions.
where \(\mathrm {H_{xj}'(\beta ,\omega )}\) is the transfer function referring to the top terminal point.
6.3.2. LFResponses
\(\mathrm {S_{LFxj}(\omega )}\) is given as direct input and subsequently used for time series generation.
6.3.3. HF time series generation
Input to the HF time series generation is the complex motion transfer function \(\mathrm {H_j(\beta _k,\omega _l)}\) , and the harmonic wave components, \(\mathrm {Z_{kl}}\) .
The harmonic components of the motion responses are simply:
The transfer functions \(\mathrm {H}\) are stored on a tabular form for given arrays of directions and frequencies different from \(\mathrm {\beta _\textrm{k}}\) and \(\mathrm {\omega _l}\) , intermediate values are obtained by interpolation (3rd order spline interpolation). Symmetry properties of the platform are utilized.
Time series are obtained by adding the harmonic components by means of an FFT algorithm.
For regular wave analysis, the HFtime series may optionally be given as direct input of amplitudes and phase angles of the top terminal point of the structure.
6.3.4. LF time series generation
Input to the LF time series generation are the complex motion response spectra\(\mathrm {S_{LF,xj}(\omega ),j=1,2,6}\)
The harmonic components of the motion responses are:
The random phase angles, \(\mathrm {\phi _{xjk}}\) , are sampled from a uniform distribution over \(\mathrm {[\pi ,\pi ]}\) .
Addition of the harmonic components to obtain time series is performed by a Fourier Transform algorithm.
The relations between time increment, \(\mathrm {\Delta t_{LF}}\) , frequency increment \(\mathrm {\Delta \omega _{LF}}\) and number of time steps for the LFsimulations is equivalent to those described in Equation (20) except that, \(\mathrm {t_{\textrm{lim}}}\) , the minimum interesting low frequency response period is greater, typically 2030 sec.
6.3.5. Total Responses
Time series of total responses for the surge, sway and yaw motions are obtained by adding the HF and LF contributions directly.
If LFresponses have been generated with larger time increment than the HFresponses, a 3rd order spline interpolation is applied to obtain identical time increments.
7. Effects from Internal Fluid Flow
The riser system model includes specification of internal fluid flow by:
Parameter  Symbol  Dimension 

Cross section area of internal volume 
\(\mathrm {A_i}\) 
\(\mathrm {L^2}\) 
Density 
\(\mathrm {\rho _i}\) 
\(\mathrm {M/L^2}\) 
Velocity 
\(\mathrm {v_i}\) 
\(\mathrm {L/T}\) 
Pressure gradient due to flow resistance 
\(\mathrm {\rho d_i}\) 
\(\mathrm {F/L^3}\) 
Pressure at end 1, including the effect of flow density 
\(\mathrm {\rho _{1i}}\) 
\(\mathrm {F/L^2}\) 
The following simplifications are made:

The fluid is incompressible.

The velocity and pressure distribution are assumed to be unaffected by the riser motion.

The pressure is included only in order to account for possible future pressuredependent stiffness of nonbounded flexible pipes and for axial wall force evaluation.
The internal flow as modelled will affect the riser configuration and behaviour in the following ways:

The effect of weight and mass
Vertical force per unit length due to internal fluid:
This weight will affect the effective tension and thus the configuration of the riser. The weight will also affect the internal pressure distribution.

The effect of internal pressure
The bending stiffness of the nonbounded flexible pipes is assumed to be dependent upon the pressure differential,
where \(\mathrm {\rho _w,\rho _i}\) = density of water and internal fluid, respectively
\(\mathrm {p_{di}}\) = pressure drop per unit length due to flow resistance
\(\mathrm {s}\) = coordinate along the riser
\(\mathrm {p_{1i}}\) = internal pressure at lower end (end 1)
This information is not utilized in the present version.

The static effect of the impulse and centrifugal force
The centrifugal force effect may be included in the riser loads. This force depends upon local curvature, velocity and mass of fluid flow.
The centrifugal force per unit length can be expressed by:
This force will act radially on any curved part of the riser. However, the force will also contribute to increase the effective riser tension:
Comparing the two cases a) and b) in Figure 11 for a case with zero bending stiffness, the radii of curvature can be expressed by:
Thus, the centrifugal force will not contribute to change the static configuration. The only effect is the increase of wall tension due to the impulse force, according to Equation (93) written:
where \(\mathrm {\boldsymbol{\dot v}}\) is the acceleration vector of the fluid. Assuming stationary, incompressible flow, the acceleration can be written:
where
\(\mathrm {\frac{\partial \boldsymbol{e}}{\partial s}}\) = curvature vector
\(\mathrm {\boldsymbol{r}}\) = acceleration of pipe
\(\mathrm {\boldsymbol{\omega }}\) = rotation of pipe
\(\mathrm {\boldsymbol{e}}\) = pipe tangent unit vector
\(\mathrm {s}\) = curve coordinate along the pipe
Thus, the following force terms are derived:
The force terms accounts for:

Force due to pipe acceleration

Force due to centripetal acceleration

Force due to pipe rotation (Coriolis force)
The present program version includes only the first term. It is included directly in the mass matrix of the elements.
A special slug modelling option is included to simulate the effect of a relatively short fluid slug moving through a riser. This is done by updating the mass matrix according to a specified slug speed in the pipe.
The second term is identical with the static centrifugal force discussed previously, but the curvature will be timedependent in the dynamic case.
The third term may have more complicated effects.
Both the second and the third terms may be included directly as forces in the time domain analysis in a possible future extension of the program.
8. Hydrodynamic Load on Partly Submerged, Floating Elements
This section will give an introduction to the hydrodynamic load models for floating partly submerged elements, which are relevant for analysis of floating flexible fish farm systems, floating offshore loading hoses, etc. The presentation is to a large extent based on Ormberg (1991).
A floating body of small volume may, during a wave cycle, accomplish large motions from its mean position. Thus, the hydrodynamic loads should be computed at the instantaneous position rather than at the mean position, which means that the assumptions within linear potential theory are violated. The applied hydrodynamic load model, described by Aarsnes et al. (1988), is an engineering approach to a solution of the problem. It involves hydrodynamic coefficients deduced by linear theory. In addition, an adequate description of the wave kinematics in the wave zone is required. Application of the model including comparison with model test results is described by Ormberg (1991).
The hydrodynamic forces are calculated based on twodimensional strip theory, which requires the structural members to be long and slender. Hydrodynamic interaction between structural parts is neglected. The waveinduced excitation forces (FroudeKrylov and diffraction forces) are computed by a long wavelength approximation which involves added mass and potential damping of the actual cross section together with the wave kinematics. The viscous loads are computed using the drag term in Morison’s equation.
In the following, a rigid structural member (section) with constant cross sectional dimensions, is considered. The section with local (corotating) coordinate system \(\mathrm {(x,y,z)}\) , is illustrated in Figure 12. The figure also defines the spacefixed (global) coordinate system with X and Yaxes in the still water plane and with the Zaxis pointing upwards.
As indicated in Figure 12 , the section is divided into subelements for hydrodynamic load calculation. The forces acting on a subelement are calculated based on the instantaneous orientation and position relative to the free surface. The water particle kinematics, which are assumed to be constant within the subelement, are calculated at the buoyancy center of the subelement. Correspondingly, the load per unit length is constant over the length of the subelement and attacks through the buoyancy center. The total force normal to the section is found by simply adding the force contribution on each subelement.
For the computer program implementation, a table of the twodimensional added mass and damping coefficients for different levels of immersion at the actual wave frequency is precomputed by use of the Frank close fit method, Frank (1967).
In the following, the hydrodynamic forces are understood as the total forces acting on the member from the surrounding fluid. This implies that buoyancy and current forces are also included in the hydrodynamic forces. Thus, the total hydrodynamic forces are given as
where:
\(\mathrm {\boldsymbol{F^{Pot}}}\) : the potential flow contribution to hydrodynamic forces
\(\mathrm {\boldsymbol{F^{FK}}}\) : FroudeKrylov forces including buoyancy forces
\(\mathrm {\boldsymbol{F^S}}\) : diffraction forces
\(\mathrm {\boldsymbol{F^R}}\) : added mass and damping forces
\(\mathrm {\boldsymbol{F^D}}\) : the viscous forces or drag forces
8.1. Potential Flow Contribution
8.1.1. Loads in longitudinal direction
The long and slender shape assumption implies that the forces on the end surfaces are dominated by the FroudeKrylov contribution. Hence, other contributions are neglected. The importance of including the pressure force on end plates is thoroughly discussed by Hooft (1972).
The end pressure forces are calculated based on dynamic FroudeKrylov pressure and the hydrostatic pressure
where \(\mathrm {g}\) is the acceleration due to gravity and \(\mathrm {\rho }\) is the density of water.
Thus, the total potential contribution to load in local xdirection of the structural member is found by the vector sum of the end pressure forces, which formally is written
where \(\mathrm {A_w}\) is the instantaneous immersed area of the end surface. The indexes (a) and (b) indicate the end surface number, see Figure 13 For practical purposes it is assumed that the following simplified formula is sufficient
where \(\mathrm {p_\textrm{tot}}\) is calculated at the wetted area center.
8.1.2. FroudeKrylov (FK) forces:
The FroudeKrylov forces including the buoyancy forces acting normal to a subelement may according to Kaplan and Hue (1959) and Dixon et al. (1979), be computed by
where \(\mathrm {\boldsymbol{f_n^{\textrm{Buoy}}}}\) is the buoyancy force calculated based on the instantaneous immersed cross section area. \(\mathrm {\Delta x}\) is the length of the subelement. \(\mathrm {\boldsymbol{\dot u_n}}\) is the water particle acceleration vector component normal to the subelement, i.e.
The buoyancy force referred to the global system is simply calculated by
where \(\mathrm {A_S}\) is the immersed cross section area of the subelement.
Transforming the buoyancy load to local system, neglecting the axial contribution to avoid double representation of this term, the FKforce referred to local system is written
where \(\mathrm {\boldsymbol{T_{ij}}}\) is the transformation matrix between the global and local coordinate systems, i.e.
8.1.3. Diffraction forces:
According to the long wave assumption, the diffraction forces may be expressed in terms of twodimensional added mass and damping
where \(A^{(2D)}=A^{(2D)}(Z_\mathrm {rel})\) and \(B^{(2D)}=B^{(2D)}(Z_{\mathrm {rel}})\) are positiondependent added mass and damping per unit length, respectively. \(Z_\mathrm {rel}\) is the vertical position relative the free surface.
8.1.4. Added mass and damping forces:
The added mass and damping forces including moment around local xaxis are calculated according to
where \(\mathrm {\boldsymbol{\dot v}}\) and \(\mathrm {\boldsymbol{\ddot v}}\) are the structural velocity and acceleration, respectively.
8.2. Simplified Approach
A further simplification of the hydrodynamic problem is included by introducing constant added mass coefficients, i.e. neglecting the depth dependency. This means that the added mass is taken as proportional to the submerged cross sectional area. Accordingly (Eq. 7.99) and (Eq. (7.100) simplify to
and
where \(\mathrm {A_S}\) is submerged cross section area. \(\mathrm {C_{my}}\) and \(\mathrm {C_{mz}}\) are constant added mass coefficients referred to completely submerged cross section.
8.3. Viscous Forces
The viscous forces are computed using the drag force term of Morison’s equation. The cross flow principle first introduced by Hoerner (1965), is used to determine the forces in transverse direction of the subelements. The drag coefficients are assumed to be constant and given for a submerged cross section.
The viscous forces or drag forces are calculated based on relative velocity (\(\mathrm {\boldsymbol{u_r}}\)) in the local system expressed by
where \(\mathrm {u_c}\) is the current velocity, \(\mathrm {u_w}\) is the wave particle velocity and \(\mathrm {\dot v}\) is the structural velocity.
8.3.1. Loads in longitudinal direction
A friction force contribution is included in the axial direction
where \(\mathrm {C_{Dx}}\) is the skin friction coefficient and \(\mathrm {L_W}\) is the instantaneous wetted part of the cross section circumference.
8.3.2. Transverse loads
From Equation (111), the transverse relative velocity vector is written
The transverse drag force is calculated according to
where
\(\begin{array}{l}\displaystyle h_\mathrm {rel}=\frac{A_S}{A}h\\\\\displaystyle b_\mathrm {rel}\,=\frac{A_S}{A}b\end{array}\)
\(\mathrm {C_{Dy}}\) and \(\mathrm {C_{Dz}}\) are drag coefficients in local y and z directions respectively, \(\mathrm {A_S}\) is the instantaneous submerged cross section area and \(\mathrm {A}\) is the cross section area. \(\mathrm {b}\) and \(\mathrm {h}\) are characteristic width (ydirection) and height (zdirection) of the section respectively.
8.4. Moment around Local xaxis
Assuming the transverse forces attack at the center of buoyancy, see Figure 14, the moment around the local xaxis due to wave excitation forces and drag forces is approximated by
where \(\mathrm {\boldsymbol{r_B}=y_B\boldsymbol{i_2}+z_B\boldsymbol{i_3}}\) is the position vector of the buoyancy center and the transverse hydrodynamic excitation forces are given by
8.5. Discussion
Computation of added mass and potential damping according to the Frank close fit method requires that the structural member is horizontal, i.e. the local xyplane is parallel with the still water plane, and that the cross section is symmetric about its local xz plane. During time domain simulation, the added mass and damping are taken as directiondependent coefficients connected to local (structural) axis directions, while strictly they are valid for global directions. Hence, introduced as a practical limit, the instantaneous pitch and roll angles have to be within \(\mathrm {\pm30^{\circ}}\) inclination to the surface.
The loads due to slamming or water impact are not included. These may be included as the time rate of change of momentum associated with the immersed portion of the structure. Calculating the vertical impact forces on horizontal members for the water entry case is well documented by e.g. Greenhow (1987), Kaplan and Silbert (1976), Faltinsen (1977). However, in this more general case regarding structure orientation, water exit as well as entry, further investigation is needed. Hence, this term is neglected.
The long wave length assumption is used in the calculation of potential flow acting on the subelement, which normally requires the wave length / diameter ratio to be greater than 5. It is further assumed that the local angular displacements are not too large relative to the free surface. This assumption is due to that the computation of hydrodynamic coefficients are precalculated for a horizontal structural member.
Experiments indicate that the drag coefficients for a cross section in the free surface zone are strongly dependent on the submergence, see e.g. Øritsland (1986). These effects are not included. The dependency on Reynolds number, KeuleganCarpenter number and roughness should be taken into account when specifying the constant drag coefficients as input.
8.6. Implementation in Finite Element Formulation
Different formulations might be used for the numerical implementation of hydrodynamic loads. These consider the order of the polynomial used for the displacement functions and the representation of the distribution of the hydrodynamic element loads. In the current version of the program, a consistent formulation is adopted for the numerical implementation of hydrodynamic loads, i.e. using the same displacement functions as for the element deformation description. The contributions from added mass and damping forces are included in the inertia and damping load vectors of the system, while the other hydrodynamical loads contribute to the external load vector.
Analogous to hydrodynamic strip theory, the load distribution is taken care of by dividing the element into a specified number of subelements for load calculation, see figure Figure 15. Within a subelement, the hydrodynamic loads are assumed to be constant. The corresponding hydrodynamic mass and damping matrices and nodal load vector are computed by numerical integration simply by adding up the contribution from each subelement. This is certainly not the most accurate method for numerical integration, but is found to be adequate considering the uncertainties in the hydrodynamic load computation.
Based on this formulation, the element length may be chosen with regard to structural deformations of the actual problem, while the lengths of the subelements are chosen with regard to hydrodynamic load distribution independently of the element length.
The hydrodynamic loads may in general be described as nonconservative loads i.e. dependent on position, direction and/or loaded area. By linearizing the incremental dynamic equation, tangential mass, damping and stiffness matrices are introduced, see chapters Finite Element Formulation and Dynamic Time Domain Analysis. However, the linearization also introduces an additional stiffness due to the displacement dependent loads. This stiffness is called load correction stiffness, are described in section Incremental Equilibrium Iterations. The inclusion of the load correction stiffness only has a consequence for the rate of convergence in a numerical solution procedure. Provided uniqueness of the solution, it is the equilibrium equation alone that governs the final solution. A thorough discussion regarding type of loads and the consequence of implementing the corresponding stiffness matrix is presented by Mathisen (1990).
In the current version of the program, a load correction stiffness matrix associated with the water plane stiffness (buoyancy springs) is included.
In the following, the formulation of hydrodynamic loads are described as developed for beam elements. This formulation applies for the bar elements by use of corresponding displacement functions and by omitting rotational degrees of freedom.
8.6.1. Hydrodynamic mass and damping matrices
The hydrodynamic (added) mass and damping matrices are established based on the hydrodynamic mass and damping calculated along the secant length of the deformed element.
The components of these matrices are computed according to
where \(\mathrm {x_i}\) is taken at the buoyancy center of the subelement. \(\mathrm {\Delta m_{\alpha \gamma }^h}\) and \(\mathrm {\Delta b_{\alpha \gamma }}\) are added mass and damping for the subelement. These terms are calculated according to the applied hydrodynamic load model.
If the added mass and damping forces are established according to Equation (105) the added mass for the subelement is computed by
and the damping computed by
where \(\mathrm {A_{ij}^{(2D)}}\) and \(\mathrm {B_{ij}^{(2D)}}\) are 2dimensional added mass and damping per unit length dependent on the instantaneous immersion. \(\mathrm {\Delta x_0}\) is the length of the subelement in the initial condition.
If the simplified approach is applied for the hydrodynamic load calculations, Equation (110), no hydrodynamic damping is considered. The added mass for the subelement is established according to
where \(\mathrm {C_{my}}\) and \(\mathrm {C_{mz}}\) are directional dependent added mass coefficients, and \(\mathrm {A_S}\) is submerged cross section area.
8.6.2. External hydrodynamic load vector
The hydrodynamic loads are computed along the secant length of the deformed element.
Considering the transverse load contributions the hydrodynamic load on each subelement is written
where
\(\mathrm {\boldsymbol{f_n^{FK}}}\) : FroudeKrylov force including buoyancy, computed according to Equation (105).
\(\mathrm {\boldsymbol{f_n^S}}\) : diffraction force computed according to (Eq. 7.99) or Equation (109).
\(\mathrm {\boldsymbol{f_n^D}}\) : drag force computed according to Equation (114).
The external moment around local xaxis is approximated according to
where \(\mathrm {\boldsymbol{r^B}}\) is the vector from the principal axis to the buoyancy center, see Figure 15. In the axial direction, \(\mathrm {F_x^{FK}}\) and \(\mathrm {f_x^D}\) are computed according to Equation (101) and Equation (112), respectively.
The element end forces, \(\mathrm {F_x^{FK}}\) , FroudeKrylov forces including hydrostatic pressure, are assumed to be constant over the element length. Hence, the cross sectional axial force (stress resultant) is expressed in terms of effective axial force, see section Section 1. The effect axial force is used directly in the expression for the geometric stiffness.
Based on these considerations, the components of the external hydrodynamic load vector is computed according to
where \(\mathrm {N}\) denotes the actual displacement function. \(\mathrm {N}\) and \(\mathrm {f}\) are taken at the buoyancy center \(\mathrm {x_i}\) of each subelement. \(\mathrm {\Delta x_0}\) is the initial length of the subelement.
8.6.3. Hydrodynamic stiffness matrix
The stiffness associated with water plane stiffness (buoyancy springs) is included in the stiffness matrix. The stiffness terms are limited to contributions due to translational displacement in global Zdirection and torsional rotation around local xaxis. For the computation of the hydrodynamic stiffness matrix, linear displacement functions are applied.
Thus, the element hydrodynamic stiffness matrix with regard to Zdisplacement, \(\mathrm {\boldsymbol{\tilde{k}_{ww}^{E}}}\) , is calculated directly in global system according to
where \(\mathrm {\frac{\partial \,\tilde{f}_Z^{~\textrm{Buoy}}}{\partial Z}}\) is the water plane stiffness of the subelement based on the water plane area in the global XYplane
The local torsional hydrodynamic stiffness matrix is calculated according to
Note that this formulation results in a symmetric hydrodynamic stiffness matrix.
9. Aerodynamic loads
9.1. Lift forces, drag forces, and torsional moments on elements
Airfoil sections which are not a part of a wind turbine may be subjected to aerodynamic loads based on interpolation from the static lift, drag, and moment curves, or Morisontype quadratic drag loads.
9.1.1. Morisontype aerodynamic drag loads (CHTYPE=MORI)
The drag load calculation for Morisontype aerodynamic drag is based on the relative wind velocity in the local system:
where \(\mathrm {\boldsymbol{u_r}}\) is the relative velocity at the element center, \(\mathrm {\boldsymbol{u_{wind}}}\) is the undisturbed incoming wind velocity at the element center, and \(\mathrm {\boldsymbol{\dot {v}}}\) is the structural velocity at the element center. In the local system, the xaxis is along the element.
The magnitude of the wind velocity normal to the elements is used in the calculation of the drag load:
Using the dimensional drag coefficients (\(\mathrm {C_{Dx}}\), \(\mathrm {C_{Dy}}\), \(\mathrm {C_{Dz}}\)), the aerodynamic load per unit length is calculated as:
The force per unit length is multiplied by the dry length (following the same water level integration as the wave loads) and applied in a lumped manner.
For an axisymmetric cross section, it is possible to give nondimensional aerodynamic coefficients (\(\mathrm {C_{dx}}\), \(\mathrm {C_{dy}}\), \(\mathrm {C_{dz}}\)) as input rather than dimensional coefficients:
where \(\mathrm {S_{2D}}\) is the crosssectional dry surface, \(\mathrm {B_y}\) is the projected area per unit length in the local \(\mathrm {y}\) direction, and \(\mathrm {B_z}\) is the projected area per unit length in the local \(\mathrm {z}\) direction. For a circular cross section with diameter \(\mathrm {D}\),
9.1.2. Lift forces, drag forces, and torsional moments on elements (CHTYPE=AIRF)
For airfoil sections which are not a part of a wind turbine, aerodynamic loads are based on interpolation from the static lift, drag, and moment curves.The wind velocity and structural velocity are computed at the center of each airfoil element based on the undisturbed wind speed and structural velocity at the nodes. All computations are carried out in the local airfoil coordinate system. The lift force per unit length is denoted \(\mathrm {F_L}\), the drag force per unit length is \(\mathrm {F_D}\), and the aerodynamic moment per unit length is \(\mathrm {M}\).
The relative velocity in the airfoil coordinate system is computed as:
The spanwise (\(\mathrm {V_{r_{z}}}\)) component of the velocity is assumed to be small, and an error is returned if the spanwise component is large compared to the velocity magnitude. The angle of attack (\(\mathrm {\alpha }\)) is computed based on the local \(\mathrm {V_{r_{x}}}\) and \(\mathrm {V_{r_{y}}}\) components of the velocity. The Reynolds number (Re) is computed as
where \(\mathrm {\rho _{air}}\) is the air density, \(\mathrm {c}\) is the chord length, and \(\mathrm {\nu_{air}}\) is the viscosity of air.
Lift (\(\mathrm {C_L}\)), drag (\(\mathrm {C_D}\)), and moment (\(\mathrm {C_M}\)) coefficients are then obtained by interpolation in the airfoil tables based on Re and \(\mathrm {\alpha }\). The forces and moment (about the airfoil coordinate system) per unit length are:
Then, \(\mathrm {F_L}\), \(\mathrm {F_D}\), and \(\mathrm {M}\) are multiplied by the dry length of the element, transformed into the appropriate coordinate system, and are applied to the nodes in a lumped manner. The dry length of the element is determined consistent with the wave load integration definition.
9.2. Blade element/momentum (BEM) theory for wind turbines
The aerodynamic load model used for wind turbines in RIFLEX and SIMO is based on the blade element momentum (BEM) theory. This theory combines momentum theory and blade element theory. The general principle of the theory is that forces developed locally at the airfoil, based upon empirical lift and drag coefficients, are balanced with the change in momentum of the air flowing through the rotor disk.
9.2.1. Outline of BEM calculation
The balance of airfoil forces with the change in momentum of the air is the basic step in the calculation. In the implemented calculation, the dynamic wake method is used. Only one solution sequence is performed for each time step, and time is used a relaxation parameter. In a quasistatic analysis, an iterative calculation would be performed for each time step.
An outline of the basic force balance is given below. Each wind turbine blade is divided into blade elements. For each blade element, aerodynamic airfoil coefficients (lift, drag and moment) are defined.

The calculation begins with a suitable guess for the induced velocity at each blade element. Induced velocity is defined as the difference between the velocity of air at the rotor plane and the velocity of air far upstream of the rotor plane.

The relative air velocity with respect to each blade element is calculated.

The angleofattack with respect to the airfoil is calculated.

Using airfoil coefficient data, the aerodynamic lift force, drag force, and moment are calculated using the Øye dynamic stall model.

Airfoil forces are rotated into the rotor plane coordinate system.

Annulusaverage forces and velocities are computed, for use in the momentum balance. Rotor planeaverage remote and induced velocities are computed, for use in computing the Prandtl factor.

The Prandtl factor is calculated. This factor modifies the basic momentum balance to account for a finite number of blades.

Induced velocity (the average over each annulus) is updated based upon momentum balance, using annulusaverage forces.

In cases where the rotor is yawed with respect to the remote wind, the wake skew angle is computed, and the variation of induced velocity with blade azimuth angle is computed. Based upon the newly calculated induced velocity, a guess can be made to feed the next iteration.

The calculation is so far quasistatic. The final step is the Øye dynamic wake model, which acts as a filter for the induced velocity. The inputs are the quasi static induced velocity and the induced velocity from the last time step.
9.2.2. Correction factors
A number of correction factors are applied in the BEM method.
Glauert correction
BEM theory is not valid for induction factors greater than 0.4. The Glauert correction is used for large induction factors: for \(\mathrm {a>0.4}\), the empirical thrust curve recommended by Burton et al (2011, p67) is used:
where \(\mathrm {a_2=1.0}\), \(\mathrm {C_{T2}=1.82}\), \(\mathrm {a_1=1.00.5\sqrt{C_{T2}}}\), \(\mathrm {C_{T1}=4a_1(1a_1)}\), and \(\mathrm {F}\) is the Prandtl factor Equation (141). For \(\mathrm {F=1}\), the thrust curve used in the calculations is compared to the momentum theory result shown in the figure Figure 17.
Prandtl factor
To account for the tip and hub loss on a blade due to a finite number of blades, the Prandtl factor is implemented. The Prandtl factor \(\mathrm {F}\) as a function of the radial location \(\mathrm {r}\) is computed as:
where
\(\mathrm {B}\) is the number of blades, \(\mathrm {R}\) is the outer radius of the rotor, and \(\mathrm {\phi }\) is the angle that the trailing vorticies make with the rotorplane. \(\mathrm {\phi }\) is assumed to be the same as the angle of incoming flow, neglecting tangential induced velocity. The default program behavior is to correct \(\mathrm {\phi }\) for yawed inflow.
Dynamic wake
The dynamic wake effect is the time lag in induced velocities due to the shedding and downstream convection of vorticity. Dynamic wake effects are most pronounced for heavily loaded rotors, corresponding to high induction factors (low wind speeds). This effect can be modeled by the Stig Øye dynamic inflow model, which acts as a filter for induced velocities, as shown below (Snel and Schepers, 1995). \(\mathrm {\boldsymbol{W_{qs}}}\) is the quasi static induced velocity vector, and \(\mathrm {\boldsymbol{W}}\) is the new induced velocity vector. \(\mathrm {\tau_1}\) and \(\mathrm {\tau_2}\) are time constants.
Dynamic stall
Dynamic stall refers to timedependent variations in the lift and drag coefficients of an airfoil due to changes in the angle of attack. The Stig Øye model is implemented and gives unsteady lift by filtering the trailing edge separation point with an empirical time constant (see also Hansen, 2008).
The degree of stall (\(\mathrm {f_s}\)) affects the lift as
where \(C_{L,\mathrm {inv}}\) is the lift coefficient without separation, and \(\mathrm {C_{L,fs}}\) is the fully separated lift coefficient. A static value \(\mathrm {f_{s}^{st}}\) can be found to reproduce the static airfoil data. The degree of stall is assumed to chase the static value:
where \(\mathrm {\tau}\) is a time constant. By integration:
For flow that is not fully stalled, the lift coefficient is computed as:
where \(\mathrm {\frac{dC_L}{d\alpha }}\) is computed at the fullstall limit (depending on the sign of \(\mathrm {\alpha }\)), and \(\mathrm {\alpha _0}\) is the angle of attack where the lift is zero.
For fully stalled flow:
where \(\mathrm {C_{L,qs}}\) is the quasistatic lift.
In order to initialize the dynamic stall calculation, the program identifies several important points in the lift curve: the angle of attack where there is zero lift (\(\mathrm {\alpha _0}\)), the maximum slope of the lift curve in the linear region for both positive and negative lift values (\([dC_L/d\alpha ]_{\mathrm {max}}(1)\) and \([dC_L/d\alpha ]_{\mathrm {max}}(2)\)), and the angle of attack corresponding to full separation for both positive and negative lift values (\(\mathrm {\alpha _{fs}(1)}\) and (\(\mathrm {\alpha _{fs}(2)}\))) For typical airfoils, these points are relatively easy to identify. For more complex airfoils, the user may input these values directly.
Skewed inflow
Skewed inflow may occur due to rotor tilt or a yaw angle between the rotor and the oncoming wind. Glauert developed the basic formulation for correcting the induction factor due to skewed inflow: \(\mathrm {a}\) is the induction factor, \(\mathrm {\chi}\) is the wake skew angle, and \(\mathrm {\Psi}\) is the azimuth angle that is zero at the most downwind position of the rotor.
Upwind tower influence
The tower influence (shadow) effect is implemented using potential flow theory. Since the incoming wind has to travel around the tower, the tower has an effect on the local inflow, even for an upwind turbine. The velocity deficit due to the tower is based on the 2D potential solution for constant flow around a circle, with the possibility to include corrections for the tower drag and the Bak coefficient. The coordinate system for the tower influence calculation is shown in figure Figure 19.
For an upwind wind turbine, the nondimensional influence at a location \(\mathrm {(x,y,z)}\) is evaluated as:
where \(D_{\mathrm {tow}}\) is the local tower diameter for the given \(\mathrm {z}\) level, \(\mathrm {b}\) is the Bak coefficient, and \(\mathrm {C_D}\) is the local tower drag coefficient for the given \(\mathrm {z}\) level. If the effect of tower drag is not considered, the Bak factor is set to zero and potential flow solution is applied.
In practice, the horizontal wind direction and speed are computed at the location \(\mathrm {(x,y,z)}\). This speed is multiplied by the velocity factors \(x_{\mathrm {infl}}\) and \(y_{\mathrm {infl}}\), and the modified speed is then transposed back to the initial wind direction (Moriarty and Hansen, 2005).
Note that the \(\mathrm {x}\)coordinate here is upwind and the velocity is set to zero for any point inside the tower. The tower is assumed to be approximately vertical and no effect on the vertical wind speed is included.
When the blade segment is above the tower height, the velocity is smoothened to match the free wind. The resulting velocity deficit factors is shown in the figures below. The blade is pointing upwards when the azimuth angle is 0 deg, and is closest to the tower at 180 deg. In these figures the distance from the hub to the center of the segment is 90 m, the tower top is located 6 m below the hub and the tower diameter is ranging from 9 m (azimuth=180 deg) to 6.54 m (tower top). There is no cone on the blade line and no tilt of the shaft.
Downwind tower influence
In order to model downwind wind turbines, a different method for computing the velocity deficit due to the presence of the tower is needed. The most common models, which are implemented in AeroDyn and Bladed, for example, are based on work by Powles. Powles’ model, as implemented in RIFLEX, assumes a cosinesquared shape for the tower shadow.
The cosinesquared model for the wake factor \(u_{\mathrm {wake}}\) is: $ u_\{}=
$ where \(\mathrm {d=\sqrt{x^2+y^2}}\). Together with the potential flow model, the local wind velocity is found as: $ U_\{}=(uu_\{})U_$ $ V_\{}=(vu_\{})U_$
where \(\mathrm {U_\infty}\) is the total instantaneous horizontal velocity, and \(\mathrm {u=1}\) when the cosinesquared model is in effect. In Powles’ model, the flow is assumed to not reverse: for \(uu_{\mathrm {wake}}<0\), we assume \(U_{\mathrm {local}}=0\). The computed \(U_{\mathrm {local}}\) and \(V_{\mathrm {local}}\) are then transformed back to the local global reference frame. Thus, the wake is assumed to align with the instantaneous horizontal wind vector.
The resulting velocity deficit for a downwind turbine is shown in figure Figure 22.
9.2.3. BEM Limitations
The BEM aerodynamic load module has the following limitations: 1. The code does not truly account for blade cone or large deflections 2. The vortical wake structure must be preserved. For this to be the case, the ratio of blade tip speed to wind speed must not be too low, and the majority of the blade must not be stalled. 3. The cone angle should be small. If the rotor is coned: for purposes of load calculation, define the aerodynamic analysis as if the cone angle were zero; for a lower bound on power output, define the aerodynamic analysis with the actual cone angle.
10. References
Burton, T., Jenkins, N., Sharpe, D. and Bossanyi, E. (2011): Wind Energy Handbook (Wiley)
Dalzell, J.F. (1999): A note on finite depth second order wavewave interactions. Applied Ocean Research 21 105.
Dean, R. G. & Dalrymple, R. A. (1991): Water Wave Mechanics for Engineers and Scientists (World Scientific)
Dixon, A.G., Greasted, C.A. and Salter, S.H. (1979): Wave Forces on Partially Submerged Cylinders ASCE, Journal of the Waterway, Port, Coastal and Ocean Division, Nov.
Engseth, A., Bech, A. and Larsen, C.M. (1988): Efficient Method for Analysis of Flexible Risers. Proc. of the Int. Conf. on Behaviour of Offshore Structures (BOSS), Vol. 3, pp. 13571371, Tapir, Trondheim.
Faltinsen, O.M. (1977): Water impact loads and dynamic response of horizontal cylinders in offshore structures. OTCpaper No. 2741, Texas.
Forristall, G.Z. (2000): Wave Crest Distributions: Observations and SecondOrder Theory. Journal of Physical Oceanography 30 1931.
Frank, W. (1967): Oscillation of Cylinders in or below the Free Surface of Deep Fluids Report 2375. Washington DC: Naval Ship Research and Development Center.
Fylling, I.J., Sødahl, N. and Bech, A. (1986): On the Effects of Slug Flow on the Dynamic Response of Flexible Risers and Submerged Hoses. IBC Technical Services, London.
Greenhow, M. (1987): Water entry and exit of horizontal circular cylinder. Technical Report, Marintek, Trondheim.
Hansen, M.O.L. (2008): Aerodynamics of Wind Turbines. 2nd Edition.
Hoerner, S.F. (1965): Fluiddynamic Drag. Published by the author, Washington D.C.
Hooft, J.P. (1972): Hydrodynamic Aspects of Semisubmersible platform in waves. H. Veeman on Zonen N.V., Wagningen.
Kaplan, P. and Hue, P.N. (1959): Virtual mass and slender body theory for bodies in waves. Proc. Sixth An. Conference on Fluid Mechanics, University of Texas, Texas.
Kaplan, P. and Silbert, M.N. (1976): Impact forces on platform horizontal members in the splash zone. OTCpaper No. 2498, Texas.
MacCamy, R. C. & Fuchs, R. A. (1954): Wave Diffraction on Piles: A Diffraction Theory. Beach Erosion Board; Corps of Engineers.
Marthinsen, T. and Winterstein, S.R. (1992): On the skewness of random surface waves. The Second International Offshore and Polar Engineering Conference, International Society of Offshore and Polar Engineers.
Mathisen, K.M. (1990): Large Displacement Analysis of Flexible and Rigid Systems considering DisplacementDependent Loads and Nonlinear Constraints. Doctoral Thesis, Div. of Structural Mechanics, The Norwegian Institute of Technology, Trondheim.
Mathisen, K.M. and Bergan, P.G. (1986): Nonlinear static analysis of flexible risers. The International Symposium & Exhibition, OMAE 1986, Tokyo.
Mei, C.C. (1989): The Applied Dynamics of Ocean Surface Waves. World Scientific.
Morgan, G.W. and Peret, J.W. (1975): Applied Mechanics of Marine Riser Systems. Petroleum Engineer, Oct. 1974  Oct. 1975.
Moriarty, P.J. and Hansen, A.C. (2005): AeroDyn Theory Manual. NREL/TP50036881
Ormberg, H. (1991): NonLinear Response Analysis of Floating Fish Farm Systems. Doctoral Thesis, Div. of Marine Structures, The Norwegian Institute of Technology, Trondheim.
Powles, S.R.J. (1983): The effects of tower shadow on the dynamics of a HAWT. Wind Engineering, 7:26–42, 1983.
Sharma, J. and Dean, R.G. (1981): SecondOrder Directional Seas and Associated Wave Forces. Publication of: Offshore Technology Conference.
Skjelbreia, L. and Hendrickson, J. (1960): Fifth order gravity wave theory. Proceedings of 7th Conference on Coastal Engineering, Vol. 1, The Netherlands.
Snel H., and Schepers J.G. (1995): Joint Investigation of Dynamic Inflow Effects and Implementation of an Engineering Method. Report ECNC–94107, Energy Research Centre of the Netherlands.
Sparks, C.P. (1984): The influence of tension, pressure and weight on pipe and riser deformations and stresses. Journal of Energy Resources Technology, Vol. 106, pp. 4654.
Stansberg, C.T., Gudmunstad, O.T. and Haver, S.K. (2008): Kinematics under extreme waves. Journal of Offshore Mechanics and Arctic Engineering 130 021010.
Øritsland O. (1986): Wave forces on a sphere in the splash zone. Technical report, MARINTEK, Trondheim.
Aarsnes, J.V., Lien E. and Oltedal, G. (1988): SIFICA  Theory Manual MARINTEK report: 513004.05.02, Trondheim.