Skip to content

SubDyn: use linear thin-rod consistent mass for rigid links and cables #3484

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
luwang00 wants to merge 3 commits into
base: rc-5.0.1
Choose a base branch
from
+45 −26
Open

SubDyn: use linear thin-rod consistent mass for rigid links and cables #3484

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
62f6f9d
SubDyn: use linear thin-rod consistent mass for rigid links and cables
luwang00 Oct 1, 2026
4628ce7
Update SubDyn user docs to clarify that cable mass is computed based …
luwang00 Oct 1, 2026
8f7b15c
Update SubDyn user docs to clarify that cable mass is computed based …
luwang00 Oct 1, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
38 changes: 26 additions & 12 deletions docs/source/user/subdyn/theory.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -633,9 +633,9 @@ Finite element formulation of a pretension cable

The rotational degrees of freedom are omitted for conciseness since
these degrees of freedom are not considered in this cable element. The
linear formulation from is applied to both nodes of a finite element,
interpreting the force at each node as the internal force that the
element exert on the nodes. Using this convention, the pretension cable
linear force--displacement relations derived above are applied to both
nodes of a finite element, interpreting the force at each node as the
internal force that the element exerts on the nodes. Using this convention, the pretension cable
element can be represented with an element stiffness matrix
:math:`\boldsymbol{K}_e` and an additional nodal load vector
:math:`\boldsymbol{f}_{e,0}` such that the static equilibrium equation
Expand Down Expand Up @@ -691,23 +691,24 @@ equations leads to the formulation of a truss element. The linear model
above is only valid for :math:`L_d-L_0>0`, that is
:math:`(L_e-L_0+u_{z,2}-u_{z,1})>0`, and the implementation should abort
if this condition is not reached at a given time. If the cable has a
positive mass density :math:`\rho`, the mass matrix of the element is
given by:
positive mass density :math:`\rho` (mass per unit length), the mass
matrix of the element is the *linear* consistent mass matrix, obtained
from the same linear interpolation of the displacement field used above:

.. math::

\begin{aligned}
\boldsymbol{M}_e = \rho L_e
\boldsymbol{M}_e = \rho L_0
\left[
\begin{array}{*{12}c}
13/35 & 0 & 0 & & & & 9/70 & 0 & 0 & & & \\
0 & 13/35 & 0 & & \boldsymbol{0}_3 & & 0 & 9/70 & 0 & & \boldsymbol{0}_3 & \\
1/3 & 0 & 0 & & & & 1/6 & 0 & 0 & & & \\
0 & 1/3 & 0 & & \boldsymbol{0}_3 & & 0 & 1/6 & 0 & & \boldsymbol{0}_3 & \\
0 & 0 & 1/3 & & & & 0 & 0 & 1/6 & & & \\
& & & & & & & & & & & \\
& \boldsymbol{0}_3 & & & \boldsymbol{0}_3 & & & \boldsymbol{0}_3 & & & \boldsymbol{0}_3 & \\
& & & & & & & & & & & \\
9/70 & 0 & 0 & & & & 13/35 & 0 & 0 & & & \\
0 & 9/70 & 0 & & \boldsymbol{0}_3 & & 0 & 13/35 & 0 & & \boldsymbol{0}_3 & \\
1/6 & 0 & 0 & & & & 1/3 & 0 & 0 & & & \\
0 & 1/6 & 0 & & \boldsymbol{0}_3 & & 0 & 1/3 & 0 & & \boldsymbol{0}_3 & \\
0 & 0 & 1/6 & & & & 0 & 0 & 1/3 & & & \\
& & & & & & & & & & & \\
& \boldsymbol{0}_3 & & & \boldsymbol{0}_3 & & & \boldsymbol{0}_3 & & & \boldsymbol{0}_3 & \\
Expand All @@ -716,8 +717,21 @@ given by:
\right]
\label{eq:MassMatrixPreTension}\end{aligned}

with :math:`L_e` the *undisplaced* length of the element (not
:math:`L_0`).
with :math:`L_0` the *rest* (unstretched) length of the cable,
:math:`L_0=L_e/(1+\epsilon_0)`.

Because the displacement varies linearly along the element (the cable is
never subdivided and carries no bending shape), the translational mass is
isotropic: the diagonal terms are :math:`1/3` and the terms coupling the
two nodes are :math:`1/6`, identical in the axial and transverse
directions. No inertia is associated with the rotational degrees of
freedom, consistent with the slender-rod idealization. A purely collinear
assembly therefore carries no inertia about its own axis. This element mass
reproduces the rigid-body inertia of a uniform thin rod exactly: a total
mass :math:`m=\rho L_0`, an inertia :math:`m L_e^2/3` about an end, and
:math:`m L_e^2/12` about its centre. The same mass matrix is also used for
*rigid-link* members, which share this straight two-node formulation. In
the case of rigid-link members, :math:`L_0=L_e`.

