4.3.7. AeroDyn Theory
This theory manual is work in progress, please refer to the AeroDyn 14 manual for more details [ad-MH05]. Many changes have occurred since AeroDyn 14 (e.g. BEM formulation, coordinate system used in the BEM equations, dynamic stall, dynamic BEM), but these changes are not yet documented here.
4.3.7.1. Steady BEM
The steady blade element momentum (BEM) equations are solved as a constrained equation, and the formulation follows the description from Ning [ad-Nin14].
4.3.7.2. Dynamic BEM Theory (DBEMT)
Two equivalent versions of Oye’s dynamic inflow model are implemented in AeroDyn.
The first one uses discrete time, it can be used with the constant-tau1 model
(DBEMT_Mod=1) or the varying-tau1 model (DBEMT_Mod=2), but it cannot be used for linearization.
The second version uses a continuous-time state-space formulation (DBEMT_Mod=1), it assumes a constant-tau1, and can be used for linearization.
For a same value of \(\tau_1\), the discrete-time and continuous-time formulations returns exactly the same results.
Oye’s dynamic inflow model consists of two first-order differential equations (see [ad-Bra17]):
where \(\boldsymbol{W}\) is the dynamic induction vector at the rotor (at a given blade position and radial position), \(\boldsymbol{W}_\text{qs}\) is the quasi-steady induction, \(\boldsymbol{W}_\text{int}\) is an intermediate value coupling the quasi-steady and the actual inductions (may be discontinuous if the quasi-steady indution is discontinuous). and \((\dot{\ })\) represents the time derivative. The coupling constant \(k\), with values between 0 and 1, is usually chosen as \(k=0.6\). Oye’s dynamic inflow model relies on two time constants, \(\tau_1\) and \(\tau_2\) :
where \(R\) is the rotor radius, \(\overline{U}_0\) is the average wind speed over the rotor, \(\overline{a}\) is the average axial induction over the rotor, and \(r\) is the radial position along the blade.
For DBEMT_Mod=1 or DBEMT_Mod=3, the user needs to provide the value of \(\tau_1\).
The continuous-time state-space formulation of the dynamic inflow model (DBEMT_Mod=3) was derived in [ad-BJP+22].
where \(\boldsymbol{I}_2\) is the 2x2 identity matrix, \(\boldsymbol{W}_\text{red}\) is the reduced induction which is a continuous, scaled, and lagged version of the quasi-steady induction, defined as:
The discrete-time version of the model is documented in the unpublished manual of DBEMT. The current discrete-time formulation is complex and in the future it can be simplified by using \(\boldsymbol{W}_\text{red}\).
4.3.7.3. Tower influence models
AeroDyn can model the influence of the tower on the flow reaching the blades through two
superimposable effects: a potential-flow disturbance of the flow around the tower
(TwrPotent) and a downstream shadow (wake) velocity deficit (TwrShadow). Both are
evaluated at each blade node from the position of the node relative to the nearest point on
the tower, expressed in a local tower reference frame centered on the tower axis:
\(\overline{x}\) is the downstream (streamwise) coordinate, \(\overline{y}\) the
lateral coordinate, and \(\overline{z}\) the coordinate along the tower axis, each
normalized by the local tower radius. The resulting disturbance-velocity fractions are scaled by
\(W_\text{tower}\), the local incoming wind speed component normal to the tower axis, and
applied in the plane normal to the tower. \(C_d\) is the local tower drag coefficient and
\(TI\) (Eames model only) is the local turbulence intensity at the tower node. It is
convenient to define \(\overline{r} = \sqrt{ \overline{x}^2 + \overline{y}^2 }\).
4.3.7.3.1. Tower potential-flow model
The baseline potential-flow model (TwrPotent=1) represents the tower cross section as a two-dimensional cylinder in potential flow (a doublet), giving the streamwise (\(u_{TwrPotent}\)) and lateral (\(v_{TwrPotent}\)) disturbance-velocity fractions:
The Bak correction (TwrPotent=2), following Bak, Madsen, and Johansen (2001), adds a drag-induced source term and offsets the streamwise coordinate by \(+0.1\) (a fixed empirical constant) to better represent the near wake. Writing \(\overline{x}' = \overline{x} + 0.1\) and \(\overline{r}' = \sqrt{ \overline{x}'^2 + \overline{y}^2 }\),
Alongside the tower, the potential-flow disturbance velocity is applied in full. Beyond the tower ends the influence is tapered and cut off; this end treatment, which is shared with the shadow model, is described in Section 4.3.7.3.3.
4.3.7.3.2. Tower shadow models
Powles tower shadow model (TwrShadow=1) is given by:
where \(\overline{r} = \sqrt{ \overline{x}^2 + \overline{y}^2 }\).
Eames tower shadow model (TwrShadow=2) is given by:
where \(TI\) is the turbulence intensity at the tower node.
To avoid excessive flow reversal behind the tower, the shadow deficit fraction is limited to \(u_{TwrShadow} \ge -0.5\).
The potential-flow and shadow contributions are superimposed and scaled by \(W_\text{tower}\) to give the velocity perturbation at the blade node in the local tower frame:
This perturbation \((v_x, v_y, 0)\) is rotated into the earth-fixed frame and added to the free-stream (undisturbed) inflow velocity to obtain the disturbed inflow velocity used by the blade-element calculations.
4.3.7.3.3. Tower ends and exclusion zones
Both tower-influence effects use the position of the blade node relative to the nearest point on the tower. Alongside the tower this nearest point is the orthogonal projection of the blade node onto the tower axis, for which the axial coordinate is \(\overline{z} = 0\). When a blade node lies axially beyond a tower end, the nearest point becomes the tower end itself, and \(\overline{z}\) is then the axial distance from that end, normalized by the local tower radius.
For each blade node, the tower clearance is computed as \(c = \lVert \mathbf{r} \rVert - R\),
where \(\mathbf{r}\) is the vector from the nearest tower point to the blade node and
\(R = \tfrac{1}{2}\) TwrDiam is the local tower radius. The disturbance is suppressed both
very close to the tower (\(c \le 0.01\,\) TwrDiam) and far from it (\(c > 20\,\) TwrDiam,
a far-field cutoff). Because the clearance beyond a tower end is measured to the end point, the
surface of minimum clearance for the tower influence models at member ends becomes a hemispherical cap
closing off the cylinder. The near-tower exclusion zone is therefore a capsule — a cylinder capped
by hemispheres at both ends — rather than a bare cylinder.
Beyond a tower end, the disturbance is faded out over one tower radius using a cosine-squared axial taper. For \(|\overline{z}| < 1\) the in-plane coordinates are divided by \(\cos\!\left(\tfrac{\pi}{2}\overline{z}\right)\),
The taper is applied purely through this coordinate scaling; how it attenuates each effect depends on how that effect depends on \(\overline{x}\) and \(\overline{y}\). Because the potential-flow disturbance velocity scales as \(1/\overline{r}^{\,2} = 1/(\overline{x}^2 + \overline{y}^2)\), dividing \(\overline{x}\) and \(\overline{y}\) in this way multiplies the potential-flow velocity by exactly \(\cos^2\!\left(\tfrac{\pi}{2}\overline{z}\right)\): it decreases smoothly from its full value at the end plane (\(\overline{z} = 0\)) to zero at \(\overline{z} = 1\) with zero slope there (\(C^1\)-continuous). The shadow deficit is tapered by the same coordinate scaling, but its dependence on \(\overline{x}\) and \(\overline{y}\) is different, so the attenuation is not a clean \(\cos^2\) factor (for the Powles model the deficit scales roughly as \(\sqrt{\cos}\) with an additional shift in its lateral argument, and for the Eames model it scales linearly in \(\cos\) with the lateral profile unchanged). Both shadow models nonetheless fade to zero as \(\overline{z} \rightarrow 1\). No tower influence is applied for \(|\overline{z}| \ge 1\).
Note
The tower shadow deficit is currently tapered near the tower ends only as a byproduct of the coordinate scaling above, which does not give a principled axial wake taper. The end treatment of the tower shadow is expected to be improved in a future release.
Furthermore, the far-field cut-off beyond 20 tower diameters of clearance is expected to be revised for the tower shadow effect in a future release. This revision would likely be made in conjunction with the addition of a low-pass filter for the tower inflow velocity used to compute the shadow deficit and direction. This is to prevent erratic behavior of the shadow region in response to instantaneous inflow velocity fluctuations.
4.3.7.4. Tower drag loads
AeroDyn can apply an aerodynamic drag load to the tower itself when tower aerodynamics are enabled (TwrAero=True). The load is a cross-flow (Morison-type) drag evaluated independently at each tower node.
At tower node \(j\) the relative wind is \(\mathbf{V}_\text{rel} = \mathbf{V}_\text{inflow} - \mathbf{V}_\text{motion}\), the difference between the local unperturbed inflow velocity and the tower structural velocity of the node. Only the component of \(\mathbf{V}_\text{rel}\) in the plane normal to the tower axis produces drag; denote this transverse relative-wind vector \(\mathbf{V}_\perp\) and its magnitude \(W_\text{tower} = \lVert \mathbf{V}_\perp \rVert\). The drag force per unit length is
where \(\rho\) is the air (or water, for MHK) density, \(C_d\) = TwrCd is the local
tower drag coefficient, and \(D\) = TwrDiam is the local tower diameter. The force acts
in the direction of the transverse relative wind, the axial (along-tower) component is zero, and
no moment is applied. This per-unit-length load is distributed along the tower line mesh and
later mapped to the ElastoDyn tower structural mesh in a coupled simulation.
4.3.7.5. Generalized support-structure influence models
The generalized support structure (GS) extends the tower influence models of
Section 4.3.7.3 to an arbitrary assembly of slender cylindrical members (for
example the columns and braces of a jacket, tripod, or floating platform). Each member is
treated exactly like the tower: at each blade node the potential-flow disturbance
(GSPotent) and the downstream shadow deficit (GSShadow) are evaluated in a local
member frame from the position of the blade node relative to the nearest point on that
member, using the same normalized coordinates \((\overline{x}, \overline{y},
\overline{z})\), the same baseline and Bak potential-flow expressions (see
Section 4.3.7.3.1), the same Powles and Eames shadow expressions (see
Section 4.3.7.3.2), and the same end handling — clearance-based exclusion capsule and
axial taper (see Section 4.3.7.3.3). The member diameter and drag
coefficient play the roles that TwrDiam and \(C_d\) play for the tower. Because the
support members generally meet at joints, the axial taper is applied at every member end
(free tip or junction), whereas for the single tower it is only ever needed at the two free
ends.
The distinguishing feature of the GS model is how the contributions of the individual members are combined at a blade node. For each blade node the potential-flow and shadow contributions are accumulated separately in the earth-fixed frame and combined by different rules, and the resulting GS disturbance is then superimposed on the tower disturbance by simple addition (there is no cross-blend between the tower field and the GS field).
Note
The generalized support-structure influence on the inflow (GSPotent and GSShadow)
is applied to the disturbed inflow at the blade nodes — and hence to the unsteady airfoil
aerodynamics — for every wake model, including the free-vortex-wake model OLAF
(Wake_Mod = 3). The one exception is the wake-convection velocity used to transport the
OLAF free vortex wake: only the tower influence is applied there, and the GS influence on
the convected wake is not yet included. This is planned to be addressed in a future release.
4.3.7.5.1. Combining the potential-flow contributions
A naive superposition (summation) of the per-member potential-flow solutions is not appropriate. Each member’s field is the potential-flow solution for that member in isolation, which already enforces the non-penetration boundary condition on that member’s surface. Adding two such fields violates the combined boundary condition and over-counts the disturbance where members are close together — most visibly at a joint, where two members meeting at a point would each contribute a full near-field doublet and roughly double the true disturbance of the single connected body. Simple approaches to remove this double count (for example detecting collinear members that continue through a joint) were found to be fragile and not general: real assemblies present an unbounded variety of geometries (slightly angled continuations, a single member branching into two, and so on) that no finite set of topology rules covers robustly.
The adopted combination is instead an influence-weighted partition of unity. Let \(\mathbf{v}_i\) be the potential-flow contribution of member \(i\) at the blade node (expressed in the earth-fixed frame). The combined potential-flow disturbance is
with blend exponent \(p = 6\). Because the weights are non-negative and sum to one, the blended magnitude never exceeds \(\max_i \lVert \mathbf{v}_i \rVert\), so the combination can never inflate the disturbance and the joint double-counting is eliminated by construction. The scheme has several convenient properties:
It reduces exactly to the single-member (and hence tower) result when only one member contributes.
It requires no knowledge of the structure topology: collinear, angled, and branching configurations are all handled through the field magnitudes alone. At a collinear joint, the two members carry almost identical fields, so any split of the weights returns essentially the single-cylinder value; at a corner the blend hands off smoothly between the members.
Weighting by the influence magnitude \(\lVert \mathbf{v}_i \rVert\) rather than by proximity ensures that a member whose field has tapered to zero (for instance at a blade node axially above the end of the member by just over one member radius) receives essentially zero weight and cannot blank out the field of a nearby member slightly further away.
The exponent \(p\) controls the sharpness of the handoff: \(p \rightarrow \infty\) recovers a hard “nearest/strongest body only” selection, while a finite \(p\) smooths the transition. An even integer is used so that \(w_i = (\mathbf{v}_i \cdot \mathbf{v}_i)^{p/2}\) is a polynomial in the velocity components and therefore smooth everywhere. The value \(p = 6\) is the smallest even integer that keeps the small residual dip at a collinear same-diameter joint (an artifact of blending a full field against its cosine-tapered neighbor) below about 5 %, while keeping the handoff gradients modest.
The per-member end taper of Section 4.3.7.3.3 is still applied before the blend: the taper makes each finite member’s field die away beyond its physical extent, and the blend only decides which member dominates where several overlap.
4.3.7.5.2. Combining the shadow contributions
The partition of unity is deliberately not used for the shadow deficit. Unlike the potential-flow disturbance, wake deficits physically stack to some extent: two overlapping wakes remove more momentum than one. The correct combined deficit therefore lies somewhere between the single-member value (which the partition of unity would return) and the linear sum (which over-counts). The shadow contributions are combined by a root-sum-square of the individual deficit magnitudes, applied along the direction of their vector sum. Writing \(\mathbf{s}_i\) for member \(i\)’s shadow contribution in the earth-fixed frame and \(\mathbf{s}_\text{sum} = \sum_i \mathbf{s}_i\),
When only one member contributes, this reduces exactly to that member’s deficit.
Note
As with the tower shadow, the GS shadow wake is directed along the wind projected into the plane normal to the member axis rather than along the true (earth-fixed) wind. For a member strongly raked into or away from the wind, this tilts the modeled wake up into the sky or down toward the ground instead of keeping it aligned with the incoming flow, which is unphysical. The current formulation is adequate for near-vertical members and mirrors the established tower model, but the shadow model is expected to be improved in a future release to advect the wake along the (ideally low-pass filtered) incident wind direction. The potential-flow part is a near-field kinematic effect and is correctly resolved in the member-normal plane, so this change would affect only the shadow model.
4.3.7.6. Generalized support-structure drag loads
The generalized support structure carries the same cross-flow drag load as the tower (Section 4.3.7.4), applied member by member when GS aerodynamics are enabled (GSAero=True). For a member element with unit axial vector \(\hat{\mathbf{k}}\), the transverse relative wind at a node is obtained by removing the along-member component of the relative wind,
and the drag force per unit length takes the same form as for the tower, with the member
diameter \(D = 2R\) and member drag coefficient \(C_d\) in place of TwrDiam and
TwrCd:
The only procedural difference from the tower is the mesh on which the load is returned. The tower load is a distributed (per-unit-length) load on a line mesh, whereas the GS load mesh is a point mesh. The distributed member drag is therefore lumped to the element end nodes: each element of length \(\Delta l\) contributes half of its integrated drag to each of its two end nodes, and the contributions of the elements meeting at a shared node are summed. As with the tower, no moment is applied. The lumped nodal forces are mapped to the SubDyn structural mesh in a coupled simulation.
Note
The GS model provides only the aerodynamic/hydrodynamic drag load on the support members.
The other load components relevant to MHK simulations — buoyancy, added mass, and fluid
inertia (see Section 4.3.7.7 and Section 4.3.7.8) — are not computed for
the generalized support structure and should instead be modeled in HydroDyn. To avoid
double-counting, the support-structure drag should be modeled in either AeroDyn (via the GS
drag load) or HydroDyn, but not both. The GS influence on the rotor inflow
(Section 4.3.7.5) is a separate, flow-disturbance effect and can always be included in
AeroDyn regardless of where the support-structure drag is modeled. If the user chooses to model
the drag force on the support structure in HydroDyn, the GS drag in AeroDyn should be disabled
by setting GSAero=False. However, in this case, the user might still want to set the GS drag
coefficients appropriately for the GS shadow model or the GS potential-flow model with Bak
correction if enabled.
4.3.7.7. Buoyancy
When a solid object is submerged in a fluid, it experiences a net force, buoyancy, from the hydrostatic fluid pressure acting on its surface. This force can often be neglected in less dense fluids, such as air, but can be significant in denser fluids, such as water. To capture the effects of this force on MHK turbines, buoyant loads are calculated for the turbine blades, tower, hub, and nacelle. Marine growth is neglected for all components. Section 4.3.7.7.1 - Section 4.3.7.7.3 detail the coordinate systems and blade, tower, hub, and nacelle buoyancy calculations.
4.3.7.7.1. Coordinate Systems
The buoyant force acting on an element depends on its instantaneous orientation and depth. The orientation is defined by heading and inclination angles, which are calculated for each element at every time step. Total water depth is defined by the user, relative to the still water level (or relative to the mean sea level when running AeroDyn in standalone mode with the AeroDyn driver). The instantaneous depth of each element is based on its position in global coordinates at each time step.
4.3.7.7.2. Blades and Tower
To allow for an efficient analytical solution, the blades and tower are modeled as tapered cylinders. The cross-sectional area of the tapered cylinders is set equal to the blade or tower cross-sectional area. Loads are estimated by breaking the blade or tower into elements of a given length and integrating the hydrostatic pressure over the wetted area of each element. For the blades, loads are applied at a user-specified center of buoyancy. For the tower, loads are applied at the centerline. When applicable, end effects are accounted for by calculating the fluid pressure on the exposed axial face of the element. The tower is assumed to be either embedded into the seabed or attached to another support structure member, such that no end effects at the tower base are needed. For MHK turbines with a support structure (i.e., any structure other than a simple tower embedded in the seabed), it is currently recommended to model the entire support structure, including the tower, in HydroDyn. Future releases will include the ability to neglect fluid loads at the interface between a tower modeled in AeroDyn and a platform modeled in HydroDyn.
The buoyancy calculation for the blades and tower is completed according to the following steps:
Calculate parameters related to element geometry that do not change with time
Check that no elements cross the free surface or go beneath the seabed
Calculate the instantaneous orientation and depth of each element
Integrate hydrostatic fluid pressure over the wetted surface of each element and express as a force acting at the center of buoyancy
For blades, calculate the buoyant force on the axial face of the blade root and tip; add the tip force to the adjacent element and store the root force
For the tower, calculate and store the buoyant force on the axial face of the tower top
Move buoyant loads from the center of buoyancy to the aerodynamic center
Express buoyant loads in the form expected by OpenFAST
Add buoyant loads to aerodynamic loads
Although the blade and tower buoyant loads are not based on volume, the volumes of these components are written to the AeroDyn summary file for reference. The blade and tower volumes are calculated by summing the volume of each element, assumed to be a tapered cylinder. The volume of a single element \((V_{elem})\) is given by:
where \(r_i\) is the element radius at node \(i\), \(r_{i+1}\) is the element radius at node \(i+1\), and \(dl\) is the element length.
4.3.7.7.3. Hub and Nacelle
The hub and nacelle are treated as separate components. The buoyant force is determined by the volume of either the hub or nacelle and applied at its user-specified center of buoyancy. Corrections are made to account for the joints between the hub and blades and the nacelle and tower, as the joint locations are not exposed to fluid pressure. No correction is made for the joint between the hub and nacelle.
The buoyancy calculation for the hub and nacelle is completed according to the following steps:
Check that the component does not cross the free surface or go beneath the seabed
Calculate the instantaneous depth of the component
Calculate the buoyant force from the volume of the component
Move buoyant loads from the center of buoyancy to the aerodynamic center
For the hub, correct loads to account for the joints with each blade
For the nacelle, correct loads to account for the joint with the tower
4.3.7.8. Added Mass and Fluid Inertia
Added mass loads are caused by body and fluid accelerations. These forces can often be neglected in less dense fluids, such as air, but can be significant in denser fluids, such as water. To capture the effects of these forces on MHK turbines, added mass and fluid inertia loads are calculated for the turbine blades and tower. Per-unit-length loads are estimated at each blade or tower node by calculating the added mass and fluid inertia forces according to the appropriate terms from Morison’s equation. The resulting loads are summed with the previously calculated hydrodynamic and/or buoyant per-unit-length loads. Loads for the blades are applied at the aerodynamic center. Loads for the tower are applied at the centerline. Marine growth and end effects are neglected, and members are not allowed to cross the free surface (i.e., members are always fully submerged). Ballast is not considered. Nodes do not need to be uniformly spaced, and axial loads are neglected. The tower is assumed to be axisymmetric (with the same coefficients used in both transverse directions), but the blade is not (with different coefficients normal and tangential to the chord, as well as an added mass coefficient for pitch).
4.3.7.8.1. Morison’s Equation
Added mass and fluid inertia loads are calculated according to the appropriate terms from Morison’s equation. The added mass force is given as
where \(\rho\) is the fluid density, \(C_a\) is the added mass coefficient, \(V\) is the element volume, \(\dot{u}\) is the fluid acceleration, and \(\dot{v}\) is the body acceleration.
The fluid inertia force is given as
where \(C_p\) is the dynamic pressure coefficient.
The fluid density and added mass and dynamic pressure coefficients are user-specified. Added mass and fluid
inertia loads can be turned off by setting the relevant coefficients to zero. Additional information about calculating added mass coefficients can be
found in Section 4.3 (“Determination of Added Mass Coefficients for Floating Hydrokinetic Turbine Blades using Computational Fluid Dynamics”).
The body and fluid accelerations are calculated internally and passed to AeroDyn. Body accelerations are available from the structural solver (or driver),
and fluid accelerations are calculated based on the inflow velocity time series. Added mass and fluid inertia loads are calculated as per-unit-length within
AeroDyn. Therefore, \(V\) is taken as the cross-sectional area at the node of interest. For the blades, the reference cross-sectional area for the normal
and tangential terms is chord*thickness (\(ct\)). This is expressed as \((c^2)(t/c)\), where \(t/c\) (i.e., t_c) is specified
in the AeroDyn blade input file and cannot be less than 0. For the tower, the reference cross-sectional area is \(\pi r^2\) where \(r\)
is calculated as (0.5 TwrDiam). The normalization for the BlCpn, BlCpt, BlCan, and BlCat coefficients should be \(\rho ct\);
the normalization for the BlCam coefficient should be \((1/12)\rho ct(c^2+t^2)\); and the normalization for the TwrCp and TwrCa coefficients should
be \(\rho\pi(0.5\) TwrDiam) \(^2\).
4.3.7.8.2. Blade Added Mass and Fluid Inertia
Added mass and fluid inertia loads are calculated for the normal-to-chord, tangential-to-chord, and pitch directions in the blade coordinate system. The following coefficients are defined by the user in the AeroDyn blade input file:
BlCpnspecifies the blade normal-to-chord dynamic pressure coefficient; to neglect normal-to-chord fluid inertia loads on the blade, setBlCpnto 0BlCptspecifies the blade tangential-to-chord dynamic pressure coefficient; to neglect tangential-to-chord fluid inertia loads on the blade, setBlCptto 0BlCanspecifies the blade normal-to-chord added mass coefficient, cannot be less than 0; to neglect normal-to-chord added mass loads on the blade, setBlCanto 0BlCatspecifies the blade tangential-to-chord added mass coefficient, cannot be less than 0; to neglect tangential-to-chord added mass loads on the blade, setBlCatto 0BlCamspecifies the blade pitch added mass coefficient, cannot be less than 0; to neglect pitch added mass loads on the blade, setBlCamto 0
4.3.7.8.3. Tower Added Mass and Fluid Inertia
Added mass and fluid inertia loads are calculated for the transverse direction in the tower coordinate system. The following coefficients are defined by the user in the AeroDyn primary input file:
TwrCpspecifies the tower transverse dynamic pressure coefficient; to neglect fluid inertia loads on the tower, setTwrCpto 0TwrCaspecifies the tower transverse added mass coefficient, cannot be less than 0; to neglect added mass loads on the tower, setTwrCato 0