AeroDyn: Using interpolated inputs causes inaccurate results

[dustinjamescondon]

AD_UpdateStates(...) calls AD_Input_ExtrapInterp(...), which interpolates/extrapolates the inputs held in AD%u at 't' and 't + dt'. However, in AeroDyn_Driver, the inputs are always given at 't' and 't +dt', so no actual interpolation takes place. We can make AD_UpdateStates actually use interpolated inputs by offsetting the time of the simulation by dt/2:

[From AeroDyn_Driver.f90 (Only change is addition of time = AD%InputTime(2) + 0.5 * dT_Dvr to force AD_UpdateStates to actually interpolate the inputs)]

      do nt = 1, numSteps
         
         !...............................
         ! set AD inputs for nt (and keep values at nt-1 as well)
         !...............................
         
         call Set_AD_Inputs(iCase,nt,DvrData,AD,errStat,errMsg)
            call CheckError()
   
            
         ! Offset the time to be simulated to by half a time-step causing AD_UpdateStates to use 
         ! interpolated inputs
         time = AD%InputTime(2) + 0.5 * dT_Dvr
            
            ! Calculate outputs at nt - 1

         call AD_CalcOutput( time, AD%u(2), AD%p, AD%x, AD%xd, AD%z, AD%OtherState, AD%y, AD%m, errStat, errMsg )
            call CheckError()
   
         call Dvr_WriteOutputLine(DvrData%OutFileData, time, AD%y%WriteOutput, errStat, errMsg)
            call CheckError()
            
            
            ! Get state variables at next step: INPUT at step nt - 1, OUTPUT at step nt

         call AD_UpdateStates( time, nt-1, AD%u, AD%InputTime, AD%p, AD%x, AD%xd, AD%z, AD%OtherState, AD%m, errStat, errMsg )
            call CheckError()
      
                  
      end do !nt=1,numSteps

And when AD_UpdateStates(...) actually uses interpolated inputs,
the non-axial forces and moments are thrown off: particularly RtAeroFzh, RtAeroFyh , RtAeroMyh, and RtAeroMzh.

I found what seems to be causing it and how to fix it. AD_Input_ExtrapInterp(...) linearly interpolates all the orientation matrices (held in HubMotion, BladeRootMotion, and BladeMotion), and the blade node positions. If the BladeMotion orientation matrices are instead recalculated using the interpolated HubMotion and BladeRootMotion matrices, and then the blade node positions recalculated using that BladeMotion orientation matrix, the non-axial forces are accurate once again.

[From AeroDyn.f90's subroutine AD_UpdateStates (Only change is call RecalcBladeMotion(uInterp%HubMotion, uInterp%BladeRootMotion, uInterp%BladeMotion) after each interpolation call)]

call AD_CopyInput( u(1), uInterp, MESH_NEWCOPY, errStat2, errMsg2)
call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
if (ErrStat >= AbortErrLev) then
   call Cleanup()
   return
end if

! set values of m%BEMT_u(2) from inputs interpolated at t+dt:
call AD_Input_ExtrapInterp(u,utimes,uInterp,t+p%DT, errStat2, errMsg2)
call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
    
call RecalcBladeMotion(uInterp%HubMotion, uInterp%BladeRootMotion, uInterp%BladeMotion)

call SetInputs(p, uInterp, m, 2, errStat2, errMsg2)
call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)

! set values of m%BEMT_u(1) from inputs (uInterp) interpolated at t:
! I'm doing this second in case we want the other misc vars at t as before, but I don't think it matters
call AD_Input_ExtrapInterp(u,utimes,uInterp, t, errStat2, errMsg2)
call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
    
call RecalcBladeMotion(uInterp%HubMotion, uInterp%BladeRootMotion, uInterp%BladeMotion)

call SetInputs(p, uInterp, m, 1, errStat2, errMsg2)
call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)


