Verification of BeamDyn for pure bending problems

Min.Li · March 13, 2025, 9:00am

Dear everyone,

I am verifying BeamDyn’s accuracy by simulating a cantilever beam’s pure bending, but the results are wrong. Did anyone else have the same problem?

Here is my case setup, references, and results.
case setup

references
https://www.nrel.gov/docs/fy15osti/63165.pdf

result

Even when I use BeamDyn to test small deformations of a cantilever beam, applying a distributed load or a concentrated force at the end also causes similar issues, resulting in strange curvatures. For example, the test results with a distributed load show the following figure.

The simulation’s boundary conditions might not be fixed-free. I only found the BD_BoundaryGA2 subroutine in the code related to boundary conditions, which sets the root’s velocity, displacement, and acceleration. Shouldn’t constraints be applied in the global stiffness matrix?

Where in the code are the boundary conditions set? Are they fixed-free?

If the boundary conditions aren’t the issue, what else could it be? Do you have any suggestions?

Many thanks in advance for any help.

Best regards,

Jason.Jonkman · March 13, 2025, 5:36pm

Dear @Min.Li,

The BeamDyn driver assumes fixed-free boundary conditions.

I don’t see this exact test in the paper you reference; did you cite the correct paper?

I do see that your beam has some twist at the tip (initial_twist); is this correct? Perhaps that twist is explaining why your results differ from what you expect.

Best regards

Min.Li · March 19, 2025, 1:53pm

Dear Dr. Jonkman,

Thank you for your reply.

I’m sorry the cited article is indeed incorrect. I cited https://www.osti.gov/servlets/purl/1371521).（Example1）

I changed the British units in the stiffness matrix to the SI. In the first post, I used one element and the lambda was 0.5 and 1.0.

The initial_twist setting was to test its impact on results. I’ve also tested the case where the initial_twist is 0, which doesn’t change the result. I’ve attached my case and results for your reference.

case

------- BEAMDYN V1.00.00 Driver INPUT FILE -------------------------------------
Static analysis of NREL 5MW blade under distributed force
---------------------- SIMULATION CONTROL --------------------------------------
0         t_initial     - Starting time of simulation (s)
30.0       t_final       - Ending time of simulation   (s)
0.002     dt            - Time increment size         (s)
---------------------- GRAVITY PARAMETER --------------------------------------
0.0        Gx            - Component of gravity vector along X direction (m/s^2)
0.0        Gy            - Component of gravity vector along Y direction (m/s^2)
0.0        Gz            - Component of gravity vector along Z direction (m/s^2)
---------------------- FRAME PARAMETER --------------------------------------
0.0        GlbPos(1)     - Component of position vector of the reference blade frame along X direction (m)
0.0        GlbPos(2)     - Component of position vector of the reference blade frame along Y direction (m)
0.0        GlbPos(3)     - Component of position vector of the reference blade frame along Z direction (m)
---The following 3 by 3 matrix is the direction cosine matirx ,GlbDCM(3,3),
---relates global frame to reference blade frame
1.000E+00    0.000E+00   0.000E+00
0.000E+00    1.000E+00   0.000E+00
0.000E+00    0.000E+00   1.000E+00
---------------------- ROOT VELOCITY PARAMETER ----------------------------------
0.0        RootVel(4)     - Component of angular velocity vector of the beam root about X axis (rad/s)
0.0        RootVel(5)     - Component of angular velocity vector of the beam root about Y axis (rad/s)
0.0        RootVel(6)     - Component of angular velocity vector of the beam root about Z axis (rad/s)
---------------------- APPLIED FORCE ----------------------------------
0.0        DistrLoad(1)  - Component of distributed force vector along X direction (N/m)
0.0        DistrLoad(2)  - Component of distributed force vector along Y direction (N/m)
0.0        DistrLoad(3)  - Component of distributed force vector along Z direction (N/m)
0.0        DistrLoad(4)  - Component of distributed moment vector along X direction (N-m/m)
0.0        DistrLoad(5)  - Component of distributed moment vector along Y direction (N-m/m)
0.0        DistrLoad(6)  - Component of distributed moment vector along Z direction (N-m/m)
0.0        TipLoad(1)    - Component of concentrated force vector at blade tip along X direction (N)
0.0        TipLoad(2)    - Component of concentrated force vector at blade tip along Y direction (N)
0.0        TipLoad(3)    - Component of concentrated force vector at blade tip along Z direction (N)
3084.5196	   TipLoad(4)    - Component of concentrated moment vector at blade tip along X direction (N-m)
0.0        TipLoad(5)    - Component of concentrated moment vector at blade tip along Y direction (N-m)
0.0        TipLoad(6)    - Component of concentrated moment vector at blade tip along Z direction (N-m)  
---------------------- PRIMARY INPUT FILE --------------------------------------
"bd_primary.inp"    InputFile - Name of the primary input file

