Configurational forces in variable-length beams for flexible multibody dynamics

This paper addresses the problem of sliding beams or beams with sliding appendages, with special emphasis on situations where the axial motion of the beam or appendage is not prescribed a priori. Hamilton’s variational principle is used to derive the weak and strong forms of governing equations based on the systematic use of Reynolds’ transport theorem. The strong form of the governing equations involve the mechanical and configurational momentum equations, together with the proper boundary conditions. It is shown that the configurational momentum equations are linear combinations of their mechanical counterparts and hence, are redundant. A weak form of the same equations is also developed; the configurational and mechanical momentum equations become independent because they combine in an integral form the strong mechanical and configurational momentum equations with their respective natural boundary conditions. The domain- and boundary-based formulations, stemming from these two forms of the governing equations, respectively, are proposed, and numerical examples are presented to contrast their relative performances. The predictions of both formulations are found to be in good agreement with those obtained from an ABAQUS model using contact pairs. The domain-based formulation presents a higher convergence rate than the boundary-based formulation. Clearly, the proper treatment of the configurational forces impacts the accuracy of the model significantly.

as joints, concentrated loads, or rigid bodies move along beams or cables. For simplicity of the exposition, the term "sliding beams" will be used to refer to both problems: sliding beams and beams with sliding appendages.
The dynamic behavior of these structures, often treated as geometrically nonlinear beams or cables, has attracted the attention of numerous researchers in recent years. For simplicity, the axial motion of the beams is often prescribed a priori. In reality, this axial motion results from the dynamic behavior of the system and hence, must be determined as part of the solution process. The present paper focuses on this problem.
Hamilton's variational principle is used to derive the weak and strong forms of the equations of motion for sliding beams. Due to the presence of sliding motion, the action integral is now defined over a time-varying material domain and the process of taking variations of this action integral must be handled carefully. The variation of functionals defined over time-varying domains appears in many engineering applications and numerous researchers have tackled this problem. For instance, variation of the cost functional in free-endpoint optimal control problems leads to the transversality condition at the free end time [1], variation of the potential energy functional of a beam contacting a rigid obstacle yields an additional boundary conditions determining the contact region [2], and first-and second-order variations of functionals defined over time-varying material domains are central to shape optimization [3][4][5].
Theoretically, action integrals defined over time-varying material domains can be transformed into action integrals over time-invariant reference domains through proper coordinate transformations. The variation process then becomes straightforward and the results can be transformed back to time-varying material domains. Reynolds' transport theorem [6] expresses this procedure mathematically: it relates the time derivative of an integral over a given control volume to the time derivative of the integrand and boundary fluxes. Clearly, Reynolds' transport theorem still holds when replacing the time derivative with a variation.
In general, structural dynamics problems are investigated using the Lagrangian formulation. Typically, for beam problems, the "material" or "Lagrangian coordinate" is selected as the arc-length coordinate of a material particle of the beam's axial line in its undeformed configuration. This coordinate is used to track down the beam's motion as it deforms in space. For sliding beams, this collection of material coordinates is time varying: for instance, an elevator cable comprises a time-varying collection of material particles as portions of the cable are reeled in or out of the capstan.
In finite element implementations, it is convenient to mesh the structure with a constant number of elements to avoid the costly re-meshing of the entire problem at each time step. The implementation of this approach requires the introduction of a "mesh coordinate." For instance, it is expeditious to model an elevator cable with a fixed number of elements, although its length varies in time. Mesh coordinate ξ ∈ [ξ 0 , ξ f ] is introduced, where coordinates ξ 0 and ξ f identify the root and end points of the variable-length cable and if the cable is modeled with N elements of equal length, the length of each element is (ξ 0 − ξ f )/N . Mesh coordinate ξ varies along the length of the cable, but no longer identifies a specific material particle of the cable. Of course, a one-to-one map relates the material and mesh coordinates.
This approach, known as the Arbitrary Lagrangian-Eulerian (ALE) formulation, has been developed by numerous authors for geometrically nonlinear sliding beams [7][8][9][10][11][12][13][14][15][16]. As will be shown in this paper, configurational forces arise at the sliding boundaries of the structure, affecting the boundary conditions at these points. Failure to identify these configurational forces properly can result prescribing erroneous boundary conditions.
Configurational forces were first introduced by Eshelby [17] to describe the driving forces acting on defects or singularities in solids; these forces arise from variations of the material configuration. A number of authors have investigated configurational forces in beams and cables both theoretically [18][19][20][21][22] and experientially [23][24][25]. Recently, a novel ALE formulation for dynamics problems of geometrically exact sliding beams [26] was developed by the first author based on dual quaternion kinematics. Material motions of the rigid-sections are decomposed into mesh motions (with mesh coordinates held fixed) and convective terms that take into account the relative motion of the material and mesh coordinates. Integration by parts then yields both strong and weak forms of the mechanical and configurational momentum equations together with the natural boundary conditions for mechanical and configurational forces at the endpoints. Although the aforementioned process is clear, integration by parts for terms involving time and spatial derivatives is complex and error prone.
The first objective of this paper is to streamline the derivation of both mechanical and configurational momentum equations for geometrically exact sliding beams. It appears that the procedure previously developed by the author [26] is needlessly complex and an approach based on Reynolds' transport theorem yields the same results far more expeditiously. This approach echoes that of Singh and Hanna [27], who develop a variational framework to derive both mechanical and configurational momentum equations together with the associate boundary conditions for a variety of continua, including three-dimensional solids, Euler's elastica, and ideal fluids.
In their strong forms, the mechanical and configurational momentum equations are not independent of each other: the configurational momentum equations form a linear combination of their mechanical counterparts. In contrast, the configurational boundary conditions are independent of their mechanical counterparts. In their weak or integral forms, the mechanical and configurational momentum equations are combined with the natural boundary conditions for mechanical and configurational forces, respectively, rendering the configurational equations independent of the mechanical counterparts. The second objective of this paper is to show that sliding beam problems can be formulated in two alternative manners: (1) domain-based approaches that solve the weak forms of mechanical and configurational momentum equations and (2) boundary-based approaches that solve the weak form of mechanical momentum equations and the strong form of configurational boundary conditions. Furthermore, the accuracies of the two approaches will be contrasted. This paper is organized as follows. Section 2 describes the mathematical relationship between the material and mesh coordinates, leading to the presentation Reynolds' transport theorem, the workhorse of the subsequent developments. Section 3 describes the kinematics of geometrically exact beams based on dual quaternions. The strong and weak forms of governing equations are derived in Sect. 4 based on Hamilton's variational principle. Section 5 presents the finite element formulations of the domain-based and boundary-based approaches. A set of numerical examples are presented in Sect. 6. Conclusions are presented at the end of the paper.