! Call into the BEMT update states    NOTE:  This is a non-standard framework interface!!!!!  GJH
call BEMT_UpdateStates(t, n, m%BEMT_u(1), m%BEMT_u(2),  p%BEMT, x%BEMT, xd%BEMT, z%BEMT, OtherState%BEMT, p%AFI%AFInfo, m%BEMT, errStat2, errMsg2)
call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) 

call Cleanup()

contains

subroutine RecalcBladeMotion(hubMotion, bladeRootMotion, bladeMotion)
   type(MeshType),                    intent(in) :: hubMotion
   type(MeshType), dimension(:),      intent(in) :: bladeRootMotion
   type(MeshType), dimension(:),   intent(inout) :: bladeMotion
       
   real(DbKi), dimension(1:3, 1:3) :: rotateMat, orientation
   real(ReKi), dimension(1:3)      :: position
   integer(intKi) :: k,j ! loop counters

   do k=1, size(bladeRootMotion)
      rotateMat = transpose( BladeRootMotion(k)%Orientation(  :,:,1) )
      rotateMat = matmul( rotateMat, BladeRootMotion(k)%RefOrientation(  :,:,1) )
      orientation = transpose(rotateMat)

      rotateMat(1,1) = rotateMat(1,1) - 1.0_ReKi
      rotateMat(2,2) = rotateMat(2,2) - 1.0_ReKi
      rotateMat(3,3) = rotateMat(3,3) - 1.0_ReKi

      do j=1,BladeMotion(k)%nnodes
         position = BladeMotion(k)%Position(:,j) - hubMotion%Position(:,1)
         BladeMotion(k)%TranslationDisp(:,j) =  hubMotion%TranslationDisp(:,1) + matmul( rotateMat, position )

         BladeMotion(k)%Orientation(  :,:,j) = matmul( BladeMotion(k)%RefOrientation(:,:,j), orientation )

         position =  BladeMotion(k)%Position(:,j) + BladeMotion(k)%TranslationDisp(:,j) - hubMotion%Position(:,1) - hubMotion%TranslationDisp(:,1)
         BladeMotion(k)%TranslationVel( :,j) = cross_product( HubMotion%RotationVel(:,1), position )! + HubMotion%TranslationVel(:,1) ! added hub vel for PDS coupling

      end do !j=nnodes

   end do !k=numBlades
       
           
end subroutine RecalcBladeMotion

Although this additional code of recalculating BladeMotion after each interpolation appears to correct the forces when actually interpolating inputs, it isn't a perfect fix, because AD_Input_ExtrapInterp(...) is still doing calculations which are just being overwritten by the additional code. But AD_Input_ExtrapInterp(...) is generated by the registry system. So I'm hoping someone can think of a better way addressing the issue.

In the attached folder:
InterpolationIssue.zip

• Test01.1(normal).out contains the results from unedited AeroDyn_Driver
• Test01.1(offsetTime).out contains results when running the version of AeroDyn_Driver when AeroDyn_Driver.f90 has the additional line of code that offsets the simulation to force it to use interpolated inputs.
• Test01.1(offsetTimeAndFix).out contains the results when running the version of AeroDyn_Driver with same addition to AeroDyn_Driver.f90 and the additional BladeMotion recalcultion code in AeroDyn.f90 in the subroutine AD_UpdateStates(...)
• ad_driver_example2.inp is the AeroDyn_Driver input file (needs to be placed in the correct directory relative to the input files it references, but I believe any AeroDyn_Driver input file will expose the issue).

[andrew-platt]

I think this idea is addressed in PR #164, which hasn't been merged in yet.

[dustinjamescondon]

@andrew-platt It looks like the issue was addressed within FAST_Solver.f90 by overwriting predicted/extrapolated AeroDyn states with calculated ones in the subroutine SolveOption2b_Inp2IfW. However the issue is still unaddressed for anybody just using the AeroDyn module. I believe the issue I've come across and the one addressed in PR #164 (DBEMT not liking extrapolated inputs) may be caused by the same thing: AD_InputExtrapInterp giving inaccurate results. I'm wondering if fixing how AD inputs are extrapolated by AD_Input_ExtraInterp would be another approach fixing the issue addressed in PR #164.