--------- BEAMDYN V1.01.* INPUT FILE -------------------------------------------
NREL 5MW blade primary input file
---------------------- SIMULATION CONTROL --------------------------------------
True     Echo           - Echo input data to "<RootName>.ech" (flag)
1         analysis_type  - 1: Static analysis; 2: Dynamic analysis (switch)
0.0       rhoinf         - Numerical Damping Parameter for Generalized-alpha integrator
1         quadrature     - 1: Gauss; 2: Trapezoidal (switch)
DEFAULT   refine         - Refinement factor for quadrature 2 (-). DEFAULT = 1
DEFAULT   n_fact         - Factorization frequency (-). DEFAULT = 5
DEFAULT   DTBeam         - Time step size (s).
DEFAULT   NRMax          - Max number of iterations in Newton-Ralphson algorithm (-). DEFAULT = 10
DEFAULT   stop_tol       - Tolerance for stopping criterion (-)
---------------------- GEOMETRY PARAMETER --------------------------------------
2     member_total       - Total number of members (-)
11    kp_total           - Total number of key points (-)
1	6                  - Member number; Number of key points in this member
2	6 
kp_xr kp_yr   kp_zr    initial_twist
(m)    (m)    (m)     (deg)
  0.000000   0.000000   0.000000   0.000000
  0.000000   0.000000   0.025400   0.000000
  0.000000   0.000000   0.050800   0.000000
  0.000000   0.000000   0.076200   0.000000
  0.000000   0.000000   0.101600   0.000000  
  0.000000   0.000000   0.127000   0.000000
  0.000000   0.000000   0.152400   0.000000
  0.000000   0.000000   0.177800   0.000000
  0.000000   0.000000   0.203200   0.000000
  0.000000   0.000000   0.228600   0.000000
  0.000000   0.000000   0.254000   0.000000  
---------------------- MESH PARAMETER ------------------------------------------
5         order_elem     - Order of interpolation (basis) function (-)
---------------------- MATERIAL PARAMETER --------------------------------------
"beam_props.inp"  BldFile - Name of file containing properties for blade (quoted string)
---------------------- PITCH ACTUATOR PARAMETERS -------------------------------
False     UsePitchAct - Whether a pitch actuator should be used (flag)
200       PitchJ      - Pitch actuator inertia (kg-m^2) [used only when UsePitchAct is true]
2.0E+7    PitchK      - Pitch actuator stiffness (kg-m^2/s^2) [used only when UsePitchAct is true]
5.0E+5    PitchC      - Pitch actuator damping (kg-m^2/s) [used only when UsePitchAct is true]
---------------------- OUTPUTS -------------------------------------------------
True      SumPrint       - Print summary data to "<RootName>.sum" (flag)
"ES10.3E2"    OutFmt          - Format used for text tabular output, excluding the time channel.
 5        NNodeOuts      - Number of nodes to output to file [0 - 9] (-)
6,7,8,9,10    OutNd          - Nodes whose values will be output  (-)
          OutList        - The next line(s) contains a list of output parameters.  See OutListParameters.xlsx.

"N1TDxr,N1TDyr,N1TDzr,N1RDxr,N1RDxr,N1RDxr"
"N2TDxr,N2TDyr,N2TDzr,N2RDxr,N2RDxr,N2RDxr"
"N3TDxr,N3TDyr,N3TDzr,N3RDxr,N3RDxr,N3RDxr"
"N4TDxr,N4TDyr,N4TDzr,N4RDxr,N4RDxr,N4RDxr"
"N5TDxr,N5TDyr,N5TDzr,N5RDxr,N5RDxr,N5RDxr"
"TipTDxr, TipTDyr, TipTDzr"
"TipRDxr, TipRDyr, TipRDzr"
END of input file (the word "END" must appear in the first 3 columns of this last OutList line)
--------------------------------------------------------------------------------

 ------- BEAMDYN V1.00.* INDIVIDUAL BLADE INPUT FILE --------------------------
 NREL5MW Blade
 ---------------------- BLADE PARAMETERS --------------------------------------
 3   station_total    - Number of blade input stations (-)
 0   damp_type        - Damping type (switch): 0: no damping; 1: viscous damping
  ---------------------- DAMPING COEFFICIENT------------------------------------
   mu1        mu2        mu3        mu4        mu5        mu6
   (s)        (s)        (s)        (s)        (s)        (s)
