BeamDyn questions

Wind & Water Structural Analysis
1

Dear @Jason.Jonkman

I have been using BeamDyn for structural analysis of blades and have observed some phenomena. The details are as follows:

Under gravity-free conditions, applying only a RootVel to the blade results in very small deformation.

图片
图片508×289 30.4 KB
图片
图片508×289 32.4 KB

图片
图片508×289 16.1 KB

Under conditions of no rotation and no gravity, applying a distributed moment in the z-direction of the local coordinate system produces zero translational deformation in the x- and y-directions, and only a torsional deformation in the z-direction that decays over time.

图片
图片508×289 3.91 KB
图片
图片508×289 3.92 KB

图片
图片508×289 13.7 KB

However, when both a root rotational velocity and the same z-direction distributed moment are applied simultaneously, a significant periodic deformation occurs, with a period matching the rotational period.

图片
图片508×289 23.8 KB
图片
图片508×289 26.6 KB

图片
图片497×288 18.8 KB

Is the structural solution computed in the blade’s local coordinate system? If so, the externally applied loads should theoretically be transformed from the global coordinate system to the local coordinate system. However, while reviewing the code, I noticed that the SetRefFramefunction does not seem to function as expected in calculating the rotation matrix—the rotation matrix remains consistent with the settings in the input file throughout the computation, without any changes. Does this imply that the solution is actually computed in the fixed initial blade local coordinate system?

I observed that the external loads are assembled in the blade’s local coordinate system at the initial time step. What variables are affected by rotation that cause this periodic response in the results? Could it be that I haven’t identified the correct location where the load transformation takes place?

Is the rotational period observed in the deformation results due to the loads being dead loads , leading to periodic variations in the rotating local coordinate system)?

I aim to use BeamDyn to input external loads in the blade’s body-fixed coordinate system and solve for the resulting deformations. I would greatly appreciate your suggestions or guidance on how to correctly achieve this objective and how to interpret the phenomena described above.

Thank you for taking the time to read this email. I look forward to your response.

Best regards,

2

Dear @Min.Li,

Can you clarify which version of OpenFAST you are running the BeamDyn driver from? There have been many improvements and bug fixes made to BeamDyn over the years. If you are not already, I suggest upgrading to the newest version of OpenFAST (currently v4.1.2).

Are you using a BeamDyn model provided by NLR or is this a custom-made model?

Regarding your case with blade rotation without applied loads, which vector are you rotating about?

Regarding your case with blade rotation and distributed loads, are you specifying the load through the driver input file or by customization of the source code? The BeamDyn driver documentation (4.6.3. Input Files — OpenFAST v4.1.2 documentation) mentions that the loads specified through the input file do not follow the blade rotation, and summarizes how the source code can be customized for your application.

Another option to address your use case is to simulate BeamDyn using OpenFAST as the driver (with all DOFs disabled in ElastoDyn) and apply prescribed loads that follow the blade as it rotates and deflects through the Structural Control (StC) submodel of ServoDyn.

Best regards,

3

Dear @Jason.Jonkman,

A master student Muhammad Talha Sajjad and myself are currently working on the last version of BeamDyn (the one included in OpenFast V5.0).

We are trying to run a simple case in standalone version of BeamDynv5.0.0.
A dynamic simulation is performed where a constant root rotational velocity (omega = -0.628 rad/s) is imposed on the beam with a square cross section with the properties defined in the attached input file. The beam is oriented vertically with its longitudinal axis aligned with the global Z axis. It was observed that the Root Reaction Fzr does not correspond to the analytical solution (F = 1/2*rho*A*omega^2*L^2). We modified the source code and switched the ChangeRefFrame in BeamDyn.f90 file from True to False, and with this configuration, the result is close to the analytical one as shown in figure. The input and output files are attached in the zip.

Note that older version of BeamDyn give consistant results too, similar to the one with ChangeRefFrame=False.

Result
Result547×453 24 KB

The input files we used can be downloaded here:

Best regards,

Alban

PS: independently of the issue, I noticed a bug in the NWTC_Num.f90 subroutine, which is not fixed yet. In the function Quaternion_to_DCM (around line 4809), q3q3 is defined as q%v(2)**2 whereas it should be q%v(3)**2.

4

Dear @Alban.Leroyer,

