Skip to content

Commit

Permalink
Merge pull request #1889 from RBergua/RBergua-SubDyn_spring_element
Browse files Browse the repository at this point in the history 
New spring element in SubDyn: 6 by 6 stiffness matrix
  • Loading branch information
@andrew-platt
andrew-platt authored Nov 29, 2023
2 parents ae56ab1 + 1b5ff78 commit f620731
Show file tree
Hide file tree
Showing 12 changed files with 13,808 additions and 13,166 deletions.
7 changes: 7 additions & 0 deletions docs/source/user/api_change.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@ The changes are tabulated according to the module input file, line number, and f
The line number corresponds to the resulting line number after all changes are implemented.
Thus, be sure to implement each in order so that subsequent line numbers are correct.


OpenFAST v3.5.1 to OpenFAST dev
----------------------------------

Expand All @@ -23,8 +24,14 @@ Module Line Flag Name Example
============================================= ==== ==================== ========================================================================================================================================================================================================
HydroDyn all Complete restructuring of input file
SeaState all New module (split from HydroDyn, so contains some inputs previously found in HydroDyn)
SubDyn 56\* ----------------------- SPRING ELEMENT PROPERTIES -------------------------------------
SubDyn 57\* NSpringPropSets 0 - Number of spring properties
SubDyn 58\* PropSetID k11 k12 k13 k14 k15 k16 k22 k23 k24 k25 k26 k33 k34 k35 k36 k44 k45 k46 k55 k56 k66
SubDyn 59\* (-) (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)
============================================= ==== ==================== ========================================================================================================================================================================================================

\*Exact line number depends on number of entries in various preceeding tables.



OpenFAST v3.5.0 to OpenFAST v3.5.1
Expand Down
6 changes: 5 additions & 1 deletion docs/source/user/subdyn/examples/OC4_Jacket_SD_Input.dat
View file
Original file line number Diff line number Diff line change
Expand Up @@ -111,7 +111,7 @@ IJointID ItfTDXss ItfTDYss ItfTDZss ItfRDXss ItfRDYss ItfRDZss
56 1 1 1 1 1 1
----------------------------------- MEMBERS -------------------------------------------
112 NMembers - Number of members (-)
MemberID MJointID1 MJointID2 MPropSetID1 MPropSetID2 MType COSMID ![MType={1:beam circ., 2:cable, 3:rigid, 4:beam arb.}. COMSID={-1:none}]
MemberID MJointID1 MJointID2 MPropSetID1 MPropSetID2 MType COSMID ![MType={1:beam circ., 2:cable, 3:rigid, 4:beam arb., 5:spring}. COMSID={-1:none}]
(-) (-) (-) (-) (-) (-) (-)
1 1 2 2 2 1 -1
2 2 3 2 2 1 -1
Expand Down Expand Up @@ -247,6 +247,10 @@ PropSetID EA MatDens T0 CtrlChannel
0 NRigidPropSets - Number of rigid link properties
PropSetID MatDens
(-) (kg/m)
----------------------- 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)
---------------------- MEMBER COSINE MATRICES COSM(i,j) -------------------------------
0 NCOSMs - Number of unique cosine matrices (i.e., of unique member alignments including principal axis rotations); ignored if NXPropSets=0 or 9999 in any element below
COSMID COSM11 COSM12 COSM13 COSM21 COSM22 COSM23 COSM31 COSM32 COSM33
Expand Down
47 changes: 27 additions & 20 deletions docs/source/user/subdyn/input_files.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -179,11 +179,10 @@ properties, the material properties are not allowed to change within a
single member.

Future releases will allow for members of different cross-sections,
i.e., noncircular members. For this reason, the input file has
(currently unused) sections dedicated to the identification of direction
cosines that in the future will allow the module to identify the correct
orientation of noncircular members. The current release only accepts
tubular (circular) members.
i.e., noncircular members. For this reason, the input file has sections
dedicated to the identification of direction cosines that in the future
will allow the module to identify the correct orientation of noncircular
members. The current release only accepts tubular (circular) members.

The file is organized into several functional sections. Each section
corresponds to an aspect of the SubDyn model and substructure.
Expand Down Expand Up @@ -419,17 +418,20 @@ MEMBER X-SECTION PROPERTY table (discussed next) for starting
cross-section properties and **MPropSetID2** specifies the identifier
for ending cross-section properties, allowing for tapered members.
The sixth column specify the member type **MType**.
A member is one of the three following types (see :numref:`SD_FEM`):
A member is one of the four following types (see :numref:`SD_FEM`):

