diff --git a/docs/source/user/subdyn/theory.rst b/docs/source/user/subdyn/theory.rst index 36c3d932be..f401048a49 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,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 & \\ @@ -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: 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