.. _SD_ControlCable:

Expand Down
27 changes: 16 additions & 11 deletions modules/subdyn/src/FEM.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1242,8 +1242,13 @@ SUBROUTINE ElemM_Beam(A, L, Ixx, Iyy, Jzz, rho, DirCos, M)

END SUBROUTINE ElemM_Beam
!------------------------------------------------------------------------------------------------------
!> Element stiffness matrix for pretension cable
SUBROUTINE ElemM_Cable(A, L, rho, DirCos, M)
!> Consistent mass for a straight 2-node bar/rod, used for both cables and rigid links.
!! Displacement is linearly interpolated in all three directions (a single, non-subdivided,
!! straight element has no bending shape), giving the linear consistent mass t/3 on the
!! translational diagonals and t/6 on the node-to-node coupling. This reproduces a uniform
!! thin rod's rigid-body inertia exactly (m*L^2/3 about an end, m*L^2/12 about the center).
!! Rotational DOF carry no inertia: a thin rod has no polar or section rotary inertia.
SUBROUTINE ElemM_Rod(A, L, rho, DirCos, M)
REAL(ReKi), INTENT( IN) :: A,rho
REAL(FEKi), INTENT( IN) :: L
REAL(FEKi), INTENT( IN) :: DirCos(3,3) !< From element to global: xg = DC.xe, Kg = DC.Ke.DC^t
Expand All @@ -1256,20 +1261,20 @@ SUBROUTINE ElemM_Cable(A, L, rho, DirCos, M)

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

M( 1, 1) = 13._FEKi/35._FEKi * t
M( 2, 2) = 13._FEKi/35._FEKi * t
M( 1, 1) = t/3.0_FEKi
M( 2, 2) = t/3.0_FEKi
M( 3, 3) = t/3.0_FEKi

M( 7, 7) = 13._FEKi/35._FEKi * t
M( 8, 8) = 13._FEKi/35._FEKi * t
M( 7, 7) = t/3.0_FEKi
M( 8, 8) = t/3.0_FEKi
M( 9, 9) = t/3.0_FEKi

M( 1, 7) = 9._FEKi/70._FEKi * t
M( 2, 8) = 9._FEKi/70._FEKi * t
M( 1, 7) = t/6.0_FEKi
M( 2, 8) = t/6.0_FEKi
M( 3, 9) = t/6.0_FEKi

M( 7, 1) = 9._FEKi/70._FEKi * t
M( 8, 2) = 9._FEKi/70._FEKi * t
M( 7, 1) = t/6.0_FEKi
M( 8, 2) = t/6.0_FEKi
M( 9, 3) = t/6.0_FEKi

DC = 0.0_FEKi
Expand All @@ -1279,7 +1284,7 @@ SUBROUTINE ElemM_Cable(A, L, rho, DirCos, M)
DC(10:12, 10:12) = DirCos

M = MATMUL( MATMUL(DC, M), TRANSPOSE(DC) ) ! TODO: change me if DirCos convention is transposed
END SUBROUTINE ElemM_Cable
END SUBROUTINE ElemM_Rod
!------------------------------------------------------------------------------------------------------
!> calculates the lumped forces and moments due to gravity on a given element:
!! the element has two nodes, with the loads for both elements stored in array F. Indexing of F is:
Expand Down
6 changes: 3 additions & 3 deletions modules/subdyn/src/SD_FEM.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2603,14 +2603,14 @@ SUBROUTINE ElemM(ep, Me)
else if (ep%eType==idMemberCable) then
Eps0 = ep%T0/(ep%YoungE*ep%Area)
L0 = ep%Length/(1+Eps0) ! "rest length" for which pretension would be 0
CALL ElemM_Cable(ep%Area, L0, ep%rho, ep%DirCos, Me)
CALL ElemM_Rod(ep%Area, L0, ep%rho, ep%DirCos, Me)

else if (ep%eType==idMemberRigid) then
if ( EqualRealNos(eP%rho, 0.0_ReKi) ) then
Me=0.0_FEKi
else
CALL ElemM_Cable(ep%Area, real(ep%Length,FEKi), ep%rho, ep%DirCos, Me)
!CALL ElemM_(A, L, rho, DirCos, Me)
! Straight 2-node bar/rod consistent mass (thin-rod rigid-body inertia); shared with cables
CALL ElemM_Rod(ep%Area, real(ep%Length,FEKi), ep%rho, ep%DirCos, Me)
endif

else if (ep%eType==idMemberSpring) then
Expand Down
Loading