- Beams (*MType=1*), Euler-Bernoulli (*FEMMod=1*) or Timoshenko (*FEMMod=3*)

- Pretension cables (*MType=2*)

- Rigid link (*MType=3*)

- Spring element (*MType=5*)

**COSMID** refers to the IDs of the members' cosine matrices for noncircular
members; the current release uses SubDyn's default direction cosine convention
if it's not present or when COSMID values are -1.
members and spring elements; the current release uses SubDyn's default direction cosine convention
if it's not present or when COSMID values are -1. Spring elements are defined between joints that
are coincident in the space and the direction cosine must be provided.


An example of member table is given below
Expand Down Expand Up @@ -525,22 +527,27 @@ An example of rigid link properties table is given below
(-) (kg/m)
12 7850.0
3 7000.0
Spring Properties
~~~~~~~~~~~~~~~~
Members that are specified as spring elements (**MType=5**),
have their properties defined in the spring element properties table.
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**). The stiffness matrix is considered symmetric and
includes diagonal (kii) and cross-coupling (kij) coefficients.
The FEM representation of the spring element is given in :numref:`SD_SpringElement`.

An example of spring properties table is given below:

.. code::
-------------------------- SPRING ELEMENT PROPERTIES ----------------------------
1 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)
2 2E7 0 0 0 0 0 1E12 0 0 0 0 1E12 0 0 0 1E12 0 0 1E8 0 1E12
Member Cosine Matrices COSM (i,j)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Expand Down
51 changes: 47 additions & 4 deletions docs/source/user/subdyn/theory.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -279,19 +279,21 @@ The following joints are supported:

- Ball joint (*JointType=4*)

A member is one of the three following types:
A member is one of the four following types:

- Beams (*MType=1*), Euler-Bernoulli (*FEMMod=1*) or Timoshenko (*FEMMod=3*)

- Pretension cables (*MType=2*)

- Rigid link (*MType=3*)

- Spring element (*MType=5*)

Beam members may be split into several elements to increase the accuracy of the model (using
the input parameter *NDiv*). Member of other types (rigid links and
pretension cables) are not split. In this document, the term *element*
the input parameter *NDiv*). Member of other types (rigid links, pretension cables and springs)
are not split. In this document, the term *element*
refers to: a sub-division of a beam member or a member of another type
than beam (rigid-link or pretension cable). The term *joints* refers to
than beam (rigid-link, pretension cable or spring). The term *joints* refers to
the points defining the extremities of the members. Some joints are
defined in the input file, while others arise from the subdivision of
beam members. The end points of an elements are called nodes and each
Expand Down Expand Up @@ -1287,9 +1289,50 @@ The constraint are applied after the full system has been assembled.



.. _SD_SpringElement:


Spring Elements
~~~~~~~~~~~

Do not confuse the spring member with the springs defined as
a boundary condition in land-based systems. The spring element
relates two joints by means of a 6 by 6 stiffness matrix that
is assumed symmetric (k_ij = k_ji).

.. math::
\begin{aligned}
K=
\begin{bmatrix}
k_{11} & k_{12} & k_{13} & k_{14} & k_{15} & k_{16} \\
k_{21} & k_{22} & k_{23} & k_{24} & k_{25} & k_{26} \\
k_{31} & k_{32} & k_{33} & k_{34} & k_{35} & k_{36} \\
k_{41} & k_{42} & k_{43} & k_{44} & k_{45} & k_{46} \\
k_{51} & k_{52} & k_{53} & k_{54} & k_{55} & k_{56} \\
k_{61} & k_{62} & k_{63} & k_{64} & k_{65} & k_{66} \\
\end{bmatrix}

The spring element does not have a mass associated. However, if desired, a lumped mass can be
defined at the joints.

Since each joint has 6 DOFs (3 translations and 3 rotations), mathematically, the
spring element has a 12 by 12 dimension.
.. math::
\begin{aligned}
K_e=
\begin{bmatrix}
k_{6x6} & -k_{6x6} \\
-k_{6x6} & k_{6x6} \\
\end{bmatrix}

The spring element must be defined between two coincident joints and the orientation has to be
provided by means of the direction cosine. This allows the assembly of the spring element in the
global system stiffness matrix.