Material and mesh coordinates
To model sliding structures with a constant mesh configuration within the framework of the finite element method, thereby avoiding the costly re-meshing of the entire structure at each time step, the concepts of material and mesh coordinates are introduced. The distinction between the two types of coordinates is key to the developments below and is illustrated in Fig. 1 that shows the motions of the material and mesh coordinates of a beam sliding along a rigid or flexible sleeve, labeled cases (a) and (b), respectively. For case (a), interface point P is a material point of the rigid sleeve. The location of material point P of the beam coincides with that of point P in the reference configuration. The material and mesh coordinates of P are denoted s b and ξ b = const, respectively, and hence, point P (ξ b ) can also be viewed as a mesh point of the beam. As the beam slides along the rigid sleeve and deforms, material point P(s b ) leaves the interface point, while mesh point P (ξ b ) remains at the interface. The corresponding material and mesh motions are shown in Fig. 1.
For case (b), the sleeve is flexible and two interface points, denoted P L and P R , located at the left end of the beam and at the right end of the sleeve, respectively, must to be considered. Points P L and P R are material points of the beam and sleeve, respectively. In the reference configuration, the location of material point P L of the sleeve coincides with that of point P L , and its material and mesh coordinates are denoted s s and ξ s = const, respectively. Similarly, the location of material point P R of the beam coincides with that of point P R , and its material and mesh coordinates are denoted s b and ξ b = const, respectively. As the beam slides, both beam and sleeve deform. The locations of material points P L (s s ) and P R (s b ) no longer coincide with those of interface mesh points P L (ξ s ) and P R (ξ b ), respectively, as shown in Fig. 1.
For case (a), the mesh point on the left end of the beam is a material point of the rigid sleeve and hence, mesh motion can be determined easily. In contrast, for case (b), the determination of the mesh motion in the beam and sleeve is more involved because both motions depend on the deformed configurations of both beam and sleeve.
The configuration of the system depicted in Fig. 1 was simplified intentionally to ease the above discussion. More realistically, it could represent a hydraulic actuator, where the "sleeve" is the hydraulic cylinder tube and the "beam" is the piston/piston rod assembly. With the present formulation, both cylinder tube and piston rod can be modeled as interacting flexible components. The pressurized hydraulic fluid in the cylinder tube then applies equal and opposite axial forces to the ends of cylinder tube and piston rod. The modeling of the hydraulic components is beyond the scope of this paper. Of course, the hydraulic actuator could be part of a more complex multibody system.

Relating material and mesh coordinates
For a more rigorous description of the kinematics, let s ∈ [s 0 , s f ] denote the material coordinate of the beam's reference axis, where s 0 and s f denote the material coordinates of the endpoints. Let ξ ∈ [ξ 0 , ξ f ] denote the mesh coordinate. The material and mesh coordinates are related by a one-to-one map Consider a kinematic quantity, k , that is a function of space and time t . Mapping (1) now implies k (s, t) = k (s(ξ, t), t), i.e., the spatial dependency can be expressed through the material or mesh coordinates interchangeably. Let L denote a function of space, time, and kinematic quantity k , i.e., L = L(k (s, t), s, t). Sectional velocities and curvature are typical kinematic quantities, whereas a Lagrangian or a Hamiltonian are typical functions of space, time, and kinematic quantities.
The material and mesh time derivatives of function L, denotedL andL, respectively, represent derivatives with respect to time t while the material and mesh coordinates are held fixed, respectively,L where chain rule is used to obtain the indicated results. Clearly, the material and mesh time derivatives are related by the chain ruleL where L ,s = dL/ds = ∂L/∂s + (∂L/∂k )(∂k /∂s). Clearly, the material time derivative of material coordinate s vanishes:ṡ =s − s ,ss = 0. Likewise, the material and mesh variations, denoted δ andδ, respectively, are introduced to represent variations of quantities while the material and mesh coordinates are kept fixed, respectively, The inverse of the Jacobian is j −1 = ∂ξ/∂s. Taking a mesh time derivative and mesh variation of identity jj −1 = 1 yieldsj

Reynolds' transport theorem
Consider next the mesh time derivative of the following integral d dt In the first line, a change of variables is performed from the material to the mesh coordinate system. Because the mesh coordinate is time independent, the time derivative now commutes with the integral over the mesh domain. In the second line, time derivative d(j L)/dt is expanded with the help of identity (2b). Identity (6a) then yields the third line. Taking a variation of the same integral and using the same identities leads to a similar result Equations (8) and (9) are statements of Reynolds' transport theorem [6] that are well adapted to the developments to appear below.

Kinematics of geometrically exact beams
Geometrically exact beams are modeled as one-dimensional Cosserat continua. The kinematics of the beam is described by the motion of its cross-sections, assumed to remain rigid, at all points along its reference axis. Figure 2 depicts the configuration of a beam, where frames define the cross-sections in their reference and deformed configurations, respectively. Points B 0 and B are located at the intersection of the beam's reference axis with planar cross-sections in the two configurations, respectively. The plane of the cross-section is defined by two mutually orthogonal unit vectors, b β0 and b β , β = 2, 3, in the reference and deformed configurations, respectively. The sectional motion is defined by unit dual-quaternion [28] In the reference configuration, the corresponding motion is denoted q 0 (s, 0).
In the Lagrangian formulation, the kinematic quantities defined in the previous paragraph are functions of material coordinate s that measures length along the axis of the beam in its reference configuration. In this work, mesh coordinate ξ is also defined and one-to-one map (1) now implies q(ξ, t) = q(s(ξ, t), t).
The mesh and material curvatures in the reference configuration are defined ask 0 = 2L(q 0 )q 0,ξ and k 0 = 2L(q 0 )q 0,s respectively; the corresponding curvatures in the deformed where subscripts 0 and 1 : 3 indicate the scalar and vector parts of dual quaternions, respectively; notation( ·) represent the skew-symmetric matrix associated with cross-product of a vector of dimension three.
The chain rule impliesk 0 = j k 0 andk = j k. The mesh and material velocities of the cross-section arev = 2L(q)q and v = 2L(q)q respectively. Finally, the mesh and material variation of the motion of cross-section areδu = 2L(q)δq and δu = 2L(q)δq respectively. Using the chain rule, identities v =v − ks and δu =δu − kδs are proved easily and echo relationship (3).
Because motion forms a nonlinear manifold, compatibility equations [29] play an important role in the formulation. Compatibility equations arise from the fact that motion is a continuous function of its spatial and temporal variables, and hence, derivatives with respect to these two sets of variables commute. First, motion is considered as a function of material coordinate s and time, leading to where Eq. (10a) results from the commutativity of the spatial derivative and variation, Eq. (10b) from that of spatial and time derivatives, and Eq. (10c) from that of time derivative and variation. Second, motion is considered to be a function of mesh coordinate ξ and time, leading to similar relationshipsδk Taking a mesh time derivative of the material curvature leads to where compatibility equation (11b) and identity (7a) were used to obtain the result. Finally, taking a mesh time derivative of the material velocity v =v − ks and introducing Eq.

Equations of motion
The equations of motion of the problem will be derived from Hamilton's variational principle. Section 4.1 focuses on the variation of the action integral that is required to develop the desired statement of Hamilton's variational principle. The strong and weak forms of the equations of motion, and natural boundary conditions appear in Sect. 4.2, 4.5, and 4.3, respectively. The relationship between the energy stored in the beam and the work done by the externally applied loads is investigated in Sect. 4.4.

Hamilton's variational principle
The beam's Lagrangian is defined as and D ∈ R 6×6 are the symmetric sectional mass and stiffness matrices, respectively. Clearly, the Lagrangian depends on s and t implicitly through the sectional velocities v(s, t) or curvatures k(s, t). The Lagrangian depends on material coordinate s explicitly if the beam is heterogeneous along its reference axis, i.e., k 0 (s), M(s), and D(s) are explicit functions of the material coordinate. The action integral for a beam segment is defined as where t 0 and t f are the initial and final times, and s 0 and s f the root and tip material coordinates of the segment. In view of Reynolds' transport theorem (9), a mesh variation of action integral (14) becomesδ where p = Mv is the sectional momentum vector and f ela = D(k − k 0 ) the sectional force vector. For clarity of the exposition, the terms stemming from the strain and kinetic energies in Eq. (15) will be treated separately.
Introducing compatibility equation (10a), the strain energy component of Eq. (15) becomes where the second line results from integration by parts. Next, attention turns to the kinetic energy component of Eq. (15) and introducing compatibility equation (10c) leads to Integration by parts of the first term on the second line is not straightforward because the order of the temporal and spatial integral cannot be inverted. Consequently, Reynolds' transport theorem (8) is applied to this first term, leading to Combining results (16) and (18) provides the variation of the action as Hamilton's variational principle [29][30][31] states The first term is the variation of the action integral given by Eq. (19). The second term represents the virtual work done by the distributed forces applied along the beam, denoted l ext mat , which are applied at material points of the beam. The third term takes into account the concentrated forces applied at the beam's endpoints; a distinction must be made between end forces applied to material and mesh coordinates, denoted f ext mat and f ext ref , respectively. Symbol ∓ indicates that negative and positive signs must be used at s 0 and s f , respectively.
Introducing variation of the action integral (19) into Hamilton's variational principle (20) then yields where symbol ± indicates that positive and negative signs must be used at s 0 and s f , respectively.

The equations of motion in strong form
Introducing identity δu =δu − kδs into weak form (21), collecting the terms multiplied byδu andδs, and invoking the arbitrariness of these variations yields the strong form of governing equations in domain (s 0 , s f ) aṡ respectively. In the original work of Eshelby [17], the configurational forces were indeed force vectors. In the present development, the term "configurational forces" is used, although  respectively. In summary, the equations of motions for sliding beams in the strong form consist of (1) mechanical momentum equations (22a), (2) mechanical boundary conditions (23a), and (3) configurational boundary conditions (25b). Because configurational momentum equations (22b) are redundant, they should not be added to the set of strong equations of motion. If the axial motion of the beam is prescribed a priori, configurational boundary conditions (25b) are not required and should not be included in the set of strong equations.

The configurational boundary conditions
If the axial motion of the beam is not prescribed a priori, configurational boundary conditions (25b) provide two additional scalar equations that can be solved to determine coordinates s 0 and s f of the material particles at the endpoints of the beam. These boundary conditions are unique to sliding beam problems and must be investigated carefully. Figure 3 shows a beam extending from material coordinate s 0 to s 2 while exhibiting a weak discontinuity at interface s 1 (s 0 ≤ s 1 ≤ s 2 ), i.e., dual quaternions q(s) are continuous at s 1 but curvature k(s) presents a discontinuity at the same point. This situation occurs, for example, if beam segment in [s 0 , s 1 ] is sliding inside a sleeve and that in [s 1 , s 2 ] is unconstrained.
Assume no external load is applied at material point s 1 , f ext mat (s 1 ) = 0. Configurational boundary condition (25b) now reduces to bs − c s 1 = 0, where notation • = •| s + 1 − •| s − 1 denotes the jump of a quantity at s 1 . This condition now expands to In practical applications, the mesh velocity of interface,s 1 , is far smaller than that of elastic waves traveling along the beam. Hence, the inertial terms can be ignored, leading to k T f ela − (k − k 0 ) T f ela /2 s 1 ≈ 0 and a simple algebraic manipulation then yields For an initially straight beam, the curvature vector is k T 0 = {1, 0, 0, 0, 0, 0} and the first term becomes k T 0 f ela = N , where N is the axial force in the beam; the second term is the strain energy density, E strain = (k − k 0 ) T f ela /2. With this definition, the jump condition at s 1 becomes Clearly, the curvature discontinuity generates a discontinuity of the axial force in the beam; an asymptotic analysis of the local equilibrium equation in a small region near the interface [24,32] leads to the same conclusion. The present analysis shows that the curvature discontinuity also generates a discontinuity of the axial force. The additional term E strain s 1 was ignored by previous researchers [10,12,16] although it can be significant, see examples below.

The work and energy principle
The jump conditions generated by a discontinuous curvature vector were investigated in the previous section and have been shown to generate a jump in the axial force, see Eq. (26). This section shows that the well-known work and energy principle is still valid for sliding beam problems, where the system comprises a time varying set of material particles.
where Reynolds' transport theorem (8) is used to yield the first equality, mechanical momentum equation (22a) leads to the second equality, and compatibility equation (10b) yields the last equality. The boundary terms can be recast as where identity E = v T p − L yields the second equality and the configurational momentum b and force c defined by Eqs. (24a), (24b) are introduced to yield the last equality.
Introducing identity (28) where identity v =v −sk is introduced to yield the second equality. As expected, identity (29) indicates that the rate of change of the total mechanical energy stored in the beam segment equals the power of externally applied loads.

The equations of motion in weak form
Introducing identity δu =δu − kδs into Hamilton's variational principle (21) yields where the configurational momentum and force defined by Eqs. (24a) and (24b), respectively, were introduced. In the present derivation, this weak form of the equations of motion is obtained directly from Reynold's transport theorem, Eqs. (8) and (9).
Because the finite element implementation is formulated in the mesh domain, it is preferable to express the material time derivative of the momentum vector in terms of its mesh time derivative counterpart, using identityṗ =p − p ,ss . This transformation gives rise to a spatial derivative of the momentum vector, p ,s = [M(v − ks)] ,s , which involves the spatial derivative of the curvature vector, k ,s . The presence of such term requires the use of an interpolation scheme presenting C 1 − continuity for the dual quaternions to satisfy the compatibility requirements of the finite element method [33]. Because such interpolation is complex, it is preferable to perform a spatial integration by parts of this inertial term, leading to where the natural boundary conditions are now enforced in a weak sense. Clearly, weak form (31a) is equivalent to mechanical momentum equation (22a) and mechanical boundary condition (23a), while weak form (31b) is equivalent to configurational momentum equation (25a) and configurational boundary condition (25b). While the strong and weak form of the governing equations are equivalent, they do not yield identical results within the framework of the finite element method that introduces approximations. The mechanical and configurational momentum equilibrium equations are linearly dependent in their strong form, see Eqs. (22a), (22b). On the other hand, the same equations become independent in their weak form (31a), (31b), which will be shown to provide a sound basis for finite element approximations.
A number of researchers [10][11][12]15] did not perform the spatial integration by parts leading to weak form (31a), (31b). Consequently, their formulations are not compatible and convergence of the finite element approximation cannot be guaranteed [33].

Finite element implementation
The following notation is defined for the finite element implementation of the proposed approach: (·) e I indicates quantities at the I th node in the eth element, array l(ξ ) = {l 1 (ξ ), l 2 (ξ ), . . . , l N (ξ )} stores Lagrange's polynomials of degree N − 1, and notation (·) T = {(·) T 1 , (·) T 2 , . . . , (·) T N } stores the nodal variables. The dual quaternion and mesh velocity fields are interpolated using the schemes developed by Bauchau et al. [34][35][36] where I k indicates the identity matrix of size k and symbol ⊗ denotes the Kronecker matrix product. The material coordinate is interpolated as where where w =vs ,s −v ,ss −kvs and ξ e 1 and ξ e N are the starting and ending mesh coordinates of the eth element. Only the first-order spatial derivative of the shape functions appears in the discretized equations because integration by parts of the inertial forces was performed, as discussed above. This is expected for second-order partial differential equations. In contrast, the formulations presented in the literature by numerous researchers [10][11][12]15] involve second-order spatial derivatives of the shape functions. These formulations are needlessly complicated, reduce the accuracy of the predictions, and are incompatible, as discussed above.
The finite element implementation presented thus far is incomplete. Indeed, system (35) provides 6n node equations for the 6 displacement and rotation components at the n node nodes of the beam. This system, however, also involves two additional unknowns, s 0 and s f , the material coordinates of the endpoints of the beam. Clearly, two additional equations are required to complete the formulation.
The weak and strong forms of the governing equations presented in Sects. 4.5 and 4.2, respectively, can both be used to complete the implementation and lead to two alternative formulations, called the domain-based and boundary-based formulations, respectively. While both formulations use the weak form of the mechanical momentum equations, leading to discrete equations (35), they differ in their treatment of the configurational momentum equations: the domain-based and boundary-based formulations stem from the weak form of configurational momentum equation and boundary conditions (31b) and from the strong form of natural boundary conditions for the configurational forces (25b), respectively.

The domain-based formulation
Introducing the interpolated quantities into the weak form of configurational momentum equation and boundary conditions (31b) leads to  Systems (35) and (36) now provide the 6n node + 2 equations required for the solution of sliding beam problems.

The boundary-based formulation
The strong form of the natural boundary conditions (25b) for the configurational momentum equations can be stated as where subscripts 0 and f indicate quantities evaluated at s 0 and s f , respectively. Equations (37a), (37b) require material curvatures k 0 and k at the endpoints of the beam. An accurate estimate of these curvatures can be obtained by extrapolation of the corresponding quantities at the Gauss points of the element. Figure 4 shows the nodes and Gauss points of the element located at the endpoint of the beam; a four-node, three-Gauss-point element is shown in the figure. Let k n * 0j and k n * j , j = 1, 2, . . . , N − 1, denote the Gauss points material curvatures in the reference and deformed configurations, respectively. Extrapolation based on Gauss point material curvatures yields the corresponding quantities at the endpoint node as whereľ j are the Lagrange polynomials with abscissae located at the Gauss points. Systems (35) and (37a), (37b) now provide the 6n node + 2 equations required for the solution of sliding beam problems.

Numerical examples
Two sliding beam dynamic problems are investigated using the proposed domain-and boundary-based formulations. A stiff beam problem is considered in the first example, while a very flexible cable is investigated in the second example. The generalized-α scheme developed by Arnold and Brüls [37,38] is used to integrate the equations of motion. Clearly, if the beam slides completely in or out of the sleeve, the governing equations of the problem need to be switched from those obtained from the proposed formulation to those of the traditional beam formulation. For clarity, the numerical examples focus on dynamic behavior before the switch happens. To validate the proposed approach, an ABAQUS model of this problem was developed featuring surface-to-surface contact pairs to simulate the interaction between the beam and sleeve. In the ABAQUS model, the two rigid surfaces of the sleeve are placed at a horizontal distance 2d = 0.02 mm, their tips present a 90 degree arc with radii of curvature r = 0.01 mm, and the beam is meshed with 400 quadratic elements. In the present analysis, the beam segments inside and outside the sleeve are both meshed to four four-node beam elements. The beam segment inside the sleeve is constrained to translate along the direction of i 2 only. For both approaches, the simulation was run for 0.5 s using a constant time step of 0.5 ms. Figure 6a shows the components of displacement for the material particle of the beam at coordinate s 2 , denoted u 1 and u 2 along unit vectors i 1 and i 2 , respectively. Figure 6b shows the corresponding velocity components, denoted v 1 and v 2 . The predictions of proposed domain-and boundary-based formulations are presented together with those of the ABAQUS model; the three predictions are found to be in close agreement.
At the beginning of the simulation, the beam moves downward under the effect of gravity and bend to the right due to the applied transverse force. The resulting curvature field generates configurational forces at point T, which in turn, cause an upwards motion of the beam during time interval t ∈ [0.24, 0.37] s. These configurational forces are associated with the jump of strain energy density at point T shown in Fig. 7. Equation (26) shows that the jumps in strain energy density and axial force must be of equal magnitudes at point T. Finally, the material coordinate of the beam located at point T, s 1 , is shown in Fig. 8.
Although the domain-and boundary-based formulations give nearly identical results, as shown in the figures above, their accuracy was assessed by carrying out a convergence study. The tip displacement and velocity vectors were computed for various mesh sizes and relative , and where the reference solutions, u ref tip and v ref tip were obtained using a mesh of 64, three-node beam elements. Figure 9 shows the displacement and velocity error measures as the number of elements increases from 2 to 32. As expected, the domain-based formulation is third-order accurate, while the boundary-based formulation is second-order accurate only.
The inertial, gravity, and external forces applied to the upper portion of the beam are equilibrated by reaction forces at the tip of the sleeve. In the ABAQUS contact model, these reaction components arise from the integrated effect of the contact pressures acting between the beam and the two side walls. These small zones are located near the tip of the sleeve. Consequently, the ABAQUS model requires 400 elements to resolve accurately the location of these contact zones and the associated pressure distributions. In the present approach, the same reaction forces arise from the configurational boundary conditions discussed in Sect. 4.3. These conditions involve discontinuities in the axial force and strain energy density; accurate predictions are obtained using only four four-node beam elements to model each of the two portions of the beam. Clearly, the proposed approach achieves great computational efficiency at the expense of increased complexity of the formulation.  The cable is of length 2L + πR, where L = 1.0 m. At the initial time, the material coordinates at the end points of the cable are denoted s 0 = 0 and s 3 = 2L + πR whereas those located where the cable enters and exits the tube are denoted s 1 = L and s 2 = L + πR, respectively, as shown in Fig. 10. At the initial time, the cable is in static equilibrium under the gravity forces, g = 9.8 m s −2 . At time t = 0, external forces F 1 = 0.1 sin(3πt) N and F 2 = 0.04 sin(6πt) N, acting along unit vectors i 1 and i 2 , respectively, are applied at material coordinate s 3 , causing the cable to sway and move through the semicircular tube.

Cable sliding in semicircular tube
The cable, modeled as a beam, consists of three segments, one in contact with the tube and the other two outside the tube; each segment is meshed with eight four-node beam elements. The portion of the cable located inside the tube is constrained to translate along the direction tangent of the semicircle only. The simulation was run for 0.23 s using a constant time step of 0.2 ms. Figure 11a shows the components of displacement for the material particle of the cable at coordinate s 3 , denoted u 1 , u 2 , and u 3 along unit vectors i 1 , i 2 , and i 3 , respectively. Figure 11b shows the corresponding velocity components, denoted v 1 , v 2 , and v 3 . The predictions of domain-and boundary-based formulations are in excellent agreement.
The time history of the material coordinate of the beam's particle located at point A, denoted s 1 (t), is shown in Fig. 12. Due to the transverse forces applied at the end of the cable, it is pulled out of the circular tube at point B: at the end of the simulation, s 1 (0.23 s) ≈ 0, i.e., the portion of the cable initially hanging out of the tube has entered it through point A and is now inside the circular tube. Jumps in strain energy density are expected to arise at points A and B, and are shown in Fig. 13. Note that the jumps of strain energy density and the resulting configurational forces are of the same order of magnitude as gravity forces; clearly, they are not negligible.
Here again, the relative accuracy of the domain-and boundary-based formulations will be assessed by means of a convergence study. Figure 14 shows the errors in tip displacement and velocity vectors, as defined by Eq. (39), for meshes of increasing size. The number of four-node beam elements in each segment was increased from 8 to 64 and the 128-element mesh was selected as the reference solution. The domain-based formulation is fourth-order accurate, while the boundary-based formulation is second-order accurate only.

Conclusions
This paper addressed the problem of sliding beams with special emphasis on situations where the axial motion of the beam is not prescribed a priori. Hamilton's variational principle was used to derive the weak and strong forms of governing equations for sliding beams. Because the action integral is defined on a non-material domain, Reynolds' transport theorem was applied to transform the variation of the action integral into the integral of the variation of the integrand, and boundary flux terms then appear.
Consequently, the variation of the integrand yields the strong form of the mechanical and configurational momentum equations and the boundary flux terms yield the natural boundary conditions for the mechanical and configurational forces. The configurational momentum equations were shown to be a linear combination of their mechanical counterparts and are therefore redundant.
A weak form of the same equations was also developed by performing spatial integration by parts. In this weak form, the configurational and mechanical momentum equations become independent because they combine in an integral sense the strong mechanical and configurational momentum equations with their respective natural boundary conditions.
The configurational boundary conditions were scrutinized closely and shown to involve discontinuities in the axial force and strain energy density. While the first discontinuity was identified by numerous researchers, the second appears to have been ignored by most, although both were shown to be of the same order of magnitude.
In view of these observations, two formulations were proposed: (1) the domain-based formulation that combines the weak forms of mechanical and configurational momentum equations, and (2) the boundary-based formulation that combines the weak form of mechanical momentum equations with the natural boundary conditions for the configurational forces.
Numerical examples have been presented to validate and compare the two proposed formulations. The predictions of both formulations are found to be in good agreement with those obtained from an ABAQUS model using contact pairs. The domain-based formulation has a convergence rate of three or four for three-or four-node beam elements, respectively, while that of boundary-based formulation is one order lower. Clearly, the proper treatment of the configurational forces at the boundaries impacts the accuracy of the model significantly.

Author contributions All authors wrote the main manuscript.
Funding Supported by the National Natural Science Foundation of China (Grant No. 12172042).

Availability of data and materials
Requests for data and materials should be made to the corresponding authors.

Declarations
Ethical approval Not applicable.

Competing interests
The authors declare no competing interests.