1.0E-03    1.0E-03    1.0E-03    1.0E-03    1.0E-03    1.0E-03
 ---------------------- DISTRIBUTED PROPERTIES---------------------------------
  0.00
   7.873300E+06    0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00
   0.000000E+00    7.873300E+06    0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00
   0.000000E+00    0.000000E+00    7.873300E+06    0.000000E+00    0.000000E+00    0.000000E+00
   0.000000E+00    0.000000E+00    0.000000E+00    2.493856E+02    0.000000E+00    0.000000E+00
   0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00    6.170070E+02    0.000000E+00
   0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00    2.341760E+01

   1.975000E+01    0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00
   0.000000E+00    1.975000E+01    0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00
   0.000000E+00    0.000000E+00    1.975000E+01    0.000000E+00    0.000000E+00    0.000000E+00
   0.000000E+00    0.000000E+00    0.000000E+00    6.583333E+00    0.000000E+00    0.000000E+00
   0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00    6.583333E+00    0.000000E+00
   0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00    1.316667E+01

  0.50
   7.873300E+06    0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00
   0.000000E+00    7.873300E+06    0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00
   0.000000E+00    0.000000E+00    7.873300E+06    0.000000E+00    0.000000E+00    0.000000E+00
   0.000000E+00    0.000000E+00    0.000000E+00    2.493856E+02    0.000000E+00    0.000000E+00
   0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00    6.170070E+02    0.000000E+00
   0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00    2.341760E+01

   1.975000E+01    0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00
   0.000000E+00    1.975000E+01    0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00
   0.000000E+00    0.000000E+00    1.975000E+01    0.000000E+00    0.000000E+00    0.000000E+00
   0.000000E+00    0.000000E+00    0.000000E+00    6.583333E+00    0.000000E+00    0.000000E+00
   0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00    6.583333E+00    0.000000E+00
   0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00    1.316667E+01

  1.00
   7.873300E+06    0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00
   0.000000E+00    7.873300E+06    0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00
   0.000000E+00    0.000000E+00    7.873300E+06    0.000000E+00    0.000000E+00    0.000000E+00
   0.000000E+00    0.000000E+00    0.000000E+00    2.493856E+02    0.000000E+00    0.000000E+00
   0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00    6.170070E+02    0.000000E+00
   0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00    2.341760E+01

   1.975000E+01    0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00
   0.000000E+00    1.975000E+01    0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00
   0.000000E+00    0.000000E+00    1.975000E+01    0.000000E+00    0.000000E+00    0.000000E+00
   0.000000E+00    0.000000E+00    0.000000E+00    6.583333E+00    0.000000E+00    0.000000E+00
   0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00    6.583333E+00    0.000000E+00
   0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00    1.316667E+01

I still used lambda equals 1.0 as I can conveniently plot the theoretical solution.

I think there’s still a strange curvature at the last node of the first element, which is causing the same issue at the first node of the second element.

result

I would be extremely grateful if you could offer some guidance and suggestions.

Best regards，

Jason.Jonkman · March 19, 2025, 2:44pm

Dear @Min.Li,

Thanks for clarifying.

Can you also clarify what you are plotting? It likes like you are using Gaussian quadrature, so, the output nodes are coincident with the FE nodes, but it looks like you are only outputting the displacements at nodes 6,7,8,9,10. Are you rerunning BModes multiple times to output the displacements at all 11 FE nodes?

How does the solution compare if you use 1 member and increase the order of the element instead of the number of elements (order_elem = 10 would use 11 FE nodes with only 1 element).

Best regards,

Topic		Replies	Views
Verification of a new tool chain that uses BeamDyn Structural Analysis	8	4883	January 22, 2019
Static structural test on rotor blade comparing BeamDyn and ElastoDyn Testing	4	527	September 29, 2023
BeamDyn effect on structral loading measurements in region 3/ 5 MW wind turbine Structural Analysis	4	480	April 22, 2022
Eigenanalysis of an isotropic beam using BeamDyn Structural Analysis	1	1776	December 2, 2019
Simulations using BeamDyn Structural Analysis	5	4595	July 20, 2020

Verification of BeamDyn for pure bending problems

Related topics