I am attempting to compute the OpenFAST outputs (specifically tower-top acceleration YawBrTAxp) from the internal states. From FASTMotions.doc, the tower-top acceleration is found as dot(EaO,b1), where B is the tower-top coordinate system in the earth reference frame (E). I compute EaO per FASTKinematics.doc using the partial linear velocities and accelerations. Attached is the Matlab code that performs this for me.
I have run OpenFAST with the OC3 spar and 5MW turbine to output all internal states to compare. The YawBrTAxp output from OpenFAST and the one that I compute using the code do not match at all. I have attached 2 plots: the direct output from OpenFAST has high initial transients around 1.4 m/s^2 and settles to a steady-state of 0. The second plot is computing YawBrTAxp using the kinematics and transforms as outlined in FASTMotions and FASTKinematics documents. I have checked and re-checked that the equations match -
Matlab Code Notes:
qddot is the vector of [qdot,qddot]. qd represents qdot. S## represents tower properties at the tower top/baseplate.
%% Tower-Top Accelerations
XvO_1 = qd_TFA1 + qd_TFA2;
XvO_2 = -(S11_TFA*q_TFA1*qd_TFA1 + S22_TFA*q_TFA2*qd_TFA2 + S12_TFA*(qd_TFA1*q_TFA2+q_TFA1*qd_TFA2) ...
+ S11_TSS*q_TSS1*qd_TSS1 + S22_TSS*q_TSS2*qd_TSS2 + S12_TSS*(qd_TSS1*q_TSS2+q_TSS1*qd_TSS2));
XvO_3 = qd_TSS1 + qd_TSS2;
XvO = XvO_1*A(1,:) + XvO_2*A(2,:) + XvO_3*A(3,:);
% Linear velocity of platform reference point in earth frame
EvZ = qd_Sg*Z(1,:) + qd_Hv*Z(2,:) - qd_Sw*Z(3,:);
% Angular velocity of tower base in earth frame
EwX = qd_R*Z(1,:) + qd_Y*Z(2,:) - qd_P*Z(3,:);
% Angular velocity of tower-top in earth frame
EwB = EwX + (dO1_TSS*qd_TSS1 + dO2_TSS*qd_TSS2)*A(1,:) - (dO1_TFA*qd_TFA1 + dO2_TFA*qd_TFA2)*A(3,:);
% Partial angular velocities of tower base in earth frame
EwX_R = Z(1,:);
EwX_P = -Z(3,:);
EwX_Y = Z(2,:);
% Linear velocity of tower top in earth frame
EvO = EvZ + XvO + cross(EwX,rZO);
% Partial linear velocities of tower-top in earth frame
EvO_Sg = Z(1,:);
EvO_Sw = -Z(3,:);
EvO_Hv = Z(2,:);
EvO_R = cross(EwX_R,rZO);
EvO_P = cross(EwX_P,rZO);
EvO_Y = cross(EwX_Y,rZO);
EvO_TFA1 = A(1,:) - (S11_TFA*q_TFA1 + S12_TFA*q_TFA2)*A(2,:);
EvO_TSS1 = A(3,:) - (S11_TSS*q_TSS1 + S12_TSS*q_TSS2)*A(2,:);
EvO_TFA2 = A(1,:) - (S22_TFA*q_TFA2 + S12_TFA*q_TFA1)*A(2,:);
EvO_TSS2 = A(3,:) - (S22_TSS*q_TSS2 + S12_TSS*q_TSS1)*A(2,:);
% Partial linear accelerations of tower top in earth frame
dEvO_R = cross(EwX_R, XvO + cross(EwX,rZO));
dEvO_P = cross(EwX_P, XvO + cross(EwX,rZO));
dEvO_Y = cross(EwX_Y, XvO + cross(EwX,rZO));
dEvO_TFA1 = -(S11_TFA*qd_TFA1 + S12_TFA*qd_TFA1)*A(2,:) + cross(EwX, EvO_TFA1);
dEvO_TSS1 = -(S11_TSS*qd_TSS1 + S12_TSS*qd_TSS2)*A(2,:) + cross(EwX, EvO_TSS1);
dEvO_TFA2 = -(S22_TFA*qd_TFA2 + S12_TFA*qd_TFA1)*A(2,:) + cross(EwX, EvO_TFA2);
dEvO_TSS2 = -(S22_TSS*qd_TSS2 + S12_TSS*qd_TSS1)*A(2,:) + cross(EwX, EvO_TSS2);
% % Term 1
term_1 = EvO_Sg*qddot_Sg + EvO_Sw*qddot_Sw + EvO_Hv*qddot_Hv + ...
EvO_R*qddot_R + EvO_P*qddot_P + EvO_Y*qddot_Y + ...
EvO_TFA1*qddot_TFA1 + EvO_TFA2*qddot_TFA2 + EvO_TSS1*qddot_TSS1 + EvO_TSS2*qddot_TSS2;
% Term 2
term_2 = dEvO_R*qd_R + dEvO_P*qd_P + dEvO_Y*qd_Y + ...
dEvO_TFA1*qd_TFA1 + dEvO_TSS1*qd_TSS1 + dEvO_TFA2*qd_TFA2 + dEvO_TSS2*qd_TSS2;
% Tower top acceleration in earth frame
EaO(i,:) = term_1 + term_2;
% YawBrTAxp(i,:) = term_1 + term_2;
YawBrTAxp(i) = dot(EaO(i,:),B(1,:));
YawBrTAyp(i) = dot(-EaO(i,:),B(3,:));
YawBrTAzp(i) = dot(EaO(i,:),B(2,:));