By fixing I mean having AD_InputExtrapInterp only interpolate/extrapolate the HubMotion orientation, and using it to calculate what BladeRootMotion and BladeMotion's matrices should be, along with calculating where the blade node positions are given all of those previous things. Except for recalculating BladeRootMotion, this is what the subroutine RecalcBladeMotion in my above posted code is doing. The thing is it's only intercepting the AD inputs after they have all been interpolated, and then overwriting BladeMotion with its results because AD_InputExtrapInterp is registry generated.

[bjonkman]

@dustinjamescondon: I'm not sure these changes are necessary after PR #164 is merged. With AeroDyn, the extrap/interp routines are just a way of getting us inputs at t and t+dt. We should always have calculated results at t, and with PR #164, we will also have calculated results at t+dt. In those cases, extrap/interp would return the calculated values, and I wouldn't want to recalculate any of the inputs that the structural code gave me.

The AeroDyn driver should also be taking steps of dt (without the change you introduced, above), and it would be using calculated (not extrapolated) inputs there, too.

[dustinjamescondon]

@bjonkman But whether or not it's used in the FAST glue-code or AeroDyn driver, wouldn't it be good to have AD_Input_ExtrapInterp give a usable result? There may be users who are just looking to use the AeroDyn module for coupling to some other piece of simulation software and may have use for the existing input interpolation/extrapolation given by AD_Input_ExtrapInterp. But right now it's not really usable when it's called with a time that isn't the same as an already given input.

In those cases, extrap/interp would return the calculated values, and I wouldn't want to recalculate any of the inputs that the structural code gave me.

The code I posted is just to show the recalculation that seems to bring the interpolated/extrapolated inputs back to being usable again. I'm not sure exactly where this calculation should take place, but wherever that is, there could just be an additional check to see if any of the input times are close enough to the interpolation/extrapolation time to warrant just returning that input.

It does seem like a finicky thing to solve though because of the fact that the interpolation code is registry generated. I haven't seen the other modules in any detail, but I'm guessing the reason the boiler-plate interpolation code doesn't quite work very well in AeroDyn is because of how the different orientation matrices in AeroDyn depend on one another: BladeMotion depends on BladeRootMotion depends on HubMotion. What may also be a cause of the problem is that linearly interpolating/extrapolating the blade node positions puts the nodes farther out (extrapolation) or closer in (interpolation) relative to the hub center.

In any case having the issue posted will let people know a problem exists before using the results of AD_Input_ExtrapInterp.

[jjonkman]

Dear @dustinjamescondon,

The AD_Input_ExtrapInterp routine is registry generated and a similar routine is registry-generated for all OpenFAST modules. I'm not sure I fully understand why you say that the routine does not produce a "useable result". The routine does exactly what it is intended to do...interpolate or extrapolate the inputs using a linear fit (if inputs at two times are provided) or using a quadratic fit (if inputs at three times are provided). If the resulting interpolation or extrapolation is "not useable", this likely means that the time step is too large such that the linear or quadratic fit to the curve is not a good fit; in this case, I would lower the time step.

Regardless, as Bonnie and Andy said, PR #164 avoids the issue all together because the input to AeroDyn at t+dt is not extrapolated; instead, it is computed directly by calling the other modules calculating inputs for AeroDyn beforehand. When running AeroDyn from the standalone driver, it is also possible to prescribe the exact inputs at t+dt. If AeroDyn was coupled to another structural solver, then I would suggest coupling AeroDyn using a similar scheme to what is proposed in PR #164. Otherwise, a correction step or a smaller time step is likely necessary.

Also, please note that the time-integration scheme used by AeroDyn assumes a fixed time step of dt, so, while inputs are used by AeroDyn at t and t+dt, inputs are not used at t+dt/2.

Best regards,

[andrew-platt]

Closing this issue now that PR #164 has been merged and there hasn't been any further discussion. @dustinjamescondon, please reopen the issue if you find the latest release (https://github.com/OpenFAST/openfast/tree/v2.2.0) does not adequately address your concerns.