The spring element can be connected to beams, kinematic joints (e.g., revolute joint, universal joint,
and spherical joint), the interface joint and rigid links. However, it cannot be connected to a cable.

.. _GenericCBReduction:

Expand Down
167 changes: 164 additions & 3 deletions modules/subdyn/src/FEM.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -951,8 +951,9 @@ END SUBROUTINE GetRigidTransformation
!!
!! bjj: note that this is the transpose of what is normally considered the Direction Cosine Matrix
!! in the FAST framework.
SUBROUTINE GetDirCos(P1, P2, DirCos, L_out, ErrStat, ErrMsg)
SUBROUTINE GetDirCos(P1, P2, eType, DirCos, L_out, ErrStat, ErrMsg)
REAL(ReKi) , INTENT(IN ) :: P1(3), P2(3) ! (x,y,z) global positions of two nodes making up an element
INTEGER(IntKi), INTENT(IN ) :: eType ! element type (1:beam circ., 2:cable, 3:rigid, 4:beam arb., 5:spring)
REAL(FEKi) , INTENT( OUT) :: DirCos(3, 3) ! calculated direction cosine matrix
REAL(ReKi) , INTENT( OUT) :: L_out ! length of element
INTEGER(IntKi), INTENT( OUT) :: ErrStat ! Error status of the operation
Expand All @@ -966,9 +967,16 @@ SUBROUTINE GetDirCos(P1, P2, DirCos, L_out, ErrStat, ErrMsg)
Dz=P2(3)-P1(3)
Dxy = sqrt( Dx**2 + Dy**2 )
L = sqrt( Dx**2 + Dy**2 + Dz**2)

! The spring element should have the same starting and ending location. P1 and P2 must be coincident (L must be 0).
IF ( .not. EqualRealNos(L, 0.0_FEKi) .and. eType == 5) THEN
ErrMsg = ' Spring(s) must be defined with the same starting and ending locations in the element.'
ErrStat = ErrID_Fatal
RETURN
ENDIF

IF ( EqualRealNos(L, 0.0_FEKi) ) THEN
ErrMsg = ' Same starting and ending location in the element.'
IF ( EqualRealNos(L, 0.0_FEKi) .and. eType/= 5) THEN
ErrMsg = ' Same starting and ending location in a beam, cable or rigid element.'
ErrStat = ErrID_Fatal
RETURN
ENDIF
Expand Down Expand Up @@ -1157,6 +1165,159 @@ SUBROUTINE ElemK_Cable(A, L, E, T0, DirCos, K)
K = MATMUL( MATMUL(DC, K), TRANSPOSE(DC) ) ! TODO: change me if DirCos convention is transposed
END SUBROUTINE ElemK_Cable
!------------------------------------------------------------------------------------------------------
!> Element stiffness matrix for spring
!! The spring element can include diagnal and cross-coupling positions.
!! Assuming that the stiffness is symmetric (21 stiffness coefficients). The stiffness matrix could also be non-symmetric, if desired.
SUBROUTINE ElemK_Spring(k11, k12, k13, k14, k15, k16, k22, k23, k24, k25, k26, k33, k34, k35, k36, k44, k45, k46, k55, k56, k66, DirCos, K)
REAL(ReKi), INTENT( IN) :: k11, k12, k13, k14, k15, k16, k22, k23, k24, k25, k26, k33, k34, k35, k36, k44, k45, k46, k55, k56, k66
REAL(FEKi), INTENT( IN) :: DirCos(3,3) !< From element to global: xg = DC.xe, Kg = DC.Ke.DC^t
REAL(FEKi), INTENT(OUT) :: K(12, 12)
! Local variables
REAL(FEKi) :: DC(12, 12)

K(1:12,1:12) = 0.0_FEKi

K( 1, 1) = k11
K( 1, 7) = -K(1,1)
K( 7, 1) = -K(1,1)
K( 7, 7) = K(1,1)

K( 1, 2) = k12
K( 1, 8) = -K(1,2)
K( 7, 2) = -K(1,2)
K( 7, 8) = K(1,2)

K( 1, 3) = k13
K( 1, 9) = -K(1,3)
K( 7, 3) = -K(1,3)
K( 7, 9) = K(1,3)

K( 1, 4) = k14
K( 1, 10) = -K(1,4)
K( 7, 4) = -K(1,4)
K( 7, 10) = K(1,4)

