Multibody Dynamics Using Quasi Energy and Momentum Conservative Algorithm Implemented with a Novel Four-Node Co-Rotational Quadrilateral Shell Element

1  Introduction

The co-rotational framework, first introduced in [1–3], is widely recognized as an efficient strategy for geometrically nonlinear analysis, particularly for problems involving large displacements and rotations [4–9]. A significant advancement in this field was the Element-Independent Co-rotational (EICR) method [10,11], which facilitates the formulation of new co-rotational elements while preserving the vector transformation matrix. This advancement broadens the applicability of the approach to various element types. Crisfield and Moita [12] incorporated the concept of enhanced strain into the traditional co-rotational quadrilateral element, applying it to the analysis of two-dimensional solid beams, and subsequently extended the approach to three-dimensional solid elements [13–15]. Masud et al. [16,17] introduced an eight-node hexahedral element based on the EICR method, with the objective of modeling the bending-dominated geometrically nonlinear behavior of multilayer shells. Utilizing the natural mode method, Argyris et al. [18] developed a triangular shell element, in which the geometric stiffness matrix is derived by introducing disturbances into the equilibrium equations. Kan et al. [19] employed a streamlined two-step methodology: it first simplifies internal force calculations by omitting projector matrix terms, and then enforces element equilibrium via a global-frame correction, thereby eliminating complex local-global rotation coupling while providing a unified, simplified framework for diverse element types.

Despite its demonstrated element generality and high accuracy [20–24], the practical application of the co-rotational formulation is hindered by prohibitive computational costs, particularly in nonlinear transient analysis. This limitation primarily stems from the asymmetric tangent stiffness matrix inherent to most co-rotational formulations. Since conventional rotational degrees of freedom (DOF) are non-additive and non-commutative, the tangent stiffness matrix cannot be derived directly from the second-order partial derivatives of the strain energy functional with respect to nodal variables. This asymmetry not only increases the complexity of numerical computation but also consumes more computer resources, significantly reducing solution efficiency and thus limiting its engineering application to large-scale or complex dynamic problems.

To address the aforementioned issues, Li and Izzuddin [25,26] introduced vectorial rotational variables as nodal rotational DOFs, developing a series of novel co-rotational beam elements [27] and curved shell elements [28–32]. In computing the second-order partial derivatives of the strain energy functional, the differentiation order with respect to any nodal variables is commutative. This commutativity ensures that the derived element’s tangent stiffness matrix exhibits symmetry under both local and global coordinate systems, thereby improving computational performance and reducing storage needs. Furthermore, during the incremental solution procedure, the increments of the vectorial rotational variables can be updated by simple addition, eliminating the need to construct rotation matrices for vector updates and leading to more efficient computation. With the introduction of vectorial rotational variables, the potential of the co-rotational method for high-fidelity simulations and advanced engineering applications is further unlocked, promising broader adoption in the future.

Traditional time integration algorithms, particularly the Newmark family of methods [33–35], while unconditionally stable for linear problems, often exhibit energy drift and numerical instability in nonlinear regimes [36,37]. This is primarily due to significant geometric nonlinearities and strong coupling between rigid body motion and elastic deformation in such systems. These two difficulties challenge the linear assumptions within time steps and lead to instability and/or inaccuracies in long-duration simulations [38,39].

To address convergence issues in nonlinear dynamic analysis, Belytschko and Schoeberle [40] proposed using energy conservation as a stability criterion for time integration algorithms. Hughes et al. [41] introduced a constrained energy method, incorporating Lagrange multipliers into the trapezoidal integration scheme. Building on this work, Kuhl and Ramm [42] added momentum conservation as an additional constraint. However, this approach of enforcing energy conservation via Lagrange multipliers still suffers from convergence difficulties during computation. In contrast, the Energy–Momentum Conservation (EMC) algorithm proposed by Simo [36,43–45] exhibits unconditional stability for both linear and nonlinear dynamic problems. This method was later extended to beams [46–48], shells [49–52], and multibody systems [53–55], demonstrating remarkable stability even in long-duration simulations.

