Journal of the Mechanics and Physics of Solids

We develop a comprehensive, geometrically-exact theory for an end-loaded heavy rod constrained to deform on a cylindrical surface. The cylinder can have arbitrary orientation relative to the direction of gravity. By viewing the rod-cylinder system as a special case of an elastic braid, we are able to obtain all forces and moments imparted by the deforming rod to the cylinder as well as all contact reactions. This framework allows for the monitoring of stresses to ascertain whether the cylinder, along with its end supports, is able to sustain the rod deformations. As an application of the theory we study buckling of the constrained rod under compressive and torsional loads, as well as the tendency of the rod to lift off the cylinder under further loading. The cases of a horizontal and vertical cylinder, with gravity having only a lateral or axial component, are amenable to exact analysis, while numerical results map out the transition in buckling mechanism between the two extremes. Weight has a stabilising effect for near-horizontal cylinders, while for near-vertical cylinders it introduces the possibility of buckling purely due to self-weight. Our results are relevant for many engineering and medical applications in which a slender structure is inserted into a cylindrical cavity.


Introduction
The deformation of slender rod-like structures constrained to a rigid cylindrical surface is encountered widely in science, engineering and medicine.Applications include the buckling of drill strings within boreholes used in offshore oil-drilling operations (Hajianmaleki and Daily, 2014), the insertion of stents within arteries and veins in endoscopic surgery (Schneider, 2003) and the soft robotic inspection of pipes in gas, oil or water supply systems (Rashid et al., 2020).
Buckling of such constrained rods inside boreholes has been studied extensively in the drilling industry.For convenient access to oil reserves, deviated and even horizontal wellbores are frequently used.The production of shale gas usually involves a combination of horizontal drilling and hydraulic fracturing of the shale stratum (Bahrami et al., 2013).Compressive as well as torsional buckling have been studied under these inclined conditions, with weight of the rod, generally having both axial and lateral components, taken into account (Paslay and Bogy, 1964).Weight is generally considered to be important in the buckling process.It is common practice that drill strings are compressed only in the lower sections and buckle under their own weight.Bucking is also a known problem, and weight a design concern, in flexible wormlike climbing robots (Wang and Yamamoto, 2017).
Critical buckling loads can be obtained from straightforward linear analysis.Nevertheless, there is a bewildering amount of approximate analyses in the drilling literature; see Cunha (2004) for a review of the sometimes conflicting results on critical loads.The study of post-buckling behaviour requires nonlinear analysis.Whenever nonlinearity is considered in the literature it usually arises in an ad hoc way by combining linear beam-column equations with the nonlinear cylindrical constraint (Mitchell, 1988(Mitchell, , 2002)), rather than performing a systematic nonlinear analysis.
Here we develop a general geometrically exact (hence nonlinear) theory for a rod constrained to deform inside a cylinder of arbitrary orientation relative to the direction of gravity.Our only assumption is that the rod is at all times in contact with the cylinder.Buckling loads depend on boundary conditions.We formulate boundary conditions for the axial and torsional loading of a clamped-clamped rod and find exact critical loads (up to negligible numerical error).Other boundary conditions could easily be applied instead.
Our approach is to view the rod-cylinder system as a braid of two strands winding around each other at constant distance.The exact equilibrium equations for such a 2-braid have recently been derived (Starostin and van der Heijden, 2014).We just need to specialise the general theory to the case where one of the strands is rigid, which we do by imposing kinematic constraints.The advantage of this approach is that we have access to all forces and moments in both the rod and the cylinder, including the contact reactions between the two.Exact expressions for all these individual strand forces and moments can be obtained explicitly in terms of the overall braid forces and moments, which have to be computed numerically by solving the nonlinear (quasi-rod) overall braid equilibrium equations.
The rod-cylinder contact problem is found to be 3-fold statically indeterminate, which means that for a full determination of all strand forces and moments three constitutive conditions need to be specified that characterise the precise nature of the rod-cylinder contact.In previous work (Shah and van der Heijden, 2023) we have shown how static friction can be rigorously incorporated into the theory, but here we make the common assumption of hard frictionless contact, in which contact between rod and cylinder is maintained solely by a normal contact pressure.Friction does not normally play an important role in drilling, especially in the case of vertical wellbores, and is ignored in most cited works.
In this paper we focus on the effect of gravity.As may be expected, for cylinders sufficiently close to horizontal, gravity, like friction (Shah and van der Heijden, 2023), has a stabilising effect in the sense that it delays compressive buckling.Once buckling has occurred, gravity has a further stabilising effect in that it helps maintain the cylindrical constraint.However, eventually, under sufficiently high loads, the rod will lift off the cylinder, an effect not usually considered in the literature.By monitoring the normal contact pressure we are able to accurately determine the lift-off load (similar to the determination of the critical load for a heavy rod to lift off a horizontal plane in van der Heijden et al. (1999)).Other forces and moments can similarly be monitored to ascertain whether the cylinder is able to sustain the required reactions under the applied loads.For instance, the cylinder (and its end supports) has to resist the twisting moment induced by the rod deforming on it.Being able to compute such forces and moments will likely be more important for free-standing tube or pipe systems than for holes bored in a wider solid.
We find that for near-horizontal cylinders gravity introduces mode-switching at critical weights, with heavier rods gradually buckling into more oscillatory modes.Smaller inclinations from the vertical lead to smaller critical compressive loads, as expected, and a gradual transition to tension-dominated buckling in which buckling under self-weight is prevented by applying a tension rather than induced by applying a compression.In both limiting cases of a horizontal and a vertical cylinder, in which gravity only has a lateral or axial component, the critical buckling conditions can be obtained analytically.Intermediate inclinations require numerical solution.
The paper is organised as follows.Section 2 sets up the braid theory and lists the equilibrium equations.In Section 3 these are specialised to the case of a rod on a cylinder and then slightly extended with the effect of gravity, resulting in a 15-dimensional system of ordinary differential equations (ODEs).Individual strand equations are introduced and explicit expressions for all forces and moments, including the normal contact pressure, are obtained.In the process all Lagrange multipliers introduced to impose the constraints are given a natural physical interpretation.In Section 4 the boundary-value problem for compressive buckling is formulated, with boundary conditions carefully designed to load the rod and not the cylinder.Critical buckling conditions are obtained and parameter continuation is used to compute post-buckling solutions and bifurcation diagrams, as well as to map out the transition from horizontal to vertical buckling.By monitoring the normal contact pressure the tendency of the rod, if it was merely lying in one-sided contact, to lift off the surface is investigated.Torsional buckling is briefly discussed as well.Finally, Section 5 ends the paper with a discussion.
Along each curve we define two moving orthonormal frames as follows (see Fig. 1).We first introduce the unit chord vector   () = (1∕)[  () −   (𝑠)].Along the first strand the so-called braid frame is then defined by {  ,   ,   }, where   ∶=   ×   , while the material frame is defined by the vectors {  ,   ,   }, where   and   lie in the plane orthogonal to   , with   ∶=   ×   .These two frames are linked by the first strand's twist angle  1 that is measured from   to   and describes a rotation about   .Similarly, along the second strand, its braid frame is defined by {  ,   ,   }, with   ∶=   ×   , and its material frame is defined by {  ,   ,   }, where   and   lie in the plane orthogonal to   , with   ∶=   ×   .These two frames differ by a rotation about   , through the The consecutive rotation sequences about   ,   and   through the angles  1 ,  and  2 respectively, relating all four orthonormal frames to each other, are as follows: where for  = 1, 2 are orthogonal rotation matrices parametrised by the angles  1 ,  2 and  respectively.
In order to represent the braid and material curvature components for each strand, four skew-symmetric 3 × 3 matrices ω, ω, Ω and Ω are defined for the four reference frames through (2) Here, represent the skew-symmetric matrix ŵ in so(3) and axial rotation vector  in R 3 respectively.A prime denotes differentiation with respect to the first-strand arclength .Thus (2) defines the braid and material curvature components for each strand in the form of four axial vectors  = ( 1 ,  2 ,  3 ) T ,  = ( 1 ,  2 ,  3 ) T , ω = ( ω1 , ω2 , ω3 ) T and Ω = ( Ω1 , Ω2 , Ω3 ) T , which respectively denote the rotation vectors of the braid frames {  ,   ,   }, {  ,   ,   } and material frames {  ,   ,   }, {  ,   ,   }.The braid and material curvature components may then be inter-related with each other by combining (2) with (1) to obtain the relationships (5) The position vector of the second strand is given by   =   +    .Differentiating this with respect to  and using the inextensibility conditions for each strand, given by   =  ′  and  ′   =  ′  , allows the following relationships to be obtained: The second of these reveals that the arclength parametrisation of the braid is regular provided  ∈ (−∕2, ∕2) and  3 < 1∕, conditions that will be satisfied by all solutions we consider (see Starostin and van der Heijden (2014) for a discussion of these conditions).

Overall braid equilibrium equations
The energy functional  for the two-strand braid may be constructed as The first term represents the total elastic strain energy density   for the braid expressed purely in terms of the first-strand braid curvatures  = ( 1 ,  2 ,  3 ) T and the strand twist angles  1 and  2 .The second term incorporates the work done on the braid by the applied end load .The final term arises due to the inextensibility constraint for the second strand, which is enforced with the help of the constant Lagrange multiplier ℎ.
In order to eliminate  1 in favour of the braid angle  = ( 1 ,  3 ), which has a clearer physical meaning, the Lagrangian function  (,  ′ ,  1 ,  ′ 1 ,  2 ,  ′ 2 , ) may be transformed into a function  ( 2 ,  3 , ,  ′ ,  1 ,  ′ 1 ,  2 ,  ′ 2 ,  1 ), thereby modifying the integrand of (7) into by using (6) to eliminate  1 and  ′ and replace the derivatives  ′ 1 and  ′ 3 by  ′ .This in turn transforms the constitutive relations (9) into where (6) 1 has again been used to obtain explicit expressions for ∕ 1 in (11) and ∕ 3 in (13).(c) Phase equations for the strand twist angles  1 and  2 , with respect to the transformed Lagrangian function  : which may be rewritten as a set of first-order equations where the new variables  1 and  2 represent the internal strand torques within the first and second strands respectively.It is worth noting from (15) that  1 and  2 become first integrals of the system in cases where the Lagrangian function  does not explicitly depend on the angles  1 and  2 .
If the braid is subject to further constraints and  depends on corresponding Lagrange multiplier variables (such as ,  and  2 below), then those variables are also subject to standard Euler-Lagrange equations analogous to (14).

Kinematic constraints
Consider a transversely isotropic, inextensible, unshearable, intrinsically straight and linearly elastic rod of length  that is constrained to lie on a rigid cylinder of radius ; see Fig. 2. Rod and cylinder are viewed as a two-strand braid with the rod as the first strand and the cylinder as the second.The overall braid is subjected to an axially applied end force  (taken positive for tension and negative for compression) and a twisting moment , as shown in Fig. 2. The centreline of the rod has constant distance  from the axis of the cylinder.This situation is therefore described by the braid kinematics of Section 2.2.
Since the second strand's material frame is constant along the strand (cylinder), we can take it as our fixed Cartesian reference frame, so we set {  ,   ,   } = { ̂, ̂, ĵ}.Using this, along with (2) 4 and (3)-( 5), the braid and material curvatures for each strand can then be written as from which the following relationships can be deduced: Eq. ( 19) 2 then enables the first-and second-strand braid curvatures to be expressed purely in terms of the braid angle  as

Braid equilibrium equations
To incorporate the effect of gravity into the theory we let the cylinder have an arbitrary angle  with the direction of gravity (the vertical).Since the cylinder is assumed rigid, it is not allowed to sag; its weight is balanced by equal and opposite reactions set up in the cylinder (and its supports).The rod is taken to have mass per unit length .The acceleration due to gravity is denoted by .The rod's weight has components in both the axial ( ̂) and lateral ( ̂) directions to the cylinder (see Fig. 2).The axial component is treated as a force external to the braid and therefore added to the overall braid force balance equation ( 8) 1 .The lateral component is treated as a force internal to the braid and is accounted for by means of a potential energy term in  .The moment balance equation ( 8) 2 is unaltered by gravity.An independent self-contained variational derivation of the equations that follow is given in Appendix.
The first of our complete set of governing equations for the heavy rod on a cylinder are therefore the overall braid force and moment balance equations, which read in componential form For easy reference, a governing differential equation for a specific variable is here labelled by this variable within square brackets [ ].
The remaining set of equilibrium equations are dependent on the precise form of the transformed Lagrangian function  from (10), in which the linearly elastic strain energy   for the overall braid can be expressed as where the general relationship ω2 2 + ω2 3 =  2 2 +  2 3 obtained by combining (4) 1 and ( 5) 1 is used, along with (3) 1 and ( 6) 1 , to re-express the energy density function in terms of the first-strand braid curvature components  2 ,  3 and braid angle , after elimination of  1 . and  are the bending and torsional stiffnesses of the rod respectively.
We need to add to  the gravitational potential energy density   of the distributed weight  of the rod in the ̂ direction: where the position vector of the rod's centreline on the cylinder is corresponding to the bottom of the cylinder (see Fig. 2).Note that cos  =   ⋅   =  ′  ⋅ ̂ =  ′ .Insertion of ( 28) into (10), along with the inclusion of ( 29), then gives the final Lagrangian  ( 2 ,  3 , ,  ′ ,  ′ 1 ,  2 ,  ′ 2 , , ,  2 ,  1 ) for the heavy rod on a cylinder: Here the kinematic constraints (20) 2 and (20) 3 for the variables  2 and  3 are appended with the help of the Lagrange multipliers  and .Furthermore, due to the Lagrangian's dependence on the second strand's twist angle  2 as a result of gravity, the final term containing the Lagrange multiplier  2 needs to be introduced to enforce the kinematic constraint (19) 2 for  2 .The multiplier is called  2 in anticipation of its role as a twisting moment associated with the angle  2 , by analogy with its definition (15) in Section 2.
The first set of remaining equations follows from the phase-angle equations ( 15) with respect to the variable  1 , which yield i.e., the first strand's internal torque  1 is constant, and we deduce, with the help of (6) 1 , the following equation for  1 : Due to the rigidity of the second strand, Eq. ( 15) does not apply for  2 .Instead,  2 is a Lagrange multiplier whose equation is obtained from the Euler-Lagrange equation for  2 : The next set of equations arises from the application of the standard Euler-Lagrange equations to (30) for the three Lagrange multipliers ,  and  2 , which lead to reproducing the constraint equations (20) 2 , (20) 3 and ( 19) 2 for  2 ,  3 and  2 .
The final set of remaining equations consists of the constitutive relations ( 11)-( 13), which after making use of (31) 1 and (20) 3 , provide expressions for the overall braid moment components  1 ,  2 and  3 as

R. Shah and G.H.M. van der Heijden
These are turned into differential equations for the remaining variables  2 ,  and  as follows.First, multiplying (37) by sin  and (39) by cos  and adding them gives an expression for the constant ℎ, ℎ = (  3 −  3 ) cos  − ( 1 +  1 ) sin  +  cos , (40) which in itself represents an additional conserved quantity for the system of equations.Subsequently differentiating (40) with respect to  and using Eqs.( 35), ( 34), ( 25) and ( 27) for the derivatives  ′ 3 ,  ′ ,  ′ 1 and  ′ 3 respectively, along with the constancy of  1 from (31) and the constitutive relation for  2 given by ( 38) results in a differential equation for  given by ) .
[] (41) An analogous governing equation for  can be obtained by substituting the expression for ℎ from (40) directly into (37) and resolving for  ′ , which yields The final equation for  2 is acquired by differentiating (38) and using the known expressions for the derivatives  ′ 2 and  ′ from ( 26) and ( 42) respectively, along with the algebraic expressions (20) 1 and (20) 3 for  1 and  3 , which results in To be able to plot rod configurations in three-dimensional space we append to the system of equations the equation for , This gives a total of 15 ODEs for the variables We observe that in addition to the overall braid equations ( 22)-( 24), gravity enters the system of equations only through (33).The equation for  2 is the equation of a circular shaft with distributed torque  at circumferential angle  2 .We therefore see that the lateral gravity component acts as an internal (to the braid) distributed torque from the rod onto the cylinder.

Intrastrand forces and moments
Once the solution for the overall braid is acquired from its equilibrium equations in Section 3.2, the individual strand equations that follow become particularly relevant for the determination of the contact force and moment reactions experienced within each of the two strands.
Projection along the chord vector direction   results in The first strand's constitutive relations for its moments  (1) 1 ,  (1)  2 and  (1)  3 are the standard ones for a hyperelastic rod and can be obtained from (9) if  is replaced by the strain energy function for the first strand (equal to   as given in (28), in which (6) 1 can be used to reintroduce  1 ) and the resulting first equation is combined with (32): Since the second strand represents a rigid cylinder, free of elastic strains, there are no second-strand constitutive relations for its moment components  (2) 1 ,  (2)  2 and  (2) 3 ; instead, they are reactions maintaining the rigidity constraints. (2)  2 is given by (46) 2 , but  (2)  1 and  (2)  3 are undetermined quantities, in addition to one of the pair  (1) 2 ∕ (2) 2 from (46) 1 .Indeed, given a braid solution (, ), we have nine equations for twelve strand force and moment components.The problem is thus 3-fold statically indeterminate.It is necessary to make three constitutive assumptions characterising the specific form of inter-strand contact in order to obtain a fully closed solution for all braid and strand variables.
Finally, from force and moment balance we also have the following relationships between the contact reactions (deducible from Eqs. ( 45), ( 48), ( 49) and the braid balance equations ( 22)-( 24), which in vectorial form read  ′ =  cos  ̂,  ′ +   ×  = ): () + ( () +    ×  () ) ′ = . (51) These equations therefore serve as an efficient means of directly obtaining the contact loads in the second strand once those in the first strand have been determined.Before specifying the three constitutive contact conditions, we note that independent of the contact model employed we can already conclude from (49) 1 that  (1)  1 = 0, after directly substituting in the constitutive relations for  (1)  1 ,  (1) 2 and  (1)  3 from ( 47).The vanishing of the contact moment  (1)  1 is a consequence of the special form of the strain energy density   in (28).It need not hold for more general   in (10).For instance, it does not hold for transversely anisotropic rods with different bending stiffnesses in two principal directions of the cross-section.For such rods the torque  1 is not constant and (47) 1 would lead to a non-zero term on the left-hand side of the equation for  (1)  1 in (49) 1 .(Note that for consideration of such anisotropic rods within the braid theory they would be required nevertheless to have circular cross-section in order to satisfy the constant-distance assumption of the theory.) We now impose frictionless contact by specifying  (1) 2 = 0,  (1) 3 = 0 and  (1) 3 = 0 as the three contact conditions.From this choice all force and moment components can be obtained.The calculation is lengthy but analogous to that for the weightless case in Shah and van der Heijden (2023) to which we refer for details.For the remaining first-strand forces we find and while  (1) 1 = 0, confirming that this choice of constitutive contact conditions indeed gives the frictionless case, with the normal contact pressure  (1)  2 the only non-zero contact load acting on the rod.The analysis is continued by evaluating the internal and distributed contact loads within the second strand, i.e., the cylinder, for which we find and With these expressions all Lagrange multipliers are given natural physical meanings.The constant ℎ enforcing inextensibility of the second strand (cylinder) becomes tension  (2) 1 in the cylinder; the multiplier  2 for the twist angle  2 , as anticipated, becomes the twisting moment  (2)  1 in the cylinder;  for the  2 constraint becomes the bending moment  (2) 2 ; and finally, by using (1) 2 to write the first term in (59) 3 as    ⋅   , we see that  for the  3 constraint becomes the magnitude of the bending moment about the circumferential axis   of the rod's first braid frame contributing a bending moment about the circumferential axis   of the cylinder.
For the distributed contact loads we find: R. Shah and G.H.M. van der Heijden in particular confirming that the distributed twisting moment  (2)  1 acting on the cylinder corresponds to the twisting moment  2 induced by gravity.
In the weightless limit ( = 0,  2 = 0) all strand force and moment expressions above reduce to those in Shah and van der Heijden (2023) for the case of hard frictionless contact, with  (1)  2 and  (2) 2 in the chord direction the only non-zero contact loads.From this complete instrastrand analysis the advantages offered by the braid theory over traditional single-rod theories become clearly apparent, particularly its ability of enabling the computation of all the reactions imparted to the cylinder by the rod.
For unilateral contact, physically realistic solutions require the normal contact pressure  (1) 2 to have uniform sign.When  (1) 2 < 0 the rod is deemed to be lying on the inside of the cylinder, while  (1)  2 > 0 has it lying on the outside of the cylinder, with the cylinder providing the required pressure, pointing inward to the cylinder in the former case and outward in the latter.We will only deal with the former case in this paper.Where  (1)  2 changes sign, the rod has to be deemed to lift off the cylinder unless prevented by a normal force, for instance adhesion or a contact force in case the rod is placed between two narrowly spaced concentric cylinders providing bilateral support.

Boundary-value problem for a clamped heavy rod in a cylinder
As an application of the equations we formulate and solve the two-point boundary-value problem for a clamped-clamped heavy rod lying on a (inclined) cylinder.We consider both compressive and torsional loading.We use  as the unit of length and ∕ 2 as the unit of force, so that we can set  = 1 (hence  ∈ [0, 1]),  = 1, ∕ = 1∕(1 + ), where  is Poisson's ratio.To get an impression of a typical magnitude of the dimensionless weight parameter  =  3 ∕ we consider realistic physical parameter values for a steel drill string of circular cross-section with radius  (Anon, 2004(Anon, , 2008(Anon, , 2016(Anon, , 2020)): where  is Young's modulus and  is the mass density.With  = , cross-sectional area  =  2 ,  = , second moment of area  =  4 ∕4, we find but note that  is sensitive to  and for long drill strings the value could very easily be much larger.
We take ∕ = 0.01 and  = 0.3 in all our numerical calculations.

Compressive loading of a clamped heavy rod
Here we consider the case of a straight rod clamped at both ends and lying at the bottom of the cylinder while being compressed by a force  < 0 applied axially to the rod.The 15 boundary conditions for this problem are as follows: The kinematical conditions (B2) and (B3) fix the end point (at  = 0) of the rod's centreline at (−, 0, 0) relative to that of the cylinder's at the origin (0, 0, 0).The angles ,  1 ,  2 form a set of Euler angles relating the rods's body frame {  ,   ,   } to the fixed inertial (and cylinder) frame { ̂, ̂, ĵ}.Conditions (B1), (B2) and (B4) enforce complete alignment of these frames at  = 0, while conditions (B7) and (B8) enforce the end tangents to be coaxial while still allowing for a twist ( 1 ) about the end tangent   (1) through angle  1 .
Conditions (B5) and (B6) are necessary to fix the values for the constants of integration from the differential forms of the algebraic relations ( 35) and (38) for  3 and  2 respectively.The subsequent force and moment conditions (B9)-(B14) at  = 1 specify that the second strand, i.e., the cylinder, carries no internal loads at its ends, with all its force and moment components made to vanish by ensuring  () =  () =  from Eqs. ( 56)-( 59).This means that any end loads applied to the overall braid are considered to simply act only on the first strand, i.e., the rod and not the second strand, i.e., the cylinder.Finally, condition (B15) sets the end value for the axial component of the braid force  ⋅   (1) = ( 1 cos  −  3 sin ) (1) =  1 (1) =  to the applied end compressive force  , after making use of (B8).This leaves the end shear force component  ⋅ ĵ (1) =  ⋅   (1) = ( 1 sin  +  3 cos ) (1) =  3 (1) unspecified, free to take on any required value depending on the solution obtained, as usual in clamped loading.

Linearisation
To find the critical buckling loads we solve the linearised boundary-value problem about the trivial straight solution.The trivial untwisted solution with  1 = 0 is given by:  1 =  + ( − 1) cos ,  2 =  3 = 0,  1 =  2 =  3 = 0,  =  2 =  3 = 0,  1 =  2 = 0,  = 0,  = 0 and  = ,  ∈ [0, 1].Using an overbar to indicate the linearised variables, the linearised problem can be reduced to the following fourth-order equation in terms of ȳ: subject to the four boundary conditions which directly follow from (B2), (B4), (B7) and (B8) upon noting that ȳ =  ξ2 and ȳ′ = η.Note that (64) is the equation for a beam-column on a linear foundation, with the effective foundation stiffness provided by the weight of the rod.The equation can be solved analytically in the case of a horizontal cylinder ( = ∕2) so that gravity only has a lateral component (see Section 4.1.2).It can also be solved in the case of a vertical cylinder ( = 0) so that gravity only has an axial component (see Section 4.1.4).For intermediate orientations of the cylinder (64) has to be solved numerically (see Section 4.1.3).
Having found ȳ, and hence ξ2 = ȳ∕, the variables η and ω2 follow directly as Other equations and boundary conditions yield while the equations can be solved by successive quadrature applying the boundary conditions The shear force F3 can then be evaluated as after which the final equation can be solved by quadrature using the last remaining boundary condition This proves that the boundary-value problem with conditions (B1)-(B15) is well-posed.
The mode switching involves the first pair of modes if ( − )∕2 = 1, the second pair if ( − )∕2 = 2, etc. Switching of the first two modes is of particular interest because the first mode (by which we will always mean the mode with lowest critical load   ) is expected to be stable, while higher modes are unstable.Switching of the first mode occurs at critical values  =  2 ( + 2) 2  4 ,  = 1, 2, 3, …, giving  = 9 4  = 8.7668, 64 4  = 62.3418, 225 4  = 219.1705, . . .(see last column in Table 1).Mode-switching causes a pairwise interweaving of buckling curves coming out of the free-column critical loads at  = 0, as can be seen by looking ahead to Fig. 5.
The analytical results are confirmed by the bifurcation diagrams in Fig. 3 in which the end geodesic curvature  2 (1) = − ′ (1) is plotted against the applied end compression − for various values of mg.The critical loads  are seen to correspond to supercritical pitchfork bifurcations from which symmetric pairs of post-buckling solution branches emanate.The diagrams are obtained by numerically solving the full 15D boundary-value problem with the help of the continuation code AUTO (Doedel and Oldeman, 2007) using  as bifurcation parameter.
At critical values   = , listed in Table 1, two pitchforks coincide and modes switch, as illustrated by the insets in Fig. 3 showing the first four modes.For instance, the first and second modes have switched between Fig. 3(b) and (c), which are for mg values straddling the first critical value in Table 1 with mode numbers (, ) = (1, 3).Associated with the codimension-two modeswitching event is a complicated interaction between the bifurcating curves that depends on the precise nonlinearity of the problem.We find that immediately after mode-switching along the first bifurcating branch there is a further, secondary, symmetry-breaking (pitchfork) bifurcation that, upon further increase of mg, quickly moves to higher loads − (indicated by the triangles in Fig. 3(c)).Out of this bifurcation comes a branch of solutions towards smaller mg values (orange in the figure) that intersects the second-mode branch in another secondary bifurcation before eventually connecting with the symmetrically related first-mode branch (second triangle).Along this branch connecting the first-and second-mode branches the solution is non-symmetric and gradually morphs from a first to a second mode (and then to a first mode again).At the intersection the solution has zero end curvature  2 (1), so in our bifurcation diagrams in Fig. 3 such modal intersections occur along the horizontal axis.At the subcritical symmetry-breaking bifurcation along the first-mode branch stability is lost and a dynamic snap will occur (assuming no lift-off has occurred earlier).This bifurcation pattern repeats itself at successive mode-switchings, as can be seen in Fig. 3(d), for  = 30, just after the second switching, with mode numbers (, ) = (1, 5).Here the third-mode branch has a symmetry-breaking bifurcation at  = −2799.71,well outside the plotted range, and the emanating (cyan) branch of non-symmetric solutions re-enters the plot and intersects the fourth-mode branch in a secondary bifurcation on the horizontal axis.
Fig. 4 shows solutions along the symmetry-breaking branch for  = 15 as they evolve from a first to a second mode (both in thick lines).The figure illustrates the wider symmetry properties of (non-symmetry-broken) modes:  is symmetric about the midpoint of the rod (and  2 anti-symmetric) in one of the modes in an interacting pair of modes, and anti-symmetric (and  2 symmetric) in the other mode of the pair.The order of modes of course changes at a mode-switching.Modes with symmetric  2 have zero end twisting moments  2 (0) and  2 (1) set up in the cylinder.These solutions are therefore self-balancing and do not require a torque to be provided by the end supports.Modes with anti-symmetric  2 , by contrast, require a torque to be provided by the end supports.
Fig. 5 shows the stability diagram.The thick solid curve labelled  ,1 indicates where the straight solution goes unstable (buckles) under compression − (assuming that this happens when the first critical load is reached at the given value of mg ).Because of the mode switching this curve consists of alternating pieces of the interweaving critical curves (dotted) coming out of the classical first and second critical loads at  = 0 (marked by triangles).A similar pair of interweaving curves comes out of the third and fourth critical loads at  = 0 and is included in the figure to illustrate the nature of the mode switching.
Also included in Fig. 5 is the lift-off curve labelled  lift where  ∶=  (1) 2 = 0 at some arclength point  along the solution (see  panels below) and lift-off of the rod from the cylinder is initiated.As expected, higher values of mg give larger ranges of post-buckling loads − before lift-off occurs.In the weightless case ( = 0) the rod immediately lifts off upon buckling.The panels below the stability diagram show solutions under increasing compressive loads with solutions at lift-off highlighted and their 3D shapes shown.
The interaction between (pairs of) modal branches in Fig. 5 makes the present case of a beam on an effective linear foundation under fixed-fixed boundary conditions markedly different from the more frequently treated case of such a beam under pinnedpinned boundary conditions (Atanackovic, 1997).In the latter case the critical loads are given by  1 =  and  =  2 1 + ∕ 2 1 , or − =  2  2 + ∕( 2  2 ), i.e., separate straight lines for modes  = 1, 2, 3, ….

Inclined cylinder (0 ≤ 𝛼 < 𝜋∕2)
Fig. 6 shows stability diagrams for inclination angles  = ∕3, ∕6, ∕18, ∕36, ∕90 and 0, corresponding to 30 • , 60 • , 80 • , 85 • , 88 • and 90 • inclinations from the horizontal.Both primary buckling (solid) and subsequent lift-off (dashed) curves are drawn.As expected, the stable region for post-buckling solutions, between these two curves, shrinks as the cylinder approaches the vertical orientation, where lift-off occurs immediately upon buckling.We also observe that under decreasing , the stability curves gradually turn downwards in a clockwise rotation, reflecting the increasing importance of the axial weight component described by the fourth term in (64), often neglected in the literature (e.g., (Mitchell, 2002)).This rotation eventually leads to intersections with the horizontal axis, corresponding to critical weights mg at which the rod buckles purely under self-weight with  = 0.For larger mg the rod has to be held in tension  > 0 to prevent buckling.If the end loading device cannot provide a tension, then a discontinuous jump will occur when mg is gradually increased past the critical values.The critical weights mg are very sensitive to variations in : a change in inclination from 88 • to 90 • leads to a decrease in critical mg from 125.95 to 74.63.
In the vertical limit ( = 0) the stability curves are very nearly straight lines (close inspection reveals a slight curvature) given by 0.5290  −  = (2) 2 (first mode), 0.5143  −  = 8.9868 2 (second mode). ( Note that these limiting critical lines do not depend on  (this is an immediate consequence of (64) when  = 0) and therefore also apply in cases without cylindrical constraint; for instance, when in a slender and flexible climbing robot (with clamped ends) adhesion in the front foot is deactivated and the robot is momentarily a free-standing column (Wang and Yamamoto, 2017).We also note that the inclination angle  behaves as an imperfection under which the interweaving curves at  = ∕2 become veering curves that avoid each other rather than intersect.So for  ≠ ∕2 there are no longer critical mg values with double roots.Such curve veering is a well-known phenomenon in degenerate eigenvalue problems, both in the context of buckling and in small vibrations (Perkins and Mote, 1986;Pierre and Plaut, 1989).
In Fig. 7 solutions at different inclinations are compared.It is seen that inclination breaks symmetry of solutions about their midpoint: they gradually become more sagged towards the  = 0 end.The effect is not strong, however, even for weights as large as mg = 140.Fig. 8 shows 3D shapes at three different inclinations for mg = 140.
Fig. 9 gives the bifurcation diagram for the vertical column at mg = 100 with solution curves presented for the first four modes.The first mode arises in tension, the other modes are all compressive for this value of mg.There is no interaction between modes; there are no secondary bifurcations.

Vertical cylinder (𝛼 = 0)
The case of a vertical cylinder, like that of a horizontal one, is amenable to exact analysis.With  = 0, the linearised equation ( 64) can be integrated once and written as This is the linearised equation for the classical problem of a column under self-weight in addition to a compressive force − (Timoshenko and Gere, 1985).The integration constant  = η′′ (1) in ( 80) is essentially the end shear force  Applying the boundary conditions ( 65) to (82) leads to the following characteristic equation to be satisfied for non-trivial solutions ( 1 ,  2 , c): where .
The integrals of the Airy functions Ai Bi can be evaluated in terms of generalised hypergeometric functions 1  2 ( 1 ;  1 ,  2 ; ⋅), but the Scorer functions Gi and Hi do not seem to have primitives in terms of standard or known special functions.However, Wolfram Mathematica 13.0 (Wolfram Research, Inc., 2023) knows the Scorer functions and a graph of (0, ) is given in Fig. 10.The roots can be calculated numerically and give the following critical buckling weights: = 74.628569, 157.032783, 325.513452, ...

𝑚𝑔
in agreement with the results in Fig. 6(f) and ( 79).Replacing Gi by Hi in (86) simply reflects the graph in Fig. 10 about the horizontal axis and therefore gives the same critical values.
The complication of having to deal with Scorer functions because of the inhomogeneous term in (81) may explain why the case of a column with coaxial clamped ends is not normally considered in textbooks (Timoshenko and Gere, 1985;Atanackovic, 1997), even though this case is one of the standard cases in the buckling of weightless columns.The case is however treated in Engelhardt (1954), where the critical weight is estimated as 74.6285 ± 0.0005 in a careful analysis of series expansions.other direction will destabilise it. 1 > 0 is found to stabilise the solution but only up to a certain amount of twist, beyond which the rod again lifts off.Fig. 11 shows a plot of this twist-induced lift-off for the first mode in the case of a horizontal cylinder.

Rod on top of the cylinder
Our equilibrium equations have a second trivial solution corresponding to the rod lying on rather than in the cylinder.This solution will be unstable but could be stabilised by suitable interaction between rod and cylinder.To study buckling of this solution we only need to replace boundary conditions (B2) and (B7) above by  2 (0) =  =  2 (1).The frame alignment is then changed to t = ̂, d = − ̂, û = − ĵ and we now have η = − ȳ′ .The resulting difference in the linearised equation ( 64) is that the coefficient of the ȳ term changes sign, corresponding to a negative foundation stiffness, which makes sense in this unstable problem as the foundation now 'helps' buckling.We now have two real and two imaginary eigenvalues: 1  2 3 .Numerical results for the case of a horizontal cylinder ( = ∕2) are displayed in Fig. 12.There is no interaction between bifurcating branches in Fig. 12(a) (no mode-switching by interweaving curves).The stability diagram in Fig. 12(b) has the same features as that for the vertical column in Fig. 6(f), with tensile buckling initiated (or prevented) at critical mg values along the horizontal axis, reflecting the role of self-weight in the buckling of a rod lying on top of the cylinder.The critical values occur when  1 =  3 =  1∕4 , where  1 solves cos  cosh  = 1, giving the successive critical values mg = 5.0056, 38.0354, 146.1763, . . .seen in the figure.The curves in Fig. 12(b) are again very close to, but not exactly, straight lines.

Discussion
We have investigated the problem of an end-loaded, heavy, linearly elastic, isotropic rod constrained to deform on a rigid cylindrical surface by employing a recently developed geometrically-exact theory of elastic two-stranded braids (Starostin and van der Heijden, 2014), slightly extended to allow for an external gravity effect.The cylinder is here one of the strands of the braid, constrained to be rigid, while the rod is the other strand.The (two) Lagrange multipliers ( and ) introduced to impose the rigidity constraints on the second strand are shown to have natural physical interpretations as reactive bending moments in the cylinder.
Since the cylinder is assumed to be rigid, gravity acting on it is balanced by equal and opposite reactions set up in the cylinder.Axial gravity of the rod is included as an external effect in our braid formulation and added to the force balance equation for the overall braid, while lateral gravity enters as an internal force to the braid, pulling down the rod but not the cylinder, thereby acting as a distributed torque on the cylinder creating a twisting moment in the cylinder, the third reactive moment,  2 .
A consequence of the rigidity of the second strand of the braid is that the rod-cylinder contact problem is 3-fold statically indeterminate (the elastic braid is 1-fold statically indeterminate (Starostin and van der Heijden, 2014)).This means that for a complete analysis of the problem three constitutive contact conditions have to be specified that characterise the nature of the rod-cylinder interaction.As shown previously (Shah and van der Heijden, 2023), this framework allows for a completely general and rigorous treatment of static friction, but here we make the assumption of hard frictionless contact in which the normal contact pressure  ∶=  (1)  2 is the only non-zero contact load.After specification of the contact conditions all strand forces and moments, including the six contact reactions, can be obtained explicitly in terms of the overall braid variables.Our formulation thus gives access to all reactions in the supporting cylinder, which can for instance be compared against critical strengths to ascertain whether the cylinder, along with its end supports, is able to sustain the required reactions under compression of the rod into the cylinder (in Fig. 4(b), for instance, we have plotted  2 ).As an application of our braid model we have studied buckling of a clamped straight rod under compressive and torsional loading.Torsional buckling is suppressed on a cylinder, while in the absence of self-weight the compressive critical loads are identical to those of an unconstrained column.Weight has a stabilising effect for near-horizontal cylinders, but gives rise to buckling purely caused by self-weight for near-vertical cylinders.By monitoring the contact pressure  between rod and cylinder we have obtained critical lift-off loads as a function of weight as well as inclination.
For the two extreme cases of a horizontal and a vertical cylinder, buckling conditions can be derived analytically.The horizontal case defines a problem of a beam on a linear foundation whose stiffness is given by the distributed weight of the rod.Compressive buckling features mode-switching, known from beam-on-foundation problems governed by a fourth-order equation, in which heavier rods buckle into increasingly more oscillatory modes.The vertical case, on the other hand, is governed by a second-order equation and has (Sturmian) separate non-interacting modes with stability of heavier rods increasingly dominated by (supercritical) buckling under self-weight.Our numerical results for intermediate inclinations map out this gradual transition from interweaving buckling curves in the horizontal case to separate curves (in fact, nearly straight and parallel lines) in the vertical case through the process of curve-veering.
In our numerical calculations we have to fix two dimensionless parameters for which we choose Poisson's ratio  and the slenderness ratio ∕.Results can therefore be interpreted physically most readily for situations in which  is constant.It is good to realise though that fixing ∕ does not fix the relative size of rod and cylinder.Since  is the distance between the axes of rod and cylinder, our numerical results can be interpreted for a one-parameter family of ratios of rod-to-cylinder radii, ∕, such that  −  =  (or  +  =  when the rod is lying on the outside of the cylinder).In the case of a vertical cylinder, primary buckling results are independent of .Critical mg conditions such as in (87) can therefore also be interpreted in dimensional terms as critical length conditions: for given mass per unit length m, the length of the tallest stable fixed-fixed column, for instance a leg of a light raised platform, is given by the critical value  cr = ( 74.628569 ) 1 3 , a factor of 1.58 taller than the critical height of a column whose upper end is free to move sideways.
Our results could be used to extract design formulae for critical buckling loads.Fig. 13 shows comparisons of our numerical results with the empirical fit of Dellinger et al. (1983) (see also (Tan and Digby, 1993;Hajianmaleki and Daily, 2014)) based on experimental data: − ,1 = 2.93 () 0.479 () 0.521  ( sin The agreement for a horizontal cylinder ( = ∕2) is reasonably good.For  = ∕3 the fit deviates from the exact numerical results for large values of mg.Part of the explanation for this will be the frictional resistance (ignored in our study) encountered in practice as a result of which a larger compressive force is needed for buckling.Fits for smaller angles  quickly deteriorate: (89) does not capture the clockwise rotation of the stability diagrams seen in Fig. 6.The formula is not meant to be valid for  ≈ 0, i.e., near-vertical cylinders, where self-weight may become important.Looking at Fig. 6, for such orientations straight-line fits could be used instead, such as in (79) for  = 0. Our theory can be straightforwardly extended to allow for nonuniform bending and torsional stiffnesses  and , as well as nonuniformly distributed mass .Intrinsic curvature can also readily be included (see Shah and van der Heijden (2023) for an example of a helical rod), while the composite rod-cylinder modelling means that any interaction between rod and cylinder, such as adhesion or electrostatics (Starostin and van der Heijden, 2014), can be incorporated.Other boundary conditions could also be considered (several examples are discussed in Shah and van der Heijden (2023)).Conditions could for instance be formulated to study coiled buckling, by allowing one end of the rod to rotate around the edge of the cylinder.So-called helical buckling is widely studied in the drilling literature (Mitchell, 1988;Tan and Digby, 1993;Huang and Pattillo, 2000;Cunha, 2004), but there is a lack of explanation how these 'helices' (i.e., coiled configurations) are supposed to arise in an initially straight rod under the proposed boundary and loading conditions.
Our theory makes no assumption about how the cylindrical constraint is maintained and therefore applies also in situations where the rod winds on the outside of the cylinder, although in that case we require  (1)  2 > 0 if the constraint is to be maintained by hard rod-cylinder contact only.This situation arises in non-buckling problems in which a rod is wrapped around a cylinder and then pulled, such as in the filament winding process under tensile loads used in the manufacturing of composite materials (Sofi et al., 2018).

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.92) reveal that the   are precisely the force components in the fixed frame:  ′ =  cos  ̂,  =  1 ̂ +  2 ĵ +  3 ̂.
Next, we make the following identifications: Differentiating these six relations and using ( 91), ( 92) and ( 93) then produces the balance Eqs. ( 22)-( 27).It may finally be verified that with the braid moment identification in (95) the Euler-Lagrange equations for  2 ,  3 and  give, after some rearrangement, the constitutive relations (37)-( 39), from which the equilibrium equations for ,  and  2 can be derived by differentiation exactly as in Section 3.2.
Had we chosen to include the lateral weight term in the form  sin  in , instead of expressed in terms of  2 , then the  1 equation would have become  ′ 1 =  sin , which would have given additional terms on the right-hand side of all three Eqs.( 22)-( 24).We would effectively have treated lateral weight as external to the braid as well.Moreover, the equation  ′ 2 = (sin )∕ could then simply have been appended to the system of equations without the need for the multiplier  2 .The twisting moment in the cylinder would have been computed as  (2)  1 .All this would not have changed the physics of the problem: critical loads and post-buckling rod configurations would have been the same.

Fig. 1 .
Fig. 1.A braid consisting of two rods at constant distance and their respective reference frames.Centrelines are drawn as thick blue curves while red curves show how the material twists.From: Starostin and van der Heijden (2014).

Fig. 2 .
Fig. 2. Schematic representation of an initially straight, end-loaded rod lying on an inclined cylinder, as described by the braid frame vectors and rotation angles.(left) Front view.(right) Cross-sectional end view at  = .

Fig. 5 .Fig. 6 .Fig. 7 .Fig. 8 .Fig. 9 .
Fig. 5. Stability diagram for the first buckling mode.The critical buckling curve, labelled  ,1 , is made up of alternating pieces of interweaving curves (light dotted lines) coming out of the first and second critical loads on the vertical axis (marked by triangles) owing to the mode switching at mg = 8.7668, 62.3418 and 219.1705.Also shown by dotted lines are similar interweaving curves coming out of the third and fourth critical loads on the vertical axis.The dashed line labelled  lift represents the lift-off curve, where the normal contact pressure  locally vanishes.Panels below display solutions for a sequence of compressions − , with solutions at lift-off highlighted and their 3D shapes shown, at (from top to bottom) mg = 5, 40, 140 and 280 (indicated by diamonds in the stability diagram).All solutions are taken along the downward branches in the bifurcation diagram.( = 0.01,  = ∕2.)

Fig. 12 .
Fig. 12. Rod on top of the cylinder ( = ∕2,  = 0.01.)(a) Bifurcation diagram for mg = 65 showing curves for the first four modes with mode shapes in insets (all taken along the downward branches, i.e.,  2 (1) < 0).(b) Stability diagram showing critical branches for the first four buckling modes.