New spring element in SubDyn: 6 by 6 stiffness matrix

[RBergua]

I believe this pull request is ready to be merged.

Feature or improvement description
One new spring element has been included in SubDyn.
The spring is a new member type. Accordingly, SubDyn now supports as members: beams, cables, rigid elements and springs.
Previously, SubDyn could only include a stiffness matrix as a boundary condition.

The spring element has the next properties:

1. It provides a 6 by 6 stiffness matrix (including diagonal and cross-coupling terms). The stiffness matrix is assumed to be symmetric (Kij = Kji).
   
2. It doesn't have a mass/inertia associated. If interested, the user may include a lumped mass at a given JointID.
3. The spring element must be defined between nodes (i.e., JointIDs) that are coincident. The code checks if this condition is true and aborts the program if the condition is not met.
4. Since the spring element is defined between two coincident nodes, the user must provide the spring orientation by means of the direction cosine.
5. When defining a spring element, the user must provide 22 inputs (the PropSetID and 21 coefficients to define the upper half of the stiffness matrix).
6. The spring is not subdivided by means of NDiv.
7. The sensors definition for the spring element (e.g., force, moment, displacement, rotation...) are the same as for the beams.
8. The spring can be connected to beams, kinematic joints (e.g., revolute joint, universal joint, and spherical joint), the interface joint and rigid links. Avoid connecting one spring directly to one cable.

The new SubDyn input file includes a new section for the spring element. See below for reference:

Internally, the spring element relates 2 nodes with 6 DOFs each. Accordingly, mathematically the spring element has a 12 by 12 dimension.

Impacted areas of the software
Only SubDyn is impacted. Since this is a new element, the verification tests should pass as the legacy elements in the code have not been modified.

Additional supporting information
Having the capability to include a stiffness matrix is very helpful for modelers. This stiffness matrix can be used for example to model:

