## Abstract

In numerous biological, medical and engineering applications, elastic rods are constrained to deform inside or around tube-like surfaces. To solve efficiently this class of problems, the equations governing the deflection of elastic rods are reformulated within the Eulerian framework of this generic tubular constraint defined as a perfectly stiff normal ringed surface. This reformulation hinges on describing the rod-deformed configuration by means of its relative position with respect to a reference curve, defined as the axis or spine curve of the constraint, and on restating the rod local equilibrium in terms of the curvilinear coordinate parametrizing this curve. Associated with a segmentation strategy, which partitions the global problem into a sequence of rod segments either in continuous contact with the constraint or free of contact (except for their extremities), this re-parametrization not only trivializes the detection of new contacts but also transforms these free boundary problems into classic two-points boundary-value problems and suppresses the isoperimetric constraints resulting from the imposition of the rod position at the extremities of each rod segment.

## 1. Introduction

A wide scope of disciplines ranging from medical science to engineering is concerned with the constrained deformation of a slender elastic body. With the nature of the constraint being contingent on the context, the elastic body, which is generally referred to as a *rod*, may be coiled/wound around another slender body or restrained inside a straight or sinuous conduit.

Applications of this general problem include the insertion of miniature instruments into blood vessels for the treatment of vascular and cerebrovascular pathologies (e.g. abdominal and thoracic aortic aneurysms, stroke) as well as the endoscopic examination of internal organs [1,2]. Another application concerns the petroleum industry, which relies on kilometres-long drillstrings to transmit the axial force and torque necessary to drill the rock formations and reach deep hydrocarbon reservoirs [3,4]. At very small lengthscales, the wrapping of DNA around histone cores [5–8] or carbon nanotubes [9] has recently been the subject of intensive research in biophysics. Finally, twining plants, which achieve vertical growth through the friction generated by revolving around poles or other supports [10,11], constitute a further natural manifestation of this relatively broad class of problems.

As discrete and/or continuous contacts develop between the rod and the constraint, modelling of these phenomena usually implies the partitioning of the global problem into a set of boundary-value problems [12,13] specifying both the position and inclination of the rod at its extremities. Each elementary problem corresponds to a segment of rod either in continuous contact with the constraint or free of contact (between its extremities.) The length of the rod spanning each elementary problem and satisfying the associated boundary conditions is, however, *a priori* unknown and, therefore, constitutes an inherent part of the solution. This feature, specific to (one-dimensional) free boundary problems [14], is also encountered in other applications, e.g. the continuous casting process and the extrusion of plastic through a draw plate [15], the laying of subsea cables and pipelines [16,17] or, in slightly different configurations, the elastically deformable arm scale [18] and the self-encapsulation of elastic rods [19]. Resorting to a conventional Lagrangian approach, the axially unconstrained nature of these problems, referred to as *self-feeding* in the following, combined with the boundary conditions specifying the rod location at both extremities lead to the establishment of integral constraints (namely isoperimetric constraints) on the unknown length of the rod.

The rod is further compelled to satisfy a non-penetration requirement ensuring that, depending on the context, the rod remains either inside or outside the constraint. The assessment of this unilateral contact condition requires, in principle, the comparison of two curves parametrized by distinct curvilinear coordinates (i.e. the rod and the constraint axis), an intensive computational task contributing to the numerical burden associated with conventional approaches.

For planar elastica, the above-mentioned difficulties have been circumvented by adopting a Eulerian description of the rod in terms of the curvilinear coordinate associated with the constraint [20,21] and by describing its deformed configuration by means of the signed distance between the rod and the axis of the constraint [22]. This approach has been shown to drastically simplify the resolution of the constrained problem by trivializing the assessment of the unilateral contact condition and limiting the degeneracy of the governing equations with diminishing clearance or bending stiffness. Generalizing these concepts, this paper aims to extend this formulation to three-dimensional configurations. In the following, the special Cosserat rod theory is reformulated within the Eulerian framework of a generic constraint defined as a perfectly stiff normal ringed surface. The axis or spine curve associated therewith, referred to as *reference curve* independently of the constraint nature, is further assumed to be known and regular (i.e. its tangent vector being everywhere well defined.) The proposed Eulerian reformulation hinges on describing the rod-deformed configuration by means of its relative position with respect to the reference curve and expressing the kinematical as well as mechanical quantities pertaining to the rod in terms of the curvilinear coordinate associated with the reference curve, rather than the rod. As the *a priori* unknown domain, *viz.* the rod length, is substituted for the known reference curve, the free boundary problem and the associated isoperimetric constraints are converted into a classical two-point boundary-value problem for which the rod eccentricity with respect to the reference curve is prescribed at both extremities of the domain.

### (a) Nomenclature

Lower case letters (Latin and Greek) are used to denote kinematic quantities pertaining to the rod while upper case is preferred for variables describing the reference curve and the constraint surface. For instance, we denote by *s* the Lagrangian coordinate, identifying a rod cross section along its centreline, and by *S* the Eulerian coordinate, identifying a section along the reference curve. Additionally, vectors of the Euclidean space ** u**, and calligraphic symbols (bold and plain) are used to denote dimensionless quantities, e.g.

*i*and

*j*have ranges over 1,2 and 1,2,3, respectively.

## 2. Lagrangian formulation