Thanks for notifying us about this. It will take some investigation to figure out where the error is in BeamDyn. I will look into it and post a fix for it.

Also, thank you for alerting us to the typo in the Quaternion_to_DCM routine. This routine was added in 2014, but never actually used. I have posted a pull request to fix it here: [BugFix] Error in `Quaternion_to_DCM` calculation by andrew-platt · Pull Request #3355 · OpenFAST/openfast · GitHub

Regards,

5

Dear @Andy.Platt , dear @Jason.Jonkman,

Thank you for your feedback. We hope that the reported issue will be fixed soon.

We are continuing our validation tests with BeamDyn and are facing a new issue, unfortunately.

We are trying to test the case of a quarter circular beam in global YZ plane. We are interested to get the orientation matrices of each node (BD_Output%BldMotion%Orientation). From the BeamDyn output, the orientation matrix of the third last node (BD_Output%BldMotion%NNodes - 2), does not seem correct (the sign of rotation matrix element [2,3] and [3,2] is not correct, see Node 28 of the Output_BeamDyn.txt in the given URL below).

We checked the .sum.yaml and found that the WM parameter corresponding to the first line of the last element (Element numer 10) has a flipped sign (highlighted below), see curved_beam_bd_driver.BD.sum.yaml.

# Element number: 9
Init_Nodes_E9: # 3 x 6

  • [ 0.00000E+00, 6.90983E-01, 9.51057E-01, -1.25283E+00, 0.00000E+00, 0.00000E+00]
  • [ 0.00000E+00, 7.53706E-01, 9.69372E-01, -1.38709E+00, 0.00000E+00, 0.00000E+00]
  • [ 0.00000E+00, 8.43566E-01, 9.87688E-01, -1.45505E+00**,** 0.00000E+00, 0.00000E+00]

Element number: 10

Init_Nodes_E10: # 3 x 6

  • [ 0.00000E+00, 8.43566E-01, 9.87688E-01, 1.39473E+00, 0.00000E+00, 0.00000E+00]
  • [ 0.00000E+00, 8.69457E-01, 9.93844E-01, -1.56557E+00, 0.00000E+00, 0.00000E+00]
  • [ 0.00000E+00, 1.00000E+00, 1.00000E+00, -1.61757E+00, 0.00000E+00, 0.00000E+00]

We do not expect an issue due to a singularity for a quarter circular beam.

The input files and the result files can be downloaded at the following URL:

Regards,

Alban Leroyer and Muhammad Talha Sajjad

6

Dear @Andy.Platt,

For a case with large motion of the root, I had to enable the ChangeRefFrame option to avoid divergence issues, even if this option seems to have a bug.

After some investigations, I found that In BD_BoundaryGA2, x%q(1) seems to be set to the global displacement respect to the inertial frame even if ChangeRefFrame is true

x%q(1:3,1) = u%RootMotion%TranslationDisp(1:3,1) + &
matmul(u%RootMotion%Position(:,1) - OtherState%GlbPos, OtherState%GlbRot)

whereas when ChangeRefFrame is true, in BD_UpdateGlobalRef , x%q(1:3, temp_id) is set to the relative motion with respect to the root frame. It seems also that x%dqdt(1:3,1) and x%dqqt(4:6,1) are not consistent too.

In addition to that, there is a bug in this routine BD_UpdateGlobalRef since applying the two loops indexed i and j , temp_id is called twice for some nodes. As a consequence, x%q(1:3, temp_id) is wrong for these nodes as the formulation depends on x%q(1:3, temp_id) itself

do i = 1, p%elem_total
do j = 1, p%nodes_per_elem
! The last node of the first element is used as the first node in the second element.
     temp_id = (i - 1)*(p%nodes_per_elem - 1) + j

     ! Calculate displacement in terms of new root motion mesh position
     x%q(1:3, temp_id) = ... ... + x%q(1:3, temp_id)

 ! The last node of the first element is used as the first node in the second element.     
temp_id = (i - 1)*(p%nodes_per_elem - 1) + j     
! Calculate displacement in terms of new root motion mesh position     
x%q(1:3, temp_id) = ... ... + x%q(1:3, temp_id)

I hope this can help to fix the bug as soon as possible.

Regards,

Alban Leroyer