1. Elastic couplings.
2. Bearings. Some OpenFAST users include part of the wind turbine drivetrain in SubDyn. The loads could be transferred to the supporting structure through one or more spring elements. Tapered roller bearings could also be represented taking advantage of the cross-coupling terms.
3. Sensors with compliance for models at scale.
4. The spring element could also be used to approximate revolute joints, universal joints or spherical joint. For example, one revolute joint could be considered by defining the diagonal coefficients as stiff (e.g., 1E12, don't provide very high values or the system will be ill-conditioned), except the direction with the DOF of interest (e.g., K44). It's important to remember at this point that SubDyn assumes that structural motions within the body are small.
5. Many more!

Test results, if applicable
Several test cases have been reviewed for the new spring element in SubDyn. Below some representative test cases are included. Results compare the output from SubDyn vs reference results (computed analytically or by means of different codes like HOTINT, SAMCEF and OpenModelica).

Checklist

• Consider adding a new subdyn driver r-test (with minimal outputs, few time steps, etc. for limited run time and disk usage, but still testing the new feature)
• Document the new feature in the documentation
• Add the API changes to the documenetation
• Update the input files in the r-test
• Final review

[RBergua]

Test Case 1: Spring and lumped mass.
One spring with a vertical stiffness of 3947.8 N/m and a lumped mass of 100 kg is considered. For the spring, one side is connected to the interface joint and the other side is attached to the lumped mass. The gravity acceleration is 9.81 m/s^2. For this condition, the system should oscillate around -0.2485 m (-9.81*100/3947.8) with a frequency of 1 Hz.

Test Case 2: Same spring and lumped mass as test case 1, but 5% damping considered.

Test Case 3: Comparison between SubDyn revolute joint and spring element.
One pendulum-type system is evaluated. See schematic representation below for reference. One revolute joint or a spring element is defined at point B. The rotational stiffness is 1E6 Nm/rad. The AB beam is 1 m long and the BC beam is 4 m long. Both beams have the same sectional properties (2 m diameter for a full circle beam) and material properties (Young Modulus = 210E9 N/m^2, Poisson ratio = 0.3 and density = 7860 kg/m^3).
The SubDyn interface joint is located at point A. Accoding to SubDyn’s methodology, it’s a floating case (there are no reaction joints). The Guyan modes correspond to rigid body modes and the gravitational loads are applied to the Guyan modes. This means that the system will only experience the elastic restoring from the revolute joint or the spring element, but not the gravitational restoring.

One force to reproduce a free-decay is applied at point C. This nonlinear force is as follows:
#Time_[s] , Fx_[N] , Fy_[N] , Fz_[N] , Mx_[Nm] , My_[Nm] , Mz_[Nm]
0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0
25 , 2.5E4 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0
50 , 2.5E4 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0
50.01 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0
100.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0

The angle rotated at the point B is computed by looking at the displacement of point C with regarding to point B.

Test Case 4: Beams and springs.
One system with 2 beams is considered. See the below schematic representation for reference. The 2 beams use a Timoshenko formulation and have the same dimensions (5 m long each, 1 m external diameter and 0.025 m thickness) and material properties (Young Modulus = 210E9 N/m^2, Poisson ratio = 0.3 and density = 7860 kg/m^3). The node 1 is clamped (interface joint for the SubDyn standalone model).
The first 4 eigenfrequies of the system are: f1 = 9.8 Hz, f2 = 9.8 Hz, f3 = 56.5 Hz and f4 = 56.5 Hz.
Gravity acceleration is not considered and a steady force along the x direction of 2E6 N is applied at node 3.

The deflection of the system looks as follows:

As it can be observed, the deflection along the x-direction does not involve a deflection in the z-direction duet to the linear assumption used in SubDyn.

This time, one spring element is introduced between the two beams. See the below representation for reference. The spring is defined between nodes 2 and 3 (coincident). The spring provides some compliance in the x-direction (k = 2E7 N/m), the rotation around the y-axis (k = 1E8 Nm/rad) and it’s relatively stiff (e.g., 1E12 N/m or Nm/rad) in the other directions.

The new eigenfrequencies for the system are: f1 = 6.3 Hz, f2 = 9.8 Hz, f3 = 27.9 Hz and f4 = 48.5 Hz.
The deflection of the system looks as follows:

As expected, the spring experiences a deflection in the x-direction of 0.1 m (2E6N/(2E7N/m)) and a rotation around the y-axis that brings the node 4 at the free-end close to a x-position of 0.95 m.

[ebranlard]

I've reviewed the code and it looks good to me. The r-test needs to be updated it seems

[ebranlard]

The changes to the SubDyn input files should be updated in the r-test repository.
And also documented here:
https://openfast.readthedocs.io/en/dev/source/user/api_change.html

And likely some small documentation changes are needed.

[RBergua]

I have updated the documentation for this new spring element. I guess the current r-test will run just modifying the SubDyn input file accordingly. This is adding, below the RIGID LINK PROPERTIES section, the new entries:

----------------------- SPRING ELEMENT PROPERTIES -------------------------------------
             0   NSpringPropSets - Number of spring properties
PropSetID   k11     k12     k13     k14     k15     k16     k22     k23     k24     k25     k26     k33     k34     k35     k36     k44      k45      k46      k55      k56      k66    
  (-)      (N/m)   (N/m)   (N/m)  (N/rad) (N/rad) (N/rad)  (N/m)   (N/m)  (N/rad) (N/rad) (N/rad)  (N/m)  (N/rad) (N/rad) (N/rad) (Nm/rad) (Nm/rad) (Nm/rad) (Nm/rad) (Nm/rad) (Nm/rad)

[deslaughter]

Found a few minor things. Overall it looks really good.

[deslaughter]

It's a little confusing that the order of the coefficients listed here doesn't match the order in which they're specified in the table.

[RBergua]

The description could be simplified by stating:
The table lists for each spring property: the property ID (**PropSetID**) and the stiffness coefficients (**K11**, **K12**, **K13**, **K14**, **K15**, **K16**, **K22**, **K23**, **K24**, **K25**, **K26**, **K33**, **K34**, **K35**, **K36**, **K44**, **K45**, **K46**, **K55**, **K56**, **K66**).

Like this it would appear in the same order as in the input file.

[andrew-platt]

See PR RBergua#2 for updates (brings this PR up to dev -- only necessary to get r-test to the right point).

[verlivkra]

Will the spring elements be included in the next release of OpenFAST? And when is that coming?

[RBergua]

The spring elements for SubDyn have been merged into OpenFAST dev branch and will be available in OpenFAST v.4.0.0. Such OpenFAST version should be released relatively soon.