### (a) Problem definition

Let us define the right-handed orthonormal basis {*e*_{j}} for the Euclidean space ** R**(

*S*)=

*X*

_{j}

*e*_{j}, the parametrization of the reference curve

*S*, referred to as the

*Eulerian coordinate*, identifies a section along

*S*. Without loss of generality, the parametrization of the reference curve is assumed to be natural or arc-length, so that its tangent is a unit vector. The constraint surface, considered frictionless and undeformable, is the normal ringed surface [23] generated by sweeping a circle of radius

*Q*(

*S*) centred on the reference curve and in the normal plane to

The elementary problem considered in this paper consists of a segment of rod, either in continuous contact with the constraint or free of contact between its extremities, that is forced to go through two fixed points in space while being subjected to both a distributed body force and a body couple. As already noted, the unknown length of the rod is free to adjust to both the external loading and the boundary conditions acting at the rod extremities. Prescribing the rod position and inclination at these points, this self-feeding feature may be obtained by imposing one extremity of the rod to be clamped while considering the other to be moving freely through a frictionless sliding sleeve. The rod, initially straight, is assumed to behave as a one-dimensional elastic body according to the special Cosserat theory of rods. Although readily extendable to encompass effects of shear and extensibility, the present reformulation is, for the sake of brevity, restricted to rods that can undergo large deformations in space by experiencing bending and torsion.

The continuous contact and free of contact configurations essentially differ by the nature of the body force acting on the rod. While solely subjected to its own weight along a contact-free problem, the rod is additionally compelled to lie on the constraint surface along a continuous contact. Assuming frictionless interactions between the rod and the constraint surface, the resulting contact pressure acts as an additional distributed body force oriented normally to the constraint surface, its magnitude being, however, *a priori* unknown.

### (b) Governing equations

In the global {*e*_{j}}-basis, the parametrization of the rod centreline is denoted ** r**(

*s*)=

*x*

_{j}

*e*_{j}, where

*s*is referred to as the

*Lagrangian coordinate*in contrast with the Eulerian coordinate

*S*. The inextensibility assumption implies that the deformation of the rod leaves invariant the distance between points measured along its centreline. The

*material*coordinate

*s*, therefore, identifies a rod cross section independently of the rod configuration. To fully characterize the spatial configuration of the rod, one has to additionally supply this space curve

*d*_{1}(

*s*),

*d*_{2}(

*s*) along two fixed material lines of the cross section, the rod is entirely defined by the two vector-valued functions

*s*

_{a}<

*s*

_{b}corresponding to the extremities

*S*

_{a}<

*S*

_{b}of the canonical problem (figure 1) and such that the rod length is ℓ=

*s*

_{b}−

*s*

_{a}. We will henceforth assume that the parametrizations of both the rod centreline and the reference curve are at least

*C*

^{4}-continuous; localized external loads and discontinuous body forces being, therefore, disregarded.

Setting *d*_{3}(*s*)=*d*_{1}×*d*_{2}, the triplet of *directors* {*d*_{j}(*s*)} constitutes, for each cross section *s*, a material frame in which the deformed state of the elastic rod is naturally described. While unshearbility implies that the normal to the rod cross section is everywhere aligned with the tangent to the rod centreline, *viz.* cross sections remain perpendicular to ** r**(

*s*) of the rod centreline is arc-length. In other words, the unshearable and inextensible assumptions mean that the director

*d*_{3}coincides with the tangent vector to the rod centreline

**(**

*u**s*)=

*u*

_{j}

*d*_{j}relates the rotation of the directors to the strain variables that enter into the constitutive equations:

*u*

_{1}and

*u*

_{2}characterize the local material curvatures and

*u*

_{3}is a kinematic measure of the twist density.

Considering only local interactions between adjacent cross sections of the rod and defining the internal force ** F**(

*s*)=

*F*

_{j}

*d*_{j}and moment

**(**

*M**s*)=

*M*

_{j}

*d*_{j}, the conservation of linear and angular momenta read [24,25, pp. 273–274]

**(**

*f**s*) and

**(**

*m**s*) are the body force and body couple per unit reference length at

*s*, respectively. The potential reaction pressure

**(**

*p**s*) acting along continuous contacts is treated as an equivalent body force.

For rods with a circular cross section made of a linear elastic material and initially straight, the components of the internal moment are related to the strain variables by means of the following constitutive equation [26, §2.4]
*B* and the torsional stiffness *C*. Although the underlying assumption of linear axisymmetric elastic homogeneous constitutive laws recently gained interest in the context of DNA mechanics [27,28], it provides a reasonably good estimate for a wide range of applications without overcomplicating the presentation of the proposed formulation.

### (c) Discussion: boundary conditions and unilateral contact condition

Classic boundary-value problems for {** r**(

*s*),

*d*_{1}(

*s*)} consist of the kinematic equations (2.2)–(2.3), the equilibrium equations (2.4) and (2.5), the constitutive relation (2.6) and a combination of kinematic and mechanical boundary conditions [25, pp. 322–328]. Considering the orthonormality of the directors, the resulting system consists of 15 first-order differential equations which requires a set of 15 boundary conditions. For an exhaustive description of the various families and possible combinations of boundary conditions associated with elastic rods, see [29, pp. 297–300]. The regular boundary-value problems associated with a fixed length ℓ=