K( 1, 5) = k15
K( 1, 11) = -K(1,5)
K( 7, 5) = -K(1,5)
K( 7, 11) = K(1,5)

K( 1, 6) = k16
K( 1, 12) = -K(1,6)
K( 7, 6) = -K(1,6)
K( 7, 12) = K(1,6)

K( 2, 2) = k22
K( 2, 8) = -K(2,2)
K( 8, 2) = -K(2,2)
K( 8, 8) = K(2,2)

K( 2, 3) = k23
K( 2, 9) = -K(2,3)
K( 8, 3) = -K(2,3)
K( 8, 9) = K(2,3)

K( 2, 4) = k24
K( 2, 10) = -K(2,4)
K( 8, 4) = -K(2,4)
K( 8, 10) = K(2,4)

K( 2, 5) = k25
K( 2, 11) = -K(2,5)
K( 8, 5) = -K(2,5)
K( 8, 11) = K(2,5)

K( 2, 6) = k26
K( 2, 12) = -K(2,6)
K( 8, 6) = -K(2,6)
K( 8, 12) = K(2,6)

K( 3, 3) = k33
K( 3, 9) = -K(3,3)
K( 9, 3) = -K(3,3)
K( 9, 9) = K(3,3)

K( 3, 4) = k34
K( 3, 10) = -K(3,4)
K( 9, 4) = -K(3,4)
K( 9, 10) = K(3,4)

K( 3, 5) = k35
K( 3, 11) = -K(3,5)
K( 9, 5) = -K(3,5)
K( 9, 11) = K(3,5)

K( 3, 6) = k36
K( 3, 12) = -K(3,6)
K( 9, 6) = -K(3,6)
K( 9, 12) = K(3,6)

K( 4, 4) = k44
K( 4, 10) = -K(4,4)
K(10, 4) = -K(4,4)
K(10, 10) = K(4,4)

K( 4, 5) = k45
K( 4, 11) = -K(4,5)
K(10, 5) = -K(4,5)
K(10, 11) = K(4,5)

K( 4, 6) = k46
K( 4, 12) = -K(4,6)
K(10, 6) = -K(4,6)
K(10, 12) = K(4,6)

K( 5, 5) = k55
K( 5, 11) = -K(5,5)
K(11, 5) = -K(5,5)
K(11, 11) = K(5,5)

K( 5, 6) = k56
K( 5, 12) = -K(5,6)
K(11, 6) = -K(5,6)
K(11, 12) = K(5,6)

K( 6, 6) = k66
K( 6, 12) = -K(6,6)
K(12, 6) = -K(6,6)
K(12, 12) = K(6,6)

! Stiffness matrix symmetry:
K(2:6, 1) = K(1,2:6)
K(2:6, 7) = K(1,8:12)
K(8:12, 1) = K(7,2:6)
K(8:12, 7) = K(7,8:12)

K(3:6, 2) = K(2,3:6)
K(3:6, 8) = K(2,9:12)
K(9:12, 2) = K(8,3:6)
K(9:12, 8) = K(8,9:12)

K(4:6, 3) = K(3,4:6)
K(4:6, 9) = K(3,10:12)
K(10:12, 3) = K(9,4:6)
K(10:12, 9) = K(9,10:12)

K(5:6, 4) = K(4,5:6)
K(5:6, 10) = K(4,11:12)
K(11:12, 4) = K(10,5:6)
K(11:12, 10) = K(10,11:12)

K(6, 5) = K(5,6)
K(6, 11) = K(5,12)
K(12, 5) = K(11,6)
K(12, 11) = K(11,12)

DC = 0.0_FEKi
DC( 1: 3, 1: 3) = DirCos
DC( 4: 6, 4: 6) = DirCos
DC( 7: 9, 7: 9) = DirCos
DC(10:12, 10:12) = DirCos

K = MATMUL( MATMUL(DC, K), TRANSPOSE(DC) ) ! TODO: change me if DirCos convention is transposed

END SUBROUTINE ElemK_Spring
!------------------------------------------------------------------------------------------------------
!> Element mass matrix for classical beam elements
SUBROUTINE ElemM_Beam(A, L, Ixx, Iyy, Jzz, rho, DirCos, M)
REAL(ReKi), INTENT( IN) :: A, L, Ixx, Iyy, Jzz, rho
Expand Down
Loading

0 comments on commit f620731

Please sign in to comment.