Nonlinear time integration methods that preserve energy conservation continue to evolve. Magisano et al. [49] improved upon Simo and Tarnow’s approach by redefining the midpoint internal force vector using (1) averaged endpoint stresses and (2) a time-averaged strain-displacement tangent operator. By retaining the Hamiltonian structure present in linear elastodynamic equations, Sánchez et al. [56] developed a symplectic Hamiltonian finite element method capable of conserving both energy and momentum. Schiebl and Romero [57] addressed the discretization of energy and momentum conserving time integration algorithms in atomic particle dynamics, with a specific focus on three key aspects: handling periodic boundary conditions, approximating three-body interactions, and extending the framework to functional potentials. Zhang et al. [58–60] enforced energy conservation as a constraint and enhanced the accuracy of the formulation by selecting three or more interpolation points for computing internal forces. To address the fundamental challenge of defining a reference frame in the deformed configuration for 3D Euler–Bernoulli beam theory, Santana et al. [61] proposed two innovative approaches based on Gram–Schmidt orthogonalization and displacement derivatives. These approaches were subsequently extended by Chhang et al. [62] to dynamic problems, yielding a stable time integration algorithm which strictly preserves both energy and momentum.

The solution scheme combining the co-rotational framework with the energy conserving algorithm can stably solve dynamic problems involving large displacements and rotations. Crisfield and Shi [20], by modifying the equations of motion and introducing correction coefficients, developed an approximately energy-conserving co-rotational algorithm for beam and rod structures, which significantly improves computational efficiency. Zhong and Crisfield [63] extended this algorithm to shell elements. Based on the co-rotational formulation, Almeida and Awruch [64] proposed an implicit time integration algorithm, which approximately conserves system energy and was successfully applied to analyze the nonlinear dynamic behavior of laminated composite shells undergoing large rotations. Yang and Xia [65,66] innovatively adopted a predictor-corrector iterative strategy combined with the generalized energy-momentum method, achieving efficient and stable solutions for co-rotational shell elements.

This paper proposes a new computational theory for multibody dynamics by employing a novel co-rotational formulation of the four-node quadrilateral shell element [29,31] and introducing a quasi-energy–momentum conserving algorithm [67]. In the proposed co-rotational formulation, the rotational variables at each node are represented by vector components. Because the sequence of differentiation is commutative when computing the second-order partial derivatives of the energy functional with respect to the nodal variables, the tangent stiffness matrix of the element exhibits symmetry. Assumed membrane and shear strains [29,68–70] are employed to mitigate the computational inaccuracies and inefficiencies arising from membrane and shear locking. Different from the exact energy–momentum conserving (EMC) algorithm [44], which constructs the strain energy functional using both strain and stress quantities, and replaces the mid-step strain values with the arithmetic mean of the endpoint strains while evaluates the stresses at the midpoint, our quasi-energy–momentum conserving algorithm directly evaluates the strain energy functional using the strain vectors. In our approach, all strain vectors, as well as their first and second partial derivatives with respect to the nodal variables, are computed via the arithmetic mean of the corresponding endpoint values. Furthermore, after discretizing the dynamic equilibrium differential equations, the term containing the first-order time derivatives of the rotational variables is incorporated into the equivalent load vector. This yields a symmetric mass matrix, consequently improving computational speed and reducing storage requirements. The proposed methodology ensures the nearly exact conservation of total energy, linear momentum, and angular momentum, while demonstrating excellent numerical stability throughout the computation.

The structure of this paper is organized as follows: Section 2 presents the co-rotational framework of the four-node quadrilateral shell element, where the element formulation is derived under local coordinate system. Section 3 establishes the dynamic equilibrium differential equations and employs a generalized midpoint interpolation method for incremental iterative solution. Section 4 validates the accuracy and stability of the developed multibody dynamics approach through four classic shell-structure examples. Section 5 provides the conclusions.

2  Co-Rotational Formulation