*s*

_{b}−

*s*

_{a}have been investigated from various perspectives and closed form solutions have been obtained under particular sets of boundary conditions and loadings [30–32]. These configurations can, however, be regarded as exceptional and one must generally resort to numerical techniques [28,33,34] to solve the general problem.

In the context under consideration, one extremity of the rod is clamped while the other moves freely through a frictionless sliding sleeve under the combined action of known axial force and twisting moment. One may assume, without loss of generality, the extremity *s*_{a} to be clamped such that the triplet
*s*_{b} of the boundary-value problem, the position of the rod and the unit tangent are imposed together with the axial force and twisting moment, i.e.
*s*_{b}−*s*_{a} being unknown, the resulting free boundary problem is well posed. The boundary conditions specifying the locations ** r**(

*s*

_{a}) and

**(**

*r**s*

_{b}) of the rod centreline at both extremities of the domain, however, lead to the establishment of isoperimetric constraints. Indeed, although the rod spatial configuration has been specified through the position vector of the rod centreline and the director

*d*_{1}(

*s*), the parametrization

**(**

*r**s*) of the space curve

*x*

_{j}(

*s*) of the rod centreline in the absolute reference frame {

*e*_{j}}, therefore, require the integration of the kinematic relation (2.2), that is

*s*

_{a},

*s*

_{b}] requires the use of involved numerical techniques, e.g. adapted shooting methods [35] or the tau method [36], which contribute to the numerical burden of this approach.

As noted in Introduction, the assessment of the unilateral contact condition constitutes an additional source of difficulties. To ensure that the rod remains either in continuous contact or free of contact along the whole domain of the problem under consideration, this supplementary constraint necessitates the evaluation of the distance between the rod centreline *s* and *S*, respectively, the evaluation of this distance reveals to be computationally intensive as it requires the identification of the mapping *s*↦*S*(*s*), from Lagrangian to Eulerian coordinates.

Finally, the ill-conditioning of the governing equations for small bending stiffness, the existence of spurious solutions associated with curling of the rod as well as the demand for increasing accuracy of the solutions with decreasing clearance, are sources of various additional complications or other bottlenecks. Consequently, the Lagrangian formulation of the problem is ineffective and laborious to solve in some circumstances [22].

## 3. Eulerian formulation

The proposed reformulation hinges on describing the rod-deformed configuration by means of its relative position about the reference curve and restating the rod local equilibrium in terms of the Eulerian curvilinear coordinate associated with

By analogy with the director basis, one may arbitrarily define a triplet {*D*_{j}(*S*)} constituting a right-handed orthonormal basis for each cross section *S* along the reference curve and such that *D*_{3}=d** R**/d

*S*is the unit vector tangent to

**(**

*U**S*)=

*U*

_{j}

*D*_{j}[37] as

*D*_{j}}-basis such that it constitutes a Bishop frame [38] is advocated as it results in substantially simpler expressions. Also referred to as the parallel transport frame, this adapted frame constitutes an alternative approach to defining a moving frame that is everywhere well defined. Its construction is based on the concept of relatively parallel fields and the observation that, while the tangent vector

*D*_{3}is unique for any given curve, the pair of orthonormal vectors {

*D*_{1},

*D*_{2}} may be chosen arbitrarily provided it remains perpendicular to the tangent vector. In particular, the pair {

*D*_{1},

*D*_{2}} may be defined such that the {

*D*_{j}}-basis smoothly varies along

*U*

_{3}(

*S*)=0, by requiring d

*D*_{1}/d

*S*⋅

*D*_{2}=0 and d

*D*_{2}/d

*S*⋅

*D*_{1}=0.

As presented in figure 1*b*, the position vector for the rod cross section centroid ** r**(

*s*)=

*x*

_{j}

*e*_{j}can, naturally, be expressed as

*eccentricity vector*

**Δ**(

*S*)=Δ

_{1}

*D*_{1}+Δ

_{2}

*D*_{2}is a measure of the rod relative position in the cross section of abscissa

*S*. Besides describing the space curve

**(**

*r**s*) connects the Eulerian and Lagrangian formulations through the mapping

*S*↦

*s*(

*S*).

In the following, derivatives of scalar- and vector-valued functions with respect to the Eulerian coordinate *S* will be denoted by the apposition of a prime while derivatives with respect to the Lagrangian coordinate *s* are explicitly specified.

### (a) Mapping and Jacobian

The reformulation of the local equilibrium (2.4)–(2.5) within the Eulerian framework requires to express the natural derivatives d⋅/d*s* in terms of Eulerian derivatives d⋅/d*S*. The Jacobian of the mapping *S*↦*s*(*S*), from Eulerian to Lagrangian coordinates, is obtained by inserting the decomposition (3.2) of the rod centreline in equation (2.2) and applying the chain rule differentiation
*D*_{3} is a unit vector, the drift between the two curvilinear coordinates is caused by the eccentricity **Δ** between the rod and the reference curve. Opting for either the positive or negative sign in this expression, *viz.* imposing the mapping *s*(*S*) to be injective, prevents the appearance of solutions associated with curling, writhing or other configurations involving self-contact. Considered unlikely to occur in the present context, these solutions are disregarded. Also the positive sign in expression (3.3) is selected.

