diff --git a/docs/source/user/aerodyn/examples/ad_primary_example.dat b/docs/source/user/aerodyn/examples/ad_primary_example.dat index 9199db451f..76f7cf97a9 100644 --- a/docs/source/user/aerodyn/examples/ad_primary_example.dat +++ b/docs/source/user/aerodyn/examples/ad_primary_example.dat @@ -7,6 +7,9 @@ False Echo - Echo the input to ".AD.ech"? (flag 1 TwrPotent - Type tower influence on wind based on potential flow around the tower (switch) {0=none, 1=baseline potential flow, 2=potential flow with Bak correction} 0 TwrShadow - Calculate tower influence on wind based on downstream tower shadow (switch) {0=none, 1=Powles model, 2=Eames model} False TwrAero - Calculate tower aerodynamic loads? (flag) + 0 GSPotent - Type generalized support (GS) structure influence on wind based on potential flow around the GS members (switch) {0=none, 1=baseline potential flow, 2=potential flow with Bak correction} + 0 GSShadow - Calculate GS influence on wind based on downstream shadow (switch) {0=none, 1=Powles model, 2=Eames model} +False GSAero - Calculate GS aerodynamic loads? (flag) False CavitCheck - Perform cavitation check? (flag) [UA_Mod must be 0 when CavitCheck=true] False NacelleDrag - Include Nacelle Drag effects? (flag) False CompAA - Flag to compute AeroAcoustics calculation [used only when Wake_Mod = 1 or 2] @@ -94,6 +97,14 @@ TwrElev TwrDiam TwrCd TwrTI TwrCb TwrCp 4.0000000E+01 5.0000000E+00 0.0000000E+00 1.0000000E-01 0.0000000E+00 0.0000000E+00 0.0000000E+00 6.0000000E+01 4.5000000E+00 0.0000000E+00 1.0000000E-01 0.0000000E+00 0.0000000E+00 0.0000000E+00 8.0000000E+01 4.0000000E+00 0.0000000E+00 1.0000000E-01 0.0000000E+00 0.0000000E+00 0.0000000E+00 +====== General support structure joints =========================================================== +0 NumGSJoints +GSJointID GSJointxi GSJointyi GSJointzi +(-) (m) (m) (m) +====== General support structure members ========================================================== +0 NumGSMembers +GSMemberID GSMJointID1 GSMJointID2 GSMDia1 GSMDia2 GSMCd1 GSMCd2 GSMTI1 GSMTI2 GSMDiv +(-) (-) (-) (m) (m) (-) (-) (-) (-) (m) ====== Outputs ==================================================================================== True SumPrint - Generate a summary file listing input options and interpolated properties to ".AD.sum"? (flag) 4 NBlOuts - Number of blade node outputs [0 - 9] (-) diff --git a/docs/source/user/aerodyn/input.rst b/docs/source/user/aerodyn/input.rst index 70b2565ba4..892e6c86d7 100644 --- a/docs/source/user/aerodyn/input.rst +++ b/docs/source/user/aerodyn/input.rst @@ -60,6 +60,8 @@ your use, but is not used by the software. +.. _ad_general_options: + General Options ~~~~~~~~~~~~~~~ @@ -115,6 +117,20 @@ enabled, the two influences will be superimposed. Set the ``TwrAero`` flag to TRUE to calculate fluid drag loads on the tower or FALSE to disable these effects. +``GSPotent``, ``GSShadow``, and ``GSAero`` are the counterparts of +``TwrPotent``, ``TwrShadow``, and ``TwrAero`` for the generalized support +structure (see :numref:`ad_gen_support`). Set ``GSPotent`` to 0 to disable the +potential-flow influence of the support members on the flow local to the blade, +1 to enable the standard potential-flow model, or 2 to include the Bak +correction. Set ``GSShadow`` to 0 to disable the downstream shadow model, 1 to +enable the Powles model, or 2 to use the Eames model. When both potential flow +and shadow are enabled, the two influences are superimposed, and the combined +support-structure influence is superimposed on the tower influence. Set the +``GSAero`` flag to TRUE to calculate fluid drag loads on the support members or +FALSE to disable them. The joint and member geometry used by these models is +defined in the GENERAL SUPPORT STRUCTURE sections described in +:numref:`ad_gen_support`. + During linearization analyses with AeroDyn coupled OpenFAST and BEM enabled (``Wake_Mod = 1``), set the ``DBEMT_Mod=-1`` to employ frozen-wake assumptions @@ -270,8 +286,10 @@ The velocity is averaged within this sector by attributing different weighting a **SectAvgPsiFwd** Forward azimuth (in degrees) relative to the blade azimuth where the sector ends. Must be positive. [used only when SectAvg=True]. Default is 60 deg. - - +.. note:: + The tower influence (potential flow and shadow) is included in the sector-averaged + inflow velocity. The generalized support structure influence is currently not accounted + for in the sector average. Dynamic Wake / Dynamic inflow model @@ -520,6 +538,100 @@ tower, set ``TwrCb`` to 0. To neglect added mass loads on the tower, set ``TwrCa`` to 0. To neglect fluid inertia loads on the tower, set ``TwrCp`` to 0. See :numref:`ad_tower_geom`. +.. _ad_gen_support: + +Generalized Support Structure +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The generalized support structure (GS) extends the tower-influence, tower-shadow, +and tower-drag models to an arbitrary assembly of slender cylindrical members +(for example a jacket, tripod, or truss). Each member disturbs the flow reaching +the blades in the same way the tower does, and, when ``GSAero = TRUE``, the members +also experience fluid drag loads. The GS model is entirely independent of the +single tower defined in the TOWER INFLUENCE AND AERODYNAMICS section: a model may +use the tower, the generalized support structure, both, or neither. In a coupled +OpenFAST simulation that includes SubDyn, the AeroDyn generalized support structure +is coupled to the SubDyn structural model. In contrast, the AeroDyn tower is coupled +to the vertical tower modeled in ElastoDyn or the Simplified ElastoDyn module. With +multirotor systems, AeroDyn can include multiple vertical towers with each tower +linked to a separate ElastoDyn instance; however, there can only be one shared +generalized support structure interfaced with the SubDyn model of the shared substructure. +If the coupled OpenFAST simulation does not include SubDyn, the generalized support +structure is simply attached to the ElastoDyn platform point. Note that in this case, +the generalized support structure is effectively rigid, and it is not possible to have +multiple rotors, which always require SubDyn. + +The two GS tables in the AeroDyn primary input file are always read, but their +contents are only used when at least one of ``GSPotent`` > 0, ``GSShadow`` > 0, +or ``GSAero = TRUE`` (see :numref:`ad_general_options`). When the feature is not +used, set ``NumGSJoints`` and ``NumGSMembers`` to 0. All rotors will share the same +generalized support structure definition. Therefore, the generalized support structure +input tables described below should appear exactly once in the AeroDyn input file. +They should not be duplicated for each rotor. + +.. note:: + In a coupled OpenFAST simulation, the motion of the generalized support structure + is provided by SubDyn (through its ``Y3Mesh``). The GS joints and members defined + here therefore represent members that are also modeled structurally in SubDyn, and + their geometry should ideally coincide with the corresponding SubDyn nodes so that + the motion and load transfer are meaningful. In the standalone AeroDyn driver, there + is no SubDyn coupling; the generalized support structure is held fixed at the joint + positions entered below, so the driver effectively supports a rigid (non-moving) + support structure. A flexible or moving support structure requires a coupled OpenFAST + simulation with SubDyn. + +**Coordinate convention.** GS joint coordinates are given in the AeroDyn global inertial +frame (z pointing up). For wind turbines, enter z-coordinates measured from the ground +(land-based) or the mean sea level (offshore). The convention for MHK turbines depends on +whether the turbine is fixed-bottom or floating. For floating MHK turbines (``MHK = 2``), +enter z-coordinates measured from the mean sea level (MSL). For fixed-bottom MHK turbines +(``MHK = 1``), enter z-coordinates measured from the seabed (``z = 0`` at the seabed); +AeroDyn internally shifts these by the water depth so that all internal calculations use +MSL as the datum. Therefore, for wind and floating MHK turbines, the GS joint coordinates +are the same as the SubDyn node coordinates, while for fixed-bottom MHK turbines, the GS +joint coordinates are the SubDyn node coordinates summed with the (positive) water depth. + +Generalized support structure joints +------------------------------------ + +``NumGSJoints`` is the number of joints that define the endpoints of the support +members and must be greater than or equal to two (or 0 to disable the feature). +It determines the number of rows in the joint table that follows (after two header +lines). For each joint, ``GSJointID`` is a unique user-specified integer identifier, +and ``GSJointXi``, ``GSJointYi``, and ``GSJointZi`` are the x-, y-, and z-coordinates +of the joint in the global inertial frame, following the coordinate convention above. +Every joint listed here must be referenced by at least one member; otherwise AeroDyn +aborts with an error. + +Generalized support structure members +------------------------------------- + +``NumGSMembers`` is the number of members and determines the number of rows in the +member table that follows (after two header lines). Each member is a straight, +tapered cylinder connecting two joints, with properties that vary linearly between +its two ends: + +- ``GSMemberID`` — unique user-specified integer identifier for the member. +- ``GSMJointID1`` and ``GSMJointID2`` — the ``GSJointID`` values of the member's two + end joints. +- ``GSMDia1`` and ``GSMDia2`` — member diameter (m) at joint 1 and joint 2. Used by + the potential-flow and shadow models. +- ``GSMCd1`` and ``GSMCd2`` — drag coefficient (-) at joint 1 and joint 2. Used by the + shadow models (``GSShadow`` > 0), by the Bak potential-flow correction + (``GSPotent = 2``), and by the drag calculation (``GSAero = TRUE``). +- ``GSMTI1`` and ``GSMTI2`` — turbulence intensity (-) at joint 1 and joint 2, entered + as a fraction (not a percentage). Used only by the Eames shadow model + (``GSShadow = 2``). +- ``GSMDiv`` — target element (division) length (m) along the member. AeroDyn divides + the member into ``ceiling(memberLength / GSMDiv)`` equal elements and inserts the + corresponding interior analysis nodes. Smaller values give finer resolution and + higher computational cost. It is functionally consistent with ``MDivSize`` in SubDyn + and HydroDyn. + +Diameter, drag coefficient, and turbulence intensity are linearly interpolated from +the joint-1 value to the joint-2 value along each member, and the analysis nodes are +placed at the joints and at the interior division points. + .. _AD-Outputs: Outputs @@ -536,6 +648,13 @@ used in the airfoil tables. This allows the user to check what values are being used in case the code has computed the parameters without user input. +When the generalized support structure is active, AeroDyn also writes a +``.GS.sum`` file listing the selected influence models and the +generated node table (joints and interior member nodes with their interpolated +radius, drag coefficient, and turbulence intensity). The blade-node clearance +outputs (``B#N#Clrnc``) report the minimum distance from each blade node to the +nearest tower or generalized-support member, whichever is closer. + AeroDyn can output aerodynamic and kinematic quantities at up to nine nodes specified along the tower and up to nine nodes along each blade. For outputs at every blade node, see :numref:`AD-Nodal-Outputs`. diff --git a/docs/source/user/aerodyn/theory.rst b/docs/source/user/aerodyn/theory.rst index 606804a236..1a939dbec9 100644 --- a/docs/source/user/aerodyn/theory.rst +++ b/docs/source/user/aerodyn/theory.rst @@ -102,16 +102,65 @@ The current discrete-time formulation is complex and in the future it can be sim +.. _AD_twr_influence: + +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: +:math:`\overline{x}` is the downstream (streamwise) coordinate, :math:`\overline{y}` the +lateral coordinate, and :math:`\overline{z}` the coordinate along the tower axis, each +normalized by the local tower radius. The resulting disturbance-velocity fractions are scaled by +:math:`W_\text{tower}`, the local incoming wind speed component normal to the tower axis, and +applied in the plane normal to the tower. :math:`C_d` is the local tower drag coefficient and +:math:`TI` (Eames model only) is the local turbulence intensity at the tower node. It is +convenient to define :math:`\overline{r} = \sqrt{ \overline{x}^2 + \overline{y}^2 }`. + +.. _AD_twr_potent: + +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 +(:math:`u_{TwrPotent}`) and lateral (:math:`v_{TwrPotent}`) disturbance-velocity fractions: + +.. math:: + u_{TwrPotent} = \frac{-\overline{x}^2 + \overline{y}^2}{\overline{r}^4}, + \qquad + v_{TwrPotent} = \frac{-2\,\overline{x}\,\overline{y}}{\overline{r}^4} + +The Bak correction (**TwrPotent=2**), following Bak, Madsen, and Johansen (2001), adds a +drag-induced source term and offsets the streamwise coordinate by :math:`+0.1` (a fixed +empirical constant) to better represent the near wake. Writing +:math:`\overline{x}' = \overline{x} + 0.1` and +:math:`\overline{r}' = \sqrt{ \overline{x}'^2 + \overline{y}^2 }`, + +.. math:: + u_{TwrPotent} = \frac{-\overline{x}'^2 + \overline{y}^2}{\overline{r}'^4} + + \frac{C_d\,\overline{x}'}{2\pi\,\overline{r}'^2}, + \qquad + v_{TwrPotent} = \frac{-2\,\overline{x}'\,\overline{y}}{\overline{r}'^4} + + \frac{C_d\,\overline{y}}{2\pi\,\overline{r}'^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 :numref:`AD_twr_ends`. + .. _AD_twr_shadow: Tower shadow models -~~~~~~~~~~~~~~~~~~~ +---------------------- Powles tower shadow model (**TwrShadow=1**) is given by: .. math:: u_{TwrShadow} = - \frac{C_d}{ \sqrt{\overline{r}} } - \cos\left( \frac{\pi/2 \overline{y}}{\sqrt{\overline{r}}}\right)^2 + \cos\left( \frac{\pi}{2}\,\frac{\overline{y}}{\sqrt{\overline{r}}}\right)^2 where :math:`\overline{r} = \sqrt{ \overline{x}^2 + \overline{y}^2 }`. @@ -124,6 +173,269 @@ Eames tower shadow model (**TwrShadow=2**) is given by: where :math:`TI` is the turbulence intensity at the tower node. +To avoid excessive flow reversal behind the tower, the shadow +deficit fraction is limited to :math:`u_{TwrShadow} \ge -0.5`. + +The potential-flow and shadow contributions are superimposed and scaled by +:math:`W_\text{tower}` to give the velocity perturbation at the blade node in the local tower +frame: + +.. math:: + v_x = \left( u_{TwrPotent} + u_{TwrShadow} \right) W_\text{tower}, + \qquad + v_y = v_{TwrPotent}\, W_\text{tower} + +This perturbation :math:`(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. + +.. _AD_twr_ends: + +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 :math:`\overline{z} = 0`. When a blade +node lies axially beyond a tower end, the nearest point becomes the tower end itself, and +:math:`\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 :math:`c = \lVert \mathbf{r} \rVert - R`, +where :math:`\mathbf{r}` is the vector from the nearest tower point to the blade node and +:math:`R = \tfrac{1}{2}` ``TwrDiam`` is the local tower radius. The disturbance is suppressed both +very close to the tower (:math:`c \le 0.01\,` ``TwrDiam``) and far from it (:math:`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 :math:`|\overline{z}| < 1` the in-plane coordinates are divided by +:math:`\cos\!\left(\tfrac{\pi}{2}\overline{z}\right)`, + +.. math:: + \overline{x} \;\rightarrow\; \frac{\overline{x}}{\cos\!\left(\tfrac{\pi}{2}\overline{z}\right)}, + \qquad + \overline{y} \;\rightarrow\; \frac{\overline{y}}{\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 :math:`\overline{x}` and :math:`\overline{y}`. Because the +potential-flow disturbance velocity scales as +:math:`1/\overline{r}^{\,2} = 1/(\overline{x}^2 + \overline{y}^2)`, dividing :math:`\overline{x}` +and :math:`\overline{y}` in this way multiplies the potential-flow velocity by exactly +:math:`\cos^2\!\left(\tfrac{\pi}{2}\overline{z}\right)`: it decreases smoothly from its full value +at the end plane (:math:`\overline{z} = 0`) to zero at :math:`\overline{z} = 1` with zero slope +there (:math:`C^1`-continuous). The shadow deficit is tapered by the same coordinate scaling, but +its dependence on :math:`\overline{x}` and :math:`\overline{y}` is different, so the attenuation is +not a clean :math:`\cos^2` factor (for the Powles model the deficit scales roughly as +:math:`\sqrt{\cos}` with an additional shift in its lateral argument, and for the Eames model it +scales linearly in :math:`\cos` with the lateral profile unchanged). Both shadow models +nonetheless fade to zero as :math:`\overline{z} \rightarrow 1`. No tower influence is applied for +:math:`|\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. + +.. _AD_twr_drag: + +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 :math:`j` the relative wind is +:math:`\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 :math:`\mathbf{V}_\text{rel}` in the plane normal to the +tower axis produces drag; denote this transverse relative-wind vector :math:`\mathbf{V}_\perp` +and its magnitude :math:`W_\text{tower} = \lVert \mathbf{V}_\perp \rVert`. The drag force per +unit length is + +.. math:: + \mathbf{f}_\text{drag} = \tfrac{1}{2}\, \rho\, C_d\, D\, W_\text{tower}\, \mathbf{V}_\perp + +where :math:`\rho` is the air (or water, for MHK) density, :math:`C_d` = ``TwrCd`` is the local +tower drag coefficient, and :math:`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. + +.. _AD_gs_influence: + +Generalized support-structure influence models +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The generalized support structure (GS) extends the tower influence models of +:numref:`AD_twr_influence` 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 :math:`(\overline{x}, \overline{y}, +\overline{z})`, the same baseline and Bak potential-flow expressions (see +:numref:`AD_twr_potent`), the same Powles and Eames shadow expressions (see +:numref:`AD_twr_shadow`), and the same end handling --- clearance-based exclusion capsule and +axial taper (see :numref:`AD_twr_ends`). The member diameter and drag +coefficient play the roles that ``TwrDiam`` and :math:`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. + +.. _AD_gs_potent_combine: + +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 +:math:`\mathbf{v}_i` be the potential-flow contribution of member :math:`i` at the blade node +(expressed in the earth-fixed frame). The combined potential-flow disturbance is + +.. math:: + \mathbf{v}_\text{pot} = \frac{\sum_i w_i\, \mathbf{v}_i}{\sum_i w_i}, + \qquad + w_i = \lVert \mathbf{v}_i \rVert^{\,p} + +with blend exponent :math:`p = 6`. Because the weights are non-negative and sum to one, the +blended magnitude never exceeds :math:`\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* :math:`\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 :math:`p` controls the sharpness of the handoff: :math:`p \rightarrow \infty` + recovers a hard "nearest/strongest body only" selection, while a finite :math:`p` smooths the + transition. An even integer is used so that :math:`w_i = (\mathbf{v}_i \cdot + \mathbf{v}_i)^{p/2}` is a polynomial in the velocity components and therefore smooth + everywhere. The value :math:`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 :numref:`AD_twr_ends` 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. + +.. _AD_gs_shadow_combine: + +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 +:math:`\mathbf{s}_i` for member :math:`i`'s shadow contribution in the earth-fixed frame and +:math:`\mathbf{s}_\text{sum} = \sum_i \mathbf{s}_i`, + +.. math:: + \mathbf{v}_\text{shad} + = \left( \sqrt{ \sum_i \lVert \mathbf{s}_i \rVert^2 } \right) + \frac{\mathbf{s}_\text{sum}}{\lVert \mathbf{s}_\text{sum} \rVert} + +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. + +.. _AD_gs_drag: + +Generalized support-structure drag loads +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The generalized support structure carries the same cross-flow drag load as the tower +(:numref:`AD_twr_drag`), applied member by member when GS aerodynamics are enabled +(**GSAero=True**). For a member element with unit axial vector :math:`\hat{\mathbf{k}}`, the +transverse relative wind at a node is obtained by removing the along-member component of the +relative wind, + +.. math:: + \mathbf{V}_\perp = \mathbf{V}_\text{rel} - (\mathbf{V}_\text{rel}\cdot\hat{\mathbf{k}})\,\hat{\mathbf{k}}, + \qquad + \mathbf{V}_\text{rel} = \mathbf{V}_\text{inflow} - \mathbf{V}_\text{motion}, + +and the drag force per unit length takes the same form as for the tower, with the member +diameter :math:`D = 2R` and member drag coefficient :math:`C_d` in place of ``TwrDiam`` and +``TwrCd``: + +.. math:: + \mathbf{f}_\text{drag} = \tfrac{1}{2}\, \rho\, C_d\, D\, W_\text{GS}\, \mathbf{V}_\perp, + \qquad W_\text{GS} = \lVert \mathbf{V}_\perp \rVert + +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 :math:`\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 :numref:`AD_buoyancy` and :numref:`AD_addedmass_inertia`) --- 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 + (:numref:`AD_gs_influence`) 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. .. _AD_buoyancy: diff --git a/docs/source/user/api_change.rst b/docs/source/user/api_change.rst index 399fedf6ef..f877be8949 100644 --- a/docs/source/user/api_change.rst +++ b/docs/source/user/api_change.rst @@ -14,6 +14,8 @@ OpenFAST v5.0.x to OpenFAST v5.1.0 Under-relaxation is introduced for the tight-coupling iterative solver to improve numerical stability, requiring two new inputs in the main OpenFAST input file. +The generalized support-structure (GS) influence model was added to AeroDyn. This introduces three new switches (``GSPotent``, ``GSShadow``, ``GSAero``) after ``TwrAero`` in the AeroDyn primary input file, and two new sections (``General support structure joints`` and ``General support structure members``) after the ``Tower Influence and Aerodynamics`` section. The two new sections are required even when the GS model is disabled (set ``NumGSJoints`` and ``NumGSMembers`` to 0, keeping the table header lines). + ============================================= ======== ==================== ========================================================================================================================================================================================================================================================================================================== Added in OpenFAST `5.1.0` -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- @@ -23,6 +25,17 @@ OpenFAST 14 AutoRelax defa OpenFAST 15 RelaxFactor default RelaxFactor - Constant or initial (if AutoRelax) under-relaxation factor for the tight-coupling iterative solver (-) [>0 and <=1; default=0.7 if AutoRelax=false; default=0.3 if AutoRelax=true] HydroDyn 28 FKMod 0 FKMod - Nonlinear Froude-Krylov and hydrostatic load model {0: none, 1: nonlinear F-K and hydrostatics} (switch) [1 to NBody if NBodyMod>1; only used when PotMod=1; must provide GeoFile if FKMod=1] HydroDyn 29 GeoFile "unused" GeoFile - Full name(s) of geometry file (ASCII STL format) for nonlinear Froude-Krylov and hydrostatic load integration [1 to NBody; only used when PotMod=1 and when the corresponding FKMod=1] +AeroDyn 10 GSPotent 0 GSPotent - Type generalized support (GS) structure influence on wind based on potential flow around the GS members (switch) {0=none, 1=baseline potential flow, 2=potential flow with Bak correction} +AeroDyn 11 GSShadow 0 GSShadow - Calculate GS influence on wind based on downstream shadow (switch) {0=none, 1=Powles model, 2=Eames model} +AeroDyn 12 GSAero False GSAero - Calculate GS aerodynamic loads? (flag) +AeroDyn \* ====== General support structure joints =========================================================== [new section, added after the "Tower Influence and Aerodynamics" section] +AeroDyn \* NumGSJoints 0 NumGSJoints - Number of general support joints (-) +AeroDyn \* GSJointID GSJointxi GSJointyi GSJointzi +AeroDyn \* (-) (m) (m) (m) +AeroDyn \* ====== General support structure members ========================================================== [new section, added after the joints section] +AeroDyn \* NumGSMembers 0 NumGSMembers - Number of general support members (-) +AeroDyn \* GSMemberID GSMJointID1 GSMJointID2 GSMDia1 GSMDia2 GSMCd1 GSMCd2 GSMTI1 GSMTI2 GSMDiv +AeroDyn \* (-) (-) (-) (m) (m) (-) (-) (-) (-) (m) ============================================= ======== ==================== ========================================================================================================================================================================================================================================================================================================== OpenFAST v4.2.x to OpenFAST v5.0.0 diff --git a/modules/aerodyn/src/AeroDyn.f90 b/modules/aerodyn/src/AeroDyn.f90 index e6d6da99c7..9da727a7bd 100644 --- a/modules/aerodyn/src/AeroDyn.f90 +++ b/modules/aerodyn/src/AeroDyn.f90 @@ -63,6 +63,7 @@ module AeroDyn ! states(z) PUBLIC :: AD_VarsPackExtInput !< Routine pack extended inputs public :: AD_CalcWind_Rotor !< Routine to calculate rotor wind inputs + public :: AD_CalcWind_GS !< Routine to calculate general support structure wind inputs contains !---------------------------------------------------------------------------------------------------------------------------------- @@ -384,6 +385,9 @@ subroutine AD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut InputFileData%TwrPotent = TwrPotent_none InputFileData%TwrShadow = TwrShadow_none InputFileData%TwrAero = TwrAero_none + InputFileData%GS%GSPotent = GSPotent_none + InputFileData%GS%GSShadow = GSShadow_none + InputFileData%GS%GSAero = GSAero_none !InputFileData%CavitCheck = .false. !InputFileData%TFinAero = .false. ! not sure if this needs to be set or not InputFileData%DBEMT_Mod = DBEMT_none @@ -442,7 +446,22 @@ subroutine AD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut ! Set pointer to FlowField data if (associated(InitInp%FlowField)) p%FlowField => InitInp%FlowField - + ! Initialize general support structure parameters + call Init_GSParam(InputFileData%GS, p%GS, InputFileData%AirDens, InitInp%MHK, InitInp%WtrDpth, errStat2, errMsg2 ) + if (Failed()) return; + if (nRotors >= 1_IntKi) then + p%rotors( 1)%hasGSMod = p%GS%hasGSMod + p%rotors( 1)%GSPotent = p%GS%GSPotent + p%rotors( 1)%GSShadow = p%GS%GSShadow + p%rotors( 1)%GSAero = p%GS%GSAero ! Only rotor 1 handles GS drag force + end if + do iR = 2, nRotors + p%rotors(iR)%hasGSMod = p%GS%hasGSMod + p%rotors(iR)%GSPotent = p%GS%GSPotent + p%rotors(iR)%GSShadow = p%GS%GSShadow + p%rotors(iR)%GSAero = GSAero_none ! Only rotor 1 handles GS drag force + enddo + !............................................................................................ ! Define and initialize inputs here !............................................................................................ @@ -546,6 +565,13 @@ subroutine AD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut if (Failed()) return enddo + !............................................................................................ + ! Initialize m%Inflow%GSInflow for tracking wind inflow + !............................................................................................ + call AllocAry( m%Inflow(1)%GSInflow%InflowVel, 3_IntKi, p%GS%NNodes, 'GSInflow%InflowVel', ErrStat2, ErrMsg2 ) + if (Failed()) return + m%Inflow(1)%GSInflow%InflowVel = 0.0_ReKi + ! Duplicte Inflow(1) (must be done after Init_OLAF) call AD_CopyInflowType(m%Inflow(1), m%Inflow(2), MESH_NEWCOPY, ErrStat2, ErrMsg2) if (Failed()) return @@ -611,8 +637,10 @@ subroutine AD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut call AD_PrintSum( InputFileData, p%rotors(iR), p, u, y, NumBlades(iR), InputFileData%rotors(iR)%BladeProps(:), ErrStat2, ErrMsg2 ) if (Failed()) return; enddo + if ( p%GS%hasGSMod ) then + call AD_PrintSum_GS( p%GS, p, u%rotors(1), ErrStat2, ErrMsg2 ); if (Failed()) return + end if end if - !............................................................................................ ! If you want to choose your own rate instead of using what the glue code suggests, tell the glue code the rate at which ! this module must be called here: @@ -760,7 +788,10 @@ subroutine Init_MiscVars(m, p, p_AD, u, y, errStat, errMsg) call AllocAry( m%TwrClrnc, p%NumBlNds, p%NumBlades, 'm%TwrClrnc', ErrStat2, ErrMsg2 ) call SetErrStat( errStat2, errMsg2, errStat, errMsg, RoutineName ) end if - + if (p%GSPotent /= GSPotent_none .or. p%GSShadow /= GSShadow_none) then + call AllocAry( m%GSClrnc, p%NumBlNds, p%NumBlades, 'm%GSClrnc', ErrStat2, ErrMsg2 ) + call SetErrStat( errStat2, errMsg2, errStat, errMsg, RoutineName ) + end if call AllocAry( m%Cant, p%NumBlNds, p%NumBlades, 'm%Cant', ErrStat2, ErrMsg2 ) call SetErrStat( errStat2, errMsg2, errStat, errMsg, RoutineName ) call AllocAry( m%Toe, p%NumBlNds, p%NumBlades, 'm%Toe', ErrStat2, ErrMsg2 ) @@ -1062,6 +1093,18 @@ subroutine Init_y(y, u, p, errStat, errMsg) errStat = ErrID_None errMsg = "" + if (p%GSAero/=GSAero_none) then + call MeshCopy ( SrcMesh = u%GSMotion & + , DestMesh = y%GSLoad & + , CtrlCode = MESH_SIBLING & + , IOS = COMPONENT_OUTPUT & + , force = .TRUE. & + , moment = .TRUE. & + , ErrStat = ErrStat2 & + , ErrMess = ErrMsg2 ) + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + if (ErrStat >= AbortErrLev) RETURN + end if if (p%NumTwrNds > 0 .and. (p%TwrAero /= TwrAero_None .or. p%MHK /= MHK_None)) then @@ -1172,12 +1215,17 @@ subroutine Init_u( u, p, p_AD, InputFileData, MHK, WtrDpth, InitInp, errStat, er ! Local variables real(reKi) :: position(3) ! node reference position real(reKi) :: positionL(3) ! node local position + real(ReKi) :: p1(3) ! node reference position + real(ReKi) :: p2(3) ! node reference position real(R8Ki) :: theta(3) ! Euler angles real(R8Ki) :: orientation(3,3) ! node reference orientation real(R8Ki) :: orientationL(3,3) ! node local orientation + real(ReKi) :: s ! normalized distance along a GS member integer(intKi) :: j ! counter for nodes integer(intKi) :: k ! counter for blades + integer(intKi) :: nnodes ! number of nodes of a GS member + integer(intKi) :: nodeCntr ! node counter integer(intKi) :: ErrStat2 ! temporary Error status character(ErrMsgLen) :: ErrMsg2 ! temporary Error message @@ -1192,11 +1240,72 @@ subroutine Init_u( u, p, p_AD, InputFileData, MHK, WtrDpth, InitInp, errStat, er call SetErrStat( errStat2, errMsg2, errStat, errMsg, RoutineName ) if (errStat >= AbortErrLev) return - - + u%UserProp = 0.0_ReKi - - ! Meshes for motion inputs (ElastoDyn and/or BeamDyn) + + ! Meshes for motion inputs (ElastoDyn, SubDyn, and/or BeamDyn) + !................ + ! general support structure + !................ + if (p%hasGSMod) then + + call MeshCreate ( BlankMesh = u%GSMotion & + ,IOS = COMPONENT_INPUT & + ,Nnodes = p_AD%GS%NNodes & + ,ErrStat = ErrStat2 & + ,ErrMess = ErrMsg2 & + ,TranslationDisp = .true. & ! Orientation is not needed: the GS models only read position, displacement and velocity + ,TranslationVel = .true. & + ) + if (failed()) return + + ! Add GS joints to the mesh first + nodeCntr = 0_IntKi + do j=1,p_AD%GS%NJoints + position = p_AD%GS%Joints(j)%position ! z offset for MHK_FixedBottom was already done in Init_GSParam + nodeCntr = nodeCntr + 1_IntKi + call MeshPositionNode(u%GSMotion, nodeCntr, position, errStat2, errMsg2) ! orientation is identity by default + if (failed()) return + end do !j + + ! Add member interior nodes to the mesh next + do j=1,p_AD%GS%NMembers + nnodes = p_AD%GS%Members(j)%NElements + 1_IntKi + p1 = p_AD%GS%Joints(p_AD%GS%Members(j)%NodeIndx( 1))%position + p2 = p_AD%GS%Joints(p_AD%GS%Members(j)%NodeIndx(nnodes))%position + do k=2_IntKi,nnodes-1_IntKi + s = real(k-1_IntKi,ReKi)/real(nnodes-1_IntKi,ReKi) + position = p1 * (1.0_ReKi-s) + p2 * s + nodeCntr = nodeCntr + 1_IntKi + call MeshPositionNode(u%GSMotion, nodeCntr, position, errStat2, errMsg2) ! orientation is identity by default + if (failed()) return + end do + ! ! create line2 elements + ! do k=1,nnodes-1_IntKi + ! call MeshConstructElement( u%GSMotion, ELEMENT_LINE2, errStat2, errMsg2, & + ! p1=p_AD%GS%Members(j)%NodeIndx( k), & + ! p2=p_AD%GS%Members(j)%NodeIndx(k+1) ) + ! call SetErrStat( errStat2, errMsg2, errStat, errMsg, RoutineName ) + ! end do + end do + + do j=1,nodeCntr + call MeshConstructElement( u%GSMotion & + ,ELEMENT_POINT & + ,errStat2 & + ,errMsg2 & + ,j & + ) + if (failed()) return + end do + + call MeshCommit(u%GSMotion, errStat2, errMsg2 ); if (failed()) return + + u%GSMotion%TranslationDisp = 0.0_R8Ki + u%GSMotion%TranslationVel = 0.0_ReKi + + end if + !................ ! tower !................ @@ -1482,6 +1591,9 @@ subroutine SetParameters( InitInp, InputFileData, RotData, p, p_AD, ErrStat, Err p%TwrPotent = InputFileData%TwrPotent p%TwrShadow = InputFileData%TwrShadow p%TwrAero = InputFileData%TwrAero + p%GSPotent = InputFileData%GS%GSPotent + p%GSShadow = InputFileData%GS%GSShadow + p%GSAero = InputFileData%GS%GSAero p%CavitCheck = InputFileData%CavitCheck p%NacelleDrag = InputFileData%NacelleDrag @@ -1885,7 +1997,7 @@ subroutine AD_UpdateStates( t, n, u, utimes, p, x, xd, z, OtherState, m, errStat if (Failed()) return do iR = 1,size(p%rotors) - call SetInputs(t, p%rotors(iR), p, uInterp%rotors(iR), InflowInterp%RotInflow(iR), m%rotors(iR), i, errStat2, errMsg2) + call SetInputs(t, p%rotors(iR), p, uInterp%rotors(iR), InflowInterp%RotInflow(iR), InflowInterp%GSInflow, m%rotors(iR), i, errStat2, errMsg2) if (Failed()) return enddo end do @@ -1973,6 +2085,11 @@ subroutine AD_CalcWind(t, u, FLowField, p, m, o, Inflow, ErrStat, ErrMsg) if(Failed()) return enddo + ! The GS is a single shared structure: every rotor's GSMotion mesh is driven from the same + ! substructure/platform source, so the GS inflow is sampled once from rotor 1 and reused for all rotors. + call AD_CalcWind_GS(t, u%rotors(1), FlowField, p%GS, p, m, Inflow%GSInflow, StartNode, ErrStat2, ErrMsg2) + if(Failed()) return + ! OLAF points if (allocated(o%WakeLocationPoints) .and. allocated(Inflow%InflowWakeVel)) then ! If rotor is MHK, add water depth to z coordinate @@ -1983,11 +2100,14 @@ subroutine AD_CalcWind(t, u, FLowField, p, m, o, Inflow, ErrStat, ErrMsg) end if if (p%FVW%MHK /= MHK_None .and. p%CompSeaSt) then ! MHK turbines with waves + ! BoxExceedAllow=.true. matches the InflowWind branch below: far-wake OLAF particles + ! that drift out of the current/wave grid extrapolate instead of aborting the run. call WaveField_GetWaveVelAcc_AD(p%WaveField, m%WaveField_m, & StartNode, t, & o%WakeLocationPoints, & Inflow%InflowWakeVel, & - NoAcc, ErrStat2, ErrMsg2) + NoAcc, ErrStat2, ErrMsg2, & + BoxExceedAllow=.true.) if(Failed()) return else call IfW_FlowField_GetVelAcc(FlowField, StartNode, t, & @@ -2143,11 +2263,63 @@ subroutine AD_CalcWind_Rotor(t, u, FlowField, p, p_AD, m, RotInflow, StartNode, contains logical function Failed() - call SetErrStat(errStat2, errMsg2, errStat, errMsg, 'AD_CalcWindRotor') + call SetErrStat(errStat2, errMsg2, errStat, errMsg, 'AD_CalcWind_Rotor') Failed = errStat >= AbortErrLev end function Failed end subroutine +subroutine AD_CalcWind_GS(t, u, FlowField, p, p_AD, m, GSInflow, StartNode, ErrStat, ErrMsg) + real(DbKi), intent(in ) :: t !< Current simulation time in seconds + type(RotInputType), intent(in ) :: u !< Inputs at Time t + type(FlowFieldType),pointer, intent(in ) :: FlowField !< Pointer to IfW flowfield + type(GSParameterType), intent(in ) :: p !< Parameters + type(AD_ParameterType), intent(in ) :: p_AD !< AD parameters + type(AD_MiscVarType), intent(inout) :: m !< Misc/optimization variables + type(ElemInflowType),target, intent(inout) :: GSInflow !< calculated inflow + integer(IntKi), intent(inout) :: StartNode !< starting node for rotor wind + integer(IntKi), intent( out) :: ErrStat !< Error status of the operation + character(*), intent( out) :: ErrMsg !< Error message if ErrStat /= ErrID_None + + integer(intKi) :: ErrStat2 + character(ErrMsgLen) :: ErrMsg2 + real(ReKi) :: PosOffset(3) + real(ReKi), allocatable :: NoAcc(:,:) + + ErrStat = ErrID_None + ErrMsg = "" + + if (.not.p%hasGSMod) return + + ! If rotor is MHK, add water depth to z coordinate + if (p%MHK /= MHK_None) then + PosOffset = [0.0_ReKi, 0.0_ReKi, p%WtrDpth] + else + PosOffset = 0.0_ReKi + end if + if (p%MHK /= MHK_None .and. p_AD%CompSeaSt) then ! MHK turbines with waves + call WaveField_GetWaveVelAcc_AD(p_AD%WaveField, m%WaveField_m, & + StartNode, t, & + real(u%GSMotion%Position + u%GSMotion%TranslationDisp, ReKi), & + GSInflow%InflowVel, & + NoAcc, ErrStat2, ErrMsg2, & + BoxExceedAllow=.true.) + if(Failed()) return + else + call IfW_FlowField_GetVelAcc(FlowField, StartNode, t, & + real(u%GSMotion%Position + u%GSMotion%TranslationDisp, ReKi), & + GSInflow%InflowVel, & + NoAcc, ErrStat2, ErrMsg2, & + BoxExceedAllow=.true., PosOffset=PosOffset) + if(Failed()) return + end if + StartNode = StartNode + p%NNodes + +contains + logical function Failed() + call SetErrStat(errStat2, errMsg2, errStat, errMsg, 'AD_CalcWind_GS') + Failed = errStat >= AbortErrLev + end function Failed +end subroutine !---------------------------------------------------------------------------------------------------------------------------------- !> Routine for computing outputs, used in both loose and tight coupling. @@ -2195,7 +2367,7 @@ subroutine AD_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, ! SetInputs, Calc BEM Outputs and Twr Outputs do iR=1,size(p%rotors) - call RotCalcOutput(t, u%rotors(iR), m%Inflow(1)%RotInflow(iR), p%rotors(iR), p, x%rotors(iR), & + call RotCalcOutput(t, u%rotors(iR), m%Inflow(1)%RotInflow(iR), m%Inflow(1)%GSInflow, p%rotors(iR), p, x%rotors(iR), & xd%rotors(iR), z%rotors(iR), OtherState%rotors(iR), & y%rotors(iR), m%rotors(iR), m, iR, ErrStat2, ErrMsg2, .false.) if(Failed()) return @@ -2234,7 +2406,7 @@ logical function Failed() end function Failed end subroutine AD_CalcOutput !---------------------------------------------------------------------------------------------------------------------------------- -subroutine RotCalcOutput( t, u, RotInflow, p, p_AD, x, xd, z, OtherState, y, m, m_AD, iRot, ErrStat, ErrMsg, NeedWriteOutput) +subroutine RotCalcOutput( t, u, RotInflow, GSInflow, p, p_AD, x, xd, z, OtherState, y, m, m_AD, iRot, ErrStat, ErrMsg, NeedWriteOutput) ! NOTE: no matter how many channels are selected for output, all of the outputs are calculated ! All of the calculated output channels are placed into the m%AllOuts(:), while the channels selected for outputs are ! placed in the y%WriteOutput(:) array. @@ -2243,6 +2415,7 @@ subroutine RotCalcOutput( t, u, RotInflow, p, p_AD, x, xd, z, OtherState, y, m, REAL(DbKi), INTENT(IN ) :: t !< Current simulation time in seconds TYPE(RotInputType), INTENT(IN ) :: u !< Inputs at Time t TYPE(RotInflowType), INTENT(IN ) :: RotInflow !< Rotor Inflow at Time t + TYPE(ElemInflowType), INTENT(IN ) :: GSInflow !< General support structure Inflow at Time t TYPE(RotParameterType), INTENT(IN ) :: p !< Parameters TYPE(AD_ParameterType), INTENT(IN ) :: p_AD !< Parameters TYPE(RotContinuousStateType), INTENT(IN ) :: x !< Continuous states at t @@ -2276,7 +2449,7 @@ subroutine RotCalcOutput( t, u, RotInflow, p, p_AD, x, xd, z, OtherState, y, m, CalcWriteOutput = .true. ! by default, calculate WriteOutput unless told that we do not need it end if - call SetInputs(t, p, p_AD, u, RotInflow, m, indx, errStat2, errMsg2) + call SetInputs(t, p, p_AD, u, RotInflow, GSInflow, m, indx, errStat2, errMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (p_AD%Wake_Mod /= WakeMod_FVW) then @@ -2299,12 +2472,16 @@ subroutine RotCalcOutput( t, u, RotInflow, p, p_AD, x, xd, z, OtherState, y, m, end if endif - if ( p%TwrAero /= TwrAero_none ) then call ADTwr_CalcOutput(p, u, RotInflow, m, y, ErrStat2, ErrMsg2 ) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) endif + if (p%GSAero/=GSAero_none) then + call ADGS_CalcOutput(p_AD%GS, u, GSInflow, m, y, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + endif + ! initialize nacelle mesh loads y%NacelleLoad%Force = 0.0_ReKi y%NacelleLoad%Moment = 0.0_ReKi @@ -3050,19 +3227,20 @@ subroutine AD_CalcConstrStateResidual( Time, u, p, x, xd, z, OtherState, m, z_re do iR=1, size(p%rotors) - call RotCalcConstrStateResidual( Time, u%rotors(iR), m%Inflow(1)%RotInflow(iR), p%rotors(iR), p, x%rotors(iR), xd%rotors(iR), z%rotors(iR), OtherState%rotors(iR), m%rotors(iR), z_residual%rotors(iR), ErrStat2, ErrMsg2 ) + call RotCalcConstrStateResidual( Time, u%rotors(iR), m%Inflow(1)%RotInflow(iR), m%Inflow(1)%GSInflow, p%rotors(iR), p, x%rotors(iR), xd%rotors(iR), z%rotors(iR), OtherState%rotors(iR), m%rotors(iR), z_residual%rotors(iR), ErrStat2, ErrMsg2 ) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) enddo end subroutine AD_CalcConstrStateResidual !---------------------------------------------------------------------------------------------------------------------------------- !> Tight coupling routine for solving for the residual of the constraint state equations -subroutine RotCalcConstrStateResidual( Time, u, RotInflow, p, p_AD, x, xd, z, OtherState, m, z_residual, ErrStat, ErrMsg ) +subroutine RotCalcConstrStateResidual( Time, u, RotInflow, GSInflow, p, p_AD, x, xd, z, OtherState, m, z_residual, ErrStat, ErrMsg ) !.................................................................................................................................. REAL(DbKi), INTENT(IN ) :: Time !< Current simulation time in seconds TYPE(RotInputType), INTENT(IN ) :: u !< Inputs at Time TYPE(RotInflowType), INTENT(IN ) :: RotInflow !< rotor inflow at Time + TYPE(ElemInflowType), INTENT(IN ) :: GSInflow !< General support structure inflow at Time TYPE(RotParameterType), INTENT(IN ) :: p !< Parameters TYPE(AD_ParameterType), INTENT(IN ) :: p_AD !< Parameters TYPE(RotContinuousStateType), INTENT(IN ) :: x !< Continuous states at Time @@ -3090,7 +3268,7 @@ subroutine RotCalcConstrStateResidual( Time, u, RotInflow, p, p_AD, x, xd, z, Ot end if - call SetInputs(Time, p, p_AD, u, RotInflow, m, indx, errStat2, errMsg2) + call SetInputs(Time, p, p_AD, u, RotInflow, GSInflow, m, indx, errStat2, errMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) @@ -3101,13 +3279,14 @@ subroutine RotCalcConstrStateResidual( Time, u, RotInflow, p, p_AD, x, xd, z, Ot end subroutine RotCalcConstrStateResidual !---------------------------------------------------------------------------------------------------------------------------------- -subroutine RotCalcContStateDeriv( t, u, RotInflow, p, p_AD, x, xd, z, OtherState, m, dxdt, ErrStat, ErrMsg ) +subroutine RotCalcContStateDeriv( t, u, RotInflow, GSInflow, p, p_AD, x, xd, z, OtherState, m, dxdt, ErrStat, ErrMsg ) ! Tight coupling routine for computing derivatives of continuous states !.................................................................................................................................. REAL(DbKi), INTENT(IN ) :: t ! Current simulation time in seconds TYPE(RotInputType), INTENT(IN ) :: u ! Inputs at t TYPE(RotInflowType), INTENT(IN ) :: RotInflow !< Rotor inflow Inputs at Time + TYPE(ElemInflowType), INTENT(IN ) :: GSInflow !< General support structure inflow at Time TYPE(RotParameterType), INTENT(IN ) :: p ! Parameters TYPE(AD_ParameterType), INTENT(IN ) :: p_AD ! Parameters TYPE(RotContinuousStateType), INTENT(IN ) :: x ! Continuous states at t @@ -3131,7 +3310,7 @@ subroutine RotCalcContStateDeriv( t, u, RotInflow, p, p_AD, x, xd, z, OtherState ErrStat = ErrID_None ErrMsg = "" - call SetInputs(t, p, p_AD, u, RotInflow, m, InputIndex, ErrStat2, ErrMsg2) + call SetInputs(t, p, p_AD, u, RotInflow, GSInflow, m, InputIndex, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) call BEMT_CalcContStateDeriv( t, m%BEMT_u(InputIndex), p%BEMT, x%BEMT, xd%BEMT, z%BEMT, OtherState%BEMT, m%BEMT, dxdt%BEMT, p_AD%AFI, ErrStat2, ErrMsg2 ) @@ -3141,17 +3320,17 @@ END SUBROUTINE RotCalcContStateDeriv !---------------------------------------------------------------------------------------------------------------------------------- !> This subroutine converts the AeroDyn inputs into values that can be used for its submodules. It calculates the disturbed inflow !! on the blade if tower shadow or tower influence are enabled, then uses these values to set m%BEMT_u(indx). -subroutine SetInputs(t, p, p_AD, u, RotInflow, m, indx, errStat, errMsg) +subroutine SetInputs(t, p, p_AD, u, RotInflow, GSInflow, m, indx, errStat, errMsg) real(DbKi), intent(in ) :: t !< Current simulation time in seconds type(RotParameterType), intent(in ) :: p !< AD parameters type(AD_ParameterType), intent(in ) :: p_AD !< AD parameters type(RotInputType), intent(in ) :: u !< AD Inputs at Time type(RotInflowType), intent(in ) :: RotInflow !< Rotor inflow Inputs at Time + type(ElemInflowType), intent(in ) :: GSInflow !< Inflow on the general support structure at Time t type(RotMiscVarType), intent(inout) :: m !< Misc/optimization variables integer, intent(in ) :: indx !< index into m%BEMT_u(indx) array; 1=t and 2=t+dt (but not checked here) integer(IntKi), intent( out) :: ErrStat !< Error status of the operation character(*), intent( out) :: ErrMsg !< Error message if ErrStat /= ErrID_None - ! local variables integer(intKi) :: ErrStat2 character(ErrMsgLen) :: ErrMsg2 @@ -3159,8 +3338,8 @@ subroutine SetInputs(t, p, p_AD, u, RotInflow, m, indx, errStat, errMsg) ErrStat = ErrID_None ErrMsg = "" - ! Disturbed inflow on blade (if tower shadow present) - call SetDisturbedInflow(p, p_AD, u, RotInflow, m, errStat2, errMsg2); call SetErrStat(errStat2, errMsg2, errStat, errMsg, RoutineName) + ! Disturbed inflow on blade (if tower shadow present, and if general support structure influence present) + call SetDisturbedInflow(p, p_AD, u, RotInflow, GSInflow, m, errStat2, errMsg2); call SetErrStat(errStat2, errMsg2, errStat, errMsg, RoutineName) if (p_AD%Wake_Mod /= WakeMod_FVW) then @@ -3176,11 +3355,12 @@ end subroutine SetInputs !---------------------------------------------------------------------------------------------------------------------------------- !> Disturbed inflow on the blade if tower shadow or tower influence are enabled -subroutine SetDisturbedInflow(p, p_AD, u, RotInflow, m, errStat, errMsg) +subroutine SetDisturbedInflow(p, p_AD, u, RotInflow, GSInflow, m, errStat, errMsg) type(RotParameterType), intent(in ) :: p !< AD parameters type(AD_ParameterType), intent(in ) :: p_AD !< AD parameters type(RotInputType), intent(in ) :: u !< AD Inputs at Time type(RotInflowType), intent(in ) :: RotInflow !< Rotor inflow at Time + type(ElemInflowType), intent(in ) :: GSInflow !< Inflow on the general support structure at Time t type(RotMiscVarType), intent(inout) :: m !< Misc/optimization variables integer(IntKi), intent( out) :: errStat !< Error status of the operation character(*), intent( out) :: errMsg !< Error message if ErrStat /= ErrID_None @@ -3201,6 +3381,12 @@ subroutine SetDisturbedInflow(p, p_AD, u, RotInflow, m, errStat, errMsg) end do end if + ! Generalized support structure (GS) influence, added on top of the tower (or undisturbed) inflow above + if (p%hasGSMod) then + call GSInfl( p, p_AD%GS, u, GSInflow, m, errStat2, errMsg2 ) + call SetErrStat(errStat2, errMsg2, errStat, errMsg, RoutineName) + end if + if (p_AD%Skew_Mod == Skew_Mod_Orthogonal) then x_hat_disk = u%HubMotion%Orientation(1,:,1) @@ -3333,6 +3519,9 @@ subroutine SetSectAvgInflow(t, p, p_AD, u, RotInflow, m, errStat, errMsg) ! TODO use a "scalar" function or change the interface of TwrInfl. Waiting for Wind Inputs of AD to be removed from AD call TwrInflArray( p, u, RotInflow, m, reshape(r_A, (/3,1/)), m%SectAvgInflow(:, j:j, k), errStat2, errMsg2); if(Failed()) return endif + ! REMINDER (future work): generalized support structure (GS) influence is NOT applied to the sector-averaged + ! inflow. GSInfl has no array/scalar form (no GSInflArray analogous to TwrInflArray), so with SectAvg enabled + ! the averaged inflow omits GS while the per-node inflow (via SetDisturbedInflow->GSInfl) includes it. enddo enddo @@ -4031,6 +4220,9 @@ subroutine SetInputsForFVW(p, u, tIndx, m, errStat, errMsg) m%FVW_u(tIndx)%V_wind = m%Inflow(tIndx)%InflowWakeVel ! Applying tower shadow to V_wind based on r_wind positions ! NOTE: m%DisturbedInflow also contains tower shadow and we need it for CalcOutput + ! REMINDER (future work): only tower influence is applied to the OLAF/FVW wake-convection wind here; + ! generalized support structure (GS) influence is not (no GSInflArray equivalent to TwrInflArray). + ! GS influence IS applied to the lifting-line disturbed inflow used for UA below (SetDisturbedInflow). if (p%FVW%TwrShadowOnWake) then do iR =1, size(p%rotors) if (p%rotors(iR)%TwrPotent /= TwrPotent_none .or. p%rotors(iR)%TwrShadow /= TwrShadow_none) then @@ -4043,7 +4235,8 @@ subroutine SetInputsForFVW(p, u, tIndx, m, errStat, errMsg) endif do iR =1, size(p%rotors) ! Disturbed inflow for UA on Lifting line Mesh Points - call SetDisturbedInflow(p%rotors(iR), p, u%rotors(iR), m%Inflow(tIndx)%RotInflow(iR), m%rotors(iR), errStat2, errMsg2) + ! GS influence is included here (regardless of Wake_Mod) so unsteady aero sees the GS-disturbed inflow + call SetDisturbedInflow(p%rotors(iR), p, u%rotors(iR), m%Inflow(tIndx)%RotInflow(iR), m%Inflow(tIndx)%GSInflow, m%rotors(iR), errStat2, errMsg2) call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) do k=1,p%rotors(iR)%NumBlades iW=p%FVW%Bld2Wings(iR,k) @@ -4401,9 +4594,20 @@ SUBROUTINE ValidateInputData( InitInp, InputFileData, NumBl, calcCrvAngle, ErrSt if (InputFileData%TwrAero /= TwrAero_none .and. InputFileData%TwrAero /= TwrAero_noVIV) then call SetErrStat ( ErrID_Fatal, 'TwrAero must be 0 (none) or 1 (Tower aero on).', ErrStat, ErrMsg, RoutineName ) end if + if (Failed()) return + + if (InputFileData%GS%GSAero /= GSAero_none .and. InputFileData%GS%GSAero /= GSAero_noVIV) then + call SetErrStat ( ErrID_Fatal, 'GSAero must be 0 (none) or 1 (multi-member generalized tower aero/hydro on).', ErrStat, ErrMsg, RoutineName ) + end if + if (InputFileData%GS%GSPotent /= GSPotent_none .and. InputFileData%GS%GSPotent /= GSPotent_baseline .and. InputFileData%GS%GSPotent /= GSPotent_Bak) then + call SetErrStat ( ErrID_Fatal, 'GSPotent must be 0 (none), 1 (baseline potential flow), or 2 (potential flow with Bak correction).', ErrStat, ErrMsg, RoutineName ) + end if + if (InputFileData%GS%GSShadow /= GSShadow_none .and. InputFileData%GS%GSShadow /= GSShadow_Powles .and. InputFileData%GS%GSShadow /= GSShadow_Eames) then + call SetErrStat ( ErrID_Fatal, 'GSShadow must be 0 (none), 1 (Powles model), or 2 (Eames model).', ErrStat, ErrMsg, RoutineName ) + end if if (Failed()) return - if (InitInp%MHK == MHK_None .and. InputFileData%CavitCheck) call SetErrStat ( ErrID_Fatal, 'A cavitation check can only be performed for an MHK turbine.', ErrStat, ErrMsg, RoutineName ) + if (InitInp%MHK == MHK_None .and. InputFileData%CavitCheck) call SetErrStat ( ErrID_Fatal, 'A cavitation check can only be performed for an MHK turbine.', ErrStat, ErrMsg, RoutineName ) if (InitInp%MHK /= MHK_None .and. InputFileData%CompAA ) call SetErrStat ( ErrID_Fatal, 'The aeroacoustics module cannot be used with an MHK turbine.', ErrStat, ErrMsg, RoutineName ) do iR = 1,size(NumBl) if (InitInp%MHK /= MHK_None .and. InputFileData%rotors(iR)%TFinAero) call SetErrStat ( ErrID_Fatal, 'A tail fin cannot be modeled for an MHK turbine.', ErrStat, ErrMsg, RoutineName ) @@ -4720,6 +4924,91 @@ SUBROUTINE ValidateInputData( InitInp, InputFileData, NumBl, calcCrvAngle, ErrSt end do end if + !.................. + ! check for generalized support structure + !.................. + ! Joints + if (InputFileData%GS%NJoints<0_IntKi) then + call SetErrStat( ErrID_Fatal, 'NumGSJoints cannot be negative.', ErrStat, ErrMsg, RoutineName ); if(Failed()) return + end if + if (InputFileData%GS%NJoints==1_IntKi) then + call SetErrStat( ErrID_Fatal, 'Generalized support structure requires at least two joints (NumGSJoints>=2).', ErrStat, ErrMsg, RoutineName ); if(Failed()) return + end if + do j = 1,InputFileData%GS%NJoints + if (InputFileData%GS%InpJoints(j)%JointID<0_IntKi) then + call SetErrStat( ErrID_Fatal, 'GSJointID cannot be negative; check row '//trim(num2lstr(j)), ErrStat, ErrMsg, RoutineName ) + if(Failed()) return + end if + do k = 1,InputFileData%GS%NJoints + if (j==k) cycle + if (InputFileData%GS%InpJoints(j)%JointID == InputFileData%GS%InpJoints(k)%JointID) then + call SetErrStat( ErrID_Fatal, 'GSJointID must be unique for each joint; check row '//trim(num2lstr(j))//' and row '//trim(num2lstr(k)), ErrStat, ErrMsg, RoutineName ) + if(Failed()) return + end if + end do + end do + ! Members + if (InputFileData%GS%NMembers<0_IntKi) then + call SetErrStat( ErrID_Fatal, 'NumGSMembers cannot be negative.', ErrStat, ErrMsg, RoutineName ); if(Failed()) return + end if + if (InputFileData%GS%NMembers==0_IntKi .and. (InputFileData%GS%GSPotent/=GSPotent_none .or. & + InputFileData%GS%GSShadow/=GSShadow_none .or. InputFileData%GS%GSAero/=GSAero_none)) then + call SetErrStat( ErrID_Warn, 'A generalized support (GS) model is enabled (GSPotent, GSShadow, or GSAero), but no GS members are defined (NumGSMembers=0); the GS model will be disabled.', ErrStat, ErrMsg, RoutineName ) + end if + do j = 1,InputFileData%GS%NMembers + if (InputFileData%GS%InpMembers(j)%MemberID<0_IntKi) then + call SetErrStat( ErrID_Fatal, 'GSMemberID cannot be negative; check row '//trim(num2lstr(j)), ErrStat, ErrMsg, RoutineName ); if(Failed()) return + end if + do k = 1,InputFileData%GS%NMembers + if (j==k) cycle + if (InputFileData%GS%InpMembers(j)%MemberID == InputFileData%GS%InpMembers(k)%MemberID) then + call SetErrStat( ErrID_Fatal, 'GSMemberID must be unique for each member; check row '//trim(num2lstr(j))//' and row '//trim(num2lstr(k)), ErrStat, ErrMsg, RoutineName ); if(Failed()) return + end if + end do + if (InputFileData%GS%InpMembers(j)%MJointID1<0_IntKi) then + call SetErrStat( ErrID_Fatal, 'GSMJointID1 cannot be negative; check row '//trim(num2lstr(j)), ErrStat, ErrMsg, RoutineName ); if(Failed()) return + end if + if (InputFileData%GS%InpMembers(j)%MJointID2<0_IntKi) then + call SetErrStat( ErrID_Fatal, 'GSMJointID2 cannot be negative; check row '//trim(num2lstr(j)), ErrStat, ErrMsg, RoutineName ); if(Failed()) return + end if + if (InputFileData%GS%InpMembers(j)%MJointID1 == InputFileData%GS%InpMembers(j)%MJointID2) then + call SetErrStat( ErrID_Fatal, 'GSMJointID1 and GSMJointID2 must be distinct for each member; check row '//trim(num2lstr(j)), ErrStat, ErrMsg, RoutineName ); if(Failed()) return + end if + if (InputFileData%GS%InpMembers(j)%MDiam1<0.0_ReKi) then + call SetErrStat( ErrID_Fatal, 'GSMDia1 cannot be negative; check row '//trim(num2lstr(j)), ErrStat, ErrMsg, RoutineName ); if(Failed()) return + elseif ((InputFileData%GS%GSPotent/=GSPotent_none .or. InputFileData%GS%GSShadow/=GSShadow_none) .and. EqualRealNos(InputFileData%GS%InpMembers(j)%MDiam1, 0.0_ReKi)) then + call SetErrStat( ErrID_Fatal, 'GSMDia1 must be greater than zero when GSPotent or GSShadow is enabled; check row '//trim(num2lstr(j)), ErrStat, ErrMsg, RoutineName ); if(Failed()) return + end if + if (InputFileData%GS%InpMembers(j)%MDiam2<0.0_ReKi) then + call SetErrStat( ErrID_Fatal, 'GSMDia2 cannot be negative; check row '//trim(num2lstr(j)), ErrStat, ErrMsg, RoutineName ); if(Failed()) return + elseif ((InputFileData%GS%GSPotent/=GSPotent_none .or. InputFileData%GS%GSShadow/=GSShadow_none) .and. EqualRealNos(InputFileData%GS%InpMembers(j)%MDiam2, 0.0_ReKi)) then + call SetErrStat( ErrID_Fatal, 'GSMDia2 must be greater than zero when GSPotent or GSShadow is enabled; check row '//trim(num2lstr(j)), ErrStat, ErrMsg, RoutineName ); if(Failed()) return + end if + if (InputFileData%GS%InpMembers(j)%MCd1<0.0_ReKi) then + call SetErrStat( ErrID_Fatal, 'GSMCd1 cannot be negative; check row '//trim(num2lstr(j)), ErrStat, ErrMsg, RoutineName ); if(Failed()) return + end if + if (InputFileData%GS%InpMembers(j)%MCd2<0.0_ReKi) then + call SetErrStat( ErrID_Fatal, 'GSMCd2 cannot be negative; check row '//trim(num2lstr(j)), ErrStat, ErrMsg, RoutineName ); if(Failed()) return + end if + if (InputFileData%GS%InpMembers(j)%MTI1<0.0_ReKi) then + call SetErrStat( ErrID_Fatal, 'GSMTI1 cannot be negative; check row '//trim(num2lstr(j)), ErrStat, ErrMsg, RoutineName ); if(Failed()) return + end if + if (InputFileData%GS%InpMembers(j)%MTI2<0.0_ReKi) then + call SetErrStat( ErrID_Fatal, 'GSMTI2 cannot be negative; check row '//trim(num2lstr(j)), ErrStat, ErrMsg, RoutineName ); if(Failed()) return + end if + ! The Eames GS shadow model divides by GSMTI, so it must be strictly positive (mirrors the Eames tower-shadow limits). + if (InputFileData%GS%GSShadow == GSShadow_Eames) then + if ( InputFileData%GS%InpMembers(j)%MTI1<=0.05_ReKi .or. InputFileData%GS%InpMembers(j)%MTI1>=1.0_ReKi .or. & + InputFileData%GS%InpMembers(j)%MTI2<=0.05_ReKi .or. InputFileData%GS%InpMembers(j)%MTI2>=1.0_ReKi ) then + call SetErrStat( ErrID_Fatal, 'The turbulence intensity (GSMTI1/GSMTI2) for the Eames GS shadow model must be greater than 0.05 and less than 1; check row '//trim(num2lstr(j)), ErrStat, ErrMsg, RoutineName ); if(Failed()) return + end if + if ( (InputFileData%GS%InpMembers(j)%MTI1>0.4_ReKi .and. InputFileData%GS%InpMembers(j)%MTI1<1.0_ReKi) .or. & + (InputFileData%GS%InpMembers(j)%MTI2>0.4_ReKi .and. InputFileData%GS%InpMembers(j)%MTI2<1.0_ReKi) ) then + call SetErrStat( ErrID_Warn, 'The turbulence intensity (GSMTI1/GSMTI2) for the Eames GS shadow model above 0.4 may return unphysical results. Interpret with caution; check row '//trim(num2lstr(j)), ErrStat, ErrMsg, RoutineName ) + end if + end if + end do + contains SUBROUTINE Fatal(ErrMsg_in) @@ -5123,7 +5412,173 @@ subroutine Cleanup() end subroutine Cleanup END SUBROUTINE Init_BEMTmodule +!---------------------------------------------------------------------------------------------------------------------------------- +!> This routine initializes the FVW module from within AeroDyn. +SUBROUTINE Init_GSParam( InputFileData, p, Density, MHK, WtrDpth, ErrStat, ErrMsg ) + type(GSInputFile), intent(in ) :: InputFileData + type(GSParameterType), intent(inout) :: p + real(ReKi), intent(in ) :: Density + integer(IntKi), intent(in ) :: MHK + real(ReKi), intent(in ) :: WtrDpth + integer(IntKi), intent( out) :: errStat !< Error status of the operation + character(*), intent( out) :: errMsg !< Error message if ErrStat /= ErrID_None + integer(IntKi) :: iMem, iJoint, iNode + integer(IntKi) :: j1, j2, j1Indx, j2Indx + integer(IntKi) :: memID + integer(IntKi) :: numDiv + real(ReKi) :: j1Pos(3), j2Pos(3) + real(ReKi) :: memLength, s + logical, allocatable :: JointUsed(:) + INTEGER(IntKi) :: ErrStat2 + CHARACTER(ErrMsgLen) :: ErrMsg2 + character(*), parameter :: RoutineName = 'Init_GSParam' + + ErrStat = ErrID_None + ErrMsg = "" + + if (InputFileData%NMembers==0_IntKi) then + p%GSPotent = GSPotent_none + p%GSShadow = GSShadow_none + p%GSAero = GSAero_none + p%hasGSMod = .false. + p%Density = 0.0_ReKi + p%MHK = MHK + p%WtrDpth = WtrDpth + p%NJoints = 0_IntKi + p%NNodes = 0_IntKi + p%NMembers = 0_IntKi + return + end if + + p%GSPotent = InputFileData%GSPotent + p%GSShadow = InputFileData%GSShadow + p%GSAero = InputFileData%GSAero + if ((p%GSPotent/=GSPotent_none) .or. (p%GSShadow/=GSShadow_none) .or. (p%GSAero/=GSAero_none)) then + p%hasGSMod = .TRUE. + else + p%hasGSMod = .FALSE. + endif + p%Density = Density + p%NMembers = InputFileData%NMembers + p%NJoints = InputFileData%NJoints + p%NNodes = InputFileData%NJoints ! Total number of nodes (joints + interior nodes); incremented below + p%MHK = MHK + p%WtrDpth = WtrDpth + + allocate( p%Members(InputFileData%NMembers), STAT = ErrStat2) + if ( ErrStat2 /= 0 ) then + errMsg = RoutineName//': Error allocating space for the members array. ' + errStat = ErrID_Fatal + call Cleanup() + return + end if + allocate( p%Joints(InputFileData%NJoints), STAT = ErrStat2) + if ( ErrStat2 /= 0 ) then + errMsg = RoutineName//': Error allocating space for the joints array. ' + errStat = ErrID_Fatal + call Cleanup() + return + end if + call AllocAry( JointUsed, InputFileData%NJoints, 'temporary logical array to track which joint is used', ErrStat2, ErrMsg2 ) + if (Failed()) return + JointUsed = .false. + + do iJoint = 1,InputFileData%NJoints + p%Joints(iJoint)%JointID = InputFileData%InpJoints(iJoint)%JointID + p%Joints(iJoint)%position = InputFileData%InpJoints(iJoint)%position + if (MHK == MHK_FixedBottom) then + p%Joints(iJoint)%position(3) = p%Joints(iJoint)%position(3) - WtrDpth + end if + end do + + do iMem = 1,InputFileData%NMembers + memID = InputFileData%InpMembers(iMem)%MemberID + p%Members(iMem)%MemberID = memID + + j1 = InputFileData%InpMembers(iMem)%MJointID1 + j2 = InputFileData%InpMembers(iMem)%MJointID2 + j1Indx = -1_IntKi + j2Indx = -1_IntKi + do iJoint = 1,InputFileData%NJoints + if (InputFileData%InpJoints(iJoint)%JointID == j1) j1Indx = iJoint + if (InputFileData%InpJoints(iJoint)%JointID == j2) j2Indx = iJoint + enddo + if (j1Indx<0_IntKi) then + ErrStat = ErrID_Fatal + ErrMsg = RoutineName//': MJointID1 of member with ID '//trim(num2lstr(memID))//' not found in the joint table. ' + call Cleanup() + return + endif + JointUsed(j1Indx) = .true. + if (j2Indx<0_IntKi) then + ErrStat = ErrID_Fatal + ErrMsg = RoutineName//': MJointID2 of member with ID '//trim(num2lstr(memID))//' not found in the joint table. ' + call Cleanup() + return + endif + JointUsed(j2Indx) = .true. + + j1Pos = InputFileData%InpJoints(j1Indx)%Position + j2Pos = InputFileData%InpJoints(j2Indx)%Position + memLength = TwoNorm(j2Pos-j1Pos) + if (memLength==0.0_ReKi) then + ErrStat = ErrID_Fatal + ErrMsg = RoutineName//': Member with ID '//trim(num2lstr(memID))//' has zero length. ' + call Cleanup() + return + endif + if (InputFileData%InpMembers(iMem)%MDivSize <= 0.0_ReKi) then + ErrStat = ErrID_Fatal + ErrMsg = RoutineName//': GSMDiv (MDivSize) must be > 0 for member with ID '//trim(num2lstr(memID))//'. ' + call Cleanup() + return + endif + numDiv = ceiling( memLength / InputFileData%InpMembers(iMem)%MDivSize ) + p%members(iMem)%RefLength = memLength + p%members(iMem)%NElements = numDiv + p%members(iMem)%dl = memLength/numDiv + + call AllocAry( p%members(iMem)%NodeIndx, numDiv+1_IntKi, 'Index of each of the member nodes in the global node list', ErrStat2, ErrMsg2 ) + if (Failed()) return + p%members(iMem)%NodeIndx(1) = j1Indx + p%members(iMem)%NodeIndx(numDiv+1_IntKi) = j2Indx + do iNode = 2,numDiv + p%members(iMem)%NodeIndx(iNode) = p%NNodes + iNode - 1_IntKi + enddo + p%NNodes = p%NNodes + numDiv - 1_IntKi + + call AllocAry( p%members(iMem)%R, numDiv+1_IntKi, 'Radius of each member nodes', ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( p%members(iMem)%Cd, numDiv+1_IntKi, 'Drag coefficient of each member nodes', ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( p%members(iMem)%TI, numDiv+1_IntKi, 'Turbulence intensity of each member nodes', ErrStat2, ErrMsg2 ); if (Failed()) return + do iNode = 1,numDiv+1_IntKi + s = real(iNode-1_IntKi,ReKi) / real(numDiv,ReKi) + p%members(iMem)%R(iNode) = (InputFileData%InpMembers(iMem)%MDiam1 * (1.0_ReKi-s) + InputFileData%InpMembers(iMem)%MDiam2 * s) * 0.5_ReKi + p%members(iMem)%Cd(iNode) = InputFileData%InpMembers(iMem)%MCd1 * (1.0_ReKi-s) + InputFileData%InpMembers(iMem)%MCd2 * s + p%members(iMem)%TI(iNode) = InputFileData%InpMembers(iMem)%MTI1 * (1.0_ReKi-s) + InputFileData%InpMembers(iMem)%MTI2 * s + enddo + end do + + if ( any(.not.JointUsed) ) then + ErrStat = ErrID_Fatal + ErrMsg = RoutineName//': All joints in the general support structure joint table must be referenced by at least one member. ' + call Cleanup() + return + end if + + call Cleanup() + +contains + subroutine Cleanup() + if (allocated(JointUsed)) deallocate(JointUsed) + end subroutine Cleanup + + logical function Failed() + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + Failed = ErrStat >= AbortErrLev + if (Failed) call CleanUp() + end function Failed +END SUBROUTINE !---------------------------------------------------------------------------------------------------------------------------------- !> This routine initializes the FVW module from within AeroDyn. SUBROUTINE Init_OLAF( InputFileData, u_AD, u, p, x, xd, z, OtherState, m, ErrStat, ErrMsg ) @@ -5463,6 +5918,65 @@ SUBROUTINE ADTwr_CalcOutput(p, u, RotInflow, m, y, ErrStat, ErrMsg ) END SUBROUTINE ADTwr_CalcOutput !---------------------------------------------------------------------------------------------------------------------------------- +!> This subroutine calculates the general support structure loads for the AeroDyn GSLoad output mesh. +SUBROUTINE ADGS_CalcOutput(p, u, GSInflow, m, y, ErrStat, ErrMsg ) + + TYPE(RotInputType), INTENT(IN ) :: u !< Inputs at Time t + TYPE(ElemInflowType), INTENT(IN ) :: GSInflow !< Inflow at Time t + TYPE(GSParameterType), INTENT(IN ) :: p !< Parameters + TYPE(RotMiscVarType), INTENT(INOUT) :: m !< Misc/optimization variables + TYPE(RotOutputType), INTENT(INOUT) :: y !< Outputs computed at t + INTEGER(IntKi), INTENT( OUT) :: ErrStat !< Error status of the operation + CHARACTER(*), INTENT( OUT) :: ErrMsg !< Error message if ErrStat /= ErrID_None + + INTEGER(IntKi) :: j, j1, j2, iMem + real(ReKi) :: q1(3), q2(3) + real(ReKi) :: V_rel(3), V_mag ! relative wind speed on a tower node + real(ReKi) :: pos1(3) + real(ReKi) :: pos2(3) + real(ReKi) :: k(3) + real(ReKi) :: dl + character(*), parameter :: RoutineName = 'ADGS_CalcOutput' + + ErrStat = ErrID_None + ErrMsg = "" + + IF (p%GSAero==GSAero_none) RETURN + + y%GSLoad%Force = 0.0_ReKi + y%GSLoad%moment = 0.0_ReKi + do iMem = 1,p%NMembers + associate( mem => p%Members(iMem)) + do j=1,mem%NElements + + j1 = mem%NodeIndx(j ) + j2 = mem%NodeIndx(j+1) + pos1 = u%GSMotion%position(:,j1) + u%GSMotion%TranslationDisp(:,j1) + pos2 = u%GSMotion%position(:,j2) + u%GSMotion%TranslationDisp(:,j2) + k = pos2-pos1 + dl = TwoNorm(k) + if (dl==0.0_ReKi) cycle ! Zero length element; shouldn't happen, but just in case. + k = k/dl + + V_rel = GSInflow%InflowVel(:,j1) - u%GSMotion%TranslationVel(:,j1) ! Total relative wind speed at GS node + V_rel = V_rel - dot_product(V_rel,k)*k ! Transverse component of relative wind speed + V_mag = TwoNorm(V_rel) + q1 = mem%Cd(j ) * p%Density * mem%R(j ) * 0.5_ReKi * dl * V_mag * V_rel + + V_rel = GSInflow%InflowVel(:,j2) - u%GSMotion%TranslationVel(:,j2) ! Total relative wind speed at GS node + V_rel = V_rel - dot_product(V_rel,k)*k ! Transverse component of relative wind speed + V_mag = TwoNorm(V_rel) + q2 = mem%Cd(j+1) * p%Density * mem%R(j+1) * 0.5_ReKi * dl * V_mag * V_rel + + y%GSLoad%force(:,j1) = y%GSLoad%force(:,j1) + q1 + y%GSLoad%force(:,j2) = y%GSLoad%force(:,j2) + q2 + + end do + end associate + end do + +END SUBROUTINE ADGS_CalcOutput +!---------------------------------------------------------------------------------------------------------------------------------- !> This routine checks for invalid inputs to the tower influence models. SUBROUTINE CheckTwrInfl(u, ErrStat, ErrMsg ) @@ -5755,6 +6269,10 @@ SUBROUTINE getLocalTowerProps(p, u, RotInflow, BladeNodePosition, theta_tower_tr end if + ! NOTE (potential fix): the 20-diameter far-field limit is a numerical safeguard, not physics. The + ! underlying shadow is a frozen 2-D wake carried in the member-normal plane, so at many diameters it is + ! unphysical (no 3-D recovery/advection) and the hard cutoff adds a small step in induced velocity. + ! Revisit: wind-aligned wake advection and/or a smooth far-field taper instead of the hard 20*Diam step. if ( TwrClrnc>20.0_ReKi*TwrDiam) then ! Far away, we skip the computation and keep undisturbed inflow DisturbInflow = .false. @@ -5850,7 +6368,7 @@ SUBROUTINE TwrInfl_NearestLine2Element(p, u, RotInflow, BladeNodePosition, r_Tow else !we're not on the element on_element = .false. end if - + end if if (on_element) then @@ -5960,7 +6478,7 @@ SUBROUTINE TwrInfl_NearestPoint(p, u, RotInflow, BladeNodePosition, r_TowerBlade ! calculate distance between points (note: this is actually the distance squared); ! will only store information once we have determined the closest node - r_TowerBlade = BladeNodePosition - p1 + r_TowerBlade = BladeNodePosition - p1 dist = dot_product( r_TowerBlade, r_TowerBlade ) if (dist .lt. min_dist) then @@ -6036,7 +6554,642 @@ SUBROUTINE TwrInfl_NearestPoint(p, u, RotInflow, BladeNodePosition, r_TowerBlade END SUBROUTINE TwrInfl_NearestPoint !---------------------------------------------------------------------------------------------------------------------------------- +!---------------------------------------------------------------------------------------------------------------------------------- +!> This routine checks for invalid inputs to the generalized support structure (GS) influence models. +!! Mirrors CheckTwrInfl, but the GS structure is a collection of members (not a single Line2 mesh), so +!! consecutive-node colocation is checked member-by-member. +SUBROUTINE CheckGSInfl(p_GS, u, ErrStat, ErrMsg ) + TYPE(GSParameterType), INTENT(IN ) :: p_GS !< General support structure parameters + TYPE(RotInputType), INTENT(IN ) :: u !< Inputs at Time t + INTEGER(IntKi), INTENT( OUT) :: ErrStat !< Error status of the operation + CHARACTER(*), INTENT( OUT) :: ErrMsg !< Error message if ErrStat /= ErrID_None + + ! local variables + real(reKi) :: ElemSize + real(reKi) :: tmp(3) + integer(intKi) :: iMem, j, n1, n2 + character(*), parameter :: RoutineName = 'CheckGSInfl' + + + ErrStat = ErrID_None + ErrMsg = "" + + !! the GS-influence models (potential flow and shadow) are valid only for small deflections; + !! so, first throw an error to avoid a division-by-zero error if any two consecutive nodes of a member are colocated. + + do iMem = 1, p_GS%NMembers + associate( mem => p_GS%Members(iMem) ) + do j = 2, mem%NElements+1 + n1 = mem%NodeIndx(j-1) + n2 = mem%NodeIndx(j ) + tmp = u%GSMotion%Position(:,n2) + u%GSMotion%TranslationDisp(:,n2) & + - u%GSMotion%Position(:,n1) - u%GSMotion%TranslationDisp(:,n1) + ElemSize = TwoNorm(tmp) + if ( EqualRealNos(ElemSize,0.0_ReKi) ) then + call SetErrStat(ErrID_Fatal, "Division by zero: nodes "//trim(num2lstr(n1))//' and '//trim(num2lstr(n2))// & + ' of member '//trim(num2lstr(mem%MemberID))//' are colocated.', ErrStat, ErrMsg, RoutineName ) + exit + end if + end do + end associate + end do + + +END SUBROUTINE CheckGSInfl +!---------------------------------------------------------------------------------------------------------------------------------- +!> This routine calculates the influence of the generalized support structure (GS) on the blade inflow, +!! following the same modeling approach as TwrInfl. The disturbance computed here is ADDED to whatever is +!! already stored in m%DisturbedInflow (e.g., tower influence), rather than overwriting it. +SUBROUTINE GSInfl( p, p_GS, u, GSInflow, m, ErrStat, ErrMsg ) +!.................................................................................................................................. + + TYPE(RotInputType), INTENT(IN ) :: u !< Inputs at Time t + TYPE(RotParameterType), INTENT(IN ) :: p !< Parameters (rotor-level GSPotent/GSShadow/hasGSMod copies) + TYPE(GSParameterType), INTENT(IN ) :: p_GS !< General support structure parameters + TYPE(ElemInflowType), INTENT(IN ) :: GSInflow !< Inflow on the GS structure at Time t + type(RotMiscVarType), intent(inout) :: m !< Misc/optimization variables + INTEGER(IntKi), INTENT( OUT) :: ErrStat !< Error status of the operation + CHARACTER(*), INTENT( OUT) :: ErrMsg !< Error message if ErrStat /= ErrID_None + + ! local variables + real(ReKi) :: xbar ! local x^ component of r_GSBlade (distance from GS structure to blade) normalized by GS diameter + real(ReKi) :: ybar ! local y^ component of r_GSBlade (distance from GS structure to blade) normalized by GS diameter + real(ReKi) :: zbar ! local z^ component of r_GSBlade (distance from GS structure to blade) normalized by GS diameter + real(ReKi) :: theta_GS_trans(3,3) ! transpose of local GS-member orientation expressed as a DCM + real(ReKi) :: GSCd ! local GS-member drag coefficient + real(ReKi) :: GSTI ! local GS-member TI (for Eames shadow model) + real(ReKi) :: W_GS ! local relative wind speed normal to the GS member + + real(ReKi) :: BladeNodePosition(3) ! local blade node position + + real(ReKi) :: v_pot(3) ! potential-flow contribution (member-local frame) for the current member + real(ReKi) :: v_shad(3) ! shadow contribution (member-local frame) for the current member + real(ReKi) :: v_pot_wsum(3) ! influence-weighted vector accumulator across members (global frame): sum_i w_i * v_i, w_i = |v_i|^p + real(ReKi) :: v_pot_wtot ! influence weight normalizer across members: sum_i w_i = sum_i |v_i|^p + real(ReKi) :: w_pot ! current member's potential-flow influence weight |v_pot|^p + real(ReKi) :: v_shad_sum(3) ! vector sum of shadow contributions across members (global frame; sets the combined deficit direction) + real(ReKi) :: v_shad_sumsq ! sum over members of squared shadow deficit magnitude (sets the root-sum-square magnitude) + real(ReKi) :: v_shad_norm ! magnitude of v_shad_sum + real(ReKi) :: v_result(3) ! combined (potential + shadow) disturbed velocity contribution + real(ReKi), parameter :: GSBlendExp = 6.0_ReKi ! partition-of-unity blend exponent p: p->inf recovers a hard argmax (one body only); finite p>=2 smooths the handoff. Even integer => weight (v.v)^(p/2) is a polynomial (C-infinity everywhere); p=6 is the smallest even p keeping the worst-case collinear-joint dip below 5%. + + real(ReKi) :: GSClrnc ! local GS clearance + logical :: FirstWarn_GSStrike + logical :: DisturbInflow + + integer(IntKi) :: j, k, iMem ! loop counters for elements, blades, and GS members + integer(intKi) :: ErrStat2 + character(ErrMsgLen) :: ErrMsg2 + character(*), parameter :: RoutineName = 'GSInfl' + + + ErrStat = ErrID_None + ErrMsg = "" + + if (.not. p%hasGSMod) return + if (p%GSPotent == GSPotent_none .and. p%GSShadow == GSShadow_none) return + + FirstWarn_GSStrike = .true. + + ! these models are valid for only small deflections; check for potential division-by-zero errors: + call CheckGSInfl( p_GS, u, ErrStat2, ErrMsg2 ) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + if (ErrStat >= AbortErrLev) return + + do k = 1, p%NumBlades + do j = 1, u%BladeMotion(k)%NNodes + + ! for each line2-element node of the blade mesh, a nearest-neighbor line2 element or node of the GS + ! structure is found in the deflected configuration, returning theta_GS_trans, W_GS, xbar, ybar, zbar, and GSCd: + + BladeNodePosition = u%BladeMotion(k)%Position(:,j) + u%BladeMotion(k)%TranslationDisp(:,j) + + ! reset per-member accumulators for this blade node: the potential effect is combined by an + ! influence-weighted partition of unity (each body carries its own field; the blend picks the + ! locally dominant one and never exceeds it), the shadow effect by root-sum-square + v_pot_wsum = 0.0_ReKi + v_pot_wtot = 0.0_ReKi + v_shad_sum = 0.0_ReKi + v_shad_sumsq = 0.0_ReKi + m%GSClrnc(j,k) = HUGE(1.0_ReKi) ! reduced to the nearest-member clearance in the loop below + + ! Loop through the GS members + do iMem = 1, p_GS%NMembers + call getLocalGSProps(p_GS, iMem, u, GSInflow, BladeNodePosition, theta_GS_trans, W_GS, xbar, ybar, zbar, GSCd, GSTI, GSClrnc, FirstWarn_GSStrike, DisturbInflow, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + if (.not. FirstWarn_GSStrike) call SetErrStat(ErrID_Fatal, "Generalized support structure strike.", ErrStat, ErrMsg, RoutineName ) + if (ErrStat >= AbortErrLev) return + + m%GSClrnc(j,k) = min( m%GSClrnc(j,k), GSClrnc ) ! nearest-member clearance, saved for output + + if ( DisturbInflow ) then + ! get the potential-flow and shadow contributions separately (member-local frame) + call CalculateGSInfluence(p, xbar, ybar, zbar, W_GS, GSCd, GSTI, v_pot, v_shad) + ! rotate each contribution into global coordinates and accumulate + v_pot = matmul( theta_GS_trans, v_pot ) + v_shad = matmul( theta_GS_trans, v_shad ) + ! potential flow: influence-weighted blend. Weight each body by its own field magnitude + ! |v_pot|^p (p=GSBlendExp), so the blade sees essentially the single strongest body's + ! field, smoothly handed off between bodies. |v_pot|^2 = dot(v_pot,v_pot) is C-infinity, + ! so raising it to p/2 is smooth; a body whose field has died out (|v_pot|->0) gets ~0 + ! weight and cannot blank the field. + w_pot = dot_product(v_pot,v_pot) ** (0.5_ReKi*GSBlendExp) + v_pot_wsum = v_pot_wsum + w_pot * v_pot + v_pot_wtot = v_pot_wtot + w_pot + v_shad_sum = v_shad_sum + v_shad ! shadow: vector sum sets the combined deficit direction + v_shad_sumsq = v_shad_sumsq + dot_product(v_shad,v_shad) ! shadow: accumulate squared magnitude for root-sum-square + end if + enddo !iMem + + v_result = 0.0_ReKi + if ( p%GSPotent /= GSPotent_none ) then + ! potential-flow effect combines by an influence-weighted partition of unity: + ! v = ( sum_i |v_i|^p v_i ) / ( sum_i |v_i|^p ). The weights sum to one, so the blended + ! magnitude never exceeds max_i|v_i| (no junction double-counting), it reduces to the + ! single-member result when only one body contributes, and it needs no knowledge of the + ! structure topology (collinear/angled/branching are all handled by field magnitude alone). + if ( v_pot_wtot > 0.0_ReKi ) then + v_result = v_result + v_pot_wsum / v_pot_wtot + end if + endif + if ( p%GSShadow /= GSShadow_none ) then + ! shadow effect combines by root-sum-square of each member's deficit magnitude, applied + ! along the direction of the vector sum of the member deficits. This is a pure vector + ! operation (frame independent): it makes no assumption about the member-local axes and + ! reduces to the single-member result when only one member contributes. + v_shad_norm = TwoNorm( v_shad_sum ) + if ( v_shad_norm > 0.0_ReKi ) then + v_result = v_result + sqrt( v_shad_sumsq ) * v_shad_sum / v_shad_norm + end if + endif + + ! Combine the per-member velocity contributions (via v_result above) into the disturbed inflow. + ! NOTE (document later): GS members are combined among themselves by the partition-of-unity / + ! root-sum-square blend above, but the GS effect is superimposed on the tower effect by simple + ! addition here (no cross-blend between the tower and GS fields). + m%DisturbedInflow(:,j,k) = m%DisturbedInflow(:,j,k) + v_result + + end do !j=NumBlNds + end do ! NumBlades + + +END SUBROUTINE GSInfl +!---------------------------------------------------------------------------------------------------------------------------------- +!> Mirrors CalculateTowerInfluence, computing the potential-flow / shadow velocity deficit induced by a +!! generalized support structure (GS) member using the GSPotent / GSShadow models. The potential-flow and +!! shadow contributions are returned separately (in the member-local frame) so the caller can combine the +!! potential effect across members by an influence-weighted blend and the shadow effect by root-sum-square. +SUBROUTINE CalculateGSInfluence(p, xbar_in, ybar, zbar, W_GS, GSCd, GSTI, v_pot, v_shad) + + TYPE(RotParameterType), INTENT(IN ) :: p !< Parameters + real(ReKi), intent(in) :: xbar_in ! local x^ component of r_GSBlade normalized by GS diameter + real(ReKi), intent(in) :: ybar ! local y^ component of r_GSBlade normalized by GS diameter + real(ReKi), intent(in) :: zbar ! local z^ component of r_GSBlade normalized by GS diameter + real(ReKi), intent(in) :: W_GS ! local relative wind speed normal to the GS member + real(ReKi), intent(in) :: GSCd ! local GS-member drag coefficient + real(ReKi), intent(in) :: GSTI ! local GS-member TI (for Eames shadow model) + real(ReKi), intent( out) :: v_pot(3) ! potential-flow velocity contribution (member-local frame) + real(ReKi), intent( out) :: v_shad(3) ! shadow velocity contribution (member-local frame) + + real(ReKi) :: denom ! denominator + real(ReKi) :: exponential ! exponential term + real(ReKi) :: xbar ! potentially modified version of xbar_in + real(ReKi) :: u_GSShadow ! axial velocity deficit fraction from GS shadow + real(ReKi) :: u_GSPotent ! axial velocity deficit fraction from GS potential flow + real(ReKi) :: v_GSPotent ! transverse velocity deficit fraction from GS potential flow + + + u_GSShadow = 0.0_ReKi + u_GSPotent = 0.0_ReKi + v_GSPotent = 0.0_ReKi + xbar = xbar_in + + ! calculate GS-member influence: + if ( abs(zbar) < 1.0_ReKi .and. p%GSPotent /= GSPotent_none ) then + + if ( p%GSPotent == GSPotent_baseline ) then + denom = (xbar**2 + ybar**2)**2 + u_GSPotent = ( -1.0*xbar**2 + ybar**2 ) / denom + v_GSPotent = ( -2.0*xbar * ybar ) / denom + + elseif (p%GSPotent == GSPotent_Bak) then + ! Reference: Bak, Madsen, Johansen (2001): Influence from Blade-Tower Interaction on Fatigue Loads and Dynamics (poster); + ! Proceedings: EWEC'01; Copenhagen (DK) + xbar = xbar + 0.1 ! offset added as part of the original model of Bak et al. + denom = (xbar**2 + ybar**2)**2 + u_GSPotent = ( -1.0*xbar**2 + ybar**2 ) / denom + v_GSPotent = ( -2.0*xbar * ybar ) / denom + denom = TwoPi*(xbar**2 + ybar**2) + u_GSPotent = u_GSPotent + GSCd*xbar / denom + v_GSPotent = v_GSPotent + GSCd*ybar / denom + xbar = xbar - 0.1 ! removing offset + + end if + end if + + ! NOTE (limitation): the wake is directed along x^ = the wind projected into the plane normal to the + ! member axis, not along the true earth-fixed wind. For a member inclined into or away from the wind + ! this tilts the wake up into the sky or down into the ground rather than keeping it parallel to the + ! ground, which is unphysical. Acceptable for near-vertical members; revisit for strongly inclined ones. + ! Intended fix: advect the wake along the incident wind direction (ideally low-pass filtered), since the + ! wake is advective; the potential-flow (kinematic) part above rightly stays in the member normal plane. + select case (p%GSShadow) + case (GSShadow_Powles) + if ( xbar > 0.0_ReKi .and. abs(zbar) < 1.0_ReKi) then + denom = sqrt( sqrt( xbar**2 + ybar**2 ) ) + if ( abs(ybar) < denom ) then + u_GSShadow = -GSCd / denom * cos( PiBy2*ybar / denom )**2 + end if + end if + case (GSShadow_Eames) + if ( xbar > 0.0_ReKi .and. abs(zbar) < 1.0_ReKi) then + exponential = ( ybar / (GSTI * xbar) )**2 + denom = GSTI * xbar * sqrt( TwoPi ) + u_GSShadow = -GSCd / denom * exp ( -0.5_ReKi * exponential ) + end if + end select + + ! We limit the deficit to avoid having too much flow reversal and accumulation of vorticity behind the member + ! Limit to -0.5 the wind speed at the GS member + u_GSShadow = max(u_GSShadow, -0.5_ReKi) + + + v_pot(1) = u_GSPotent*W_GS + v_pot(2) = v_GSPotent*W_GS + v_pot(3) = 0.0_ReKi + + v_shad(1) = u_GSShadow*W_GS + v_shad(2) = 0.0_ReKi + v_shad(3) = 0.0_ReKi + + +END SUBROUTINE CalculateGSInfluence +!---------------------------------------------------------------------------------------------------------------------------------- +!> This routine returns the GS-member constants necessary to compute the GS influence, mirroring getLocalTowerProps. +!! if u%GSMotion does not have any nodes there will be serious problems. I assume that has been checked earlier. +SUBROUTINE getLocalGSProps(p_GS, iMem, u, GSInflow, BladeNodePosition, theta_GS_trans, W_GS, xbar, ybar, zbar, GSCd, GSTI, GSClrnc, FirstWarn_GSStrike, DisturbInflow, ErrStat, ErrMsg) +!.................................................................................................................................. + TYPE(RotInputType), INTENT(IN ) :: u !< Inputs at Time t + TYPE(GSParameterType), INTENT(IN ) :: p_GS !< General support structure parameters + TYPE(ElemInflowType), INTENT(IN ) :: GSInflow !< Inflow on the GS structure at Time t + INTEGER(IntKi) ,INTENT( IN) :: iMem !< GS member index + REAL(ReKi) ,INTENT(IN ) :: BladeNodePosition(3) !< local blade node position + REAL(ReKi) ,INTENT( OUT) :: theta_GS_trans(3,3) !< transpose of local GS-member orientation expressed as a DCM + LOGICAL ,INTENT(INOUT) :: FirstWarn_GSStrike !< Whether we should check and warn for a GS-structure strike + LOGICAL ,INTENT( OUT) :: DisturbInflow !< Whether GS clearance is in the range of values where it should disturb the inflow + REAL(ReKi) ,INTENT( OUT) :: W_GS !< local relative wind speed normal to the GS member + REAL(ReKi) ,INTENT( OUT) :: xbar !< local x^ component of r_GSBlade normalized by GS diameter + REAL(ReKi) ,INTENT( OUT) :: ybar !< local y^ component of r_GSBlade normalized by GS diameter + REAL(ReKi) ,INTENT( OUT) :: zbar !< local z^ component of r_GSBlade normalized by GS diameter + REAL(ReKi) ,INTENT( OUT) :: GSCd !< local GS-member drag coefficient + REAL(ReKi) ,INTENT( OUT) :: GSTI !< local GS-member TI (for Eames shadow model) + REAL(ReKi) ,INTENT( OUT) :: GSClrnc !< GS clearance for potential output + INTEGER(IntKi), INTENT( OUT) :: ErrStat !< Error status of the operation + CHARACTER(*), INTENT( OUT) :: ErrMsg !< Error message if ErrStat /= ErrID_None + + ! local variables + real(ReKi) :: r_GSBlade(3) ! distance vector from GS member to blade + real(ReKi) :: GSDiam ! local GS-member diameter + logical :: found + character(*), parameter :: RoutineName = 'getLocalGSProps' + + + ErrStat = ErrID_None + ErrMsg = "" + + ! .............................................. + ! option 1: nearest line2 element + ! .............................................. + call GSInfl_NearestLine2Element(p_GS, iMem, u, GSInflow, BladeNodePosition, r_GSBlade, theta_GS_trans, W_GS, xbar, ybar, zbar, GSCd, GSTI, GSDiam, found) + + if ( .not. found) then + ! .............................................. + ! option 2: nearest node + ! .............................................. + call GSInfl_NearestPoint(p_GS, iMem, u, GSInflow, BladeNodePosition, r_GSBlade, theta_GS_trans, W_GS, xbar, ybar, zbar, GSCd, GSTI, GSDiam) + + end if + + GSClrnc = TwoNorm(r_GSBlade) - 0.5_ReKi*GSDiam + + if (FirstWarn_GSStrike) then + if ( GSClrnc <= 0.0_ReKi ) then + call WrScr( NewLine//NewLine//"** WARNING: Generalized support structure strike. ** This warning will not be repeated though the condition may persist."//NewLine//NewLine ) + FirstWarn_GSStrike = .false. + end if + end if + + + ! NOTE (potential fix): the 20-diameter far-field limit is a numerical safeguard, not physics. The + ! underlying shadow is a frozen 2-D wake carried in the member-normal plane (see limitation note in + ! CalculateGSInfluence), so at many diameters it is unphysical (no 3-D recovery/advection) and the hard + ! cutoff adds a small step in induced velocity. Revisit: wind-aligned wake advection and/or a smooth + ! far-field taper instead of the hard 20*Diam step. + if ( GSClrnc>20.0_ReKi*GSDiam) then + ! Far away, we skip the computation and keep undisturbed inflow + DisturbInflow = .false. + elseif ( GSClrnc<=0.01_ReKi*GSDiam) then + ! Inside the member, or very close, (will happen for vortex elements) we keep undisturbed inflow + ! We don't want to reach the stagnation points + DisturbInflow = .false. + else + DisturbInflow = .true. + end if + +END SUBROUTINE getLocalGSProps +!---------------------------------------------------------------------------------------------------------------------------------- +!> Option 1: Find the nearest-neighbor line2 element (within any member) of the GS structure for which the blade +!! line2-element node projects orthogonally onto the member's line2-element domain. Mirrors TwrInfl_NearestLine2Element, +!! but loops over each member's own element list (mem%NodeIndx) instead of a single Line2 ElemTable, since the GS +!! structure may consist of several members that do not form one contiguous line. +SUBROUTINE GSInfl_NearestLine2Element(p_GS, iMem, u, GSInflow, BladeNodePosition, r_GSBlade, theta_GS_trans, W_GS, xbar, ybar, zbar, GSCd, GSTI, GSDiam, found) +!.................................................................................................................................. + TYPE(RotInputType), INTENT(IN ) :: u !< Inputs at Time t + TYPE(GSParameterType), INTENT(IN ) :: p_GS !< General support structure parameters + TYPE(ElemInflowType), INTENT(IN ) :: GSInflow !< Inflow on the GS structure at Time t + REAL(ReKi) ,INTENT(IN ) :: BladeNodePosition(3) !< local blade node position + INTEGER(IntKi) ,INTENT( IN) :: iMem !< GS member index + REAL(ReKi) ,INTENT( OUT) :: r_GSBlade(3) !< distance vector from GS member to blade + REAL(ReKi) ,INTENT( OUT) :: theta_GS_trans(3,3) !< transpose of local GS-member orientation expressed as a DCM + REAL(ReKi) ,INTENT( OUT) :: W_GS !< local relative wind speed normal to the GS member + REAL(ReKi) ,INTENT( OUT) :: xbar !< local x^ component of r_GSBlade normalized by GS diameter + REAL(ReKi) ,INTENT( OUT) :: ybar !< local y^ component of r_GSBlade normalized by GS diameter + REAL(ReKi) ,INTENT( OUT) :: zbar !< local z^ component of r_GSBlade normalized by GS diameter + REAL(ReKi) ,INTENT( OUT) :: GSCd !< local GS-member drag coefficient + REAL(ReKi) ,INTENT( OUT) :: GSTI !< local GS-member TI (Eames shadow model) + REAL(ReKi) ,INTENT( OUT) :: GSDiam !< local GS-member diameter + logical ,INTENT( OUT) :: found !< whether a mapping was found with this option + + ! local variables + REAL(ReKi) :: denom + REAL(ReKi) :: dist + REAL(ReKi) :: min_dist + REAL(ReKi) :: elem_position, elem_position2 + REAL(SiKi) :: elem_position_SiKi + + REAL(ReKi) :: p1(3), p2(3) ! position vectors for nodes on a GS-member line 2 element + + REAL(ReKi) :: V_rel_GS(3) + + REAL(ReKi) :: n1_n2_vector(3) ! vector going from node 1 to node 2 in Line2 element + REAL(ReKi) :: n1_Point_vector(3) ! vector going from node 1 in Line 2 element to Destination Point + REAL(ReKi) :: tmp(3) ! temporary vector for cross product calculation + + INTEGER(IntKi) :: jElem ! do-loop counter for elements within a member + + INTEGER(IntKi) :: n1, n2 ! global node indices associated with an element + + LOGICAL :: on_element + + + found = .false. + min_dist = HUGE(min_dist) + + associate( mem => p_GS%Members(iMem) ) + do jElem = 1, mem%NElements ! number of elements on this member + ! grab global node numbers associated with the jElem_th element of this member + n1 = mem%NodeIndx(jElem ) + n2 = mem%NodeIndx(jElem+1) + + p1 = u%GSMotion%Position(:,n1) + u%GSMotion%TranslationDisp(:,n1) + p2 = u%GSMotion%Position(:,n2) + u%GSMotion%TranslationDisp(:,n2) + + ! Calculate vectors used in projection operation + n1_n2_vector = p2 - p1 + n1_Point_vector = BladeNodePosition - p1 + + denom = DOT_PRODUCT( n1_n2_vector, n1_n2_vector ) ! we've already checked that these aren't zero + + ! project point onto line defined by n1 and n2 + + elem_position = DOT_PRODUCT(n1_n2_vector,n1_Point_vector) / denom + + ! note: i forumlated it this way because Fortran doesn't necessarially do shortcutting and I don't want to call EqualRealNos if we don't need it: + if ( elem_position .ge. 0.0_ReKi .and. elem_position .le. 1.0_ReKi ) then !we're ON the element (between the two nodes) + on_element = .true. + else + elem_position_SiKi = REAL( elem_position, SiKi ) + if (EqualRealNos( elem_position_SiKi, 1.0_SiKi )) then !we're ON the element (at a node) + on_element = .true. + elem_position = 1.0_ReKi + elseif (EqualRealNos( elem_position_SiKi, 0.0_SiKi )) then !we're ON the element (at a node) + on_element = .true. + elem_position = 0.0_ReKi + else !we're not on the element + on_element = .false. + end if + + end if + + if (on_element) then + + ! calculate distance between point and line (note: this is actually the distance squared); + ! will only store information once we have determined the closest element + elem_position2 = 1.0_ReKi - elem_position + + r_GSBlade = BladeNodePosition - elem_position2*p1 - elem_position*p2 + dist = dot_product( r_GSBlade, r_GSBlade ) + + if (dist .lt. min_dist) then + found = .true. + min_dist = dist + + V_rel_GS = ( GSInflow%InflowVel(:,n1) - u%GSMotion%TranslationVel(:,n1) ) * elem_position2 & + + ( GSInflow%InflowVel(:,n2) - u%GSMotion%TranslationVel(:,n2) ) * elem_position + + GSDiam = elem_position2*2.0_ReKi*mem%R(jElem) + elem_position*2.0_ReKi*mem%R(jElem+1) + GSCd = elem_position2*mem%Cd( jElem) + elem_position*mem%Cd( jElem+1) + GSTI = elem_position2*mem%TI( jElem) + elem_position*mem%TI( jElem+1) + + + ! z_hat + theta_GS_trans(:,3) = n1_n2_vector / sqrt( denom ) ! = n1_n2_vector / twoNorm( n1_n2_vector ) + + tmp = V_rel_GS - dot_product(V_rel_GS,theta_GS_trans(:,3)) * theta_GS_trans(:,3) + denom = TwoNorm( tmp ) + if (.not. EqualRealNos( denom, 0.0_ReKi ) ) then + ! x_hat + theta_GS_trans(:,1) = tmp / denom + + ! y_hat + tmp = cross_product( theta_GS_trans(:,3), V_rel_GS ) + theta_GS_trans(:,2) = tmp / denom + + W_GS = dot_product( V_rel_GS,theta_GS_trans(:,1) ) + xbar = 2.0/GSDiam * dot_product( r_GSBlade, theta_GS_trans(:,1) ) + ybar = 2.0/GSDiam * dot_product( r_GSBlade, theta_GS_trans(:,2) ) + zbar = 0.0_ReKi + + else + ! there is no GS influence because dot_product(V_rel_GS,x_hat) = 0 + ! thus, we don't need to set the other values (except we don't want the sum of xbar^2 and ybar^2 to be 0) + theta_GS_trans = 0.0_ReKi + W_GS = 0.0_ReKi + xbar = 1.0_ReKi + ybar = 0.0_ReKi + zbar = 0.0_ReKi + end if + + + end if !the point is closest to this line2 element + + end if + + end do !jElem + end associate + +END SUBROUTINE GSInfl_NearestLine2Element +!---------------------------------------------------------------------------------------------------------------------------------- +!> Option 2: used when the blade node does not orthogonally intersect a GS-member element. +!! Find the nearest-neighbor node within the member (for a straight member this is, by construction, a member +!! end; a strongly curved/deflected member could instead make an interior node nearest, which would then also +!! receive the end taper below -- acceptable for the near-rigid members assumed here). The true +!! axial offset zbar is used to apply the same cosine-squared taper as TwrInfl_NearestPoint (dividing +!! xbar,ybar by cos(PiBy2*zbar)), so the deficit decays smoothly to zero at |zbar|=1. The taper is applied +!! at every member end (junction or free tip), not only at structural free ends as in the tower model. +SUBROUTINE GSInfl_NearestPoint(p_GS, iMem, u, GSInflow, BladeNodePosition, r_GSBlade, theta_GS_trans, W_GS, xbar, ybar, zbar, GSCd, GSTI, GSDiam) +!.................................................................................................................................. + TYPE(RotInputType), INTENT(IN ) :: u !< Inputs at Time t + TYPE(GSParameterType), INTENT(IN ) :: p_GS !< General support structure parameters + TYPE(ElemInflowType), INTENT(IN ) :: GSInflow !< Inflow on the GS structure at Time t + REAL(ReKi) ,INTENT(IN ) :: BladeNodePosition(3) !< local blade node position + INTEGER(IntKi) ,INTENT( IN) :: iMem !< GS member index + REAL(ReKi) ,INTENT( OUT) :: r_GSBlade(3) !< distance vector from GS member to blade + REAL(ReKi) ,INTENT( OUT) :: theta_GS_trans(3,3) !< transpose of local GS-member orientation expressed as a DCM + REAL(ReKi) ,INTENT( OUT) :: W_GS !< local relative wind speed normal to the GS member + REAL(ReKi) ,INTENT( OUT) :: xbar !< local x^ component of r_GSBlade normalized by GS diameter + REAL(ReKi) ,INTENT( OUT) :: ybar !< local y^ component of r_GSBlade normalized by GS diameter + REAL(ReKi) ,INTENT( OUT) :: zbar !< local z^ component of r_GSBlade normalized by GS diameter + REAL(ReKi) ,INTENT( OUT) :: GSCd !< local GS-member drag coefficient + REAL(ReKi) ,INTENT( OUT) :: GSTI !< local GS-member TI (for Eames shadow model) + REAL(ReKi) ,INTENT( OUT) :: GSDiam !< local GS-member diameter + + ! local variables + REAL(ReKi) :: denom + REAL(ReKi) :: dist + REAL(ReKi) :: min_dist + REAL(ReKi) :: cosTaper + + REAL(ReKi) :: p1(3) ! position vectors for nodes on the GS structure + REAL(ReKi) :: V_rel_GS(3) + + REAL(ReKi) :: tmp(3) ! temporary vector for cross product calculation + REAL(ReKi) :: zaxis(3) ! member axis (from the element adjacent to the nearest node) + + INTEGER(IntKi) :: j ! do-loop counters for members and local member nodes + INTEGER(IntKi) :: n1 ! global node index + INTEGER(IntKi) :: nA, nB ! global node indices of the element adjacent to the nearest node + + INTEGER(IntKi) :: j_min, n1_min ! local-node/global-node with minimum distance found so far + + + + !................. + ! find the closest node (searching over each member's own node list) + !................. + + min_dist = HUGE(min_dist) + j_min = 0 + n1_min = 0 + + associate( mem => p_GS%Members(iMem) ) + do j = 1, mem%NElements+1 ! number of nodes on this member + + n1 = mem%NodeIndx(j) + p1 = u%GSMotion%Position(:,n1) + u%GSMotion%TranslationDisp(:,n1) + + ! calculate distance between points (note: this is actually the distance squared); + ! will only store information once we have determined the closest node + r_GSBlade = BladeNodePosition - p1 + dist = dot_product( r_GSBlade, r_GSBlade ) + + if (dist .lt. min_dist) then + min_dist = dist + j_min = j + n1_min = n1 + end if !the point is (so far) closest to this GS node + + end do !j + + !................. + ! calculate the values to be returned: + !.................. + if (n1_min == 0) then + j_min = 1 + n1_min = mem%NodeIndx(1) + if (NWTC_VerboseLevel == NWTC_Verbose) call WrScr( 'AD:GSInfl_NearestPoint:Error finding minimum distance. Positions may be invalid.' ) + end if + + n1 = n1_min + + r_GSBlade = BladeNodePosition - u%GSMotion%Position(:,n1) - u%GSMotion%TranslationDisp(:,n1) + V_rel_GS = GSInflow%InflowVel(:,n1) - u%GSMotion%TranslationVel(:,n1) + GSDiam = 2.0_ReKi * mem%R( j_min) + GSCd = mem%Cd( j_min) + GSTI = mem%TI( j_min) + + ! z_hat: use the actual member axis (mirroring GSInfl_NearestLine2Element) taken from the element + ! adjacent to the nearest node. The point-mesh Orientation is identity and would give (0,0,1), + ! which is wrong for inclined members and, by making the axial offset appear as an in-plane offset, + ! would spuriously apply full influence from members the blade node is actually axially far beyond. + if ( j_min == mem%NElements+1 ) then + nA = mem%NodeIndx(j_min-1) + nB = mem%NodeIndx(j_min) + else + nA = mem%NodeIndx(j_min) + nB = mem%NodeIndx(j_min+1) + end if + zaxis = ( u%GSMotion%Position(:,nB) + u%GSMotion%TranslationDisp(:,nB) ) & + - ( u%GSMotion%Position(:,nA) + u%GSMotion%TranslationDisp(:,nA) ) + theta_GS_trans(:,3) = zaxis / TwoNorm(zaxis) + + tmp = V_rel_GS - dot_product(V_rel_GS,theta_GS_trans(:,3)) * theta_GS_trans(:,3) + denom = TwoNorm( tmp ) + + if (.not. EqualRealNos( denom, 0.0_ReKi ) ) then + + ! x_hat + theta_GS_trans(:,1) = tmp / denom + + ! y_hat + tmp = cross_product( theta_GS_trans(:,3), V_rel_GS ) + theta_GS_trans(:,2) = tmp / denom + + W_GS = dot_product( V_rel_GS,theta_GS_trans(:,1) ) + + ! The nearest node is, by construction, a member end (the blade node projects beyond the member's + ! extent). Apply the same cosine-squared taper as TwrInfl_NearestPoint (option 2b): dividing xbar and + ! ybar by cos(PiBy2*zbar) makes the 1/r^2 doublet deficit scale as cos^2(PiBy2*zbar), so the influence + ! ramps smoothly to zero with zero slope at |zbar|=1. Unlike the tower (which tapers only at its free + ! top/bottom), we taper at every member end -- junction or free tip -- because each finite GS member's + ! doublet field must decay past its own end; a neighbouring member covers the region beyond a junction. + zbar = 2.0/GSDiam * dot_product( r_GSBlade, theta_GS_trans(:,3) ) + if (abs(zbar) < 1) then + cosTaper = cos( PiBy2*zbar ) + xbar = 2.0/GSDiam * dot_product( r_GSBlade, theta_GS_trans(:,1) ) / cosTaper + ybar = 2.0/GSDiam * dot_product( r_GSBlade, theta_GS_trans(:,2) ) / cosTaper + else ! zbar is checked (<1) before xbar,ybar are used downstream, but set them to safe values anyway + xbar = 1.0_ReKi + ybar = 0.0_ReKi + end if + + else + + ! there is no GS influence because W_GS = dot_product(V_rel_GS,x_hat) = 0 + ! thus, we don't need to set the other values (except we don't want the sum of xbar^2 and ybar^2 to be 0) + W_GS = 0.0_ReKi + theta_GS_trans = 0.0_ReKi + xbar = 1.0_ReKi + ybar = 0.0_ReKi + zbar = 0.0_ReKi + + end if + + end associate + +END SUBROUTINE GSInfl_NearestPoint +!---------------------------------------------------------------------------------------------------------------------------------- subroutine AD_InitVars(iR, u, p, x, z, OtherState, y, m, InitOut, InputFileData, Linearize, CompAeroMaps, ErrStat, ErrMsg) integer(IntKi), intent(in) :: iR !< Rotor number type(RotInputType), intent(inout) :: u !< An initial guess for the input; input mesh must be defined @@ -6175,6 +7328,13 @@ subroutine AD_InitVars(iR, u, p, x, z, OtherState, y, m, InitOut, InputFileData, Mesh=u%TowerMotion, & Perturbs=[PerturbTower, Perturb, PerturbTower, PerturbTower]) + ! Add general support motion (no orientation: the GS models only use translational quantities) + call MV_AddMeshVar(InitOut%Vars%u, "General Support", [FieldTransDisp, FieldTransVel], & + DatLoc(AD_u_GSMotion), & + Mesh=u%GSMotion, & + Perturbs=[PerturbTower, PerturbTower], & + Active=p%hasGSMod) + ! Add blade root motion do j = 1, p%NumBlades call MV_AddMeshVar(InitOut%Vars%u, "Blade root "//Num2LStr(j), [FieldOrientation], & @@ -6243,6 +7403,11 @@ subroutine AD_InitVars(iR, u, p, x, z, OtherState, y, m, InitOut, InputFileData, call MV_AddMeshVar(InitOut%Vars%y, "Tower", LoadFields, DatLoc(AD_y_TowerLoad), & Mesh=y%TowerLoad) + ! Add general support load + call MV_AddMeshVar(InitOut%Vars%y, "General Support", LoadFields, DatLoc(AD_y_GSLoad), & + Mesh=y%GSLoad, & + Active=(p%GSAero/=GSAero_none)) + ! Loop through blades, add blade loads do j = 1, p%NumBlades Flags = VF_Line @@ -6358,6 +7523,7 @@ SUBROUTINE AD_JacobianPInput(Vars, iRotor, t, u_AD, p_AD, x_AD, xd_AD, z_AD, Oth type(FlowFieldType),target :: FF_perturb type(FlowFieldType),pointer :: FF_ptr ! need a pointer in the CalcWind_Rotor routine type(RotInflowType) :: RotInflow_perturb !< Rotor inflow, perturbed by FlowField extended inputs + type(ElemInflowType) :: GSInflow_perturb ErrStat = ErrID_None ErrMsg = '' @@ -6373,7 +7539,8 @@ SUBROUTINE AD_JacobianPInput(Vars, iRotor, t, u_AD, p_AD, x_AD, xd_AD, z_AD, Oth z => z_AD%rotors(iRotor), & OtherState => OtherState_AD%rotors(iRotor), & y_lin => m_AD%y_lin%rotors(iRotor), & - RotInflow => m_AD%Inflow(1)%RotInflow(iRotor)) + RotInflow => m_AD%Inflow(1)%RotInflow(iRotor), & + GSInflow => m_AD%Inflow(1)%GSInflow) ! Find indices for extended input variables iVarHWindSpeed = 0 @@ -6392,6 +7559,8 @@ SUBROUTINE AD_JacobianPInput(Vars, iRotor, t, u_AD, p_AD, x_AD, xd_AD, z_AD, Oth call AD_CalcWind_Rotor(t, u, p_AD%FlowField, p, p_AD, m_AD, RotInflow, StartNode, ErrStat, ErrMsg) if (ErrStat >= AbortErrLev) return + call AD_CalcWind_GS(t, u, p_AD%FlowField, p_AD%GS, p_AD, m_AD, GSInflow, StartNode, ErrStat, ErrMsg) + if (ErrStat >= AbortErrLev) return ! If flow field will need to be perturbed (HWindSpeed, PLexp, or PropagationDir variables) if (iVarHWindSpeed > 0 .or. iVarPLexp > 0 .or. iVarPropagationDir > 0) then @@ -6405,7 +7574,7 @@ SUBROUTINE AD_JacobianPInput(Vars, iRotor, t, u_AD, p_AD, x_AD, xd_AD, z_AD, Oth ! Get OP values here (i.e., set inputs for BEMT): if (p%DBEMT_Mod == DBEMT_frozen) then - call SetInputs(t, p, p_AD, u, RotInflow, m, indx, errStat2, errMsg2); if (Failed()) return + call SetInputs(t, p, p_AD, u, RotInflow, GSInflow, m, indx, errStat2, errMsg2); if (Failed()) return ! compare m%BEMT_y arguments with call to BEMT_CalcOutput call computeFrozenWake(m%BEMT_u(indx), p%BEMT, m%BEMT_y, m%BEMT) @@ -6422,7 +7591,7 @@ SUBROUTINE AD_JacobianPInput(Vars, iRotor, t, u_AD, p_AD, x_AD, xd_AD, z_AD, Oth ! Initialize x_init so that we get accurrate values for first step ! changes values only if states haven't been initialized if ((p_AD%Wake_Mod /= WakeMod_FVW) .and. (.not. OtherState%BEMT%nodesInitialized)) then - call SetInputs(t, p, p_AD, u, RotInflow, m, indx, errStat2, errMsg2); if (Failed()) return + call SetInputs(t, p, p_AD, u, RotInflow, GSInflow, m, indx, errStat2, errMsg2); if (Failed()) return call BEMT_InitStates(t, m%BEMT_u(indx), p%BEMT, m%x_init%BEMT, xd%BEMT, z%BEMT, & m%OtherState_init%BEMT, m%BEMT, p_AD%AFI, ErrStat2, ErrMsg2); if (Failed()) return end if @@ -6431,6 +7600,10 @@ SUBROUTINE AD_JacobianPInput(Vars, iRotor, t, u_AD, p_AD, x_AD, xd_AD, z_AD, Oth call AD_CopyRotInputType(u, u_perturb, MESH_UPDATECOPY, ErrStat2, ErrMsg2); if (Failed()) return call AD_VarsPackInput(Vars, u, m%Jac%u) + ! Copy rotor and GS inflow types for perturbation (shared by the dYdu and dXdu loops below) + call AD_CopyRotInflowType(RotInflow, RotInflow_perturb, MESH_NEWCOPY, ErrStat2, ErrMsg2); if (Failed()) return + call AD_CopyElemInflowType(GSInflow, GSInflow_perturb, MESH_NEWCOPY, ErrStat2, ErrMsg2); if (Failed()) return + ! Calculate the partial derivative of the output functions (Y) with respect to the inputs (u) here: if (present(dYdu)) then @@ -6439,9 +7612,6 @@ SUBROUTINE AD_JacobianPInput(Vars, iRotor, t, u_AD, p_AD, x_AD, xd_AD, z_AD, Oth call AllocAry(dYdu, Vars%Ny, Vars%Nu, 'dYdu', ErrStat2, ErrMsg2); if (Failed()) return end if - ! Copy rotor inflow type for perturbation - call AD_CopyRotInflowType(RotInflow, RotInflow_perturb, MESH_NEWCOPY, ErrStat2, ErrMsg2); if (Failed()) return - ! Loop through input variables do i = 1, size(Vars%u) @@ -6456,9 +7626,10 @@ SUBROUTINE AD_JacobianPInput(Vars, iRotor, t, u_AD, p_AD, x_AD, xd_AD, z_AD, Oth if (associated(FF_ptr, FF_perturb)) call PerturbFlowField(Vars%u(i), p_AD%FlowField, 1, FF_ptr) StartNode = 1 call AD_CalcWind_Rotor(t, u_perturb, FF_ptr, p, p_AD, m_AD, RotInflow_perturb, StartNode, ErrStat2, ErrMsg2); if (Failed()) return - call SetInputs(t, p, p_AD, u_perturb, RotInflow_perturb, m, indx, ErrStat2, ErrMsg2); if (Failed()) return + call AD_CalcWind_GS(t, u_perturb, FF_ptr, p_AD%GS, p_AD, m_AD, GSInflow_perturb, StartNode, ErrStat2, ErrMsg2); if (Failed()) return + call SetInputs(t, p, p_AD, u_perturb, RotInflow_perturb, GSInflow_perturb, m, indx, ErrStat2, ErrMsg2); if (Failed()) return call UpdatePhi(m%BEMT_u(indx), p%BEMT, m%z_lin%BEMT%phi, p_AD%AFI, m%BEMT, m%OtherState_jac%BEMT%ValidPhi, ErrStat2, ErrMsg2); if (Failed()) return - call RotCalcOutput(t, u_perturb, RotInflow_perturb, p, p_AD, m%x_init, xd, m%z_lin, m%OtherState_jac, y_lin, m, m_AD, iRotor, ErrStat2, ErrMsg2); if (Failed()) return + call RotCalcOutput(t, u_perturb, RotInflow_perturb, GSInflow_perturb, p, p_AD, m%x_init, xd, m%z_lin, m%OtherState_jac, y_lin, m, m_AD, iRotor, ErrStat2, ErrMsg2); if (Failed()) return if (p_AD%Wake_Mod == WakeMod_FVW) then call SetInputsForFVW(p_AD, m_AD%u_perturb, 1, m_AD, ErrStat2, ErrMsg2); if(Failed()) return call FVW_CalcOutput(t, m_AD%FVW_u(1), p_AD%FVW, x_AD%FVW, xd_AD%FVW, z_AD%FVW, OtherState_AD%FVW, m_AD%FVW_y, m_AD%FVW, ErrStat2, ErrMsg2); if(Failed()) return @@ -6474,9 +7645,10 @@ SUBROUTINE AD_JacobianPInput(Vars, iRotor, t, u_AD, p_AD, x_AD, xd_AD, z_AD, Oth if (associated(FF_ptr, FF_perturb)) call PerturbFlowField(Vars%u(i), p_AD%FlowField, -1, FF_ptr) StartNode = 1 call AD_CalcWind_Rotor(t, u_perturb, FF_ptr, p, p_AD, m_AD, RotInflow_perturb, StartNode, ErrStat2, ErrMsg2); if (Failed()) return - call SetInputs(t, p, p_AD, u_perturb, RotInflow_perturb, m, indx, ErrStat2, ErrMsg2); if (Failed()) return + call AD_CalcWind_GS(t, u_perturb, FF_ptr, p_AD%GS, p_AD, m_AD, GSInflow_perturb, StartNode, ErrStat2, ErrMsg2); if (Failed()) return + call SetInputs(t, p, p_AD, u_perturb, RotInflow_perturb, GSInflow_perturb, m, indx, ErrStat2, ErrMsg2); if (Failed()) return call UpdatePhi(m%BEMT_u(indx), p%BEMT, m%z_lin%BEMT%phi, p_AD%AFI, m%BEMT, m%OtherState_jac%BEMT%ValidPhi, ErrStat2, ErrMsg2); if (Failed()) return - call RotCalcOutput(t, u_perturb, RotInflow_perturb, p, p_AD, m%x_init, xd, m%z_lin, m%OtherState_jac, y_lin, m, m_AD, iRotor, ErrStat2, ErrMsg2); if (Failed()) return + call RotCalcOutput(t, u_perturb, RotInflow_perturb, GSInflow_perturb, p, p_AD, m%x_init, xd, m%z_lin, m%OtherState_jac, y_lin, m, m_AD, iRotor, ErrStat2, ErrMsg2); if (Failed()) return if (p_AD%Wake_Mod == WakeMod_FVW) then call SetInputsForFVW(p_AD, m_AD%u_perturb, 1, m_AD, ErrStat2, ErrMsg2); if(Failed()) return call FVW_CalcOutput(t, m_AD%FVW_u(1), p_AD%FVW, x_AD%FVW, xd_AD%FVW, z_AD%FVW, OtherState_AD%FVW, m_AD%FVW_y, m_AD%FVW, ErrStat2, ErrMsg2); if(Failed()) return @@ -6520,7 +7692,8 @@ SUBROUTINE AD_JacobianPInput(Vars, iRotor, t, u_AD, p_AD, x_AD, xd_AD, z_AD, Oth if (associated(FF_ptr, FF_perturb)) call PerturbFlowField(Vars%u(i), p_AD%FlowField, 1, FF_ptr) StartNode = 1 call AD_CalcWind_Rotor(t, u_perturb, FF_ptr, p, p_AD, m_AD, RotInflow_perturb, StartNode, ErrStat2, ErrMsg2); if (Failed()) return - call RotCalcContStateDeriv(t, u_perturb, RotInflow_perturb, p, p_AD, m%x_init, xd, z, m%OtherState_init, m, m%dxdt_lin, ErrStat2, ErrMsg2) ; if (Failed()) return + call AD_CalcWind_GS(t, u_perturb, FF_ptr, p_AD%GS, p_AD, m_AD, GSInflow_perturb, StartNode, ErrStat2, ErrMsg2); if (Failed()) return + call RotCalcContStateDeriv(t, u_perturb, RotInflow_perturb, GSInflow_perturb, p, p_AD, m%x_init, xd, z, m%OtherState_init, m, m%dxdt_lin, ErrStat2, ErrMsg2) ; if (Failed()) return call AD_VarsPackContState(Vars, m%dxdt_lin, m%Jac%x_pos) ! Calculate negative perturbation @@ -6529,7 +7702,8 @@ SUBROUTINE AD_JacobianPInput(Vars, iRotor, t, u_AD, p_AD, x_AD, xd_AD, z_AD, Oth if (associated(FF_ptr, FF_perturb)) call PerturbFlowField(Vars%u(i), p_AD%FlowField, -1, FF_ptr) StartNode = 1 call AD_CalcWind_Rotor(t, u_perturb, FF_ptr, p, p_AD, m_AD, RotInflow_perturb, StartNode, ErrStat2, ErrMsg2); if (Failed()) return - call RotCalcContStateDeriv(t, u_perturb, RotInflow_perturb, p, p_AD, m%x_init, xd, z, m%OtherState_init, m, m%dxdt_lin, ErrStat2, ErrMsg2) ; if (Failed()) return + call AD_CalcWind_GS(t, u_perturb, FF_ptr, p_AD%GS, p_AD, m_AD, GSInflow_perturb, StartNode, ErrStat2, ErrMsg2); if (Failed()) return + call RotCalcContStateDeriv(t, u_perturb, RotInflow_perturb, GSInflow_perturb, p, p_AD, m%x_init, xd, z, m%OtherState_init, m, m%dxdt_lin, ErrStat2, ErrMsg2) ; if (Failed()) return call AD_VarsPackContState(Vars, m%dxdt_lin, m%Jac%x_neg) ! Calculate column index @@ -6580,6 +7754,8 @@ logical function Failed() end function subroutine cleanup() + call AD_DestroyRotInflowType(RotInflow_perturb, ErrStat2, ErrMsg2) + call AD_DestroyElemInflowType(GSInflow_perturb, ErrStat2, ErrMsg2) m_AD%rotors(iRotor)%BEMT%UseFrozenWake = .false. end subroutine cleanup end subroutine @@ -6622,14 +7798,16 @@ SUBROUTINE AD_JacobianPContState(Vars, iRotor, t, u, p, x, xd, z, OtherState, y, StartNode = 1 call AD_CalcWind_Rotor(t, u%rotors(iRotor), p%FlowField, p%rotors(iRotor), p, m, m%Inflow(1)%RotInflow(iRotor), StartNode, ErrStat, ErrMsg) if (ErrStat >= AbortErrLev) return - call RotJacobianPContState(Vars, iRotor, t, u%rotors(iRotor), m%Inflow(1)%RotInflow(iRotor), p%rotors(iRotor), p, x%rotors(iRotor), xd%rotors(iRotor), z%rotors(iRotor), OtherState%rotors(iRotor), y%rotors(iRotor), m%rotors(iRotor), m, ErrStat, ErrMsg, dYdx, dXdx, dXddx, dZdx) + call AD_CalcWind_GS(t, u%rotors(1), p%FlowField, p%GS, p, m, m%Inflow(1)%GSInflow, StartNode, ErrStat, ErrMsg) + if (ErrStat >= AbortErrLev) return + call RotJacobianPContState(Vars, iRotor, t, u%rotors(iRotor), m%Inflow(1)%RotInflow(iRotor), m%Inflow(1)%GSInflow, p%rotors(iRotor), p, x%rotors(iRotor), xd%rotors(iRotor), z%rotors(iRotor), OtherState%rotors(iRotor), y%rotors(iRotor), m%rotors(iRotor), m, ErrStat, ErrMsg, dYdx, dXdx, dXddx, dZdx) END SUBROUTINE AD_JacobianPContState !---------------------------------------------------------------------------------------------------------------------------------- !> Routine to compute the Jacobians of the output (Y), continuous- (X), discrete- (Xd), and constraint-state (Z) functions !! with respect to the continuous states (x). The partial derivatives dY/dx, dX/dx, dXd/dx, and dZ/dx are returned. -SUBROUTINE RotJacobianPContState(Vars, iRotor, t, u, RotInflow, p, p_AD, x, xd, z, OtherState, y, m, m_AD, ErrStat, ErrMsg, dYdx, dXdx, dXddx, dZdx) +SUBROUTINE RotJacobianPContState(Vars, iRotor, t, u, RotInflow, GSInflow, p, p_AD, x, xd, z, OtherState, y, m, m_AD, ErrStat, ErrMsg, dYdx, dXdx, dXddx, dZdx) !.................................................................................................................................. TYPE(ModVarsType), INTENT(IN ) :: Vars !< Module variables for packing arrays @@ -6637,6 +7815,7 @@ SUBROUTINE RotJacobianPContState(Vars, iRotor, t, u, RotInflow, p, p_AD, x, xd, REAL(DbKi), INTENT(IN ) :: t !< Time in seconds at operating point TYPE(RotInputType), INTENT(IN ) :: u !< Inputs at operating point (may change to inout if a mesh copy is required) TYPE(RotInflowType), INTENT(IN ) :: RotInflow !< Rotor inflow + TYPE(ElemInflowType), INTENT(IN ) :: GSInflow !< General support structure inflow TYPE(RotParameterType), INTENT(IN ) :: p !< Parameters TYPE(AD_ParameterType), INTENT(IN ) :: p_AD !< Parameters TYPE(RotContinuousStateType), INTENT(IN ) :: x !< Continuous states at operating point @@ -6667,7 +7846,7 @@ SUBROUTINE RotJacobianPContState(Vars, iRotor, t, u, RotInflow, p, p_AD, x, xd, ! Get OP values here (i.e., set inputs for BEMT): if (p%DBEMT_Mod == DBEMT_frozen) then - call SetInputs(t, p, p_AD, u, RotInflow, m, indx, errStat2, errMsg2); if (Failed()) return + call SetInputs(t, p, p_AD, u, RotInflow, GSInflow, m, indx, errStat2, errMsg2); if (Failed()) return ! compare m%BEMT_y arguments with call to BEMT_CalcOutput call computeFrozenWake(m%BEMT_u(indx), p%BEMT, m%BEMT_y, m%BEMT) @@ -6681,7 +7860,7 @@ SUBROUTINE RotJacobianPContState(Vars, iRotor, t, u, RotInflow, p, p_AD, x, xd, ! Initialize x_init so that we get accurrate values for first step ! changes values only if states haven't been initialized if (.not. OtherState%BEMT%nodesInitialized) then - call SetInputs(t, p, p_AD, u, RotInflow, m, indx, errStat2, errMsg2); if (Failed()) return + call SetInputs(t, p, p_AD, u, RotInflow, GSInflow, m, indx, errStat2, errMsg2); if (Failed()) return call BEMT_InitStates(t, m%BEMT_u(indx), p%BEMT, m%x_init%BEMT, xd%BEMT, z%BEMT, & m%OtherState_init%BEMT, m%BEMT, p_AD%AFI, ErrStat2, ErrMsg2); if (Failed()) return end if @@ -6707,13 +7886,13 @@ SUBROUTINE RotJacobianPContState(Vars, iRotor, t, u, RotInflow, p, p_AD, x, xd, ! Calculate positive perturbation call MV_Perturb(Vars%x(i), j, 1, m%Jac%x, m%Jac%x_perturb) call AD_VarsUnpackContState(Vars, m%Jac%x_perturb, m%x_perturb) - call RotCalcOutput(t, u, RotInflow, p, p_AD, m%x_perturb, xd, z, m%OtherState_init, m%y_lin, m, m_AD, iRotor, ErrStat2, ErrMsg2) ; if (Failed()) return + call RotCalcOutput(t, u, RotInflow, GSInflow, p, p_AD, m%x_perturb, xd, z, m%OtherState_init, m%y_lin, m, m_AD, iRotor, ErrStat2, ErrMsg2) ; if (Failed()) return call AD_VarsPackOutput(Vars, m%y_lin, m%Jac%y_pos) ! Calculate negative perturbation call MV_Perturb(Vars%x(i), j, -1, m%Jac%x, m%Jac%x_perturb) call AD_VarsUnpackContState(Vars, m%Jac%x_perturb, m%x_perturb) - call RotCalcOutput(t, u, RotInflow, p, p_AD, m%x_perturb, xd, z, m%OtherState_init, m%y_lin, m, m_AD, iRotor, ErrStat2, ErrMsg2) ; if (Failed()) return + call RotCalcOutput(t, u, RotInflow, GSInflow, p, p_AD, m%x_perturb, xd, z, m%OtherState_init, m%y_lin, m, m_AD, iRotor, ErrStat2, ErrMsg2) ; if (Failed()) return call AD_VarsPackOutput(Vars, m%y_lin, m%Jac%y_neg) ! Calculate column index @@ -6743,13 +7922,13 @@ SUBROUTINE RotJacobianPContState(Vars, iRotor, t, u, RotInflow, p, p_AD, x, xd, ! Calculate positive perturbation call MV_Perturb(Vars%x(i), j, 1, m%Jac%x, m%Jac%x_perturb) call AD_VarsUnpackContState(Vars, m%Jac%x_perturb, m%x_perturb) - call RotCalcContStateDeriv(t, u, RotInflow, p, p_AD, m%x_perturb, xd, z, m%OtherState_init, m, m%dxdt_lin, ErrStat2, ErrMsg2); if (Failed()) return + call RotCalcContStateDeriv(t, u, RotInflow, GSInflow, p, p_AD, m%x_perturb, xd, z, m%OtherState_init, m, m%dxdt_lin, ErrStat2, ErrMsg2); if (Failed()) return call AD_VarsPackContState(Vars, m%dxdt_lin, m%Jac%x_pos) ! Calculate negative perturbation call MV_Perturb(Vars%x(i), j, -1, m%Jac%x, m%Jac%x_perturb) call AD_VarsUnpackContState(Vars, m%Jac%x_perturb, m%x_perturb) - call RotCalcContStateDeriv(t, u, RotInflow, p, p_AD, m%x_perturb, xd, z, m%OtherState_init, m, m%dxdt_lin, ErrStat2, ErrMsg2); if (Failed()) return + call RotCalcContStateDeriv(t, u, RotInflow, GSInflow, p, p_AD, m%x_perturb, xd, z, m%OtherState_init, m, m%dxdt_lin, ErrStat2, ErrMsg2); if (Failed()) return call AD_VarsPackContState(Vars, m%dxdt_lin, m%Jac%x_neg) ! Calculate column index @@ -6855,7 +8034,9 @@ SUBROUTINE AD_JacobianPConstrState(Vars, t, u, p, x, xd, z, OtherState, y, m, Er StartNode = 1 call AD_CalcWind_Rotor(t, u%rotors(iR), p%FlowField, p%rotors(iR), p, m, m%Inflow(1)%RotInflow(iR), StartNode, ErrStat, ErrMsg) if (ErrStat >= AbortErrLev) return - call RotJacobianPConstrState(t, u%rotors(iR), m%Inflow(1)%RotInflow(iR), p%rotors(iR), p, x%rotors(iR), xd%rotors(iR), z%rotors(iR), OtherState%rotors(iR), y%rotors(iR), m%rotors(iR), m, iR, errStat, errMsg, dYdz, dXdz, dXddz, dZdz) + call AD_CalcWind_GS(t, u%rotors(1), p%FlowField, p%GS, p, m, m%Inflow(1)%GSInflow, StartNode, ErrStat, ErrMsg) + if (ErrStat >= AbortErrLev) return + call RotJacobianPConstrState(t, u%rotors(iR), m%Inflow(1)%RotInflow(iR), m%Inflow(1)%GSInflow, p%rotors(iR), p, x%rotors(iR), xd%rotors(iR), z%rotors(iR), OtherState%rotors(iR), y%rotors(iR), m%rotors(iR), m, iR, errStat, errMsg, dYdz, dXdz, dXddz, dZdz) END SUBROUTINE AD_JacobianPConstrState @@ -6863,10 +8044,11 @@ END SUBROUTINE AD_JacobianPConstrState !---------------------------------------------------------------------------------------------------------------------------------- !> Routine to compute the Jacobians of the output (Y), continuous- (X), discrete- (Xd), and constraint-state (Z) functions !! with respect to the constraint states (z). The partial derivatives dY/dz, dX/dz, dXd/dz, and dZ/dz are returned. -SUBROUTINE RotJacobianPConstrState( t, u, RotInflow, p, p_AD, x, xd, z, OtherState, y, m, m_AD, iRot, ErrStat, ErrMsg, dYdz, dXdz, dXddz, dZdz ) +SUBROUTINE RotJacobianPConstrState( t, u, RotInflow, GSInflow, p, p_AD, x, xd, z, OtherState, y, m, m_AD, iRot, ErrStat, ErrMsg, dYdz, dXdz, dXddz, dZdz ) REAL(DbKi), INTENT(IN ) :: t !< Time in seconds at operating point TYPE(RotInputType), INTENT(IN ) :: u !< Inputs at operating point (may change to inout if a mesh copy is required) - TYPE(RotInflowType), INTENT(IN ) :: RotInflow !< Inflow on rotor + TYPE(RotInflowType), INTENT(IN ) :: RotInflow !< Inflow on rotor + TYPE(ElemInflowType), INTENT(IN ) :: GSInflow !< Inflow on general support structure TYPE(RotParameterType), INTENT(IN ) :: p !< Parameters TYPE(AD_ParameterType), INTENT(IN ) :: p_AD !< Parameters TYPE(RotContinuousStateType), INTENT(IN ) :: x !< Continuous states at operating point @@ -6909,7 +8091,7 @@ SUBROUTINE RotJacobianPConstrState( t, u, RotInflow, p, p_AD, x, xd, z, OtherSta ! get OP values here: !call AD_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, ErrStat2, ErrMsg2 ) ! (bjj: is this necessary? if not, still need to get BEMT inputs) - call SetInputs(t, p, p_AD, u, RotInflow, m, indx, errStat2, errMsg2); if (Failed()) return; + call SetInputs(t, p, p_AD, u, RotInflow, GSInflow, m, indx, errStat2, errMsg2); if (Failed()) return; call BEMT_CopyInput( m%BEMT_u(indx), m%BEMT_u(op_indx), MESH_UPDATECOPY, ErrStat2, ErrMsg2); if (Failed()) return; ! copy the BEMT OP inputs to a temporary location that won't be overwritten @@ -6954,13 +8136,13 @@ SUBROUTINE RotJacobianPConstrState( t, u, RotInflow, p, p_AD, x, xd, z, OtherSta z_perturb%BEMT%phi(j,k) = z%BEMT%phi(j,k) + delta_p ! compute y at z_op + delta_p z - call RotCalcOutput( t, u, RotInflow, p, p_AD, x, xd, z_perturb, OtherState, y_p, m, m_AD, iRot, ErrStat2, ErrMsg2 ) ; if (Failed()) return; + call RotCalcOutput( t, u, RotInflow, GSInflow, p, p_AD, x, xd, z_perturb, OtherState, y_p, m, m_AD, iRot, ErrStat2, ErrMsg2 ) ; if (Failed()) return; ! get z_op - delta_m z z_perturb%BEMT%phi(j,k) = z%BEMT%phi(j,k) - delta_m ! compute y at z_op - delta_m z - call RotCalcOutput( t, u, RotInflow, p, p_AD, x, xd, z_perturb, OtherState, y_m, m, m_AD, iRot, ErrStat2, ErrMsg2 ) ; if (Failed()) return; + call RotCalcOutput( t, u, RotInflow, GSInflow, p, p_AD, x, xd, z_perturb, OtherState, y_m, m, m_AD, iRot, ErrStat2, ErrMsg2 ) ; if (Failed()) return; ! get central difference: call Compute_dY( p, p_AD, y_p, y_m, delta_p, delta_m, dYdz(:,i) ) @@ -7009,13 +8191,13 @@ SUBROUTINE RotJacobianPConstrState( t, u, RotInflow, p, p_AD, x, xd, z, OtherSta z_perturb%BEMT%phi(j,k) = z%BEMT%phi(j,k) + delta_p ! compute z_p at z_op + delta_p z - call RotCalcConstrStateResidual( t, u, RotInflow, p, p_AD, x, xd, z_perturb, OtherState, m, z_p, ErrStat2, ErrMsg2 ) ; if (Failed()) return; + call RotCalcConstrStateResidual( t, u, RotInflow, GSInflow, p, p_AD, x, xd, z_perturb, OtherState, m, z_p, ErrStat2, ErrMsg2 ) ; if (Failed()) return; ! get z_op - delta_m z z_perturb%BEMT%phi(j,k) = z%BEMT%phi(j,k) - delta_m ! compute z_m at u_op - delta_m u - call RotCalcConstrStateResidual( t, u, RotInflow, p, p_AD, x, xd, z_perturb, OtherState, m, z_m, ErrStat2, ErrMsg2 ) ; if (Failed()) return; + call RotCalcConstrStateResidual( t, u, RotInflow, GSInflow, p, p_AD, x, xd, z_perturb, OtherState, m, z_m, ErrStat2, ErrMsg2 ) ; if (Failed()) return; ! get central difference: do k2=1,p%NumBlades ! size(z%BEMT%Phi,2) diff --git a/modules/aerodyn/src/AeroDyn_IO.f90 b/modules/aerodyn/src/AeroDyn_IO.f90 index 5f919bb211..8493b31a57 100644 --- a/modules/aerodyn/src/AeroDyn_IO.f90 +++ b/modules/aerodyn/src/AeroDyn_IO.f90 @@ -296,17 +296,32 @@ subroutine Calc_WriteOutput_AD() end do ! blades end if - ! blade node tower clearance (requires tower influence calculation): - if (p%TwrPotent /= TwrPotent_none .or. p%TwrShadow /= TwrShadow_none) then + ! blade node clearance to tower and/or generalized support (requires influence calculation): + ! each clearance array is allocated only when its model is active; report the smallest. + if (allocated(m%TwrClrnc)) then do k=1,min(p%numBlades,AD_MaxBl_Out) do beta=1,p%NBlOuts j=p%BlOutNd(beta) m%AllOuts( BNClrnc( beta,k) ) = m%TwrClrnc(j,k) end do end do + ! reduce with the GS clearance when it is also active + if (allocated(m%GSClrnc)) then + do k=1,min(p%numBlades,AD_MaxBl_Out) + do beta=1,p%NBlOuts + j=p%BlOutNd(beta) + m%AllOuts( BNClrnc( beta,k) ) = min( m%AllOuts( BNClrnc( beta,k) ), m%GSClrnc(j,k) ) + end do + end do + end if + else if (allocated(m%GSClrnc)) then + do k=1,min(p%numBlades,AD_MaxBl_Out) + do beta=1,p%NBlOuts + j=p%BlOutNd(beta) + m%AllOuts( BNClrnc( beta,k) ) = m%GSClrnc(j,k) + end do + end do end if - - m%AllOuts( RtSpeed ) = omega*RPS2RPM @@ -702,8 +717,9 @@ SUBROUTINE ParsePrimaryFileInfo( PriPath, InitInp, InputFile, RootName, NumBlade character(ErrMsgLen) :: ErrMsg2 !< Temporary Error message character(ErrMsgLen) :: ErrMsg_NoAllBldNdOuts integer(IntKi) :: CurLine !< current entry in FileInfo_In%Lines array - real(ReKi) :: TmpRe7(7) !< temporary 8 number array for reading values in + real(ReKi) :: TmpRe10(10) !< temporary 10 number array for reading values in logical :: TwrAeroLogical !< convert TwrAero from logical (input file) to integer (new) + logical :: GSAeroLogical !< convert GSAero from logical (input file) to integer (new) character(1024) :: sDummy !< temporary string character(1024) :: tmpOutStr !< temporary string for writing to screen logical :: wakeModProvided, frozenWakeProvided, skewModProvided, AFAeroModProvided, UAModProvided, isLegalComment, firstWarn !< Temporary for legacy purposes @@ -798,7 +814,22 @@ SUBROUTINE ParsePrimaryFileInfo( PriPath, InitInp, InputFile, RootName, NumBlade else InputFileData%TwrAero = TwrAero_None end if - + + ! GSPotent - Type of general support influence on wind based on potential flow around the member (switch) + call ParseVar( FileInfo_In, CurLine, "GSPotent", InputFileData%GS%GSPotent, ErrStat2, ErrMsg2, UnEc ) + if (Failed()) return + ! GSShadow - Type of general support influence on wind based on downstream tower shadow {0=none, 1=Powles model, 2=Eames model} + call ParseVar( FileInfo_In, CurLine, "GSShadow", InputFileData%GS%GSShadow, ErrStat2, ErrMsg2, UnEc ) + if (Failed()) return + + ! GSAero - Calculate multi-member generalized tower aerodynamic loads? (flag) + call ParseVar( FileInfo_In, CurLine, "GSAero", GSAeroLogical, ErrStat2, ErrMsg2, UnEc ) + if (Failed()) return + if (GSAeroLogical) then + InputFileData%GS%GSAero = GSAero_NoVIV + else + InputFileData%GS%GSAero = GSAero_None + end if ! FrozenWake - Assume frozen wake during linearization? (flag) [used only when WakeMod=1 and when linearizing] call ParseVar( FileInfo_In, CurLine, "FrozenWake", FrozenWake_Old, ErrStat2, ErrMsg2, UnEc ) frozenWakeProvided = legacyInputPresent('FrozenWake', Curline, ErrStat2, ErrMsg2, 'DBEMTMod=-1 (FrozenWake=True) or DBEMTMod>-1 (FrozenWake=False)') @@ -1111,18 +1142,84 @@ SUBROUTINE ParsePrimaryFileInfo( PriPath, InitInp, InputFile, RootName, NumBlade CALL AllocAry( InputFileData%rotors(iR)%TwrCa, InputFileData%rotors(iR)%NumTwrNds, 'TwrCa', ErrStat2, ErrMsg2) if (Failed()) return do I=1,InputFileData%rotors(iR)%NumTwrNds - call ParseAry ( FileInfo_In, CurLine, 'Properties for tower node '//trim( Int2LStr( I ) )//'.', TmpRe7, 7, ErrStat2, ErrMsg2, UnEc ) + call ParseAry ( FileInfo_In, CurLine, 'Properties for tower node '//trim( Int2LStr( I ) )//'.', TmpRe10, 7, ErrStat2, ErrMsg2, UnEc ) if (Failed()) return; - InputFileData%rotors(iR)%TwrElev(I) = TmpRe7( 1) - InputFileData%rotors(iR)%TwrDiam(I) = TmpRe7( 2) - InputFileData%rotors(iR)%TwrCd(I) = TmpRe7( 3) - InputFileData%rotors(iR)%TwrTI(I) = TmpRe7( 4) - InputFileData%rotors(iR)%TwrCb(I) = TmpRe7( 5) - InputFileData%rotors(iR)%TwrCp(I) = TmpRe7( 6) - InputFileData%rotors(iR)%TwrCa(I) = TmpRe7( 7) + InputFileData%rotors(iR)%TwrElev(I) = TmpRe10( 1) + InputFileData%rotors(iR)%TwrDiam(I) = TmpRe10( 2) + InputFileData%rotors(iR)%TwrCd(I) = TmpRe10( 3) + InputFileData%rotors(iR)%TwrTI(I) = TmpRe10( 4) + InputFileData%rotors(iR)%TwrCb(I) = TmpRe10( 5) + InputFileData%rotors(iR)%TwrCp(I) = TmpRe10( 6) + InputFileData%rotors(iR)%TwrCa(I) = TmpRe10( 7) end do enddo + !====== General support structure joints ============================================================ + if ( InputFileData%Echo ) WRITE(UnEc, '(A)') FileInfo_In%Lines(CurLine) ! Write section break to echo + CurLine = CurLine + 1 + ! NumGSJoints - Number of general support joints used in the analysis (-) + call ParseVar( FileInfo_In, CurLine, "NumGSJoints", InputFileData%GS%NJoints, ErrStat2, ErrMsg2, UnEc ) + if (Failed()) return + !GSJointID GSJointXi GSJointYi GSJointZi + if ( InputFileData%Echo ) WRITE(UnEc, '(A)') 'GS Joint Table Header: '//FileInfo_In%Lines(CurLine) ! Write section break to echo + CurLine = CurLine + 1 + !(-) (m) (m) (m) + if ( InputFileData%Echo ) WRITE(UnEc, '(A)') 'GS Joint Table Header: '//FileInfo_In%Lines(CurLine) ! Write section break to echo + CurLine = CurLine + 1 + + ! Allocate space for general support joint table + allocate( InputFileData%GS%InpJoints(InputFileData%GS%NJoints), STAT = ErrStat2) + if ( ErrStat2 /= 0 ) then + ErrStat2 = ErrID_Fatal + ErrMsg2 = 'Error allocating space for GSInpJoints array.' + if (Failed()) return; + end if + + do I=1,InputFileData%GS%NJoints + call ParseAry ( FileInfo_In, CurLine, 'General support joint '//trim( Int2LStr( I ) )//'.', TmpRe10, 4, ErrStat2, ErrMsg2, UnEc ) + if (Failed()) return; + InputFileData%GS%InpJoints(I)%JointID = NINT(TmpRe10(1)) + InputFileData%GS%InpJoints(I)%position(1) = TmpRe10(2) + InputFileData%GS%InpJoints(I)%position(2) = TmpRe10(3) + InputFileData%GS%InpJoints(I)%position(3) = TmpRe10(4) + end do + + !====== General support structure members ============================================================ + if ( InputFileData%Echo ) WRITE(UnEc, '(A)') FileInfo_In%Lines(CurLine) ! Write section break to echo + CurLine = CurLine + 1 + ! NumGSMembers - Number of general support members used in the analysis (-) + call ParseVar( FileInfo_In, CurLine, "NumGSMembers", InputFileData%GS%NMembers, ErrStat2, ErrMsg2, UnEc ) + if (Failed()) return + !GSMemberID GSMJointID1 GSMJointID2 GSMDia1 GSMDia2 GSMCd1 GSMCd2 GSMTI1 GSMTI2 GSMDiv + if ( InputFileData%Echo ) WRITE(UnEc, '(A)') 'GS Member Table Header: '//FileInfo_In%Lines(CurLine) ! Write section break to echo + CurLine = CurLine + 1 + !(-) (-) (-) (m) (m) (-) (-) (-) (-) (-) + if ( InputFileData%Echo ) WRITE(UnEc, '(A)') 'GS Member Table Header: '//FileInfo_In%Lines(CurLine) ! Write section break to echo + CurLine = CurLine + 1 + + ! Allocate space for general support joint table + allocate( InputFileData%GS%InpMembers(InputFileData%GS%NMembers), STAT = ErrStat2) + if ( ErrStat2 /= 0 ) then + ErrStat2 = ErrID_Fatal + ErrMsg2 = 'Error allocating space for GSInpMembers array.' + if (Failed()) return; + end if + + do I=1,InputFileData%GS%NMembers + call ParseAry ( FileInfo_In, CurLine, 'General support member '//trim( Int2LStr( I ) )//'.', TmpRe10, 10, ErrStat2, ErrMsg2, UnEc ) + if (Failed()) return; + InputFileData%GS%InpMembers(I)%MemberID = NINT(TmpRe10( 1)) + InputFileData%GS%InpMembers(I)%MJointID1 = NINT(TmpRe10( 2)) + InputFileData%GS%InpMembers(I)%MJointID2 = NINT(TmpRe10( 3)) + InputFileData%GS%InpMembers(I)%MDiam1 = TmpRe10( 4) + InputFileData%GS%InpMembers(I)%MDiam2 = TmpRe10( 5) + InputFileData%GS%InpMembers(I)%MCd1 = TmpRe10( 6) + InputFileData%GS%InpMembers(I)%MCd2 = TmpRe10( 7) + InputFileData%GS%InpMembers(I)%MTI1 = TmpRe10( 8) + InputFileData%GS%InpMembers(I)%MTI2 = TmpRe10( 9) + InputFileData%GS%InpMembers(I)%MDivSize = TmpRe10(10) + end do + !====== Outputs ==================================================================================== if ( InputFileData%Echo ) WRITE(UnEc, '(A)') FileInfo_In%Lines(CurLine) ! Write section break to echo CurLine = CurLine + 1 @@ -2104,6 +2201,108 @@ END SUBROUTINE AD_PrintSum !---------------------------------------------------------------------------------------------------------------------------------- +!---------------------------------------------------------------------------------------------------------------------------------- +SUBROUTINE AD_PrintSum_GS( p, p_AD, u, ErrStat, ErrMsg ) +! This routine generates the summary file, which contains a summary of input file options. + ! passed variables + TYPE(GSParameterType), INTENT(IN) :: p ! Parameters + TYPE(AD_ParameterType), INTENT(IN) :: p_AD ! Parameters + TYPE(RotInputType), INTENT(IN) :: u ! inputs + + INTEGER(IntKi), INTENT(OUT) :: ErrStat + CHARACTER(*), INTENT(OUT) :: ErrMsg + + ! Local variables. + + INTEGER(IntKi) :: I, J ! Generic loop counter + INTEGER(IntKi) :: JointID, NodeIndx + INTEGER(IntKi) :: UnSu ! I/O unit number for the summary output file + CHARACTER(100) :: Msg ! temporary string for writing appropriate text to summary file + + ErrStat = ErrID_None + ErrMsg = "" + + if (p%NMembers<=0_IntKi) return + + ! Open the summary file and give it a heading. + + !$OMP critical(fileopen_critical) + CALL GetNewUnit( UnSu, ErrStat, ErrMsg ) + CALL OpenFOutFile ( UnSu, TRIM( p_AD%RootName )//'.GS.sum', ErrStat, ErrMsg ) + !$OMP end critical(fileopen_critical) + IF ( ErrStat >= AbortErrLev ) RETURN + + ! Heading: + WRITE (UnSu,'(/,A)') 'This summary information was generated by '//TRIM( GetNVD(AD_Ver) )// & + ' on '//CurDate()//' at '//CurTime()//'.' + + WRITE (UnSu,'(/,/,A)') '====== General support structure ==================================================================' + + select case (p%GSPotent) + case (GSPotent_none) + Msg = 'none' + case (GSPotent_baseline) + Msg = 'baseline model' + case (GSPotent_bak) + Msg = "Bak correction" + case default + Msg = 'unknown' + end select + WRITE (UnSu,'(/,A)') 'Potential-flow influence model: '//trim(num2lstr(p%GSPotent))//' ('//trim(Msg)//')' + + select case (p%GSShadow) + case (GSShadow_none) + Msg = 'none' + case (GSShadow_Powles) + Msg = 'Powles model' + case (GSShadow_Eames) + Msg = "Eames model" + case default + Msg = 'unknown' + end select + WRITE (UnSu,'(A)') 'Downstream shadow/wake model: '//trim(num2lstr(p%GSShadow))//' ('//trim(Msg)//')' + + select case (p%GSAero) + case (GSAero_none) + Msg = 'none' + case (GSAero_NoVIV) + Msg = 'drag only' + case default + Msg = 'unknown' + end select + WRITE (UnSu,'(A)') 'Flow-induced load model: '//trim(num2lstr(p%GSAero))//' ('//trim(Msg)//')' + + WRITE (UnSu,'(/,A)') 'Number of joints: '//trim(num2lstr(p%NJoints)) + WRITE (UnSu,'(A)') 'Number of members: '//trim(num2lstr(p%NMembers)) + WRITE (UnSu,'(A)') 'Number of nodes: '//trim(num2lstr(u%GSMotion%Nnodes))//' (joints + member interior nodes)' + + WRITE (UnSu,'(/,A)') 'In the table below, a member interior node will have joint ID -1.' + + ! MemberID, JointID, NodeIndx, NodeX, NodeY, NodeZ, R, Cd, TI + WRITE (UnSu,'(/,9(1x,A14))') 'MemberID', 'JointID', 'NodeIndx', 'Node_xi', 'Node_yi', 'Node_zi', 'Node_R', 'Node_Cd', 'Node_TI' + WRITE (UnSu,'(9(1x,A14))') ' (-)', ' (-)', ' (-)', ' (m)', ' (m)', ' (m)', ' (m)', ' (-)', ' (-)' + DO I=1,p%NMembers + associate(mem=>p%Members(I)) + DO J=1,p%Members(I)%NElements+1 + NodeIndx = mem%NodeIndx(J) + IF (NodeIndx<=p%NJoints) THEN + JointID = p%Joints(NodeIndx)%JointID + ELSE + JointID = -1 + END IF + WRITE(UnSu,'(3(I15),6(F15.6))') & + mem%MemberID, JointID, NodeIndx, & + u%GSMotion%position(1,NodeIndx), u%GSMotion%position(2,NodeIndx), u%GSMotion%position(3,NodeIndx), & + mem%R(J),mem%Cd(J),mem%TI(J) + END DO + end associate + END DO + CLOSE(UnSu) +RETURN +END SUBROUTINE AD_PrintSum_GS +!---------------------------------------------------------------------------------------------------------------------------------- + + !********************************************************************************************************************************** ! NOTE: The following lines of code were generated by a Matlab script called "Write_ChckOutLst.m" ! using the parameters listed in the "OutListParameters.xlsx" Excel file. Any changes to these diff --git a/modules/aerodyn/src/AeroDyn_Inflow_C_Binding.f90 b/modules/aerodyn/src/AeroDyn_Inflow_C_Binding.f90 index e97c8223e3..1b8f9041e1 100644 --- a/modules/aerodyn/src/AeroDyn_Inflow_C_Binding.f90 +++ b/modules/aerodyn/src/AeroDyn_Inflow_C_Binding.f90 @@ -686,6 +686,24 @@ SUBROUTINE ADI_C_Init( ADinputFilePassed, ADinputFileString_C, ADinputFileString call CheckNodes(iWT); if (Failed()) return enddo + ! Tower influence/drag and the generalized support structure (GS) are not wired up through the + ! C-binding interface: no tower/GS motion inputs are received and no tower/GS loads are returned + ! (only blade mesh motions and loads are transferred). Error out if any of these are enabled. + ! These flags/switches are set once in the shared AeroDyn input file (not rotor-specific), so + ! checking the first rotor is sufficient. Internally they are all integer switches where 0 means + ! "none"/disabled (TwrAero and GSAero are logical flags in the input file that map to 0 when False). + if ( ADI%p%AD%rotors(1)%TwrPotent /= 0 .or. ADI%p%AD%rotors(1)%TwrShadow /= 0 .or. & + ADI%p%AD%rotors(1)%TwrAero /= 0 .or. ADI%p%AD%rotors(1)%GSPotent /= 0 .or. & + ADI%p%AD%rotors(1)%GSShadow /= 0 .or. ADI%p%AD%rotors(1)%GSAero /= 0 ) then + ErrStat_F2 = ErrID_Fatal + ErrMsg_F2 = "Tower influence (TwrPotent, TwrShadow, TwrAero) and generalized support structure "// & + "(GSPotent, GSShadow, GSAero) options are not supported through the AeroDyn C-binding "// & + "interface; the tower and GS motions and loads are not transferred through the interface. "// & + "Set TwrPotent, TwrShadow, GSPotent, and GSShadow to 0, and TwrAero and GSAero to False, "// & + "in the AeroDyn input file." + if (Failed()) return + endif + !------------------------------------------------------------- ! Set the interface meshes for motion inputs and loads output diff --git a/modules/aerodyn/src/AeroDyn_Inflow_Types.f90 b/modules/aerodyn/src/AeroDyn_Inflow_Types.f90 index 470536defc..ec1974508e 100644 --- a/modules/aerodyn/src/AeroDyn_Inflow_Types.f90 +++ b/modules/aerodyn/src/AeroDyn_Inflow_Types.f90 @@ -212,17 +212,19 @@ MODULE AeroDyn_Inflow_Types integer(IntKi), public, parameter :: ADI_u_AD_rotors_BladeRootMotion = 15 ! ADI%AD%rotors(DL%i1)%BladeRootMotion(DL%i2) integer(IntKi), public, parameter :: ADI_u_AD_rotors_BladeMotion = 16 ! ADI%AD%rotors(DL%i1)%BladeMotion(DL%i2) integer(IntKi), public, parameter :: ADI_u_AD_rotors_TFinMotion = 17 ! ADI%AD%rotors(DL%i1)%TFinMotion - integer(IntKi), public, parameter :: ADI_u_AD_rotors_UserProp = 18 ! ADI%AD%rotors(DL%i1)%UserProp - integer(IntKi), public, parameter :: ADI_y_AD_rotors_NacelleLoad = 19 ! ADI%AD%rotors(DL%i1)%NacelleLoad - integer(IntKi), public, parameter :: ADI_y_AD_rotors_HubLoad = 20 ! ADI%AD%rotors(DL%i1)%HubLoad - integer(IntKi), public, parameter :: ADI_y_AD_rotors_TowerLoad = 21 ! ADI%AD%rotors(DL%i1)%TowerLoad - integer(IntKi), public, parameter :: ADI_y_AD_rotors_BladeLoad = 22 ! ADI%AD%rotors(DL%i1)%BladeLoad(DL%i2) - integer(IntKi), public, parameter :: ADI_y_AD_rotors_TFinLoad = 23 ! ADI%AD%rotors(DL%i1)%TFinLoad - integer(IntKi), public, parameter :: ADI_y_AD_rotors_WriteOutput = 24 ! ADI%AD%rotors(DL%i1)%WriteOutput - integer(IntKi), public, parameter :: ADI_y_HHVel = 25 ! ADI%HHVel - integer(IntKi), public, parameter :: ADI_y_PLExp = 26 ! ADI%PLExp - integer(IntKi), public, parameter :: ADI_y_IW_WriteOutput = 27 ! ADI%IW_WriteOutput - integer(IntKi), public, parameter :: ADI_y_WriteOutput = 28 ! ADI%WriteOutput + integer(IntKi), public, parameter :: ADI_u_AD_rotors_GSMotion = 18 ! ADI%AD%rotors(DL%i1)%GSMotion + integer(IntKi), public, parameter :: ADI_u_AD_rotors_UserProp = 19 ! ADI%AD%rotors(DL%i1)%UserProp + integer(IntKi), public, parameter :: ADI_y_AD_rotors_NacelleLoad = 20 ! ADI%AD%rotors(DL%i1)%NacelleLoad + integer(IntKi), public, parameter :: ADI_y_AD_rotors_HubLoad = 21 ! ADI%AD%rotors(DL%i1)%HubLoad + integer(IntKi), public, parameter :: ADI_y_AD_rotors_TowerLoad = 22 ! ADI%AD%rotors(DL%i1)%TowerLoad + integer(IntKi), public, parameter :: ADI_y_AD_rotors_BladeLoad = 23 ! ADI%AD%rotors(DL%i1)%BladeLoad(DL%i2) + integer(IntKi), public, parameter :: ADI_y_AD_rotors_TFinLoad = 24 ! ADI%AD%rotors(DL%i1)%TFinLoad + integer(IntKi), public, parameter :: ADI_y_AD_rotors_GSLoad = 25 ! ADI%AD%rotors(DL%i1)%GSLoad + integer(IntKi), public, parameter :: ADI_y_AD_rotors_WriteOutput = 26 ! ADI%AD%rotors(DL%i1)%WriteOutput + integer(IntKi), public, parameter :: ADI_y_HHVel = 27 ! ADI%HHVel + integer(IntKi), public, parameter :: ADI_y_PLExp = 28 ! ADI%PLExp + integer(IntKi), public, parameter :: ADI_y_IW_WriteOutput = 29 ! ADI%IW_WriteOutput + integer(IntKi), public, parameter :: ADI_y_WriteOutput = 30 ! ADI%WriteOutput contains @@ -1980,6 +1982,8 @@ function ADI_InputMeshPointer(u, DL) result(Mesh) Mesh => u%AD%rotors(DL%i1)%BladeMotion(DL%i2) case (ADI_u_AD_rotors_TFinMotion) Mesh => u%AD%rotors(DL%i1)%TFinMotion + case (ADI_u_AD_rotors_GSMotion) + Mesh => u%AD%rotors(DL%i1)%GSMotion end select end function @@ -1999,6 +2003,8 @@ function ADI_OutputMeshPointer(y, DL) result(Mesh) Mesh => y%AD%rotors(DL%i1)%BladeLoad(DL%i2) case (ADI_y_AD_rotors_TFinLoad) Mesh => y%AD%rotors(DL%i1)%TFinLoad + case (ADI_y_AD_rotors_GSLoad) + Mesh => y%AD%rotors(DL%i1)%GSLoad end select end function @@ -2191,6 +2197,8 @@ subroutine ADI_VarPackInput(V, u, ValAry) call MV_PackMesh(V, u%AD%rotors(DL%i1)%BladeMotion(DL%i2), ValAry) ! Mesh case (ADI_u_AD_rotors_TFinMotion) call MV_PackMesh(V, u%AD%rotors(DL%i1)%TFinMotion, ValAry) ! Mesh + case (ADI_u_AD_rotors_GSMotion) + call MV_PackMesh(V, u%AD%rotors(DL%i1)%GSMotion, ValAry) ! Mesh case (ADI_u_AD_rotors_UserProp) VarVals = u%AD%rotors(DL%i1)%UserProp(V%iLB:V%iUB,V%j) ! Rank 2 Array case default @@ -2227,6 +2235,8 @@ subroutine ADI_VarUnpackInput(V, ValAry, u) call MV_UnpackMesh(V, ValAry, u%AD%rotors(DL%i1)%BladeMotion(DL%i2)) ! Mesh case (ADI_u_AD_rotors_TFinMotion) call MV_UnpackMesh(V, ValAry, u%AD%rotors(DL%i1)%TFinMotion) ! Mesh + case (ADI_u_AD_rotors_GSMotion) + call MV_UnpackMesh(V, ValAry, u%AD%rotors(DL%i1)%GSMotion) ! Mesh case (ADI_u_AD_rotors_UserProp) u%AD%rotors(DL%i1)%UserProp(V%iLB:V%iUB, V%j) = VarVals ! Rank 2 Array end select @@ -2249,6 +2259,8 @@ function ADI_InputFieldName(DL) result(Name) Name = "u%AD%rotors("//trim(Num2LStr(DL%i1))//")%BladeMotion("//trim(Num2LStr(DL%i2))//")" case (ADI_u_AD_rotors_TFinMotion) Name = "u%AD%rotors("//trim(Num2LStr(DL%i1))//")%TFinMotion" + case (ADI_u_AD_rotors_GSMotion) + Name = "u%AD%rotors("//trim(Num2LStr(DL%i1))//")%GSMotion" case (ADI_u_AD_rotors_UserProp) Name = "u%AD%rotors("//trim(Num2LStr(DL%i1))//")%UserProp" case default @@ -2282,6 +2294,8 @@ subroutine ADI_VarPackOutput(V, y, ValAry) call MV_PackMesh(V, y%AD%rotors(DL%i1)%BladeLoad(DL%i2), ValAry) ! Mesh case (ADI_y_AD_rotors_TFinLoad) call MV_PackMesh(V, y%AD%rotors(DL%i1)%TFinLoad, ValAry) ! Mesh + case (ADI_y_AD_rotors_GSLoad) + call MV_PackMesh(V, y%AD%rotors(DL%i1)%GSLoad, ValAry) ! Mesh case (ADI_y_AD_rotors_WriteOutput) VarVals = y%AD%rotors(DL%i1)%WriteOutput(V%iLB:V%iUB) ! Rank 1 Array case (ADI_y_HHVel) @@ -2324,6 +2338,8 @@ subroutine ADI_VarUnpackOutput(V, ValAry, y) call MV_UnpackMesh(V, ValAry, y%AD%rotors(DL%i1)%BladeLoad(DL%i2)) ! Mesh case (ADI_y_AD_rotors_TFinLoad) call MV_UnpackMesh(V, ValAry, y%AD%rotors(DL%i1)%TFinLoad) ! Mesh + case (ADI_y_AD_rotors_GSLoad) + call MV_UnpackMesh(V, ValAry, y%AD%rotors(DL%i1)%GSLoad) ! Mesh case (ADI_y_AD_rotors_WriteOutput) y%AD%rotors(DL%i1)%WriteOutput(V%iLB:V%iUB) = VarVals ! Rank 1 Array case (ADI_y_HHVel) @@ -2352,6 +2368,8 @@ function ADI_OutputFieldName(DL) result(Name) Name = "y%AD%rotors("//trim(Num2LStr(DL%i1))//")%BladeLoad("//trim(Num2LStr(DL%i2))//")" case (ADI_y_AD_rotors_TFinLoad) Name = "y%AD%rotors("//trim(Num2LStr(DL%i1))//")%TFinLoad" + case (ADI_y_AD_rotors_GSLoad) + Name = "y%AD%rotors("//trim(Num2LStr(DL%i1))//")%GSLoad" case (ADI_y_AD_rotors_WriteOutput) Name = "y%AD%rotors("//trim(Num2LStr(DL%i1))//")%WriteOutput" case (ADI_y_HHVel) diff --git a/modules/aerodyn/src/AeroDyn_Registry.txt b/modules/aerodyn/src/AeroDyn_Registry.txt index 1c1f30babf..a44c0c58e8 100644 --- a/modules/aerodyn/src/AeroDyn_Registry.txt +++ b/modules/aerodyn/src/AeroDyn_Registry.txt @@ -34,9 +34,20 @@ param ^ - IntKi TwrShadow_none - 0 - "no tower s param ^ - IntKi TwrShadow_Powles - 1 - "Powles tower shadow model" - param ^ - IntKi TwrShadow_Eames - 2 - "Eames tower shadow model" - -param ^ - IntKi TwrAero_none - 0 - "no tower aero" - -param ^ - IntKi TwrAero_noVIV - 1 - "Tower aero model without VIV" - -param ^ - IntKi TwrAero_VIV - 2 - "Tower aero model with VIV" - +param ^ - IntKi TwrAero_none - 0 - "no tower aero" - +param ^ - IntKi TwrAero_noVIV - 1 - "Tower aero model without VIV" - +param ^ - IntKi TwrAero_VIV - 2 - "Tower aero model with VIV" - + +param ^ - IntKi GSPotent_none - 0 - "no general support potential flow" - +param ^ - IntKi GSPotent_baseline - 1 - "baseline general support potential flow" - +param ^ - IntKi GSPotent_Bak - 2 - "general support potential flow with Bak correction" - + +param ^ - IntKi GSShadow_none - 0 - "no general support shadow" - +param ^ - IntKi GSShadow_Powles - 1 - "Powles general support shadow model" - +param ^ - IntKi GSShadow_Eames - 2 - "Eames general support shadow model" - + +param ^ - IntKi GSAero_none - 0 - "no general support aero" - +param ^ - IntKi GSAero_noVIV - 1 - "General support aero model without VIV" - param ^ - IntKi SA_Wgt_Uniform - 1 - "Sector average weighting - Uniform" - #param ^ - IntKi SA_Wgt_Impulse - 1 - "Sector average weighting - Impulse" - @@ -92,7 +103,27 @@ typedef ^ AD_VTK_BLSurfaceType SiKi AirfoilCoords typedef ^ AD_VTK_RotSurfaceType AD_VTK_BLSurfaceType BladeShape {:} - - "AirfoilCoords for each blade" - typedef ^ ^ SiKi TowerRad {:} - - "radius of each ED tower node" m - +# General support data types +typedef ^ GSJointType IntKi JointID - - - "User-specified joint ID for the general support structure" - +typedef ^ GSJointType ReKi Position {3} - - "Undisplaced position of the joint" m +typedef ^ GSInpMemberType IntKi MemberID - - - "User-specified member ID for the general support structure" - +typedef ^ GSInpMemberType IntKi MJointID1 - - - "Starting joint ID of the member" - +typedef ^ GSInpMemberType IntKi MJointID2 - - - "Ending joint ID of the member" - +typedef ^ GSInpMemberType ReKi MDiam1 - - - "Starting diameter of the member" m +typedef ^ GSInpMemberType ReKi MDiam2 - - - "Ending diameter of the member" m +typedef ^ GSInpMemberType ReKi MCd1 - - - "Starting drag coefficient of the member" - +typedef ^ GSInpMemberType ReKi MCd2 - - - "Ending drag coefficient of the member" - +typedef ^ GSInpMemberType ReKi MTI1 - - - "Starting turbulence intensity of the member for the Eames tower shadow model" - +typedef ^ GSInpMemberType ReKi MTI2 - - - "Ending turbulence intensity of the member for the Eames tower shadow model" - +typedef ^ GSInpMemberType ReKi MDivSize - - - "Maximum member discretized element length" m +typedef ^ GSMemberType IntKi MemberID - - - "User-supplied integer ID for this member" - +typedef ^ GSMemberType IntKi NElements - - - "Number of elements in this member" - +typedef ^ GSMemberType IntKi NodeIndx {:} - - "Index of each of the member's nodes in the global node list" - +typedef ^ GSMemberType ReKi RefLength - - - "The reference total length for this member" m +typedef ^ GSMemberType ReKi dl - - - "The reference element length of this member" m +typedef ^ GSMemberType ReKi R {:} - - "Outer member radius at each node" m +typedef ^ GSMemberType ReKi Cd {:} - - "Drag coefficient at each node" - +typedef ^ GSMemberType ReKi TI {:} - - "Turbulence intensity at each node for the Eames tower shadow model" - # ..... Initialization data ....................................................................................................... # Define inputs that the initialization routine may need here: @@ -189,6 +220,14 @@ typedef ^ RotInputFile LOGICAL TFinAero - .FALSE. - "Calculate t typedef ^ RotInputFile CHARACTER(1024) TFinFile - - - "Input file for tail fin aerodynamics [used only when TFinAero=True]" - typedef ^ RotInputFile TFinInputFileType TFin - - - "Input file data for tail fin" - +typedef ^ GSInputFile IntKi NJoints - - - "Number of user-specified joints" - +typedef ^ GSInputFile GSJointType InpJoints {:} - - "Array of user-specified joints" - +typedef ^ GSInputFile IntKi NMembers - - - "Number of user-specified members" - +typedef ^ GSInputFile GSInpMemberType InpMembers {:} - - "Array of user-specified members" - +typedef ^ GSInputFile IntKi GSPotent - - - "Type of general support influence based on potential flow around the member {0=none, 1=baseline potential flow, 2=potential flow with Bak correction}" - +typedef ^ GSInputFile IntKi GSShadow - - - "Type of general support downstream shadow effect {0=none, 1=Powles model, 2=Eames model}" - +typedef ^ GSInputFile IntKi GSAero - - - "Calculate multi-member generalized tower aerodynamic loads? {0=none, 1=aero without VIV, 2=aero with VIV}" - + typedef ^ AD_InputFile Logical Echo - - - "Echo input file to echo file" - typedef ^ AD_InputFile DbKi DTAero - - - "Time interval for aerodynamic calculations {or \"default\"}" s typedef ^ AD_InputFile IntKi Wake_Mod - - - "Type of wake/induction model {0=none, 1=BEMT, 2=DBEMT, 3=FVW}" - @@ -253,7 +292,7 @@ typedef ^ AD_InputFile IntKi BldNd_BladesOut - - - "The blades to output (AD typedef ^ AD_InputFile ReKi UAStartRad - - - Starting radius for dynamic stall (fraction of rotor radius) typedef ^ AD_InputFile ReKi UAEndRad - - - Ending radius for dynamic stall (fraction of rotor radius) typedef ^ AD_InputFile RotInputFile rotors {:} - - "Rotor (blades and tower) input file data" - - +typedef ^ AD_InputFile GSInputFile GS - - - "General support input file data" - @@ -290,7 +329,7 @@ typedef ^ OtherStateType ReKi WakeLocationPoints {:}{:} - - "wake points velocit typedef ^ ElemInflowType ReKi InflowVel {:}{:} - - "U,V,W at nodes on element (note if we change the requirement that NumNodes is the same for each blade, this will need to change)" m/s typedef ^ ElemInflowType ReKi InflowAcc {:}{:} - - "Wind acceleration at nodes on element (blade or tower) (note if we change the requirement that NumNodes is the same for each blade, this will need to change)" m/s typedef ^ RotInflowType ElemInflowType Blade {:} - - "Blade wind inputs" - -typedef ^ RotInflowType ElemInflowType Tower - - - "Blade wind inputs" - +typedef ^ RotInflowType ElemInflowType Tower - - - "Tower wind inputs" - typedef ^ RotInflowType ReKi InflowOnHub {3}{1} - - "U,V,W at hub" m/s typedef ^ RotInflowType ReKi InflowOnNacelle {3}{1} - - "U,V,W at nacelle" m/s typedef ^ RotInflowType ReKi InflowOnTailFin {3}{1} - - "U,V,W at tailfin" m/s @@ -298,6 +337,7 @@ typedef ^ RotInflowType ReKi AvgDiskVel {3} - 0.0 "disk-averaged U,V,W" m/s typedef ^ AD_InflowType ReKi InflowWakeVel {:}{:} - - "U,V,W at wake points" m/s typedef ^ AD_InflowType RotInflowType RotInflow {:} - - "Inflow on rotor" - +typedef ^ AD_InflowType ElemInflowType GSInflow - - - "Inflow on general support structure" - # ..... Parameters ................................................................................................................ @@ -385,7 +425,24 @@ typedef ^ RotParameterType IntKi BldNd_NumNodesOut - - - "The blades to outpu # Tail fin parameters (per rotor) typedef ^ RotParameterType LOGICAL TFinAero - .FALSE. - "Calculate tail fin aerodynamics model (flag)" flag typedef ^ RotParameterType TFinParameterType TFin - - - "Parameters for tail fin of current rotor" - - +# General support parameters (per rotor) +typedef ^ RotParameterType LOGICAL hasGSMod - .FALSE. - "Has general support structure model (flag)" flag +typedef ^ RotParameterType IntKi GSPotent - - - "Type of general support influence based on potential flow around the member {0=none, 1=baseline potential flow, 2=potential flow with Bak correction}" - +typedef ^ RotParameterType IntKi GSShadow - - - "Type of general support downstream shadow effect {0=none, 1=Powles model, 2=Eames model}" - +typedef ^ RotParameterType IntKi GSAero - - - "Type of flow-induced loads on the general support {0=none, 1=Drag}" - +# General support parameters (shared) +typedef ^ GSParameterType LOGICAL hasGSMod - .FALSE. - "Has general support structure model (flag)" flag +typedef ^ GSParameterType IntKi GSPotent - - - "Type of general support influence based on potential flow around the member {0=none, 1=baseline potential flow, 2=potential flow with Bak correction}" - +typedef ^ GSParameterType IntKi GSShadow - - - "Type of general support downstream shadow effect {0=none, 1=Powles model, 2=Eames model}" - +typedef ^ GSParameterType IntKi GSAero - - - "Type of flow-induced loads on the general support {0=none, 1=Drag}" - +typedef ^ GSParameterType IntKi NJoints - - - "Total number of user-specified joints for the general support structure" - +typedef ^ GSParameterType IntKi NNodes - - - "Total number of general support structure nodes (joints+interior nodes)" - +typedef ^ GSParameterType IntKi NMembers - - - "Total number of general support structure members" - +typedef ^ GSParameterType GSMemberType Members {:} - - "General support structure members" - +typedef ^ GSParameterType GSJointType Joints {:} - - "General support structure joints" - +typedef ^ GSParameterType ReKi Density - - - "Fluid density" kg/m^3 +typedef ^ GSParameterType IntKi MHK - - - "MHK turbine type switch" - +typedef ^ GSParameterType ReKi WtrDpth - - - "Water depth" m # Parameters for all rotors: typedef ^ ParameterType RotParameterType rotors {:} - - "Parameter types for each rotor" - # Time step for integration of continuous states (if a fixed-step integrator is used) and update of discrete states: @@ -405,6 +462,7 @@ typedef ^ ^ ReKi SA_PsiBwd - - - "Sector Average typedef ^ ^ ReKi SA_PsiFwd - - - "Sector Average - Forward Azimuth (>0)" deg typedef ^ ^ IntKi SA_nPerSec - - - "Sector Average - Number of points per sector (>1)" - typedef ^ ParameterType SeaSt_WaveFieldType *WaveField - - - "Pointer to SeaState wave field data type" - +typedef ^ ParameterType GSParameterType GS - - - "General support structure parameters" - # ..... Inputs .................................................................................................................... # Define inputs that are contained on a mesh here: @@ -414,6 +472,7 @@ typedef ^ RotInputType MeshType HubMotion - - - "motion on the hub" - typedef ^ RotInputType MeshType BladeRootMotion {:} - - "motion on each blade root" - typedef ^ RotInputType MeshType BladeMotion {:} - - "motion on each blade" - typedef ^ RotInputType MeshType TFinMotion - - - "motion of tail fin (at tail fin ref point)" - +typedef ^ RotInputType MeshType GSMotion - - - "motion of the general support structure" - # Define inputs that are not on a mesh here: typedef ^ RotInputType ReKi UserProp {:}{:} - - "Optional user property for interpolating airfoils (per element per blade)" - @@ -428,6 +487,7 @@ typedef ^ RotOutputType MeshType HubLoad - - - "loads on the hub" - typedef ^ RotOutputType MeshType TowerLoad - - - "loads on the tower" - typedef ^ RotOutputType MeshType BladeLoad {:} - - "loads on each blade" - typedef ^ RotOutputType MeshType TFinLoad - - - "loads on tail fin (at tail fin ref point)" - +typedef ^ RotOutputType MeshType GSLoad - - - "loads on the general support structure" - # Define outputs that are not on a mesh here: typedef ^ RotOutputType ReKi WriteOutput {:} - - "Data to be written to an output file: see WriteOutputHdr for names of each variable" "see WriteOutputUnt" @@ -463,6 +523,7 @@ typedef ^ RotMiscVarType ReKi Y_Twr {:} - - "local y-component of force per unit typedef ^ RotMiscVarType ReKi Cant {:}{:} - - "curvature angle, saved for possible output to file" rad typedef ^ RotMiscVarType ReKi Toe {:}{:} - - "Toe angle, saved for possible output to file" rad typedef ^ RotMiscVarType ReKi TwrClrnc {:}{:} - - "Distance between tower (including tower radius) and blade node (not including blade width), saved for possible output to file" m +typedef ^ RotMiscVarType ReKi GSClrnc {:}{:} - - "Distance between nearest GS member (including member radius) and blade node (not including blade width), saved for possible output to file" m typedef ^ RotMiscVarType ReKi X {:}{:} - - "normal force per unit length (normal to the plane, not chord) of the jth node in the kth blade" N/m typedef ^ RotMiscVarType ReKi Y {:}{:} - - "tangential force per unit length (tangential to the plane, not chord) of the jth node in the kth blade" N/m typedef ^ RotMiscVarType ReKi Z {:}{:} - - "axial force per unit length (tangential to the plane, not chord) of the jth node in the kth blade" N/m diff --git a/modules/aerodyn/src/AeroDyn_Types.f90 b/modules/aerodyn/src/AeroDyn_Types.f90 index 01447c6f27..ff49e7ff1b 100644 --- a/modules/aerodyn/src/AeroDyn_Types.f90 +++ b/modules/aerodyn/src/AeroDyn_Types.f90 @@ -52,6 +52,14 @@ MODULE AeroDyn_Types INTEGER(IntKi), PUBLIC, PARAMETER :: TwrAero_none = 0 ! no tower aero [-] INTEGER(IntKi), PUBLIC, PARAMETER :: TwrAero_noVIV = 1 ! Tower aero model without VIV [-] INTEGER(IntKi), PUBLIC, PARAMETER :: TwrAero_VIV = 2 ! Tower aero model with VIV [-] + INTEGER(IntKi), PUBLIC, PARAMETER :: GSPotent_none = 0 ! no general support potential flow [-] + INTEGER(IntKi), PUBLIC, PARAMETER :: GSPotent_baseline = 1 ! baseline general support potential flow [-] + INTEGER(IntKi), PUBLIC, PARAMETER :: GSPotent_Bak = 2 ! general support potential flow with Bak correction [-] + INTEGER(IntKi), PUBLIC, PARAMETER :: GSShadow_none = 0 ! no general support shadow [-] + INTEGER(IntKi), PUBLIC, PARAMETER :: GSShadow_Powles = 1 ! Powles general support shadow model [-] + INTEGER(IntKi), PUBLIC, PARAMETER :: GSShadow_Eames = 2 ! Eames general support shadow model [-] + INTEGER(IntKi), PUBLIC, PARAMETER :: GSAero_none = 0 ! no general support aero [-] + INTEGER(IntKi), PUBLIC, PARAMETER :: GSAero_noVIV = 1 ! General support aero model without VIV [-] INTEGER(IntKi), PUBLIC, PARAMETER :: SA_Wgt_Uniform = 1 ! Sector average weighting - Uniform [-] INTEGER(IntKi), PUBLIC, PARAMETER :: TFinAero_none = 0 ! no tail fin aero [-] INTEGER(IntKi), PUBLIC, PARAMETER :: TFinAero_polar = 1 ! polar-based tail fin aerodynamics [-] @@ -106,6 +114,38 @@ MODULE AeroDyn_Types REAL(SiKi) , DIMENSION(:), ALLOCATABLE :: TowerRad !< radius of each ED tower node [m] END TYPE AD_VTK_RotSurfaceType ! ======================= +! ========= GSJointType ======= + TYPE, PUBLIC :: GSJointType + INTEGER(IntKi) :: JointID = 0_IntKi !< User-specified joint ID for the general support structure [-] + REAL(ReKi) , DIMENSION(1:3) :: Position = 0.0_ReKi !< Undisplaced position of the joint [m] + END TYPE GSJointType +! ======================= +! ========= GSInpMemberType ======= + TYPE, PUBLIC :: GSInpMemberType + INTEGER(IntKi) :: MemberID = 0_IntKi !< User-specified member ID for the general support structure [-] + INTEGER(IntKi) :: MJointID1 = 0_IntKi !< Starting joint ID of the member [-] + INTEGER(IntKi) :: MJointID2 = 0_IntKi !< Ending joint ID of the member [-] + REAL(ReKi) :: MDiam1 = 0.0_ReKi !< Starting diameter of the member [m] + REAL(ReKi) :: MDiam2 = 0.0_ReKi !< Ending diameter of the member [m] + REAL(ReKi) :: MCd1 = 0.0_ReKi !< Starting drag coefficient of the member [-] + REAL(ReKi) :: MCd2 = 0.0_ReKi !< Ending drag coefficient of the member [-] + REAL(ReKi) :: MTI1 = 0.0_ReKi !< Starting turbulence intensity of the member for the Eames tower shadow model [-] + REAL(ReKi) :: MTI2 = 0.0_ReKi !< Ending turbulence intensity of the member for the Eames tower shadow model [-] + REAL(ReKi) :: MDivSize = 0.0_ReKi !< Maximum member discretized element length [m] + END TYPE GSInpMemberType +! ======================= +! ========= GSMemberType ======= + TYPE, PUBLIC :: GSMemberType + INTEGER(IntKi) :: MemberID = 0_IntKi !< User-supplied integer ID for this member [-] + INTEGER(IntKi) :: NElements = 0_IntKi !< Number of elements in this member [-] + INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: NodeIndx !< Index of each of the member's nodes in the global node list [-] + REAL(ReKi) :: RefLength = 0.0_ReKi !< The reference total length for this member [m] + REAL(ReKi) :: dl = 0.0_ReKi !< The reference element length of this member [m] + REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: R !< Outer member radius at each node [m] + REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: Cd !< Drag coefficient at each node [-] + REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: TI !< Turbulence intensity at each node for the Eames tower shadow model [-] + END TYPE GSMemberType +! ======================= ! ========= RotInitInputType ======= TYPE, PUBLIC :: RotInitInputType INTEGER(IntKi) :: NumBlades = 0_IntKi !< Number of blades on the turbine [-] @@ -210,6 +250,17 @@ MODULE AeroDyn_Types TYPE(TFinInputFileType) :: TFin !< Input file data for tail fin [-] END TYPE RotInputFile ! ======================= +! ========= GSInputFile ======= + TYPE, PUBLIC :: GSInputFile + INTEGER(IntKi) :: NJoints = 0_IntKi !< Number of user-specified joints [-] + TYPE(GSJointType) , DIMENSION(:), ALLOCATABLE :: InpJoints !< Array of user-specified joints [-] + INTEGER(IntKi) :: NMembers = 0_IntKi !< Number of user-specified members [-] + TYPE(GSInpMemberType) , DIMENSION(:), ALLOCATABLE :: InpMembers !< Array of user-specified members [-] + INTEGER(IntKi) :: GSPotent = 0_IntKi !< Type of general support influence based on potential flow around the member {0=none, 1=baseline potential flow, 2=potential flow with Bak correction} [-] + INTEGER(IntKi) :: GSShadow = 0_IntKi !< Type of general support downstream shadow effect {0=none, 1=Powles model, 2=Eames model} [-] + INTEGER(IntKi) :: GSAero = 0_IntKi !< Calculate multi-member generalized tower aerodynamic loads? {0=none, 1=aero without VIV, 2=aero with VIV} [-] + END TYPE GSInputFile +! ======================= ! ========= AD_InputFile ======= TYPE, PUBLIC :: AD_InputFile LOGICAL :: Echo = .false. !< Echo input file to echo file [-] @@ -273,6 +324,7 @@ MODULE AeroDyn_Types REAL(ReKi) :: UAStartRad = 0.0_ReKi !< Starting [radius] REAL(ReKi) :: UAEndRad = 0.0_ReKi !< Ending [radius] TYPE(RotInputFile) , DIMENSION(:), ALLOCATABLE :: rotors !< Rotor (blades and tower) input file data [-] + TYPE(GSInputFile) :: GS !< General support input file data [-] END TYPE AD_InputFile ! ======================= ! ========= RotContinuousStateType ======= @@ -331,7 +383,7 @@ MODULE AeroDyn_Types ! ========= RotInflowType ======= TYPE, PUBLIC :: RotInflowType TYPE(ElemInflowType) , DIMENSION(:), ALLOCATABLE :: Blade !< Blade wind inputs [-] - TYPE(ElemInflowType) :: Tower !< Blade wind inputs [-] + TYPE(ElemInflowType) :: Tower !< Tower wind inputs [-] REAL(ReKi) , DIMENSION(1:3,1:1) :: InflowOnHub = 0.0_ReKi !< U,V,W at hub [m/s] REAL(ReKi) , DIMENSION(1:3,1:1) :: InflowOnNacelle = 0.0_ReKi !< U,V,W at nacelle [m/s] REAL(ReKi) , DIMENSION(1:3,1:1) :: InflowOnTailFin = 0.0_ReKi !< U,V,W at tailfin [m/s] @@ -342,6 +394,7 @@ MODULE AeroDyn_Types TYPE, PUBLIC :: AD_InflowType REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: InflowWakeVel !< U,V,W at wake points [m/s] TYPE(RotInflowType) , DIMENSION(:), ALLOCATABLE :: RotInflow !< Inflow on rotor [-] + TYPE(ElemInflowType) :: GSInflow !< Inflow on general support structure [-] END TYPE AD_InflowType ! ======================= ! ========= RotParameterType ======= @@ -422,8 +475,28 @@ MODULE AeroDyn_Types INTEGER(IntKi) :: BldNd_NumNodesOut = 0_IntKi !< The blades to output (AD_AllBldNdOuts) [-] LOGICAL :: TFinAero = .FALSE. !< Calculate tail fin aerodynamics model (flag) [flag] TYPE(TFinParameterType) :: TFin !< Parameters for tail fin of current rotor [-] + LOGICAL :: hasGSMod = .FALSE. !< Has general support structure model (flag) [flag] + INTEGER(IntKi) :: GSPotent = 0_IntKi !< Type of general support influence based on potential flow around the member {0=none, 1=baseline potential flow, 2=potential flow with Bak correction} [-] + INTEGER(IntKi) :: GSShadow = 0_IntKi !< Type of general support downstream shadow effect {0=none, 1=Powles model, 2=Eames model} [-] + INTEGER(IntKi) :: GSAero = 0_IntKi !< Type of flow-induced loads on the general support {0=none, 1=Drag} [-] END TYPE RotParameterType ! ======================= +! ========= GSParameterType ======= + TYPE, PUBLIC :: GSParameterType + LOGICAL :: hasGSMod = .FALSE. !< Has general support structure model (flag) [flag] + INTEGER(IntKi) :: GSPotent = 0_IntKi !< Type of general support influence based on potential flow around the member {0=none, 1=baseline potential flow, 2=potential flow with Bak correction} [-] + INTEGER(IntKi) :: GSShadow = 0_IntKi !< Type of general support downstream shadow effect {0=none, 1=Powles model, 2=Eames model} [-] + INTEGER(IntKi) :: GSAero = 0_IntKi !< Type of flow-induced loads on the general support {0=none, 1=Drag} [-] + INTEGER(IntKi) :: NJoints = 0_IntKi !< Total number of user-specified joints for the general support structure [-] + INTEGER(IntKi) :: NNodes = 0_IntKi !< Total number of general support structure nodes (joints+interior nodes) [-] + INTEGER(IntKi) :: NMembers = 0_IntKi !< Total number of general support structure members [-] + TYPE(GSMemberType) , DIMENSION(:), ALLOCATABLE :: Members !< General support structure members [-] + TYPE(GSJointType) , DIMENSION(:), ALLOCATABLE :: Joints !< General support structure joints [-] + REAL(ReKi) :: Density = 0.0_ReKi !< Fluid density [kg/m^3] + INTEGER(IntKi) :: MHK = 0_IntKi !< MHK turbine type switch [-] + REAL(ReKi) :: WtrDpth = 0.0_ReKi !< Water depth [m] + END TYPE GSParameterType +! ======================= ! ========= AD_ParameterType ======= TYPE, PUBLIC :: AD_ParameterType TYPE(RotParameterType) , DIMENSION(:), ALLOCATABLE :: rotors !< Parameter types for each rotor [-] @@ -443,6 +516,7 @@ MODULE AeroDyn_Types REAL(ReKi) :: SA_PsiFwd = 0.0_ReKi !< Sector Average - Forward Azimuth (>0) [deg] INTEGER(IntKi) :: SA_nPerSec = 0_IntKi !< Sector Average - Number of points per sector (>1) [-] TYPE(SeaSt_WaveFieldType) , POINTER :: WaveField => NULL() !< Pointer to SeaState wave field data type [-] + TYPE(GSParameterType) :: GS !< General support structure parameters [-] END TYPE AD_ParameterType ! ======================= ! ========= RotInputType ======= @@ -453,6 +527,7 @@ MODULE AeroDyn_Types TYPE(MeshType) , DIMENSION(:), ALLOCATABLE :: BladeRootMotion !< motion on each blade root [-] TYPE(MeshType) , DIMENSION(:), ALLOCATABLE :: BladeMotion !< motion on each blade [-] TYPE(MeshType) :: TFinMotion !< motion of tail fin (at tail fin ref point) [-] + TYPE(MeshType) :: GSMotion !< motion of the general support structure [-] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: UserProp !< Optional user property for interpolating airfoils (per element per blade) [-] END TYPE RotInputType ! ======================= @@ -468,6 +543,7 @@ MODULE AeroDyn_Types TYPE(MeshType) :: TowerLoad !< loads on the tower [-] TYPE(MeshType) , DIMENSION(:), ALLOCATABLE :: BladeLoad !< loads on each blade [-] TYPE(MeshType) :: TFinLoad !< loads on tail fin (at tail fin ref point) [-] + TYPE(MeshType) :: GSLoad !< loads on the general support structure [-] REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: WriteOutput !< Data to be written to an output file: see WriteOutputHdr for names of each variable [see WriteOutputUnt] END TYPE RotOutputType ! ======================= @@ -504,6 +580,7 @@ MODULE AeroDyn_Types REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: Cant !< curvature angle, saved for possible output to file [rad] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: Toe !< Toe angle, saved for possible output to file [rad] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: TwrClrnc !< Distance between tower (including tower radius) and blade node (not including blade width), saved for possible output to file [m] + REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: GSClrnc !< Distance between nearest GS member (including member radius) and blade node (not including blade width), saved for possible output to file [m] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: X !< normal force per unit length (normal to the plane, not chord) of the jth node in the kth blade [N/m] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: Y !< tangential force per unit length (tangential to the plane, not chord) of the jth node in the kth blade [N/m] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: Z !< axial force per unit length (tangential to the plane, not chord) of the jth node in the kth blade [N/m] @@ -585,13 +662,15 @@ MODULE AeroDyn_Types integer(IntKi), public, parameter :: AD_u_BladeRootMotion = 8 ! AD%BladeRootMotion(DL%i1) integer(IntKi), public, parameter :: AD_u_BladeMotion = 9 ! AD%BladeMotion(DL%i1) integer(IntKi), public, parameter :: AD_u_TFinMotion = 10 ! AD%TFinMotion - integer(IntKi), public, parameter :: AD_u_UserProp = 11 ! AD%UserProp - integer(IntKi), public, parameter :: AD_y_NacelleLoad = 12 ! AD%NacelleLoad - integer(IntKi), public, parameter :: AD_y_HubLoad = 13 ! AD%HubLoad - integer(IntKi), public, parameter :: AD_y_TowerLoad = 14 ! AD%TowerLoad - integer(IntKi), public, parameter :: AD_y_BladeLoad = 15 ! AD%BladeLoad(DL%i1) - integer(IntKi), public, parameter :: AD_y_TFinLoad = 16 ! AD%TFinLoad - integer(IntKi), public, parameter :: AD_y_WriteOutput = 17 ! AD%WriteOutput + integer(IntKi), public, parameter :: AD_u_GSMotion = 11 ! AD%GSMotion + integer(IntKi), public, parameter :: AD_u_UserProp = 12 ! AD%UserProp + integer(IntKi), public, parameter :: AD_y_NacelleLoad = 13 ! AD%NacelleLoad + integer(IntKi), public, parameter :: AD_y_HubLoad = 14 ! AD%HubLoad + integer(IntKi), public, parameter :: AD_y_TowerLoad = 15 ! AD%TowerLoad + integer(IntKi), public, parameter :: AD_y_BladeLoad = 16 ! AD%BladeLoad(DL%i1) + integer(IntKi), public, parameter :: AD_y_TFinLoad = 17 ! AD%TFinLoad + integer(IntKi), public, parameter :: AD_y_GSLoad = 18 ! AD%GSLoad + integer(IntKi), public, parameter :: AD_y_WriteOutput = 19 ! AD%WriteOutput contains @@ -901,6 +980,232 @@ subroutine AD_UnPackVTK_RotSurfaceType(RF, OutData) call RegUnpackAlloc(RF, OutData%TowerRad); if (RegCheckErr(RF, RoutineName)) return end subroutine +subroutine AD_CopyGSJointType(SrcGSJointTypeData, DstGSJointTypeData, CtrlCode, ErrStat, ErrMsg) + type(GSJointType), intent(in) :: SrcGSJointTypeData + type(GSJointType), intent(inout) :: DstGSJointTypeData + integer(IntKi), intent(in ) :: CtrlCode + integer(IntKi), intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + character(*), parameter :: RoutineName = 'AD_CopyGSJointType' + ErrStat = ErrID_None + ErrMsg = '' + DstGSJointTypeData%JointID = SrcGSJointTypeData%JointID + DstGSJointTypeData%Position = SrcGSJointTypeData%Position +end subroutine + +subroutine AD_DestroyGSJointType(GSJointTypeData, ErrStat, ErrMsg) + type(GSJointType), intent(inout) :: GSJointTypeData + integer(IntKi), intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + character(*), parameter :: RoutineName = 'AD_DestroyGSJointType' + ErrStat = ErrID_None + ErrMsg = '' +end subroutine + +subroutine AD_PackGSJointType(RF, Indata) + type(RegFile), intent(inout) :: RF + type(GSJointType), intent(in) :: InData + character(*), parameter :: RoutineName = 'AD_PackGSJointType' + if (RF%ErrStat >= AbortErrLev) return + call RegPack(RF, InData%JointID) + call RegPack(RF, InData%Position) + if (RegCheckErr(RF, RoutineName)) return +end subroutine + +subroutine AD_UnPackGSJointType(RF, OutData) + type(RegFile), intent(inout) :: RF + type(GSJointType), intent(inout) :: OutData + character(*), parameter :: RoutineName = 'AD_UnPackGSJointType' + if (RF%ErrStat /= ErrID_None) return + call RegUnpack(RF, OutData%JointID); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%Position); if (RegCheckErr(RF, RoutineName)) return +end subroutine + +subroutine AD_CopyGSInpMemberType(SrcGSInpMemberTypeData, DstGSInpMemberTypeData, CtrlCode, ErrStat, ErrMsg) + type(GSInpMemberType), intent(in) :: SrcGSInpMemberTypeData + type(GSInpMemberType), intent(inout) :: DstGSInpMemberTypeData + integer(IntKi), intent(in ) :: CtrlCode + integer(IntKi), intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + character(*), parameter :: RoutineName = 'AD_CopyGSInpMemberType' + ErrStat = ErrID_None + ErrMsg = '' + DstGSInpMemberTypeData%MemberID = SrcGSInpMemberTypeData%MemberID + DstGSInpMemberTypeData%MJointID1 = SrcGSInpMemberTypeData%MJointID1 + DstGSInpMemberTypeData%MJointID2 = SrcGSInpMemberTypeData%MJointID2 + DstGSInpMemberTypeData%MDiam1 = SrcGSInpMemberTypeData%MDiam1 + DstGSInpMemberTypeData%MDiam2 = SrcGSInpMemberTypeData%MDiam2 + DstGSInpMemberTypeData%MCd1 = SrcGSInpMemberTypeData%MCd1 + DstGSInpMemberTypeData%MCd2 = SrcGSInpMemberTypeData%MCd2 + DstGSInpMemberTypeData%MTI1 = SrcGSInpMemberTypeData%MTI1 + DstGSInpMemberTypeData%MTI2 = SrcGSInpMemberTypeData%MTI2 + DstGSInpMemberTypeData%MDivSize = SrcGSInpMemberTypeData%MDivSize +end subroutine + +subroutine AD_DestroyGSInpMemberType(GSInpMemberTypeData, ErrStat, ErrMsg) + type(GSInpMemberType), intent(inout) :: GSInpMemberTypeData + integer(IntKi), intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + character(*), parameter :: RoutineName = 'AD_DestroyGSInpMemberType' + ErrStat = ErrID_None + ErrMsg = '' +end subroutine + +subroutine AD_PackGSInpMemberType(RF, Indata) + type(RegFile), intent(inout) :: RF + type(GSInpMemberType), intent(in) :: InData + character(*), parameter :: RoutineName = 'AD_PackGSInpMemberType' + if (RF%ErrStat >= AbortErrLev) return + call RegPack(RF, InData%MemberID) + call RegPack(RF, InData%MJointID1) + call RegPack(RF, InData%MJointID2) + call RegPack(RF, InData%MDiam1) + call RegPack(RF, InData%MDiam2) + call RegPack(RF, InData%MCd1) + call RegPack(RF, InData%MCd2) + call RegPack(RF, InData%MTI1) + call RegPack(RF, InData%MTI2) + call RegPack(RF, InData%MDivSize) + if (RegCheckErr(RF, RoutineName)) return +end subroutine + +subroutine AD_UnPackGSInpMemberType(RF, OutData) + type(RegFile), intent(inout) :: RF + type(GSInpMemberType), intent(inout) :: OutData + character(*), parameter :: RoutineName = 'AD_UnPackGSInpMemberType' + if (RF%ErrStat /= ErrID_None) return + call RegUnpack(RF, OutData%MemberID); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%MJointID1); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%MJointID2); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%MDiam1); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%MDiam2); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%MCd1); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%MCd2); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%MTI1); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%MTI2); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%MDivSize); if (RegCheckErr(RF, RoutineName)) return +end subroutine + +subroutine AD_CopyGSMemberType(SrcGSMemberTypeData, DstGSMemberTypeData, CtrlCode, ErrStat, ErrMsg) + type(GSMemberType), intent(in) :: SrcGSMemberTypeData + type(GSMemberType), intent(inout) :: DstGSMemberTypeData + integer(IntKi), intent(in ) :: CtrlCode + integer(IntKi), intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + integer(B4Ki) :: LB(1), UB(1) + integer(IntKi) :: ErrStat2 + character(*), parameter :: RoutineName = 'AD_CopyGSMemberType' + ErrStat = ErrID_None + ErrMsg = '' + DstGSMemberTypeData%MemberID = SrcGSMemberTypeData%MemberID + DstGSMemberTypeData%NElements = SrcGSMemberTypeData%NElements + if (allocated(SrcGSMemberTypeData%NodeIndx)) then + LB(1:1) = lbound(SrcGSMemberTypeData%NodeIndx) + UB(1:1) = ubound(SrcGSMemberTypeData%NodeIndx) + if (.not. allocated(DstGSMemberTypeData%NodeIndx)) then + allocate(DstGSMemberTypeData%NodeIndx(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstGSMemberTypeData%NodeIndx.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstGSMemberTypeData%NodeIndx = SrcGSMemberTypeData%NodeIndx + end if + DstGSMemberTypeData%RefLength = SrcGSMemberTypeData%RefLength + DstGSMemberTypeData%dl = SrcGSMemberTypeData%dl + if (allocated(SrcGSMemberTypeData%R)) then + LB(1:1) = lbound(SrcGSMemberTypeData%R) + UB(1:1) = ubound(SrcGSMemberTypeData%R) + if (.not. allocated(DstGSMemberTypeData%R)) then + allocate(DstGSMemberTypeData%R(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstGSMemberTypeData%R.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstGSMemberTypeData%R = SrcGSMemberTypeData%R + end if + if (allocated(SrcGSMemberTypeData%Cd)) then + LB(1:1) = lbound(SrcGSMemberTypeData%Cd) + UB(1:1) = ubound(SrcGSMemberTypeData%Cd) + if (.not. allocated(DstGSMemberTypeData%Cd)) then + allocate(DstGSMemberTypeData%Cd(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstGSMemberTypeData%Cd.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstGSMemberTypeData%Cd = SrcGSMemberTypeData%Cd + end if + if (allocated(SrcGSMemberTypeData%TI)) then + LB(1:1) = lbound(SrcGSMemberTypeData%TI) + UB(1:1) = ubound(SrcGSMemberTypeData%TI) + if (.not. allocated(DstGSMemberTypeData%TI)) then + allocate(DstGSMemberTypeData%TI(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstGSMemberTypeData%TI.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstGSMemberTypeData%TI = SrcGSMemberTypeData%TI + end if +end subroutine + +subroutine AD_DestroyGSMemberType(GSMemberTypeData, ErrStat, ErrMsg) + type(GSMemberType), intent(inout) :: GSMemberTypeData + integer(IntKi), intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + character(*), parameter :: RoutineName = 'AD_DestroyGSMemberType' + ErrStat = ErrID_None + ErrMsg = '' + if (allocated(GSMemberTypeData%NodeIndx)) then + deallocate(GSMemberTypeData%NodeIndx) + end if + if (allocated(GSMemberTypeData%R)) then + deallocate(GSMemberTypeData%R) + end if + if (allocated(GSMemberTypeData%Cd)) then + deallocate(GSMemberTypeData%Cd) + end if + if (allocated(GSMemberTypeData%TI)) then + deallocate(GSMemberTypeData%TI) + end if +end subroutine + +subroutine AD_PackGSMemberType(RF, Indata) + type(RegFile), intent(inout) :: RF + type(GSMemberType), intent(in) :: InData + character(*), parameter :: RoutineName = 'AD_PackGSMemberType' + if (RF%ErrStat >= AbortErrLev) return + call RegPack(RF, InData%MemberID) + call RegPack(RF, InData%NElements) + call RegPackAlloc(RF, InData%NodeIndx) + call RegPack(RF, InData%RefLength) + call RegPack(RF, InData%dl) + call RegPackAlloc(RF, InData%R) + call RegPackAlloc(RF, InData%Cd) + call RegPackAlloc(RF, InData%TI) + if (RegCheckErr(RF, RoutineName)) return +end subroutine + +subroutine AD_UnPackGSMemberType(RF, OutData) + type(RegFile), intent(inout) :: RF + type(GSMemberType), intent(inout) :: OutData + character(*), parameter :: RoutineName = 'AD_UnPackGSMemberType' + integer(B4Ki) :: LB(1), UB(1) + integer(IntKi) :: stat + logical :: IsAllocAssoc + if (RF%ErrStat /= ErrID_None) return + call RegUnpack(RF, OutData%MemberID); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%NElements); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%NodeIndx); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%RefLength); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%dl); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%R); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%Cd); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%TI); if (RegCheckErr(RF, RoutineName)) return +end subroutine + subroutine AD_CopyRotInitInputType(SrcRotInitInputTypeData, DstRotInitInputTypeData, CtrlCode, ErrStat, ErrMsg) type(RotInitInputType), intent(in) :: SrcRotInitInputTypeData type(RotInitInputType), intent(inout) :: DstRotInitInputTypeData @@ -2128,6 +2433,164 @@ subroutine AD_UnPackRotInputFile(RF, OutData) call AD_UnpackTFinInputFileType(RF, OutData%TFin) ! TFin end subroutine +subroutine AD_CopyGSInputFile(SrcGSInputFileData, DstGSInputFileData, CtrlCode, ErrStat, ErrMsg) + type(GSInputFile), intent(in) :: SrcGSInputFileData + type(GSInputFile), intent(inout) :: DstGSInputFileData + integer(IntKi), intent(in ) :: CtrlCode + integer(IntKi), intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + integer(B4Ki) :: i1 + integer(B4Ki) :: LB(1), UB(1) + integer(IntKi) :: ErrStat2 + character(ErrMsgLen) :: ErrMsg2 + character(*), parameter :: RoutineName = 'AD_CopyGSInputFile' + ErrStat = ErrID_None + ErrMsg = '' + DstGSInputFileData%NJoints = SrcGSInputFileData%NJoints + if (allocated(SrcGSInputFileData%InpJoints)) then + LB(1:1) = lbound(SrcGSInputFileData%InpJoints) + UB(1:1) = ubound(SrcGSInputFileData%InpJoints) + if (.not. allocated(DstGSInputFileData%InpJoints)) then + allocate(DstGSInputFileData%InpJoints(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstGSInputFileData%InpJoints.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + do i1 = LB(1), UB(1) + call AD_CopyGSJointType(SrcGSInputFileData%InpJoints(i1), DstGSInputFileData%InpJoints(i1), CtrlCode, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + if (ErrStat >= AbortErrLev) return + end do + end if + DstGSInputFileData%NMembers = SrcGSInputFileData%NMembers + if (allocated(SrcGSInputFileData%InpMembers)) then + LB(1:1) = lbound(SrcGSInputFileData%InpMembers) + UB(1:1) = ubound(SrcGSInputFileData%InpMembers) + if (.not. allocated(DstGSInputFileData%InpMembers)) then + allocate(DstGSInputFileData%InpMembers(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstGSInputFileData%InpMembers.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + do i1 = LB(1), UB(1) + call AD_CopyGSInpMemberType(SrcGSInputFileData%InpMembers(i1), DstGSInputFileData%InpMembers(i1), CtrlCode, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + if (ErrStat >= AbortErrLev) return + end do + end if + DstGSInputFileData%GSPotent = SrcGSInputFileData%GSPotent + DstGSInputFileData%GSShadow = SrcGSInputFileData%GSShadow + DstGSInputFileData%GSAero = SrcGSInputFileData%GSAero +end subroutine + +subroutine AD_DestroyGSInputFile(GSInputFileData, ErrStat, ErrMsg) + type(GSInputFile), intent(inout) :: GSInputFileData + integer(IntKi), intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + integer(B4Ki) :: i1 + integer(B4Ki) :: LB(1), UB(1) + integer(IntKi) :: ErrStat2 + character(ErrMsgLen) :: ErrMsg2 + character(*), parameter :: RoutineName = 'AD_DestroyGSInputFile' + ErrStat = ErrID_None + ErrMsg = '' + if (allocated(GSInputFileData%InpJoints)) then + LB(1:1) = lbound(GSInputFileData%InpJoints) + UB(1:1) = ubound(GSInputFileData%InpJoints) + do i1 = LB(1), UB(1) + call AD_DestroyGSJointType(GSInputFileData%InpJoints(i1), ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + end do + deallocate(GSInputFileData%InpJoints) + end if + if (allocated(GSInputFileData%InpMembers)) then + LB(1:1) = lbound(GSInputFileData%InpMembers) + UB(1:1) = ubound(GSInputFileData%InpMembers) + do i1 = LB(1), UB(1) + call AD_DestroyGSInpMemberType(GSInputFileData%InpMembers(i1), ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + end do + deallocate(GSInputFileData%InpMembers) + end if +end subroutine + +subroutine AD_PackGSInputFile(RF, Indata) + type(RegFile), intent(inout) :: RF + type(GSInputFile), intent(in) :: InData + character(*), parameter :: RoutineName = 'AD_PackGSInputFile' + integer(B4Ki) :: i1 + integer(B4Ki) :: LB(1), UB(1) + if (RF%ErrStat >= AbortErrLev) return + call RegPack(RF, InData%NJoints) + call RegPack(RF, allocated(InData%InpJoints)) + if (allocated(InData%InpJoints)) then + call RegPackBounds(RF, 1, lbound(InData%InpJoints), ubound(InData%InpJoints)) + LB(1:1) = lbound(InData%InpJoints) + UB(1:1) = ubound(InData%InpJoints) + do i1 = LB(1), UB(1) + call AD_PackGSJointType(RF, InData%InpJoints(i1)) + end do + end if + call RegPack(RF, InData%NMembers) + call RegPack(RF, allocated(InData%InpMembers)) + if (allocated(InData%InpMembers)) then + call RegPackBounds(RF, 1, lbound(InData%InpMembers), ubound(InData%InpMembers)) + LB(1:1) = lbound(InData%InpMembers) + UB(1:1) = ubound(InData%InpMembers) + do i1 = LB(1), UB(1) + call AD_PackGSInpMemberType(RF, InData%InpMembers(i1)) + end do + end if + call RegPack(RF, InData%GSPotent) + call RegPack(RF, InData%GSShadow) + call RegPack(RF, InData%GSAero) + if (RegCheckErr(RF, RoutineName)) return +end subroutine + +subroutine AD_UnPackGSInputFile(RF, OutData) + type(RegFile), intent(inout) :: RF + type(GSInputFile), intent(inout) :: OutData + character(*), parameter :: RoutineName = 'AD_UnPackGSInputFile' + integer(B4Ki) :: i1 + integer(B4Ki) :: LB(1), UB(1) + integer(IntKi) :: stat + logical :: IsAllocAssoc + if (RF%ErrStat /= ErrID_None) return + call RegUnpack(RF, OutData%NJoints); if (RegCheckErr(RF, RoutineName)) return + if (allocated(OutData%InpJoints)) deallocate(OutData%InpJoints) + call RegUnpack(RF, IsAllocAssoc); if (RegCheckErr(RF, RoutineName)) return + if (IsAllocAssoc) then + call RegUnpackBounds(RF, 1, LB, UB); if (RegCheckErr(RF, RoutineName)) return + allocate(OutData%InpJoints(LB(1):UB(1)),stat=stat) + if (stat /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating OutData%InpJoints.', RF%ErrStat, RF%ErrMsg, RoutineName) + return + end if + do i1 = LB(1), UB(1) + call AD_UnpackGSJointType(RF, OutData%InpJoints(i1)) ! InpJoints + end do + end if + call RegUnpack(RF, OutData%NMembers); if (RegCheckErr(RF, RoutineName)) return + if (allocated(OutData%InpMembers)) deallocate(OutData%InpMembers) + call RegUnpack(RF, IsAllocAssoc); if (RegCheckErr(RF, RoutineName)) return + if (IsAllocAssoc) then + call RegUnpackBounds(RF, 1, LB, UB); if (RegCheckErr(RF, RoutineName)) return + allocate(OutData%InpMembers(LB(1):UB(1)),stat=stat) + if (stat /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating OutData%InpMembers.', RF%ErrStat, RF%ErrMsg, RoutineName) + return + end if + do i1 = LB(1), UB(1) + call AD_UnpackGSInpMemberType(RF, OutData%InpMembers(i1)) ! InpMembers + end do + end if + call RegUnpack(RF, OutData%GSPotent); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%GSShadow); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%GSAero); if (RegCheckErr(RF, RoutineName)) return +end subroutine + subroutine AD_CopyInputFile(SrcInputFileData, DstInputFileData, CtrlCode, ErrStat, ErrMsg) type(AD_InputFile), intent(in) :: SrcInputFileData type(AD_InputFile), intent(inout) :: DstInputFileData @@ -2263,6 +2726,9 @@ subroutine AD_CopyInputFile(SrcInputFileData, DstInputFileData, CtrlCode, ErrSta if (ErrStat >= AbortErrLev) return end do end if + call AD_CopyGSInputFile(SrcInputFileData%GS, DstInputFileData%GS, CtrlCode, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + if (ErrStat >= AbortErrLev) return end subroutine subroutine AD_DestroyInputFile(InputFileData, ErrStat, ErrMsg) @@ -2299,6 +2765,8 @@ subroutine AD_DestroyInputFile(InputFileData, ErrStat, ErrMsg) end do deallocate(InputFileData%rotors) end if + call AD_DestroyGSInputFile(InputFileData%GS, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) end subroutine subroutine AD_PackInputFile(RF, Indata) @@ -2377,6 +2845,7 @@ subroutine AD_PackInputFile(RF, Indata) call AD_PackRotInputFile(RF, InData%rotors(i1)) end do end if + call AD_PackGSInputFile(RF, InData%GS) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -2462,6 +2931,7 @@ subroutine AD_UnPackInputFile(RF, OutData) call AD_UnpackRotInputFile(RF, OutData%rotors(i1)) ! rotors end do end if + call AD_UnpackGSInputFile(RF, OutData%GS) ! GS end subroutine subroutine AD_CopyRotContinuousStateType(SrcRotContinuousStateTypeData, DstRotContinuousStateTypeData, CtrlCode, ErrStat, ErrMsg) @@ -3321,6 +3791,9 @@ subroutine AD_CopyInflowType(SrcInflowTypeData, DstInflowTypeData, CtrlCode, Err if (ErrStat >= AbortErrLev) return end do end if + call AD_CopyElemInflowType(SrcInflowTypeData%GSInflow, DstInflowTypeData%GSInflow, CtrlCode, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + if (ErrStat >= AbortErrLev) return end subroutine subroutine AD_DestroyInflowType(InflowTypeData, ErrStat, ErrMsg) @@ -3346,6 +3819,8 @@ subroutine AD_DestroyInflowType(InflowTypeData, ErrStat, ErrMsg) end do deallocate(InflowTypeData%RotInflow) end if + call AD_DestroyElemInflowType(InflowTypeData%GSInflow, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) end subroutine subroutine AD_PackInflowType(RF, Indata) @@ -3365,6 +3840,7 @@ subroutine AD_PackInflowType(RF, Indata) call AD_PackRotInflowType(RF, InData%RotInflow(i1)) end do end if + call AD_PackElemInflowType(RF, InData%GSInflow) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -3391,6 +3867,7 @@ subroutine AD_UnPackInflowType(RF, OutData) call AD_UnpackRotInflowType(RF, OutData%RotInflow(i1)) ! RotInflow end do end if + call AD_UnpackElemInflowType(RF, OutData%GSInflow) ! GSInflow end subroutine subroutine AD_CopyRotParameterType(SrcRotParameterTypeData, DstRotParameterTypeData, CtrlCode, ErrStat, ErrMsg) @@ -3815,6 +4292,10 @@ subroutine AD_CopyRotParameterType(SrcRotParameterTypeData, DstRotParameterTypeD call AD_CopyTFinParameterType(SrcRotParameterTypeData%TFin, DstRotParameterTypeData%TFin, CtrlCode, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return + DstRotParameterTypeData%hasGSMod = SrcRotParameterTypeData%hasGSMod + DstRotParameterTypeData%GSPotent = SrcRotParameterTypeData%GSPotent + DstRotParameterTypeData%GSShadow = SrcRotParameterTypeData%GSShadow + DstRotParameterTypeData%GSAero = SrcRotParameterTypeData%GSAero end subroutine subroutine AD_DestroyRotParameterType(RotParameterTypeData, ErrStat, ErrMsg) @@ -4034,6 +4515,10 @@ subroutine AD_PackRotParameterType(RF, Indata) call RegPack(RF, InData%BldNd_NumNodesOut) call RegPack(RF, InData%TFinAero) call AD_PackTFinParameterType(RF, InData%TFin) + call RegPack(RF, InData%hasGSMod) + call RegPack(RF, InData%GSPotent) + call RegPack(RF, InData%GSShadow) + call RegPack(RF, InData%GSAero) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -4146,6 +4631,183 @@ subroutine AD_UnPackRotParameterType(RF, OutData) call RegUnpack(RF, OutData%BldNd_NumNodesOut); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%TFinAero); if (RegCheckErr(RF, RoutineName)) return call AD_UnpackTFinParameterType(RF, OutData%TFin) ! TFin + call RegUnpack(RF, OutData%hasGSMod); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%GSPotent); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%GSShadow); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%GSAero); if (RegCheckErr(RF, RoutineName)) return +end subroutine + +subroutine AD_CopyGSParameterType(SrcGSParameterTypeData, DstGSParameterTypeData, CtrlCode, ErrStat, ErrMsg) + type(GSParameterType), intent(in) :: SrcGSParameterTypeData + type(GSParameterType), intent(inout) :: DstGSParameterTypeData + integer(IntKi), intent(in ) :: CtrlCode + integer(IntKi), intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + integer(B4Ki) :: i1 + integer(B4Ki) :: LB(1), UB(1) + integer(IntKi) :: ErrStat2 + character(ErrMsgLen) :: ErrMsg2 + character(*), parameter :: RoutineName = 'AD_CopyGSParameterType' + ErrStat = ErrID_None + ErrMsg = '' + DstGSParameterTypeData%hasGSMod = SrcGSParameterTypeData%hasGSMod + DstGSParameterTypeData%GSPotent = SrcGSParameterTypeData%GSPotent + DstGSParameterTypeData%GSShadow = SrcGSParameterTypeData%GSShadow + DstGSParameterTypeData%GSAero = SrcGSParameterTypeData%GSAero + DstGSParameterTypeData%NJoints = SrcGSParameterTypeData%NJoints + DstGSParameterTypeData%NNodes = SrcGSParameterTypeData%NNodes + DstGSParameterTypeData%NMembers = SrcGSParameterTypeData%NMembers + if (allocated(SrcGSParameterTypeData%Members)) then + LB(1:1) = lbound(SrcGSParameterTypeData%Members) + UB(1:1) = ubound(SrcGSParameterTypeData%Members) + if (.not. allocated(DstGSParameterTypeData%Members)) then + allocate(DstGSParameterTypeData%Members(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstGSParameterTypeData%Members.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + do i1 = LB(1), UB(1) + call AD_CopyGSMemberType(SrcGSParameterTypeData%Members(i1), DstGSParameterTypeData%Members(i1), CtrlCode, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + if (ErrStat >= AbortErrLev) return + end do + end if + if (allocated(SrcGSParameterTypeData%Joints)) then + LB(1:1) = lbound(SrcGSParameterTypeData%Joints) + UB(1:1) = ubound(SrcGSParameterTypeData%Joints) + if (.not. allocated(DstGSParameterTypeData%Joints)) then + allocate(DstGSParameterTypeData%Joints(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstGSParameterTypeData%Joints.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + do i1 = LB(1), UB(1) + call AD_CopyGSJointType(SrcGSParameterTypeData%Joints(i1), DstGSParameterTypeData%Joints(i1), CtrlCode, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + if (ErrStat >= AbortErrLev) return + end do + end if + DstGSParameterTypeData%Density = SrcGSParameterTypeData%Density + DstGSParameterTypeData%MHK = SrcGSParameterTypeData%MHK + DstGSParameterTypeData%WtrDpth = SrcGSParameterTypeData%WtrDpth +end subroutine + +subroutine AD_DestroyGSParameterType(GSParameterTypeData, ErrStat, ErrMsg) + type(GSParameterType), intent(inout) :: GSParameterTypeData + integer(IntKi), intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + integer(B4Ki) :: i1 + integer(B4Ki) :: LB(1), UB(1) + integer(IntKi) :: ErrStat2 + character(ErrMsgLen) :: ErrMsg2 + character(*), parameter :: RoutineName = 'AD_DestroyGSParameterType' + ErrStat = ErrID_None + ErrMsg = '' + if (allocated(GSParameterTypeData%Members)) then + LB(1:1) = lbound(GSParameterTypeData%Members) + UB(1:1) = ubound(GSParameterTypeData%Members) + do i1 = LB(1), UB(1) + call AD_DestroyGSMemberType(GSParameterTypeData%Members(i1), ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + end do + deallocate(GSParameterTypeData%Members) + end if + if (allocated(GSParameterTypeData%Joints)) then + LB(1:1) = lbound(GSParameterTypeData%Joints) + UB(1:1) = ubound(GSParameterTypeData%Joints) + do i1 = LB(1), UB(1) + call AD_DestroyGSJointType(GSParameterTypeData%Joints(i1), ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + end do + deallocate(GSParameterTypeData%Joints) + end if +end subroutine + +subroutine AD_PackGSParameterType(RF, Indata) + type(RegFile), intent(inout) :: RF + type(GSParameterType), intent(in) :: InData + character(*), parameter :: RoutineName = 'AD_PackGSParameterType' + integer(B4Ki) :: i1 + integer(B4Ki) :: LB(1), UB(1) + if (RF%ErrStat >= AbortErrLev) return + call RegPack(RF, InData%hasGSMod) + call RegPack(RF, InData%GSPotent) + call RegPack(RF, InData%GSShadow) + call RegPack(RF, InData%GSAero) + call RegPack(RF, InData%NJoints) + call RegPack(RF, InData%NNodes) + call RegPack(RF, InData%NMembers) + call RegPack(RF, allocated(InData%Members)) + if (allocated(InData%Members)) then + call RegPackBounds(RF, 1, lbound(InData%Members), ubound(InData%Members)) + LB(1:1) = lbound(InData%Members) + UB(1:1) = ubound(InData%Members) + do i1 = LB(1), UB(1) + call AD_PackGSMemberType(RF, InData%Members(i1)) + end do + end if + call RegPack(RF, allocated(InData%Joints)) + if (allocated(InData%Joints)) then + call RegPackBounds(RF, 1, lbound(InData%Joints), ubound(InData%Joints)) + LB(1:1) = lbound(InData%Joints) + UB(1:1) = ubound(InData%Joints) + do i1 = LB(1), UB(1) + call AD_PackGSJointType(RF, InData%Joints(i1)) + end do + end if + call RegPack(RF, InData%Density) + call RegPack(RF, InData%MHK) + call RegPack(RF, InData%WtrDpth) + if (RegCheckErr(RF, RoutineName)) return +end subroutine + +subroutine AD_UnPackGSParameterType(RF, OutData) + type(RegFile), intent(inout) :: RF + type(GSParameterType), intent(inout) :: OutData + character(*), parameter :: RoutineName = 'AD_UnPackGSParameterType' + integer(B4Ki) :: i1 + integer(B4Ki) :: LB(1), UB(1) + integer(IntKi) :: stat + logical :: IsAllocAssoc + if (RF%ErrStat /= ErrID_None) return + call RegUnpack(RF, OutData%hasGSMod); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%GSPotent); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%GSShadow); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%GSAero); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%NJoints); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%NNodes); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%NMembers); if (RegCheckErr(RF, RoutineName)) return + if (allocated(OutData%Members)) deallocate(OutData%Members) + call RegUnpack(RF, IsAllocAssoc); if (RegCheckErr(RF, RoutineName)) return + if (IsAllocAssoc) then + call RegUnpackBounds(RF, 1, LB, UB); if (RegCheckErr(RF, RoutineName)) return + allocate(OutData%Members(LB(1):UB(1)),stat=stat) + if (stat /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating OutData%Members.', RF%ErrStat, RF%ErrMsg, RoutineName) + return + end if + do i1 = LB(1), UB(1) + call AD_UnpackGSMemberType(RF, OutData%Members(i1)) ! Members + end do + end if + if (allocated(OutData%Joints)) deallocate(OutData%Joints) + call RegUnpack(RF, IsAllocAssoc); if (RegCheckErr(RF, RoutineName)) return + if (IsAllocAssoc) then + call RegUnpackBounds(RF, 1, LB, UB); if (RegCheckErr(RF, RoutineName)) return + allocate(OutData%Joints(LB(1):UB(1)),stat=stat) + if (stat /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating OutData%Joints.', RF%ErrStat, RF%ErrMsg, RoutineName) + return + end if + do i1 = LB(1), UB(1) + call AD_UnpackGSJointType(RF, OutData%Joints(i1)) ! Joints + end do + end if + call RegUnpack(RF, OutData%Density); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%MHK); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%WtrDpth); if (RegCheckErr(RF, RoutineName)) return end subroutine subroutine AD_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) @@ -4210,6 +4872,9 @@ subroutine AD_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) DstParamData%SA_PsiFwd = SrcParamData%SA_PsiFwd DstParamData%SA_nPerSec = SrcParamData%SA_nPerSec DstParamData%WaveField => SrcParamData%WaveField + call AD_CopyGSParameterType(SrcParamData%GS, DstParamData%GS, CtrlCode, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + if (ErrStat >= AbortErrLev) return end subroutine subroutine AD_DestroyParam(ParamData, ErrStat, ErrMsg) @@ -4245,6 +4910,8 @@ subroutine AD_DestroyParam(ParamData, ErrStat, ErrMsg) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) nullify(ParamData%FlowField) nullify(ParamData%WaveField) + call AD_DestroyGSParameterType(ParamData%GS, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) end subroutine subroutine AD_PackParam(RF, Indata) @@ -4300,6 +4967,7 @@ subroutine AD_PackParam(RF, Indata) call SeaSt_WaveField_PackSeaSt_WaveFieldType(RF, InData%WaveField) end if end if + call AD_PackGSParameterType(RF, InData%GS) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -4389,6 +5057,7 @@ subroutine AD_UnPackParam(RF, OutData) else OutData%WaveField => null() end if + call AD_UnpackGSParameterType(RF, OutData%GS) ! GS end subroutine subroutine AD_CopyRotInputType(SrcRotInputTypeData, DstRotInputTypeData, CtrlCode, ErrStat, ErrMsg) @@ -4448,6 +5117,9 @@ subroutine AD_CopyRotInputType(SrcRotInputTypeData, DstRotInputTypeData, CtrlCod call MeshCopy(SrcRotInputTypeData%TFinMotion, DstRotInputTypeData%TFinMotion, CtrlCode, ErrStat2, ErrMsg2 ) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return + call MeshCopy(SrcRotInputTypeData%GSMotion, DstRotInputTypeData%GSMotion, CtrlCode, ErrStat2, ErrMsg2 ) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + if (ErrStat >= AbortErrLev) return if (allocated(SrcRotInputTypeData%UserProp)) then LB(1:2) = lbound(SrcRotInputTypeData%UserProp) UB(1:2) = ubound(SrcRotInputTypeData%UserProp) @@ -4499,6 +5171,8 @@ subroutine AD_DestroyRotInputType(RotInputTypeData, ErrStat, ErrMsg) end if call MeshDestroy( RotInputTypeData%TFinMotion, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + call MeshDestroy( RotInputTypeData%GSMotion, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (allocated(RotInputTypeData%UserProp)) then deallocate(RotInputTypeData%UserProp) end if @@ -4533,6 +5207,7 @@ subroutine AD_PackRotInputType(RF, Indata) end do end if call MeshPack(RF, InData%TFinMotion) + call MeshPack(RF, InData%GSMotion) call RegPackAlloc(RF, InData%UserProp) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -4576,6 +5251,7 @@ subroutine AD_UnPackRotInputType(RF, OutData) end do end if call MeshUnpack(RF, OutData%TFinMotion) ! TFinMotion + call MeshUnpack(RF, OutData%GSMotion) ! GSMotion call RegUnpackAlloc(RF, OutData%UserProp); if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -4716,6 +5392,9 @@ subroutine AD_CopyRotOutputType(SrcRotOutputTypeData, DstRotOutputTypeData, Ctrl call MeshCopy(SrcRotOutputTypeData%TFinLoad, DstRotOutputTypeData%TFinLoad, CtrlCode, ErrStat2, ErrMsg2 ) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return + call MeshCopy(SrcRotOutputTypeData%GSLoad, DstRotOutputTypeData%GSLoad, CtrlCode, ErrStat2, ErrMsg2 ) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + if (ErrStat >= AbortErrLev) return if (allocated(SrcRotOutputTypeData%WriteOutput)) then LB(1:1) = lbound(SrcRotOutputTypeData%WriteOutput) UB(1:1) = ubound(SrcRotOutputTypeData%WriteOutput) @@ -4758,6 +5437,8 @@ subroutine AD_DestroyRotOutputType(RotOutputTypeData, ErrStat, ErrMsg) end if call MeshDestroy( RotOutputTypeData%TFinLoad, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + call MeshDestroy( RotOutputTypeData%GSLoad, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (allocated(RotOutputTypeData%WriteOutput)) then deallocate(RotOutputTypeData%WriteOutput) end if @@ -4783,6 +5464,7 @@ subroutine AD_PackRotOutputType(RF, Indata) end do end if call MeshPack(RF, InData%TFinLoad) + call MeshPack(RF, InData%GSLoad) call RegPackAlloc(RF, InData%WriteOutput) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -4813,6 +5495,7 @@ subroutine AD_UnPackRotOutputType(RF, OutData) end do end if call MeshUnpack(RF, OutData%TFinLoad) ! TFinLoad + call MeshUnpack(RF, OutData%GSLoad) ! GSLoad call RegUnpackAlloc(RF, OutData%WriteOutput); if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -5106,6 +5789,18 @@ subroutine AD_CopyRotMiscVarType(SrcRotMiscVarTypeData, DstRotMiscVarTypeData, C end if DstRotMiscVarTypeData%TwrClrnc = SrcRotMiscVarTypeData%TwrClrnc end if + if (allocated(SrcRotMiscVarTypeData%GSClrnc)) then + LB(1:2) = lbound(SrcRotMiscVarTypeData%GSClrnc) + UB(1:2) = ubound(SrcRotMiscVarTypeData%GSClrnc) + if (.not. allocated(DstRotMiscVarTypeData%GSClrnc)) then + allocate(DstRotMiscVarTypeData%GSClrnc(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstRotMiscVarTypeData%GSClrnc.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstRotMiscVarTypeData%GSClrnc = SrcRotMiscVarTypeData%GSClrnc + end if if (allocated(SrcRotMiscVarTypeData%X)) then LB(1:2) = lbound(SrcRotMiscVarTypeData%X) UB(1:2) = ubound(SrcRotMiscVarTypeData%X) @@ -5637,6 +6332,9 @@ subroutine AD_DestroyRotMiscVarType(RotMiscVarTypeData, ErrStat, ErrMsg) if (allocated(RotMiscVarTypeData%TwrClrnc)) then deallocate(RotMiscVarTypeData%TwrClrnc) end if + if (allocated(RotMiscVarTypeData%GSClrnc)) then + deallocate(RotMiscVarTypeData%GSClrnc) + end if if (allocated(RotMiscVarTypeData%X)) then deallocate(RotMiscVarTypeData%X) end if @@ -5821,6 +6519,7 @@ subroutine AD_PackRotMiscVarType(RF, Indata) call RegPackAlloc(RF, InData%Cant) call RegPackAlloc(RF, InData%Toe) call RegPackAlloc(RF, InData%TwrClrnc) + call RegPackAlloc(RF, InData%GSClrnc) call RegPackAlloc(RF, InData%X) call RegPackAlloc(RF, InData%Y) call RegPackAlloc(RF, InData%Z) @@ -5965,6 +6664,7 @@ subroutine AD_UnPackRotMiscVarType(RF, OutData) call RegUnpackAlloc(RF, OutData%Cant); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%Toe); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%TwrClrnc); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%GSClrnc); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%X); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%Y); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%Z); if (RegCheckErr(RF, RoutineName)) return @@ -6501,6 +7201,10 @@ SUBROUTINE AD_Input_ExtrapInterp1(u1, u2, tin, u_out, tin_out, ErrStat, ErrMsg ) CALL MeshExtrapInterp1(u1%rotors(i01)%TFinMotion, u2%rotors(i01)%TFinMotion, tin, u_out%rotors(i01)%TFinMotion, tin_out, ErrStat2, ErrMsg2) CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg,RoutineName) END DO + do i01 = lbound(u_out%rotors,1),ubound(u_out%rotors,1) + CALL MeshExtrapInterp1(u1%rotors(i01)%GSMotion, u2%rotors(i01)%GSMotion, tin, u_out%rotors(i01)%GSMotion, tin_out, ErrStat2, ErrMsg2) + CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg,RoutineName) + END DO do i01 = lbound(u_out%rotors,1),ubound(u_out%rotors,1) IF (ALLOCATED(u_out%rotors(i01)%UserProp) .AND. ALLOCATED(u1%rotors(i01)%UserProp)) THEN u_out%rotors(i01)%UserProp = a1*u1%rotors(i01)%UserProp + a2*u2%rotors(i01)%UserProp @@ -6599,6 +7303,10 @@ SUBROUTINE AD_Input_ExtrapInterp2(u1, u2, u3, tin, u_out, tin_out, ErrStat, ErrM CALL MeshExtrapInterp2(u1%rotors(i01)%TFinMotion, u2%rotors(i01)%TFinMotion, u3%rotors(i01)%TFinMotion, tin, u_out%rotors(i01)%TFinMotion, tin_out, ErrStat2, ErrMsg2) CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg,RoutineName) END DO + do i01 = lbound(u_out%rotors,1),ubound(u_out%rotors,1) + CALL MeshExtrapInterp2(u1%rotors(i01)%GSMotion, u2%rotors(i01)%GSMotion, u3%rotors(i01)%GSMotion, tin, u_out%rotors(i01)%GSMotion, tin_out, ErrStat2, ErrMsg2) + CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg,RoutineName) + END DO do i01 = lbound(u_out%rotors,1),ubound(u_out%rotors,1) IF (ALLOCATED(u_out%rotors(i01)%UserProp) .AND. ALLOCATED(u1%rotors(i01)%UserProp)) THEN u_out%rotors(i01)%UserProp = a1*u1%rotors(i01)%UserProp + a2*u2%rotors(i01)%UserProp + a3*u3%rotors(i01)%UserProp @@ -6729,6 +7437,10 @@ SUBROUTINE AD_Output_ExtrapInterp1(y1, y2, tin, y_out, tin_out, ErrStat, ErrMsg CALL MeshExtrapInterp1(y1%rotors(i01)%TFinLoad, y2%rotors(i01)%TFinLoad, tin, y_out%rotors(i01)%TFinLoad, tin_out, ErrStat2, ErrMsg2) CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg,RoutineName) END DO + do i01 = lbound(y_out%rotors,1),ubound(y_out%rotors,1) + CALL MeshExtrapInterp1(y1%rotors(i01)%GSLoad, y2%rotors(i01)%GSLoad, tin, y_out%rotors(i01)%GSLoad, tin_out, ErrStat2, ErrMsg2) + CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg,RoutineName) + END DO do i01 = lbound(y_out%rotors,1),ubound(y_out%rotors,1) IF (ALLOCATED(y_out%rotors(i01)%WriteOutput) .AND. ALLOCATED(y1%rotors(i01)%WriteOutput)) THEN y_out%rotors(i01)%WriteOutput = a1*y1%rotors(i01)%WriteOutput + a2*y2%rotors(i01)%WriteOutput @@ -6817,6 +7529,10 @@ SUBROUTINE AD_Output_ExtrapInterp2(y1, y2, y3, tin, y_out, tin_out, ErrStat, Err CALL MeshExtrapInterp2(y1%rotors(i01)%TFinLoad, y2%rotors(i01)%TFinLoad, y3%rotors(i01)%TFinLoad, tin, y_out%rotors(i01)%TFinLoad, tin_out, ErrStat2, ErrMsg2) CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg,RoutineName) END DO + do i01 = lbound(y_out%rotors,1),ubound(y_out%rotors,1) + CALL MeshExtrapInterp2(y1%rotors(i01)%GSLoad, y2%rotors(i01)%GSLoad, y3%rotors(i01)%GSLoad, tin, y_out%rotors(i01)%GSLoad, tin_out, ErrStat2, ErrMsg2) + CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg,RoutineName) + END DO do i01 = lbound(y_out%rotors,1),ubound(y_out%rotors,1) IF (ALLOCATED(y_out%rotors(i01)%WriteOutput) .AND. ALLOCATED(y1%rotors(i01)%WriteOutput)) THEN y_out%rotors(i01)%WriteOutput = a1*y1%rotors(i01)%WriteOutput + a2*y2%rotors(i01)%WriteOutput + a3*y3%rotors(i01)%WriteOutput @@ -6965,6 +7681,12 @@ SUBROUTINE AD_InflowType_ExtrapInterp1(u1, u2, tin, u_out, tin_out, ErrStat, Err u_out%RotInflow(i01)%AvgDiskVel = a1*u1%RotInflow(i01)%AvgDiskVel + a2*u2%RotInflow(i01)%AvgDiskVel END DO END IF ! check if allocated + IF (ALLOCATED(u_out%GSInflow%InflowVel) .AND. ALLOCATED(u1%GSInflow%InflowVel)) THEN + u_out%GSInflow%InflowVel = a1*u1%GSInflow%InflowVel + a2*u2%GSInflow%InflowVel + END IF ! check if allocated + IF (ALLOCATED(u_out%GSInflow%InflowAcc) .AND. ALLOCATED(u1%GSInflow%InflowAcc)) THEN + u_out%GSInflow%InflowAcc = a1*u1%GSInflow%InflowAcc + a2*u2%GSInflow%InflowAcc + END IF ! check if allocated END SUBROUTINE SUBROUTINE AD_InflowType_ExtrapInterp2(u1, u2, u3, tin, u_out, tin_out, ErrStat, ErrMsg ) @@ -7065,6 +7787,12 @@ SUBROUTINE AD_InflowType_ExtrapInterp2(u1, u2, u3, tin, u_out, tin_out, ErrStat, u_out%RotInflow(i01)%AvgDiskVel = a1*u1%RotInflow(i01)%AvgDiskVel + a2*u2%RotInflow(i01)%AvgDiskVel + a3*u3%RotInflow(i01)%AvgDiskVel END DO END IF ! check if allocated + IF (ALLOCATED(u_out%GSInflow%InflowVel) .AND. ALLOCATED(u1%GSInflow%InflowVel)) THEN + u_out%GSInflow%InflowVel = a1*u1%GSInflow%InflowVel + a2*u2%GSInflow%InflowVel + a3*u3%GSInflow%InflowVel + END IF ! check if allocated + IF (ALLOCATED(u_out%GSInflow%InflowAcc) .AND. ALLOCATED(u1%GSInflow%InflowAcc)) THEN + u_out%GSInflow%InflowAcc = a1*u1%GSInflow%InflowAcc + a2*u2%GSInflow%InflowAcc + a3*u3%GSInflow%InflowAcc + END IF ! check if allocated END SUBROUTINE function AD_InputMeshPointer(u, DL) result(Mesh) @@ -7085,6 +7813,8 @@ function AD_InputMeshPointer(u, DL) result(Mesh) Mesh => u%BladeMotion(DL%i1) case (AD_u_TFinMotion) Mesh => u%TFinMotion + case (AD_u_GSMotion) + Mesh => u%GSMotion end select end function @@ -7104,6 +7834,8 @@ function AD_OutputMeshPointer(y, DL) result(Mesh) Mesh => y%BladeLoad(DL%i1) case (AD_y_TFinLoad) Mesh => y%TFinLoad + case (AD_y_GSLoad) + Mesh => y%GSLoad end select end function @@ -7240,6 +7972,8 @@ subroutine AD_VarPackInput(V, u, ValAry) call MV_PackMesh(V, u%BladeMotion(DL%i1), ValAry) ! Mesh case (AD_u_TFinMotion) call MV_PackMesh(V, u%TFinMotion, ValAry) ! Mesh + case (AD_u_GSMotion) + call MV_PackMesh(V, u%GSMotion, ValAry) ! Mesh case (AD_u_UserProp) VarVals = u%UserProp(V%iLB:V%iUB,V%j) ! Rank 2 Array case default @@ -7276,6 +8010,8 @@ subroutine AD_VarUnpackInput(V, ValAry, u) call MV_UnpackMesh(V, ValAry, u%BladeMotion(DL%i1)) ! Mesh case (AD_u_TFinMotion) call MV_UnpackMesh(V, ValAry, u%TFinMotion) ! Mesh + case (AD_u_GSMotion) + call MV_UnpackMesh(V, ValAry, u%GSMotion) ! Mesh case (AD_u_UserProp) u%UserProp(V%iLB:V%iUB, V%j) = VarVals ! Rank 2 Array end select @@ -7298,6 +8034,8 @@ function AD_InputFieldName(DL) result(Name) Name = "u%BladeMotion("//trim(Num2LStr(DL%i1))//")" case (AD_u_TFinMotion) Name = "u%TFinMotion" + case (AD_u_GSMotion) + Name = "u%GSMotion" case (AD_u_UserProp) Name = "u%UserProp" case default @@ -7331,6 +8069,8 @@ subroutine AD_VarPackOutput(V, y, ValAry) call MV_PackMesh(V, y%BladeLoad(DL%i1), ValAry) ! Mesh case (AD_y_TFinLoad) call MV_PackMesh(V, y%TFinLoad, ValAry) ! Mesh + case (AD_y_GSLoad) + call MV_PackMesh(V, y%GSLoad, ValAry) ! Mesh case (AD_y_WriteOutput) VarVals = y%WriteOutput(V%iLB:V%iUB) ! Rank 1 Array case default @@ -7365,6 +8105,8 @@ subroutine AD_VarUnpackOutput(V, ValAry, y) call MV_UnpackMesh(V, ValAry, y%BladeLoad(DL%i1)) ! Mesh case (AD_y_TFinLoad) call MV_UnpackMesh(V, ValAry, y%TFinLoad) ! Mesh + case (AD_y_GSLoad) + call MV_UnpackMesh(V, ValAry, y%GSLoad) ! Mesh case (AD_y_WriteOutput) y%WriteOutput(V%iLB:V%iUB) = VarVals ! Rank 1 Array end select @@ -7385,6 +8127,8 @@ function AD_OutputFieldName(DL) result(Name) Name = "y%BladeLoad("//trim(Num2LStr(DL%i1))//")" case (AD_y_TFinLoad) Name = "y%TFinLoad" + case (AD_y_GSLoad) + Name = "y%GSLoad" case (AD_y_WriteOutput) Name = "y%WriteOutput" case default diff --git a/modules/externalinflow/src/ExternalInflow.f90 b/modules/externalinflow/src/ExternalInflow.f90 index c12b324caa..1c49f2b6c7 100644 --- a/modules/externalinflow/src/ExternalInflow.f90 +++ b/modules/externalinflow/src/ExternalInflow.f90 @@ -90,6 +90,17 @@ SUBROUTINE Init_ExtInfw( InitInp, p_FAST, AirDens, u_AD, initOut_AD, y_AD, ExtIn call SetErrStat(ErrID_Warn, ErrMsg2, ErrStat, ErrMsg, RoutineName) endif + ! ExternalInflow does not export the AeroDyn generalized-support-structure (GS) node positions + ! to CFD, nor map velocities back for them. GS nodes are appended to AeroDyn's wind-sampling + ! list, so leaving them uncounted here would silently return zero inflow on the GS (and shift + ! any OLAF wake-point indices). Reject the combination rather than run with bad inflow. + if (u_AD%rotors(1)%GSMotion%NNodes > 0) then + call SetErrStat(ErrID_Fatal, 'The AeroDyn generalized support structure (GSPotent/GSShadow/GSAero) '// & + 'is not supported with ExternalInflow (CompInflow=2). Disable the GS model to use this inflow option.', & + ErrStat, ErrMsg, RoutineName) + return + end if + !--------------------------- ! Motion points from AD15 !--------------------------- diff --git a/modules/openfast-library/src/FAST_Funcs.f90 b/modules/openfast-library/src/FAST_Funcs.f90 index ff46ce942b..995ba6fb7a 100644 --- a/modules/openfast-library/src/FAST_Funcs.f90 +++ b/modules/openfast-library/src/FAST_Funcs.f90 @@ -870,7 +870,7 @@ subroutine FAST_CalcOutput(ModData, Mappings, ThisTime, iInput, iState, T, ErrSt subroutine FAST_GetOP(ModData, ThisTime, iInput, iState, T, ErrStat, ErrMsg, & u_op, y_op, x_op, dx_op, z_op, u_glue, y_glue, x_glue, dx_glue, z_glue) - use AeroDyn, only: AD_CalcWind_Rotor + use AeroDyn, only: AD_CalcWind_Rotor, AD_CalcWind_GS type(ModDataType), intent(in) :: ModData !< Module information real(DbKi), intent(in) :: ThisTime !< Time integer(IntKi), intent(in) :: iInput !< Input index @@ -1086,8 +1086,14 @@ subroutine FAST_GetOP(ModData, ThisTime, iInput, iState, T, ErrStat, ErrMsg, & T%AD%m, T%AD%m%Inflow(iInput)%RotInflow(ModData%Ins), & i, ErrStat2, ErrMsg2) if (Failed()) return + call AD_CalcWind_GS(ThisTime, T%AD%Input(iInput)%rotors(ModData%Ins), & + T%AD%p%FlowField, T%AD%p%GS, T%AD%p, & + T%AD%m, T%AD%m%Inflow(iInput)%GSInflow, & + i, ErrStat2, ErrMsg2) + if (Failed()) return call RotCalcContStateDeriv(ThisTime, T%AD%Input(iInput)%rotors(ModData%Ins), & T%AD%m%Inflow(iInput)%RotInflow(ModData%Ins), & + T%AD%m%Inflow(iInput)%GSInflow, & T%AD%p%rotors(ModData%Ins), T%AD%p, & T%AD%x(iState)%rotors(ModData%Ins), & T%AD%xd(iState)%rotors(ModData%Ins), & diff --git a/modules/openfast-library/src/FAST_Mapping.f90 b/modules/openfast-library/src/FAST_Mapping.f90 index 0f4444c139..f11f3202c0 100644 --- a/modules/openfast-library/src/FAST_Mapping.f90 +++ b/modules/openfast-library/src/FAST_Mapping.f90 @@ -706,6 +706,14 @@ subroutine InitMappings_AD(Mappings, SrcMod, DstMod, Turbine, ErrStat, ErrMsg) Active=NotCompAeroMaps) if (Failed()) return + associate (AD_p_Rotor => Turbine%AD%p%Rotors(DstMod%Ins)) + call MapMotionMesh(Turbine, Mappings, & + SrcMod=SrcMod, SrcDL=DatLoc(ED_y_PlatformPtMesh), & ! ED%y%PlatformPtMesh + DstMod=DstMod, DstDL=DatLoc(AD_u_GSMotion), & ! AD%u%rotors(DstMod%Ins)%GSMotion + ErrStat=ErrStat2, ErrMsg=ErrMsg2, & + Active=(AD_p_Rotor%hasGSMod.and.(Turbine%p_FAST%CompSub/=Module_SD).and.NotCompAeroMaps)) + if (Failed()) return + end associate case (Module_SED) @@ -781,6 +789,15 @@ subroutine InitMappings_AD(Mappings, SrcMod, DstMod, Turbine, ErrStat, ErrMsg) call MapCustom(Mappings, Custom_SrvD_to_AD, SrcMod, DstMod) + case (Module_SD) + associate (AD_p_Rotor => Turbine%AD%p%Rotors(DstMod%Ins)) + call MapMotionMesh(Turbine, Mappings, & + SrcMod=SrcMod, SrcDL=DatLoc(SD_y_Y3Mesh), & ! SD%y(SrcMod%Ins)%Y3Mesh + DstMod=DstMod, DstDL=DatLoc(AD_u_GSMotion), & ! AD%u%rotors(DstMod%Ins)%GSMotion + ErrStat=ErrStat2, ErrMsg=ErrMsg2, & + Active=(AD_p_Rotor%hasGSMod)) + if (Failed()) return + end associate end select contains @@ -1012,6 +1029,18 @@ subroutine InitMappings_ED(Mappings, SrcMod, DstMod, Turbine, ErrStat, ErrMsg) Active=CompAeroAD .and. NotCompAeroMaps) if (Failed()) return + ! Generalized support structure + associate (AD_p_Rotor => Turbine%AD%p%Rotors(SrcMod%Ins)) + call MapLoadMesh(Turbine, Mappings, SrcMod=SrcMod, DstMod=DstMod, & + SrcDL=DatLoc(AD_y_GSLoad), & ! AD%y%rotors(SrcMod%Ins)%GSLoad + SrcDispDL=DatLoc(AD_u_GSMotion),& ! AD%u%rotors(SrcMod%Ins)%GSMotion + DstDL=DatLoc(ED_u_PlatformPtMesh), & ! ED%u%PlatformPtMesh + DstDispDL=DatLoc(ED_y_PlatformPtMesh),& ! ED%y%PlatformPtMesh + ErrStat=ErrStat2, ErrMsg=ErrMsg2, & + Active=((AD_p_Rotor%GSAero/=GSAero_none).and.(Turbine%p_FAST%CompSub/=Module_SD).and.NotCompAeroMaps)) + if (Failed()) return + end associate + case (Module_ADsk) ! Hub Loads @@ -1955,6 +1984,18 @@ subroutine InitMappings_SD(Mappings, SrcMod, DstMod, Turbine, ErrStat, ErrMsg) select case (SrcMod%ID) + case (Module_AD) + associate (AD_p_Rotor => Turbine%AD%p%Rotors(SrcMod%Ins)) + call MapLoadMesh(Turbine, Mappings, SrcMod=SrcMod, DstMod=DstMod, & + SrcDL=DatLoc(AD_y_GSLoad), & ! AD%y%rotors(SrcMod%Ins)%GSLoad + SrcDispDL=DatLoc(AD_u_GSMotion),& ! AD%u%rotors(SrcMod%Ins)%GSMotion + DstDL=DatLoc(SD_u_LMesh), & ! SD%u(DstMod%Ins)%LMesh + DstDispDL=DatLoc(SD_y_y3Mesh),& ! SD%y(DstMod%Ins)%y3Mesh + ErrStat=ErrStat2, ErrMsg=ErrMsg2, & + Active=AD_p_Rotor%GSAero/=GSAero_none) + if (Failed()) return + end associate + case (Module_ED) call MapMotionMesh(Turbine, Mappings, SrcMod=SrcMod, DstMod=DstMod, & diff --git a/modules/openfast-library/src/FAST_Subs.f90 b/modules/openfast-library/src/FAST_Subs.f90 index 1567849764..08507d8b68 100644 --- a/modules/openfast-library/src/FAST_Subs.f90 +++ b/modules/openfast-library/src/FAST_Subs.f90 @@ -6006,6 +6006,33 @@ SUBROUTINE WrVTK_AllMeshes(p_FAST, y_FAST, ED, SED, BD, AD, IfW, ExtInfw, HD, SD end do end if end do + + ! General support structure (GS): a single shared structure across all rotors, so + ! write at most one copy (no rotor suffix) to avoid duplicates in multirotor cases. + ! Prefer the rotor with a committed GSLoad (GSAero on => drag vectors), else fall + ! back to a committed GSMotion (influence-only, points). + j = 0 + do iRot = 1, p_FAST%NRotors + if (AD%y%rotors(iRot)%GSLoad%Committed) then + j = iRot + exit + end if + end do + if (j>0) then + call MeshWrVTK(p_FAST%TurbinePos, AD%y%rotors(j)%GSLoad, & + trim(p_FAST%VTK_OutFileRoot)//'.AD_GS', & + y_FAST%VTK_count, p_FAST%VTK_fields, ErrStat2, ErrMsg2, p_FAST%VTK_tWidth, & + AD%Input(INPUT_CURR)%rotors(j)%GSMotion ) + else + do iRot = 1, p_FAST%NRotors + if (AD%Input(INPUT_CURR)%rotors(iRot)%GSMotion%Committed) then + call MeshWrVTK(p_FAST%TurbinePos, AD%Input(INPUT_CURR)%rotors(iRot)%GSMotion, & + trim(p_FAST%VTK_OutFileRoot)//'.AD_GS', & + y_FAST%VTK_count, p_FAST%VTK_fields, ErrStat2, ErrMsg2, p_FAST%VTK_tWidth ) + exit + end if + end do + end if end if ! FVW submodule of AD15 diff --git a/modules/subdyn/src/SubDyn_Types.f90 b/modules/subdyn/src/SubDyn_Types.f90 index 9fee52b5f2..5bdaa7d288 100644 --- a/modules/subdyn/src/SubDyn_Types.f90 +++ b/modules/subdyn/src/SubDyn_Types.f90 @@ -170,7 +170,7 @@ MODULE SubDyn_Types REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: MemberSpin !< Member spin angle about its axis - for rectangular members [rad] REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: MemberDivSize !< Optional maximum element length for each member [m] INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: MemberNDiv !< Resolved number of finite elements per member [-] - INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: MemberElemStart !< First element index for each member in p%Elems [-] + INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: MemberElemStart !< First element index for each member in p%Elems [-] CHARACTER(ChanLen) , DIMENSION(:), ALLOCATABLE :: SSOutList !< List of Output Channels [-] LOGICAL :: OutCOSM = .false. !< Output Cos-matrices Flag [-] LOGICAL :: TabDelim = .false. !< Generate a tab-delimited output file in OutJckF-Flag [-] @@ -196,7 +196,7 @@ MODULE SubDyn_Types REAL(R8Ki) , DIMENSION(:,:), ALLOCATABLE :: K !< System stiffness matrix [-] REAL(R8Ki) , DIMENSION(:,:), ALLOCATABLE :: M !< System mass matrix [-] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: ElemProps !< Element properties(A, L, Ixx, Iyy, Jzz, Shear, Kappa, E, G, Rho, DirCos(1,1), DirCos(2, 1), ....., DirCos(3, 3) ) [-] - INTEGER(IntKi) , DIMENSION(:,:), ALLOCATABLE :: MemberNodes !< Member number and list of nodes making up a member (>2 if subdivided) [-] + INTEGER(IntKi) , DIMENSION(:,:), ALLOCATABLE :: MemberNodes !< Member number and endpoint node IDs (interior member nodes reconstructed from connectivity) [-] INTEGER(IntKi) , DIMENSION(:,:), ALLOCATABLE :: NodesConnN !< Nodes that connect to a common node [-] INTEGER(IntKi) , DIMENSION(:,:), ALLOCATABLE :: NodesConnE !< Elements that connect to a common node [-] LOGICAL :: SSSum = .false. !< SubDyn Summary File Flag [-] diff --git a/openfast_io/openfast_io/FAST_reader.py b/openfast_io/openfast_io/FAST_reader.py index 6c6c9e351a..98a85340d5 100644 --- a/openfast_io/openfast_io/FAST_reader.py +++ b/openfast_io/openfast_io/FAST_reader.py @@ -1070,6 +1070,9 @@ def read_AeroDyn(self): self.fst_vt['AeroDyn']['TwrPotent'] = int(f.readline().split()[0]) self.fst_vt['AeroDyn']['TwrShadow'] = int(f.readline().split()[0]) self.fst_vt['AeroDyn']['TwrAero'] = bool_read(f.readline().split()[0]) + self.fst_vt['AeroDyn']['GSPotent'] = int(f.readline().split()[0]) + self.fst_vt['AeroDyn']['GSShadow'] = int(f.readline().split()[0]) + self.fst_vt['AeroDyn']['GSAero'] = bool_read(f.readline().split()[0]) self.fst_vt['AeroDyn']['CavitCheck'] = bool_read(f.readline().split()[0]) self.fst_vt['AeroDyn']['NacelleDrag'] = bool_read(f.readline().split()[0]) self.fst_vt['AeroDyn']['CompAA'] = bool_read(f.readline().split()[0]) @@ -1201,6 +1204,50 @@ def read_AeroDyn(self): self.fst_vt['AeroDyn']['TwrCp'][i] = data[5] self.fst_vt['AeroDyn']['TwrCa'][i] = data[6] + # Generalized support structure joints + f.readline() + self.fst_vt['AeroDyn']['NumGSJoints'] = int(f.readline().split()[0]) + f.readline() + f.readline() + self.fst_vt['AeroDyn']['GSJointID'] = [None]*self.fst_vt['AeroDyn']['NumGSJoints'] + self.fst_vt['AeroDyn']['GSJointXi'] = [None]*self.fst_vt['AeroDyn']['NumGSJoints'] + self.fst_vt['AeroDyn']['GSJointYi'] = [None]*self.fst_vt['AeroDyn']['NumGSJoints'] + self.fst_vt['AeroDyn']['GSJointZi'] = [None]*self.fst_vt['AeroDyn']['NumGSJoints'] + for i in range(self.fst_vt['AeroDyn']['NumGSJoints']): + data = [float(val) for val in f.readline().split()] + self.fst_vt['AeroDyn']['GSJointID'][i] = int(data[0]) + self.fst_vt['AeroDyn']['GSJointXi'][i] = data[1] + self.fst_vt['AeroDyn']['GSJointYi'][i] = data[2] + self.fst_vt['AeroDyn']['GSJointZi'][i] = data[3] + + # Generalized support structure members + f.readline() + self.fst_vt['AeroDyn']['NumGSMembers'] = int(f.readline().split()[0]) + f.readline() + f.readline() + self.fst_vt['AeroDyn']['GSMemberID'] = [None]*self.fst_vt['AeroDyn']['NumGSMembers'] + self.fst_vt['AeroDyn']['GSMJointID1'] = [None]*self.fst_vt['AeroDyn']['NumGSMembers'] + self.fst_vt['AeroDyn']['GSMJointID2'] = [None]*self.fst_vt['AeroDyn']['NumGSMembers'] + self.fst_vt['AeroDyn']['GSMDia1'] = [None]*self.fst_vt['AeroDyn']['NumGSMembers'] + self.fst_vt['AeroDyn']['GSMDia2'] = [None]*self.fst_vt['AeroDyn']['NumGSMembers'] + self.fst_vt['AeroDyn']['GSMCd1'] = [None]*self.fst_vt['AeroDyn']['NumGSMembers'] + self.fst_vt['AeroDyn']['GSMCd2'] = [None]*self.fst_vt['AeroDyn']['NumGSMembers'] + self.fst_vt['AeroDyn']['GSMTI1'] = [None]*self.fst_vt['AeroDyn']['NumGSMembers'] + self.fst_vt['AeroDyn']['GSMTI2'] = [None]*self.fst_vt['AeroDyn']['NumGSMembers'] + self.fst_vt['AeroDyn']['GSMDiv'] = [None]*self.fst_vt['AeroDyn']['NumGSMembers'] + for i in range(self.fst_vt['AeroDyn']['NumGSMembers']): + data = [float(val) for val in f.readline().split()] + self.fst_vt['AeroDyn']['GSMemberID'][i] = int(data[0]) + self.fst_vt['AeroDyn']['GSMJointID1'][i] = int(data[1]) + self.fst_vt['AeroDyn']['GSMJointID2'][i] = int(data[2]) + self.fst_vt['AeroDyn']['GSMDia1'][i] = data[3] + self.fst_vt['AeroDyn']['GSMDia2'][i] = data[4] + self.fst_vt['AeroDyn']['GSMCd1'][i] = data[5] + self.fst_vt['AeroDyn']['GSMCd2'][i] = data[6] + self.fst_vt['AeroDyn']['GSMTI1'][i] = data[7] + self.fst_vt['AeroDyn']['GSMTI2'][i] = data[8] + self.fst_vt['AeroDyn']['GSMDiv'][i] = data[9] + # Outputs f.readline() self.fst_vt['AeroDyn']['SumPrint'] = bool_read(f.readline().split()[0]) diff --git a/openfast_io/openfast_io/FAST_writer.py b/openfast_io/openfast_io/FAST_writer.py index d97d1304b0..e6342ee097 100644 --- a/openfast_io/openfast_io/FAST_writer.py +++ b/openfast_io/openfast_io/FAST_writer.py @@ -979,6 +979,9 @@ def write_AeroDyn(self): f.write('{:<22d} {:<11} {:}'.format(self.fst_vt['AeroDyn']['TwrPotent'], 'TwrPotent', '- Type tower influence on wind based on potential flow around the tower (switch) {0=none, 1=baseline potential flow, 2=potential flow with Bak correction}\n')) f.write('{:<22d} {:<11} {:}'.format(self.fst_vt['AeroDyn']['TwrShadow'], 'TwrShadow', '- Calculate tower influence on wind based on downstream tower shadow (switch) {0=none, 1=Powles model, 2=Eames model}\n')) f.write('{!s:<22} {:<11} {:}'.format(self.fst_vt['AeroDyn']['TwrAero'], 'TwrAero', '- Calculate tower aerodynamic loads? (flag)\n')) + f.write('{:<22d} {:<11} {:}'.format(self.fst_vt['AeroDyn']['GSPotent'], 'GSPotent', '- Type of general support structure influence on wind based on potential flow around the member (switch) {0=none, 1=baseline potential flow, 2=potential flow with Bak correction}\n')) + f.write('{:<22d} {:<11} {:}'.format(self.fst_vt['AeroDyn']['GSShadow'], 'GSShadow', '- Calculate general support structure influence on wind based on downstream shadow (switch) {0=none, 1=Powles model, 2=Eames model}\n')) + f.write('{!s:<22} {:<11} {:}'.format(self.fst_vt['AeroDyn']['GSAero'], 'GSAero', '- Calculate general support structure aerodynamic loads? (flag)\n')) f.write('{!s:<22} {:<11} {:}'.format(self.fst_vt['AeroDyn']['CavitCheck'], 'CavitCheck', '- Perform cavitation check? (flag) [UA_Mod must be 0 when CavitCheck=true]\n')) f.write('{!s:<22} {:<11} {:}'.format(self.fst_vt['AeroDyn']['NacelleDrag'], 'NacelleDrag', '- Include Nacelle Drag effects? (flag)\n')) f.write('{!s:<22} {:<11} {:}'.format(self.fst_vt['AeroDyn']['CompAA'], 'CompAA', '- Flag to compute AeroAcoustics calculation [used only when Wake_Mod = 1 or 2]\n')) @@ -1061,6 +1064,18 @@ def write_AeroDyn(self): f.write('(m) (m) (-) (-) (-) (-) (-)\n') for TwrElev, TwrDiam, TwrCd, TwrTI, TwrCb, TwrCp, TwrCa in zip(self.fst_vt['AeroDyn']['TwrElev'], self.fst_vt['AeroDyn']['TwrDiam'], self.fst_vt['AeroDyn']['TwrCd'], self.fst_vt['AeroDyn']['TwrTI'], self.fst_vt['AeroDyn']['TwrCb'], self.fst_vt['AeroDyn']['TwrCp'], self.fst_vt['AeroDyn']['TwrCa']): f.write('{: 2.15e} {: 2.15e} {: 2.15e} {: 2.15e} {: 2.15e} {: 2.15e} {: 2.15e} \n'.format(TwrElev, TwrDiam, TwrCd, TwrTI, TwrCb, TwrCp, TwrCa)) + f.write('====== Generalized Support Structure Joints ========================================================\n') + f.write('{:<22d} {:<11} {:}'.format(self.fst_vt['AeroDyn']['NumGSJoints'], 'NumGSJoints', '- Number of general support structure joints used in the analysis (-)\n')) + f.write('GSJointID GSJointXi GSJointYi GSJointZi\n') + f.write('(-) (m) (m) (m)\n') + for GSJointID, GSJointXi, GSJointYi, GSJointZi in zip(self.fst_vt['AeroDyn']['GSJointID'], self.fst_vt['AeroDyn']['GSJointXi'], self.fst_vt['AeroDyn']['GSJointYi'], self.fst_vt['AeroDyn']['GSJointZi']): + f.write('{:<14d} {: 2.15e} {: 2.15e} {: 2.15e}\n'.format(int(GSJointID), GSJointXi, GSJointYi, GSJointZi)) + f.write('====== Generalized Support Structure Members =======================================================\n') + f.write('{:<22d} {:<11} {:}'.format(self.fst_vt['AeroDyn']['NumGSMembers'], 'NumGSMembers', '- Number of general support structure members used in the analysis (-)\n')) + f.write('GSMemberID GSMJointID1 GSMJointID2 GSMDia1 GSMDia2 GSMCd1 GSMCd2 GSMTI1 GSMTI2 GSMDiv\n') + f.write('(-) (-) (-) (m) (m) (-) (-) (-) (-) (m)\n') + for GSMemberID, GSMJointID1, GSMJointID2, GSMDia1, GSMDia2, GSMCd1, GSMCd2, GSMTI1, GSMTI2, GSMDiv in zip(self.fst_vt['AeroDyn']['GSMemberID'], self.fst_vt['AeroDyn']['GSMJointID1'], self.fst_vt['AeroDyn']['GSMJointID2'], self.fst_vt['AeroDyn']['GSMDia1'], self.fst_vt['AeroDyn']['GSMDia2'], self.fst_vt['AeroDyn']['GSMCd1'], self.fst_vt['AeroDyn']['GSMCd2'], self.fst_vt['AeroDyn']['GSMTI1'], self.fst_vt['AeroDyn']['GSMTI2'], self.fst_vt['AeroDyn']['GSMDiv']): + f.write('{:<14d} {:<15d} {:<15d} {: 2.15e} {: 2.15e} {: 2.15e} {: 2.15e} {: 2.15e} {: 2.15e} {: 2.15e}\n'.format(int(GSMemberID), int(GSMJointID1), int(GSMJointID2), GSMDia1, GSMDia2, GSMCd1, GSMCd2, GSMTI1, GSMTI2, GSMDiv)) f.write('====== Outputs ====================================================================================\n') f.write('{!s:<22} {:<11} {:}'.format(self.fst_vt['AeroDyn']['SumPrint'], 'SumPrint', '- Generate a summary file listing input options and interpolated properties to ".AD.sum"? (flag)\n')) f.write('{:<22d} {:<11} {:}'.format(self.fst_vt['AeroDyn']['NBlOuts'], 'NBlOuts', '- Number of blade node outputs [0 - 9] (-)\n')) diff --git a/reg_tests/r-test b/reg_tests/r-test index 174a8ef776..8aefcb0e74 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit 174a8ef7764a845782932122ab6ec5992c0f72e9 +Subproject commit 8aefcb0e74b711bc7eebf77fab246d2f9bb7c817