In Fig. 1, d4 and t4 denote the translational displacement vectors of Node 4 in the global (O-X-Y-Z) and local (o-x-y-z) coordinate system, respectively. p10 and p1 represent the normal vectors of the mid-surface at Node 1 under the initial and deformed configurations in the global coordinate system, while r10 and r1 are the corresponding mid-surface normal vectors under the initial and deformed configuration within the local coordinate system. U, V, and W represent the displacement components along the X, Y, and Z directions of the global coordinate system, respectively. The element’s transition from the initial to the current configuration consists of two stages. First, it rotates as a rigid body to an intermediate configuration without deformation. Then, it proceeds from this intermediate configuration to the current configuration through pure deformation, during which rigid-body translation also occurs.

Figure 1: Co-rotational framework of the element.

For the element in its initial configuration, the local coordinate system’s directional vectors are defined by the following equations:

{ex0=e130−e240|e130−e240|ey0=e130+e240|e130+e240|ez0=ex0×ey0(1a,b,c)

where the vectors e130 and e240, oriented from Node 1 toward Node 3 and from Node 2 toward Node 4, respectively, are the normalized direction vectors along the quadrilateral’s diagonals. They are calculated from the vectors v130 and v240 as

{e130=v130|v130|=X30−X10|X30−X10|e240=v240|v240|=X40−X20|X40−X20|(2a,b)

where Xi0(i=1,2,3,4) denotes the global coordinates of Node i in the initial configuration.

Similarly, for the deformed configuration, these directional vectors of the local coordinate system can be obtained from the normalized vectors e13 and e24 using:

{ex=e13−e24|e13−e24|ey=e13+e24|e13+e24|ez=ex×ey(3a,b,c)

where e13 and e24 are calculated as follows:

{e13=v13|v13|=(X30+d3)−(X10+d1)|(X30+d3)−(X10+d1)|e24=v24|v24|=(X40+d4)−(X20+d2)|(X40+d4)−(X20+d2)|(4a,b)

In the global coordinate system, this four-node quadrilateral shell element has five or six DOFs per node—three for translation and two or three for rotation. The element’s nodal variable vector is:

uGT=⟨d1Tng1Td2Tng2Td3Tng3Td4Tng4T⟩(5)

where diT=⟨UiViWi⟩ represents the three translational DOFs of Node i. For nodes of a smooth shell, or nodes away from the intersection line of a non-smooth shell, ngiT=⟨pi,npi,m⟩ represents the two rotational DOFs of Node i, pi,n and pi,m denote the two components of the mid-surface normal vector pi at Node i with smaller absolute values. For nodes located on the intersection lines of a non-smooth shell, ngiT=⟨eiy,neiy,meiz,n⟩, eiy,n, eiy,m and eiz,n denote the three components with smaller absolute values of the nodal orientation vectors eiy and eiz at Node i [32].

In the local coordinate system, the four-node quadrilateral shell element has five DOFs per node, resulting in a total of 20 nodal variables:

uLT=⟨t1Tθ1Tt2Tθ2Tt3Tθ3Tt4Tθ4T⟩(6)

where tiT=⟨uiviwi⟩ represents the three translational DOFs of Node i, and θiT=⟨rixriy⟩ represents the two rotational DOFs of Node i, rix and riydenote the components along the x and y directions of the mid-surface normal vector ri at Node i, respectively.

In the initial configuration, the nodal coordinates xi0 in the local coordinate system can be calculated from the global-coordinate ones as:

xi0={xi0yi0zi0}=R0vi0(7)

where R0T=[ex0ey0ez0] is the directional matrix composed of the normalized vectors along the three axes of the local coordinate system in the initial configuration; vi0=Xi0−X10.

The nodal translational DOFs in the local coordinate system are obtained from the global-coordinate ones as:

ti=xi−xi0=R(di+vi0)−R0vi0(8)

where RT=[exeyez] is the directional matrix of the local coordinate system in the current configuration, and xi represents the nodal coordinates in the local coordinate system at the current configuration.

For nodes of a smooth shell, or nodes of a non-smooth shell located away from the intersection, the relationship between the local and global components of the shell’s mid-surface normal vector is expressed as:

ri0=R0pi0(9)

ri=Rpi(10)

where ri0 and ri denote the mid-surface normal vectors at Node i in the local coordinate system for the initial and current configurations, respectively; pi0 and pi represent their counterparts in the global coordinate system.

For nodes on the intersection line of non-smooth shells, three global rotational DOFs are employed per node. The relationship between the local and global components of the shell’s mid-surface normal vector is then given by:

ri=RRiTRi0pi0(11)

where, Ri0T=[eix0eiy0eiz0] and RiT=[eixeiyeiz] represent the orientation matrices of Node i under the initial and current configurations, respectively.

At a node on the intersection line of a non-smooth shell, the triad of orientational vectors eix0, eiy0 and eiz0 in the initial configuration coincides with the direction vectors of the three axes of the global coordinate system. The initial orientation matrix Ri0 formed by these three vectors, is then given by:

Ri0=[eix0Teiy0Teiz0T]=[100010001](12)

In the global coordinate system, the initial value of the nodal mid-surface normal vector is calculated from the cross product of the tangent vectors along the natural coordinate ξ-axis and η-axis at the node:

p¯i0=∂[∑j=14Nj(ξ,η)Xj0]∂ξ×∂[∑j=14Nj(ξ,η)Xj0]∂η|(ξi,ηi)(13)

where “×” is the symbol of cross product, (ξi,ηi) are the natural coordinates corresponding to Node i, and Ni(ξ,η) is the shape function for Node i of the four-node quadrilateral shell element.

For a node shared by adjacent elements in a smooth shell, or by elements in the same patch of a non-smooth shell, the mid-surface normal is the average of the element normals at that node:

pi0=∑e(p¯i0/|p¯i0|)|∑e(p¯i0/|p¯i0|)|(14)

where e represents all elements containing node i.

In the local coordinate system, the quadrilateral element employs the Green-Lagrange strain measure suitable for shallow shells, which is split into the membrane strain εm, the bending strain z1χ, and the out-of-plane shear strain γ. The element’s potential energy functional is given by:

Π=12∫V[(εm+z1χ)TD1(εm+z1χ)+γTD2γ]dV−Wext(15)

where V is the element’s volume; Wext=fextTuL denotes the work done by external forces; fext represents the external force vector of the element; uL denotes the nodal variables of the element in the local coordinate system; D1 and D2 are the elasticity matrices.

Taking the stationarity of the potential energy, the element’s internal force vector in the local coordinate system is obtained:

fint=fext=∫V[(Bm+z1Bb)TD1(εm+z1χ)+BγTD2γ]dV(16)

The element’s tangent stiffness matrix is obtained by computing the partial derivative of the internal force vector with respect to the nodal variables uL in the local coordinate system:

kT=∫V[BmTD1Bm+(z1)2BbTD1Bb+∂BmT∂uLTD1εm+BγTD2Bγ]dV(17)

where Bm, z1Bb and Bγ are respectively the first-order partial derivatives of the membrane strain εm, bending strain z1χ, and shear strain γ with respect to the nodal variables uL in the local coordinate system (see Eqs. (A.1a–2b, A.4a–c, A.7a,b) on Pages 202–203 of the reference [29]); z1=ζh2, and ζ is the natural coordinate along the shell thickness; ∂BmT∂uLT denotes the second-order partial derivative of the membrane strain with respect to the nodal variables uL (see Eq. (A.8) on Page 203 of the reference [29]).

It can be seen from Eqs. (15)–(17) that the element tangent stiffness matrix kT is essentially obtained by calculating the second partial derivatives of the element strain energy functional (a scalar) with respect to the element nodal vector uL in the local coordinate system. Since the order of differentiation with respect to the nodal variables is commutative when computing these second partial derivatives, the resulting element tangent stiffness matrix kT is necessarily symmetric.

3  Computational Theory of Multibody Dynamics for Four-Node Quadrilateral Shell Element

For a conservative system, the Hamiltonian H, which represents the total mechanical energy, is defined as the summation of the kinetic energy K and the total potential energy V, and the total potential energy V comprises the internal strain energy Vint and the potential energy Vext arising from conservative external loads:

H=K+Vint−Vext=constant(18)

Differentiating the Hamiltonian H with respect to time gives:

dHdt=dKdt+dVintdt−dVextdt=0(19)

In the global coordinate system, the coordinates of any point within the element can be calculated by:

X(ξ,η,ζ,t)=X0(ξ,η,0)+d(ξ,η,t)+h2ζp(ξ,η,t)(20)

where, X0(ξ,η,0) represents the coordinates of the projection point of this point onto the shell’s mid-surface in the initial configuration. The variable d(ξ,η,t) denotes the translational displacement of this projection point, while p(ξ,η,t) represents the normal vector to the mid-surface at the same location. The symbol h indicates the thickness of the shell.

The velocity is computed by taking the time derivative of the coordinates:

X˙(ξ,η,ζ,t)=d˙(ξ,η,t)+12ζhp˙(ξ,η,t)=∑i=14Ni(d˙i+12ζhp˙i)(21)

where d˙i and p˙i represent the first-order derivatives of di and pi with respect to time, respectively.

Integrating the nodal kinetic energy over the element volume yields the kinetic energy functional of the four-node quadrilateral shell element:

Ke=12ρ∫V[∑i=14Ni(d˙i+ζhp˙i2)]T[∑j=14Nj(d˙j+ζhp˙j2)]dV(22)

The time derivative of the element kinetic energy functional is calculated as:

dKedt=ρ∫V[(∑i=14Nid˙i)T(∑j=14Njd¨j)+ζ2h24(∑i=14Nip˙i)T(∑j=14Njp¨j)]dV(23)

For Node i away from intersections of shells, the first derivative of the normal vector p˙i with respect to time is:

p˙i=∂pi∂ngiTn˙gi=[∂pi,X∂pi,n∂pi,X∂pi,m∂pi,Y∂pi,n∂pi,Y∂pi,m∂pi,Z∂pi,n∂pi,Z∂pi,m]{p˙i,np˙i,m}(24)

where m≠n≠l, and m,n,l∈{X,Y,Z}.

The second derivative of the normal vector p¨i with respect to time is:

p¨i=∂pi∂ngiTn¨gi+∂2pi∂ngiT2n˙gi⊗n˙gi(25)

where:

n¨gi={p¨i,np¨i,m}(26)

Define the matrix A1i and the vector A2i:

A1i=∂pi∂ngiT=[∂pi,X∂pi,n∂pi,X∂pi,m∂pi,Y∂pi,n∂pi,Y∂pi,m∂pi,Z∂pi,n∂pi,Z∂pi,m](27)

A2i=∂2pi∂ngiT2n˙gi⊗n˙gi={∂2pi,X∂pi,n2p˙2i,n+2∂2pi,X∂pi,n∂pi,mp˙i,np˙i,m+∂2pi,X∂pi,m2p˙2i,m∂2pi,Y∂pi,n2p˙2i,n+2∂2pi,Y∂pi,n∂pi,mp˙i,np˙i,m+∂2pi,Y∂pi,m2p˙2i,m∂2pi,Z∂pi,n2p˙2i,n+2∂2pi,Z∂pi,n∂pi,mp˙i,np˙i,m+∂2pi,Z∂pi,m2p˙2i,m}(28)

where:

{∂2pi,n∂pi,n2=∂2pi,n∂pi,m2=∂2pi,n∂pi,n∂pi,m=0∂2pi,m∂pi,n2=∂2pi,m∂pi,m2=∂2pi,m∂pi,n∂pi,m=0∂2pi,l∂pi,n2=−1pi,l−p2i,np3i,l∂2pi,l∂pi,m2=−1pi,l−p2i,mp3i,l∂2pi,l∂pi,n∂pi,m=−pi,npi,mp3i,l(29a–e)

Therefore, p˙i and p¨i, can be expressed as:

{p˙i=A1in˙gip¨i=A1in¨gi+A2i(30a,b)

For Node i located on an intersection line of non-smooth shells,

p˙i=(∂RiTRi0pi0∂ngiT)n˙gi=(∂[eixeiyeiz]∂ngiT[100010001]pi0){e˙iy,ne˙iy,me˙iz,n}(31)

p¨i=∂pi∂ng⁡iT⁢n¨g⁡i+∂2pi∂ng⁡iT⁢2⁢n˙g⁡i⊗n˙g⁡i(32)

Write Eqs. (31) and (32) in the form of Eq. (30a,b); then we have:

A1i=∂pi∂ngiT=∂[eixeiyeiz]∂ngiTpi0(33)

A2i=∂2pi∂ngiT2n˙gi⊗n˙gi=([∂2eix∂ngiT∂ngiT∂2eiy∂ngiT∂ngiT∂2eiz∂ngiT∂ngiT]pi0)n˙gi⊗n˙gi(34)

where:

{∂eiy∂ngiT=[∂eiy∂eiy,n∂eiy∂eiy,m∂eiy∂eiz,n]∂eiz∂ngiT=[∂eiz∂eiy,n∂eiz∂eiy,m∂eiz∂eiz,n]∂eix∂ngiT=eiy×∂eiz∂ngiT+∂eiy∂ngiT×eiz(35a,b,c)

∂2eix∂ngiT∂ngiT=eiy×∂2eiz∂ngiT∂ngiT+∂2eiy∂ngiT∂ngiT×eiz+∂eiy∂ngiT×∂eiz∂ngiT+∂eiy∂ngiT×∂eiz∂ngiT(36)

The expressions for the first-order partial derivatives in Eqs. (35) and (36) are given by Eqs. (B20a–k) and (B21a–c) on Pages 596–597 of the reference [32], while those for the second-order partial derivatives are given by Eq. (B47a–r) on Page 600 of the reference [32].

Substitute the time derivatives of the shell’s mid-surface normal vector at Node i into Eq. (23) and define the term fm containing A2j:

fm=ρh312{0∫−11∫−11N1∑j=14NjA11TA2jJadξdη⋮0∫−11∫−11N4∑j=14NjA14TA2jJadξdη}(37)

where Ja represents the Jacobian determinant,

Ja=|∂x∂ξ∂y∂ξ∂x∂η∂y∂η|(38)

After rearrangement, the following expression is obtained:

dKedt={d˙1n˙g1⋮d˙4n˙g4}T[C110…C1400Q11…0Q14⋮⋮⋱⋮⋮C410…C4400Q41…0Q44]{d¨1n¨g1⋮d¨4n¨g4}+{d˙1n˙g1⋮d˙4n˙g4}Tfm(39)

where:

{Cij=ρh∫−11∫−11NiNjJadξdη[100010001]Qij=ρh312∫−11∫−11NiNjA1iTA1jJadξdη(40a,b)

The definition of the element mass matrix m is as follows:

m=[C110…C1400Q11…0Q14⋮⋮⋱⋮⋮C410…C4400Q41…0Q44](41)

Thus, the first-order time derivative of the element kinetic energy is given by:

dKedt=u˙GTmu¨G+u˙GTfm(42)

The mass matrices m of all four-node quadrilateral elements in the system are assembled into the global mass matrix M, the nodal variable vectors uG of all elements are assembled into the global nodal variable vector U, and fm are assembled into Fm. The first-order time derivative of system’s total kinetic energy is then obtained:

dKdt=∑ndKedt=U˙TMU¨+U˙TFm(43)

The time derivative of the element strain energy functional is calculated as:

dVintedt=u˙GT(∂uL∂uGT)Tfint==u˙GTTTfint=u˙GTfG(44)

where T denotes the transformation matrix from local to global coordinate system (see Eq. (39) on Page 580, and Eqs. (B1)–(B21c) on Pages 595–597 of the reference [32]), fG represents the element’s internal force vector in the global coordinate system, with

{T=∂uL∂uGTfG=TTfint(45a,b)

The total linear momentum Le and the total angular momentum Je (the reference point is at the origin of global coordinate system) of a shell element can be calculated as following,

Le=∫VρX˙dV=∫Vρ(d˙+12hζp˙)dV=∫Aρd˙hdA=∫Aρh∑i=14Nid˙idA(46)

Je=∫VX×ρX˙dV=∫AX¯×ρhd˙dA+112ρh3∫Ap×p˙dA=ρh∫A∑i=14[Ni(Xi0+di)×d˙i]dA+112ρh3∫A∑i=14(Nipi×p˙i)dA(47)

where X¯ is the global coordinate of a point of shell mid-surface at Time t, with

X¯=X¯0(ξ,η,0)+d(ξ,η,t)(48)

The total linear momentum and total angular momentum of the system are respectively expressed as:

L=∑eLe=∑eρh∫A∑i=14Nid˙idA(49)

J=∑eJe=∑eρh∫A[∑i=14Ni(Xi0+di)×d˙i+∑i=14h212Ni(pi×p˙i)]dA(50)

The time derivative of the system’s total strain energy is obtained by assembling the time derivatives of the individual elemental strain energies:

dVintdt=∑ndVintedt=∑n(u˙GTfG)=U˙TFint(51)

where n represents all elements in the system; Fint is the system’s internal force vector in the global coordinate system, assembled from fG,

Fint=∑nfG=∑n(TTfint)(52)

The system’s external work accumulated from tn to tn+1 is:

Vext=∫tntn+1U˙TFextdt(53)

The first-order time derivative of the system’s external work is:

dVextdt=U˙TFext(54)

Substituting the time derivatives of the system’s kinetic energy functional dKdt, strain energy fuctional dVintdt, and external work dVextdt into Eq. (19), and eliminating the U˙Tterm yields the system’s dynamic equilibrium differential equation:

MU¨+Fm+Fint−Fext=0(55)

The term Fm containing the first-order time derivatives of the rotational variables is then incorporated into the equivalent load vector.

For [tn,tn+1] with time step length Δt=tn+1−tn, based on the generalized energy-momentum method [71,72], appropriate generalized midpoints tn+1−αf and tn+1−αm are selected:

{tn+1−αf=(1−αf)tn+1+αftntn+1−αm=(1−αm)tn+1+αmtn(56a,b)

where αf and αm are integration parameters, controlled by the spectral radius ρ∞ (see Eq. (67a,b) in the reference [67]).

Based on the Newmark-β method, the nodal velocities U˙(tn+1−αf) and accelerations U¨(tn+1−αm) at the generalized midpoints are given by:

{U˙(tn+1−αf)=(1−αf)γβΔt[U(tn+1)−U(tn)]−γ−β−αfγβU˙(tn)−(γ−2β)(1−αf)2βU¨(tn)ΔtU¨(tn+1−αm)=(1−αm)βΔt2[U(tn+1)−U(tn)]−(1−αm)βΔtU˙(tn)−(1−2β−αm)2βU¨(tn)(57a,b)

where γ and β are integration parameters, which are also controlled by the spectral radius ρ∞ (see Eq. (70a,b) in the reference [67]).

At the generalized midpoints tn+1−αf and tn+1−αm, the system’s dynamic equilibrium differential Eq. (55) can be expressed as:

M(tn+1−αf)U¨(tn+1−αm)+Fm(tn+1−αf)+Fint(tn+1−αf)−Fext(tn+1−αf)=0(58)

Considering Eq. (57), the first term of Eq. (58) can be rewritten as:

M(tn+1−αf)U¨(tn+1−αm)=M(tn+1−αf)[(1−αm)βΔt2(U(tn+1)−U(tn))−(1−αm)βΔtU˙(tn)−(1−2β−αm)2βU¨(tn)](59)