Similarly, the derivatives of the inverse mapping *s*↦*S*(*s*), from Lagrangian to Eulerian coordinates, can be explicitly written in terms of the Eulerian functions Δ_{i}(*S*) and *U*_{j}(*S*) (with *i*=1,2 and *j*=1,2,3.) Defining the Eulerian function *J*_{k}(*S*) as
*J*_{1}(*S*)=1/*s*′(*S*), the recursive relation
*k*>1.

### (b) Directors and strain variables

Through the introduction of the eccentricity vector, the space curve *D*_{j}}-basis and accounting for the definition (3.1) of the vector ** U**(

*S*), the components of

*d*_{3}in the frame attached to the reference curve are

*D*_{j}}’s be defined such that they constitute a Bishop frame, the former two expressions considerably simplify as

*U*

_{3}(

*S*)=0.

The attitude of the rod cross section being uniquely defined by the director *d*_{3}, the knowledge of the angle between either director *d*_{1} or *d*_{2} and a specified direction is sufficient to fully characterize the rod spatial configuration. Although straightforward, the characterization of this rotation with respect to the normal and binormal to the space curve *viz.* where the rod curvature *κ* vanishes. Alternatively, as depicted in figure 2, the rotation of the rod cross section about the director *d*_{3} can be described with respect to the pair {*k*_{1},*k*_{2}} as
*k*_{1} and *k*_{2} are the images of *D*_{1} and *D*_{2}, respectively, through the rotation mapping *D*_{3} on *d*_{3}. Rodrigues formula [40] yields
** ϑ**=

*D*_{3}×

*d*_{3}, with

*k*_{3}=

*d*_{3}, the triplet {

*k*_{j}} constitutes a right-handed orthonormal basis that coincides with the directors for

*D*_{j}}-basis for

**=**

*ϑ***0**. Entirely defined by the three Eulerian functions

*g*

_{j}(see appendix A), this

*intermediate*basis attached to

*D*_{j}}-basis as well as the rudiments to re-express the strain variables in terms of Eulerian functions. Henceforth, tildes will be used to distinguish components of vectorial quantities in the {

*k*_{j}}-basis from their counterpart in the directors basis, e.g.

*u*

_{j}=

**⋅**

*u*

*d*_{j}and

Analogous to relation (2.3) for the directors, the kinematic of the intermediate {*k*_{j}}-basis along the rod centreline is described by a *fictitious* twist vector ** u** as a result of the relative rotation

*φ*of the directors {

*d*_{1},

*d*_{2}} with respect to the pair {

*k*_{1},

*k*_{2}}. Hence, in view of the kinematic relations (2.3) and (3.13), differentiation of expressions (3.10)–(3.11) leads to

**=**

*u***+**

*w**J*

_{1}

*φ*′

*d*_{3}for the twist vector and

*d*_{3}/d

*s*=

*κ*

**on**

*n*

*k*_{1}and

*k*_{2}, respectively, while

*k*_{1},

*k*_{2}} about

*k*_{3}. Their expressions, in terms of the Eulerian quantities, are obtained by substituting the definition (3.12) for the

*k*_{j}’s in the projection of equation (3.13) on the {

*k*_{j}}-basis

*k*_{j}}-basis may be decomposed into the kinematics of the {

*D*_{j}}-basis, characterized by the Darboux vector

**(**

*U**S*), and its rotation with respect to this basis. The components

**(**

*U**S*) in the {

*k*_{j}}-basis are given in appendix A in terms of the

*g*

_{j}’s and

*U*

_{j}’s.

In conclusion, through substitution of equations (3.10)–(3.11) and (3.14)–(3.16), the constitutive law (2.6) reads
*g*_{j}(*S*) defined in equations (3.7)–(3.9). Provided the body force and body couple per unit reference length may be expressed as functions of the Eulerian curvilinear coordinate, *viz.* ** f**(

*S*)=

**(**

*f**S*(

*s*)) and

**(**

*m**S*)=

**(**

*m**S*(

*s*)), the balance of linear and angular momenta (2.4)–(2.5) may be expressed as

### (c) Governing equations and boundary conditions

Both the Eulerian coordinate *S* and the components of the eccentricity vector are scaled by the known length of the boundary-value problem, *L*=*S*_{b}−*S*_{a}, leading to the introduction of the dimensionless curvilinear coordinate *ξ*=(*S*−*S*_{a})/*L* and eccentricity vector ** δ**(

*ξ*)=

*δ*

_{1}

*D*_{1}+

*δ*

_{2}

*D*_{2}with

*L*, the Jacobians derivatives (3.6) read

*k*≥1. Similarly, the fictitious twist vector and the Darboux vector scale as

*ξ*.

Considering the previous definitions and denoting the characteristic force *F*^{⋆}=*B*/*L*^{2}, the scaling leads to the introduction of the following vector fields
*ν*=*B*/*C*−1 is Poisson’s ratio. Relation (3.20), therefore, becomes
*φ*(*ξ*)=*φ*(*S*(*ξ*)). The projection of the equilibrium equations (3.21)–(3.22) in the {*k*_{j}}-basis then yields
** ω**(

*ξ*) in the intermediate basis read

*d*_{3}in the {

*D*_{j}}-basis given by

The six differential equations (3.32)–(3.37) with expression (3.27) for the Jacobian *g*_{j}’s result in an 11th order system of nonlinear differential equations in the six Eulerian unknowns *φ* and two third-order, (3.35)–(3.36), in the *δ*_{i}’s (with *i*=1,2 and *j*=1,2,3.)

This mixed-order system consequently requires the specification of 11 constants of integration. The Eulerian counterpart to the boundary conditions (2.7) prescribing the clamped extremity of the rod reads
*i*=1,2. At the end *S*_{b} of the boundary-value problem, the rod remains free to move smoothly through a frictionless sliding sleeve, whose position is fixed in space. Both eccentricity and inclination of the rod with respect to the reference curve must, therefore, be imposed. Additionally, according to the constitutive relations (3.31), the known twisting moment acting at this extremity also requires to prescribe the first derivative of the angle *φ*. Hence, the boundary conditions (2.8) becomes
*i*=1,2. As an essential outcome of the proposed reformulation, the isoperimetric constraints on the unknown length of the rod, a source of difficulty in the Lagrangian formulation, disappear and the system of equations (3.32)–(3.37) with the boundary conditions (3.44))–(3.45) constitute a classical boundary-value problem (as opposed to a free boundary problem.) Note that the proposed Eulerian reformulation of the rod-governing equations could be extended to other boundary conditions.

In conclusion, it has been shown that the rod configuration, entirely defined by ** r**(

*s*) and

*d*_{1}(

*s*) in the Lagrangian formulation (2.1), reduces to the knowledge of the three Eulerian functions

*ξ*={0,1} of the canonical problem corresponding to the rod extremities

*s*

_{a}and

*s*

_{b}, respectively. While the components

*δ*

_{1},

*δ*

_{2}of the eccentricity vector describe the rod relative position with respect to the reference curve, the angle

*φ*characterizes the rotation of the rod cross section about the director

*d*_{3}. This representation of elastic rods leads to a reduced coordinate formulation involving a minimal number of degrees of freedom similar to the

*curve-angle*representations proposed in [33,41–43]. In this formulation, the centreline is explicitly represented and the material frame is characterized by a unique angle.

Besides the suppression of the isoperimetric constraints, the proposed Eulerian reformulation of the rod-governing equations drastically simplifies the assessment of the unilateral contact condition, which requires the evaluation of the distance between the two curves

#### (i) Free rod

The *a posteriori* assessment of the unilateral contact condition along contact-free problems is drastically simplified as it reduces to ensuring that a threshold on the magnitude of the eccentricity vector is not violated. This constraint, indeed, merely consists in checking that either ∥** δ**∥<

*ϵ*or ∥

**∥>**

*δ**ϵ*, with the scaled constraint radius

*ϵ*(

*ξ*)=

*Q*(

*S*(

*ξ*))/

*L*, for interior and exterior configurations, respectively.

#### (ii) Continuous contact

Alternatively, although the magnitude of the eccentricity vector is known along continuous contact problems, the magnitude of the reaction pressure ** ρ**(

*ξ*) included in the body force

**(**

*σ**ξ*) and acting normally to the constraint surface is

*a priori*unknown. The nonlinear boundary-value problem (3.32)–(3.45) is, therefore, supplemented by the relation ∥

**∥=**

*δ**ϵ*to close the formulation and ensure the continuous contact along the whole domain. As emphasized in [42,44], one may always recombine equations (3.32)–(3.33) and (3.35)–(3.36) to eliminate the algebraic component,

*viz.*the reaction pressure, from the resulting system of differential-algebraic equations. This approach, illustrated in §4b for the continuous contact with a cylindrical constraint, is equivalent to projecting the governing equations (3.21)–(3.22) along the normal and geodesic directions to the constraint surface, allowing to reduce the system to a set of ordinary differential equations.

## 4. Validation and applications

As suggested in Introduction, the range of disciplines and applications concerned with the constrained deformation of a slender elastic body is potentially very wide. Exploring all the possible features of the proposed Eulerian reformulation exceeds the scope of this paper and, merely, a few of its benefits and limitations are highlighted in the following. Considering a series of elementary problems, the validity of the method is assessed by comparison with known analytical results and numerical solutions obtained using the Lagrangian formulation.

### (a) Planar elastica

Planar deflection can either result from the rod cross section geometry and loading symmetry (these solutions do not need to be stable) or be imposed by external constraints (in which case the transversal force is not necessarily null). The purely bi-dimensional configuration corresponding to the external constraint *d*_{2}(*ξ*)=*e*_{2} is investigated here by assuming *δ*_{2}(*ξ*)=0 and a planar reference curve. The {*D*_{j}}-basis attached to *D*_{2}(*ξ*)=*e*_{2}; hence, the only non-null component of the Darboux vector

Consequently, the components of the director *d*_{3} in the frame attached to the reference curve, cf. equations (3.41)–(3.43), reduce to
*φ*′(*ξ*)=0, the Eulerian formulation proposed by [22] for a constrained elastica is recovered.

### (b) Straight conduit

The helical buckling and post-buckling behaviour of a rod within a cylindrical constraint has been studied extensively [45–50] not only due to the major interest it represents to several areas of science and engineering (e.g. mechanics of climbing in twining plants, drilling industry, fittings in overhead transmission lines or buckling of optical fibres in loose-tube packaging), but also as a result of the position this fundamental problem occupies in advanced continuum mechanics. Here, it is shown that well-known results associated with twisted elastic rods can easily be recovered from the proposed Eulerian formulation of the rod-governing equations.

Along a straight reference curve framed by a constant {*D*_{j}}-basis, the rod assumes a helical deflection with radius *δ*=∥** δ**∥ and pitch angle

*θ*provided the components of the eccentricity vector satisfy

*θ*is also the inclination of the rod on the reference curve as defined in equation (3.12). The Jacobian of the mapping from Eulerian to Lagrangian coordinates reduces then to

**(**

*ρ**ξ*) is equivalent to a body force acting normally to the constraint surface, i.e.

**⋅**

*ρ*

*D*_{3}=0.

In this context, it may be convenient to decompose the shear forces into their normal *δ*_{1} and *δ*_{2} in the equilibrium equations (3.35) and (3.36) leads to
*φ* then read
*ρ*>0 corresponding to an outward-pointing radial reaction pressure. In particular, the free helix is associated with a null reaction pressure corresponding to an axial force

### (c) Helical conduit

The progressive buckling of an initially straight elastic rod constrained inside a helical conduit is investigated in the following. Although axially unconstrained, i.e. free to flow in/out of the conduit, the rod is assumed to be transversally clamped on the conduit axis at both extremities and subjected to an end torque *a*). Assuming the rod to be weightless and maintaining the end torque constant, the magnitude of the axial force is progressively increased until the rod contacts the restraint and a discrete contact subsequently emerges. The constraint considered in the present application is defined by its helical axis, with pitch angle *Θ*=*π*/4 and radius *ϵ*=*Q*/*L* where the length *ϵ*=1/40 and 1/20 are investigated. Anticipating the discrete contact to appear at midspan *ξ*=1/2 for symmetry reasons, the attached {*D*_{j}}-basis is chosen to be the Bishop frame aligned with the Frenet–Serret apparatus at this point.

Although the present application shares multiple features with the problem of an elastica constrained inside a tube that has been investigated in [53–55], it differs by the unprescribed nature of the rod length, i.e. self-feeding. Indeed, while the distance between the extremities of a rod of finite length evolves as a result of the loading, this distance is kept constant in the present analysis and the rod length freely adapts to the loading.

Before analysing the initial contact-free and the subsequent discrete and continuous contact phases of this constrained buckling for

#### (i) Numerical implementation

Various numerical methods have been applied to two-point boundary-value problems [56]. The so-called shooting method, which combines a root-finding algorithm and a high-order differential equation solver, has been frequently applied to the Lagrangian problem associated with rods of known length [5,34,55]. Although this numerical technique has also been successfully extended to the Lagrangian free boundary-value problem associated with a planar elastica and its Eulerian reformulation [22], a collocation method was favoured to treat the mixed-order system of nonlinear ordinary differential equations (3.32)–(3.37) and the boundary conditions (3.44)–(3.45) constituting the boundary-value problem under consideration.

Defining a partition *Π* of the domain [0,1] constituted of *N* subintervals, an approximate solution *B*-splines are chosen as basis functions [58] such that
*k*≥3 is the number of collocation points per subinterval and *n* on the partition *Π*. The unknown *B*-spline coefficients are determined by solving the nonlinear system of 6 *k* *N* equations, obtained by collocating the governing equations (3.32)–(3.37) at Gaussian points, and imposing the 11 boundary conditions. This method has been implemented in Matlab by combining the *Spline Toolbox* [59] to evaluate the *B*-splines and their derivatives with a *line search* (or backtracking) algorithm [60] to solve the nonlinear system of collocation equations.

#### (ii) Free rod

The initial configuration is chosen to be the helical deformation characterized by a uniformly null eccentricity vector; this state is, however, not to be confused with the straight stress-free configuration *ϵ*≥*Σ*. Following a procedure similar to the one presented in §4b but considering solutions that are not homoclinic to the straight configuration, one may solve the governing equations (3.32)–(3.37) with the boundary conditions *ρ*(*ξ*)=0, that is

Maintaining the torque *exact* Eulerian formulation are compared with both the small inclination approximation (cf. appendix B) and the conventional Lagrangian formulation in figure 4 for *ϵ*, is expected to cause subsequent issues related to the non-penetration condition such as the detection of emerging contacts or the distinction between continuous and discrete contacts. This is further illustrated in the following.

It is also seen that the magnitude of the eccentricity vector is maximal at midspan where *δ*_{2}=0 such that the emergence of a contact between the rod and the constraint can readily be verified by ensuring that |*δ*_{1}|<*ϵ* at *ξ*=1/2. The evolution of the eccentricity *δ*_{1} at midspan is pictured in figure 3*b* as a function of the end force *self-feeding* ability characterizing the proposed method is emphasized in figure 3*d* as the scaled rod length ℓ/*L*, which typically corresponds to the domain of the classical Lagrangian formulation, is seen to freely evolve as a result of the loading.

#### (iii) Discrete contact

As the eccentricity ∥** δ**∥ reaches the clearance

*ϵ*, a discrete contact emerges between the rod and the constraint giving rise to a reaction force

*a priori*unknown, the global problem is partitioned into two elementary boundary-value problems, each corresponding to a rod segment free of contact. Although both the position and inclination of the rod are prescribed at the extremities of the global problem, only the magnitude of the eccentricity vector is known at the junction of the two elementary problems. Resorting to a shooting method [56], the missing boundary conditions are substituted with educated guesses and the rod integrity is restored by ensuring the continuity of the internal moment across the contact, i.e.

The magnitude *L* are presented in figure 3*c* and *d*, respectively, as the end force varies. For consistency purposes with the preceding phase of the loading, these quantities are plotted in the scaling associated with the global problem. Points marking the appearance of the discrete contact are identified with circles. As anticipated, the small inclination approximation is seen to slightly overestimate the magnitude of the end force at which the contact emerges and, therefore, underestimate the intensity of the reaction force (for the same end force.) Owing to the symmetry of the rod-deformed configuration, the position of the discrete contact along the global problem is, however, properly identified. Fields obtained with the Eulerian formulation and its small inclination approximation are depicted in figure 5*a* for *ϵ*=1/20 and

Further increasing the magnitude of the end force may, however, lead to a violation of the unilateral contact condition and an additional restriction on the curvatures of both the rod and the constraint surface at the contact point should be considered. This criterion, which also marks the transition form discrete to continuous contact, ensures that the normal curvature of the constraint surface is smaller than the curvature of the osculating circle to the space curve

#### (iv) Continuous contact

Along a continuous contact, the rod curvature *κ* can be decomposed between its normal *κ*_{n} and geodesic *κ*_{g} components according to [39]
** N** is the inward pointing unit normal to the constraint surface. While

**=−**

*N*

*D*_{1}for the inner (

*δ*

_{1}=

*ϵ*) and

**=**

*N*

*D*_{1}for the outer (

*δ*

_{1}=−

*ϵ*) discrete contact, it can be shown that, due to the symmetry of the rod-deformed configuration, the unit normal vector to

**=**

*n*

*D*_{1}at these points irrespectively of the contacting side. Additionally, following Euler’s theorem, the normal curvature reads

*D*_{2}and

*D*_{3}and the angle

*θ*is defined as the inclination of the director

*d*_{3}on the reference curve according to equation (3.12). Therefore, projecting equation (4.17) on the unit normal vector

**=±**

*n***, it can readily be shown that the contacts remain discrete provided the rod curvature satisfies**

*N**κ*/

