AeroDyn interface problem

Wind & Water Computer-Aided Engineering Software Tools
1

Hi,

I’m planning to couple AeroDyn with a FEM code as part of my MSc thesis research. I thought it was a good idea to write a simple piece of code to get to know AeroDyn’s interface. The structure I’ve made is completely static, it just calculates the aerodynamic loads on the structure. I initially thought it was working fine but am starting to have some doubts. I have attached the files to this post, I apologize for the sloppy fortran code, it’s the first piece of code I’ve ever written outside of matlab.

There are two weird observations I’ve made:

1.) When using steady wind flow going back to an earlier time gives the result of the previous timestep e.g.

Time is 0.0000000E+00 Force components 30.05191 1.2289569E-02 0.0000000E+00 Moment components 0.0000000E+00 0.0000000E+00 -12.46170 Time is 1.000000 Force components 30.33880 1.2406891E-02 0.0000000E+00 Moment components 0.0000000E+00 0.0000000E+00 -12.58066 Time is 2.000000 Force components 30.62705 1.2524771E-02 0.0000000E+00 Moment components 0.0000000E+00 0.0000000E+00 -12.70019 Time is 3.000000 Force components 30.91667 1.2643209E-02 0.0000000E+00 Moment components 0.0000000E+00 0.0000000E+00 -12.82029 Time is 4.000000 Force components 31.20765 1.2762201E-02 0.0000000E+00 Moment components 0.0000000E+00 0.0000000E+00 -12.94095 Time is 5.000000 Force components 31.49998 1.2881753E-02 0.0000000E+00 Moment components 0.0000000E+00 0.0000000E+00 -13.06217 Time is 0.0000000E+00 Force components 31.49998 1.2881753E-02 0.0000000E+00 Moment components 0.0000000E+00 0.0000000E+00 -13.06217
I can imagine that currtime has to be increasing, but I don’t think I’ve read this in the documentation.

2.) I’ve noticed that when using a turbulent full field wind model only the first time step works correctly, other time steps give NaN errors. I initially thought it might be an allocation problem but using a different variable to store the entries gave the same error. I am not quite sure what the origin of the problem is, I don’t think it’s related to the wind file because the first entry works ok no matter what the currtime value is for the entry.

Time is 0.0000000E+00 Force components 375.8227 0.1050251 0.0000000E+00 Moment components 0.0000000E+00 0.0000000E+00 -155.8176 Time is 1.000000 Force components NaN NaN 0.0000000E+00 Moment components 0.0000000E+00 0.0000000E+00 NaN Time is 2.000000 Force components NaN NaN 0.0000000E+00 Moment components 0.0000000E+00 0.0000000E+00 NaN Time is 3.000000 Force components NaN NaN 0.0000000E+00 Moment components 0.0000000E+00 0.0000000E+00 NaN Time is 4.000000 Force components NaN NaN 0.0000000E+00 Moment components 0.0000000E+00 0.0000000E+00 NaN Time is 5.000000 Force components NaN NaN 0.0000000E+00 Moment components 0.0000000E+00 0.0000000E+00 NaN Time is 0.0000000E+00 Force components NaN NaN 0.0000000E+00 Moment components 0.0000000E+00 0.0000000E+00 NaN

I’d appreciate any help on the matter. I am not quite sure how to proceed.

2

For some reason I can’t attach files to my posts. here’s the fortran code:

[code]program Aeroding

use AeroDyn
use sharedtypes

implicit none

INTEGER :: ErrStat, i, j
REAL(ReKi) :: time, Rotatebottom(3,3), test(3), time2
INTEGER :: Numbl, hubradius, nelem
LOGICAL :: OutputPlottingInfo

type(AeroConfig) :: ADInterfaceComponents
type(AD_InitOptions) :: ADOptions
type(AllAeroMarkers) :: ADAeroMarkers
type(AllAeroLoads) :: ADAeroLoads
type(AllAeroLoads) :: ADAeroLoads2
type(AeroLoadsOptions) :: ALOptions

CALL NWTC_Init(‘Aerotest’, ‘0.01’)
Numbl = 2
Hubradius=1
Nelem= 3

OutputPlottingInfo = .true.

ADOptions%ADInputFile = ‘aerodyn2.ipt’
ADOptions%OutRootName = ‘./out’
ADOptions%WrSumFile = .TRUE.

