From 62f6f9d44786c5dfda1610ab26535336522c6f4d Mon Sep 17 00:00:00 2001 From: Lu Wang Date: Wed, 30 Sep 2026 18:44:27 -0600 Subject: [PATCH 1/3] SubDyn: use linear thin-rod consistent mass for rigid links and cables Rigid links and cables previously reused the beam's cubic transverse mass coefficients (13/35, 9/70) while dropping all rotational-DOF inertia. For a rigid-body rotation this overstates the inertia about an offset point (13/35 m L^2 instead of the correct 1/3 m L^2, ~11.4%) and is kinematically inconsistent with the element's linear displacement field. Correct ElemM (renamed from ElemM_Cable to ElemM_Rod) to the linear consistent mass: isotropic t/3 on the translational diagonals and t/6 on the node-to-node coupling, zero rotational inertia. This reproduces a uniform thin rod's rigid-body inertia exactly (m L^2/3 about an end, m L^2/12 about the centre) and is shared by cables and rigid links. Update SubDyn theory doc mass matrix accordingly. Co-authored-by: GitHub Copilot Co-authored-by: Claude --- docs/source/user/subdyn/theory.rst | 32 +++++++++++++++++++++--------- modules/subdyn/src/FEM.f90 | 27 +++++++++++++++---------- modules/subdyn/src/SD_FEM.f90 | 6 +++--- 3 files changed, 42 insertions(+), 23 deletions(-) diff --git a/docs/source/user/subdyn/theory.rst b/docs/source/user/subdyn/theory.rst index 36c3d932be..3157dd483c 100755 --- a/docs/source/user/subdyn/theory.rst +++ b/docs/source/user/subdyn/theory.rst @@ -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 @@ -691,8 +691,9 @@ 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:: @@ -700,14 +701,14 @@ given by: \boldsymbol{M}_e = \rho L_e \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 & \\ @@ -719,6 +720,19 @@ given by: with :math:`L_e` the *undisplaced* length of the element (not :math:`L_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. This element mass +reproduces the rigid-body inertia of a uniform thin rod exactly: a total +mass :math:`m=\rho L_e`, 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 used for +*rigid-link* members, which share this straight two-node formulation. A +purely collinear assembly therefore carries no inertia about its own +axis. + .. _SD_ControlCable: Controlled pretension cable diff --git a/modules/subdyn/src/FEM.f90 b/modules/subdyn/src/FEM.f90 index f645c10c5c..780105e963 100644 --- a/modules/subdyn/src/FEM.f90 +++ b/modules/subdyn/src/FEM.f90 @@ -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 @@ -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 @@ -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: diff --git a/modules/subdyn/src/SD_FEM.f90 b/modules/subdyn/src/SD_FEM.f90 index 2b56d40019..6e83270c3d 100644 --- a/modules/subdyn/src/SD_FEM.f90 +++ b/modules/subdyn/src/SD_FEM.f90 @@ -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 From 4628ce79c85be16e7acd08e1cae18957837de96a Mon Sep 17 00:00:00 2001 From: Lu Wang Date: Thu, 1 Oct 2026 17:33:13 -0600 Subject: [PATCH 2/3] Update SubDyn user docs to clarify that cable mass is computed based on the unstretched length L_0, not the nominal length L_e --- docs/source/user/subdyn/theory.rst | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/docs/source/user/subdyn/theory.rst b/docs/source/user/subdyn/theory.rst index 3157dd483c..404aa89f6c 100755 --- a/docs/source/user/subdyn/theory.rst +++ b/docs/source/user/subdyn/theory.rst @@ -698,7 +698,7 @@ 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} 1/3 & 0 & 0 & & & & 1/6 & 0 & 0 & & & \\ @@ -717,8 +717,9 @@ from the same linear interpolation of the displacement field used above: \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)`. (Note that the above mass matrix is +also used for rigid-link members, in which case, :math:`L_0=L_e`.) Because the displacement varies linearly along the element (the cable is never subdivided and carries no bending shape), the translational mass is From 8f7b15c0925dbb23a56a1d8c199dc0d999fca7b8 Mon Sep 17 00:00:00 2001 From: Lu Wang Date: Thu, 1 Oct 2026 17:52:22 -0600 Subject: [PATCH 3/3] Update SubDyn user docs to clarify that cable mass is computed based on the unstretched length L_0, not the nominal length L_e --- docs/source/user/subdyn/theory.rst | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/docs/source/user/subdyn/theory.rst b/docs/source/user/subdyn/theory.rst index 404aa89f6c..f401048a49 100755 --- a/docs/source/user/subdyn/theory.rst +++ b/docs/source/user/subdyn/theory.rst @@ -718,21 +718,20 @@ from the same linear interpolation of the displacement field used above: \label{eq:MassMatrixPreTension}\end{aligned} with :math:`L_0` the *rest* (unstretched) length of the cable, -:math:`L_0=L_e/(1+\epsilon_0)`. (Note that the above mass matrix is -also used for rigid-link members, in which case, :math:`L_0=L_e`.) +: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. This element mass +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_e`, 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 used for -*rigid-link* members, which share this straight two-node formulation. A -purely collinear assembly therefore carries no inertia about its own -axis. +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: