Motion transfer from ElastoDyn to AeroDyn

Hello,

I am using FAST v8.15 with AeroDyn v15.
I have questions about the transfer of the motions from the mesh of ElastoDyn to the mesh used by AeroDyn.

  1. The angular velocities of the blade elements are not computed in ElastoDyn, because they are not used either in ElastoDyn nor AeroDyn. However, they are used when transferring the translation velocities in the routine Transfer_Motions_Line2_to_Point in the module ModMesh_Mapping. So not taking them into account would lead to a small error on the results ?
    Therefore I uncommented the computation of AngVelHM done in CalculateAngularPosVelPAcc in the module ElastoDyn, and transferred the array to BladeLn2Mesh%RotationVel at the end of ED_CalcOutput. Do you agree with this approach ? Is there a reason for not computing those angular velocities ?
    (note that I also need the angular velocities for the coupling of FAST with my own aerodynamic module)

  2. In the routine Transfer_Motions_Line2_to_Point, the computation of the direction cosine matrix has been replaced with a new version, apparently intended to save one matrix multiplication, and I guess, improve the performances. However, with this version, I get small oscillations in the orientation arrays every time the azimuth angle is close to 180°, which can also be observed on other variables such as the torque, even if the amplitude is very small. Have you ever observed such an issue ? Does this version of the DCM computation lead to significant improvement of the performances of the code ?

Best regards,
Bastien

Dear Bastien,

Regarding (1), good question. Actually, I thought we had already included calculations of the rotational velocity of the blade nodes in ElastoDyn, but looking at the code, you are correct that they have not been calculated. While AeroDyn doesn’t need the rotational velocity of the blade nodes, the AeroDyn nodes may be offset from the ElastoDyn nodes, so a rotational velocity of the ElastoDyn nodes should lead to a change in the translational velocity of the AeroDyn nodes. I will have to dig into this some more to see why we neglected these terms and to include them in a future release of FAST.

Regarding (2), I’m not sure I understand what you mean when you say a new version of the DCM has replaced an older version. I’m not aware of problems in the implementation around 180 degrees azimuth; please enlighten us as to what problems you are seeing.

Best regards,

Dear Jason,

thank you for your reply.

Regarding (2), I’m looking at this part of the code in the routine Transfer_Motions_Line2_to_Point :

#ifdef __ORIGINAL_LOGMAP
            ! calculate Rotation matrix for FieldValueN1 and convert to tensor:
         RotationMatrix = MATMUL( MATMUL( Dest%RefOrientation(:,:,i), TRANSPOSE( Src%RefOrientation(:,:,n1) ) )&
                                 , Src%Orientation(:,:,n1) )

         CALL DCM_logmap( RotationMatrix, FieldValue(:,1), ErrStat, ErrMsg )
         IF (ErrStat >= AbortErrLev) RETURN

            ! calculate Rotation matrix for FieldValueN2 and convert to tensor:
         RotationMatrix = MATMUL( MATMUL( Dest%RefOrientation(:,:,i), TRANSPOSE( Src%RefOrientation(:,:,n2) ) )&
                                 , Src%Orientation(:,:,n2) )

         CALL DCM_logmap( RotationMatrix, FieldValue(:,2), ErrStat, ErrMsg )
         IF (ErrStat >= AbortErrLev) RETURN

         CALL DCM_SetLogMapForInterp( FieldValue )  ! make sure we don't cross a 2pi boundary


            ! interpolate tensors:
         TmpVec =   MeshMap%MapMotions(i)%shape_fn(1)*FieldValue(:,1)  &
                  + MeshMap%MapMotions(i)%shape_fn(2)*FieldValue(:,2)

            ! convert back to DCM:
         Dest%Orientation(:,:,i) = DCM_exp( TmpVec )


#else
!this should be equivalent, with one less matrix multiply

         ! bjj: because of numerical issues when the angle of rotation is pi, (where
         ! DCM_exp( DCM_logmap (x) ) isn't quite x
      if ( EqualRealNos( MeshMap%MapMotions(i)%shape_fn(1), 1.0_ReKi ) ) then
         RotationMatrixD = MATMUL( TRANSPOSE( Src%RefOrientation(:,:,n1) ), Src%Orientation(:,:,n1) )

      elseif ( EqualRealNos( MeshMap%MapMotions(i)%shape_fn(2), 1.0_ReKi ) ) then
         RotationMatrixD = MATMUL( TRANSPOSE( Src%RefOrientation(:,:,n2) ), Src%Orientation(:,:,n2) )

      else

            ! calculate Rotation matrix for FieldValueN1 and convert to tensor:
         RotationMatrixD = MATMUL( TRANSPOSE( Src%RefOrientation(:,:,n1) ), Src%Orientation(:,:,n1) )

         CALL DCM_logmap( RotationMatrixD, FieldValue(:,1), ErrStat, ErrMsg, theta(1) )
         IF (ErrStat >= AbortErrLev) RETURN

            ! calculate Rotation matrix for FieldValueN2 and convert to tensor:
         RotationMatrixD = MATMUL( TRANSPOSE( Src%RefOrientation(:,:,n2) ), Src%Orientation(:,:,n2) )

         CALL DCM_logmap( RotationMatrixD, FieldValue(:,2), ErrStat, ErrMsg, theta(1) )
         IF (ErrStat >= AbortErrLev) RETURN

         CALL DCM_SetLogMapForInterp( FieldValue )  ! make sure we don't cross a 2pi boundary

            ! interpolate tensors:
         tensor_interp =   MeshMap%MapMotions(i)%shape_fn(1)*FieldValue(:,1)  &
                         + MeshMap%MapMotions(i)%shape_fn(2)*FieldValue(:,2)

            ! convert back to DCM:
         RotationMatrixD = DCM_exp( tensor_interp )

      end if

      RotationMatrix = REAL( RotationMatrixD, ReKi )
      Dest%Orientation(:,:,i) = MATMUL( Dest%RefOrientation(:,:,i), RotationMatrix  )

#endif

The block “ORIGINAL_LOGMAP” seems to have been replaced by the second block after the “#else” keyword.
With this second block, which is actually used in the v8.15, I have oscillations I don’t have if I use the original block instead, as you can see on the figure below, where I plotted the torque vs. the azimuth angle of the rotor, in a very simple simulation with rigid bodies and no controller.

Best regards,
Bastien

Dear Bastien,

I don’t recall the detailed history of the code development–and while I familiar with the theory, this section of code was written by Bonnie Jonkman, who now works at Envision Energy–but I noticed that this section of code in Transfer_Motions_Line2_to_Point has been updated quite a bit between FAST v8.15 and FAST v8.16 e.g. the “__ORIGINAL_LOGMAP” compiler directive no longer exists in FAST v8.16. Does upgrading from FAST v8.15 to FAST v8.16 resolve the problem you are seeing? FAST v8.16 was released over a year ago: nwtc.nrel.gov/FAST8.

Best regards,

Dear Jason,

first of all note that I can only see this issue in the particular case where all the DOFs of ElastoDyn are disabled, the gravity and the tilt angle are set to zero, the aerodynamic effects of the tower are disabled, and with a steady wind, using FAST v8.15.

Indeed, upgrading to v8.16 solves the issue. In v8.16, the last matrix multiplication used to compute Dest%Orientation is done before the result is converted back to single precision real numbers, unlike the v8.15, which solves the problem.

Thank you very much for your help,
Best regards,

Bastien

Dear Bastien,

OK, thanks for checking. I’ve flagged the missing blade-node rotational velocities in ElastoDyn as a bug that we’ll fix in a future version.

Best regards,