*κ*

_{n}

**⋅**

*n***≥1. Considering expression (4.18) for the normal curvature**

*N**κ*

_{n}and the above-mentioned definitions of the principal curvatures

*κ*/

*κ*

_{n}

**⋅**

*n***is depicted in figure 3**

*N**e*as the end force varies. Points marking the transition from discrete to continuous contact, i.e. at which the inequalities are strictly satisfied, are identified with squares.

At this stage of the loading, the method used to solve the discrete contact problem ceases to be valid and the global problem should be partitioned into three elementary boundary-value problems. As the extent and geometry of the central problem, corresponding to the rod in continuous contact with the constraint surface, are *a priori* unknown a shooting method similar to the one previously introduced can be used. Likewise, the missing boundary conditions at the connections between elementary problems are determined by ascertaining the integrity of the rod along the global problem, that is by ensuring the continuity of the internal moment

The reaction pressure along the continuous contact is indeed supplemented by two concentrated reaction forces, *κ*/*κ*_{n}** n**⋅

**=1 (figure 3**

*N**e*). The resulting two reaction forces similarly partition the global problem into three elementary boundary-value problems, each corresponding to a rod segment free of contact. The discrepancy between the solutions obtained with the Eulerian formulation and its small inclination approximation seems to only slightly affect the rod overall response. As shown in figure 5

*b*for

*ϵ*=1/20 and

As the magnitude of the axial end force is further increased, the reaction pressure along the central continuous contact becomes locally outward pointing (i.e. negative). This elementary problem must, therefore, be split to allow these parts of the rod to leave the constraint and preserve the validity of the unilateral contact condition. This forthcoming stage of the loading is, however, not analysed here as it is beyond the scope of the paper.

## 5. Conclusion and limitations

The Eulerian formulation of elastic rods proposed in this paper was motivated by the wealth of biological, medical and engineering applications concerned by the constrained deformation of elastic rods and the need to efficiently solve this class of problems. This formulation hinges on (i) describing the rod-deformed configuration

The restatement of the rod local equilibrium in terms of the Eulerian curvilinear coordinate is particularly appropriate to treat free boundary problems associated with rods forced to go through two distinct fixed points in space. This re-parametrization indeed allows to substitute the unknown domain of the problem, namely the rod length, with the length of the reference curve between the extremities of the problem and, therefore, reduce the resolution to that of a classical boundary-value problem. The isoperimetric constraints that would otherwise ensue from a conventional Lagrangian formulation also vanish. Additionally, the *a posteriori* assessment of the unilateral contact condition ensuring that, depending on the context, the rod remains either in continuous contact or free of contact along an elementary problem has been trivialized through the definition of the eccentricity vector. The detection of emerging contacts throughout contact-free problems hence reduces to the comparison of the eccentricity magnitude with the constraint radius. Alternatively, along continuous contacts, the rod is compelled to lie on the constraining surface by imposing the magnitude of this vector.

The proposed formulation reaches its limits with the single-valuedness of the mapping *S*(*s*), from Lagrangian to Eulerian coordinates, and its reciprocal *s*(*S*), from Eulerian to Lagrangian coordinates. Geometrically, these knotty situations emerge either as a section along the reference curve is crossed multiple times or if the eccentricity vectors at distinct Eulerian coordinates intersect. The former situation, which usually arises in conjunction with deformed configurations characterized by curling, writhing or other solutions involving self-contact such as folding, can be identified by the orthogonality of the director *d*_{3} and the tangent vector to the reference curve *D*_{3}, i.e. situations where *g*_{3}(*ξ*) vanishes at least once along the domain. Alternatively, the singular configuration associated with the latter situation occurs when the eccentricity vector aligns with the normal vector *D*_{3}′ to the reference curve and its norm is at least equal to the reference curve curvature *g*_{3}(*ξ*) to remain strictly positive throughout the domain. Therefore, situations in which the scalar product *d*_{3}⋅*D*_{3} changes sign cannot be captured by the proposed method.

Although the Eulerian formulation proposed in this paper is restricted to quasi-static solutions of the governing equations, it can readily be extended to incorporate rod dynamics by redefining the Jacobians and the eccentricity vector as time-dependent. Similarly, the assumption of frictionless interactions between the rod and the constraint surface can be relaxed. However, it will require to deal with evolutive problems as solutions become then history dependent. Finally, the definition of the constraint surface may be extended to include constraints with non-circular cross sections, e.g. by describing its local geometry in the {*D*_{1},*D*_{2}}-basis, and account for its deformability, which is essential to biomedical applications.

## Authors' contributions

The authors have contributed equally to this work.

## Competing interests

We have no competing interests.

## Funding

A.H. was supported by the National Fund for Scientific Research of Belgium.

## Acknowledgements

We thank the anonymous referees for their thorough reading of our manuscript and their salient comments and suggestions.

## Appendix A. Intermediate basis

The {*k*_{j}}-basis is the image of the {*D*_{j}}-basis through the rotation mapping *D*_{3} on *d*_{3}. The axis describing the rotation, therefore, reads ** ϑ**(

*ξ*)=

*g*

_{1}

*D*_{2}−

*g*

_{2}

*D*_{1}such that, upon application of Rodrigues formula [40], the

*k*_{j}’s are given by

*g*

_{j}’s are the components of the unit tangent

*d*_{3}in the basis attached to the reference curve, see (3.7)–(3.9). The components of the Darboux vector

**(**

*U**S*) in this basis read

*U*

_{j}(

*S*)=

**⋅**

*U*

*D*_{j}.

## Appendix B. Small inclination approximation

The developments presented in §3 lead to the rigorous reformulation of the special Cosserat rod theory within the Eulerian framework associated with the reference curve. Although this reformulation hinges on the description of the rod-deformed configuration by means of its relative deflection about ** δ**(

*ξ*). Along interior problems, the rod eccentricity is, however, expected to be at most equal to the constraint radius, while for many typical exterior problems, the rod remains within close range of the reference curve such that the norm ∥

**∥ be of the order of the constraint radius. The following rescaling of the eccentricity vector is, therefore, naturally proposed**

*δ**ε*a representative value of

*ϵ*(

*ξ*) along the elementary problem under consideration such that

*ε*→0. For small values of ∥

**∥, i.e. for**

*δ**ε*≪1, the space curve

For eccentricities of order *ε* such as previously defined, the inclination of the rod on the reference curve is expected to be comparably small, that is *ε*→0, and, consequently, the {*k*_{j}}-basis to be expressed as a perturbation of the {*D*_{j}}-basis attached to the reference curve. Expanding the rotation vector ** ϑ**(

*ξ*)=

*D*_{3}×

*k*_{3}in powers of

*ε*and inserting the resulting expression in the definition (3.12) of the intermediate basis leads to

*ε*→0, where the first corrections read

**(**

*ϑ**ξ*) being, by definition, orthogonal to

*D*_{3}(

*ξ*), it may be expressed as

*ε*→0, where the corrections read

*ϑ*

_{1,k}(

*ξ*) and

*ϑ*

_{2,k}(

*ξ*) functions of the

Comparison of expression (B 2) with the definition of the intermediate basis in terms of the *g*_{j}(*ξ*)’s (see equations (A 1)–(A 3)) allows to identify the perturbation expansions for the components of the director *d*_{3} in terms of the *D*_{j}}-basis, i.e. Bishop frames, the only non-zero components of the Darboux vector *D*_{1},*D*_{2}}, and the normal and binormal to *ε*→0.

Considering only the leading order terms, the Jacobian of the mapping from stretched to Eulerian coordinates reduces to
*ε*→0, which emphasizes that the drift between these two curvilinear coordinates becomes negligible either for reasonably small curvature of the reference curve or when the eccentricity vector ** δ**(

*ξ*) is orthogonal to the osculating plane to

*contribute*to its curvature. In many applications, the reference curve is only slightly tortuous and its radius of curvature is large compared with the length

*L*of the problem. Hence, assuming the product

*ε*→0, the third terms in equations (B6) and (B7) become negligible and the Jacobian reduces to

*ε*→0, leading to a substantial simplification of the governing system of equations (3.32)–(3.37)

*ε*→0 and in the absence of body couple.

- Received August 6, 2015.
- Accepted May 31, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.