From cf58377c890add573dfb2f332cbb7952ba7c0b11 Mon Sep 17 00:00:00 2001 From: Bonnie Jonkman Date: Tue, 8 Mar 2022 14:24:02 -0700 Subject: [PATCH 1/3] AD/AA: remove unused variables --- modules/aerodyn/src/AeroAcoustics.f90 | 38 +++++++++------------------ 1 file changed, 13 insertions(+), 25 deletions(-) diff --git a/modules/aerodyn/src/AeroAcoustics.f90 b/modules/aerodyn/src/AeroAcoustics.f90 index dbc0c6b214..58d11b8d34 100644 --- a/modules/aerodyn/src/AeroAcoustics.f90 +++ b/modules/aerodyn/src/AeroAcoustics.f90 @@ -1101,7 +1101,7 @@ SUBROUTINE CalcAeroAcousticsOutput(u,p,m,xd,y,errStat,errMsg) ! Amiet's Inflow Noise Model is Calculated as long as InflowNoise is On CALL InflowNoise(AlphaNoise,p%BlChord(J,I),Unoise,m%ChordAngleLE(K,J,I),m%SpanAngleLE(K,J,I),& - elementspan,m%rLEtoObserve(K,J,I),xd%MeanVxVyVz(J,I),xd%TIVx(J,I),m%LE_Location(3,J,I),0.050,p,m%SPLti,errStat2,errMsg2 ) + elementspan,m%rLEtoObserve(K,J,I),xd%TIVx(J,I),p,m%SPLti,errStat2,errMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) ! If Guidati model (simplified or full version) is also on then the 'SPL correction' to Amiet's model will be added IF ( p%IInflow .EQ. 2 ) THEN @@ -1648,30 +1648,19 @@ SUBROUTINE TIPNOIS(ALPHTIP,ALPRAT2,C,U ,THETA,PHI, R,p,SPLTIP, errStat, errMsg) ENDDO END SUBROUTINE TipNois !==================================================================================================================================! -SUBROUTINE InflowNoise(AlphaNoise,Chord,U,THETA,PHI,d,RObs,MeanVNoise,TINoise,LE_Location,dissip,p,SPLti,errStat,errMsg) -! REAL(ReKi), INTENT(IN ) :: AlphaNoise ! AOA -! REAL(ReKi), INTENT(IN ) :: Chord ! Chord Length -! REAL(ReKi), INTENT(IN ) :: U ! -! REAL(ReKi), INTENT(IN ) :: d ! element span -! REAL(ReKi), INTENT(IN ) :: RObs ! distance to observer -! REAL(ReKi), INTENT(IN ) :: THETA ! -! REAL(ReKi), INTENT(IN ) :: PHI ! Spanwise directivity angle +SUBROUTINE InflowNoise(AlphaNoise,Chord,U,THETA,PHI,d,RObs,TINoise,p,SPLti,errStat,errMsg) + REAL(ReKi), INTENT(IN ) :: AlphaNoise ! AOA + REAL(ReKi), INTENT(IN ) :: Chord ! Chord Length + REAL(ReKi), INTENT(IN ) :: U ! + REAL(ReKi), INTENT(IN ) :: THETA ! + REAL(ReKi), INTENT(IN ) :: PHI ! Spanwise directivity angle + REAL(ReKi), INTENT(IN ) :: d ! element span + REAL(ReKi), INTENT(IN ) :: RObs ! distance to observer ! REAL(ReKi), INTENT(IN ) :: MeanVNoise ! -! REAL(ReKi), INTENT(IN ) :: TINoise ! + REAL(ReKi), INTENT(IN ) :: TINoise ! ! REAL(ReKi), INTENT(IN ) :: LE_Location ! - - REAL(ReKi) :: AlphaNoise ! AOA - REAL(ReKi) :: Chord ! Chord Length - REAL(ReKi) :: U ! - REAL(ReKi) :: d ! element span - REAL(ReKi) :: RObs ! distance to observer - REAL(ReKi) :: THETA ! - REAL(ReKi) :: PHI ! Spanwise directivity angle - REAL(ReKi) :: MeanVNoise ! - REAL(ReKi) :: TINoise ! - REAL(ReKi) :: LE_Location ! - - REAL(ReKi), INTENT(IN ) :: dissip ! + +! REAL(ReKi), INTENT(IN ) :: dissip ! TYPE(AA_ParameterType), INTENT(IN ) :: p ! Parameters REAL(ReKi),DIMENSION(size(p%FreqList)), INTENT( OUT) :: SPLti ! INTEGER(IntKi), INTENT( OUT) :: errStat ! Error status of the operation @@ -2481,8 +2470,7 @@ SUBROUTINE Aero_Tests() !CALL TIPNOIS(AlphaNoise,p%ALpRAT,p%BlChord(J,I),UNoise,m%ChordAngleTE(K,J,I),m%SpanAngleTE(K,J,I), & ! m%rTEtoObserve(K,J,I), p, m%SPLTIP,ErrStat2,errMsg2) !--------Inflow Turbulence Noise ------------------------------------------------! - !CALL InflowNoise(3.0d0,0.22860d0,63.920d0,90.0d0,90.0d0,0.5090d0,1.220d0, & - ! xd%MeanVrel(J,I),0.050d0,0.050d0,p,m%SPLti,ErrStat2,errMsg2 ) + !CALL InflowNoise(3.0d0,0.22860d0,63.920d0,90.0d0,90.0d0,0.5090d0,1.220d0, xd%TIVx(J,I),0.050d0,p,m%SPLti,ErrStat2,errMsg2 ) !CALL FullGuidati(3.0d0,63.920d0,0.22860d0,0.5090d0,1.220d0,90.0d0,90.0d0,xd%MeanVrel(J,I),xd%TIVrel(J,I), & ! p,p%BlAFID(J,I),m%SPLTIGui,ErrStat2 ) !CALL Simple_Guidati(UNoise,0.22860d0,0.120d0,0.020d0,p,m%SPLTIGui,ErrStat2,errMsg2 ) From 17795d7c06361d002847fa73e0acfaa749b81e32 Mon Sep 17 00:00:00 2001 From: Bonnie Jonkman Date: Wed, 9 Mar 2022 14:35:15 -0700 Subject: [PATCH 2/3] AD sync: sort BEMT/FVW outputs in AeroDyn_IO.f90 - some outputs weren't calculated for both FVW and BEMT but could be (merged into one subroutine) - also adds a couple of variables to help sync with Envision --- modules/aerodyn/src/AeroDyn.f90 | 154 +++++++-- modules/aerodyn/src/AeroDyn_IO.f90 | 473 ++++++++++++-------------- modules/aerodyn/src/BEMT_Registry.txt | 1 + modules/aerodyn/src/BEMT_Types.f90 | 54 +++ 4 files changed, 397 insertions(+), 285 deletions(-) diff --git a/modules/aerodyn/src/AeroDyn.f90 b/modules/aerodyn/src/AeroDyn.f90 index 6a0598be7f..29f3a07cfc 100644 --- a/modules/aerodyn/src/AeroDyn.f90 +++ b/modules/aerodyn/src/AeroDyn.f90 @@ -276,7 +276,7 @@ subroutine AD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut if (Failed()) return; NumBlades(iR) = InitInp%rotors(iR)%NumBlades p%rotors(iR)%NumBlades = InitInp%rotors(iR)%NumBlades - if (nRotors>1) then + if (nRotors > 1) then p%rotors(iR)%RootName = TRIM(InitInp%RootName)//'.AD.R'//trim(num2lstr(iR)) else p%rotors(iR)%RootName = TRIM(InitInp%RootName)//'.AD' @@ -1160,7 +1160,8 @@ subroutine AD_End( u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg ) TYPE(AD_MiscVarType), INTENT(INOUT) :: m !< Misc/optimization variables INTEGER(IntKi), INTENT( OUT) :: ErrStat !< Error status of the operation CHARACTER(*), INTENT( OUT) :: ErrMsg !< Error message if ErrStat /= ErrID_None - integer :: iW + + integer :: iW @@ -1513,7 +1514,7 @@ subroutine RotWriteOutputs( t, u, p, p_AD, x, xd, z, OtherState, y, m, m_AD, iRo ! get values to output to file: !------------------------------------------------------- if (p%NumOuts > 0) then - call Calc_WriteOutput( p, p_AD, u, m, m_AD, y, OtherState, xd, indx, iRot, ErrStat2, ErrMsg2 ) + call Calc_WriteOutput( p, p_AD, u, x, m, m_AD, y, OtherState, xd, indx, iRot, ErrStat2, ErrMsg2 ) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) !............................................................................................................................... @@ -1722,6 +1723,7 @@ subroutine SetInputs(p, p_AD, u, m, indx, errStat, errMsg) endif end subroutine SetInputs +!---------------------------------------------------------------------------------------------------------------------------------- !> Disturbed inflow on the blade if tower shadow or tower influence are enabled subroutine SetDisturbedInflow(p, u, m, errStat, errMsg) type(RotParameterType), intent(in ) :: p !< AD parameters @@ -1777,14 +1779,15 @@ subroutine SetInputsForBEMT(p, u, m, indx, errStat, errMsg) ! note ErrStat and ErrMsg are set in GeomWithoutSweepPitchTwist: ! Get disk average values and orientations - call DiskAvgValues(p, u, m, x_hat_disk, y_hat_disk, z_hat_disk, Azimuth) - call GeomWithoutSweepPitchTwist(p,u,m,thetaBladeNds,ErrStat,ErrMsg) + call DiskAvgValues(p, u, m, x_hat_disk, y_hat_disk, z_hat_disk, Azimuth) ! also sets m%V_diskAvg, m%V_dot_x + call GeomWithoutSweepPitchTwist(p,u,x_hat_disk,m,thetaBladeNds,ErrStat,ErrMsg) if (ErrStat >= AbortErrLev) return ! Velocity in disk normal m%BEMT_u(indx)%Un_disk = m%V_dot_x - ! "Angular velocity of rotor" rad/s + + ! "Angular velocity of rotor" rad/s m%BEMT_u(indx)%omega = dot_product( u%HubMotion%RotationVel(:,1), x_hat_disk ) ! "Angle between the vector normal to the rotor plane and the wind vector (e.g., the yaw angle in the case of no tilt)" rad @@ -1850,8 +1853,10 @@ subroutine SetInputsForBEMT(p, u, m, indx, errStat, errMsg) m%BEMT_u(indx)%UserProp = u%UserProp - ! ................ TSR ..................... + !.......................... + ! TSR + !.......................... if ( EqualRealNos( m%V_dot_x, 0.0_ReKi ) ) then m%BEMT_u(indx)%TSR = 0.0_ReKi else @@ -1924,9 +1929,58 @@ subroutine DiskAvgValues(p, u, m, x_hat_disk, y_hat_disk, z_hat_disk, Azimuth) end subroutine DiskAvgValues !---------------------------------------------------------------------------------------------------------------------------------- -subroutine GeomWithoutSweepPitchTwist(p,u,m,thetaBladeNds,ErrStat,ErrMsg) +subroutine Calculate_MeshOrientation_Rel2Hub(Mesh1, HubMotion, x_hat_disk, orientationAnnulus, elemPosRelToHub_save, elemPosRotorProj_save) + TYPE(MeshType), intent(in) :: Mesh1 !< either BladeMotion or BladeRootMotion mesh + TYPE(MeshType), intent(in) :: HubMotion !< HubMotion mesh + REAL(R8Ki), intent(in) :: x_hat_disk(3) + REAL(ReKi), intent(out) :: orientationAnnulus(3,3,Mesh1%NNodes) + real(R8Ki), optional, intent(out) :: elemPosRelToHub_save( 3,Mesh1%NNodes) + real(R8Ki), optional, intent(out) :: elemPosRotorProj_save(3,Mesh1%NNodes) + + real(R8Ki) :: x_hat_annulus(3) ! rotor normal unit vector (local rotor reference frame) + real(R8Ki) :: y_hat_annulus(3) ! annulus tangent unit vector (local rotor reference frame) + real(R8Ki) :: z_hat_annulus(3) ! annulus radial unit vector (local rotor reference frame) +! real(R8Ki) :: chordVec(3) + + integer(intKi) :: j ! loop counter for nodes + + REAL(R8Ki) :: HubAbsPosition(3) + real(R8Ki) :: elemPosRelToHub(3) ! local copies of + real(R8Ki) :: elemPosRotorProj(3) ! local copies of + + + HubAbsPosition = HubMotion%Position(:,1) + HubMotion%TranslationDisp(:,1) + + !.......................... + ! orientation + !.......................... + + do j=1,Mesh1%NNodes + !chordVec(:,j) = Mesh1%orientation(:,2,j) + ! Project element position onto the rotor plane + elemPosRelToHub = Mesh1%Position(:,j) + Mesh1%TranslationDisp(:,j) - HubAbsPosition ! + 0.00_ReKi*chordVec(:,j)*p%BEMT%chord(j,k) + elemPosRotorProj = elemPosRelToHub - x_hat_disk * dot_product( x_hat_disk, elemPosRelToHub ) + + ! Get unit vectors of the local annulus reference frame + z_hat_annulus = elemPosRotorProj / TwoNorm( elemPosRotorProj ) + x_hat_annulus = x_hat_disk + y_hat_annulus = cross_product( z_hat_annulus, x_hat_annulus ) + + ! Form a orientation matrix for the annulus reference frame + orientationAnnulus(1,:,j) = x_hat_annulus + orientationAnnulus(2,:,j) = y_hat_annulus + orientationAnnulus(3,:,j) = z_hat_annulus + + if (present(elemPosRelToHub_save) ) elemPosRelToHub_save( :,j) = elemPosRelToHub + if (present(elemPosRotorProj_save)) elemPosRotorProj_save(:,j) = elemPosRotorProj + end do + +end subroutine Calculate_MeshOrientation_Rel2Hub +!---------------------------------------------------------------------------------------------------------------------------------- +subroutine GeomWithoutSweepPitchTwist(p,u,x_hat_disk,m,thetaBladeNds,ErrStat,ErrMsg) type(RotParameterType), intent(in ) :: p !< AD parameters type(RotInputType), intent(in ) :: u !< AD Inputs at Time + real(R8Ki), intent(in ) :: x_hat_disk(3) type(RotMiscVarType), intent(inout) :: m !< Misc/optimization variables real(R8Ki), intent( out) :: thetaBladeNds(p%NumBlNds,p%NumBlades) integer(IntKi), intent( out) :: ErrStat !< Error status of the operation @@ -1945,6 +1999,7 @@ subroutine GeomWithoutSweepPitchTwist(p,u,m,thetaBladeNds,ErrStat,ErrMsg) ErrMsg = "" if (p%AeroProjMod==0) then + ! theta, "Twist angle (includes all sources of twist)" rad ! Vx, "Local axial velocity at node" m/s ! Vy, "Local tangential velocity at node" m/s @@ -1974,22 +2029,18 @@ subroutine GeomWithoutSweepPitchTwist(p,u,m,thetaBladeNds,ErrStat,ErrMsg) call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) theta = EulerExtract( orientation ) !root(k)WithoutPitch_theta(j)_blade(k) - thetaBladeNds(j,k) = -theta(3) ! local pitch + twist (aerodyanmic + elastic) angle of the jth node in the kth blade - + m%Curve( j,k) = theta(2) ! save value for possible output later + thetaBladeNds(j,k) = -theta(3) ! local pitch + twist (aerodyanmic + elastic) angle of the jth node in the kth blade theta(1) = 0.0_ReKi theta(3) = 0.0_ReKi - m%Curve(j,k) = theta(2) ! save value for possible output later m%WithoutSweepPitchTwist(:,:,j,k) = matmul( EulerConstruct( theta ), orientation_nopitch ) ! WithoutSweepPitch+Twist_theta(j)_Blade(k) end do !j=nodes end do !k=blades + else if (p%AeroProjMod==1) then - !m%AllOuts( BPitch( k) ) = 0.0_ReKi ! save this value of pitch for potential output; ill-defined, TODO - ! - !m%hub_theta_x_root(k) = 0.0_ReKi ! ill-defined, TODO - do k=1,p%NumBlades call LAPACK_gemm( 'n', 't', 1.0_R8Ki, u%BladeRootMotion(k)%Orientation(:,:,1), u%HubMotion%Orientation(:,:,1), 0.0_R8Ki, orientation, errStat2, errMsg2) call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) @@ -2000,16 +2051,25 @@ subroutine GeomWithoutSweepPitchTwist(p,u,m,thetaBladeNds,ErrStat,ErrMsg) theta(3) = 0.0_ReKi m%hub_theta_x_root(k) = theta(1) ! save this value for FAST.Farm end do - - ! Generic blade, we don't assume where the axes are, and we keep the default orientation + + do k=1,p%NumBlades - + call Calculate_MeshOrientation_Rel2Hub(u%BladeMotion(k), u%HubMotion, x_hat_disk, m%WithoutSweepPitchTwist(:,:,:,k)) + do j=1,p%NumBlNds - thetaBladeNds(j,k) = 0.0_ReKi ! local pitch + twist (aerodyanmic + elastic) angle of the jth node in the kth blade - m%Curve(j,k) = 0.0_ReKi ! ill-defined, TODO m%WithoutSweepPitchTwist(:,:,j,k) = u%BladeMotion(k)%Orientation(:,:,j) enddo enddo + + do k=1,p%NumBlades + do j=1,p%NumBlNds + orientation = matmul( u%BladeMotion(k)%Orientation(:,:,j), transpose( m%WithoutSweepPitchTwist(:,:,j,k) ) ) + theta = EulerExtract( orientation ) + m%Curve( j,k) = theta(2) + thetaBladeNds(j,k) = -theta(3) + enddo + enddo + else ErrStat = ErrID_Fatal ErrMsg ='GeomWithoutSweepPitchTwist: AeroProjMod not supported '//trim(num2lstr(p%AeroProjMod)) @@ -2040,7 +2100,7 @@ subroutine SetInputsForFVW(p, u, m, errStat, errMsg) ! Get disk average values and orientations ! NOTE: needed because it sets m%V_diskAvg and m%V_dot_x, needed by CalcOutput.. call DiskAvgValues(p%rotors(iR), u(tIndx)%rotors(iR), m%rotors(iR), x_hat_disk) ! also sets m%V_diskAvg and m%V_dot_x - call GeomWithoutSweepPitchTwist(p%rotors(iR),u(tIndx)%rotors(iR), m%rotors(iR), thetaBladeNds,ErrStat,ErrMsg) + call GeomWithoutSweepPitchTwist(p%rotors(iR),u(tIndx)%rotors(iR), x_hat_disk, m%rotors(iR), thetaBladeNds,ErrStat,ErrMsg) if (ErrStat >= AbortErrLev) return ! Rather than use a meshcopy, we will just copy what we need to the WingsMesh @@ -2310,7 +2370,7 @@ subroutine SetOutputsFromFVW(t, u, p, OtherState, x, xd, m, y, ErrStat, ErrMsg) end subroutine SetOutputsFromFVW !---------------------------------------------------------------------------------------------------------------------------------- -!> This routine validates the inputs from the AeroDyn input files. +!> This routine validates the number of blades on each rotor. SUBROUTINE ValidateNumBlades( NumBl, ErrStat, ErrMsg ) integer(IntKi), intent(in) :: NumBl !< Number of blades integer(IntKi), intent(out) :: ErrStat !< Error status @@ -2342,6 +2402,17 @@ SUBROUTINE ValidateInputData( InitInp, InputFileData, NumBl, ErrStat, ErrMsg ) ErrStat = ErrID_None ErrMsg = "" +! do iR = 1,size(NumBl) +! if (NumBl(iR) < 1) then +! call SetErrStat( ErrID_Fatal, 'Number of blades must be at least 1.', ErrStat, ErrMsg, RoutineName ) +! return ! return early because InputFileData%BladeProps may not be allocated properly otherwise... +! else +! if (NumBl(iR) > AD_MaxBl_Out .and. InitInp%Linearize) then +! call SetErrStat( ErrID_Fatal, 'Number of blades must be no larger than '//trim(num2lstr(AD_MaxBl_Out))//' for linearizaton analysis.', ErrStat, ErrMsg, RoutineName ) +! return ! return early because InputFileData%BladeProps may not be allocated properly otherwise... +! end if +! end if +! end do if (InputFileData%DTAero <= 0.0) call SetErrStat ( ErrID_Fatal, 'DTAero must be greater than zero.', ErrStat, ErrMsg, RoutineName ) if (InputFileData%WakeMod /= WakeMod_None .and. InputFileData%WakeMod /= WakeMod_BEMT .and. InputFileData%WakeMod /= WakeMod_DBEMT .and. InputFileData%WakeMod /= WakeMod_FVW) then @@ -2748,6 +2819,7 @@ SUBROUTINE Init_BEMTmodule( InputFileData, RotInputFileData, u_AD, u, p, p_AD, x real(ReKi) :: tmp(3), tmp_sz_y, tmp_sz real(ReKi) :: y_hat_disk(3) real(ReKi) :: z_hat_disk(3) + real(ReKi) :: position(3) real(ReKi) :: rMax real(ReKi) :: frac integer(IntKi) :: ErrStat2 @@ -2786,7 +2858,7 @@ SUBROUTINE Init_BEMTmodule( InputFileData, RotInputFileData, u_AD, u, p, p_AD, x call AllocAry(InitInp%zLocal,InitInp%numBladeNodes,InitInp%numBlades,'zLocal', ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) call AllocAry(InitInp%rLocal,InitInp%numBladeNodes,InitInp%numBlades,'rLocal', ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) call AllocAry(InitInp%zTip, InitInp%numBlades,'zTip', ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) - + call AllocAry(InitInp%rTipFix, InitInp%numBlades,'rTipFix',ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) call AllocAry(InitInp%UAOff_innerNode, InitInp%numBlades,'UAOff_innerNode',ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) call AllocAry(InitInp%UAOff_outerNode, InitInp%numBlades,'UAOff_outerNode',ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) @@ -2796,7 +2868,7 @@ SUBROUTINE Init_BEMTmodule( InputFileData, RotInputFileData, u_AD, u, p, p_AD, x end if - ! Compute zLocal, zHub, zTip, rLocal, rMax + ! Compute zLocal, zHub, zTip, rLocal, rMax, rTipFix rMax = 0.0_ReKi do k=1,p%numBlades @@ -2823,7 +2895,25 @@ SUBROUTINE Init_BEMTmodule( InputFileData, RotInputFileData, u_AD, u, p, p_AD, x tmp_sz = dot_product( tmp, z_hat_disk )**2 InitInp%rLocal(j,k) = sqrt( tmp_sz + tmp_sz_y ) rMax = max(rMax, InitInp%rLocal(j,k)) - end do !j=nodes + end do !j=nodes + + + !......... + ! compute fixed rLocal at tip node (without prebend) for Bladed-like calculations: + !......... + tmp(1) = 0.0_ReKi !RotInputFile%BladeProps(k)%BlCrvAC(p%NumBlNds) + tmp(2) = 0.0_ReKi !RotInputFile%BladeProps(k)%BlSwpAC(p%NumBlNds) + tmp(3) = RotInputFileData%BladeProps(k)%BlSpn(p%NumBlNds) + position = u_AD%BladeRootMotion(k)%Position(:,1) + matmul(tmp,u_AD%BladeRootMotion(k)%RefOrientation(:,:,1)) ! note that because positionL is a 1-D array, we're doing the transpose of matmul(transpose(u%BladeRootMotion(k)%RefOrientation),positionL) + + ! position of the coned tip node in the kth blade relative to the hub: + tmp = position - u_AD%HubMotion%Position(:,1) + + ! local radius (normalized distance from rotor centerline) + tmp_sz_y = dot_product( tmp, y_hat_disk )**2 + tmp_sz = dot_product( tmp, z_hat_disk )**2 + InitInp%rTipFix(k) = sqrt( tmp_sz + tmp_sz_y ) + end do !k=blades @@ -3313,7 +3403,7 @@ FUNCTION CalculateTowerInfluence(p, xbar, ybar, zbar, W_tower, TwrCd, TwrTI) RES u_TwrShadow = -TwrCd / denom * cos( PiBy2*ybar / denom )**2 end if end if - case (TwrShadow_Eames) + case (TwrShadow_Eames) if ( xbar > 0.0_ReKi .and. abs(zbar) < 1.0_ReKi) then exponential = ( ybar / (TwrTI * xbar) )**2 denom = TwrTI * xbar * sqrt( TwoPi ) @@ -3395,6 +3485,9 @@ SUBROUTINE getLocalTowerProps(p, u, BladeNodePosition, theta_tower_trans, W_towe ! Inside the tower, or very close, (will happen for vortex elements) we keep undisturbed inflow ! We don't want to reach the stagnation points DisturbInflow = .false. + !elseif ( TwrClrnc<= 0.0_ReKi) then + ! ! Tower strike + ! DisturbInflow = .false. else DisturbInflow = .true. end if @@ -3884,7 +3977,7 @@ SUBROUTINE Rot_JacobianPInput( t, u, p, p_AD, x, xd, z, OtherState, y, m, m_AD, ! get central difference: - call Compute_dY( p, y_p, y_m, delta_p, delta_m, dYdu(:,i) ) + call Compute_dY( p, p_AD, y_p, y_m, delta_p, delta_m, dYdu(:,i) ) end do @@ -4166,7 +4259,7 @@ SUBROUTINE RotJacobianPContState( t, u, p, p_AD, x, xd, z, OtherState, y, m, m_A ! get central difference: - call Compute_dY( p, y_p, y_m, delta_p, delta_m, dYdx(:,i) ) + call Compute_dY( p, p_AD, y_p, y_m, delta_p, delta_m, dYdx(:,i) ) end do @@ -4533,7 +4626,7 @@ SUBROUTINE RotJacobianPConstrState( t, u, p, p_AD, x, xd, z, OtherState, y, m, m ! get central difference: - call Compute_dY( p, y_p, y_m, delta_p, delta_m, dYdz(:,i) ) + call Compute_dY( p, p_AD, y_p, y_m, delta_p, delta_m, dYdz(:,i) ) ! put z_perturb back (for next iteration): @@ -5645,9 +5738,10 @@ END SUBROUTINE Perturb_x !---------------------------------------------------------------------------------------------------------------------------------- !> This routine uses values of two output types to compute an array of differences. !! Do not change this packing without making sure subroutine aerodyn::init_jacobian is consistant with this routine! -SUBROUTINE Compute_dY(p, y_p, y_m, delta_p, delta_m, dY) +SUBROUTINE Compute_dY(p, p_AD, y_p, y_m, delta_p, delta_m, dY) TYPE(RotParameterType) , INTENT(IN ) :: p !< parameters + TYPE(AD_ParameterType) , INTENT(IN ) :: p_AD !< parameters TYPE(RotOutputType) , INTENT(IN ) :: y_p !< AD outputs at \f$ u + \Delta_p u \f$ or \f$ x + \Delta_p x \f$ (p=plus) TYPE(RotOutputType) , INTENT(IN ) :: y_m !< AD outputs at \f$ u - \Delta_m u \f$ or \f$ x - \Delta_m x \f$ (m=minus) REAL(R8Ki) , INTENT(IN ) :: delta_p !< difference in inputs or states \f$ delta_p = \Delta_p u \f$ or \f$ delta_p = \Delta_p x \f$ diff --git a/modules/aerodyn/src/AeroDyn_IO.f90 b/modules/aerodyn/src/AeroDyn_IO.f90 index 9273b345ec..e31a73161d 100644 --- a/modules/aerodyn/src/AeroDyn_IO.f90 +++ b/modules/aerodyn/src/AeroDyn_IO.f90 @@ -1654,89 +1654,98 @@ END FUNCTION Calc_Chi0 !---------------------------------------------------------------------------------------------------------------------------------- -SUBROUTINE Calc_WriteOutput( p, p_AD, u, m, m_AD, y, OtherState, xd, indx, iRot, ErrStat, ErrMsg ) +SUBROUTINE Calc_WriteOutput( p, p_AD, u, x, m, m_AD, y, OtherState, xd, indx, iRot, ErrStat, ErrMsg ) - TYPE(RotParameterType), INTENT(IN ) :: p ! The rotor parameters - TYPE(AD_ParameterType), INTENT(IN ) :: p_AD ! The module parameters - TYPE(RotInputType), INTENT(IN ) :: u ! inputs - TYPE(RotMiscVarType), INTENT(INOUT) :: m ! misc variables - TYPE(AD_MiscVarType), INTENT(INOUT) :: m_AD ! misc variables - TYPE(RotOutputType), INTENT(IN ) :: y ! outputs - TYPE(RotOtherStateType), INTENT(IN ) :: OtherState ! other states at t (for DBEMT and UA) - TYPE(RotDiscreteStateType),INTENT(IN ) :: xd ! Discrete states - integer, intent(in ) :: indx ! index into m%BEMT_u(indx) array; 1=t and 2=t+dt (but not checked here) - integer, intent(in ) :: iRot ! Rotor index, needed for FVW - INTEGER(IntKi), INTENT( OUT) :: ErrStat ! The error status code - CHARACTER(*), INTENT( OUT) :: ErrMsg ! The error message, if an error occurred + TYPE(RotParameterType), INTENT(IN ) :: p ! The rotor parameters + TYPE(AD_ParameterType), INTENT(IN ) :: p_AD ! The module parameters + TYPE(RotInputType), INTENT(IN ) :: u ! inputs + TYPE(RotContinuousStateType), INTENT(IN ) :: x !< Continuous states at t + TYPE(RotMiscVarType), INTENT(INOUT) :: m ! misc variables + TYPE(AD_MiscVarType), INTENT(INOUT) :: m_AD ! misc variables + TYPE(RotOutputType), INTENT(IN ) :: y ! outputs + TYPE(RotOtherStateType), INTENT(IN ) :: OtherState ! other states at t (for DBEMT and UA) + TYPE(RotDiscreteStateType), INTENT(IN ) :: xd ! Discrete states + integer, intent(in ) :: indx ! index into m%BEMT_u(indx) array; 1=t and 2=t+dt (but not checked here) + integer, intent(in ) :: iRot ! Rotor index, needed for FVW + INTEGER(IntKi), INTENT( OUT) :: ErrStat ! The error status code + CHARACTER(*), INTENT( OUT) :: ErrMsg ! The error message, if an error occurred ! local variables - CHARACTER(*), PARAMETER :: RoutineName = 'Calc_WriteOutput' - INTEGER(intKi) :: ErrStat2 - CHARACTER(ErrMsgLen) :: ErrMsg2 - - INTEGER(IntKi) :: j,k,beta - REAL(ReKi) :: tmp(3) - REAL(ReKi) :: force(3) - REAL(ReKi) :: moment(3) - REAL(ReKi) :: denom, rmax - REAL(ReKi) :: ct, st ! cosine, sine of theta - REAL(ReKi) :: cp, sp ! cosine, sine of phi - + CHARACTER(*), PARAMETER :: RoutineName = 'Calc_WriteOutput' + INTEGER(intKi) :: ErrStat2 + CHARACTER(ErrMsgLen) :: ErrMsg2 + INTEGER(IntKi) :: j,k,beta + REAL(ReKi) :: tmp(3) + REAL(ReKi) :: force(3) + REAL(ReKi) :: moment(3) + REAL(ReKi) :: denom, rmax, omega + REAL(ReKi) :: ct, st ! cosine, sine of theta + REAL(ReKi) :: cp, sp ! cosine, sine of phi + ! start routine: ErrStat = ErrID_None ErrMsg = "" - ! tower outputs - do beta=1,p%NTwOuts - j = p%TwOutNd(beta) - - tmp = matmul( u%TowerMotion%Orientation(:,:,j) , u%InflowOnTower(:,j) ) - m%AllOuts( TwNVUnd(:,beta) ) = tmp - - tmp = matmul( u%TowerMotion%Orientation(:,:,j) , u%TowerMotion%TranslationVel(:,j) ) - m%AllOuts( TwNSTV(:,beta) ) = tmp - - m%AllOuts( TwNVrel(beta) ) = m%W_Twr(j) ! relative velocity - m%AllOuts( TwNDynP(beta) ) = 0.5 * p%AirDens * m%W_Twr(j)**2 ! dynamic pressure - m%AllOuts( TwNRe( beta) ) = p%TwrDiam(j) * m%W_Twr(j) / p%KinVisc / 1.0E6 ! reynolds number (in millions) - m%AllOuts( TwNM( beta) ) = m%W_Twr(j) / p%SpdSound ! Mach number - m%AllOuts( TwNFdx( beta) ) = m%X_Twr(j) - m%AllOuts( TwNFdy( beta) ) = m%Y_Twr(j) - - end do ! out nodes + ! Compute max radius and rotor speed if (p_AD%WakeMod /= WakeMod_FVW) then - call Calc_WriteOutput_BEMT + rmax = 0.0_ReKi + do k=1,p%NumBlades + do j=1,p%NumBlNds + rmax = max(rmax, m%BEMT_u(indx)%rLocal(j,k) ) + end do !j=nodes + end do !k=blades + +! rmax = p%BEMT%rTipFixMax + omega = m%BEMT_u(indx)%omega else - call Calc_WriteOutput_FVW + rmax = Calc_MaxRadius(p, u) + omega = Calc_Omega(u) endif + + call Calc_WriteOutput_AD() ! need to call this before calling the BEMT vs FVW versions of outputs so that the integrated output quantities are known + + if (p_AD%WakeMod /= WakeMod_FVW) then + call Calc_WriteOutput_BEMT() + else + call Calc_WriteOutput_FVW() + endif - ! blade node tower clearance (requires tower influence calculation): - if (p%TwrPotent /= TwrPotent_none .or. p%TwrShadow /= TwrShadow_none) then - do k=1,p%numBlades - do beta=1,p%NBlOuts - j=p%BlOutNd(beta) - m%AllOuts( BNClrnc( beta,k) ) = m%TwrClrnc(j,k) - end do - end do - end if + ! set these for debugging +! m%AllOuts( Debug1 ) = 0.0_ReKi !TwoNorm( m%BEMT%u_SkewWake(1)%v_qsw ) +! m%AllOuts( Debug2 ) = 0.0_ReKi !TwoNorm( x%BEMT%v_w ) +! m%AllOuts( Debug3 ) = 0.0_ReKi CONTAINS !.......................................................................................... - subroutine Calc_WriteOutput_BEMT - real(ReKi) :: omega - omega = m%BEMT_u(indx)%omega + subroutine Calc_WriteOutput_AD() + + ! tower outputs + do beta=1,p%NTwOuts + j = p%TwOutNd(beta) + + tmp = matmul( u%TowerMotion%Orientation(:,:,j) , u%InflowOnTower(:,j) ) + m%AllOuts( TwNVUnd(:,beta) ) = tmp + + tmp = matmul( u%TowerMotion%Orientation(:,:,j) , u%TowerMotion%TranslationVel(:,j) ) + m%AllOuts( TwNSTV(:,beta) ) = tmp + + m%AllOuts( TwNVrel(beta) ) = m%W_Twr(j) ! relative velocity + m%AllOuts( TwNDynP(beta) ) = 0.5 * p%AirDens * m%W_Twr(j)**2 ! dynamic pressure + m%AllOuts( TwNRe( beta) ) = p%TwrDiam(j) * m%W_Twr(j) / p%KinVisc / 1.0E6 ! reynolds number (in millions) + m%AllOuts( TwNM( beta) ) = m%W_Twr(j) / p%SpdSound ! Mach number + m%AllOuts( TwNFdx( beta) ) = m%X_Twr(j) + m%AllOuts( TwNFdy( beta) ) = m%Y_Twr(j) + + end do ! out nodes + ! blade outputs do k=1,min(p%numBlades,3) ! limit this - m%AllOuts( BAzimuth(k) ) = MODULO( m%BEMT_u(indx)%psi(k)*R2D, 360.0_ReKi ) - ! m%AllOuts( BPitch( k) ) = calculated in SetInputsForBEMT - do beta=1,p%NBlOuts - j=p%BlOutNd(beta) tmp = matmul( m%WithoutSweepPitchTwist(:,:,j,k), u%InflowOnBlade(:,j,k) ) @@ -1753,72 +1762,39 @@ subroutine Calc_WriteOutput_BEMT m%AllOuts( BNSTVx( beta,k) ) = tmp(1) m%AllOuts( BNSTVy( beta,k) ) = tmp(2) m%AllOuts( BNSTVz( beta,k) ) = tmp(3) - - m%AllOuts( BNVrel( beta,k) ) = m%BEMT_y%Vrel(j,k) - m%AllOuts( BNDynP( beta,k) ) = 0.5 * p%airDens * m%BEMT_y%Vrel(j,k)**2 - m%AllOuts( BNRe( beta,k) ) = p%BEMT%chord(j,k) * m%BEMT_y%Vrel(j,k) / p%KinVisc / 1.0E6 - m%AllOuts( BNM( beta,k) ) = m%BEMT_y%Vrel(j,k) / p%SpdSound - - m%AllOuts( BNVIndx(beta,k) ) = - m%BEMT_u(indx)%Vx(j,k) * m%BEMT_y%axInduction( j,k) - m%AllOuts( BNVIndy(beta,k) ) = m%BEMT_u(indx)%Vy(j,k) * m%BEMT_y%tanInduction(j,k) - - m%AllOuts( BNAxInd(beta,k) ) = m%BEMT_y%axInduction(j,k) - m%AllOuts( BNTnInd(beta,k) ) = m%BEMT_y%tanInduction(j,k) - - m%AllOuts( BNAlpha(beta,k) ) = Rad2M180to180Deg( m%BEMT_y%phi(j,k) - m%BEMT_u(indx)%theta(j,k) ) - m%AllOuts( BNTheta(beta,k) ) = m%BEMT_u(indx)%theta(j,k)*R2D - m%AllOuts( BNPhi( beta,k) ) = m%BEMT_y%phi(j,k)*R2D + m%AllOuts( BNCurve(beta,k) ) = m%Curve(j,k)*R2D - - !m%AllOuts( BNCl( beta,k) ) = m%BEMT_y%Cl(j,k) - !m%AllOuts( BNCd( beta,k) ) = m%BEMT_y%Cd(j,k) - - m%AllOuts( BNCpmin( beta,k) ) = m%BEMT_y%Cpmin(j,k) + m%AllOuts( BNSigCr( beta,k) ) = m%SigmaCavitCrit(j,k) m%AllOuts( BNSgCav( beta,k) ) = m%SigmaCavit(j,k) - cp=cos(m%BEMT_y%phi(j,k)) - sp=sin(m%BEMT_y%phi(j,k)) - m%AllOuts( BNCl( beta,k) ) = m%BEMT_y%Cx(j,k)*cp + m%BEMT_y%Cy(j,k)*sp - m%AllOuts( BNCd( beta,k) ) = m%BEMT_y%Cx(j,k)*sp - m%BEMT_y%Cy(j,k)*cp - m%AllOuts( BNCm( beta,k) ) = m%BEMT_y%Cm(j,k) - m%AllOuts( BNCx( beta,k) ) = m%BEMT_y%Cx(j,k) - m%AllOuts( BNCy( beta,k) ) = m%BEMT_y%Cy(j,k) - - ct=cos(m%BEMT_u(indx)%theta(j,k)) - st=sin(m%BEMT_u(indx)%theta(j,k)) - m%AllOuts( BNCn( beta,k) ) = m%BEMT_y%Cx(j,k)*ct + m%BEMT_y%Cy(j,k)*st - m%AllOuts( BNCt( beta,k) ) =-m%BEMT_y%Cx(j,k)*st + m%BEMT_y%Cy(j,k)*ct - - m%AllOuts( BNFl( beta,k) ) = m%X(j,k)*cp - m%Y(j,k)*sp - m%AllOuts( BNFd( beta,k) ) = m%X(j,k)*sp + m%Y(j,k)*cp - m%AllOuts( BNMm( beta,k) ) = m%M(j,k) - m%AllOuts( BNFx( beta,k) ) = m%X(j,k) - m%AllOuts( BNFy( beta,k) ) = -m%Y(j,k) - m%AllOuts( BNFn( beta,k) ) = m%X(j,k)*ct - m%Y(j,k)*st - m%AllOuts( BNFt( beta,k) ) = -m%X(j,k)*st - m%Y(j,k)*ct - - m%AllOuts( BNGam( beta,k) ) = 0.5_ReKi * p%BEMT%chord(j,k) * m%BEMT_y%Vrel(j,k) * m%BEMT_y%Cl(j,k) ! "Gam" [m^2/s] end do ! nodes end do ! blades + - ! rotor outputs: - rmax = 0.0_ReKi - do k=1,p%NumBlades - do j=1,p%NumBlNds - rmax = max(rmax, m%BEMT_u(indx)%rLocal(j,k) ) - end do !j=nodes - end do !k=blades - m%AllOuts( RtSpeed ) = m%BEMT_u(indx)%omega*RPS2RPM - m%AllOuts( RtArea ) = pi*rmax**2 + ! blade node tower clearance (requires tower influence calculation): + if (p%TwrPotent /= TwrPotent_none .or. p%TwrShadow /= TwrShadow_none) then + do k=1,p%numBlades + do beta=1,p%NBlOuts + j=p%BlOutNd(beta) + m%AllOuts( BNClrnc( beta,k) ) = m%TwrClrnc(j,k) + end do + end do + end if + + + + m%AllOuts( RtSpeed ) = omega*RPS2RPM + m%AllOuts( RtArea ) = pi * rmax**2 + tmp = matmul( u%HubMotion%Orientation(:,:,1), m%V_DiskAvg ) m%AllOuts( RtVAvgxh ) = tmp(1) m%AllOuts( RtVAvgyh ) = tmp(2) m%AllOuts( RtVAvgzh ) = tmp(3) - m%AllOuts( RtSkew ) = m%BEMT_u(indx)%chi0*R2D + ! integrate force/moments over blades by performing mesh transfer to hub point: force = 0.0_ReKi @@ -1827,76 +1803,155 @@ subroutine Calc_WriteOutput_BEMT call Transfer_Line2_to_Point( y%BladeLoad(k), m%HubLoad, m%B_L_2_H_P(k), ErrStat2, ErrMsg2, u%BladeMotion(k), u%HubMotion ) force = force + m%HubLoad%force( :,1) moment = moment + m%HubLoad%moment(:,1) - - if (k<=4) then + + if (k<=size(BAeroFxg)) then ! Power contribution of blade wrt hub tmp = matmul( u%HubMotion%Orientation(:,:,1), m%HubLoad%moment(:,1) ) - m%AllOuts( BAeroPwr(k) ) = omega * tmp(1) - - ! In global, wrt hub! + m%AllOuts( BAeroPwr(k) ) = omega * tmp(1) + + ! In global, wrt hub! m%AllOuts( BAeroFxg(k) ) = m%HubLoad%force(1,1) m%AllOuts( BAeroFyg(k) ) = m%HubLoad%force(2,1) m%AllOuts( BAeroFzg(k) ) = m%HubLoad%force(3,1) m%AllOuts( BAeroMxg(k) ) = m%HubLoad%moment(1,1) m%AllOuts( BAeroMyg(k) ) = m%HubLoad%moment(2,1) m%AllOuts( BAeroMzg(k) ) = m%HubLoad%moment(3,1) - - endif + end if end do - ! In global + + ! In global m%AllOuts( RtAeroFxg ) = force(1) m%AllOuts( RtAeroFyg ) = force(2) m%AllOuts( RtAeroFzg ) = force(3) m%AllOuts( RtAeroMxg ) = moment(1) m%AllOuts( RtAeroMyg ) = moment(2) m%AllOuts( RtAeroMzg ) = moment(3) - ! In hub coord tmp = matmul( u%HubMotion%Orientation(:,:,1), force ) m%AllOuts( RtAeroFxh ) = tmp(1) m%AllOuts( RtAeroFyh ) = tmp(2) m%AllOuts( RtAeroFzh ) = tmp(3) - + tmp = matmul( u%HubMotion%Orientation(:,:,1), moment ) m%AllOuts( RtAeroMxh ) = tmp(1) m%AllOuts( RtAeroMyh ) = tmp(2) m%AllOuts( RtAeroMzh ) = tmp(3) - - m%AllOuts( RtAeroPwr ) = m%BEMT_u(indx)%omega * m%AllOuts( RtAeroMxh ) - - - - m%AllOuts( RtTSR ) = m%BEMT_u(indx)%TSR - - if ( EqualRealNos( m%V_dot_x, 0.0_ReKi ) ) then - m%AllOuts( RtTSR ) = 0.0_ReKi - m%AllOuts( RtAeroCp ) = 0.0_ReKi - m%AllOuts( RtAeroCq ) = 0.0_ReKi - m%AllOuts( RtAeroCt ) = 0.0_ReKi - else - denom = 0.5*p%AirDens*m%AllOuts( RtArea )*m%V_dot_x**2 - m%AllOuts( RtTSR ) = m%BEMT_u(indx)%omega * rmax / m%V_dot_x - - m%AllOuts( RtAeroCp ) = m%AllOuts( RtAeroPwr ) / (denom * m%V_dot_x) - m%AllOuts( RtAeroCq ) = m%AllOuts( RtAeroMxh ) / (denom * rmax) - m%AllOuts( RtAeroCt ) = m%AllOuts( RtAeroFxh ) / denom - end if + + m%AllOuts( RtAeroPwr ) = omega * m%AllOuts( RtAeroMxh ) + + ! Integrate force/moments over blades by performing mesh transfer to blade root points: - do k=1,min(p%NumBlades,4) ! Temporary hack for at least one more blae outputs + do k=1,p%NumBlades call Transfer_Line2_to_Point( y%BladeLoad(k), m%BladeRootLoad(k), m%B_L_2_R_P(k), ErrStat2, ErrMsg2, u%BladeMotion(k), u%BladeRootMotion(k) ) + end do + do k=1,min(p%NumBlades,size(BAeroFx)) ! Transform force vector to blade root coordinate system tmp = matmul( u%BladeRootMotion(k)%Orientation(:,:,1), m%BladeRootLoad(k)%force( :,1) ) m%AllOuts( BAeroFx(k) ) = tmp(1) m%AllOuts( BAeroFy(k) ) = tmp(2) m%AllOuts( BAeroFz(k) ) = tmp(3) + ! Transform moment vector to blade root coordinate system tmp = matmul( u%BladeRootMotion(k)%Orientation(:,:,1), m%BladeRootLoad(k)%moment( :,1) ) m%AllOuts( BAeroMx(k) ) = tmp(1) m%AllOuts( BAeroMy(k) ) = tmp(2) m%AllOuts( BAeroMz(k) ) = tmp(3) - end do ! k=blades + end do ! k=blades + + ! rotor outputs + if ( EqualRealNos( m%V_dot_x, 0.0_ReKi ) ) then + m%AllOuts( RtTSR ) = 0.0_ReKi + m%AllOuts( RtAeroCp ) = 0.0_ReKi + m%AllOuts( RtAeroCq ) = 0.0_ReKi + m%AllOuts( RtAeroCt ) = 0.0_ReKi + else + m%AllOuts( RtTSR ) = omega * rmax / m%V_dot_x + + denom = 0.5*p%AirDens*m%AllOuts( RtArea )*m%V_dot_x**2 + m%AllOuts( RtAeroCp ) = m%AllOuts( RtAeroPwr ) / (denom * m%V_dot_x) + m%AllOuts( RtAeroCq ) = m%AllOuts( RtAeroMxh ) / (denom * rmax ) + m%AllOuts( RtAeroCt ) = m%AllOuts( RtAeroFxh ) / denom + end if + + + end subroutine Calc_WriteOutput_AD + !.......................................................................................... + subroutine Calc_WriteOutput_BEMT() + REAL(R8Ki) :: orient(3,3) + REAL(R8Ki) :: theta(3) + REAL(ReKi) :: denom !, rmax + REAL(ReKi) :: ct, st ! cosine, sine of theta + REAL(ReKi) :: cp, sp ! cosine, sine of phi + + + + + ! blade outputs + do k=1,min(p%numBlades,size(BAzimuth) ) ! limit this + m%AllOuts( BAzimuth(k) ) = MODULO( m%BEMT_u(indx)%psi(k)*R2D, 360.0_ReKi ) + ! m%AllOuts( BPitch( k) ) = calculated in SetInputsForBEMT + + do beta=1,p%NBlOuts + + j=p%BlOutNd(beta) + + m%AllOuts( BNVrel( beta,k) ) = m%BEMT_y%Vrel(j,k) + m%AllOuts( BNDynP( beta,k) ) = 0.5 * p%airDens * m%BEMT_y%Vrel(j,k)**2 + m%AllOuts( BNRe( beta,k) ) = p%BEMT%chord(j,k) * m%BEMT_y%Vrel(j,k) / p%KinVisc / 1.0E6 + m%AllOuts( BNM( beta,k) ) = m%BEMT_y%Vrel(j,k) / p%SpdSound + + m%AllOuts( BNVIndx(beta,k) ) = - m%BEMT_u(indx)%Vx(j,k) * m%BEMT_y%axInduction( j,k) + m%AllOuts( BNVIndy(beta,k) ) = m%BEMT_u(indx)%Vy(j,k) * m%BEMT_y%tanInduction(j,k) + + m%AllOuts( BNAxInd(beta,k) ) = m%BEMT_y%axInduction(j,k) + m%AllOuts( BNTnInd(beta,k) ) = m%BEMT_y%tanInduction(j,k) + + m%AllOuts( BNAlpha(beta,k) ) = Rad2M180to180Deg( m%BEMT_y%phi(j,k) - m%BEMT_u(indx)%theta(j,k) ) + m%AllOuts( BNTheta(beta,k) ) = m%BEMT_u(indx)%theta(j,k)*R2D + m%AllOuts( BNPhi( beta,k) ) = m%BEMT_y%phi(j,k)*R2D + + ! m%AllOuts( BNCurve(beta,k) ) = m%Curve(j,k)*R2D + + m%AllOuts( BNCpmin( beta,k) ) = m%BEMT_y%Cpmin(j,k) + ! m%AllOuts( BNSigCr( beta,k) ) = m%SigmaCavitCrit(j,k) + ! m%AllOuts( BNSgCav( beta,k) ) = m%SigmaCavit(j,k) + + !m%AllOuts( BNCl( beta,k) ) = m%BEMT_y%Cl(j,k) + !m%AllOuts( BNCd( beta,k) ) = m%BEMT_y%Cd(j,k) + cp=cos(m%BEMT_y%phi(j,k)) + sp=sin(m%BEMT_y%phi(j,k)) + m%AllOuts( BNCl( beta,k) ) = m%BEMT_y%Cx(j,k)*cp + m%BEMT_y%Cy(j,k)*sp + m%AllOuts( BNCd( beta,k) ) = m%BEMT_y%Cx(j,k)*sp - m%BEMT_y%Cy(j,k)*cp + m%AllOuts( BNCm( beta,k) ) = m%BEMT_y%Cm(j,k) + m%AllOuts( BNCx( beta,k) ) = m%BEMT_y%Cx(j,k) + m%AllOuts( BNCy( beta,k) ) = m%BEMT_y%Cy(j,k) + + ct=cos(m%BEMT_u(indx)%theta(j,k)) + st=sin(m%BEMT_u(indx)%theta(j,k)) + m%AllOuts( BNCn( beta,k) ) = m%BEMT_y%Cx(j,k)*ct + m%BEMT_y%Cy(j,k)*st + m%AllOuts( BNCt( beta,k) ) =-m%BEMT_y%Cx(j,k)*st + m%BEMT_y%Cy(j,k)*ct + + m%AllOuts( BNFl( beta,k) ) = m%X(j,k)*cp - m%Y(j,k)*sp + m%AllOuts( BNFd( beta,k) ) = m%X(j,k)*sp + m%Y(j,k)*cp + m%AllOuts( BNMm( beta,k) ) = m%M(j,k) + m%AllOuts( BNFx( beta,k) ) = m%X(j,k) + m%AllOuts( BNFy( beta,k) ) = -m%Y(j,k) + m%AllOuts( BNFn( beta,k) ) = m%X(j,k)*ct - m%Y(j,k)*st + m%AllOuts( BNFt( beta,k) ) = -m%X(j,k)*st - m%Y(j,k)*ct + + m%AllOuts( BNGam( beta,k) ) = 0.5_ReKi * p%BEMT%chord(j,k) * m%BEMT_y%Vrel(j,k) * m%BEMT_y%Cl(j,k) ! "Gam" [m^2/s] + + end do ! nodes + end do ! blades + + + ! rotor outputs: + + m%AllOuts( RtSkew ) = m%BEMT_u(indx)%chi0*R2D +! m%AllOuts( RtTSR ) = m%BEMT_u(indx)%TSR m%AllOuts( DBEMTau1 ) = OtherState%BEMT%DBEMT%tau1 + end subroutine Calc_WriteOutput_BEMT @@ -1907,11 +1962,6 @@ end subroutine Calc_WriteOutput_BEMT !! Make sure these are set! subroutine Calc_WriteOutput_FVW integer :: iW - real(ReKi) :: rmax, omega - - ! Compute max radius and rotor speed - rmax = Calc_MaxRadius(p, u) - omega = Calc_Omega(u) ! blade outputs do k=1,min(p%numBlades,3) @@ -1920,22 +1970,6 @@ subroutine Calc_WriteOutput_FVW do beta=1,p%NBlOuts j=p%BlOutNd(beta) - ! --- Setting AD outputs - tmp = matmul( m%WithoutSweepPitchTwist(:,:,j,k), u%InflowOnBlade(:,j,k) ) - m%AllOuts( BNVUndx(beta,k) ) = tmp(1) - m%AllOuts( BNVUndy(beta,k) ) = tmp(2) - m%AllOuts( BNVUndz(beta,k) ) = tmp(3) - - tmp = matmul( m%WithoutSweepPitchTwist(:,:,j,k), m%DisturbedInflow(:,j,k) ) - m%AllOuts( BNVDisx(beta,k) ) = tmp(1) - m%AllOuts( BNVDisy(beta,k) ) = tmp(2) - m%AllOuts( BNVDisz(beta,k) ) = tmp(3) - - tmp = matmul( m%WithoutSweepPitchTwist(:,:,j,k), u%BladeMotion(k)%TranslationVel(:,j) ) - m%AllOuts( BNSTVx( beta,k) ) = tmp(1) - m%AllOuts( BNSTVy( beta,k) ) = tmp(2) - m%AllOuts( BNSTVz( beta,k) ) = tmp(3) - m%AllOuts( BNVrel( beta,k) ) = m_AD%FVW%W(iW)%BN_Vrel(j) m%AllOuts( BNDynP( beta,k) ) = 0.5 * p%airDens * m_AD%FVW%W(iW)%BN_Vrel(j)**2 m%AllOuts( BNRe( beta,k) ) = m_AD%FVW%W(iW)%BN_Re(j) / 1.0E6 @@ -1953,9 +1987,6 @@ subroutine Calc_WriteOutput_FVW ! m%AllOuts( BNCurve(beta,k) ) = m%Curve(j,k)*R2D ! TODO ! m%AllOuts( BNCpmin( beta,k) ) = m%BEMT_y%Cpmin(jk) ! TODO - m%AllOuts( BNSigCr( beta,k) ) = m%SigmaCavitCrit(j,k) - m%AllOuts( BNSgCav( beta,k) ) = m%SigmaCavit(j,k) - m%AllOuts( BNCl( beta,k) ) = m_AD%FVW%W(iW)%BN_Cl(j) m%AllOuts( BNCd( beta,k) ) = m_AD%FVW%W(iW)%BN_Cd(j) m%AllOuts( BNCm( beta,k) ) = m_AD%FVW%W(iW)%BN_Cm(j) @@ -1982,84 +2013,10 @@ subroutine Calc_WriteOutput_FVW end do ! blades - m%AllOuts( RtSpeed ) = omega*RPS2RPM - m%AllOuts( RtArea ) = pi*rmax**2 ! TODO vertical axis - - tmp = matmul( u%HubMotion%Orientation(:,:,1), m%V_DiskAvg ) - m%AllOuts( RtVAvgxh ) = tmp(1) - m%AllOuts( RtVAvgyh ) = tmp(2) - m%AllOuts( RtVAvgzh ) = tmp(3) - +! m%AllOuts( RtArea ) = pi*rmax**2 ! TODO vertical axis m%AllOuts( RtSkew ) = Calc_Chi0(m%V_diskAvg, m%V_dot_x) * R2D - ! integrate force/moments over blades by performing mesh transfer to hub point: - force = 0.0_ReKi - moment = 0.0_ReKi - do k=1,p%NumBlades - call Transfer_Line2_to_Point( y%BladeLoad(k), m%HubLoad, m%B_L_2_H_P(k), ErrStat2, ErrMsg2, u%BladeMotion(k), u%HubMotion ) - force = force + m%HubLoad%force( :,1) - moment = moment + m%HubLoad%moment(:,1) - - if (k<=4) then - ! Power contribution of blade wrt hub - tmp = matmul( u%HubMotion%Orientation(:,:,1), m%HubLoad%moment(:,1) ) - m%AllOuts( BAeroPwr(k) ) = omega * tmp(1) - ! In global, wrt hub! - m%AllOuts( BAeroFxg(k) ) = m%HubLoad%force(1,1) - m%AllOuts( BAeroFyg(k) ) = m%HubLoad%force(2,1) - m%AllOuts( BAeroFzg(k) ) = m%HubLoad%force(3,1) - m%AllOuts( BAeroMxg(k) ) = m%HubLoad%moment(1,1) - m%AllOuts( BAeroMyg(k) ) = m%HubLoad%moment(2,1) - m%AllOuts( BAeroMzg(k) ) = m%HubLoad%moment(3,1) - endif - end do - ! In global - m%AllOuts( RtAeroFxg ) = force(1) - m%AllOuts( RtAeroFyg ) = force(2) - m%AllOuts( RtAeroFzg ) = force(3) - m%AllOuts( RtAeroMxg ) = moment(1) - m%AllOuts( RtAeroMyg ) = moment(2) - m%AllOuts( RtAeroMzg ) = moment(3) - ! In hub coord - tmp = matmul( u%HubMotion%Orientation(:,:,1), force ) - m%AllOuts( RtAeroFxh ) = tmp(1) - m%AllOuts( RtAeroFyh ) = tmp(2) - m%AllOuts( RtAeroFzh ) = tmp(3) - - tmp = matmul( u%HubMotion%Orientation(:,:,1), moment ) - m%AllOuts( RtAeroMxh ) = tmp(1) - m%AllOuts( RtAeroMyh ) = tmp(2) - m%AllOuts( RtAeroMzh ) = tmp(3) - - m%AllOuts( RtAeroPwr ) = omega * m%AllOuts( RtAeroMxh ) - - if ( EqualRealNos( m%V_dot_x, 0.0_ReKi ) ) then - m%AllOuts( RtTSR ) = 0.0_ReKi - m%AllOuts( RtAeroCp ) = 0.0_ReKi - m%AllOuts( RtAeroCq ) = 0.0_ReKi - m%AllOuts( RtAeroCt ) = 0.0_ReKi - else - denom = 0.5*p%AirDens*m%AllOuts( RtArea )*m%V_dot_x**2 - m%AllOuts( RtTSR ) = omega * rmax / m%V_dot_x - m%AllOuts( RtAeroCp ) = m%AllOuts( RtAeroPwr ) / (denom * m%V_dot_x) - m%AllOuts( RtAeroCq ) = m%AllOuts( RtAeroMxh ) / (denom * rmax) - m%AllOuts( RtAeroCt ) = m%AllOuts( RtAeroFxh ) / denom - end if - - ! Integrate force/moments over blades by performing mesh transfer to blade root points: - do k=1,min(p%NumBlades,4) - call Transfer_Line2_to_Point( y%BladeLoad(k), m%BladeRootLoad(k), m%B_L_2_R_P(k), ErrStat2, ErrMsg2, u%BladeMotion(k), u%BladeRootMotion(k) ) - ! Transform force vector to blade root coordinate system - tmp = matmul( u%BladeRootMotion(k)%Orientation(:,:,1), m%BladeRootLoad(k)%force( :,1) ) - m%AllOuts( BAeroFx(k) ) = tmp(1) - m%AllOuts( BAeroFy(k) ) = tmp(2) - m%AllOuts( BAeroFz(k) ) = tmp(3) - ! Transform moment vector to blade root coordinate system - tmp = matmul( u%BladeRootMotion(k)%Orientation(:,:,1), m%BladeRootLoad(k)%moment( :,1) ) - m%AllOuts( BAeroMx(k) ) = tmp(1) - m%AllOuts( BAeroMy(k) ) = tmp(2) - m%AllOuts( BAeroMz(k) ) = tmp(3) - end do ! k=blades +! m%AllOuts( DBEMTau1 ) = 0.0_ReKi ! not valid with FVW end subroutine Calc_WriteOutput_FVW @@ -2123,7 +2080,7 @@ SUBROUTINE ReadInputFiles( InputFileName, InputFileData, Default_DT, OutFileRoot iBld = iBld+1 ! Increment blade counter END DO - ENDDO ! Loop on rotors + end do ! loop on rotors CALL Cleanup ( ) @@ -2661,8 +2618,12 @@ SUBROUTINE AD_PrintSum( InputFileData, p, p_AD, u, y, ErrStat, ErrMsg ) CHARACTER(*), PARAMETER :: FmtDatT = '(A,T35,1(:,F13.8))' ! Format for outputting time steps. CHARACTER(30) :: OutPFmt ! Format to print list of selected output channels to summary file + CHARACTER(30) :: OutPFmtS ! Format to print list of selected output channels to summary file CHARACTER(100) :: Msg ! temporary string for writing appropriate text to summary file + CHARACTER(ChanLen),PARAMETER :: TitleStr(2) = (/ 'Parameter', 'Units ' /) + CHARACTER(ChanLen),PARAMETER :: TitleStrLines(2) = (/ '---------------', '---------------' /) + ! Open the summary file and give it a heading. CALL GetNewUnit( UnSu, ErrStat, ErrMsg ) @@ -2845,7 +2806,7 @@ SUBROUTINE AD_PrintSum( InputFileData, p, p_AD, u, y, ErrStat, ErrMsg ) Msg = 'Stieg Oye dynamic stall model' case (UA_BV) Msg = 'Boeing-Vertol dynamic stall model (e.g. used in CACTUS)' - case default + case default Msg = 'unknown' end select WRITE (UnSu,Ec_IntFrmt) InputFileData%UAMod, 'UAMod', 'Unsteady Aero Model: '//TRIM(Msg) @@ -2887,11 +2848,12 @@ SUBROUTINE AD_PrintSum( InputFileData, p, p_AD, u, y, ErrStat, ErrMsg ) end if - OutPFmt = '( 15x, I4, 2X, A '//TRIM(Num2LStr(ChanLen))//',1 X, A'//TRIM(Num2LStr(ChanLen))//' )' + OutPFmt = '( 15x, I4, 3X,A '//TRIM(Num2LStr(ChanLen))//',1 X, A'//TRIM(Num2LStr(ChanLen))//' )' + OutPFmtS = '( 15x, A4, 3X,A '//TRIM(Num2LStr(ChanLen))//',1 X, A'//TRIM(Num2LStr(ChanLen))//' )' + WRITE (UnSu,'(15x,A)') WRITE (UnSu,'(15x,A)') 'Requested Output Channels:' - WRITE (UnSu,'(15x,A)') 'Col Parameter Units' - WRITE (UnSu,'(15x,A)') '---- -------------- -----' - + WRITE (UnSu,OutPFmtS ) "Col", TitleStr + WRITE (UnSu,OutPFmtS ) "---", TitleStrLines DO I = 0,p%NumOuts WRITE (UnSu,OutPFmt) I, p%OutParam(I)%Name, p%OutParam(I)%Units END DO @@ -2899,11 +2861,12 @@ SUBROUTINE AD_PrintSum( InputFileData, p, p_AD, u, y, ErrStat, ErrMsg ) WRITE (UnSu,'(15x,A)') WRITE (UnSu,'(15x,A)') WRITE (UnSu,'(15x,A)') 'Requested Output Channels at each blade station:' - WRITE (UnSu,'(15x,A)') 'Col Parameter Units' - WRITE (UnSu,'(15x,A)') '---- -------------- -----' + WRITE (UnSu,OutPFmtS ) "Col", TitleStr + WRITE (UnSu,OutPFmtS ) "---", TitleStrLines DO I = 1,p%BldNd_NumOuts WRITE (UnSu,OutPFmt) I, p%BldNd_OutParam(I)%Name, p%BldNd_OutParam(I)%Units - END DO + END DO + CLOSE(UnSu) @@ -3445,7 +3408,7 @@ SUBROUTINE SetOutParam(OutList, p, p_AD, ErrStat, ErrMsg ) ! BNClrnc is set only when we're computing the tower influence do i = 1,size(BNClrnc,2) ! all blades (need to do this in a loop because we need the index of InvalidOutput to be an array of rank one) - InvalidOutput( BNClrnc(:,i) ) = .true. + InvalidOutput( BNClrnc(:,i) ) = .true. end do end if diff --git a/modules/aerodyn/src/BEMT_Registry.txt b/modules/aerodyn/src/BEMT_Registry.txt index 73b2d2cd4d..e3b28be8f5 100644 --- a/modules/aerodyn/src/BEMT_Registry.txt +++ b/modules/aerodyn/src/BEMT_Registry.txt @@ -47,6 +47,7 @@ typedef ^ ^ ReKi typedef ^ ^ ReKi zLocal {:}{:} - - "Distance to blade node, measured along the blade" m typedef ^ ^ ReKi zTip {:} - - "Distance to blade tip, measured along the blade" m typedef ^ ^ ReKi rLocal {:}{:} - - "Radial distance to blade node from the center of rotation, measured in the rotor plane, needed for DBEMT" m +typedef ^ ^ ReKi rTipFix {:} - - "Nominally the coned rotor diameter (without prebend)" m typedef ^ ^ INTEGER UAMod - - - "Model for the dynamic stall equations [1 = Leishman/Beddoes, 2 = Gonzalez, 3 = Minnema]" - typedef ^ ^ LOGICAL UA_Flag - - - "logical flag indicating whether to use UnsteadyAero" - typedef ^ ^ LOGICAL Flookup - - - "Use table lookup for f' and f'' " - diff --git a/modules/aerodyn/src/BEMT_Types.f90 b/modules/aerodyn/src/BEMT_Types.f90 index d7dd816ca9..bb5752eff9 100644 --- a/modules/aerodyn/src/BEMT_Types.f90 +++ b/modules/aerodyn/src/BEMT_Types.f90 @@ -57,6 +57,7 @@ MODULE BEMT_Types REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: zLocal !< Distance to blade node, measured along the blade [m] REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: zTip !< Distance to blade tip, measured along the blade [m] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: rLocal !< Radial distance to blade node from the center of rotation, measured in the rotor plane, needed for DBEMT [m] + REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: rTipFix !< Nominally the coned rotor diameter (without prebend) [m] INTEGER(IntKi) :: UAMod !< Model for the dynamic stall equations [1 = Leishman/Beddoes, 2 = Gonzalez, 3 = Minnema] [-] LOGICAL :: UA_Flag !< logical flag indicating whether to use UnsteadyAero [-] LOGICAL :: Flookup !< Use table lookup for f' and f'' [-] @@ -296,6 +297,18 @@ SUBROUTINE BEMT_CopyInitInput( SrcInitInputData, DstInitInputData, CtrlCode, Err END IF END IF DstInitInputData%rLocal = SrcInitInputData%rLocal +ENDIF +IF (ALLOCATED(SrcInitInputData%rTipFix)) THEN + i1_l = LBOUND(SrcInitInputData%rTipFix,1) + i1_u = UBOUND(SrcInitInputData%rTipFix,1) + IF (.NOT. ALLOCATED(DstInitInputData%rTipFix)) THEN + ALLOCATE(DstInitInputData%rTipFix(i1_l:i1_u),STAT=ErrStat2) + IF (ErrStat2 /= 0) THEN + CALL SetErrStat(ErrID_Fatal, 'Error allocating DstInitInputData%rTipFix.', ErrStat, ErrMsg,RoutineName) + RETURN + END IF + END IF + DstInitInputData%rTipFix = SrcInitInputData%rTipFix ENDIF DstInitInputData%UAMod = SrcInitInputData%UAMod DstInitInputData%UA_Flag = SrcInitInputData%UA_Flag @@ -359,6 +372,9 @@ SUBROUTINE BEMT_DestroyInitInput( InitInputData, ErrStat, ErrMsg ) IF (ALLOCATED(InitInputData%rLocal)) THEN DEALLOCATE(InitInputData%rLocal) ENDIF +IF (ALLOCATED(InitInputData%rTipFix)) THEN + DEALLOCATE(InitInputData%rTipFix) +ENDIF IF (ALLOCATED(InitInputData%UAOff_innerNode)) THEN DEALLOCATE(InitInputData%UAOff_innerNode) ENDIF @@ -445,6 +461,11 @@ SUBROUTINE BEMT_PackInitInput( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrM IF ( ALLOCATED(InData%rLocal) ) THEN Int_BufSz = Int_BufSz + 2*2 ! rLocal upper/lower bounds for each dimension Re_BufSz = Re_BufSz + SIZE(InData%rLocal) ! rLocal + END IF + Int_BufSz = Int_BufSz + 1 ! rTipFix allocated yes/no + IF ( ALLOCATED(InData%rTipFix) ) THEN + Int_BufSz = Int_BufSz + 2*1 ! rTipFix upper/lower bounds for each dimension + Re_BufSz = Re_BufSz + SIZE(InData%rTipFix) ! rTipFix END IF Int_BufSz = Int_BufSz + 1 ! UAMod Int_BufSz = Int_BufSz + 1 ! UA_Flag @@ -629,6 +650,21 @@ SUBROUTINE BEMT_PackInitInput( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrM Re_Xferred = Re_Xferred + 1 END DO END DO + END IF + IF ( .NOT. ALLOCATED(InData%rTipFix) ) THEN + IntKiBuf( Int_Xferred ) = 0 + Int_Xferred = Int_Xferred + 1 + ELSE + IntKiBuf( Int_Xferred ) = 1 + Int_Xferred = Int_Xferred + 1 + IntKiBuf( Int_Xferred ) = LBOUND(InData%rTipFix,1) + IntKiBuf( Int_Xferred + 1) = UBOUND(InData%rTipFix,1) + Int_Xferred = Int_Xferred + 2 + + DO i1 = LBOUND(InData%rTipFix,1), UBOUND(InData%rTipFix,1) + ReKiBuf(Re_Xferred) = InData%rTipFix(i1) + Re_Xferred = Re_Xferred + 1 + END DO END IF IntKiBuf(Int_Xferred) = InData%UAMod Int_Xferred = Int_Xferred + 1 @@ -866,6 +902,24 @@ SUBROUTINE BEMT_UnPackInitInput( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrStat, E Re_Xferred = Re_Xferred + 1 END DO END DO + END IF + IF ( IntKiBuf( Int_Xferred ) == 0 ) THEN ! rTipFix not allocated + Int_Xferred = Int_Xferred + 1 + ELSE + Int_Xferred = Int_Xferred + 1 + i1_l = IntKiBuf( Int_Xferred ) + i1_u = IntKiBuf( Int_Xferred + 1) + Int_Xferred = Int_Xferred + 2 + IF (ALLOCATED(OutData%rTipFix)) DEALLOCATE(OutData%rTipFix) + ALLOCATE(OutData%rTipFix(i1_l:i1_u),STAT=ErrStat2) + IF (ErrStat2 /= 0) THEN + CALL SetErrStat(ErrID_Fatal, 'Error allocating OutData%rTipFix.', ErrStat, ErrMsg,RoutineName) + RETURN + END IF + DO i1 = LBOUND(OutData%rTipFix,1), UBOUND(OutData%rTipFix,1) + OutData%rTipFix(i1) = ReKiBuf(Re_Xferred) + Re_Xferred = Re_Xferred + 1 + END DO END IF OutData%UAMod = IntKiBuf(Int_Xferred) Int_Xferred = Int_Xferred + 1 From 8156e6c2cd80f6b4f9fd28ab1edb2d9c5365b00e Mon Sep 17 00:00:00 2001 From: Bonnie Jonkman Date: Mon, 21 Mar 2022 14:54:48 -0600 Subject: [PATCH 3/3] AD: fix some typos + convert `pi/180` to `D2R` --- modules/aerodyn/src/AeroDyn_AllBldNdOuts_IO.f90 | 2 +- modules/aerodyn/src/AeroDyn_Driver_Subs.f90 | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/modules/aerodyn/src/AeroDyn_AllBldNdOuts_IO.f90 b/modules/aerodyn/src/AeroDyn_AllBldNdOuts_IO.f90 index 3255a0475e..5569b3a641 100644 --- a/modules/aerodyn/src/AeroDyn_AllBldNdOuts_IO.f90 +++ b/modules/aerodyn/src/AeroDyn_AllBldNdOuts_IO.f90 @@ -1253,7 +1253,7 @@ SUBROUTINE AllBldNdOuts_SetParameters( InputFileData, p, p_AD, ErrStat, ErrMsg ) IF ( (InputFileData%BldNd_BladesOut < 0_IntKi) ) then p%BldNd_BladesOut = 0_IntKi ELSE IF ((InputFileData%BldNd_BladesOut > p%NumBlades) ) THEN - CALL SetErrStat( ErrID_Warn, " Number of blades to output data at all blade nodes (BldNd_BladesOut) must be less than "//TRIM(Num2LStr(p%NumBlades))//".", ErrStat, ErrMsg, RoutineName) + CALL SetErrStat( ErrID_Warn, " Number of blades to output data at all blade nodes (BldNd_BladesOut) must be no more than the total number of blades, "//TRIM(Num2LStr(p%NumBlades))//".", ErrStat, ErrMsg, RoutineName) p%BldNd_BladesOut = p%NumBlades ! NOTE: we are forgiving and plateau to numBlades ELSE p%BldNd_BladesOut = InputFileData%BldNd_BladesOut diff --git a/modules/aerodyn/src/AeroDyn_Driver_Subs.f90 b/modules/aerodyn/src/AeroDyn_Driver_Subs.f90 index ce50e3f811..163e4a7cba 100644 --- a/modules/aerodyn/src/AeroDyn_Driver_Subs.f90 +++ b/modules/aerodyn/src/AeroDyn_Driver_Subs.f90 @@ -1392,8 +1392,8 @@ subroutine Dvr_ReadInputFile(fileName, dvr, errStat, errMsg ) call ParseAry(FileInfo_In, CurLine, 'nacOrigin_t'//sWT , wt%nac%origin_t, 3 , errStat2, errMsg2, unEc); if(Failed()) return call ParseAry(FileInfo_In, CurLine, 'hubOrigin_n'//sWT , wt%hub%origin_n, 3 , errStat2, errMsg2, unEc); if(Failed()) return call ParseAry(FileInfo_In, CurLine, 'hubOrientation_n'//sWT , wt%hub%orientation_n, 3 , errStat2, errMsg2, unEc); if(Failed()) return - wt%hub%orientation_n = wt%hub%orientation_n*Pi/180_ReKi - wt%orientationInit = wt%orientationInit*Pi/180_ReKi + wt%hub%orientation_n = wt%hub%orientation_n*D2R + wt%orientationInit = wt%orientationInit*D2R ! Blades call ParseCom(FileInfo_In, CurLine, Line, errStat2, errMsg2, unEc); if(Failed()) return call ParseVar(FileInfo_In, CurLine, 'numBlades'//sWT , wt%numBlades, errStat2, errMsg2, unEc); if(Failed()) return @@ -1837,7 +1837,7 @@ subroutine Dvr_InitializeDriverOutputs(dvr, errStat, errMsg) dvr%out%WriteOutputHdr(j) = 'ShearExp' if (dvr%CompInflow==1) then - dvr%out%WriteOutputUnt(j) = '(NVALID)'; j=j+1 + dvr%out%WriteOutputUnt(j) = '(INVALID)'; j=j+1 else dvr%out%WriteOutputUnt(j) = '(-)'; j=j+1 endif