ADInterfaceComponents%Hub%Position(:slight_smile: = 0.0
ADInterfaceComponents%Hub%Position(3) = 55

ADInterfaceComponents%Hub%Orientation(1,1) = 0.0
ADInterfaceComponents%Hub%Orientation(1,1) = 1.0
ADInterfaceComponents%Hub%Orientation(2,2) = 1.0
ADInterfaceComponents%Hub%Orientation(3,3) = 1.0

IF (.NOT. ALLOCATED( ADInterfaceComponents%Blade ) ) THEN
ALLOCATE( ADInterfaceComponents%Blade( NumBl ), STAT = ErrStat )
IF ( ErrStat /= 0 ) THEN
CALL ProgAbort( ’ Error allocating space for ADInterfaceComponents%Blade.’ )
END IF
END IF

IF (.NOT. ALLOCATED( ADAeroLoads%Blade ) ) THEN
ALLOCATE( ADAeroLoads%Blade(Nelem, NumBl ), STAT = ErrStat )
IF ( ErrStat /= 0 ) THEN
CALL ProgAbort( ’ Error allocating space for ADAeroLoads%Blade.’ )
END IF
END IF

IF (.NOT. ALLOCATED( ADAeroLoads2%Blade ) ) THEN
ALLOCATE( ADAeroLoads2%Blade(Nelem, NumBl ), STAT = ErrStat )
IF ( ErrStat /= 0 ) THEN
CALL ProgAbort( ’ Error allocating space for ADAeroLoads2%Blade.’ )
END IF
END IF

ADInterfaceComponents%Blade(1)%Position(1) = 0.0
ADInterfaceComponents%Blade(1)%Position(2) = 0.0
ADInterfaceComponents%Blade(1)%Position(3) = 56

ADInterfaceComponents%Blade(1)%TranslationVel(:slight_smile: = 0.0
ADInterfaceComponents%Blade(1)%RotationVel(:slight_smile: = 0.0
ADInterfaceComponents%Blade(1)%Orientation(1,1) = 0
ADInterfaceComponents%Blade(1)%Orientation(2,1) = 0.0
ADInterfaceComponents%Blade(1)%Orientation(3,1) = 1.0

ADInterfaceComponents%Blade(1)%Orientation(1,2) = 0.0
ADInterfaceComponents%Blade(1)%Orientation(2,2) = 1.0
ADInterfaceComponents%Blade(1)%Orientation(3,2) = 0.0

ADInterfaceComponents%Blade(1)%Orientation(1,3) = 0
ADInterfaceComponents%Blade(1)%Orientation(2,3) = 0.0
ADInterfaceComponents%Blade(1)%Orientation(3,3) = 1

ADInterfaceComponents%Blade(2)%Position(1) = 0
ADInterfaceComponents%Blade(2)%Position(2) = 0.0
ADInterfaceComponents%Blade(2)%Position(3) = 54

ADInterfaceComponents%Blade(2)%TranslationVel(:slight_smile: = 0.0
ADInterfaceComponents%Blade(2)%RotationVel(:slight_smile: = 0.0
ADInterfaceComponents%Blade(2)%Orientation(1,1) = 0
ADInterfaceComponents%Blade(2)%Orientation(2,1) = 0.0
ADInterfaceComponents%Blade(2)%Orientation(3,1) = 1.0

ADInterfaceComponents%Blade(2)%Orientation(1,2) = 0.0
ADInterfaceComponents%Blade(2)%Orientation(2,2) = -1.0
ADInterfaceComponents%Blade(2)%Orientation(3,2) = 0.0

ADInterfaceComponents%Blade(2)%Orientation(1,3) = 0
ADInterfaceComponents%Blade(2)%Orientation(2,3) = 0.0
ADInterfaceComponents%Blade(2)%Orientation(3,3) = -1

ADInterfaceComponents%BladeLength = 15

ADInterfaceComponents%RotorFurl%Position(1) = 0.0
ADInterfaceComponents%RotorFurl%Position(2) = 0.0
ADInterfaceComponents%RotorFurl%Position(3) = 54.5

ADInterfaceComponents%RotorFurl%Orientation(1,1) = 1
ADInterfaceComponents%RotorFurl%Orientation(2,1) = 0.0
ADInterfaceComponents%RotorFurl%Orientation(3,1) = 0.0

ADInterfaceComponents%RotorFurl%Orientation(1,2) = 0.0
ADInterfaceComponents%RotorFurl%Orientation(2,2) = 1.0
ADInterfaceComponents%RotorFurl%Orientation(3,2) = 0.0

ADInterfaceComponents%RotorFurl%Orientation(1,3) = 0.0
ADInterfaceComponents%RotorFurl%Orientation(2,3) = 0.0
ADInterfaceComponents%RotorFurl%Orientation(3,3) = 1

ADInterfaceComponents%RotorFurl%TranslationVel(:slight_smile: = 0.0
ADInterfaceComponents%RotorFurl%RotationVel(:slight_smile: = 0.0

ADInterfaceComponents%Nacelle%Position(1) = -1
ADInterfaceComponents%Nacelle%Position(2) = 0.0
ADInterfaceComponents%Nacelle%Position(3) = 0

ADInterfaceComponents%Nacelle%Orientation(1,1) = 1
ADInterfaceComponents%Nacelle%Orientation(2,1) = 0.0
ADInterfaceComponents%Nacelle%Orientation(3,1) = 0.0

ADInterfaceComponents%Nacelle%Orientation(1,2) = 0.0
ADInterfaceComponents%Nacelle%Orientation(2,2) = 1.0
ADInterfaceComponents%Nacelle%Orientation(3,2) = 0.0

ADInterfaceComponents%Nacelle%Orientation(1,3) = 0.0
ADInterfaceComponents%Nacelle%Orientation(2,3) = 0.0
ADInterfaceComponents%Nacelle%Orientation(3,3) = 1

ADInterfaceComponents%Nacelle%TranslationVel(:slight_smile: = 0.0
ADInterfaceComponents%Nacelle%RotationVel(:slight_smile: = 0.0

ADInterfaceComponents%Tower%Position(1) = 0
ADInterfaceComponents%Tower%Position(2) = 0.0
ADInterfaceComponents%Tower%Position(3) = 0

ADInterfaceComponents%Tower%Orientation(1,1) = 0
ADInterfaceComponents%Tower%Orientation(2,1) = 0.0
ADInterfaceComponents%Tower%Orientation(3,1) = 1.0

ADInterfaceComponents%Tower%Orientation(1,2) = 0.0
ADInterfaceComponents%Tower%Orientation(2,2) = 1.0
ADInterfaceComponents%Tower%Orientation(3,2) = 0.0

ADInterfaceComponents%Tower%Orientation(1,3) = 1.0
ADInterfaceComponents%Tower%Orientation(2,3) = 0.0
ADInterfaceComponents%Tower%Orientation(3,3) = 0

ADInterfaceComponents%Tower%TranslationVel(:slight_smile: = 0.0
ADInterfaceComponents%Tower%RotationVel(:slight_smile: = 0.0

ALLOCATE(ALOptions%SetMulTabLoc(3,2))
ALLOCATE(ALOptions%MulTabLoc(3,2))
ALOptions%SetMulTabLoc(:,:slight_smile: = .FALSE.
ALOptions%LinearizeFlag = .FALSE.

ErrStat = 0

time = 0

ADAeroMarkers = AD_Init(ADOptions, ADInterfaceComponents, ErrStat)

write(,) 'Hub radius ’ , DOT_PRODUCT(ADInterfaceComponents%Blade(2)%Position(:slight_smile: - ADInterfaceComponents%Hub%Position(:), ADInterfaceComponents%Blade(2)%Orientation(3,:slight_smile: )

!write(,) 'Tip radius ’ , ADInterfaceComponents%BladeLength + HubRadius

!! AD_Init gives BLADE FIXED coordinates of the aerodynamic markers, AD_CaculateLoads needs GLOBAL coordinates to function.

Rotatebottom(1,1)=1
Rotatebottom(1,2)=0
Rotatebottom(1,3)=0
Rotatebottom(2,1)=0
Rotatebottom(2,2)=-1
Rotatebottom(2,3)=0
Rotatebottom(3,1)=0
Rotatebottom(3,2)=0
Rotatebottom(3,3)=-1
!Change orientation on bottom blade
do i=1, nelem
ADAeromarkers%Blade(i,2)%Orientation(:,:)=matmul(ADAeromarkers%Blade(i,2)%Orientation(:,:),RotateBottom(:,:))
enddo

do i=1, Nelem
ADAeromarkers%Blade(i,2)%Position=Matmul(Rotatebottom,ADAeromarkers%Blade(i,2)%Position)
enddo

! Offset position for GLOBAL reference frame
do j=1, 3
ADAeroMarkers%Blade(j,1)%Position(3)=ADAeroMarkers%Blade(j,1)%Position(3)+50+hubradius
ADAeroMarkers%Blade(j,2)%Position(3)=ADAeroMarkers%Blade(j,2)%Position(3)+50-hubradius
enddo

!write(,) ‘Rotorfurl’, ADAeroMarkers%RotorFurl(1)%Position(:slight_smile:
ADAeroMarkers%Hub%Position(3)=ADAeroMarkers%Hub%Position(3)+50
ADAeroMarkers%RotorFurl%Position(3)=ADAeroMarkers%RotorFurl%Position(3)+50
ADAeroMarkers%Nacelle%Position(3)=ADAeroMarkers%Nacelle%Position(3)+50
ADAeroMarkers%Tower%Position(3)=ADAeroMarkers%Tower%Position(3)+50
ADAeroMarkers%Tail%Position(3)=ADAeroMarkers%Tail%Position(3)+50
!write(,) ADAeromarkers%Blade(3,2)%Position

! Read a value if interested.
!do, i=1,3
!write(,) ( Rotatebottom(i,j), j=1,3 )
!enddo

! Read a value if interested.
!DO i = 1, 3
!PRINT*,(ADAeroMarkers%Blade(3,2)%Position(i))
!END DO
!write(,) ADInterfaceComponents%Rotorfurl%RotationVel
!do, i=1,3
!write(,) ( ADAeromarkers%Blade(1,1)%Orientation(i,j), j=1,3 )
!enddo
do time=0,1,0.02
ADAeroLoads = AD_CalculateLoads(time,ADAeroMarkers,ADInterfaceComponents,ALOptions,ErrStat)
write(,) 'Time is ', time
write(,) ‘Force components’, ADAeroLoads%Blade(1,1)%Force(:slight_smile:
write(,) ‘Moment components’, ADAeroLoads%Blade(1,1)%Moment(:slight_smile:
!write(,) 'Orientation of el1 bl1 ', ADAeromarkers%Blade(2,2)%Position(:,:slight_smile:
end do

time2 = 0
ADAeroLoads2 = AD_CalculateLoads(time2,ADAeroMarkers,ADInterfaceComponents,ALOptions,ErrStat)
write(,) ‘Time is’, time2
write(,) ‘Force components’, ADAeroLoads2%Blade(1,1)%Force(:slight_smile:
write(,) ‘Moment components’, ADAeroLoads2%Blade(1,1)%Moment(:slight_smile:

CALL AD_Terminate(ErrStat)

end program Aeroding[/code]

Steady wind input file:

! Steady wind file created 7/22/98 7:58:37 PM by YawDynVB Version 2.0 ! Time Wind Wind Vert. Horiz. Vert. LinV Gust ! Speed Dir Speed Shear Shear Shear Speed 0.0 5 0 0 0 0 0 0 630.0 20 0 0 0 0 0 0
And finally the AeroDyn input file:

Combined Experiment Baseline for YawDyn version 12.5 SI SysUnits - System of units for used for input and output [must be SI for FAST] (unquoted string) STEADY StallMod - Dynamic stall included [BEDDOES or STEADY] (unquoted string) USE_CM UseCm - Use aerodynamic pitching moment model? [USE_CM or NO_CM] (unquoted string) EQUIL InfModel - Inflow model [DYNIN or EQUIL] (unquoted string) SWIRL IndModel - Induction-factor model [NONE or WAKE or SWIRL] (unquoted string) 5.0000E-03 AToler - Induction-factor tolerance (convergence criteria) (-) PRANDtl TLModel - Tip-loss model (EQUIL only) [PRANDtl, GTECH, or NONE] (unquoted string) PRANDtl HLModel - Hub-loss model (EQUIL only) [PRANDtl or NONE] (unquoted string) "TurbSim.wnd" WindFile - Name of file containing wind data (quoted string) 55.0 HH - Wind reference (hub) height [TowerHt+Twr2Shft+OverHang*SIN(ShftTilt)] (m) 0 TwrShad - Tower-shadow velocity deficit (-) 9999.9 ShadHWid - Tower-shadow half width (m) 9999.9 T_Shad_Refpt - Tower-shadow reference point (m) 1.225 AirDens - Air density (kg/m^3) 1.625e-5 KinVisc - Kinematic air viscosity (m^2/sec) 0.025 DTAero - Time interval for aerodynamic calculations (sec) 1 NumFoil - Number of airfoil files (-) "S809_cln.dat" FoilNm - Names of the airfoil files [NumFoil lines] (quoted strings) 3 BldNodes - Number of blade nodes used for analysis (-) RELM Twist DR Chord NFoil ElPrList { Twist ignored by ADAMS (but placeholders must be present) } 3.5 5.0000 5 1.5000 1 PRINT 8.5 5.0000 5 1.5000 1 PRINT 13.5 5.0000 5 1.5000 1 PRINT

3

I also noticed that the output file generated by AeroDyn is corrupted?

out.opt

This file was generated by AeroDyn (v13.00.01a-bjj, 16-Feb-2012) on 15-Jan-2013 at 01:52:28.

Inputs read in from the AeroDyn input file:
Combined Experiment Baseline for YawDyn version 12.5

SI	Units for input and output         
STEADY	Dynamic stall model                 [NO Dynamic stall]
USE_CM	Aerodynamic pitching moment model   [Pitching Moments calculated]
EQUIL	Inflow model                        [Equilibrium]
SWIRL	Induction factor model              [Normal and Radial flow induction factors calculated]
5.00000E-03	Convergence tolerance for induction factor
PRAND	Tip-loss model                      [Prandtl model]
PRAND	Hub-loss model                      [Prandtl model]
"./TurbSim.wnd"	  Wind file name
55	Wind reference (hub) height, m
0	Tower shadow centerline velocity deficit
9999.9	Tower shadow half width, m
9999.9	Tower shadow reference point, m
1.225	Air density, kg/m^3
1.62500E-05	Kinematic air viscosity, m^2/sec
2.50000E-02	Time interval for aerodynamic calculations, sec
1	Number of airfoil files used. Files listed below:
"./S809_cln.dat"
3	Number of blade elements per blade

     Element     RELM       Twist       DR        Chord      NFoil     Print?    Tip-loss   Hub-loss 
       (-)        (m)       (deg)       (m)        (m)        (-)     (Yes/No)   constant   constant 
   ---------- ---------- ---------- ---------- ---------- ---------- ---------- ---------- ----------
            1    3.50000    5.00000    5.00000    1.50000          1        Yes    3.57143    0.00000
            2    8.50000    5.00000    5.00000    1.50000          1        Yes    0.88235    0.00000
            3   13.50000    5.00000    5.00000    1.50000          1        Yes    0.18519    0.00000


 Rotor radius     =  -0.000 m
 Hub radius       =  -0.000 m
 Number of blades =   2

Blade element aerodynamic time series data written to file.

I am not sure why the rotor and hub radius are at 0, I double checked it even added a piece of code to calculate it to my fortran file (which gave the correct answer)

out.elm

This file was generated by \00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00\00 on 15-Jan-2013 at 01:52:28. Time VX03 VY03 VZ03 Alpha01 DynPres01 CLift01 CDrag01 CNorm01 CTang01 CMomt01 Pitch01 AxInd01 TanInd01 ForcN01 ForcT01 Pmomt01 ReNum01 Alpha02 DynPres02 CLift02 CDrag02 CNorm02 CTang02 CMomt02 Pitch02 AxInd02 TanInd02 ForcN02 ForcT02 Pmomt02 ReNum02 Alpha03 DynPres03 CLift03 CDrag03 CNorm03 CTang03 CMomt03 Pitch03 AxInd03 TanInd03 ForcN03 ForcT03 Pmomt03 ReNum03 (sec) (m/sec) (m/sec) (m/sec) (deg) (Pa) (-) (-) (-) (-) (-) (deg) (-) (-) (N) (N) (N-m) (x10^6) (deg) (Pa) (-) (-) (-) (-) (-) (deg) (-) (-) (N) (N) (N-m) (x10^6) (deg) (Pa) (-) (-) (-) (-) (-) (deg) (-) (-) (N) (N) (N-m) (x10^6)

The \00\00… part is corrupted I believe?

4

Hi,

I’m not sure I’m the best person to answer your questions, but here are my two cents:

AeroDyn has some logic in it (see the first few lines of AD_CalculateLoads) so that if it hasn’t advanced at least DTAero seconds, it returns the previously calculated values.

I am also pretty sure that you should be calling it at least every DTAero seconds (don’t skip steps), and that it needs to calculate loads at least 4 times before Time = 1 second. (The integration method needs to be initialized with at least 4 time steps, but it checks that time > 1 second instead of checking that the number of steps > 4.)

These are things that will be fixed once we convert to the new modularization framework, but we aren’t there, yet…

The turbulent wind will produce errors if you’ve asked for wind outside its bounds (i.e., the FF grid isn’t big enough or long enough), whereas the HH wind files will always return some value. When the turbulent wind routines return errors, it generally indicate a problem computing positions elsewhere in AeroDyn.

As for the file corruption part, perhaps the ProgName and ProgVer variables from the NWTC Library haven’t been initialized? I would also recommend that you compile AeroDyn with the /QSave and /QZero compiler options.

5

Hi Bonnie,

Thanks for your input.

I think that might explain the result of the steady wind problem.

I set TurbSim to a grid 80 m wide, 80 m high and a hub height of 55 m. The fictitious structure has a hub height of 55 m and blades of 15 m, so should stay within the 80 m height. Here’s the input I used with TurbSim:

[code]TurbSim Input File. Valid for TurbSim v1.06.00, 21-Sep-2012

---------Runtime Options-----------------------------------
2318573 RandSeed1 - First random seed (-2147483648 to 2147483647)
RANLUX RandSeed2 - Second random seed (-2147483648 to 2147483647) for intrinsic pRNG, or an alternative pRNG: “RanLux” or “RNSNLW”
False WrBHHTP - Output hub-height turbulence parameters in binary form? (Generates RootName.bin)
False WrFHHTP - Output hub-height turbulence parameters in formatted form? (Generates RootName.dat)
False WrADHH - Output hub-height time-series data in AeroDyn form? (Generates RootName.hh)
True WrADFF - Output full-field time-series data in TurbSim/AeroDyn form? (Generates Rootname.bts)
True WrBLFF - Output full-field time-series data in BLADED/AeroDyn form? (Generates RootName.wnd)
False WrADTWR - Output tower time-series data? (Generates RootName.twr)
False WrFMTFF - Output full-field time-series data in formatted (readable) form? (Generates RootName.u, RootName.v, RootName.w)
True WrACT - Output coherent turbulence time steps in AeroDyn form? (Generates RootName.cts)
True Clockwise - Clockwise rotation looking downwind? (used only for full-field binary files - not necessary for AeroDyn)
0 ScaleIEC - Scale IEC turbulence models to exact target standard deviation? [0=no additional scaling; 1=use hub scale uniformly; 2=use individual scales]

--------Turbine/Model Specifications-----------------------
13 NumGrid_Z - Vertical grid-point matrix dimension
13 NumGrid_Y - Horizontal grid-point matrix dimension
0.05 TimeStep - Time step [seconds]
1000 AnalysisTime - Length of analysis time series [seconds] (program will add time if necessary: AnalysisTime = MAX(AnalysisTime, UsableTime+GridWidth/MeanHHWS) )
600 UsableTime - Usable length of output time series [seconds] (program will add GridWidth/MeanHHWS seconds)
55.00 HubHt - Hub height [m] (should be > 0.5GridHeight)
80.00 GridHeight - Grid height [m]
80.00 GridWidth - Grid width [m] (should be >= 2(RotorRadius+ShaftLength))
0 VFlowAng - Vertical mean flow (uptilt) angle [degrees]
0 HFlowAng - Horizontal mean flow (skew) angle [degrees]

--------Meteorological Boundary Conditions-------------------
“IECVKM” TurbModel - Turbulence model (“IECKAI”=Kaimal, “IECVKM”=von Karman, “GP_LLJ”, “NWTCUP”, “SMOOTH”, “WF_UPW”, “WF_07D”, “WF_14D”, “TIDAL”, or “NONE”)
“1-ED3” IECstandard - Number of IEC 61400-x standard (x=1,2, or 3 with optional 61400-1 edition number (i.e. “1-Ed2”) )
“A” IECturbc - IEC turbulence characteristic (“A”, “B”, “C” or the turbulence intensity in percent) (“KHTEST” option with NWTCUP model, not used for other models)
“NTM” IEC_WindType - IEC turbulence type (“NTM”=normal, “xETM”=extreme turbulence, “xEWM1”=extreme 1-year wind, “xEWM50”=extreme 50-year wind, where x=wind turbine class 1, 2, or 3)
default ETMc - IEC Extreme Turbulence Model “c” parameter [m/s]
default WindProfileType - Wind profile type (“JET”;“LOG”=logarithmic;“PL”=power law;“H2L”=Log law for TIDAL spectral model;“IEC”=PL on rotor disk, LOG elsewhere; or “default”)
55 RefHt - Height of the reference wind speed [m]
10 URef - Mean (total) wind speed at the reference height [m/s] (or “default” for JET wind profile)
default ZJetMax - Jet height [m] (used only for JET wind profile, valid 70-490 m)
default PLExp - Power law exponent [-] (or “default”)
default Z0 - Surface roughness length [m] (or “default”)

--------Non-IEC Meteorological Boundary Conditions------------
default Latitude - Site latitude [degrees] (or “default”)
0.05 RICH_NO - Gradient Richardson number
default UStar - Friction or shear velocity [m/s] (or “default”)
default ZI - Mixing layer depth [m] (or “default”)
default PC_UW - Hub mean u’w’ Reynolds stress (or “default”)
default PC_UV - Hub mean u’v’ Reynolds stress (or “default”)
default PC_VW - Hub mean v’w’ Reynolds stress (or “default”)
default IncDec1 - u-component coherence parameters (e.g. “10.0 0.3e-3” in quotes) (or “default”)
default IncDec2 - v-component coherence parameters (e.g. “10.0 0.3e-3” in quotes) (or “default”)
default IncDec3 - w-component coherence parameters (e.g. “10.0 0.3e-3” in quotes) (or “default”)
default CohExp - Coherence exponent (or “default”)

--------Coherent Turbulence Scaling Parameters-------------------
“c:\aero\turbsim\EventData” CTEventPath - Name of the path where event data files are located
“Random” CTEventFile - Type of event files (“LES”, “DNS”, or “RANDOM”)
true Randomize - Randomize the disturbance scale and locations? (true/false)
1.0 DistScl - Disturbance scale (ratio of wave height to rotor disk). (Ignored when Randomize = true.)
0.5 CTLy - Fractional location of tower centerline from right (looking downwind) to left side of the dataset. (Ignored when Randomize = true.)
0.5 CTLz - Fractional location of hub height from the bottom of the dataset. (Ignored when Randomize = true.)
30.0 CTStartTime - Minimum start time for coherent structures in RootName.cts [seconds]

==================================================
NOTE: Do not add or remove any lines in this file!
==================================================[/code]

I tried using smaller steps, but even with a dt=0.05 the second calculation still gets NaN errors.

[code]Using NWTC Subroutine Library (v1.04.01, 21-Feb-2012).

Aerodynamic loads calculated using AeroDyn (v13.00.01a-bjj, 16-Feb-2012).
Heading of the AeroDyn input file: Combined Experiment Baseline for YawDyn version 12.5
Using NWTC Subroutine Library (v1.04.01, 21-Feb-2012).

Using InflowWind (v1.00.01b-bjj, 14-Nov-2011)

Assuming ./TurbSim.wnd is a binary FF wind file.
Using NWTC Subroutine Library (v1.04.01, 21-Feb-2012).

Reading a 13x13 grid (80 m wide, 15 m to 95 m above ground) with a characterstic wind speed of
15 m/s.
Hub radius 1.000000 of 20-Hz full-field data (605.3 seconds).
Time is 0.0000000E+00
Force components 375.8227 0.1050251 0.0000000E+00
Moment components 0.0000000E+00 0.0000000E+00 -155.8176
Time is 5.0000001E-02
Force components NaN NaN 0.0000000E+00
Moment components 0.0000000E+00 0.0000000E+00 NaN
Time is 0.1000000
Force components NaN NaN 0.0000000E+00
Moment components 0.0000000E+00 0.0000000E+00 NaN
Time is 0.1500000
Force components NaN NaN 0.0000000E+00
Moment components 0.0000000E+00 0.0000000E+00 NaN
Time is 0.2000000
Force components NaN NaN 0.0000000E+00
Moment components 0.0000000E+00 0.0000000E+00 NaN
[/code]

What I think is interesting is that no matter what the currtime is set to at the first calculation, it always returns an answer. e.g.

Reading a 13x13 grid (80 m wide, 15 m to 95 m above ground) with a characterstic wind speed of 15 m/s. Hub radius 1.000000 of 20-Hz full-field data (605.3 seconds). Time is 10.00000 Force components 221.2051 0.9392253 0.0000000E+00 Moment components 0.0000000E+00 0.0000000E+00 -91.68458 Time is 10.05000 Force components NaN NaN 0.0000000E+00 Moment components 0.0000000E+00 0.0000000E+00 NaN

Would the structural code have to initialise NWTC subroutines for AeroDyn output files to work (I thought AeroDyn did all of that internally).

6

Hi Sander,

AeroDyn does initialize the library internally, but the library’s global variables ProgName and ProgVer (which will be removed from future versions) are set by the structural code and AeroDyn prints them in its output file. I don’t know why it didn’t just print spaces instead of the \00 characters, but that is something you can check.

I would always start AeroDyn simulations at time = 0. AeroDyn tends to be pretty “unforgiving” in the sense that if you give it positions/orientations that aren’t realistic, it blows up quickly. I’m not an aerodynamic expert, so I can’t say exactly why that is… You might want to plot the markers you’ve set and see if anything looks strange. The AeroDyn code has a parameter called “OutputPlottingInfo” that you can set to .true. and then it will print all the markers to a file that you can read and plot.

Maybe someone who knows something about aerodynamics can chime in here… :slight_smile:

7

Hi Bonnie,

Adding a call to NWTC seemed to fix the 00/00 problem. The output file still doesn’t list any results though (I enabled print for all the elements in the AeroDyn input file).

I was not aware of a plotting parameter in AeroDyn, how do I enable it?

Thanks,

Sander

8

Another little quirk in AeroDyn is that the driver program has to call AeroDyn’s routine to write to the element file. If you add CALL ElemOut() after you calculate the aerodynamic loads, you should get some results.

I put an undocumented feature in AeroDyn at one point to help visualize the markers. Line 69 in AeroDyn.f90 says:

LOGICAL, PRIVATE, PARAMETER :: OutputPlottingInfo = .false. !BJJ SET THIS = .FALSE. !!!! To enable the feature, just change it to .true. and recompile.

I think it prints every time step that loads are calculated, so you would never want this enabled during a real simualation (or for more than a couple of timesteps or two). The trick is really in plotting the information. I’m not prepared to support the plotting efforts, but I do have some Matlab scripts that you are welcome to use if that helps you. (email me)

9

I believe I’ve fixed the errors I got in the turbulent flow simulation. I think there was a mistake in one of the orientation matrices.

I’ve included the fortran code in my post, in case someone ever wants to have a look at it.

Thanks for the help,

Sander

[code]program Aeroding

use AeroDyn
use sharedtypes
use aerogensubs

implicit none

INTEGER :: ErrStat, i, j
REAL(ReKi) :: time, Rotatebottom(3,3), test(3), time2
INTEGER :: Numbl, hubradius, nelem
LOGICAL :: OutputPlottingInfo

type(AeroConfig) :: ADInterfaceComponents
type(AD_InitOptions) :: ADOptions
type(AllAeroMarkers) :: ADAeroMarkers
type(AllAeroLoads) :: ADAeroLoads
type(AllAeroLoads) :: ADAeroLoads2
type(AeroLoadsOptions) :: ALOptions

CALL NWTC_Init(‘Aerotest’, ‘0.01’)
Numbl = 2
Hubradius=1
Nelem= 3

OutputPlottingInfo = .true.

ADOptions%ADInputFile = ‘aerodyn2.ipt’
ADOptions%OutRootName = ‘./out’
ADOptions%WrSumFile = .TRUE.

ADInterfaceComponents%Hub%Position(:slight_smile: = 0.0
ADInterfaceComponents%Hub%Position(3) = 55

ADInterfaceComponents%Hub%Orientation(1,1) = 0.0
ADInterfaceComponents%Hub%Orientation(1,1) = 1.0
ADInterfaceComponents%Hub%Orientation(2,2) = 1.0
ADInterfaceComponents%Hub%Orientation(3,3) = 1.0

IF (.NOT. ALLOCATED( ADInterfaceComponents%Blade ) ) THEN
ALLOCATE( ADInterfaceComponents%Blade( NumBl ), STAT = ErrStat )
IF ( ErrStat /= 0 ) THEN
CALL ProgAbort( ’ Error allocating space for ADInterfaceComponents%Blade.’ )
END IF
END IF

IF (.NOT. ALLOCATED( ADAeroLoads%Blade ) ) THEN
ALLOCATE( ADAeroLoads%Blade(Nelem, NumBl ), STAT = ErrStat )
IF ( ErrStat /= 0 ) THEN
CALL ProgAbort( ’ Error allocating space for ADAeroLoads%Blade.’ )
END IF
END IF

IF (.NOT. ALLOCATED( ADAeroLoads2%Blade ) ) THEN
ALLOCATE( ADAeroLoads2%Blade(Nelem, NumBl ), STAT = ErrStat )
IF ( ErrStat /= 0 ) THEN
CALL ProgAbort( ’ Error allocating space for ADAeroLoads2%Blade.’ )
END IF
END IF

ADInterfaceComponents%Blade(1)%Position(1) = 0.0
ADInterfaceComponents%Blade(1)%Position(2) = 0.0
ADInterfaceComponents%Blade(1)%Position(3) = 56

ADInterfaceComponents%Blade(1)%TranslationVel(:slight_smile: = 0.0
ADInterfaceComponents%Blade(1)%RotationVel(:slight_smile: = 0.0
ADInterfaceComponents%Blade(1)%Orientation(1,1) = 1
ADInterfaceComponents%Blade(1)%Orientation(2,1) = 0.0
ADInterfaceComponents%Blade(1)%Orientation(3,1) = 0.0

ADInterfaceComponents%Blade(1)%Orientation(1,2) = 0.0
ADInterfaceComponents%Blade(1)%Orientation(2,2) = 1.0
ADInterfaceComponents%Blade(1)%Orientation(3,2) = 0.0

ADInterfaceComponents%Blade(1)%Orientation(1,3) = 0
ADInterfaceComponents%Blade(1)%Orientation(2,3) = 0.0
ADInterfaceComponents%Blade(1)%Orientation(3,3) = 1

ADInterfaceComponents%Blade(2)%Position(1) = 0
ADInterfaceComponents%Blade(2)%Position(2) = 0.0
ADInterfaceComponents%Blade(2)%Position(3) = 54

ADInterfaceComponents%Blade(2)%TranslationVel(:slight_smile: = 0.0
ADInterfaceComponents%Blade(2)%RotationVel(:slight_smile: = 0.0
ADInterfaceComponents%Blade(2)%Orientation(1,1) = 1
ADInterfaceComponents%Blade(2)%Orientation(2,1) = 0.0
ADInterfaceComponents%Blade(2)%Orientation(3,1) = 0.0

ADInterfaceComponents%Blade(2)%Orientation(1,2) = 0.0
ADInterfaceComponents%Blade(2)%Orientation(2,2) = -1.0
ADInterfaceComponents%Blade(2)%Orientation(3,2) = 0.0

ADInterfaceComponents%Blade(2)%Orientation(1,3) = 0
ADInterfaceComponents%Blade(2)%Orientation(2,3) = 0.0
ADInterfaceComponents%Blade(2)%Orientation(3,3) = -1

ADInterfaceComponents%BladeLength = 15

ADInterfaceComponents%RotorFurl%Position(1) = 0.0
ADInterfaceComponents%RotorFurl%Position(2) = 0.0
ADInterfaceComponents%RotorFurl%Position(3) = 54.5

ADInterfaceComponents%RotorFurl%Orientation(1,1) = 1
ADInterfaceComponents%RotorFurl%Orientation(2,1) = 0.0
ADInterfaceComponents%RotorFurl%Orientation(3,1) = 0.0

ADInterfaceComponents%RotorFurl%Orientation(1,2) = 0.0
ADInterfaceComponents%RotorFurl%Orientation(2,2) = 1.0
ADInterfaceComponents%RotorFurl%Orientation(3,2) = 0.0

ADInterfaceComponents%RotorFurl%Orientation(1,3) = 0.0
ADInterfaceComponents%RotorFurl%Orientation(2,3) = 0.0
ADInterfaceComponents%RotorFurl%Orientation(3,3) = 1

ADInterfaceComponents%RotorFurl%TranslationVel(:slight_smile: = 0.0
ADInterfaceComponents%RotorFurl%RotationVel(:slight_smile: = 0.0

ADInterfaceComponents%Nacelle%Position(1) = 54
ADInterfaceComponents%Nacelle%Position(2) = 0.0
ADInterfaceComponents%Nacelle%Position(3) = 0

ADInterfaceComponents%Nacelle%Orientation(1,1) = 1
ADInterfaceComponents%Nacelle%Orientation(2,1) = 0.0
ADInterfaceComponents%Nacelle%Orientation(3,1) = 0.0

ADInterfaceComponents%Nacelle%Orientation(1,2) = 0.0
ADInterfaceComponents%Nacelle%Orientation(2,2) = 1.0
ADInterfaceComponents%Nacelle%Orientation(3,2) = 0.0

ADInterfaceComponents%Nacelle%Orientation(1,3) = 0.0
ADInterfaceComponents%Nacelle%Orientation(2,3) = 0.0
ADInterfaceComponents%Nacelle%Orientation(3,3) = 1

ADInterfaceComponents%Nacelle%TranslationVel(:slight_smile: = 0.0
ADInterfaceComponents%Nacelle%RotationVel(:slight_smile: = 0.0

ADInterfaceComponents%Tower%Position(1) = 0
ADInterfaceComponents%Tower%Position(2) = 0.0
ADInterfaceComponents%Tower%Position(3) = 0

ADInterfaceComponents%Tower%Orientation(1,1) = 0
ADInterfaceComponents%Tower%Orientation(2,1) = 0.0
ADInterfaceComponents%Tower%Orientation(3,1) = 1.0

ADInterfaceComponents%Tower%Orientation(1,2) = 0.0
ADInterfaceComponents%Tower%Orientation(2,2) = 1.0
ADInterfaceComponents%Tower%Orientation(3,2) = 0.0

ADInterfaceComponents%Tower%Orientation(1,3) = 1.0
ADInterfaceComponents%Tower%Orientation(2,3) = 0.0
ADInterfaceComponents%Tower%Orientation(3,3) = 0

ADInterfaceComponents%Tower%TranslationVel(:slight_smile: = 0.0
ADInterfaceComponents%Tower%RotationVel(:slight_smile: = 0.0

ALLOCATE(ALOptions%SetMulTabLoc(3,2))
ALLOCATE(ALOptions%MulTabLoc(3,2))
ALOptions%SetMulTabLoc(:,:slight_smile: = .FALSE.
ALOptions%LinearizeFlag = .FALSE.

ErrStat = 0

time = 0

write(,) ‘Hub radius blade 1’ , DOT_PRODUCT(ADInterfaceComponents%Blade(1)%Position(:slight_smile: - ADInterfaceComponents%Hub%Position(:), ADInterfaceComponents%Blade(1)%Orientation(3,:slight_smile: )
write(,) ‘Orientation’, ADInterfaceComponents%Blade(1)%Orientation(3,:slight_smile:
ADAeroMarkers = AD_Init(ADOptions, ADInterfaceComponents, ErrStat)

write(,) ‘Hub radius blade 2’ , DOT_PRODUCT(ADInterfaceComponents%Blade(2)%Position(:slight_smile: - ADInterfaceComponents%Hub%Position(:), ADInterfaceComponents%Blade(2)%Orientation(3,:slight_smile: )

!write(,) 'Tip radius ’ , ADInterfaceComponents%BladeLength + HubRadius

!! AD_Init gives BLADE FIXED coordinates of the aerodynamic markers, AD_CaculateLoads needs GLOBAL coordinates to function.

Rotatebottom(1,1)=1
Rotatebottom(1,2)=0
Rotatebottom(1,3)=0
Rotatebottom(2,1)=0
Rotatebottom(2,2)=-1
Rotatebottom(2,3)=0
Rotatebottom(3,1)=0
Rotatebottom(3,2)=0
Rotatebottom(3,3)=-1
!Change orientation on bottom blade
do i=1, nelem
ADAeromarkers%Blade(i,2)%Orientation(:,:)=matmul(ADAeromarkers%Blade(i,2)%Orientation(:,:),RotateBottom(:,:))
enddo

do i=1, Nelem
ADAeromarkers%Blade(i,2)%Position=Matmul(Rotatebottom,ADAeromarkers%Blade(i,2)%Position)
enddo

! Offset position for GLOBAL reference frame
do j=1, 3
ADAeroMarkers%Blade(j,1)%Position(3)=ADAeroMarkers%Blade(j,1)%Position(3)+50+hubradius
ADAeroMarkers%Blade(j,2)%Position(3)=ADAeroMarkers%Blade(j,2)%Position(3)+50-hubradius
enddo

!write(,) ‘Rotorfurl’, ADAeroMarkers%RotorFurl(1)%Position(:slight_smile:
ADAeroMarkers%Hub%Position(3)=ADAeroMarkers%Hub%Position(3)+50
ADAeroMarkers%RotorFurl%Position(3)=ADAeroMarkers%RotorFurl%Position(3)+50
ADAeroMarkers%Nacelle%Position(3)=ADAeroMarkers%Nacelle%Position(3)+50
ADAeroMarkers%Tower%Position(3)=ADAeroMarkers%Tower%Position(3)+50
ADAeroMarkers%Tail%Position(3)=ADAeroMarkers%Tail%Position(3)+50
!write(,) ADAeromarkers%Blade(3,2)%Position

! Read a value if interested.
!do, i=1,3
!write(,) ( Rotatebottom(i,j), j=1,3 )
!enddo

! Read a value if interested.
!DO i = 1, 3
!PRINT*,(ADAeroMarkers%Blade(3,2)%Position(i))
!END DO
!write(,) ADInterfaceComponents%Rotorfurl%RotationVel
!do, i=1,3
!write(,) ( ADAeromarkers%Blade(1,1)%Orientation(i,j), j=1,3 )
!enddo
do time=0,1,0.1
ADAeroLoads = AD_CalculateLoads(time,ADAeroMarkers,ADInterfaceComponents,ALOptions,ErrStat)
write(,) 'Time is ', time
write(,) ‘Force components’, ADAeroLoads%Blade(1,1)%Force(:slight_smile:
write(,) ‘Moment components’, ADAeroLoads%Blade(1,1)%Moment(:slight_smile:
!write(,) 'Orientation of el1 bl1 ', ADAeromarkers%Blade(2,2)%Position(:,:slight_smile:
CALL ElemOut()
end do

CALL AD_Terminate(ErrStat)

end program Aeroding[/code]

10

On a related note, does anyone have experience with wrapping/calling AeroDyn in/from python?

The structural code I am going to couple AeroDyn with is written is C++ and has a python wrapper. The developer of the code recommended doing the coupling on the python level. I read about f2py but am already experiencing errors using f2py on the NWTC library.

I am still trying to get it working, but any existing knowledge on this topic could come in handy.

Sander

11

Sander,

I’m not aware on anyone who has called AeroDyn from Python. We do have people here who call Fortran routines from Python, but not AeroDyn specifically.

Marshall

12

Hi Sander,

A DTU student I met a couple of years ago told me he’d created a Python wrapper for AeroDyn. I don’t recall many details other than that he was waiting for a newer version of Python to deal with Fortran’s user-defined TYPES.