From 18d85d7e9157e6a6e97d50ff123355db954a0514 Mon Sep 17 00:00:00 2001 From: Hannah Ross Date: Mon, 28 Sep 2026 12:39:52 -0600 Subject: [PATCH 01/12] Add UA_Mod=9 (IAG) registry definitions: new model enum, Ka/Kv/dCNdA/CnMax/CnMin airfoil parameters, IAG OtherState vortex variables, and updated continuous-state documentation --- modules/aerodyn/src/AirfoilInfo_Registry.txt | 11 +++ modules/aerodyn/src/AirfoilInfo_Types.f90 | 43 +++++++++++ modules/aerodyn/src/UnsteadyAero.f90 | 6 +- modules/aerodyn/src/UnsteadyAero_Registry.txt | 6 +- modules/aerodyn/src/UnsteadyAero_Types.f90 | 74 ++++++++++++++++++- 5 files changed, 136 insertions(+), 4 deletions(-) diff --git a/modules/aerodyn/src/AirfoilInfo_Registry.txt b/modules/aerodyn/src/AirfoilInfo_Registry.txt index c4fff142f8..690a5595b4 100644 --- a/modules/aerodyn/src/AirfoilInfo_Registry.txt +++ b/modules/aerodyn/src/AirfoilInfo_Registry.txt @@ -31,6 +31,7 @@ param AirfoilInfo/AFI - INTEGER UA_HGMV param AirfoilInfo/AFI - INTEGER UA_Oye - 6 - "Stieg Oye dynamic stall model" - param AirfoilInfo/AFI - INTEGER UA_BV - 7 - "Boeing-Vertol dynamic stall model (e.g. used in CACTUS)" - param AirfoilInfo/AFI - INTEGER UA_HGMV360 - 8 - "continuous variant of HGM (Hansen) model with vortex modifications modified for 360-deg" - +param AirfoilInfo/AFI - INTEGER UA_IAG - 9 - "IAG dynamic stall model (first-order, state-space)" - # ..... Airfoil data ............................................................................................................... @@ -82,6 +83,12 @@ typedef ^ ^ ReKi CnBreakUppe typedef ^ ^ ReKi alphaBreakLower - - 2pi "(calculated) Angle of attack where normal and reverse flow CnAttached intersect; between -pi and 0; will be near -pi/2 deg in most cases" "rad" typedef ^ ^ ReKi CnBreakLower - - - "(calculated) CnAttached value at alphaBreakLower where normal and reverse flow CnAttached intersect; will be negative" "-" +typedef ^ ^ ReKi Ka - - - "IAG model: impulsive (non-circulatory) normal-force gain [default=0.75]" - +typedef ^ ^ ReKi Kv - - - "IAG model: center-of-pressure amplitude for vortex moment [default=0.2]" - +typedef ^ ^ ReKi dCNdA - - - "IAG model: static normal-force curve slope from linear-fit approach" 1/rad +typedef ^ ^ ReKi CnMax - - - "IAG model: (calculated) maximum static Cn, used as positive CN_CRIT" - +typedef ^ ^ ReKi CnMin - - - "IAG model: (calculated) minimum static Cn, used as negative CN_CRIT" - + typedef AirfoilInfo/AFI AFI_UA_BL_Default_Type LOGICAL alpha0 - .true. - "Calculate value for this input?" - typedef ^ ^ LOGICAL alpha1 - .true. - "Calculate value for this input?" - typedef ^ ^ LOGICAL alpha2 - .true. - "Calculate value for this input?" - @@ -120,6 +127,10 @@ typedef ^ ^ LOGICAL filtCutOff typedef ^ ^ LOGICAL alphaUpper - .true. - "Calculate value for this input?" - typedef ^ ^ LOGICAL alphaLower - .true. - "Calculate value for this input?" - +typedef ^ ^ LOGICAL Ka - .true. - "Calculate value for this input?" - +typedef ^ ^ LOGICAL Kv - .true. - "Calculate value for this input?" - +typedef ^ ^ LOGICAL dCNdA - .true. - "Calculate value for this input?" - + # The following derived type stores data for an airfoil at a single combination of Re and control setting. typedef ^ AFI_Table_Type ReKi Alpha {:} - - "Angle-of-attack vector that matches the Coefs matrix" rad diff --git a/modules/aerodyn/src/AirfoilInfo_Types.f90 b/modules/aerodyn/src/AirfoilInfo_Types.f90 index ae02f958b3..daac5fe418 100644 --- a/modules/aerodyn/src/AirfoilInfo_Types.f90 +++ b/modules/aerodyn/src/AirfoilInfo_Types.f90 @@ -45,6 +45,7 @@ MODULE AirfoilInfo_Types INTEGER(IntKi), PUBLIC, PARAMETER :: UA_Oye = 6 ! Stieg Oye dynamic stall model [-] INTEGER(IntKi), PUBLIC, PARAMETER :: UA_BV = 7 ! Boeing-Vertol dynamic stall model (e.g. used in CACTUS) [-] INTEGER(IntKi), PUBLIC, PARAMETER :: UA_HGMV360 = 8 ! continuous variant of HGM (Hansen) model with vortex modifications modified for 360-deg [-] + INTEGER(IntKi), PUBLIC, PARAMETER :: UA_IAG = 9 ! IAG dynamic stall model (first-order, state-space) [-] ! ========= AFI_UA_BL_Type ======= TYPE, PUBLIC :: AFI_UA_BL_Type REAL(ReKi) :: alpha0 = 0.0_ReKi !< Angle of attack for zero lift (also used in HGM) [input in degrees; stored as radians] @@ -91,6 +92,11 @@ MODULE AirfoilInfo_Types REAL(ReKi) :: CnBreakUpper = 0.0_ReKi !< (calculated) CnAttached value at alphaBreakUpper where normal and reverse flow CnAttached intersect; will be positive [-] REAL(ReKi) :: alphaBreakLower = 0.0_ReKi !< (calculated) Angle of attack where normal and reverse flow CnAttached intersect; between -pi and 0; will be near -pi/2 deg in most cases [rad] REAL(ReKi) :: CnBreakLower = 0.0_ReKi !< (calculated) CnAttached value at alphaBreakLower where normal and reverse flow CnAttached intersect; will be negative [-] + REAL(ReKi) :: Ka = 0.0_ReKi !< IAG model: impulsive (non-circulatory) normal-force gain [default=0.75] [-] + REAL(ReKi) :: Kv = 0.0_ReKi !< IAG model: center-of-pressure amplitude for vortex moment [default=0.2] [-] + REAL(ReKi) :: dCNdA = 0.0_ReKi !< IAG model: static normal-force curve slope from linear-fit approach [1/rad] + REAL(ReKi) :: CnMax = 0.0_ReKi !< IAG model: (calculated) maximum static Cn, used as positive CN_CRIT [-] + REAL(ReKi) :: CnMin = 0.0_ReKi !< IAG model: (calculated) minimum static Cn, used as negative CN_CRIT [-] END TYPE AFI_UA_BL_Type ! ======================= ! ========= AFI_UA_BL_Default_Type ======= @@ -131,6 +137,9 @@ MODULE AirfoilInfo_Types LOGICAL :: filtCutOff = .true. !< Calculate value for this input? [-] LOGICAL :: alphaUpper = .true. !< Calculate value for this input? [-] LOGICAL :: alphaLower = .true. !< Calculate value for this input? [-] + LOGICAL :: Ka = .true. !< Calculate value for this input? [-] + LOGICAL :: Kv = .true. !< Calculate value for this input? [-] + LOGICAL :: dCNdA = .true. !< Calculate value for this input? [-] END TYPE AFI_UA_BL_Default_Type ! ======================= ! ========= AFI_Table_Type ======= @@ -272,6 +281,11 @@ subroutine AFI_CopyUA_BL_Type(SrcUA_BL_TypeData, DstUA_BL_TypeData, CtrlCode, Er DstUA_BL_TypeData%CnBreakUpper = SrcUA_BL_TypeData%CnBreakUpper DstUA_BL_TypeData%alphaBreakLower = SrcUA_BL_TypeData%alphaBreakLower DstUA_BL_TypeData%CnBreakLower = SrcUA_BL_TypeData%CnBreakLower + DstUA_BL_TypeData%Ka = SrcUA_BL_TypeData%Ka + DstUA_BL_TypeData%Kv = SrcUA_BL_TypeData%Kv + DstUA_BL_TypeData%dCNdA = SrcUA_BL_TypeData%dCNdA + DstUA_BL_TypeData%CnMax = SrcUA_BL_TypeData%CnMax + DstUA_BL_TypeData%CnMin = SrcUA_BL_TypeData%CnMin end subroutine subroutine AFI_DestroyUA_BL_Type(UA_BL_TypeData, ErrStat, ErrMsg) @@ -332,6 +346,11 @@ subroutine AFI_PackUA_BL_Type(RF, Indata) call RegPack(RF, InData%CnBreakUpper) call RegPack(RF, InData%alphaBreakLower) call RegPack(RF, InData%CnBreakLower) + call RegPack(RF, InData%Ka) + call RegPack(RF, InData%Kv) + call RegPack(RF, InData%dCNdA) + call RegPack(RF, InData%CnMax) + call RegPack(RF, InData%CnMin) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -384,6 +403,11 @@ subroutine AFI_UnPackUA_BL_Type(RF, OutData) call RegUnpack(RF, OutData%CnBreakUpper); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%alphaBreakLower); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%CnBreakLower); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%Ka); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%Kv); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%dCNdA); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%CnMax); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%CnMin); if (RegCheckErr(RF, RoutineName)) return end subroutine subroutine AFI_CopyUA_BL_Default_Type(SrcUA_BL_Default_TypeData, DstUA_BL_Default_TypeData, CtrlCode, ErrStat, ErrMsg) @@ -431,6 +455,9 @@ subroutine AFI_CopyUA_BL_Default_Type(SrcUA_BL_Default_TypeData, DstUA_BL_Defaul DstUA_BL_Default_TypeData%filtCutOff = SrcUA_BL_Default_TypeData%filtCutOff DstUA_BL_Default_TypeData%alphaUpper = SrcUA_BL_Default_TypeData%alphaUpper DstUA_BL_Default_TypeData%alphaLower = SrcUA_BL_Default_TypeData%alphaLower + DstUA_BL_Default_TypeData%Ka = SrcUA_BL_Default_TypeData%Ka + DstUA_BL_Default_TypeData%Kv = SrcUA_BL_Default_TypeData%Kv + DstUA_BL_Default_TypeData%dCNdA = SrcUA_BL_Default_TypeData%dCNdA end subroutine subroutine AFI_DestroyUA_BL_Default_Type(UA_BL_Default_TypeData, ErrStat, ErrMsg) @@ -483,6 +510,9 @@ subroutine AFI_PackUA_BL_Default_Type(RF, Indata) call RegPack(RF, InData%filtCutOff) call RegPack(RF, InData%alphaUpper) call RegPack(RF, InData%alphaLower) + call RegPack(RF, InData%Ka) + call RegPack(RF, InData%Kv) + call RegPack(RF, InData%dCNdA) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -527,6 +557,9 @@ subroutine AFI_UnPackUA_BL_Default_Type(RF, OutData) call RegUnpack(RF, OutData%filtCutOff); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%alphaUpper); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%alphaLower); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%Ka); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%Kv); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%dCNdA); if (RegCheckErr(RF, RoutineName)) return end subroutine subroutine AFI_CopyTable_Type(SrcTable_TypeData, DstTable_TypeData, CtrlCode, ErrStat, ErrMsg) @@ -1351,6 +1384,11 @@ SUBROUTINE AFI_UA_BL_Type_ExtrapInterp1(u1, u2, tin, u_out, tin_out, ErrStat, Er u_out%CnBreakUpper = a1*u1%CnBreakUpper + a2*u2%CnBreakUpper CALL Angles_ExtrapInterp( u1%alphaBreakLower, u2%alphaBreakLower, tin, u_out%alphaBreakLower, tin_out ) u_out%CnBreakLower = a1*u1%CnBreakLower + a2*u2%CnBreakLower + u_out%Ka = a1*u1%Ka + a2*u2%Ka + u_out%Kv = a1*u1%Kv + a2*u2%Kv + u_out%dCNdA = a1*u1%dCNdA + a2*u2%dCNdA + u_out%CnMax = a1*u1%CnMax + a2*u2%CnMax + u_out%CnMin = a1*u1%CnMin + a2*u2%CnMin END SUBROUTINE SUBROUTINE AFI_UA_BL_Type_ExtrapInterp2(u1, u2, u3, tin, u_out, tin_out, ErrStat, ErrMsg ) @@ -1450,6 +1488,11 @@ SUBROUTINE AFI_UA_BL_Type_ExtrapInterp2(u1, u2, u3, tin, u_out, tin_out, ErrStat u_out%CnBreakUpper = a1*u1%CnBreakUpper + a2*u2%CnBreakUpper + a3*u3%CnBreakUpper CALL Angles_ExtrapInterp( u1%alphaBreakLower, u2%alphaBreakLower, u3%alphaBreakLower, tin, u_out%alphaBreakLower, tin_out ) u_out%CnBreakLower = a1*u1%CnBreakLower + a2*u2%CnBreakLower + a3*u3%CnBreakLower + u_out%Ka = a1*u1%Ka + a2*u2%Ka + a3*u3%Ka + u_out%Kv = a1*u1%Kv + a2*u2%Kv + a3*u3%Kv + u_out%dCNdA = a1*u1%dCNdA + a2*u2%dCNdA + a3*u3%dCNdA + u_out%CnMax = a1*u1%CnMax + a2*u2%CnMax + a3*u3%CnMax + u_out%CnMin = a1*u1%CnMin + a2*u2%CnMin + a3*u3%CnMin END SUBROUTINE function AFI_InputMeshPointer(u, DL) result(Mesh) diff --git a/modules/aerodyn/src/UnsteadyAero.f90 b/modules/aerodyn/src/UnsteadyAero.f90 index 396cc74197..5d3cc1285f 100644 --- a/modules/aerodyn/src/UnsteadyAero.f90 +++ b/modules/aerodyn/src/UnsteadyAero.f90 @@ -2601,8 +2601,9 @@ SUBROUTINE HGM_Steady( i, j, u, p, x, AFInfo, ErrStat, ErrMsg ) ! States !x1: Downwash memory term 1 (rad) !x2: Downwash memory term 2 (rad) - !x3: Clp', Lift coefficient with a time lag to the attached lift coeff + !x3: lagged attached-flow coefficient: Clp' (lift) for UA_Mod=4; Cnp (normal force) for UA_Mod=5,8,9 !x4: f'' , Final separation point function + !x5: vortex-lift normal force (UA_Mod=5,9) ! Steady states @@ -2737,8 +2738,9 @@ subroutine UA_CalcContStateDeriv( i, j, t, u_in, p, x, OtherState, AFInfo, m, dx ! States !x1: Downwash memory term 1 (rad) !x2: Downwash memory term 2 (rad) - !x3: Clp', Lift coefficient with a time lag to the attached lift coeff + !x3: lagged attached-flow coefficient: Clp' (lift) for UA_Mod=4; Cnp (normal force) for UA_Mod=5,8,9 !x4: f'' , Final separation point function + !x5: vortex-lift normal force (UA_Mod=5,9) ! Constraining x4 between 0 and 1 increases numerical stability (should be done elsewhere, but we'll double check here in case there were perturbations on the state value) x4 = max( min( x%x(4), 1.0_R8Ki ), 0.0_R8Ki ) diff --git a/modules/aerodyn/src/UnsteadyAero_Registry.txt b/modules/aerodyn/src/UnsteadyAero_Registry.txt index 3b3a30df66..d5a6be2370 100644 --- a/modules/aerodyn/src/UnsteadyAero_Registry.txt +++ b/modules/aerodyn/src/UnsteadyAero_Registry.txt @@ -108,7 +108,7 @@ typedef ^ UA_KelvinChainType ReKi # ..... States .................................................................................................................... # Define continuous (differentiable) states here: -typedef ^ UA_ElementContinuousStateType R8Ki x 7 - - "continuous states when UA_Mod=4 (x1 and x2:Downwash memory terms; x3:Clp', Lift coefficient with a time lag to the attached lift coeff; x4: f'' , Final separation point function)" "{rad, rad, - -}" +typedef ^ UA_ElementContinuousStateType R8Ki x 7 - - "continuous states for UA_Mod=4,5,6,8,9 (x1,x2: downwash memory terms; x3: lagged attached-flow coefficient -- Clp' for UA_Mod=4, Cnp for UA_Mod=5,8,9; x4: f'', final separation point function; x5: vortex-lift normal force, UA_Mod=5,9)" "{rad, rad, - -}" typedef ^ ContinuousStateType UA_ElementContinuousStateType element {:}{:} - - "continuous states when UA_Mod=4 for each blade/node" "-" # Define discrete (non-differentiable) states here: @@ -168,6 +168,10 @@ typedef ^ OtherStateType ReKi typedef ^ OtherStateType LOGICAL PositivePressure {:}{:} - - "HGMV model: logical flag indicating if the vortex lift became active because of positive pressure (or negative)" - typedef ^ OtherStateType LOGICAL vortexOn {:}{:} - - "HGMV model: logical flag indicating if the vortex lift term is active" - typedef ^ OtherStateType LOGICAL BelowThreshold {:}{:} - - "HGMV model: logical flag indicating if cn fell below threshold to form another vortex" - +typedef ^ OtherStateType ReKi tau_v_IAG {:}{:} - - "IAG model: non-dimensional vortex time" - +typedef ^ OtherStateType LOGICAL VortexOn_IAG {:}{:} - - "IAG model: logical flag indicating if the vortex lift term is active" - +typedef ^ OtherStateType LOGICAL PositiveStall {:}{:} - - "IAG model: sign of the stall state (x3>=0) latched at vortex initiation" - +typedef ^ OtherStateType ReKi alpha_minus1_IAG {:}{:} - - "IAG model: angle of attack at the previous step, used for the upstroke test" rad typedef ^ OtherStateType LOGICAL activeL {:}{:} - - "BV model: logical flag indicating if the lift stall is active" - typedef ^ OtherStateType LOGICAL activeD {:}{:} - - "BV model: logical flag indicating if the drag stall is active" - diff --git a/modules/aerodyn/src/UnsteadyAero_Types.f90 b/modules/aerodyn/src/UnsteadyAero_Types.f90 index 41a52d6d4c..9bf4591904 100644 --- a/modules/aerodyn/src/UnsteadyAero_Types.f90 +++ b/modules/aerodyn/src/UnsteadyAero_Types.f90 @@ -121,7 +121,7 @@ MODULE UnsteadyAero_Types ! ======================= ! ========= UA_ElementContinuousStateType ======= TYPE, PUBLIC :: UA_ElementContinuousStateType - REAL(R8Ki) , DIMENSION(1:7) :: x = 0.0_R8Ki !< continuous states when UA_Mod=4 (x1 and x2:Downwash memory terms; x3:Clp', Lift coefficient with a time lag to the attached lift coeff; x4: f'' , Final separation point function) [{rad, rad, - -}] + REAL(R8Ki) , DIMENSION(1:7) :: x = 0.0_R8Ki !< continuous states for UA_Mod=4,5,6,8,9 (x1,x2: downwash memory terms; x3: lagged attached-flow coefficient -- Clp' for UA_Mod=4, Cnp for UA_Mod=5,8,9; x4: f'', final separation point function; x5: vortex-lift normal force, UA_Mod=5,9) [{rad, rad, - -}] END TYPE UA_ElementContinuousStateType ! ======================= ! ========= UA_ContinuousStateType ======= @@ -187,6 +187,10 @@ MODULE UnsteadyAero_Types LOGICAL , DIMENSION(:,:), ALLOCATABLE :: PositivePressure !< HGMV model: logical flag indicating if the vortex lift became active because of positive pressure (or negative) [-] LOGICAL , DIMENSION(:,:), ALLOCATABLE :: vortexOn !< HGMV model: logical flag indicating if the vortex lift term is active [-] LOGICAL , DIMENSION(:,:), ALLOCATABLE :: BelowThreshold !< HGMV model: logical flag indicating if cn fell below threshold to form another vortex [-] + REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: tau_v_IAG !< IAG model: non-dimensional vortex time [-] + LOGICAL , DIMENSION(:,:), ALLOCATABLE :: VortexOn_IAG !< IAG model: logical flag indicating if the vortex lift term is active [-] + LOGICAL , DIMENSION(:,:), ALLOCATABLE :: PositiveStall !< IAG model: sign of the stall state (x3>=0) latched at vortex initiation [-] + REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: alpha_minus1_IAG !< IAG model: angle of attack at the previous step, used for the upstroke test [rad] LOGICAL , DIMENSION(:,:), ALLOCATABLE :: activeL !< BV model: logical flag indicating if the lift stall is active [-] LOGICAL , DIMENSION(:,:), ALLOCATABLE :: activeD !< BV model: logical flag indicating if the drag stall is active [-] END TYPE UA_OtherStateType @@ -1621,6 +1625,54 @@ subroutine UA_CopyOtherState(SrcOtherStateData, DstOtherStateData, CtrlCode, Err end if DstOtherStateData%BelowThreshold = SrcOtherStateData%BelowThreshold end if + if (allocated(SrcOtherStateData%tau_v_IAG)) then + LB(1:2) = lbound(SrcOtherStateData%tau_v_IAG) + UB(1:2) = ubound(SrcOtherStateData%tau_v_IAG) + if (.not. allocated(DstOtherStateData%tau_v_IAG)) then + allocate(DstOtherStateData%tau_v_IAG(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstOtherStateData%tau_v_IAG.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstOtherStateData%tau_v_IAG = SrcOtherStateData%tau_v_IAG + end if + if (allocated(SrcOtherStateData%VortexOn_IAG)) then + LB(1:2) = lbound(SrcOtherStateData%VortexOn_IAG) + UB(1:2) = ubound(SrcOtherStateData%VortexOn_IAG) + if (.not. allocated(DstOtherStateData%VortexOn_IAG)) then + allocate(DstOtherStateData%VortexOn_IAG(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstOtherStateData%VortexOn_IAG.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstOtherStateData%VortexOn_IAG = SrcOtherStateData%VortexOn_IAG + end if + if (allocated(SrcOtherStateData%PositiveStall)) then + LB(1:2) = lbound(SrcOtherStateData%PositiveStall) + UB(1:2) = ubound(SrcOtherStateData%PositiveStall) + if (.not. allocated(DstOtherStateData%PositiveStall)) then + allocate(DstOtherStateData%PositiveStall(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstOtherStateData%PositiveStall.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstOtherStateData%PositiveStall = SrcOtherStateData%PositiveStall + end if + if (allocated(SrcOtherStateData%alpha_minus1_IAG)) then + LB(1:2) = lbound(SrcOtherStateData%alpha_minus1_IAG) + UB(1:2) = ubound(SrcOtherStateData%alpha_minus1_IAG) + if (.not. allocated(DstOtherStateData%alpha_minus1_IAG)) then + allocate(DstOtherStateData%alpha_minus1_IAG(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstOtherStateData%alpha_minus1_IAG.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstOtherStateData%alpha_minus1_IAG = SrcOtherStateData%alpha_minus1_IAG + end if if (allocated(SrcOtherStateData%activeL)) then LB(1:2) = lbound(SrcOtherStateData%activeL) UB(1:2) = ubound(SrcOtherStateData%activeL) @@ -1703,6 +1755,18 @@ subroutine UA_DestroyOtherState(OtherStateData, ErrStat, ErrMsg) if (allocated(OtherStateData%BelowThreshold)) then deallocate(OtherStateData%BelowThreshold) end if + if (allocated(OtherStateData%tau_v_IAG)) then + deallocate(OtherStateData%tau_v_IAG) + end if + if (allocated(OtherStateData%VortexOn_IAG)) then + deallocate(OtherStateData%VortexOn_IAG) + end if + if (allocated(OtherStateData%PositiveStall)) then + deallocate(OtherStateData%PositiveStall) + end if + if (allocated(OtherStateData%alpha_minus1_IAG)) then + deallocate(OtherStateData%alpha_minus1_IAG) + end if if (allocated(OtherStateData%activeL)) then deallocate(OtherStateData%activeL) end if @@ -1739,6 +1803,10 @@ subroutine UA_PackOtherState(RF, Indata) call RegPackAlloc(RF, InData%PositivePressure) call RegPackAlloc(RF, InData%vortexOn) call RegPackAlloc(RF, InData%BelowThreshold) + call RegPackAlloc(RF, InData%tau_v_IAG) + call RegPackAlloc(RF, InData%VortexOn_IAG) + call RegPackAlloc(RF, InData%PositiveStall) + call RegPackAlloc(RF, InData%alpha_minus1_IAG) call RegPackAlloc(RF, InData%activeL) call RegPackAlloc(RF, InData%activeD) if (RegCheckErr(RF, RoutineName)) return @@ -1774,6 +1842,10 @@ subroutine UA_UnPackOtherState(RF, OutData) call RegUnpackAlloc(RF, OutData%PositivePressure); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%vortexOn); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%BelowThreshold); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%tau_v_IAG); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%VortexOn_IAG); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%PositiveStall); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%alpha_minus1_IAG); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%activeL); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%activeD); if (RegCheckErr(RF, RoutineName)) return end subroutine From 9ec93d17cfd680e2467c132b7d9a8c3e31a20f7e Mon Sep 17 00:00:00 2001 From: Hannah Ross Date: Mon, 28 Sep 2026 14:32:12 -0600 Subject: [PATCH 02/12] Add IAG dynamic stall airfoil preprocessing (dCNdA fit, separation functions, CnMax/CnMin) to AirfoilInfo with summary-table output --- modules/aerodyn/src/AirfoilInfo.f90 | 424 +++++++++++++++++++++++++++- 1 file changed, 415 insertions(+), 9 deletions(-) diff --git a/modules/aerodyn/src/AirfoilInfo.f90 b/modules/aerodyn/src/AirfoilInfo.f90 index bd4409c627..9177ad97ac 100644 --- a/modules/aerodyn/src/AirfoilInfo.f90 +++ b/modules/aerodyn/src/AirfoilInfo.f90 @@ -43,6 +43,23 @@ MODULE AirfoilInfo integer, parameter :: MaxNumAFCoeffs = 7 !cl,cd,cm,cpMin, UA:f_st, FullySeparate, FullyAttached + !> Sentinel written into UA_BL%CnMax / CnMin when the IAG (UAMod=9) vortex logic must stay + !! dormant: non-IAG models, and degenerate (cylinder-like) polars with no stall peak. + !! + !! DO NOT change this to huge() or a similarly extreme value. Every UA_BL field is blended + !! across airfoil tables by the registry-generated AFI_UA_BL_Type_ExtrapInterp routines + !! (AirfoilInfo_Types.f90) as u_out%CnMax = a1*u1%CnMax + a2*u2%CnMax. With huge() that + !! arithmetic misbehaves in three silent ways: + !! - one table degenerate + one normal -> ~0.5*huge, not a usable CnMax; + !! - extrapolation weights (a1>1, a2<0) -> overflow to +/-Infinity; + !! - the value also overflows the F17.5 summary-table format into '*****'. + !! None of these trap, because OpenFAST does not build with -ffpe-trap. + !! + !! 999 is ~500x any physically meaningful |Cn|, so the abs(x3) > CN_CRIT vortex test can + !! still never be satisfied, while the value stays finite, interpolates harmlessly and + !! prints cleanly. + real(ReKi), parameter :: IAG_CnCrit_Disabled = 999.0_ReKi + CONTAINS @@ -657,6 +674,16 @@ SUBROUTINE ReadAFfile ( InitInp, NumCoefsIn, p, ErrStat, ErrMsg, UnEc ) CALL ParseVarWDefault ( FileInfo, CurLine, 'x_cp_bar', p%Table(iTable)%UA_BL%x_cp_bar, .2_ReKi, ErrStat2, ErrMsg2, UnEc ) CalcDefaults(iTable)%x_cp_bar = ErrStat2 >= AbortErrLev + + ! IAG model (UAMod=9) parameters: + CALL ParseVarWDefault ( FileInfo, CurLine, 'Ka', p%Table(iTable)%UA_BL%Ka, 0.75_ReKi, ErrStat2, ErrMsg2, UnEc ) + CalcDefaults(iTable)%Ka = ErrStat2 >= AbortErrLev + + CALL ParseVarWDefault ( FileInfo, CurLine, 'Kv', p%Table(iTable)%UA_BL%Kv, 0.2_ReKi, ErrStat2, ErrMsg2, UnEc ) + CalcDefaults(iTable)%Kv = ErrStat2 >= AbortErrLev + + CALL ParseVarWDefault ( FileInfo, CurLine, 'dCNdA', p%Table(iTable)%UA_BL%dCNdA, TwoPi, ErrStat2, ErrMsg2, UnEc ) + CalcDefaults(iTable)%dCNdA = ErrStat2 >= AbortErrLev CALL ParseVarWDefault ( FileInfo, CurLine, 'UACutout', p%Table(iTable)%UA_BL%UACutout, 45.0_ReKi, ErrStat2, ErrMsg2, UnEc ) CalcDefaults(iTable)%UACutout = ErrStat2 >= AbortErrLev @@ -757,7 +784,12 @@ SUBROUTINE ReadAFfile ( InitInp, NumCoefsIn, p, ErrStat, ErrMsg, UnEc ) if ( p%Table(iTable)%ConstData ) then p%Table(iTable)%InclUAdata = .false. else - call CalculateUACoeffs(CalcDefaults(iTable), p%Table(iTable), p%ColCl, p%ColCd, p%ColCm, p%ColUAf, InitInp%UAMod) + call CalculateUACoeffs(CalcDefaults(iTable), p%Table(iTable), p%ColCl, p%ColCd, p%ColCm, p%ColUAf, InitInp%UAMod, ErrStat2, ErrMsg2) + CALL SetErrStat( ErrStat2, trim(ErrMsg2)//' (airfoil table '//trim(num2lstr(iTable))//' of "'//TRIM( InitInp%FileName )//'")', ErrStat, ErrMsg, RoutineName ) + IF (ErrStat >= AbortErrLev) THEN + CALL Cleanup() + RETURN + END IF end if ! Let's make sure that the data go from -Pi to Pi and that the values are the same for both @@ -824,7 +856,7 @@ END SUBROUTINE Cleanup END SUBROUTINE ReadAFfile !---------------------------------------------------------------------------------------------------------------------------------- - SUBROUTINE CalculateUACoeffs(CalcDefaults,p,ColCl,ColCd,ColCm,ColUAf,UAMod) + SUBROUTINE CalculateUACoeffs(CalcDefaults,p,ColCl,ColCd,ColCm,ColUAf,UAMod,ErrStat,ErrMsg) TYPE (AFI_UA_BL_Default_Type),intent(in):: CalcDefaults TYPE (AFI_Table_Type), intent(inout) :: p ! This structure stores all the module parameters that are set by AirfoilInfo during the initialization phase. integer(IntKi), intent(in ) :: ColCl ! column for cl @@ -832,6 +864,8 @@ SUBROUTINE CalculateUACoeffs(CalcDefaults,p,ColCl,ColCd,ColCm,ColUAf,UAMod) integer(IntKi), intent(in ) :: ColCm ! column for cm integer(IntKi), intent(in ) :: ColUAf ! column for UA f_st (based on Cl or cn) integer(IntKi), intent(in ) :: UAMod ! UA model; determines how to compute f_st? + integer(IntKi), intent( out) :: ErrStat ! Error status of the operation + character(*), intent( out) :: ErrMsg ! Error message if ErrStat /= ErrID_None INTEGER(IntKi) :: Row ! The row of a table to be parsed in the FileInfo structure. INTEGER(IntKi) :: col_fs ! column for UA cn/cl_fs (fully separated cn or cl) @@ -872,6 +906,9 @@ SUBROUTINE CalculateUACoeffs(CalcDefaults,p,ColCl,ColCd,ColCm,ColUAf,UAMod) CHARACTER(ErrMsgLen) :: ErrMsg2 CHARACTER(*), PARAMETER :: RoutineName = 'CalculateUACoeffs' + ErrStat = ErrID_None + ErrMsg = "" + if ( UAMod == UA_HGMV360 ) then LimitAlphaRange = TwoPi ! range we're limiting our equations to (in radians) else @@ -911,6 +948,31 @@ SUBROUTINE CalculateUACoeffs(CalcDefaults,p,ColCl,ColCd,ColCm,ColUAf,UAMod) if (CalcDefaults%x_cp_bar ) p%UA_BL%x_cp_bar = 0.20_ReKi if (CalcDefaults%filtCutOff ) p%UA_BL%filtCutOff = 0.50_ReKi + if (UAMod == UA_IAG) then + ! IAG model defaults from Bangga, Parkinson & Collier (2023), Table 1. + ! NOTE: b1 and T_VL differ from OpenFAST's B-L defaults set above, so they are + ! overridden here (only when the user did not supply a value). + if (CalcDefaults%Ka ) p%UA_BL%Ka = 0.75_ReKi + if (CalcDefaults%Kv ) p%UA_BL%Kv = 0.20_ReKi + if (CalcDefaults%A1 ) p%UA_BL%A1 = 0.30_ReKi + if (CalcDefaults%A2 ) p%UA_BL%A2 = 0.70_ReKi + if (CalcDefaults%b1 ) p%UA_BL%b1 = 0.70_ReKi ! IAG value (OpenFAST B-L default is 0.14) + if (CalcDefaults%b2 ) p%UA_BL%b2 = 0.53_ReKi + if (CalcDefaults%T_p ) p%UA_BL%T_p = 1.70_ReKi + if (CalcDefaults%T_f0 ) p%UA_BL%T_f0 = 3.00_ReKi + if (CalcDefaults%T_V0 ) p%UA_BL%T_V0 = 6.00_ReKi + if (CalcDefaults%T_VL ) p%UA_BL%T_VL = 6.00_ReKi ! IAG value (OpenFAST B-L default is 11.0) + else + p%UA_BL%Ka = 0.0_ReKi + p%UA_BL%Kv = 0.0_ReKi + p%UA_BL%dCNdA = 0.0_ReKi + end if + + ! these are only meaningful for UAMod=UA_IAG; initialize so that the vortex logic + ! can never trigger for other models (or for degenerate polars, see below) + p%UA_BL%CnMax = IAG_CnCrit_Disabled + p%UA_BL%CnMin = -IAG_CnCrit_Disabled + if (UAMod == UA_HGMV360) then ! set defaults for this model (note: we don't turn off UA) if (CalcDefaults%St_sh ) p%UA_BL%St_sh = 0.14_ReKi if (CalcDefaults%UACutout ) p%UA_BL%UACutout = TwoPi*2 ! don't turn off UA for this model @@ -955,8 +1017,19 @@ SUBROUTINE CalculateUACoeffs(CalcDefaults,p,ColCl,ColCd,ColCm,ColUAf,UAMod) end if Cn = Calculate_Cn(alpha=p%alpha, cl=p%Coefs(:,ColCl), cd=p%Coefs(:,ColCd), cd0=p%UA_BL%Cd0) call ComputeUA360_AttachedFlow(p, ColUAf, Cn, iLower, iUpper) - call ComputeUA360_updateSeparationF( p, ColUAf, Cn, iLower, iUpper ) - call ComputeUA360_updateCnSeparated( p, ColUAf, Cn, iLower ) + + if ( UAMod == UA_IAG ) then + ! Degenerate (cylinder-like) polar: there is no meaningful lift slope or stall + ! peak. Use a nominal slope so the Eq. 42 denominator is not singular, and leave + ! CnMax/CnMin at their +/-huge sentinels so the vortex logic stays dormant. + ! The Kirchhoff consistency check is meaningless on such a polar, so it is skipped. + if (CalcDefaults%dCNdA) p%UA_BL%dCNdA = TwoPi + call ComputeIAG_SeparationFunctions( p, ColUAf, Cn, iLower, iUpper, CheckConsistency=.false., ErrStat=ErrStat2, ErrMsg=ErrMsg2 ) + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + else + call ComputeUA360_updateSeparationF( p, ColUAf, Cn, iLower, iUpper ) + call ComputeUA360_updateCnSeparated( p, ColUAf, Cn, iLower ) + end if end if else @@ -1085,8 +1158,33 @@ SUBROUTINE CalculateUACoeffs(CalcDefaults,p,ColCl,ColCd,ColCm,ColUAf,UAMod) call Compute_iLoweriUpper(p, iLower, iUpper) ! calculating iLower and iUpper here (for alpha1 and alpha2) else call ComputeUA360_AttachedFlow(p, ColUAf, Cn, iLower, iUpper) - call ComputeUA360_updateSeparationF( p, ColUAf, Cn, iLower, iUpper ) - call ComputeUA360_updateCnSeparated( p, ColUAf, Cn, iLower ) + + if ( UAMod == UA_IAG ) then + !------------------------------------------------------------------------ + ! IAG model. The separation-function columns are built from IAG's own + ! sinusoidal attached-flow curve dCNdA*sin(alpha-alpha0) (paper Eq. 37), + ! so dCNdA must be resolved FIRST -- including any user override, or the + ! columns and the run-time CN_f (Eq. 44) would use different slopes. + !------------------------------------------------------------------------ + if (CalcDefaults%dCNdA) then + call ComputeIAG_dCNdA( p%alpha, cn, p%UA_BL%alpha0, p%UA_BL%C_nalpha, p%UA_BL%dCNdA ) + end if + + if (p%UA_BL%dCNdA <= 0.0_ReKi) then + ! A zero or negative slope makes the Eq. 42 denominator degenerate for every + ! row of the table. UA_ValidateAFI issues the user-facing error; fall back to + ! the thin-airfoil value here so the columns below remain well defined. + p%UA_BL%dCNdA = TwoPi + end if + + call ComputeIAG_SeparationFunctions( p, ColUAf, Cn, iLower, iUpper, CheckConsistency=.true., ErrStat=ErrStat2, ErrMsg=ErrMsg2 ) + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + if (ErrStat >= AbortErrLev) return + call ComputeIAG_CnMaxCnMin( p, cn, iLower, iUpper ) + else + call ComputeUA360_updateSeparationF( p, ColUAf, Cn, iLower, iUpper ) + call ComputeUA360_updateCnSeparated( p, ColUAf, Cn, iLower ) + end if end if @@ -1578,6 +1676,302 @@ REAL(ReKi) FUNCTION ComputeUA360_CnOffset(p, cn_cl, Row, iLower) RESULT(offset) SlopeScale = 0.1_ReKi*R2D offset = CnOffset * ( tanh(SlopeScale*(p%alpha(Row)+PiBy2)) - tanh(SlopeScale*(p%alpha(Row)-PiBy2)) ) / 2.0_ReKi; !Only apply Cn offset in vicinity of AoA 0 deg END FUNCTION ComputeUA360_CnOffset +!---------------------------------------------------------------------------------------------------------------------------------- + SUBROUTINE ComputeIAG_dCNdA( alpha, cn, alpha0, C_nalpha, dCNdA ) + ! IAG "linear fit approach" for the static normal-force curve slope. + ! Bangga, Parkinson & Collier (2023), Energies 16:3994, section 2.4, p. 11: + ! (i) compute the local gradient of cn(alpha) at every polar point between alpha0 and 7 deg; + ! (ii) retain only points whose local gradient lies between 1.8*pi and 2.5*pi; + ! (iii) least-squares fit cn = dCNdA*(alpha - alpha0) through the retained points. + ! If too few points survive the filter, fall back to the B-L C_nalpha. + REAL(ReKi), intent(in ) :: alpha(:) ! angle of attack table (rad) + REAL(ReKi), intent(in ) :: cn(:) ! static normal force coefficient + REAL(ReKi), intent(in ) :: alpha0 ! zero-lift angle of attack (rad) + REAL(ReKi), intent(in ) :: C_nalpha ! fallback slope + REAL(ReKi), intent( out) :: dCNdA ! resulting slope (1/rad) + + REAL(ReKi) :: SlopeMin ! lower gradient filter: 1.8*pi + REAL(ReKi) :: SlopeMax ! upper gradient filter: 2.5*pi + REAL(ReKi) :: alphaFitMax ! upper edge of the fit window: 7 deg + + REAL(ReKi) :: LocalSlope + REAL(ReKi) :: da, dcn + REAL(ReKi) :: SumXX, SumXY + INTEGER(IntKi) :: Row, NumAlf, nKept + + SlopeMin = 1.8_ReKi * Pi + SlopeMax = 2.5_ReKi * Pi + alphaFitMax = 7.0_ReKi * D2R + + NumAlf = size(alpha) + SumXX = 0.0_ReKi + SumXY = 0.0_ReKi + nKept = 0 + + do Row = 2, NumAlf-1 + if (alpha(Row) < alpha0 .or. alpha(Row) > alphaFitMax) cycle + + ! central-difference local gradient + da = alpha(Row+1) - alpha(Row-1) + if (EqualRealNos(da, 0.0_ReKi)) cycle + LocalSlope = ( cn(Row+1) - cn(Row-1) ) / da + if (LocalSlope < SlopeMin .or. LocalSlope > SlopeMax) cycle + + ! accumulate the (zero-intercept) least-squares fit about alpha0 + da = alpha(Row) - alpha0 + dcn = cn(Row) + SumXX = SumXX + da*da + SumXY = SumXY + da*dcn + nKept = nKept + 1 + end do + + if (nKept >= 2 .and. SumXX > 0.0_ReKi) then + dCNdA = SumXY / SumXX + else + dCNdA = C_nalpha ! fallback: too few points passed the gradient filter + end if + + END SUBROUTINE ComputeIAG_dCNdA +!---------------------------------------------------------------------------------------------------------------------------------- + SUBROUTINE ComputeIAG_SeparationFunctions( p, ColUAf, cn, iLower, iUpper, CheckConsistency, ErrStat, ErrMsg ) + ! Build the IAG-specific FullyAttached / f_st / FullySeparate columns. + ! + ! Unlike the HGM/HGMV/HGMV360 path, IAG's attached-flow curve is the sinusoidal + ! CN_C = dCNdA*sin(alphaE - alpha0) (paper Eq. 37), and its separated normal force is the + ! squared-Kirchhoff blend CN_f = dCNdA*((1+sqrt(x4))/2)^2 * sin(alphaE-alpha0) + CN_I (Eq. 44). + ! The separation function must therefore be inverted from that same curve (Eq. 42), or f and + ! CN_f are inconsistent and the model fails to reduce to the static polar at zero pitch rate. + ! + ! Note the FullySeparate column is NOT used by IAG at run time (CN_f reads only dCNdA, alpha0 + ! and x4). It is filled anyway so AFI_WrTables writes something meaningful and the array is + ! never left uninitialized. + TYPE (AFI_Table_Type), intent(inout) :: p ! airfoil table + integer(IntKi), intent(in ) :: ColUAf ! column for UA f_st + REAL(ReKi), intent(in ) :: cn(:) ! static normal force coefficient + INTEGER(IntKi), intent(in ) :: iLower ! lower index of the fully attached region + INTEGER(IntKi), intent(in ) :: iUpper ! upper index of the fully attached region + LOGICAL, intent(in ) :: CheckConsistency ! verify the Kirchhoff identity (skip for degenerate polars) + integer(IntKi), intent( out) :: ErrStat ! Error status of the operation + character(*), intent( out) :: ErrMsg ! Error message if ErrStat /= ErrID_None + + REAL(ReKi) :: cn_fa ! IAG fully-attached cn at this alpha + REAL(ReKi) :: CnRatio + REAL(ReKi) :: CnFaTol ! tolerance for the cn_fa zero crossings + REAL(ReKi) :: f_st_(p%NumAlf) ! temporary for the smoothing pass + LOGICAL :: RowIsExact(p%NumAlf) ! .true. where the Kirchhoff identity must hold exactly + REAL(ReKi) :: cn_check ! cn reconstructed from the columns + REAL(ReKi) :: Residual, MaxResidual + REAL(ReKi) :: ResidualTol + INTEGER(IntKi) :: Row, iWorst, nChecked + INTEGER(IntKi) :: col_fs, col_fa + + ErrStat = ErrID_None + ErrMsg = "" + + col_fs = ColUAf + 1 ! fully separated + col_fa = ColUAf + 2 ! fully attached + + ! cn produced by roughly 0.6 deg of incidence; small enough not to perturb the polar, + ! large enough that we never divide by a near-zero cn_fa. + CnFaTol = 0.01_ReKi * p%UA_BL%dCNdA + + RowIsExact = .true. + + do Row = 1, p%NumAlf + cn_fa = p%UA_BL%dCNdA * sin( p%alpha(Row) - p%UA_BL%alpha0 ) ! Eq. 37 curve + p%Coefs(Row,col_fa) = cn_fa + + ! cn_fa -> 0 at alpha0 and again at alpha0 +/- 180 deg. These two cases are NOT the + ! same: near alpha0 the flow is attached (cn -> 0 too, so the ratio is finite), but + ! near +/-180 deg the section is fully stalled. Do not collapse them. + if ( abs(cn_fa) < CnFaTol ) then + if ( abs(p%alpha(Row) - p%UA_BL%alpha0) < PiBy2 ) then + CnRatio = 1.0_ReKi ! attached branch: f = 1 + else + CnRatio = 0.0_ReKi ! reverse-flow branch: clipped to f = 0 below + end if + RowIsExact(Row) = .false. ! CnRatio substituted, not computed from the polar + else + CnRatio = cn(Row) / cn_fa + end if + + if ( CnRatio < 0.25_ReKi ) then + CnRatio = 0.25_ReKi ! below 1/4 => fully separated + RowIsExact(Row) = .false. ! clipped + end if + + p%Coefs(Row,ColUAf) = ( 2.0_ReKi * sqrt( CnRatio ) - 1.0_ReKi )**2 + + if ( p%Coefs(Row,ColUAf) > 1.0_ReKi ) then + p%Coefs(Row,ColUAf) = 1.0_ReKi ! f <= 1 + RowIsExact(Row) = .false. ! clipped + end if + p%Coefs(Row,ColUAf) = max( 0.0_ReKi, p%Coefs(Row,ColUAf) ) ! f >= 0 + end do + + ! Where the IAG attached curve matches the polar by construction, set f = 1. + ! This is imposed, not derived, so these rows cannot be held to the identity. + do Row = iLower, iUpper + p%Coefs(Row,ColUAf) = 1.0_ReKi + RowIsExact(Row) = .false. + end do + + ! light smoothing, mirroring the conditioning applied on the HGMV360 path. + ! Any row the smoother actually moves is likewise no longer exact. + f_st_ = p%Coefs(:,ColUAf) + do Row = iUpper+1, p%NumAlf-1 + if (EqualRealNos(f_st_(Row),0.0_ReKi)) f_st_(Row+1) = 0.0_ReKi + if ( f_st_(Row+1) > f_st_(Row) ) f_st_(Row) = 0.5_ReKi * (f_st_(Row+1) + f_st_(Row-1)) + end do + do Row = iLower-1, 2, -1 + if (EqualRealNos(f_st_(Row),0.0_ReKi)) f_st_(Row-1) = 0.0_ReKi + if ( f_st_(Row-1) > f_st_(Row) ) f_st_(Row) = 0.5_ReKi * (f_st_(Row+1) + f_st_(Row-1)) + end do + + do Row = 1, p%NumAlf + if ( .not. EqualRealNos( f_st_(Row), p%Coefs(Row,ColUAf) ) ) RowIsExact(Row) = .false. + end do + p%Coefs(:,ColUAf) = f_st_ + + ! fully-separated column: invert the linear blend so the static polar is reproduced. + ! (Written for output only -- see the routine header.) + do Row = 1, p%NumAlf + if ( EqualRealNos( p%Coefs(Row,ColUAf), 1.0_ReKi ) ) then + p%Coefs(Row,col_fs) = 0.5_ReKi * cn(Row) + else + p%Coefs(Row,col_fs) = ( cn(Row) - p%Coefs(Row,col_fa) * p%Coefs(Row,ColUAf) ) & + / ( 1.0_ReKi - p%Coefs(Row,ColUAf) ) + end if + end do + + !------------------------------------------------------------------------------------- + ! Verify the Kirchhoff identity that Eq. 44 actually evaluates: + ! + ! cn = cn_fa * ( (1 + sqrt(f)) / 2 )**2 + ! + ! SCOPE AND LIMITS OF THIS CHECK -- read before relying on it. + ! + ! On any row where f was computed from the polar and NOT clipped, this identity is an + ! algebraic tautology, because f is *defined* just above as (2*sqrt(cn/cn_fa) - 1)**2: + ! + ! ((1 + sqrt(f))/2)**2 = ((1 + 2*sqrt(cn/cn_fa) - 1)/2)**2 = cn/cn_fa + ! + ! It therefore holds to round-off for ANY value of cn_fa, including a badly wrong one. + ! This check consequently does NOT validate dCNdA or alpha0 -- a 40% error in dCNdA + ! passes it cleanly (verified numerically). The implementation plan's claim that it + ! "fails immediately if the wrong attached curve is used" is only true for a curve + ! substituted *after* f was built, not for one used to build f. + ! + ! What it does still catch, and why it is kept: + ! - col_fa and ColUAf getting out of step (wrong column indices, a partial overwrite + ! by another routine, or a future edit that rebuilds one column but not the other); + ! - f being recomputed or smoothed without the corresponding cn_fa update. + ! These are real and plausible regressions, and this is a cheap guard against them. + ! + ! It is NOT a substitute for validating dCNdA/alpha0. That requires comparing the + ! reconstructed *dynamic* polar against the static one, which cannot be done here + ! because it needs the run-time CN_f path (Task 5). + ! + ! Do NOT check the linear blend f*cn_fa + (1-f)*cn_fs == cn here. That is the + ! HGM/HGMV/Oye form; IAG never evaluates it. + !------------------------------------------------------------------------------------- + if ( CheckConsistency ) then + MaxResidual = 0.0_ReKi + iWorst = 0 + nChecked = 0 + + ! scale the tolerance to the size of the polar so the test means the same thing for + ! a lightly loaded section as for a highly loaded one + ResidualTol = 1.0e-4_ReKi * max( 1.0_ReKi, maxval(abs(cn)) ) + + do Row = 1, p%NumAlf + if ( .not. RowIsExact(Row) ) cycle + + nChecked = nChecked + 1 + cn_check = p%Coefs(Row,col_fa) & + * ( 0.5_ReKi * ( 1.0_ReKi + sqrt( p%Coefs(Row,ColUAf) ) ) )**2 + Residual = abs( cn_check - cn(Row) ) + + if ( Residual > MaxResidual ) then + MaxResidual = Residual + iWorst = Row + end if + end do + + if ( iWorst > 0 .and. MaxResidual > ResidualTol ) then + ErrStat = ErrID_Fatal + ErrMsg = 'IAG model (UAMod=9): the f_st and FullyAttached columns are mutually '// & + 'inconsistent. The Kirchhoff identity cn = cn_fa*((1+sqrt(f))/2)^2 is '// & + 'violated by '//trim(num2lstr(MaxResidual))//' (tolerance '// & + trim(num2lstr(ResidualTol))//') at alpha = '// & + trim(num2lstr(p%alpha(iWorst)*R2D))//' deg, where cn = '// & + trim(num2lstr(cn(iWorst)))//', cn_fa = '//trim(num2lstr(p%Coefs(iWorst,col_fa)))// & + ', f_st = '//trim(num2lstr(p%Coefs(iWorst,ColUAf)))// & + '. One of the two columns has been overwritten or rebuilt independently '// & + 'of the other; this is an internal error, not a bad airfoil file.' + else if ( nChecked == 0 ) then + ErrStat = ErrID_Warn + ErrMsg = 'IAG model (UAMod=9): every row of the separation function was clipped or '// & + 'imposed, so the f_st/FullyAttached columns could not be cross-checked. '// & + 'This often means dCNdA or alpha0 is badly wrong for this airfoil, or the '// & + 'attached-flow region (alphaLower..alphaUpper) spans the whole table.' + end if + end if + + END SUBROUTINE ComputeIAG_SeparationFunctions +!---------------------------------------------------------------------------------------------------------------------------------- + SUBROUTINE ComputeIAG_CnMaxCnMin( p, cn, iLower, iUpper ) + ! Compute the IAG critical normal-force coefficients CnMax / CnMin: the static stall peaks + ! either side of alpha0. These are NOT Cn1/Cn2 (which are cn at f = 0.7). + ! + ! The search window matters. A global maxval over a 360-deg table can return the deep-stall + ! secondary peak near +/-45 deg, which on thin sections exceeds the stall peak; CN_CRIT would + ! then be unreachable and vortex shedding silently disabled. A window that is too narrow + ! (e.g. the +/-20 deg LimitAlphaRange) clips the genuine peak on thick sections and the vortex + ! fires far too early. Both failures are silent, so both bounds are checked below. + TYPE (AFI_Table_Type), intent(inout) :: p ! airfoil table + REAL(ReKi), intent(in ) :: cn(:) ! static normal force coefficient + INTEGER(IntKi), intent(in ) :: iLower ! lower index of the fully attached region + INTEGER(IntKi), intent(in ) :: iUpper ! upper index of the fully attached region + + REAL(ReKi) :: SearchMargin ! how far past alphaUpper/alphaLower to look + REAL(ReKi) :: PeakEdgeTol ! how close to the window edge counts as "clipped" + REAL(ReKi) :: aMaxSearch, aMinSearch + INTEGER(IntKi) :: iPeak + + PeakEdgeTol = 0.5_ReKi*D2R + + SearchMargin = 20.0_ReKi*D2R + + aMaxSearch = p%UA_BL%alphaUpper + SearchMargin + aMinSearch = p%UA_BL%alphaLower - SearchMargin + + p%UA_BL%CnMax = maxval( cn, MASK = p%alpha >= p%UA_BL%alpha0 .and. p%alpha <= aMaxSearch ) + p%UA_BL%CnMin = minval( cn, MASK = p%alpha <= p%UA_BL%alpha0 .and. p%alpha >= aMinSearch ) + + ! Sanity: detect a search window that CLIPPED the stall peak rather than containing it. + ! The symptom of clipping is that the extremum lies at the outer edge of the window, so + ! that is what we test. Do NOT test CnMax > cn(iUpper): the stall peak legitimately + ! coincides with alphaUpper on many airfoils (it is roughly where the attached region is + ! defined to end), and that comparison then falls back to Cn1 on perfectly good polars. + iPeak = maxloc( cn, DIM=1, MASK = p%alpha >= p%UA_BL%alpha0 .and. p%alpha <= aMaxSearch ) + if ( iPeak >= 1 ) then + if ( p%alpha(iPeak) >= aMaxSearch - PeakEdgeTol ) p%UA_BL%CnMax = p%UA_BL%Cn1 ! peak at window edge => clipped + end if + + iPeak = minloc( cn, DIM=1, MASK = p%alpha <= p%UA_BL%alpha0 .and. p%alpha >= aMinSearch ) + if ( iPeak >= 1 ) then + if ( p%alpha(iPeak) <= aMinSearch + PeakEdgeTol ) p%UA_BL%CnMin = p%UA_BL%Cn2 ! peak at window edge => clipped + end if + + ! Last resort: if the result is still degenerate, disable the vortex logic entirely + ! rather than leaving CnMax near zero, which would chatter the vortex on and off. + if ( p%UA_BL%CnMin >= p%UA_BL%CnMax ) then + p%UA_BL%CnMax = IAG_CnCrit_Disabled + p%UA_BL%CnMin = -IAG_CnCrit_Disabled + end if + + END SUBROUTINE ComputeIAG_CnMaxCnMin !---------------------------------------------------------------------------------------------------------------------------------- subroutine FindBoundingTables(p, secondaryDepVal, lowerTable, upperTable, xVals) @@ -1874,7 +2268,7 @@ subroutine AFI_WrHeader(delim, FileName, unOutFile, ErrStat, ErrMsg) character(*), parameter :: RoutineName = 'AFI_WrHeader' integer, parameter :: MaxLen = 17 - integer, parameter :: NumChans = 46 + integer, parameter :: NumChans = 51 character(MaxLen) :: ChanName( NumChans) character(MaxLen) :: ChanUnit( NumChans) @@ -1926,6 +2320,11 @@ subroutine AFI_WrHeader(delim, FileName, unOutFile, ErrStat, ErrMsg) ChanName(i) = 'CnBreakUpper'; ChanUnit(i) = '(-)'; i = i+1; ChanName(i) = 'alphaBreakLower'; ChanUnit(i) = '(deg)'; i = i+1; ChanName(i) = 'CnBreakLower'; ChanUnit(i) = '(-)'; i = i+1; + ChanName(i) = 'Ka'; ChanUnit(i) = '(-)'; i = i+1; + ChanName(i) = 'Kv'; ChanUnit(i) = '(-)'; i = i+1; + ChanName(i) = 'dCNdA'; ChanUnit(i) = '(-/rad)'; i = i+1; + ChanName(i) = 'CnMax'; ChanUnit(i) = '(-)'; i = i+1; + ChanName(i) = 'CnMin'; ChanUnit(i) = '(-)'; i = i+1; !$OMP critical(fileopen_critical) CALL GetNewUnit( unOutFile, ErrStat, ErrMsg ) @@ -1972,7 +2371,7 @@ subroutine AFI_WrData(k, unOutFile, delim, AFInfo) integer(IntKi) :: i integer, parameter :: MaxLen = 17 - integer, parameter :: NumChans = 46 + integer, parameter :: NumChans = 51 real(ReKi) :: TmpValues(NumChans) character(3) :: MaxLenStr character(80) :: Fmt @@ -2028,7 +2427,12 @@ subroutine AFI_WrData(k, unOutFile, delim, AFInfo) AFInfo%Table(i)%UA_BL%alphaBreakUpper*R2D, & AFInfo%Table(i)%UA_BL%CnBreakUpper , & AFInfo%Table(i)%UA_BL%alphaBreakLower*R2D, & - AFInfo%Table(i)%UA_BL%CnBreakLower + AFInfo%Table(i)%UA_BL%CnBreakLower , & + AFInfo%Table(i)%UA_BL%Ka , & + AFInfo%Table(i)%UA_BL%Kv , & + AFInfo%Table(i)%UA_BL%dCNdA , & + AFInfo%Table(i)%UA_BL%CnMax , & + AFInfo%Table(i)%UA_BL%CnMin ELSE WRITE(unOutFile, Fmt) k, i, TmpValues(3:) @@ -2117,6 +2521,8 @@ subroutine AFI_WrTables(AFI_Params,UAMod,OutRootName) end if WRITE (unOutFile,'(/,A/)') 'These predictions were generated by AirfoilInfo on '//CurDate()//' at '//CurTime()//'.' + WRITE (unOutFile,'(A)') 'The f_st, '//Prefix//'FullySep and '//trim(sFullyAtt)//' columns are built for UAMod = ' & + //trim(num2lstr(UAMod))//' and are NOT comparable across UA models.' WRITE (unOutFile,'(/,A/)') ' ' if (AFI_Params%ColUAf > 0) then From 6d27e42572abe1e70c004e15f35560a0066b42d9 Mon Sep 17 00:00:00 2001 From: Hannah Ross Date: Mon, 28 Sep 2026 15:14:24 -0600 Subject: [PATCH 03/12] Implement IAG dynamic stall ODEs and steady-state init for UA_Mod=9 --- modules/aerodyn/src/UnsteadyAero.f90 | 150 ++++++++++++++++++++++++--- 1 file changed, 137 insertions(+), 13 deletions(-) diff --git a/modules/aerodyn/src/UnsteadyAero.f90 b/modules/aerodyn/src/UnsteadyAero.f90 index 5d3cc1285f..14928a6d7e 100644 --- a/modules/aerodyn/src/UnsteadyAero.f90 +++ b/modules/aerodyn/src/UnsteadyAero.f90 @@ -750,7 +750,7 @@ subroutine UA_SetParameters( dt, InitInp, p, AFInfo, AFIndx, ErrStat, ErrMsg ) p%ShedEffect = InitInp%ShedEffect p%UA_OUTS = InitInp%UA_OUTS - if (p%UAMod==UA_HGM .or. p%UAMod==UA_HGMV .or. p%UAMod==UA_HGMV360) then + if (p%UAMod==UA_HGM .or. p%UAMod==UA_HGMV .or. p%UAMod==UA_HGMV360 .or. p%UAMod==UA_IAG) then UA_NumLinStates = 4 ! set the maximum number of states ! note: we will subtract states for nodes where UA is off for good, below @@ -889,7 +889,7 @@ subroutine UA_InitStates_Misc( p, x, xd, OtherState, m, ErrStat, ErrMsg ) ! allocate all the state arrays - if (p%UAMod == UA_HGM .or. p%UAMod == UA_HGMV .or. p%UAMod == UA_OYE .or. p%UAMod==UA_HGMV360) then + if (p%UAMod == UA_HGM .or. p%UAMod == UA_HGMV .or. p%UAMod == UA_OYE .or. p%UAMod==UA_HGMV360 .or. p%UAMod == UA_IAG) then allocate( x%element( p%nNodesPerBlade, p%numBlades ), stat=ErrStat2 ) if (ErrStat2 /= 0) call SetErrStat(ErrID_Fatal,"Cannot allocate x%x.",ErrStat,ErrMsg,RoutineName) @@ -897,6 +897,11 @@ subroutine UA_InitStates_Misc( p, x, xd, OtherState, m, ErrStat, ErrMsg ) allocate( OtherState%n(p%nNodesPerBlade, p%numBlades), stat=ErrStat2) if (ErrStat2 /= 0 ) call SetErrStat( ErrID_Fatal, " Error allocating OtherState%n.", ErrStat, ErrMsg, RoutineName) + if (p%UAMod == UA_IAG) then + allocate( OtherState%VortexOn_IAG(p%nNodesPerBlade, p%numBlades), stat=ErrStat2) + if (ErrStat2 /= 0 ) call SetErrStat( ErrID_Fatal, " Error allocating OtherState%VortexOn_IAG.", ErrStat, ErrMsg, RoutineName) + end if + if (p%UAMod == UA_HGMV) then allocate( OtherState%t_vortexBegin(p%nNodesPerBlade, p%numBlades), stat=ErrStat2) if (ErrStat2 /= 0 ) call SetErrStat( ErrID_Fatal, " Error allocating OtherState%t_vortexBegin.", ErrStat, ErrMsg, RoutineName) @@ -1017,7 +1022,7 @@ subroutine UA_ReInit( p, x, xd, OtherState, m, ErrStat, ErrMsg ) end do end do - if ( p%UAMod == UA_HGM .or. p%UAMod == UA_HGMV .or. p%UAMod == UA_OYE .or. p%UAMod==UA_HGMV360) then + if ( p%UAMod == UA_HGM .or. p%UAMod == UA_HGMV .or. p%UAMod == UA_OYE .or. p%UAMod==UA_HGMV360 .or. p%UAMod == UA_IAG) then OtherState%n = -1 ! we haven't updated OtherState%xdot, yet @@ -1042,6 +1047,10 @@ subroutine UA_ReInit( p, x, xd, OtherState, m, ErrStat, ErrMsg ) OtherState%BelowThreshold = .true. end if + if (p%UAMod == UA_IAG) then + OtherState%VortexOn_IAG = .false. + end if + elseif (p%UAMod == UA_BV) then xd%alpha_minus1 = 0.0_ReKi @@ -1489,7 +1498,7 @@ subroutine UA_ValidateInput(InitInp, ErrStat, ErrMsg) type(UA_InitInputType), intent(in ) :: InitInp ! Input data for initialization routine integer(IntKi), intent( out) :: ErrStat ! Error status of the operation character(*), intent( out) :: ErrMsg ! Error message if ErrStat /= ErrID_None - integer, parameter :: UA_VALID(8) = (/UA_None, UA_Gonzalez, UA_MinnemaPierce, UA_HGM, UA_HGMV, UA_Oye, UA_BV, UA_HGMV360/) + integer, parameter :: UA_VALID(9) = (/UA_None, UA_Gonzalez, UA_MinnemaPierce, UA_HGM, UA_HGMV, UA_Oye, UA_BV, UA_HGMV360, UA_IAG/) character(*), parameter :: RoutineName = 'UA_ValidateInput' @@ -1498,13 +1507,13 @@ subroutine UA_ValidateInput(InitInp, ErrStat, ErrMsg) if (.not.(any(InitInp%UAMod==UA_VALID))) call SetErrStat( ErrID_Fatal, & "In this version, UAMod must be 0 (None), 2 (Gonzalez's variant), 3 (Minnema/Pierce variant), 4 (continuous HGM model), 5 (HGM with vortex), & - &6 (Oye), 7 (Boeing-Vertol), or 8 (HGM-360)", ErrStat, ErrMsg, RoutineName ) ! NOTE: for later- 1 (baseline/original) + &6 (Oye), 7 (Boeing-Vertol), 8 (HGM-360), or 9 (IAG)", ErrStat, ErrMsg, RoutineName ) ! NOTE: for later- 1 (baseline/original) if (.not. InitInp%FLookUp ) call SetErrStat( ErrID_Fatal, 'FLookUp must be TRUE for this version.', ErrStat, ErrMsg, RoutineName ) if (InitInp%a_s <= 0.0) call SetErrStat ( ErrID_Fatal, 'The speed of sound (SpdSound) must be greater than zero.', ErrStat, ErrMsg, RoutineName ) - if (InitInp%UAMod == UA_HGM .or. InitInp%UAMod == UA_HGMV .or. InitInp%UAMod == UA_OYE .or. InitInp%UAMod == UA_HGMV360) then ! these are the continuous methods that integrate states + if (InitInp%UAMod == UA_HGM .or. InitInp%UAMod == UA_HGMV .or. InitInp%UAMod == UA_OYE .or. InitInp%UAMod == UA_HGMV360 .or. InitInp%UAMod == UA_IAG) then ! these are the continuous methods that integrate states if ( InitInp%IntegrationMethod /= UA_Method_RK4 & .and. InitInp%IntegrationMethod /= UA_Method_AB4 & .and. InitInp%IntegrationMethod /= UA_Method_ABM4 & @@ -1702,7 +1711,7 @@ subroutine UA_TurnOff_param(p, AFInfo, ErrStat, ErrMsg) ErrStat = ErrID_Fatal ErrMsg = 'UA parameters are not included in airfoil.' return - else if ( (p%UAMod == UA_HGM .or. p%UAMod == UA_OYE .or. p%UAMod == UA_HGMV .or. p%UAMod==UA_HGMV360) .and. & + else if ( (p%UAMod == UA_HGM .or. p%UAMod == UA_OYE .or. p%UAMod == UA_HGMV .or. p%UAMod==UA_HGMV360 .or. p%UAMod == UA_IAG) .and. & (maxval( AFInfo%Table(j)%Coefs(:, AFInfo%ColUAf) ) == 0.0_ReKi ) ) then ErrStat = ErrID_Fatal ErrMsg = 'separation function is 0 at all values.' @@ -1736,6 +1745,27 @@ subroutine UA_TurnOff_param(p, AFInfo, ErrStat, ErrMsg) elseif (p%UAMod == UA_HGMV .or. p%UAMod==UA_HGMV360) then ! pass + elseif (p%UAMod == UA_IAG) then + ! Get_alphaF's UA_IAG branch divides by dCNdA (Eq. 41), so a zero slope must turn UA off, + ! mirroring the C_lalpha treatment for UA_HGM above. + do j=1, AFInfo%NumTabs + if ( EqualRealNos(AFInfo%Table(j)%UA_BL%dCNdA, 0.0_ReKi) ) then + ErrStat = ErrID_Fatal + ErrMsg = 'dCNdA is 0.' + return + end if + end do + + ! now check about interpolated values: + do j=2, AFInfo%NumTabs + if ( sign( 1.0_ReKi, AFInfo%Table(j)%UA_BL%dCNdA) /= & + sign( 1.0_ReKi, AFInfo%Table(1)%UA_BL%dCNdA) ) then + ErrStat = ErrID_Fatal + ErrMsg = 'dCNdA (interpolated value) could be 0.' + return + end if + end do + elseif (p%UAMod == UA_Baseline .or. p%UAMod == UA_Gonzalez .or. p%UAMod == UA_MinnemaPierce) then ! unsteady aerodynamics will be turned off is Cn,alpha =0 do j=1, AFInfo%NumTabs @@ -2383,7 +2413,7 @@ subroutine UA_UpdateStates( i, j, t, n, u, uTimes, p, x, xd, OtherState, AFInfo, call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) - if (p%UAMod == UA_HGM .or. p%UAMod == UA_HGMV .or. p%UAMod == UA_OYE .or. p%UAMod == UA_HGMV360) then + if (p%UAMod == UA_HGM .or. p%UAMod == UA_HGMV .or. p%UAMod == UA_OYE .or. p%UAMod == UA_HGMV360 .or. p%UAMod == UA_IAG) then ! initialize states to steady-state values: if (OtherState%FirstPass(i,j)) then @@ -2534,7 +2564,7 @@ subroutine UA_InitStates_AllNodes( u, p, x, OtherState, AFInfo, AFIndx ) !............................................................................................................................... ! compute UA states at t=0 (with known inputs) !............................................................................................................................... - if (p%UAMod == UA_HGM .or. p%UAMod == UA_HGMV .or. p%UAMod == UA_OYE .or. p%UAMod == UA_HGMV360) then + if (p%UAMod == UA_HGM .or. p%UAMod == UA_HGMV .or. p%UAMod == UA_OYE .or. p%UAMod == UA_HGMV360 .or. p%UAMod == UA_IAG) then do j = 1,size(p%UA_off_forGood,2) ! blades do i = 1,size(p%UA_off_forGood,1) ! nodes @@ -2573,6 +2603,7 @@ SUBROUTINE HGM_Steady( i, j, u, p, x, AFInfo, ErrStat, ErrMsg ) type(AFI_UA_BL_Type) :: BL_p ! potentially interpolated UA parameters type(AFI_OutputType) :: AFI_Interp + type(AFI_OutputType) :: AFI_Interp_F ! interpolated values at alphaF (UA_IAG only; alphaF /= alphaE there) character(ErrMsgLen) :: errMsg2 integer(IntKi) :: errStat2 character(*), parameter :: RoutineName = 'HGM_Steady' @@ -2644,6 +2675,27 @@ SUBROUTINE HGM_Steady( i, j, u, p, x, AFInfo, ErrStat, ErrMsg ) ! calculate x%x(4) = fs_aF = f_st(alphaF): !call AFI_ComputeAirfoilCoefs( alphaF, u%Re, u%UserProp, AFInfo, AFI_interp, ErrStat, ErrMsg) !x%x(4) = AFI_interp%f_st + elseif (p%UAMod==UA_IAG) then + ! IAG Eq. 40 with x3dot = 0 gives x3 = CN^P. At steady state alphadot = 0, so the impulsive + ! term CN^I (Eq. 38) vanishes and CN^P = CN^C = dCNdA*sin(alphaE-alpha0) (Eq. 37). + call AddOrSub2Pi(BL_p%alpha0, alphaE) ! wrap before the sin() + x%x(3) = BL_p%dCNdA * sin(alphaE - BL_p%alpha0) ! Eqs. 37/39/40 + + ! Unlike HGM, alphaF does NOT collapse to alphaE for IAG: inverting the sinusoidal CN^C + ! through the linearized Eq. 41 gives alphaF = alpha0 + sin(alphaE-alpha0). The offset is + ! small (~0.4 deg at 20 deg from alpha0) but f_st is steep near stall, so evaluating the + ! separation function at alphaE would produce a visible start-up transient. + ! NOTE: x%x(3) MUST already be assigned -- Get_alphaF's UA_IAG branch reads it from x. + alphaF = Get_alphaF(p, u, x, BL_p, alpha_34, alphaE) ! Eq. 41 + call AFI_ComputeAirfoilCoefs( alphaF, u%Re, u%UserProp, AFInfo, AFI_interp_F, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) + if (ErrStat >= AbortErrLev) return + + ! Eq. 43 with x4dot = 0 gives x4 = f_st(alphaF). Set it here rather than falling through to + ! the shared assignment below, which uses AFI_interp (populated at alphaE) -- the wrong angle. + x%x(4) = AFI_interp_F%f_st + x%x(5) = 0.0_R8Ki ! Eq. 45, no active vortex + return else call WrScr('>>> HGM_steady logic error: should never happen.') call SetErrStat(ErrID_FATAL,"Programming error.",ErrStat,ErrMsg,RoutineName) @@ -2692,8 +2744,11 @@ subroutine UA_CalcContStateDeriv( i, j, t, u_in, p, x, OtherState, AFInfo, m, dx real(ReKi) :: alpha_34 real(ReKi) :: TuOmega real(R8Ki), parameter :: U_dot = 0.0_R8Ki ! at some point we may add this term + real(R8Ki) :: UdotTerm ! velocity-transient coefficient added to b1/b2; model-dependent, see below TYPE(UA_InputType) :: u ! Inputs at t real(R8Ki) :: CnC_dot, One_Plus_Sqrt_x4, cv_dot, CnC + real(R8Ki) :: CN_I ! IAG: impulsive (non-circulatory) normal force, Eq. 38 + real(R8Ki) :: alphaE_dot ! IAG: d(alphaE)/dt, used in the Eq. 46 vortex term ! Initialize ErrStat @@ -2745,12 +2800,31 @@ subroutine UA_CalcContStateDeriv( i, j, t, u_in, p, x, OtherState, AFInfo, m, dx ! Constraining x4 between 0 and 1 increases numerical stability (should be done elsewhere, but we'll double check here in case there were perturbations on the state value) x4 = max( min( x%x(4), 1.0_R8Ki ), 0.0_R8Ki ) + ! Velocity-transient ("added mass") coefficient added to b1/b2 in the x1/x2 ODEs below. + ! The HGM and IAG references disagree by a factor of 2 on this term: + ! HGM (Hansen et al. [40], Eqs. 8-9): c*U_dot/(2*U**2) + ! IAG (Bangga et al. 2023, Eqs. 34-35): c*V_dot/( V**2) + ! Both cite the same lineage, so one of them carries a factor-of-2 error. It cannot be + ! adjudicated from the OpenFAST source alone, so each model is implemented to match its + ! OWN reference: that keeps HGM/HGMV bit-for-bit unchanged and makes UA_Mod=9 faithful to + ! the IAG paper. Note c*V_dot/(2*V**2) = -d(Tu)/dt exactly, which makes the HGM form the + ! clean product-rule term for a time-varying Tu -- suggestive, but not decisive. + ! + ! This is currently inert because U_dot is hard-coded to 0.0 (parameter, declared above), + ! but it is written to be correct if U_dot is ever made a real input. Do NOT collapse this + ! back to a single shared expression on the grounds that "the term is zero anyway". + if (p%UAMod == UA_IAG) then + UdotTerm = p%c(i,j) * U_dot / (u%u**2) ! IAG Eqs. 34-35 + else + UdotTerm = p%c(i,j) * U_dot / (2.0_R8Ki*u%u**2) ! HGM Eqs. 8-9 [40] + end if + if (p%ShedEffect) then - if (.NOT. EqualRealNos(BL_p%A1,0.0_ReKi)) call AddOrSub2Pi(real(x%x(1)/BL_p%A1,ReKi), alpha_34) ! beause U_dot == 0, dx%x1 is A1*b1/Tu*(alpha_34 - x1/A1), we want the angle difference to be calculated correctly - dxdt%x(1) = -1.0_R8Ki / Tu * (BL_p%b1 + p%c(i,j) * U_dot/(2*u%u**2)) * x%x(1) + BL_p%b1 * BL_p%A1 / Tu * alpha_34 ! Eq. 8 [40] + if (.NOT. EqualRealNos(BL_p%A1,0.0_ReKi)) call AddOrSub2Pi(real(x%x(1)/BL_p%A1,ReKi), alpha_34) ! when U_dot == 0, dx%x1 is A1*b1/Tu*(alpha_34 - x1/A1), we want the angle difference to be calculated correctly + dxdt%x(1) = -1.0_R8Ki / Tu * (BL_p%b1 + UdotTerm) * x%x(1) + BL_p%b1 * BL_p%A1 / Tu * alpha_34 ! Eq. 8 [40] / IAG Eq. 34 - if (.NOT. EqualRealNos(BL_p%A2,0.0_ReKi)) call AddOrSub2Pi(real(x%x(2)/BL_p%A2,ReKi), alpha_34) ! beause U_dot == 0, dx%x2 is A2*b2/Tu*(alpha_34 - x2/A2), we want the angle difference to be calculated correctly - dxdt%x(2) = -1.0_R8Ki / Tu * (BL_p%b2 + p%c(i,j) * U_dot/(2*u%u**2)) * x%x(2) + BL_p%b2 * BL_p%A2 / Tu * alpha_34 ! Eq. 9 [40] + if (.NOT. EqualRealNos(BL_p%A2,0.0_ReKi)) call AddOrSub2Pi(real(x%x(2)/BL_p%A2,ReKi), alpha_34) ! when U_dot == 0, dx%x2 is A2*b2/Tu*(alpha_34 - x2/A2), we want the angle difference to be calculated correctly + dxdt%x(2) = -1.0_R8Ki / Tu * (BL_p%b2 + UdotTerm) * x%x(2) + BL_p%b2 * BL_p%A2 / Tu * alpha_34 ! Eq. 9 [40] / IAG Eq. 35 else dxdt%x(1) = 0.0_R8Ki dxdt%x(2) = 0.0_R8Ki @@ -2808,6 +2882,47 @@ subroutine UA_CalcContStateDeriv( i, j, t, u_in, p, x, OtherState, AFInfo, m, dx dxdt%x(5) = cv_dot - x%x(5)/(BL_p%T_V0 * Tu) end if + elseif (p%UAMod == UA_IAG) then + + ! x1/x2 (Eqs. 34-35) are integrated by the shared `if (p%ShedEffect)` block above -- + ! including the AddOrSub2Pi angle fix-ups and the ShedEffect=.false. zeroing -- so nothing + ! is written for them here. The one place the IAG equations differ from HGM's Eqs. 8/9 is + ! the velocity-transient coefficient (factor of 2); that is handled by UdotTerm above, + ! which selects the IAG form for UA_IAG. See the comment there. + ! + ! alphaE (Eq. 36) likewise comes from Get_HGM_constants; do not recompute it. + + ! Impulsive (non-circulatory) normal force, Eq. 38: CN^I = 4*Ka*(c/V)*alphadot. + ! c/V = 2*Tu, so this is 8*Ka*Tu*omega. TuOmega is the pre-clamped Tu*u%omega. + CN_I = 8.0_R8Ki * BL_p%Ka * TuOmega ! Eq. 38 + + call AddOrSub2Pi(BL_p%alpha0, alphaE) ! wrap before the sin() + CnC = BL_p%dCNdA * sin(alphaE - BL_p%alpha0) ! Eq. 37 + Clp = CnC + CN_I ! Eq. 39 (this is CN^P) + + ! Separated flow. NOTE: BL_p%T_p and BL_p%T_f0 were already multiplied by Tu above, so they + ! are in seconds here. T_V0 below is NOT pre-scaled and genuinely needs the *Tu. + dxdt%x(3) = ( Clp - x%x(3) ) / BL_p%T_p ! Eq. 40 + dxdt%x(4) = ( AFI_AlphaF%f_st - x4 ) / BL_p%T_f0 ! Eq. 43 + ! AFI_AlphaF was computed above at alphaF = Get_alphaF(...), i.e. Eq. 41, and f_st is the + ! IAG-consistent separation column built in AirfoilInfo::CalculateUACoeffs, clamped to [0,1]. + + ! Vortex lift, Eqs. 45-46. + if (OtherState%VortexOn_IAG(i,j)) then + One_Plus_Sqrt_x4 = 1.0_R8Ki + sqrt(x4) + + ! d(alphaE)/dt from Eq. 36 term by term, with d(alpha_34)/dt = u%omega. + ! Structurally identical to the HGMV CnC_dot construction above. + alphaE_dot = u%omega * (1.0_R8Ki - BL_p%A1 - BL_p%A2) + dxdt%x(1) + dxdt%x(2) + CnC_dot = BL_p%dCNdA * cos(alphaE - BL_p%alpha0) * alphaE_dot ! d/dt of Eq. 37 + + cv_dot = CnC_dot*(1.0_R8Ki - 0.25_R8Ki*(One_Plus_Sqrt_x4)**2) + cv_dot = cv_dot - CnC*0.25_R8Ki*One_Plus_Sqrt_x4/sqrt(max(0.0001_R8Ki,x4))*dxdt%x(4) + else + cv_dot = 0.0_R8Ki + end if + + dxdt%x(5) = cv_dot - x%x(5)/(BL_p%T_V0 * Tu) ! Eq. 45 else call WrScr('>>> UA_CalcContStateDeriv logic error: should never happen.') call SetErrStat(ErrID_FATAL,"Programming error.",ErrStat,ErrMsg,RoutineName) @@ -2878,6 +2993,15 @@ FUNCTION Get_alphaF(p, u, x, BL_p, alpha_34, alphaE_in) RESULT(alphaF) !note: BL_p%c_lalpha cannot be zero. UA is turned off at initialization if this occurs. alphaF = x%x(3)/BL_p%c_lalpha + BL_p%alpha0 ! Eq. 15 [40] + elseif (p%UAMod == UA_IAG) then + + ! IAG Eq. 41: invert the sinusoidal attached-flow relation CN^C = dCNdA*sin(alphaF-alpha0). + ! The paper linearizes the inversion, so this is algebraically the same form as the UA_HGM + ! branch above with dCNdA in place of c_lalpha -- it is NOT alpha0 + asin(x3/dCNdA). + !note: BL_p%dCNdA cannot be zero. UA is turned off at initialization if this occurs + ! (AirfoilInfo.f90 falls back to C_nalpha and then to 2*pi if the linear fit degenerates). + alphaF = x%x(3)/BL_p%dCNdA + BL_p%alpha0 ! Eq. 41 + elseif (p%UAMod==UA_HGMV360 .or. p%UAMod == UA_HGMV) then call MPi2Pi(alphaE) From e31669da86d8ddb67073734bf5632757b2956f3b Mon Sep 17 00:00:00 2001 From: Hannah Ross Date: Mon, 28 Sep 2026 16:02:58 -0600 Subject: [PATCH 04/12] Add IAG vortex shedding logic and tau_v recurrence to UA_UpdateStates --- modules/aerodyn/src/UnsteadyAero.f90 | 142 +++++++++++++++++++++++++++ 1 file changed, 142 insertions(+) diff --git a/modules/aerodyn/src/UnsteadyAero.f90 b/modules/aerodyn/src/UnsteadyAero.f90 index 14928a6d7e..52fb8fb3fe 100644 --- a/modules/aerodyn/src/UnsteadyAero.f90 +++ b/modules/aerodyn/src/UnsteadyAero.f90 @@ -900,6 +900,15 @@ subroutine UA_InitStates_Misc( p, x, xd, OtherState, m, ErrStat, ErrMsg ) if (p%UAMod == UA_IAG) then allocate( OtherState%VortexOn_IAG(p%nNodesPerBlade, p%numBlades), stat=ErrStat2) if (ErrStat2 /= 0 ) call SetErrStat( ErrID_Fatal, " Error allocating OtherState%VortexOn_IAG.", ErrStat, ErrMsg, RoutineName) + + allocate( OtherState%tau_v_IAG(p%nNodesPerBlade, p%numBlades), stat=ErrStat2) + if (ErrStat2 /= 0 ) call SetErrStat( ErrID_Fatal, " Error allocating OtherState%tau_v_IAG.", ErrStat, ErrMsg, RoutineName) + + allocate( OtherState%PositiveStall(p%nNodesPerBlade, p%numBlades), stat=ErrStat2) + if (ErrStat2 /= 0 ) call SetErrStat( ErrID_Fatal, " Error allocating OtherState%PositiveStall.", ErrStat, ErrMsg, RoutineName) + + allocate( OtherState%alpha_minus1_IAG(p%nNodesPerBlade, p%numBlades), stat=ErrStat2) + if (ErrStat2 /= 0 ) call SetErrStat( ErrID_Fatal, " Error allocating OtherState%alpha_minus1_IAG.", ErrStat, ErrMsg, RoutineName) end if if (p%UAMod == UA_HGMV) then @@ -1049,6 +1058,14 @@ subroutine UA_ReInit( p, x, xd, OtherState, m, ErrStat, ErrMsg ) if (p%UAMod == UA_IAG) then OtherState%VortexOn_IAG = .false. + OtherState%tau_v_IAG = 0.0_ReKi + OtherState%PositiveStall = .true. + ! Sentinel, NOT 0 and NOT u%alpha: UA_ReInit has no access to the inputs, so there is + ! no previous angle yet. The n > 0 guard in UA_UpdateStates skips the dAlpha test on + ! the first step and the unconditional write-back there seeds the real value. + ! Initializing to 0 instead would make the first dAlpha the full angle of attack + ! rather than a rate, which can spuriously trip the upstroke test. + OtherState%alpha_minus1_IAG = huge(1.0_ReKi) end if elseif (p%UAMod == UA_BV) then @@ -2389,8 +2406,13 @@ subroutine UA_UpdateStates( i, j, t, n, u, uTimes, p, x, xd, OtherState, AFInfo, character(*), parameter :: RoutineName = 'UA_UpdateStates' type(UA_InputType) :: u_interp_raw ! Input at current timestep, t and t+dt type(UA_InputType) :: u_interp ! Input at current timestep, t and t+dt + type(UA_InputType) :: u_vortex ! Input at t+dt, used only by the UA_IAG vortex block (see note there) type(AFI_UA_BL_Type) :: BL_p ! airfoil UA parameters retrieved in Kelvin Chain real(R8Ki) :: Tu + real(R8Ki) :: x3_IAG ! IAG: lagged normal force at t+dt + real(ReKi) :: alpha_w, alpha_p, dAlpha ! IAG: wrapped alpha, previous alpha, and their difference + real(ReKi) :: tau_v, tau_v_n1 ! IAG: non-dimensional vortex time at steps n and n-1 + logical :: PosStall, overCrit, upstroke ! IAG: vortex initiation/termination tests ! Initialize variables @@ -2522,6 +2544,103 @@ subroutine UA_UpdateStates( i, j, t, n, u, uTimes, p, x, xd, OtherState, AFInfo, end if ! p%UAMod == UA_HGMV + if (p%UAMod == UA_IAG) then + + ! Vortex initiation/termination and the non-dimensional vortex clock tau_v (Eq. 49). + ! This block MUST run after the integration above, because it tests x3 at t+dt. + ! + ! Time level: u_interp above is at t for RK4/AB4/ABM4, but the BDF2 case + ! re-extrapolates it to t+dt and does not restore it. Reusing it here would make + ! dAlpha -- which is the entire upstroke test -- depend on p%integrationMethod, an + ! integrator-dependent change in when the vortex fires. So interpolate a dedicated + ! input at an explicit t+dt, consistent with the x3(t+dt) we test against. + CALL UA_Input_ExtrapInterp( u, utimes, u_interp_raw, t+p%dt, ErrStat2, ErrMsg2 ) + CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) + IF ( ErrStat >= AbortErrLev ) RETURN + call UA_fixInputs(u_interp_raw, u_vortex, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + + ! BL_p and Tu are not otherwise in scope in this routine (they are locals of + ! UA_CalcContStateDeriv), so compute them here as the HGMV block above does. + call AFI_ComputeUACoefs( AFInfo, u_vortex%Re, u_vortex%UserProp, BL_p, ErrMsg2, ErrStat2 ) + call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) + if (ErrStat >= AbortErrLev) return + Tu = Get_Tu( u_vortex%u, p%c(i,j) ) + + x3_IAG = x%element(i,j)%x(3) ! x3 at t+dt + + ! Wrap alpha before BOTH the |alpha| > pi/2 test and the dAlpha difference: a wrap + ! event would otherwise fabricate a ~2*pi dAlpha and invert the upstroke flag. + alpha_w = u_vortex%alpha + call MPi2Pi(alpha_w) + alpha_p = OtherState%alpha_minus1_IAG(i,j) + call AddOrSub2Pi(alpha_w, alpha_p) ! put the previous angle on the same branch + dAlpha = alpha_w - alpha_p + + ! Eq. 49 is a recurrence tau_v(n) = F(tau_v(n-1)). Load the previous value once, + ! use ONLY that local on every right-hand side, and write back once at the end. + tau_v_n1 = OtherState%tau_v_IAG(i,j) + + if (.not. OtherState%VortexOn_IAG(i,j)) then + + ! No vortex running: test for initiation using the instantaneous sign of x3. + PosStall = x3_IAG >= 0.0_R8Ki + ! CN_CRIT is the static polar extremum (CnMax/CnMin), not the Cn1/Cn2 the + ! Beddoes-Leishman models use. Strict inequality: x3 lags the quasi-steady + ! normal force, and it is that lag which lets it transiently exceed the static + ! peak at all. + if (PosStall) then + overCrit = x3_IAG > BL_p%CnMax + else + overCrit = x3_IAG < BL_p%CnMin + end if + upstroke = IAG_Upstroke( alpha_w, dAlpha, PosStall ) + + ! Do not start a vortex while UA is being blended off (mirrors HGMV), and never + ! on the first step, when alpha_minus1_IAG is still the huge() sentinel from + ! UA_ReInit and dAlpha is therefore meaningless. + if (upstroke .and. overCrit .and. EqualRealNos(m%weight(i,j),1.0_ReKi) .and. n > 0) then + OtherState%VortexOn_IAG(i,j) = .true. + OtherState%PositiveStall(i,j) = PosStall ! latch for this vortex's lifetime + tau_v_n1 = 0.0_ReKi ! restart the convection clock + end if + + else + + ! Vortex running: test for termination using the LATCHED sign. Using the + ! instantaneous sign here would let CN_CRIT flip from CnMax to CnMin as x3 + ! crosses zero, inverting overCrit and killing the vortex on the spot. + PosStall = OtherState%PositiveStall(i,j) + if (PosStall) then + overCrit = x3_IAG > BL_p%CnMax + else + overCrit = x3_IAG < BL_p%CnMin + end if + upstroke = IAG_Upstroke( alpha_w, dAlpha, PosStall ) + + ! tau_v_n1 is the previous step's value; the update happens below. The vortex + ! therefore survives the step that carries tau_v past T_VL and terminates on the + ! next -- a one-step lag, consistent with HGMV's elapsed-time test above. + ! NOTE: tau_v is dimensionless and T_VL is compared to it directly. Do NOT + ! multiply T_VL by Tu here; HGMV does, because it stores a wall-clock time. + OtherState%VortexOn_IAG(i,j) = upstroke .and. overCrit .and. tau_v_n1 < BL_p%T_VL + + end if + + ! Eq. 49, advancing n-1 -> n. Get_Tu clamps Tu to [0.001,50] s, so 2*Tu is never 0. + ! 2*Tu = c/V exactly, so the paper's 0.45*dt*V/c is 0.45*dt/(2*Tu). + if (OtherState%VortexOn_IAG(i,j)) then + tau_v = tau_v_n1 + 0.45_ReKi * real(p%dt,ReKi) / real(2.0_R8Ki*Tu,ReKi) + else + tau_v = tau_v_n1 * exp( -real(p%dt,ReKi) / real(2.0_R8Ki*Tu,ReKi) ) + end if + + OtherState%tau_v_IAG(i,j) = tau_v + ! Unconditional: this also seeds the sentinel on the first step. + OtherState%alpha_minus1_IAG(i,j) = alpha_w + + end if ! p%UAMod == UA_IAG + elseif (p%UAMod == UA_BV) then ! Integrate discrete states (alpha_dot, alpha_filt_minus1) call UA_UpdateDiscOtherState_BV( i, j, u_interp, p, xd, OtherState, AFInfo, m, ErrStat2, ErrMsg2 ) @@ -3037,6 +3156,29 @@ FUNCTION Get_alphaF(p, u, x, BL_p, alpha_34, alphaE_in) RESULT(alphaF) END FUNCTION Get_alphaF !--------------------------------------------------------------------------------- +!> IAG model: is the airfoil moving further into stall of the given sign? +!! Factored out because the caller must evaluate it twice -- once against the instantaneous +!! sign of x3 when testing for vortex initiation, and once against the sign LATCHED at +!! initiation when testing for termination. The two cannot be hoisted into a single call: +!! using the instantaneous sign during a vortex's lifetime would let CN_CRIT flip from CnMax +!! to CnMin as x3 crosses zero, inverting the test and terminating the vortex immediately. +!! +!! Beyond +/-90 deg the airfoil is in reverse flow and the sense of "increasing incidence" +!! inverts, hence the branch on abs(alpha). alpha is expected to be wrapped to [-pi,pi] +!! already. +pure LOGICAL FUNCTION IAG_Upstroke( alpha, dAlpha, PosStall ) + REAL(ReKi), INTENT(IN ) :: alpha !< angle of attack, already wrapped to [-pi,pi] + REAL(ReKi), INTENT(IN ) :: dAlpha !< change in alpha over the last step + LOGICAL, INTENT(IN ) :: PosStall !< .true. for positive stall (Cn >= 0) + + if (abs(alpha) <= PiBy2) then + IAG_Upstroke = ( PosStall .eqv. (dAlpha > 0.0_ReKi) ) + else + IAG_Upstroke = ( PosStall .neqv. (dAlpha > 0.0_ReKi) ) + end if + +END FUNCTION IAG_Upstroke +!--------------------------------------------------------------------------------- !> Compute angle of attack at 3/4 chord point based on values at Aerodynamic center real(ReKi) function Get_Alpha34(v_ac, omega, d_34_to_ac) real(ReKi), intent(in) :: v_ac(2) !< Velocity at aerodynamic center (AC) From e8b4a6b6b85c2f0a8ce69a9abbb4bc39bdc14345 Mon Sep 17 00:00:00 2001 From: Hannah Ross Date: Mon, 28 Sep 2026 16:47:15 -0600 Subject: [PATCH 05/12] Add IAG (UA_Mod=9) output branch to UnsteadyAero --- modules/aerodyn/src/UnsteadyAero.f90 | 184 ++++++++++++++++++++++++++- 1 file changed, 180 insertions(+), 4 deletions(-) diff --git a/modules/aerodyn/src/UnsteadyAero.f90 b/modules/aerodyn/src/UnsteadyAero.f90 index 52fb8fb3fe..02b39334b7 100644 --- a/modules/aerodyn/src/UnsteadyAero.f90 +++ b/modules/aerodyn/src/UnsteadyAero.f90 @@ -63,6 +63,14 @@ module UnsteadyAero real(ReKi), parameter, public :: UA_u_min = 0.01_ReKi ! m/s; used to provide a minimum value so UA equations don't blow up (this should be much lower than range where UA is turned off) real(ReKi), parameter :: K1pos=1.0_ReKi, K1neg=0.5_ReKi ! K1 coefficients for BV model real(ReKi), parameter :: MaxTuOmega = 1.5_ReKi ! adding a little safety factor for UA models + ! IAG model blend limits, in degrees of |alpha|. The dynamic Cd/Cm are faded to static + ! between BlendLo and BlendHi, and the vortex state x5 is faded out between x5FadeLo and + ! x5FadeHi. Both come from the IAG theory text rather than the paper's equations, so they + ! are named constants here rather than user inputs -- see the notes on Task 5. + real(ReKi), parameter :: IAG_BlendLo = 30.0_ReKi ! deg; start of the dynamic->static blend + real(ReKi), parameter :: IAG_BlendHi = 45.0_ReKi ! deg; end of the dynamic->static blend + real(ReKi), parameter :: IAG_x5FadeLo = 45.0_ReKi ! deg; start of the x5 fade-out + real(ReKi), parameter :: IAG_x5FadeHi = 75.0_ReKi ! deg; end of the x5 fade-out contains @@ -1226,6 +1234,8 @@ subroutine UA_Init_Outputs(InitInp, p, y, InitOut, errStat, errMsg) p%NumOuts = 21 elseif(p%UAMod == UA_HGMV360) then p%NumOuts = 22 + elseif(p%UAMod == UA_IAG) then + p%NumOuts = 24 ! channels 1-20 as HGM, plus x5, tau_v, alphaF, VortexOn elseif(p%UAMod == UA_BV) then p%NumOuts = 26 else @@ -1287,7 +1297,7 @@ subroutine UA_Init_Outputs(InitInp, p, y, InitOut, errStat, errMsg) iOffAcc = iOffset+11 - elseif (p%UAmod == UA_HGM .or. p%UAMod == UA_HGMV .or. p%UAMod == UA_OYE .or. p%UAMod == UA_HGMV360) then + elseif (p%UAmod == UA_HGM .or. p%UAMod == UA_HGMV .or. p%UAMod == UA_OYE .or. p%UAMod == UA_HGMV360 .or. p%UAMod == UA_IAG) then InitOut%WriteOutputHdr(iOffset+ 8) = trim(chanPrefix)//'omega' InitOut%WriteOutputHdr(iOffset+ 9) = trim(chanPrefix)//'alphaE' @@ -1331,6 +1341,16 @@ subroutine UA_Init_Outputs(InitInp, p, y, InitOut, errStat, errMsg) InitOut%WriteOutputUnt(iOffset+21) = '(-)' InitOut%WriteOutputUnt(iOffset+22) = '(-)' iOffAcc = iOffset+22 + else if (p%UAmod == UA_IAG) then + InitOut%WriteOutputHdr(iOffset+21) = trim(chanPrefix)//'x5' + InitOut%WriteOutputHdr(iOffset+22) = trim(chanPrefix)//'tau_v' + InitOut%WriteOutputHdr(iOffset+23) = trim(chanPrefix)//'alphaF' + InitOut%WriteOutputHdr(iOffset+24) = trim(chanPrefix)//'VortexOn' + InitOut%WriteOutputUnt(iOffset+21) = '(-)' + InitOut%WriteOutputUnt(iOffset+22) = '(-)' + InitOut%WriteOutputUnt(iOffset+23) = '(deg)' + InitOut%WriteOutputUnt(iOffset+24) = '(-)' + iOffAcc = iOffset+24 end if elseif(p%UAMod == UA_BV) then @@ -2623,7 +2643,11 @@ subroutine UA_UpdateStates( i, j, t, n, u, uTimes, p, x, xd, OtherState, AFInfo, ! next -- a one-step lag, consistent with HGMV's elapsed-time test above. ! NOTE: tau_v is dimensionless and T_VL is compared to it directly. Do NOT ! multiply T_VL by Tu here; HGMV does, because it stores a wall-clock time. - OtherState%VortexOn_IAG(i,j) = upstroke .and. overCrit .and. tau_v_n1 < BL_p%T_VL + ! The weight test also terminates a running vortex once UA starts blending off: + ! a vortex left running through the cutout blend would keep tau_v advancing + ! against loads that are being faded to static. + OtherState%VortexOn_IAG(i,j) = upstroke .and. overCrit .and. tau_v_n1 < BL_p%T_VL & + .and. EqualRealNos(m%weight(i,j),1.0_ReKi) end if @@ -3179,6 +3203,60 @@ pure LOGICAL FUNCTION IAG_Upstroke( alpha, dAlpha, PosStall ) END FUNCTION IAG_Upstroke !--------------------------------------------------------------------------------- +!> Normal-force coefficient from an airfoil-table interpolation. Consistent with +!! AirfoilInfo::Calculate_Cn and with the inlined expression in Get_f_from_Lookup. +!! Cd0 is subtracted so that Cn -> 0 at alpha0. +!! alpha need not be wrapped: cos/sin are 2*pi-periodic and the table lookup wraps +!! internally. Do not add a wrap here. +!! NOTE: AFI_interp%Cd0 is forced to 0 for tables without UA data, so this is only +!! meaningful for UA-enabled tables -- always true on the IAG code path. +pure real(ReKi) function Get_Cn( AFI_interp, alpha ) + type(AFI_OutputType), intent(in) :: AFI_interp + real(ReKi), intent(in) :: alpha + Get_Cn = AFI_interp%Cl*cos(alpha) + (AFI_interp%Cd - AFI_interp%Cd0)*sin(alpha) +end function Get_Cn +!--------------------------------------------------------------------------------- +!> Chordwise-force coefficient from an airfoil-table interpolation, consistent with the +!! inlined expression in Get_f_c_from_Lookup. Cd0 is subtracted for the same reason as +!! in Get_Cn. +!! NOTE the module is not self-consistent about this: the HGM/HGMV output paths compute +!! y%Cc = y%Cl*SinAlpha - y%Cd*CosAlpha *without* subtracting Cd0. IAG's Eq. 51 defines +!! CT_D as the viscous static chordwise force, which is the Cd0-subtracted quantity. +pure real(ReKi) function Get_Cc( AFI_interp, alpha ) + type(AFI_OutputType), intent(in) :: AFI_interp + real(ReKi), intent(in) :: alpha + Get_Cc = AFI_interp%Cl*sin(alpha) - (AFI_interp%Cd - AFI_interp%Cd0)*cos(alpha) +end function Get_Cc +!--------------------------------------------------------------------------------- +!> IAG model: linearly blend the dynamic Cd and Cm back to their static values between +!! IAG_BlendLo and IAG_BlendHi degrees of incidence. +!! Past ~30 deg the model's separated-flow construction loses validity, so the dynamic +!! corrections are faded out rather than extrapolated. Cl/Cn/Cc are deliberately NOT +!! blended here: Cl is reconstructed from Cn and Cc by Eq. 52, so blending it separately +!! would make the three mutually inconsistent. +!! This is applied BEFORE the shared UA_BlendSteady so that the two blends compose: this +!! one fades dynamic->static, that one fades UA->steady at the cutout. +subroutine IAG_BlendStatic( alpha_w, AFI_static, AFInfo, y ) + real(ReKi), intent(in ) :: alpha_w !< angle of attack, wrapped to [-pi,pi] + type(AFI_OutputType), intent(in ) :: AFI_static !< static interpolation at u%alpha + type(AFI_ParameterType), intent(in ) :: AFInfo !< airfoil parameters (for ColCm) + type(UA_OutputType), intent(inout) :: y !< outputs, modified in place + + real(ReKi) :: w + + ! w = 1 fully dynamic, w = 0 fully static + w = 1.0_ReKi - min( max( (abs(alpha_w)*R2D - IAG_BlendLo)/(IAG_BlendHi - IAG_BlendLo), & + 0.0_ReKi ), 1.0_ReKi ) + + if (w < 1.0_ReKi) then + y%Cd = w*y%Cd + (1.0_ReKi - w)*AFI_static%Cd + if (AFInfo%ColCm /= 0) then + y%Cm = w*y%Cm + (1.0_ReKi - w)*AFI_static%Cm + end if + end if + +end subroutine IAG_BlendStatic +!--------------------------------------------------------------------------------- !> Compute angle of attack at 3/4 chord point based on values at Aerodynamic center real(ReKi) function Get_Alpha34(v_ac, omega, d_34_to_ac) real(ReKi), intent(in) :: v_ac(2) !< Velocity at aerodynamic center (AC) @@ -3773,6 +3851,14 @@ subroutine UA_CalcOutput( i, j, t, u_in, p, x, xd, OtherState, AFInfo, y, misc, type(AFI_OutputType) :: AFI_interp type(AFI_OutputType) :: AFI_interpE type(AFI_OutputType) :: AFI_interpF + ! for UA_IAG + type(AFI_OutputType) :: AFI_interpA ! static interpolation at u%alpha + real(ReKi) :: alpha_w ! u%alpha wrapped to [-pi,pi] + real(ReKi) :: tau_v ! non-dimensional vortex time + real(ReKi) :: CN_I, CN_C, CN_f, CN_D, CT_D, CM_C, CPv + real(ReKi) :: One_Plus_Sqrt_x4 + real(ReKi) :: fs_alpha ! separation function at u%alpha + real(ReKi) :: x5_fade, x5_eff ErrStat = ErrID_None ! no error has occurred @@ -3859,7 +3945,7 @@ subroutine UA_CalcOutput( i, j, t, u_in, p, x, xd, OtherState, AFInfo, y, misc, call BV_CalcOutput() if (ErrStat >= AbortErrLev) return - elseif (p%UAMod == UA_HGM .or. p%UAMod == UA_HGMV .or. p%UAMod == UA_OYE .or. p%UAMod == UA_HGMV360) then + elseif (p%UAMod == UA_HGM .or. p%UAMod == UA_HGMV .or. p%UAMod == UA_OYE .or. p%UAMod == UA_HGMV360 .or. p%UAMod == UA_IAG) then ! --- CalcOutput State Space models x_in = x%element(i,j) @@ -3879,6 +3965,16 @@ subroutine UA_CalcOutput( i, j, t, u_in, p, x, xd, OtherState, AFInfo, y, misc, call AFI_ComputeAirfoilCoefs( alphaE, u%Re, u%UserProp, AFInfo, AFI_interpE, ErrStat2, ErrMsg2 ) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + ! IAG's Eq. 53 and its static blend both need the STATIC values at u%alpha, which the + ! other HGM-family models never look up (they work from alphaE). Only done for IAG so + ! the other models keep their exact current cost. + if (p%UAMod == UA_IAG) then + call AFI_ComputeAirfoilCoefs( u%alpha, u%Re, u%UserProp, AFInfo, AFI_interpA, ErrStat2, ErrMsg2 ) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + if (ErrStat >= AbortErrLev) return + fs_alpha = AFI_interpA%f_st + end if + ! Constraining x4 between 0 and 1 increases numerical stability (should be done elsewhere, but we'll double check here in case there were perturbations on the state value) x4 = max( min( x_in%x(4), 1.0_R8Ki ), 0.0_R8Ki ) @@ -3973,6 +4069,79 @@ subroutine UA_CalcOutput( i, j, t, u_in, p, x, xd, OtherState, AFInfo, y, misc, y%Cm = AFI_interpE%Cm + cn_circ * delta_c_mf_primeprime - 0.0_ReKi * piBy2 * TuOmega - 0.25_ReKi*(1.0_ReKi - cos(pi * tV_ratio ))*x5 end if + elseif (p%UAMod == UA_IAG) then + + ! None of alphaF, AFI_interpF, x5, CN_I or CN_C is in scope on this path: alphaF is + ! only computed inside the HGMV360 sub-branch, x5 only inside the HGMV sub-branch, + ! and CN_I/CN_C are UA_CalcContStateDeriv locals. All are formed here, using the + ! SAME expressions as the state-derivative branch so the outputs cannot drift from + ! the ODEs. BL_p and the clamped TuOmega are already available above. + x5 = x_in%x(5) ! vortex normal force, Eq. 45 + tau_v = OtherState%tau_v_IAG(i,j) ! written by UA_UpdateStates + + alpha_w = u%alpha + call MPi2Pi(alpha_w) ! for the CM_C and blend tests + + alphaF = Get_alphaF(p, u, x_in, BL_p, alpha_34, alphaE) ! Eq. 41 + call AFI_ComputeAirfoilCoefs( alphaF, u%Re, u%UserProp, AFInfo, AFI_interpF, ErrStat2, ErrMsg2 ) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + if (ErrStat >= AbortErrLev) return + + call AddOrSub2Pi(BL_p%alpha0, alphaE) ! wrap before the sin() + CN_I = 8.0_ReKi * BL_p%Ka * real(TuOmega,ReKi) ! Eq. 38 (clamped TuOmega) + CN_C = BL_p%dCNdA * sin(alphaE - BL_p%alpha0) ! Eq. 37 + + ! Kirchhoff blend toward the SINUSOIDAL attached curve -- the same curve the + ! IAG separation column was built against in AirfoilInfo (plan section 3.2). + One_Plus_Sqrt_x4 = 1.0_ReKi + sqrt(x4) + CN_f = BL_p%dCNdA * (0.5_ReKi*One_Plus_Sqrt_x4)**2 * sin(alphaE - BL_p%alpha0) + CN_I ! Eq. 44 + + ! x5 is faded out between 45 and 75 deg (theory text): past deep stall the vortex + ! construction is no longer meaningful. Applied to every x5 contribution below. + x5_fade = 1.0_ReKi - BlendCosine( abs(alpha_w)*R2D, IAG_x5FadeLo, IAG_x5FadeHi ) + x5_eff = x5 * x5_fade + + CN_D = CN_f + x5_eff ! Eq. 50 (CN_V = x5) + CT_D = Get_Cc( AFI_interpF, alphaF ) ! Eq. 51, viscous Cc at alphaF + + y%Cn = CN_D + y%Cc = CT_D + y%Cl = CN_D*CosAlpha - CT_D*SinAlpha ! Eq. 52 + + ! Eq. 53. AFI_interp is the static interpolation at u%alpha, made further below for + ! the other models; the IAG branch needs it here, so it is computed above. + delta_c_df_primeprime = (0.5_ReKi*(1.0_ReKi - sqrt(x4)))**2 & + - (0.5_ReKi*(1.0_ReKi - sqrt(fs_alpha)))**2 + call AddOrSub2Pi(u%alpha, alphaE) + y%Cd = AFI_interpA%Cd + (u%alpha - alphaE)*CN_C & + + (AFI_interpA%Cd - BL_p%Cd0)*delta_c_df_primeprime & + + x5_eff*SinAlpha + + if (AFInfo%ColCm == 0) then ! we don't have a cm column, so make everything 0 + y%Cm = 0.0_ReKi + else + ! Eq. 55. The sign flips when backwinded. Written as an explicit test on the + ! WRAPPED alpha rather than SIGN(): Fortran's SIGN(x,0.0) returns +|x|, and an + ! unwrapped angle would make abs(alpha) > PiBy2 true for e.g. 2*pi-0.1, which is + ! physically a small positive incidence. + ! The non-backwinded branch matches the -piBy2*TuOmega added-mass term used by + ! HGM and HGMV above; the sign flip is the IAG-specific part. + if (abs(alpha_w) > PiBy2) then + CM_C = PiBy2 * real(TuOmega,ReKi) + else + CM_C = -PiBy2 * real(TuOmega,ReKi) + end if + + ! Eq. 31: vortex center-of-pressure travel. tau_v is dimensionless, as is T_VL. + ! Clamped so the cosine cannot wrap back down if tau_v overshoots T_VL. + CPv = BL_p%Kv * (1.0_ReKi - cos( pi * min(tau_v/BL_p%T_VL, 1.0_ReKi) )) + y%Cm = AFI_interpF%Cm - CPv*x5_eff + CM_C ! Eqs. 54/29/30 + end if + + ! IAG-specific linear blend back to static Cd/Cm between 30 and 45 deg, applied + ! BEFORE the shared UA_BlendSteady below so the two blends compose. + call IAG_BlendStatic( alpha_w, AFI_interpA, AFInfo, y ) + else call SetErrStat(ErrID_Fatal, "Programming error, UAMod continuous model not accounted for", ErrStat, ErrMsg, RoutineName) @@ -4163,7 +4332,7 @@ subroutine CalcWriteOutputs() y%WriteOutput(iOffset+11) = alpha_34*R2D iOffAcc = iOffset+11 - elseif (p%UAMod == UA_HGM .or. p%UAMod == UA_HGMV .or. p%UAMod == UA_OYE .or. p%UAMod == UA_HGMV360) then + elseif (p%UAMod == UA_HGM .or. p%UAMod == UA_HGMV .or. p%UAMod == UA_OYE .or. p%UAMod == UA_HGMV360 .or. p%UAMod == UA_IAG) then y%WriteOutput(iOffset+ 8) = u%omega*R2D y%WriteOutput(iOffset+ 9) = alphaE*R2D y%WriteOutput(iOffset+10) = Tu @@ -4187,6 +4356,13 @@ subroutine CalcWriteOutputs() y%WriteOutput(iOffset+21) = x_in%x(6) !x%element(i,j)%x(6) y%WriteOutput(iOffset+22) = x_in%x(7) !x%element(i,j)%x(7) iOffAcc = iOffset+22 + else if (p%UAMod == UA_IAG) then + y%WriteOutput(iOffset+21) = x_in%x(5) + y%WriteOutput(iOffset+22) = OtherState%tau_v_IAG(i,j) + y%WriteOutput(iOffset+23) = alphaF*R2D + ! LOGICAL has no implicit conversion to real in Fortran; use merge. + y%WriteOutput(iOffset+24) = merge(1.0_ReKi, 0.0_ReKi, OtherState%VortexOn_IAG(i,j)) + iOffAcc = iOffset+24 end if elseif(p%UAMod == UA_BV) then From 4f2ad58564c8c2d530192dc11f23a281b24e8a0b Mon Sep 17 00:00:00 2001 From: Hannah Ross Date: Mon, 28 Sep 2026 17:16:50 -0600 Subject: [PATCH 06/12] Add IAG (UA_Mod=9) validation checks to UA_ValidateAFI --- modules/aerodyn/src/UnsteadyAero.f90 | 74 +++++++++++++++++++++------- 1 file changed, 56 insertions(+), 18 deletions(-) diff --git a/modules/aerodyn/src/UnsteadyAero.f90 b/modules/aerodyn/src/UnsteadyAero.f90 index 02b39334b7..317ed00d5d 100644 --- a/modules/aerodyn/src/UnsteadyAero.f90 +++ b/modules/aerodyn/src/UnsteadyAero.f90 @@ -1590,9 +1590,12 @@ subroutine UA_ValidateAFI(UAMod, FLookup, AFInfo, ErrStat, ErrMsg) if ( tab%InclUAdata ) then ! parameters used only for UAMod/=UA_HGM) - if (UAMod == UA_Baseline .or. UAMod == UA_Gonzalez .or. UAMod == UA_MinnemaPierce .or. UAMod == UA_HGMV) then + ! NOTE: UA_IAG is included here for the T_VL/T_V0 checks below, which it needs + ! (both appear in Eqs. 45 and 31). The inner blocks that do NOT apply to it are + ! excluded individually, the same way UA_HGMV is. + if (UAMod == UA_Baseline .or. UAMod == UA_Gonzalez .or. UAMod == UA_MinnemaPierce .or. UAMod == UA_HGMV .or. UAMod == UA_IAG) then - if (UAMod /= UA_HGMV) then + if (UAMod /= UA_HGMV .and. UAMod /= UA_IAG) then if ( EqualRealNos(tab%UA_BL%St_sh, 0.0_ReKi) ) then call SetErrStat(ErrID_Fatal, 'UA St_sh parameter must not be 0.', ErrStat_tab, ErrMsg_tab, "" ) end if @@ -1632,36 +1635,71 @@ subroutine UA_ValidateAFI(UAMod, FLookup, AFInfo, ErrStat, ErrMsg) call SetErrStat(ErrID_Fatal, 'UA T_V0 parameter must be greater than 0.', ErrStat_tab, ErrMsg_tab, "" ) end if - if (tab%UA_BL%Cn2 >= tab%UA_BL%Cn1) call SetErrStat(ErrID_Fatal, 'Cn2 must be less than Cn1.', ErrStat_tab, ErrMsg_tab, "" ) + ! IAG uses CnMax/CnMin as its vortex-onset thresholds and never reads Cn1/Cn2, + ! which are only computed on the non-circular-polar path and may legitimately + ! be absent. Requiring Cn2 < Cn1 would reject otherwise-valid tables. + if (UAMod == UA_IAG) then + if (tab%UA_BL%CnMin >= tab%UA_BL%CnMax) call SetErrStat(ErrID_Fatal, 'CnMin must be less than CnMax.', ErrStat_tab, ErrMsg_tab, "" ) + else + if (tab%UA_BL%Cn2 >= tab%UA_BL%Cn1) call SetErrStat(ErrID_Fatal, 'Cn2 must be less than Cn1.', ErrStat_tab, ErrMsg_tab, "" ) + end if end if + ! IAG-specific parameters (Eqs. 38, 31 and 37/41/44 respectively). + if (UAMod == UA_IAG) then + ! Ka scales the impulsive normal force CN^I; a negative value would invert it. + if ( tab%UA_BL%Ka < 0.0_ReKi ) then + call SetErrStat(ErrID_Fatal, 'UA Ka parameter must not be negative.', ErrStat_tab, ErrMsg_tab, "" ) + end if + ! Kv scales the vortex center-of-pressure travel CPv. + if ( tab%UA_BL%Kv < 0.0_ReKi ) then + call SetErrStat(ErrID_Fatal, 'UA Kv parameter must not be negative.', ErrStat_tab, ErrMsg_tab, "" ) + end if + ! Get_alphaF's UA_IAG branch divides by dCNdA (Eq. 41), and Eq. 44 multiplies + ! by it. AirfoilInfo falls back to 2*pi if its own fit degenerates, so a value + ! of 0 here can only come from the user's table. + if ( tab%UA_BL%dCNdA <= 0.0_ReKi ) then + call SetErrStat(ErrID_Fatal, 'UA dCNdA parameter must be greater than 0.', ErrStat_tab, ErrMsg_tab, "" ) + end if + end if + if (UAMod /= UA_HGMV) then if ( tab%UA_BL%alpha0 > pi .or. tab%UA_BL%alpha0 < -pi ) then call SetErrStat(ErrID_Fatal, 'UA alpha0 parameter must be between -180 and 180 degrees.', ErrStat_tab, ErrMsg_tab, "" ) end if end if ! Not UA_HGM - if ( tab%UA_BL%UACutout < Pi .and. (UAMod == UA_HGM .or. UAMod == UA_HGMV .or. UAMod == UA_OYE .or. UAMod == UA_HGMV360) ) then + if ( tab%UA_BL%UACutout < Pi .and. (UAMod == UA_HGM .or. UAMod == UA_HGMV .or. UAMod == UA_OYE .or. UAMod == UA_HGMV360 .or. UAMod == UA_IAG) ) then cl_fs = InterpStp( tab%UA_BL%UACutout, tab%alpha, tab%Coefs(:,AFInfo%ColUAf), indx, tab%NumAlf ) if (.not. EqualRealNos( cl_fs, 0.0_ReKi ) ) then call SetErrStat(ErrID_Severe, 'UA cutout parameter should be at a value where the separation function is 0;'// & ' separation function is '//trim(num2lstr(cl_fs))//'.' , ErrStat_tab, ErrMsg_tab, "" ) end if - ! C_alpha should have a reasonable value - if (abs(tab%UA_BL%C_lalpha)>9.11_ReKi) then ! 45% above 2*pi, arbitrary.. - call SetErrStat(ErrID_Severe, 'Large value of C_lalpha.'// & - " C_lalpha="//trim(num2lstr(tab%UA_BL%C_lalpha))//& - ". We advise to check this value or provide it in the input file.", ErrStat_tab, ErrMsg_tab, "" ) - endif - ! NOTE: check if C_nalpha is alwasy defined - ! C_lalpha and C_nalpha should be in the same ballpark - if (abs(tab%UA_BL%C_nalpha-tab%UA_BL%C_lalpha)>3.0_ReKi) then ! arbitrary criteria.. - call SetErrStat(ErrID_Severe, 'Large difference between C_lalpha and C_nalpha.'// & - " C_lalpha="//trim(num2lstr(tab%UA_BL%C_lalpha))//& - " C_nalpha="//trim(num2lstr(tab%UA_BL%C_nalpha))//& - ". We advise to check these values or provide them in the input file.", ErrStat_tab, ErrMsg_tab, "" ) - endif + + ! The C_lalpha/C_nalpha sanity checks below do NOT apply to IAG: it blends on + ! dCNdA (Eq. 44) and never reads either of them, so a table tuned for IAG can + ! legitimately carry values that trip these warnings. dCNdA is range-checked + ! separately above. The f_st checks that follow DO apply -- they validate the + ! ColUAf column, which IAG builds itself in AirfoilInfo -- so they are outside + ! this exclusion. + if (UAMod /= UA_IAG) then + ! C_alpha should have a reasonable value + if (abs(tab%UA_BL%C_lalpha)>9.11_ReKi) then ! 45% above 2*pi, arbitrary.. + call SetErrStat(ErrID_Severe, 'Large value of C_lalpha.'// & + " C_lalpha="//trim(num2lstr(tab%UA_BL%C_lalpha))//& + ". We advise to check this value or provide it in the input file.", ErrStat_tab, ErrMsg_tab, "" ) + endif + ! NOTE: check if C_nalpha is alwasy defined + ! C_lalpha and C_nalpha should be in the same ballpark + if (abs(tab%UA_BL%C_nalpha-tab%UA_BL%C_lalpha)>3.0_ReKi) then ! arbitrary criteria.. + call SetErrStat(ErrID_Severe, 'Large difference between C_lalpha and C_nalpha.'// & + " C_lalpha="//trim(num2lstr(tab%UA_BL%C_lalpha))//& + " C_nalpha="//trim(num2lstr(tab%UA_BL%C_nalpha))//& + ". We advise to check these values or provide them in the input file.", ErrStat_tab, ErrMsg_tab, "" ) + endif + end if + vmax = maxval(tab%Coefs(:,AFInfo%ColUAf)) if (vmax>1.00_ReKi) then call SetErrStat(ErrID_Severe, 'The separation function f_st exceeds 1;'// & From fe6f8616d3d58012d1e40dddfaa699d82d27c6b4 Mon Sep 17 00:00:00 2001 From: Hannah Ross Date: Tue, 29 Sep 2026 09:28:30 -0600 Subject: [PATCH 07/12] Enable UA_Mod=9 (IAG) for linearization and add it to AD_PrintSum --- modules/aerodyn/src/AeroDyn.f90 | 9 +++++++-- modules/aerodyn/src/AeroDyn_IO.f90 | 4 +++- 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/modules/aerodyn/src/AeroDyn.f90 b/modules/aerodyn/src/AeroDyn.f90 index 9da727a7bd..aebac233cd 100644 --- a/modules/aerodyn/src/AeroDyn.f90 +++ b/modules/aerodyn/src/AeroDyn.f90 @@ -4894,8 +4894,13 @@ SUBROUTINE ValidateInputData( InitInp, InputFileData, NumBl, calcCrvAngle, ErrSt call SetErrStat( ErrID_Fatal, 'Wake_Mod must be 0 or 1 for linearization.', ErrStat, ErrMsg, RoutineName ) endif - if (InputFileData%UA_Init%UAMod /= UA_None .and. InputFileData%UA_Init%UAMod /= UA_HGM .and. InputFileData%UA_Init%UAMod /= UA_HGMV .and. InputFileData%UA_Init%UAMod /= UA_OYE) then - call SetErrStat( ErrID_Fatal, 'UA_Mod must be 0, 4, 5, or 6 for linearization.', ErrStat, ErrMsg, RoutineName ) + ! NOTE: this is a chain of exclusions, not an allowed-list, so a new model is + ! REJECTED for linearization unless it is named here. UA_IAG is admitted on the + ! same basis as UA_HGMV: both carry a 5th (vortex) state that is not linearized, + ! and both linearize x1-x4 only (see p%lin_nx in UA_SetParameters). UA_HGMV360 is + ! deliberately still absent. + if (InputFileData%UA_Init%UAMod /= UA_None .and. InputFileData%UA_Init%UAMod /= UA_HGM .and. InputFileData%UA_Init%UAMod /= UA_HGMV .and. InputFileData%UA_Init%UAMod /= UA_OYE .and. InputFileData%UA_Init%UAMod /= UA_IAG) then + call SetErrStat( ErrID_Fatal, 'UA_Mod must be 0, 4, 5, 6, or 9 for linearization.', ErrStat, ErrMsg, RoutineName ) end if select case(InputFileData%DBEMT_Mod) diff --git a/modules/aerodyn/src/AeroDyn_IO.f90 b/modules/aerodyn/src/AeroDyn_IO.f90 index 8493b31a57..b98bead7a8 100644 --- a/modules/aerodyn/src/AeroDyn_IO.f90 +++ b/modules/aerodyn/src/AeroDyn_IO.f90 @@ -981,7 +981,7 @@ SUBROUTINE ParsePrimaryFileInfo( PriPath, InitInp, InputFile, RootName, NumBlade ! UAMod (Legacy) call ParseVar( FileInfo_In, CurLine, "UAMod", UAMod_Old, ErrStat2, ErrMsg2, UnEc ) UAModProvided = legacyInputPresent('UAMod', CurLine, ErrStat2, ErrMsg2, 'UA_Mod=0 (AFAeroMod=1), UA_Mod>1 (AFAeroMod=2 and UA_Mod=UAMod') - ! UA_Mod - Unsteady Aero Model Switch (switch) {0=Quasi-steady (no UA), 2=Gonzalez's variant (changes in Cn,Cc,Cm), 3=Minnema/Pierce variant (changes in Cc and Cm)} + ! UA_Mod - Unsteady Aero Model Switch (switch) {0=Quasi-steady (no UA), 2=Gonzalez's variant (changes in Cn,Cc,Cm), 3=Minnema/Pierce variant (changes in Cc and Cm), 4=HGM, 5=HGM+vortex, 6=Oye, 7=Boeing-Vertol, 9=IAG} call ParseVar( FileInfo_In, CurLine, "UA_Mod", InputFileData%UA_Init%UAMod, ErrStat2, ErrMsg2, UnEc ) if (newInputMissing('UA_Mod', CurLine, errStat2, errMsg2)) then ! We'll deal with it when we deal with AFAeroMod @@ -2089,6 +2089,8 @@ SUBROUTINE AD_PrintSum( InputFileData, p, p_AD, u, y, NumBlades, BladeInputFileD Msg = 'Stieg Oye dynamic stall model' case (UA_BV) Msg = 'Boeing-Vertol dynamic stall model (e.g. used in CACTUS)' + case (UA_IAG) + Msg = 'IAG dynamic stall model (first-order, state-space)' case default Msg = 'unknown' end select From c438b386af2da2828030c096a6b103b88b45ce8e Mon Sep 17 00:00:00 2001 From: Hannah Ross Date: Tue, 29 Sep 2026 09:29:59 -0600 Subject: [PATCH 08/12] Update documentation for the IAG (UA_Mod=9) dynamic stall model --- docs/source/user/aerodyn/bibliography.bib | 22 ++++++++ .../aerodyn/examples/ad_primary_example.dat | 2 +- docs/source/user/aerodyn/input.rst | 26 +++++++++- docs/source/user/aerodyn/theory_ua.rst | 50 +++++++++++++++++++ docs/source/user/api_change.rst | 2 + 5 files changed, 100 insertions(+), 2 deletions(-) diff --git a/docs/source/user/aerodyn/bibliography.bib b/docs/source/user/aerodyn/bibliography.bib index 56ea901810..4360599d4e 100644 --- a/docs/source/user/aerodyn/bibliography.bib +++ b/docs/source/user/aerodyn/bibliography.bib @@ -107,6 +107,28 @@ @techreport{ad-Murray:2011 institution={49th AIAA Aerospace Sciences Meeting, Orlando, Florida} } +@article{ad-Bangga:2020, + author = {Galih Bangga and Thorsten Lutz and Matthias Arnold}, + title = {An improved second-order dynamic stall model for wind turbine airfoils}, + year = {2020}, + journal = {Wind Energy Science}, + volume = {5}, + number = {3}, + pages = {1037--1058}, + doi = {10.5194/wes-5-1037-2020} +} + +@article{ad-Bangga:2023, + author = {Galih Bangga and Jason Parkinson and William Collier}, + title = {Development and Validation of the IAG Dynamic Stall Model in State-Space Representation for Wind Turbine Airfoils}, + year = {2023}, + journal = {Energies}, + volume = {16}, + number = {10}, + pages = {3994}, + doi = {10.3390/en16103994} +} + @article{ad-hammam2022, author = {Mohamed M. Hammam and David H. Wood}, diff --git a/docs/source/user/aerodyn/examples/ad_primary_example.dat b/docs/source/user/aerodyn/examples/ad_primary_example.dat index 76f7cf97a9..50086aade5 100644 --- a/docs/source/user/aerodyn/examples/ad_primary_example.dat +++ b/docs/source/user/aerodyn/examples/ad_primary_example.dat @@ -48,7 +48,7 @@ False SectAvg - Use sector averaging (flag) "unused" OLAFInputFileName - Input file for OLAF [used only when Wake_Mod=3] ====== Unsteady Airfoil Aerodynamics Options ==================================================== True AoA34 - Sample the angle of attack (AoA) at the 3/4 chord or the AC point {default=True} [always used] - 3 UA_Mod - Unsteady Aero Model Switch (switch) {0=Quasi-steady (no UA), 2=B-L Gonzalez, 3=B-L Minnema/Pierce, 4=B-L HGM 4-states, 5=B-L HGM+vortex 5 states, 6=Oye, 7=Boeing-Vertol} + 3 UA_Mod - Unsteady Aero Model Switch (switch) {0=Quasi-steady (no UA), 2=B-L Gonzalez, 3=B-L Minnema/Pierce, 4=B-L HGM 4-states, 5=B-L HGM+vortex 5 states, 6=Oye, 7=Boeing-Vertol, 9=IAG} True FLookup - Flag to indicate whether a lookup for f' will be calculated (TRUE) or whether best-fit exponential equations will be used (FALSE); if FALSE S1-S4 must be provided in airfoil input files (flag) [used only when UA_Mod=2 or UA_Mod=3] 3 IntegrationMethod - Switch to indicate which integration method UA uses (1=RK4, 2=AB4, 3=ABM4, 4=BDF2) 0 UAStartRad - Starting radius for dynamic stall (fraction of rotor radius [0.0,1.0]) [used only when UA_Mod>0; if line is missing UAStartRad=0] diff --git a/docs/source/user/aerodyn/input.rst b/docs/source/user/aerodyn/input.rst index 892e6c86d7..d1fd9317c3 100644 --- a/docs/source/user/aerodyn/input.rst +++ b/docs/source/user/aerodyn/input.rst @@ -358,8 +358,13 @@ Most ``UA_Mod`` will require `AoA34` to be set to true. But when using quasi-ste - ``5``: 5-states continuous-time B-L model similar to HGM with an additional state for vortex generation - ``6``: 1-state continuous-time developed by Oye - ``7``: discrete-time Boeing-Vertol (BV) model +- ``9``: 5-states continuous-time IAG model (first-order, state-space), with a vortex state -Linearization is supported with ``UA_Mod=4,5,6`` (which use continuous-time states) but not with the other models. The different models are described in :numref:`AD_UA`. +Linearization is supported with ``UA_Mod=4,5,6,9`` (which use continuous-time states) but not with the other models. The different models are described in :numref:`AD_UA`. + +.. note:: + For ``UA_Mod=9``, only the first four states are linearized; the fifth (vortex) state + is excluded, exactly as for ``UA_Mod=5``. .. note:: Link to old inputs: If `UA_Mod>0`, then this is equivalent to the old `AFAeroMod=2`. @@ -841,6 +846,25 @@ or calculating it based on the polar coefficient data in the airfoil table: - ``C_lalpha`` is the slope of the 2D normal lift coefficient curve in the linear region; Used for ``UA_Mod=4,6``. +- ``Ka`` is the impulsive (non-circulatory) normal-force gain of the IAG + model; used only when ``UA_Mod=9``. If the keyword ``DEFAULT`` is entered + in place of a numerical value, ``Ka`` is set to 0.75. + +- ``Kv`` is the amplitude of the vortex center-of-pressure travel in the + IAG model; used only when ``UA_Mod=9``. If the keyword ``DEFAULT`` is + entered in place of a numerical value, ``Kv`` is set to 0.2. + +- ``dCNdA`` is the slope of the static :math:`C_n` curve used by the IAG + model; used only when ``UA_Mod=9``. If the keyword ``DEFAULT`` is entered + in place of a numerical value, ``dCNdA`` is obtained from a linear fit to + the attached-flow region of the supplied polar. Note that ``UA_Mod=9`` + uses ``dCNdA`` rather than ``C_nalpha`` or ``C_lalpha``, and that it + builds its **own** separation function by tabulating + :math:`dCNdA\,\sin(\alpha-\alpha_0)`, rather than using the + piecewise-linear fully-attached curve shared by the HGM/HGMV models. The + ``f_st`` values written for ``UA_Mod=9`` are therefore **not** on a common + basis with those written for ``UA_Mod=5``. + - ``T_f0`` is the initial value of the time constant associated with *Df* in the expressions of *Df* and *f’*; if the keyword ``DEFAULT`` is entered in place of a numerical value, ``T_f0`` is set to 3.0; diff --git a/docs/source/user/aerodyn/theory_ua.rst b/docs/source/user/aerodyn/theory_ua.rst index 2a67b8fdd5..0928896f96 100644 --- a/docs/source/user/aerodyn/theory_ua.rst +++ b/docs/source/user/aerodyn/theory_ua.rst @@ -494,6 +494,56 @@ The moment coefficient is calculated based on values at the aerodynamic center a where :math:`\alpha_{50}` is computed the same way as :math:`\alpha_{34}` (using the velocity at the aerodynamic center and the rotational rate of the airfoil) but using the distance from the aerodynamic center to the mid-chord (see :numref:`ua_notations`). +IAG model (UAMod=9) +~~~~~~~~~~~~~~~~~~~ + +The IAG model :cite:`ad-Bangga:2020,ad-Bangga:2023` is a five-state, continuous-time +Beddoes-Leishman variant. Like the HGMV model (``UA_Mod=5``) it carries two downwash memory +states :math:`x_1,x_2`, a lagged attached-flow state :math:`x_3`, a separation state +:math:`x_4`, and a vortex state :math:`x_5`. Linearization is supported, but only +:math:`x_1`-:math:`x_4` are linearized; the vortex state is excluded. + +The model differs from HGM/HGMV in three ways that matter when comparing output: + +**1. It uses a sinusoidal attached-flow curve, not a piecewise-linear one.** The circulatory +normal force is + +.. math:: + C_N^C = \frac{\mathrm{d}C_N}{\mathrm{d}\alpha}\,\sin(\alpha_E - \alpha_0) + +so the model is driven by the airfoil input ``dCNdA`` rather than by ``C_nalpha`` or +``C_lalpha``, neither of which it reads. ``dCNdA`` is obtained by a linear fit to the +attached-flow region of the supplied polar unless the user overrides it. + +**2. It builds its own separation function.** The ``f_st`` column used by ``UA_Mod=9`` is +tabulated by inverting the squared-Kirchhoff relation against the sinusoidal curve above, + +.. math:: + C_N^f = \frac{\mathrm{d}C_N}{\mathrm{d}\alpha} + \left(\frac{1+\sqrt{x_4}}{2}\right)^{2}\sin(\alpha_F-\alpha_0) + C_N^I + +rather than against the piecewise-linear ``FullyAttached`` curve shared by the HGM/HGMV +models. **The** ``f_st`` **values written for** ``UA_Mod=9`` **are therefore not on a common +basis with those written for** ``UA_Mod=5``, and the two should not be compared directly. + +**3. The effective and separation angles do not coincide.** Because the attached-flow +relation is sinusoidal, inverting it for :math:`\alpha_F` gives +:math:`\alpha_F = \alpha_0 + \sin(\alpha_E-\alpha_0)` under the model's linearized +inversion, so :math:`\alpha_F \neq \alpha_E` in general. For HGM the two collapse to the +same value. + +The impulsive (non-circulatory) normal force is scaled by the airfoil input ``Ka``, and the +vortex center-of-pressure travel by ``Kv``. Vortex shedding is triggered on the calculated +``CnMax``/``CnMin`` thresholds rather than on ``Cn1``/``Cn2``, which the model does not use. + +Beyond roughly 30 degrees of incidence the separated-flow construction loses validity, so +the dynamic :math:`C_d` and :math:`C_m` are faded linearly back to their static values +between 30 and 45 degrees. This is applied before the shared UA cutout blend, so the two +compose. :math:`C_l`, :math:`C_n` and :math:`C_c` are deliberately not faded separately, +since :math:`C_l` is reconstructed from :math:`C_n` and :math:`C_c` and blending it +independently would make the three mutually inconsistent. + + diff --git a/docs/source/user/api_change.rst b/docs/source/user/api_change.rst index f877be8949..fc653d0f6f 100644 --- a/docs/source/user/api_change.rst +++ b/docs/source/user/api_change.rst @@ -16,6 +16,8 @@ Under-relaxation is introduced for the tight-coupling iterative solver to improv 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). +The IAG dynamic stall model was added to AeroDyn as a new ``UA_Mod=9`` option. This change is **backwards compatible**: no lines are added to or removed from the AeroDyn primary input file (only the comment text of the existing ``UA_Mod`` line changes), and the three new airfoil-file inputs it introduces (``Ka``, ``Kv``, ``dCNdA``) are all optional and accept ``DEFAULT``, so existing airfoil files continue to work unchanged. The new airfoil inputs are placed after ``x_cp_bar`` and before ``UACutout``, and are read only when ``UA_Mod=9``. + ============================================= ======== ==================== ========================================================================================================================================================================================================================================================================================================== Added in OpenFAST `5.1.0` -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- From d027f6eae6ba8a0b8b9259b21ee9c9fe676248d9 Mon Sep 17 00:00:00 2001 From: Hannah Ross Date: Tue, 29 Sep 2026 11:31:27 -0600 Subject: [PATCH 09/12] Share the AFI summary-table channel count between AFI_WrHeader/AFI_WrData and udpate documentation --- docs/source/user/aerodyn/input.rst | 12 ++++++++ docs/source/user/aerodyn/theory_ua.rst | 4 +++ modules/aerodyn/src/AirfoilInfo.f90 | 42 +++++++++++++++++++++++--- 3 files changed, 54 insertions(+), 4 deletions(-) diff --git a/docs/source/user/aerodyn/input.rst b/docs/source/user/aerodyn/input.rst index d1fd9317c3..fd94a320aa 100644 --- a/docs/source/user/aerodyn/input.rst +++ b/docs/source/user/aerodyn/input.rst @@ -865,6 +865,18 @@ or calculating it based on the polar coefficient data in the airfoil table: ``f_st`` values written for ``UA_Mod=9`` are therefore **not** on a common basis with those written for ``UA_Mod=5``. +- ``CnMax`` and ``CnMin`` are not airfoil-file inputs. They are the maximum + and minimum of the static :math:`C_n` polar, calculated at initialization, + and they set the positive and negative vortex-shedding thresholds of the + IAG model. They are reported in the unsteady-aero summary table written + when ``UA_Mod=9`` and ``SumPrint = TRUE``. A reported value of + ``999.00000`` for ``CnMax`` (and ``-999.00000`` for ``CnMin``) is a + sentinel, not a physical coefficient: it means the vortex logic has been + disabled for that table, either because the model is not IAG or because + the polar has no identifiable stall peak, as for a cylinder-like table at + the blade root. Because :math:`|C_n|` can never reach 999, the shedding + trigger simply never fires and the model runs without vortex lift. + - ``T_f0`` is the initial value of the time constant associated with *Df* in the expressions of *Df* and *f’*; if the keyword ``DEFAULT`` is entered in place of a numerical value, ``T_f0`` is set to 3.0; diff --git a/docs/source/user/aerodyn/theory_ua.rst b/docs/source/user/aerodyn/theory_ua.rst index 0928896f96..731adea891 100644 --- a/docs/source/user/aerodyn/theory_ua.rst +++ b/docs/source/user/aerodyn/theory_ua.rst @@ -535,6 +535,10 @@ same value. The impulsive (non-circulatory) normal force is scaled by the airfoil input ``Ka``, and the vortex center-of-pressure travel by ``Kv``. Vortex shedding is triggered on the calculated ``CnMax``/``CnMin`` thresholds rather than on ``Cn1``/``Cn2``, which the model does not use. +``CnMax`` and ``CnMin`` are the extrema of the static :math:`C_n` polar and are reported in +the unsteady-aero summary table; values of :math:`\pm 999` there are a sentinel indicating +that no stall peak could be identified for that table and that vortex shedding is therefore +disabled for it (see :numref:`airfoil_data_input_file`). Beyond roughly 30 degrees of incidence the separated-flow construction loses validity, so the dynamic :math:`C_d` and :math:`C_m` are faded linearly back to their static values diff --git a/modules/aerodyn/src/AirfoilInfo.f90 b/modules/aerodyn/src/AirfoilInfo.f90 index 9177ad97ac..c5df26bdfc 100644 --- a/modules/aerodyn/src/AirfoilInfo.f90 +++ b/modules/aerodyn/src/AirfoilInfo.f90 @@ -43,6 +43,24 @@ MODULE AirfoilInfo integer, parameter :: MaxNumAFCoeffs = 7 !cl,cd,cm,cpMin, UA:f_st, FullySeparate, FullyAttached + !> Column count and field width of the UA summary table written by AFI_WrHeader / AFI_WrData. + !! + !! These MUST be shared: the two routines write the header line and the data lines of the same + !! file, and nothing else ties them together. They were previously declared as separate local + !! parameters in each routine, so adding a channel to one and not the other silently misaligned + !! every column of the table -- a failure with no error message and no obvious symptom. + !! + !! AFI_NumSumChans counts ALL columns, including the two leading integers (AirfoilNumber, + !! TableNumber). When adding a channel you must therefore touch three places: + !! 1. bump AFI_NumSumChans here; + !! 2. add the ChanName/ChanUnit pair in AFI_WrHeader; + !! 3. add the matching UA_BL field to the write list in AFI_WrData. + !! Step 1 vs 2 is checked at run time by the assertion at the end of the AFI_WrHeader + !! assignment block. Step 3 cannot be checked automatically -- a Fortran output list has no + !! introspectable length -- so it remains the one manual step. + integer, parameter :: AFI_NumSumChans = 51 + integer, parameter :: AFI_SumChanWidth = 17 + !> Sentinel written into UA_BL%CnMax / CnMin when the IAG (UAMod=9) vortex logic must stay !! dormant: non-IAG models, and degenerate (cylinder-like) polars with no stall peak. !! @@ -2267,12 +2285,16 @@ subroutine AFI_WrHeader(delim, FileName, unOutFile, ErrStat, ErrMsg) character(ErrMsgLen) :: ErrMsg2 character(*), parameter :: RoutineName = 'AFI_WrHeader' - integer, parameter :: MaxLen = 17 - integer, parameter :: NumChans = 51 + integer, parameter :: MaxLen = AFI_SumChanWidth + integer, parameter :: NumChans = AFI_NumSumChans character(MaxLen) :: ChanName( NumChans) character(MaxLen) :: ChanUnit( NumChans) + ErrStat = ErrID_None + ErrMsg = '' + unOutFile = -1 ! defined on every path, including the early return below + i=1 ChanName(i) = 'AirfoilNumber'; ChanUnit(i) = '(-)'; i = i+1; ChanName(i) = 'TableNumber'; ChanUnit(i) = '(-)'; i = i+1; @@ -2326,6 +2348,18 @@ subroutine AFI_WrHeader(delim, FileName, unOutFile, ErrStat, ErrMsg) ChanName(i) = 'CnMax'; ChanUnit(i) = '(-)'; i = i+1; ChanName(i) = 'CnMin'; ChanUnit(i) = '(-)'; i = i+1; + ! Catch a channel added here without bumping AFI_NumSumChans (or vice versa). Without this, + ! too few assignments leave trailing columns filled with uninitialized garbage, and too many + ! overrun the array -- neither of which is reported unless bounds checking happens to be on. + ! This cannot check AFI_WrData's write list; see the AFI_NumSumChans comment. + if ( i-1 /= NumChans ) then + call SetErrStat( ErrID_Fatal, 'Programming error: AFI_WrHeader assigned '//trim(num2lstr(i-1))// & + ' channels but AFI_NumSumChans is '//trim(num2lstr(NumChans))// & + '. Update AFI_NumSumChans and the AFI_WrData write list to match.', & + ErrStat, ErrMsg, RoutineName ) + return + end if + !$OMP critical(fileopen_critical) CALL GetNewUnit( unOutFile, ErrStat, ErrMsg ) if (ErrStat < AbortErrLev) then @@ -2370,8 +2404,8 @@ subroutine AFI_WrData(k, unOutFile, delim, AFInfo) integer(IntKi) :: i - integer, parameter :: MaxLen = 17 - integer, parameter :: NumChans = 51 + integer, parameter :: MaxLen = AFI_SumChanWidth + integer, parameter :: NumChans = AFI_NumSumChans real(ReKi) :: TmpValues(NumChans) character(3) :: MaxLenStr character(80) :: Fmt From f56e9e0b18adcc1a62deb533b2f54be26479d438 Mon Sep 17 00:00:00 2001 From: Hannah Ross Date: Wed, 30 Sep 2026 09:50:56 -0600 Subject: [PATCH 10/12] Add UA_SetIAGBlendBounds verification hook and document that the IAG Cc output channel is not comparable with UA_Mod=4/5 --- docs/source/user/aerodyn/theory_ua.rst | 17 +++- modules/aerodyn/src/UnsteadyAero.f90 | 123 ++++++++++++++++++++++++- 2 files changed, 134 insertions(+), 6 deletions(-) diff --git a/docs/source/user/aerodyn/theory_ua.rst b/docs/source/user/aerodyn/theory_ua.rst index 731adea891..2f27beca7d 100644 --- a/docs/source/user/aerodyn/theory_ua.rst +++ b/docs/source/user/aerodyn/theory_ua.rst @@ -503,7 +503,7 @@ states :math:`x_1,x_2`, a lagged attached-flow state :math:`x_3`, a separation s :math:`x_4`, and a vortex state :math:`x_5`. Linearization is supported, but only :math:`x_1`-:math:`x_4` are linearized; the vortex state is excluded. -The model differs from HGM/HGMV in three ways that matter when comparing output: +The model differs from HGM/HGMV in four ways that matter when comparing output: **1. It uses a sinusoidal attached-flow curve, not a piecewise-linear one.** The circulatory normal force is @@ -532,6 +532,21 @@ relation is sinusoidal, inverting it for :math:`\alpha_F` gives inversion, so :math:`\alpha_F \neq \alpha_E` in general. For HGM the two collapse to the same value. +**4. The** ``Cc`` **output channel is a different quantity than for the other models.** +When unsteady-aero outputs are enabled, ``UA_Mod=9`` writes the *viscous* chordwise force +evaluated at :math:`\alpha_F` with :math:`C_{d0}` removed, + +.. math:: + C_c = C_l^{st}(\alpha_F)\sin\alpha_F - \left[C_d^{st}(\alpha_F) - C_{d0}\right]\cos\alpha_F + +which is the :math:`C_T^D` of the IAG formulation. The HGM and HGMV models instead write +:math:`C_c = C_l\sin\alpha - C_d\cos\alpha` evaluated at the instantaneous :math:`\alpha` +and **without** subtracting :math:`C_{d0}`. The two are not the same quantity and **the** +``Cc`` **column should not be compared between** ``UA_Mod=9`` **and the other models.** +This affects the reported channel only: :math:`C_l`, :math:`C_d` and :math:`C_m` are the +quantities passed to the rest of AeroDyn, and for the IAG model those are reconstructed +from :math:`C_N^D` and :math:`C_T^D` before they are returned, so loads are unaffected. + The impulsive (non-circulatory) normal force is scaled by the airfoil input ``Ka``, and the vortex center-of-pressure travel by ``Kv``. Vortex shedding is triggered on the calculated ``CnMax``/``CnMin`` thresholds rather than on ``Cn1``/``Cn2``, which the model does not use. diff --git a/modules/aerodyn/src/UnsteadyAero.f90 b/modules/aerodyn/src/UnsteadyAero.f90 index 317ed00d5d..3b8ca50c24 100644 --- a/modules/aerodyn/src/UnsteadyAero.f90 +++ b/modules/aerodyn/src/UnsteadyAero.f90 @@ -59,6 +59,9 @@ module UnsteadyAero public :: UA_ReInit public :: UA_InitStates_AllNodes ! used for AD linearization initialization + public :: UA_SetIAGBlendBounds ! verification hook; see the declarations below + public :: UA_GetIAGBlendBounds + real(ReKi), parameter :: Gonzalez_factor = 0.2_ReKi ! this factor, proposed by Gonzalez (for "all" models) is used to modify Cc to account for negative values seen at f=0 (see Eqn 1.40) real(ReKi), parameter, public :: UA_u_min = 0.01_ReKi ! m/s; used to provide a minimum value so UA equations don't blow up (this should be much lower than range where UA is turned off) real(ReKi), parameter :: K1pos=1.0_ReKi, K1neg=0.5_ReKi ! K1 coefficients for BV model @@ -66,11 +69,26 @@ module UnsteadyAero ! IAG model blend limits, in degrees of |alpha|. The dynamic Cd/Cm are faded to static ! between BlendLo and BlendHi, and the vortex state x5 is faded out between x5FadeLo and ! x5FadeHi. Both come from the IAG theory text rather than the paper's equations, so they - ! are named constants here rather than user inputs -- see the notes on Task 5. - real(ReKi), parameter :: IAG_BlendLo = 30.0_ReKi ! deg; start of the dynamic->static blend - real(ReKi), parameter :: IAG_BlendHi = 45.0_ReKi ! deg; end of the dynamic->static blend - real(ReKi), parameter :: IAG_x5FadeLo = 45.0_ReKi ! deg; start of the x5 fade-out - real(ReKi), parameter :: IAG_x5FadeHi = 75.0_ReKi ! deg; end of the x5 fade-out + ! are not exposed in the user input files -- see the notes on Task 5. + ! + ! They are module variables rather than PARAMETERs so that verification drivers can widen + ! the bands and reach model branches that the default schedule makes unreachable. In + ! particular the backwinded (|alpha| > 90 deg) CM_C sign flip in the IAG CalcOutput branch + ! is multiplied by a blend weight that is identically zero for |alpha| >= IAG_BlendHi, so + ! with the defaults that branch can never influence y%Cm. See UA_SetIAGBlendBounds. + ! + ! These are PROTECTED: only UA_SetIAGBlendBounds may modify them, and it is intended for + ! test/verification use. They are global (not per-instance), so a driver that changes them + ! changes them for every UA instance in the process. + real(ReKi), parameter :: IAG_BlendLo_Def = 30.0_ReKi ! deg; default start of the dynamic->static blend + real(ReKi), parameter :: IAG_BlendHi_Def = 45.0_ReKi ! deg; default end of the dynamic->static blend + real(ReKi), parameter :: IAG_x5FadeLo_Def = 45.0_ReKi ! deg; default start of the x5 fade-out + real(ReKi), parameter :: IAG_x5FadeHi_Def = 75.0_ReKi ! deg; default end of the x5 fade-out + + real(ReKi), protected :: IAG_BlendLo = IAG_BlendLo_Def + real(ReKi), protected :: IAG_BlendHi = IAG_BlendHi_Def + real(ReKi), protected :: IAG_x5FadeLo = IAG_x5FadeLo_Def + real(ReKi), protected :: IAG_x5FadeHi = IAG_x5FadeHi_Def contains @@ -3266,6 +3284,95 @@ pure real(ReKi) function Get_Cc( AFI_interp, alpha ) Get_Cc = AFI_interp%Cl*sin(alpha) - (AFI_interp%Cd - AFI_interp%Cd0)*cos(alpha) end function Get_Cc !--------------------------------------------------------------------------------- +!> Override the IAG blend/fade-out bounds, in degrees of |alpha|. +!! +!! *** THIS IS A VERIFICATION HOOK, NOT A USER INPUT. *** +!! +!! The defaults (30/45 deg for the Cd/Cm dynamic->static blend, 45/75 deg for the x5 +!! vortex fade) reproduce the IAG theory text and are what production runs must use. +!! They are exposed here only so that unit-test drivers can widen the bands and reach +!! branches of the IAG output path that the default schedule renders unreachable -- +!! most notably the backwinded (|alpha| > 90 deg) CM_C sign flip, whose contribution to +!! y%Cm is multiplied by a weight that is identically zero for |alpha| >= IAG_BlendHi. +!! +!! Passing a bound as absent leaves that bound at its current value. Call with no +!! optional arguments plus Reset=.true. to restore all four defaults. +!! +!! NOTE: these are module-global, not per-instance. Changing them affects every UA +!! instance in the process, and they are not saved/restored across UA_Init. +subroutine UA_SetIAGBlendBounds( BlendLo, BlendHi, x5FadeLo, x5FadeHi, Reset, ErrStat, ErrMsg ) + real(ReKi), optional, intent(in ) :: BlendLo !< deg; start of the dynamic->static Cd/Cm blend + real(ReKi), optional, intent(in ) :: BlendHi !< deg; end of the dynamic->static Cd/Cm blend + real(ReKi), optional, intent(in ) :: x5FadeLo !< deg; start of the x5 vortex fade-out + real(ReKi), optional, intent(in ) :: x5FadeHi !< deg; end of the x5 vortex fade-out + logical, optional, intent(in ) :: Reset !< if .true., restore defaults first + integer(IntKi), intent( out) :: ErrStat !< error status + character(*), intent( out) :: ErrMsg !< error message + + character(*), parameter :: RoutineName = 'UA_SetIAGBlendBounds' + real(ReKi) :: bLo, bHi, xLo, xHi + + ErrStat = ErrID_None + ErrMsg = "" + + if (present(Reset)) then + if (Reset) then + IAG_BlendLo = IAG_BlendLo_Def + IAG_BlendHi = IAG_BlendHi_Def + IAG_x5FadeLo = IAG_x5FadeLo_Def + IAG_x5FadeHi = IAG_x5FadeHi_Def + end if + end if + + ! stage into locals so a rejected set leaves the module state untouched + bLo = IAG_BlendLo; bHi = IAG_BlendHi + xLo = IAG_x5FadeLo; xHi = IAG_x5FadeHi + + if (present(BlendLo)) bLo = BlendLo + if (present(BlendHi)) bHi = BlendHi + if (present(x5FadeLo)) xLo = x5FadeLo + if (present(x5FadeHi)) xHi = x5FadeHi + + ! The blends are written as (|alpha|-Lo)/(Hi-Lo), so Hi must be strictly greater than Lo + ! or the weight is a divide-by-zero. |alpha| is wrapped to [-pi,pi], so 180 deg is the + ! largest meaningful bound and a bound at/above it simply disables the blend. + if (bHi <= bLo) then + call SetErrStat(ErrID_Fatal, 'IAG BlendHi ('//trim(num2lstr(bHi))//' deg) must be greater '// & + 'than BlendLo ('//trim(num2lstr(bLo))//' deg).', ErrStat, ErrMsg, RoutineName) + end if + if (xHi <= xLo) then + call SetErrStat(ErrID_Fatal, 'IAG x5FadeHi ('//trim(num2lstr(xHi))//' deg) must be greater '// & + 'than x5FadeLo ('//trim(num2lstr(xLo))//' deg).', ErrStat, ErrMsg, RoutineName) + end if + if (bLo < 0.0_ReKi .or. xLo < 0.0_ReKi) then + call SetErrStat(ErrID_Fatal, 'IAG blend lower bounds must be non-negative (degrees of |alpha|).', & + ErrStat, ErrMsg, RoutineName) + end if + if (ErrStat >= AbortErrLev) return + + IAG_BlendLo = bLo + IAG_BlendHi = bHi + IAG_x5FadeLo = xLo + IAG_x5FadeHi = xHi + +end subroutine UA_SetIAGBlendBounds +!--------------------------------------------------------------------------------- +!> Report the IAG blend/fade-out bounds currently in effect, in degrees of |alpha|. +!! Lets a driver record what it actually ran with, and lets a test assert that the +!! defaults are unchanged. +subroutine UA_GetIAGBlendBounds( BlendLo, BlendHi, x5FadeLo, x5FadeHi ) + real(ReKi), optional, intent( out) :: BlendLo + real(ReKi), optional, intent( out) :: BlendHi + real(ReKi), optional, intent( out) :: x5FadeLo + real(ReKi), optional, intent( out) :: x5FadeHi + + if (present(BlendLo)) BlendLo = IAG_BlendLo + if (present(BlendHi)) BlendHi = IAG_BlendHi + if (present(x5FadeLo)) x5FadeLo = IAG_x5FadeLo + if (present(x5FadeHi)) x5FadeHi = IAG_x5FadeHi + +end subroutine UA_GetIAGBlendBounds +!--------------------------------------------------------------------------------- !> IAG model: linearly blend the dynamic Cd and Cm back to their static values between !! IAG_BlendLo and IAG_BlendHi degrees of incidence. !! Past ~30 deg the model's separated-flow construction loses validity, so the dynamic @@ -4143,6 +4250,12 @@ subroutine UA_CalcOutput( i, j, t, u_in, p, x, xd, OtherState, AFInfo, y, misc, CT_D = Get_Cc( AFI_interpF, alphaF ) ! Eq. 51, viscous Cc at alphaF y%Cn = CN_D + ! NOTE y%Cc is NOT the same quantity here as on the HGM/HGMV paths, which set + ! y%Cc = y%Cl*SinAlpha - y%Cd*CosAlpha at u%alpha and do NOT subtract Cd0. + ! This is CT_D: the viscous chordwise force at alphaF, Cd0 removed (Eq. 51). + ! Nothing downstream reads y%Cc -- BEMT, AeroDyn and FVW take only Cl/Cd/Cm -- so + ! this affects the UA_OUTS 'Cc' channel and nothing else. It does mean the Cc + ! column is not comparable between UA_Mod=9 and UA_Mod=4/5; see theory_ua.rst. y%Cc = CT_D y%Cl = CN_D*CosAlpha - CT_D*SinAlpha ! Eq. 52 From 693614f905255729438739171110f71770f0db91 Mon Sep 17 00:00:00 2001 From: Hannah Ross Date: Wed, 30 Sep 2026 13:11:16 -0600 Subject: [PATCH 11/12] Add mandatory ErrStat/ErrMsg to Get_alphaF so an unhandled UAMod aborts instead of silently returning alphaF=0 --- modules/aerodyn/src/UnsteadyAero.f90 | 38 +++++++++++++++++++++++----- 1 file changed, 31 insertions(+), 7 deletions(-) diff --git a/modules/aerodyn/src/UnsteadyAero.f90 b/modules/aerodyn/src/UnsteadyAero.f90 index 3b8ca50c24..418259bc22 100644 --- a/modules/aerodyn/src/UnsteadyAero.f90 +++ b/modules/aerodyn/src/UnsteadyAero.f90 @@ -2885,7 +2885,9 @@ SUBROUTINE HGM_Steady( i, j, u, p, x, AFInfo, ErrStat, ErrMsg ) ! small (~0.4 deg at 20 deg from alpha0) but f_st is steep near stall, so evaluating the ! separation function at alphaE would produce a visible start-up transient. ! NOTE: x%x(3) MUST already be assigned -- Get_alphaF's UA_IAG branch reads it from x. - alphaF = Get_alphaF(p, u, x, BL_p, alpha_34, alphaE) ! Eq. 41 + alphaF = Get_alphaF(p, u, x, BL_p, alpha_34, alphaE, ErrStat2, ErrMsg2) ! Eq. 41 + call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) + if (ErrStat >= AbortErrLev) return call AFI_ComputeAirfoilCoefs( alphaF, u%Re, u%UserProp, AFInfo, AFI_interp_F, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) if (ErrStat >= AbortErrLev) return @@ -2982,7 +2984,9 @@ subroutine UA_CalcContStateDeriv( i, j, t, u_in, p, x, OtherState, AFInfo, m, dx ! calculate fs_aF (stored in AFI_interp%f_st): ! find alphaF where FullyAttached(alphaF) = x(3) - alphaF = Get_alphaF(p, u, x, BL_p, alpha_34, alphaE) + alphaF = Get_alphaF(p, u, x, BL_p, alpha_34, alphaE, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) + if (ErrStat >= AbortErrLev) return call AFI_ComputeAirfoilCoefs( alphaF, u%Re, u%UserProp, AFInfo, AFI_AlphaF, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) @@ -3165,7 +3169,15 @@ SUBROUTINE Get_HGM_constants(i, j, p, u, x, BL_p, Tu, alpha_34, alphaE) END SUBROUTINE Get_HGM_constants !--------------------------------------------------------------------------------- -FUNCTION Get_alphaF(p, u, x, BL_p, alpha_34, alphaE_in) RESULT(alphaF) +!> Invert the fully-attached normal-force (or lift) curve to find the angle alphaF at which the +!! attached-flow coefficient equals the lagged state x3. Dispatches on p%UAMod. +!! +!! ErrStat/ErrMsg are mandatory: the final ELSE of the dispatcher is reachable only if a new +!! UAMod is added without a branch here, and the failure mode is silent rather than loud +!! (alphaF = 0 places the separation lookup at the wrong angle, typically giving f_st ~ 1, i.e. +!! stall quietly switched off, and the run completes with plausible-looking but wrong loads). +!! Callers must therefore check ErrStat and abort. +FUNCTION Get_alphaF(p, u, x, BL_p, alpha_34, alphaE_in, ErrStat, ErrMsg) RESULT(alphaF) TYPE(UA_InputType), INTENT(IN ) :: u ! Inputs at t TYPE(UA_ParameterType), INTENT(IN ) :: p ! Parameters TYPE(UA_ElementContinuousStateType), INTENT(IN ) :: x ! Continuous states at t @@ -3173,13 +3185,18 @@ FUNCTION Get_alphaF(p, u, x, BL_p, alpha_34, alphaE_in) RESULT(alphaF) REAL(ReKi), INTENT(IN ) :: alpha_34 REAL(ReKi), INTENT(IN ) :: alphaE_in + INTEGER(IntKi), INTENT( OUT) :: ErrStat ! Error status of the operation + CHARACTER(*), INTENT( OUT) :: ErrMsg ! Error message if ErrStat /= ErrID_None REAL(ReKi) :: alphaF ! function result REAL(ReKi) :: alphaE ! value that can be changed (+/- 2pi) REAL(ReKi) :: alpha_(2), c_(2) REAL(ReKi) :: alphaN_(4), cN_(4) integer(IntKi) :: Indx + character(*), parameter :: RoutineName = 'Get_alphaF' + ErrStat = ErrID_None + ErrMsg = "" alphaE = alphaE_in @@ -3229,8 +3246,11 @@ FUNCTION Get_alphaF(p, u, x, BL_p, alpha_34, alphaE_in) RESULT(alphaF) end if else - !PROGRAMMING ERROR IF WE GET TO THIS PART OF THE IF STATEMENT! - alphaF = 0 + ! PROGRAMMING ERROR: a UAMod reached this routine without an inversion branch. Returning + ! 0 silently would disable stall rather than stop, so raise a fatal error instead. + alphaF = 0.0_ReKi + call SetErrStat(ErrID_Fatal, 'Programming error: Get_alphaF has no branch for UAMod='// & + trim(num2lstr(p%UAMod))//'.', ErrStat, ErrMsg, RoutineName) end if @@ -4151,7 +4171,9 @@ subroutine UA_CalcOutput( i, j, t, u_in, p, x, xd, OtherState, AFInfo, y, misc, y%Cl = y%Cn * CosAlpha + y%Cc * SinAlpha; - alphaF = Get_alphaF(p, u, x_in, BL_p, alpha_34, alphaE) + alphaF = Get_alphaF(p, u, x_in, BL_p, alpha_34, alphaE, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + if (ErrStat >= AbortErrLev) return call AFI_ComputeAirfoilCoefs( alphaF, u%Re, u%UserProp, AFInfo, AFI_interpF, ErrStat2, ErrMsg2 ) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) @@ -4227,7 +4249,9 @@ subroutine UA_CalcOutput( i, j, t, u_in, p, x, xd, OtherState, AFInfo, y, misc, alpha_w = u%alpha call MPi2Pi(alpha_w) ! for the CM_C and blend tests - alphaF = Get_alphaF(p, u, x_in, BL_p, alpha_34, alphaE) ! Eq. 41 + alphaF = Get_alphaF(p, u, x_in, BL_p, alpha_34, alphaE, ErrStat2, ErrMsg2) ! Eq. 41 + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + if (ErrStat >= AbortErrLev) return call AFI_ComputeAirfoilCoefs( alphaF, u%Re, u%UserProp, AFInfo, AFI_interpF, ErrStat2, ErrMsg2 ) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return From 896bb193fd690bf5fd0d3b2ab61a3a16b0f2465d Mon Sep 17 00:00:00 2001 From: Hannah Ross Date: Wed, 30 Sep 2026 16:46:54 -0600 Subject: [PATCH 12/12] Update r-tests --- reg_tests/r-test | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/reg_tests/r-test b/reg_tests/r-test index 8aefcb0e74..bc1bf601e4 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit 8aefcb0e74b711bc7eebf77fab246d2f9bb7c817 +Subproject commit bc1bf601e45c9d2245ba681355fd7e0575b179bf