Rotational Displacement Output Bug in BeamDyn

[jjonkman]

While the translational displacement outputs from BeamDyn are correct, the corresponding rotational displacement outputs from BeamDyn are not. More information is discussed in the following forum topic: https://wind.nrel.gov/forum/wind/viewtopic.php?f=4&t=1564&p=7780.

[michaelasprague]

@jjonkman, do you know if there's a straightforward way to expose this bug in the stand-alone beamdyn?

[jjonkman]

@michaelasprague, I have not set-up a BeamDyn standalone case to test this problem myself.

The forum topic linked above mentions that converting the BeamDyn rotational parameter outputs to a rotation matrix following the equations in the manual does not result in a rotation matrix with a determinant of one. This could be one clue.

From what I've seen examining the BeamDyn rotational parameter outputs from Test26 in the FAST CertTest, the outputs (1) did not seem properly consistent between the 3 blades and (2) did not match well with the tip rotations output by ElastoDyn from Test18.

[michaelasprague]

@jjonkman, I ran some tests and the rotation parameters correctly correspond to rotation matrices, i.e., det(R)=1 and R^T R = I. Please let me know if you still see discrepancies.

[drWindstrom]

@michaelasprague, I just got an email from you about my post "BeamDyn: Convert Wiener-Milenkovic to rotation matrix" a while ago. I ended up using this Matlab Script for converting:

function [R, Xphi, Ytheta, Zpsi] = wieMilToR(rx, ry, rz)
% ---------------------------------------------------------------------
% WIEMILTOR.M Converts Wiener-Milenkovic parameters to rotation matrix
%
% Version 1.0    
% Last modified  Winstroth 08 July 2016 
% Created        Winstroth 08 July 2016   
%
% This function converts the Wiener-Milenkovic rotation parameters used
% by NREL FAST v8 to a rotation matrix and also to Tait-Bryan angles
% with the rotation sequence x-y'-z''.
% ---------------------------------------------------------------------
% Input
% rx  x-component of the Wiener-Milenkovic parameter
%
% ry  y-component of the Wiener-Milenkovic parameter
%
% rz  z-component of the Wiener-Milenkovic parameter
% ---------------------------------------------------------------------
% Output
% R  3x3 Rotation matrix that corresponds to the Wiener-Milenkovic 
%    parameters [rx; ry; rz]
%
% Xphi  Rotation angle of the first rotation about the x-axis
%
% Ytheta  Rotation angle of the second rotation about the y'-axis
%
% Zpsi  Rotation angle of the third rotation about the z''-axis
% ---------------------------------------------------------------------

% Create Wiener-Milenkovic vector
c = [rx; ry; rz];

% Convert to rotation matrix
c0 = 2.0 - 1.0/8.0*(c')*c;

R11 = (c0^2 + rx^2 - ry^2 - rz^2);
R12 = (2.0*(rx*ry - c0*rz));
R13 = (2.0*(rx*rz + c0*ry));

R21 = (2.0*(rx*ry + c0*rz));
R22 = (c0^2 - rx^2 + ry^2 - rz^2);
R23 = (2.0*(ry*rz - c0*rx));

R31 = (2.0*(rx*rz - c0*ry));
R32 = (2.0*(ry*rz + c0*rx));
R33 = (c0^2 - rx^2 - ry^2 + rz^2);

R = zeros(3, 3);

R(1,1,:) = R11;
R(1,2,:) = R12;
R(1,3,:) = R13;

R(2,1,:) = R21;
R(2,2,:) = R22;
R(2,3,:) = R23;

R(3,1,:) = R31;
R(3,2,:) = R32;
R(3,3,:) = R33;

scale_factor = 1.0/(4.0 - c0)^2;
R = scale_factor * R;

Xphi = atan2d(-R(2,3), R(3,3));
Ytheta = asind(R(1,3));
Zpsi = atan2d(-R(1,2), R(1,1));

end

I don't remember where I got the equations from. Possibly the BeamDyn manual and the book "Flexible Multibody Dynamics" by O. A. Bauchau. Using the script, the resulting angles appeared fine to me, but I did not do any further checking.

Hopefully this helps. Regards

Jan

[michaelasprague]

Jan, Thanks. I believe those are indeed the equations that are found in the BD manual.

[michaelasprague]

@jjonkman, can we close this issue, or are you still seeing discrepancies in Test 26?

[jjonkman]

@michaelasprague -- The bug is still present. I've attached results from three models:

1. Test18 - Run without modification except for adding the tip-rotations for all three blades to the OutList,
2. Test26 - Run without modification, and
3. Test26 with Upwind Precurve - Test26 with the BeamDyn primary and AeroDyn blade input files modified with a small upwind precurve--a quadratic function from 0 m at the root to -1 m at the tip. These files are available from the following forum topic: https://wind.nrel.gov/forum/wind/viewtopic.php?f=4&t=1451&p=6563.

The results show that Test18 and Test26 give quite consistent results. While the tip rotation outputs from Test26 have a different convention from Test18, the trends from Test26 are consistent with those from Test18 and the results are as expected. However, the results from Test26 with Upwind Precurve are not as expected. In this test, blade 1 has a near zero mean tip twist, blade 2 has a negative mean tip twist, and blade 1 has a positive mean tip twist, despite all blades being identical. The mean RDyr rotations also look inconsistent between the three blades. The bug seems to be tied to blades with initial curvature.

TipRotations.pdf

[bjonkman]

line 1624 in BeamDyn_IO.f90 (the first CrvCompose in the NRDr output computation) is currently this:
CALL BD_CrvCompose(temp_vec2,temp_vec,temp_glb,FLAG_R1R2)
@andrew-platt noticed this is backwards from other similar places in the code. He changed it to
CALL BD_CrvCompose(temp_vec2,temp_glb,temp_vec,FLAG_R1R2) and it didn't seem to change results in the models we are running. However, this might be related to problems with the precurve....

[jjonkman]

As stated above, I just committed a fix to this bug to jjonkman/OpenFAST. I've made a pull request to OpenFAST/OpenFAST. @bjonkman was correct; line 1544 was also incorrect. The attached results show the effect of the bug fix on Test26 and Test26 with Upwind Precurve.
TipRotations.pdf

[michaelasprague]

Closing issue per @jjonkman commit.