Exact Solutions of Nonlinear Equation of Rod Deflections Involving the Lauricella Hypergeometric Functions

Abstract

The stress induced in a loaded beam will not exceed some threshold, but also its maximum deflection, as for all the elastic systems, will be controlled. Nevertheless, the linear beam theory fails to describe the large deflections; highly flexible linear elements, namely, rods, typically found in aerospace or oil applications, may experience large displacements—but small strains, for not leaving the field of linear elasticity—so that geometric nonlinearities become significant. In this article, we provide analytical solutions to large deflections problem of a straight, cantilevered rod under different coplanar loadings. Our researches are led by means of the elliptic integrals, but the main achievement concerns the Lauricella 𝐹_𝐷⁽³⁾ hypergeometric functions use for solving elasticity problems. Each of our analytic solutions has been individually validated by comparison with other tools, so that it can be used in turn as a benchmark, that is, for testing other methods based on the finite elements approximation.

1. Introduction

1.1. A Literature's Overview

As early as 1691, Jakob Bernoulli proposed to find out the deformed centerline (planar “elastica”) of a thin, homogeneous, straight, and flexible rod under a force applied at its end, and ignoring its weight. Studying the geometry of its deflection, he established, Curvatura laminae elasticae (1694), the differential equation of the centerline, say 𝑦=𝑦(𝑥), of the rod ab appenso pondere curvata to be 𝑥d𝑦=²d𝑥√1−𝑥⁴. (1.1) Not succeeding in solving it, he was led to look for its arclength, arriving at the same elliptic integral 𝑠=^𝑥₀d𝜉1−𝜉⁴, (1.2) of lemniscate rectification. Daniel Bernoulli, Johann's son and Jakob's nephew, obtained in 1738 subspecie integralis the potential energy stored in a bent rod; in October 1742 he wrote to Leonard Euler the rod should attain such a shape as to minimize the function of squared curvature. Accordingly, Euler dealt elastica as an isoperimetric problem of calculus of variations, De curvis elasticis (1744), arriving at elastica's differential equation and identifying nine shapes of the curve, all equilibrium states of minimal energy. Anyway, Colin Maclaurin had realized the elastica's connection to elliptic integrals, see A treatise on fluxions (1742), § 927 entitled “The construction of the elastic curve, and of other Figures, by the rectification of the conic sections”. In 1757, Euler authored a paper, Sur la force des colonnes, concerning the buckling of columns again, where the critical load is approached through a simplified expression to the elastica's curvature, coming again over the problem again in 1770 and 1775; for all this historical concern, see [1, 2]. Some analytic solutions to elastic planar curves through elliptic integrals of I and II kind can be read at [3–5]. Almost half a century ago, the analytic treatment of flexible rods inspired a nice book, [6], which collects many problems ordered according to the constraint type, and all solved by the elliptic integrals of I and II kind, while the III kind and Theta functions appear there marginally in the 6th Chapter concerning the three dimensional deformations. In order to fix just right from the start, the true concept of rod will be recalled from [7] as follows:A rod is something of a hybrid, that is, a mathematical curve made up of material points. It has no cross section, yet it has stiffness. It can weigh nothing or it can be heavy. We can twist it, bend it, stretch it and shake it, but we cannot break it, it is completely elastic. Of course, it is a mathematical object, and does not exist in the physical world; yet it finds wide application in structural (and biological) mechanics: columns, struts, cables, thread, and DNA have all been modeled with rod theory.A rod which offers no resistance to bending is called a string: it is completely flexible. The string model finds applications in the study of long slender structural members such as cables, where an external load dominates the mechanics and the rods flexural stiffness is negligible.

The link between elliptic functions and elastica was always understood so close, to induce the elliptic functions to be plotted through suitable curved rubber rods (G. Greenhill: Graphical representation of the elliptic functions by means of a bent elastic beam, 1876), while the paper [8] is founded on the use of the Weierstrass ℘ function. For recent developments, see: [8–15]. Inextensibility and circular shape of undeformed rod are (numerically) analyzed in [16, 17] where the load is a uniform centrally directed force.

Our object is the analytical solutions to large-deflection problems of a planar slender cantilever under coplanar loads. The slender structures are those much longer in one direction than in the other two, and their deflection will be described by elastic rod theory if one is interested in phenomena on length scales much larger than the lateral dimensions of the structure. For Dirichlet problems relevant to large deformations, see [18, 19]. A simplified version of this theory, beam theory, suffices if only small deflections need to be considered; the majority of structural engineering applications are adequately modelled by a beam; bridges and buildings are not designed to suffer large deformations. But in other applications of structural mechanics, deflections may be large, meaning that nonlinear geometric effects must be taken into account; the full (geometrically nonlinear) rod theory is then necessary. Mathematically, a rod is modelled as a curve with effective mechanical properties such as bending and torsional stiffnesses, see [20]. Although such problems belong to statics, they have a strong relation with dynamics; arclength along the rod plays a role similar to time in a dynamical system. In its simplest form, this analogy (Kirchhoff's Dynamic Analogy: Über das Gleichgewicht und die Bewegung eines unendlich dűnnen elastischen Stabes, 1859, see [21]) is seen between the deformations of an elastic strut and the motions of a swinging pendulum. A less known analogy links a twisted rod and a spinning top, and the precession of a top does correspond to a rod helical deformation. Engineering problems where rod theory is applied include the looping of ocean cables [22] and the buckling of oil well drilling strings, which nowadays can be over 5 km long. Furthermore, there is a deep reason for a great interest in applying the elastic rod theory to the supercoiling and packing of DNA molecules and proteins (such as collagen): as a matter of fact, due to its double helical nature, DNA is unusually stiff (in both, bending and torsion) for a polymer and, hence, well described by a rod. Finally, the study of strings and rods is of interest to electrodynamic space tethers¹, generally credited of great potential for future space missions in order to bring down space debris (derelict spacecraft, etc.). Unlike the majority of tethers, which can be treated as cables, Space European Tether (SET, see [23]) is a prototype designed to withstand torsional and bending forces so that it needs to be modelled as a spinning rod (and not a cable), namely, a 3D-elastica². All above provide enough motivations to implement exact solutions within the rod theory, see [10].The exact governing equations for the finite displacement rod theory become highly nonlinear, and, hence, it is very difficult to solve these equations analytically. Thus, for practical purposes, problems of this kind are usually solved by the methods based on the finite elements approximations. However, the closed form solutions for these problems are still important to the practical point that the accuracy of the approximate methods can be precisely evaluated by these solutions, to say nothing of the mathematical importance.

1.2. Aim of the Paper

It is well known why the elastic curve of a rod is essential, in fact not only the stresses, but also its maximum deflection will not exceed a fixed amount. Furthermore, (see Mattiasson [24]).When developing finite element methods, there is also a need for checking the results against an exact solution. In the case of finite element methods for geometrically nonlinear beam, frame and shell structures, elliptic integral solutions of large deflection beam and frame problems offer such exact solution.

Nevertheless, his approach is a bit obsolete nowadays because he computes the elliptic integrals numerically through the arithmetic-geometrical means.

We are going to analyze an elastic, thin, flexible rod, clamped at one end (cantilever) bent under different loadings whose plane coincides with that of the strained centerline. Equating the exact curvature in cartesian coordinates to the bending moment divided into the flexural rigidity, one is led to a second-order nonlinear and nonautonomous differential equation. Boundary conditions are given both on the clamped end (Cauchy problem), so that, for example, if the constraint is perfect, displacement and rotation have to be zero there. We discarded the illusion of a general system of loads and rods and restricted ourselves to only one configuration, charged by a single load and treating each single problem by means of the special functions of mathematical physics. The chosen loads are constant couple, evenly distributed load, shear at the tip, hydrostatic load, and sinusoidal load. All our force loads have to be understood as “dead load”, namely, always equal to itself and then not “follower” as are currently named the loads keeping a fixed inclination to the rod before and after the strain.

For one of the above loads, based on the rod's unextensibility, we computed in closed form the Ω free tip coordinates, parametrized to the load itself. We use everywhere nondimensional quantities; the reader will see how advantageous is such a choice in producing “universal” relationships which (in elliptical and/or hypergeometric cloth) keep their formal validity in several situations.

1.3. Special Functions Tools

We assume that the reader knows the elliptic functions (see [25] for details), but let us recall hereinafter the less known Lauricella hypergeometric ones.

The hypergeometric series first appeared in the Wallis's Arithmetica infinitorum (1656) ₂𝐹₁(𝑎,𝑏;𝑐;𝑥)=1+𝑎⋅𝑏()()1⋅𝑐𝑥+𝑎⋅𝑎+1⋅𝑏⋅𝑏+1()𝑥1⋅2⋅𝑐⋅𝑐+1²+⋯, (1.3) for |𝑥|<1 and real parameters 𝑎,𝑏,𝑐. The product of 𝑛 factors (𝜆)_𝑛=𝜆(𝜆+1)⋯(𝜆+𝑛−1), (1.4) called Pochhammer symbol (or truncated factorial) allows to write₂𝐹₁ as ₂𝐹₁(𝑎,𝑏;𝑐;𝑥)=^∞_𝑛=0(𝑎)_𝑛(𝑏)_𝑛(𝑐)_𝑛𝑥^𝑛𝑛!. (1.5) The first meaningful contributions about various ₂𝐹₁ representations are due to Euler³; nevertheless, he does not seem to have known the integral representation ₂𝐹₁()=Γ(𝑐)𝑎,𝑏;𝑐;𝑥Γ(𝑎)Γ(𝑐−𝑎)¹₀𝑢^𝑎−1(1−𝑢)^{𝑐−𝑎−1}(1−𝑥𝑢)^𝑏d𝑢, (1.6) traditionally ascribed to him, but really due to A. M. Legendre⁴. For all this and for the Stirling contributions, see [26]. The above integral relationship is true if 𝑐>𝑎>0 and for |𝑥|<1, even if this bound can be discarded due to analytic continuation.

Several functions have been introduced in the 19th century for generalizing the hypergeometric functions to multiple variables, but we will mention only those introduced and investigated by Lauricella (1893) and Saran (1954). Among them our interest is focused on the hypergeometric function 𝐹_𝐷^(𝑛) of 𝑛∈ℕ⁺ variables (and 𝑛+2 parameters), see [27, 28], defined as 𝐹_𝐷^(𝑛)𝑎,𝑏₁,…,𝑏_𝑛;𝑐;𝑥₁,…,𝑥_𝑛∶=_{𝑚1,…,𝑚𝑛∈ℕ}(𝑎)_{𝑚1+⋯+𝑚𝑛}𝑏_1𝑚1⋯𝑏_𝑛𝑚𝑛(𝑐)_{𝑚1+⋯+𝑚𝑛}𝑚₁!⋯𝑚_𝑛!𝑥^𝑚1₁⋯𝑥^𝑚𝑚_𝑛, (1.7) with the hypergeometric series usual convergence requirements |𝑥₁|<1,…,|𝑥_𝑛|<1. If Re𝑐>𝑅𝑒𝑎>0, the relevant integral representation theorem (IRT) provides 𝐹_𝐷^(𝑛)𝑎,𝑏₁,…,𝑏_𝑛;𝑐;𝑥₁,…,𝑥_𝑛=Γ(𝑐)Γ(𝑎)Γ(𝑐−𝑎)¹₀𝑢^𝑎−1(1−𝑢)^{𝑐−𝑎−1}1−𝑥₁𝑢^𝑏1⋯1−𝑥_𝑛𝑢^𝑏𝑛d𝑢, (1.8) allowing the analytic continuation to ℂ^𝑛 deprived of the cartesian 𝑛-dimensional product of the interval ]1,∞[ with itself.

2. Bent Slender Rod: A Large Deflections 2D-Model

Let us tackle a 𝐿-long, thin, and weightless rod whose end A is clamped and B free. We put at B the origin O of a (𝑥,𝑦) cartesian reference frame with 𝑥 along the undeformed rod, from B to A; and 𝑦 downwards, normal at B to 𝑥. Our basic assumptions are(A₁)the rod is thin, initially straight, homogeneous, with uniform cross-section and uniform flexural stiffness 𝐸𝐽, where 𝐸 is the Young modulus, and 𝐽 the cross-section moment of inertia about a “neutral” axis normal to the plane of bending and passing through the central line;(A₂)the slender rod is charged by coplanar dead loads;(A₃)linear constitutive load: the induced curvature is proportional via 1/(𝐸𝐽) to the bending moment. The rod undergoes large deflections but small strains;(A₄)the shear transverse deformation is ignored, and the rod is assumed to keep constant its previous length (inextensibility);(A₅)a stationary strain field, by isostatic equilibrium of active loads and reactive forces, takes place, and, due to rest of static equilibrium, no rod element undergoes acceleration.

The key Bernoulli-Euler linear constitutive law connects flexural stiffness 𝐸𝐽 and bending moment 𝑀 to the consequent curvature 𝜘 of the deflected centerline ±𝐸𝐽𝜘=𝑀(𝑥). With our (traditional) reference frame (𝑥,𝑦), the curvature has always opposite sign to the bending moment 𝑦𝐸𝐽^𝑥𝑥1+𝑦²_𝑥^3/2=−𝑀(𝑥), (2.1) where 𝑦=𝑦(𝑥) is the unknown elastica's cartesian equation and 𝑦_𝑥=d𝑦d𝑥,𝑦_𝑥𝑥=d²𝑦d𝑥². (2.2) To (2.1) a useful dress can be given d𝑦d𝑥1+^𝑥𝑦²_𝑥1=−𝑀(𝑥)𝐸𝐽, (2.3) to be solved together to the perfect clamping boundary conditions𝑦(𝐿)=0,𝑦^′(𝐿)=0. (2.4) Throughout all the paper, it will be minded that the Cauchy problem d𝑦d𝑥1+^𝑥𝑦²_𝑥1=−𝑀(𝑥)𝑦(𝐿)𝑦𝐸𝐽,0≤𝑥≤𝐿,=0,_𝑥(𝐿)=0, (2.5) can be solved in such a way. Putting 𝐻(𝑥)∶=^𝐿_𝑥𝑀(𝜉)d𝜉, (2.6) the solution of (2.5) is𝑦(𝑥)=−^𝐿_𝑥𝐻(𝜉)(𝐸𝐽)²−𝐻²(𝜉)d𝜉. (2.7)

In the sequel, we are going to use dimensionless variables 𝜉=𝑥/𝐿,𝜂=𝑦/𝐿, so that the previous relationship becomes𝜂_𝜉=𝐻(𝐿𝜉)(𝐸𝐽)²−𝐻²(),𝜂(1)𝐿𝜉=0, (2.8) whose solution is𝜂(𝜉)=^𝜉₁𝐻(𝐿𝜁)(𝐸𝐽)²−𝐻²(𝐿𝜁)d𝜁. (2.9)

2.1. Rod Inflected by a Constant Bending Couple at Its Tip

The clamped section is at 𝑥=𝐿, see Figure 1. In such a case, the equation to be integrated is (2.5), with 𝑀(𝑥)=−ℳ₀, and (2.9) provides, in nondimensional form, 𝜂(𝜉)=1𝜒−1𝜒²−(𝜉−1)², (2.10) where 𝜒=ℳ₀𝐿/𝐸𝐽>0. The elastica is more readable reverting at a dimensional form 𝑥²+𝑦²−2𝐿𝑥−2𝜒𝐿𝑦+𝐿²=0, (2.11) half circle of centre 𝐶(𝐿,𝐿𝜒) and radius 𝐿𝜒, see Figure 1, as it is well known.

Figure 1 

Horizontal cantilever rod inflected by a constant bending couple at its tip.

2.2. The Horizontal Heavy Rod

The classic elastica theory started by the “lamina”, namely, a thin rod without weight but capable to withstand the flexure and loaded at its free end only. Chapter 5 of [6] on the specific problem of the horizontal cantilever under uniform load gives (page 173):This problem… has shown considerable resistance to a closed form solution and all investigations dealing with it have employed approximate methods only.

Accordingly, he obtains elastica coordinates as parametric functions of the arclength 𝑠, stopping their Maclaurin expansions to 𝑠¹⁰.

Let us consider at Figure 2 a horizontal, flexible, homogeneous, heavy and cantilever; it will be inflected by the uniform evenly continuous load 𝑞 (weight for unity span length). Its elastica, according to the approximate beam theory, would consist of a fourth-order rational function of 𝑥. Being 𝑀(𝑥)1=−2𝑞𝑥²,𝜈=𝑞𝐿³6𝐸𝐽, (2.12) we obtain, in nondimensional coordinates, using (2.9) 𝜂(𝜉)=𝜈3^𝜉3₀⁻¹𝜎^1/2(1+𝜎)^2/3(1−𝜈𝜎)^1/2(1+𝜈𝜎)d𝜎, (2.13) and finally, passing to the variable 𝜆=𝜎/(𝜉³−1): 𝜂(𝜉)=𝜈𝜉³−1²3¹₀𝜆1−𝜆1−𝜉^32/31−𝜆𝜈𝜉³−11−𝜆𝜈1−𝜉³d𝜆. (2.14)

Figure 2 

Horizontal cantilever rod bent by a uniform load.

By identification, we obtain the solution to our elastica through a Lauricella function of 3 variables 𝑥₁=1−𝜉³,𝑥₂=𝜈(𝜉³−1),𝑥₃=𝜈(1−𝜉³), so that 𝜂(𝜉)=𝜈6𝜉³−1²𝐹_𝐷⁽³⁾22;3,12,12;3;1−𝜉³,𝜈𝜉³−1,𝜈1−𝜉³. (2.15) For being compliant with the IRT requirements, it shall be |𝑥_𝑘|<1, and then 𝜈<1. Otherwise, it will be necessary another method.

We got there three 𝜈-ranges. The first is the very narrow one, where the small deflections assumption is true, first-order theory. By |𝜂^′(𝜉)|<ℎ²,ℎ∈(0,1), we get 𝜈𝜉³−1<ℎ²1−𝜈²𝜉³−1², (2.16) which at 𝜉=0, gives 𝜈<ℎ²/(1+ℎ⁴). For instance, if ℎ=1/10, then(i)if 0<𝜈<0.0099, the first-order approximation can be accepted;(ii)if 0.0099<𝜈<1 by (IRT), the hypergeometric function 𝐹_𝐷⁽³⁾ is capable of describing our solution (2.14);(iii)if 𝜈≥1, the solution (2.14) will be numerically computed.

Formula (2.15), founded on the Lauricella 𝐹_𝐷⁽³⁾ hypergeometric functions, provides a very compact explicit solution to the elastica of a thin clamped heavy rod. We cannot give now the due importance to (2.15) due to the lack of software packages whose quickness and precision of computing are comparable to those available on Mathematica for the Gauss or Appell hypergeometric functions. We did invite Wolfram to consider such a necessity, but no forecasting is possible. As far as we are concerned, the integral (2.14) cannot be led to any other known special functions or to any combination of them.

2.3. A Tip-Sheared Horizontal Rod

We refer to Figure 3, where the force 𝑅 is transversal (in [29] is studied the inclined end load) but not a follower one; we have 𝑀(𝑥)=𝑅𝑥, so that, after defining the nondimensional 𝜇=𝑅𝐿²/2𝐸𝐽, we obtain 𝜂(𝜉)=𝜇¹_𝜉1−𝜁²1−𝜇1−𝜁²1+𝜇1−𝜁²d𝜁, (2.17) being the bending moment always greater than zero, so will be 𝐻(𝑥), and then it will be 𝜂_𝜉>0. For this load too, the hypergeometric treatment is possible. Passing from 𝜁 to 𝜆, with 𝜆=(1−𝜁²)/(1−𝜉²), we obtain𝜂(𝜉)=𝜇21−𝜉²²¹₀𝜆1−𝜆1−𝜉²1−𝜆𝜇1−𝜉²1−𝜆𝜇𝜉²−1d𝜆. (2.18)

Figure 3 

Horizontal cantilever rod bent under tip-concentrated upwards load.

We can see that, recognizing the Lauricella parameters 𝑎=2,𝑏_𝑖=1/2,𝑐=3𝜂(𝜉)=𝜇41−𝜉²²𝐹_𝐷⁽³⁾12;2,12,12;3;1−𝜉²,𝜇𝜉²−1,𝜇1−𝜉², (2.19) compliant to the IRT requirements if 𝜇<1. For this load, we have to the Lauricella's alternative of the elliptic integrals (not applicable to the heavy nor to the submerged rod), so that we will cover with much more details this approach. Investigating the elastica of a thin cantilever, the majority of authors (Landau and Lifchitz [30], Love [21], Mathieu [31], Burgatti [4], Panayotounakos and Theocaris [11, 12, 14], Chucheepsakul and Phungpaigram [9], Basoco [8], Frisch-Fay [6]) use the slope-arclength formulation (𝜗,𝑠) they integrate a first time for getting 𝜗 as a function of 𝑠. But d𝑥=d𝑠cos𝜗(𝑠) and d𝑦=d𝑠sin𝜗(𝑠), so that they get, through a second integration, the elastica's 𝑥 and 𝑦 coordinates parametrized on 𝑠. Only few develop a different (𝑥,𝑦) parametric equations, for example, Saalschütz, in [32]. Anyway a treatment leading to explicit sole cartesian coordinates does not seem to be ever established, and it is exactly what we are doing. An old reference, rather fit to Jakob Bernoulli (1694) and Euler (1744) formulae, is that of Kiepert (1912), who⁵ stops without reducing his elliptic integral to the Legendre standard form. It will be just the elliptic integral (2.17) our startup. First let us change variables in (2.17) passing from 𝜁 to 𝑡=1−𝜁², obtaining 𝜂(𝜉)=12^1−𝜉2₀𝑡1/𝜇−𝑡(1−𝑡)𝑡+1/𝜇d𝑡. (2.20) For the reader convenience, let us fix some notation. Let 𝑎>𝑏>𝑥>𝑦>𝑐. Consider the definite elliptic integral ℐ𝑎,𝑏,𝑐,𝑥,𝑦=^𝑥_𝑦𝑡(𝑎−𝑡)(𝑏−𝑡)(𝑡−𝑐)d𝑡. (2.21) Combining formulae given in [25]: integral 234.16 page 76 and integral 340.01 page 205, introducing 𝑘²=𝑏−𝑐(𝑥)𝑎−𝑐,𝜑=arcsin(𝑎−𝑐)(𝑏−𝑥)(𝑏−𝑐)(𝑎−𝑥), (2.22) we obtain: ℐ=2𝑎,𝑏,𝑐,𝑥,𝑦√𝑎𝐹𝜑𝑦𝜑(𝑥)−(Π𝜑𝑦𝑎−𝑐,𝑘−𝐹,𝑘𝑎−𝑏),𝑘²,𝑘−Π𝜑(𝑥),𝑘²,𝑘, (2.23) where 𝐹(𝜑,𝑘) and Π(𝜑,𝑘²,𝑘) are the I and III kind elliptic integrals. Now assuming that 𝜇<1 we have the identification 1𝑎𝑏𝑥𝑦𝑐↕↕↕↕↕𝜇11−𝜉²0−1𝜇 (2.24) which allows to give (2.17) via (2.20) by means of the integration formula 𝜂𝜉;𝜇=𝜇21𝜇ℱ+1𝜉,𝜇1−𝜇𝒫𝜉,𝜇, (2.25) where: ℱ𝜉,𝜇=𝐹𝜑𝜇,1+𝜇2−𝐹𝜓𝜇,1+𝜇2,𝒫𝜉,𝜇=Π𝜑𝜇,1+𝜇2,1+𝜇2−Π𝜓𝜇,1+𝜇2,1+𝜇2, (2.26) being 𝜑𝜇=arcsin2𝜇𝜇1+𝜇,𝜓=arcsin2𝜇𝜉1+𝜇1−𝜇1²−𝜉². (2.27) Formula (2.25) provides the dimensionless deflection 𝜂 to the position 𝜉 and to the load 𝜇: its 3D-plot, via Mathematica, as shown in Figure 4.

Figure 4 

The sheared cantilever: a 3D dimensionless plot of deflection 𝜂 versus position 𝜉 and load intensity 𝜇.

From it, one could obtain the relevant contour plots cutting by planes of 𝜇= constant. The thing is shown in Figure 5, where the clamped end is on the left. The deflection curves are given for some values of 𝜇∶0.2;0.4;0.6;0.8. Each of them is compared with the dot curve coming from the approximate theory⁶. Mind that 𝜉 is ranged as 𝜉_Ω<𝜉<1, and 𝜉_Ω has to be found from (3.12). Accordingly, each single 𝜇-curve has been plotted between 𝜉=1 and 𝜉_Ω(𝜇). From the physical viewpoint, the lines'⁷ comparison shows at the low loads an agreement which is becoming worst for increasing loads when the approximate theory is progressively failing.

Figure 5 

The sheared cantilever: solution to (2.17) with some 𝜇 values compared to the dotted lines of the linearized theory.

It is worth noting that the power series expansion of (2.25) drives back to the approximate expression for 𝜂(𝜉) used in literature. Of course, we deeply exploited the computer algebra power for this task. First we compute the Maclaurin third-order expansion of (2.25) with respect to 𝜉 obtaining 𝜂(𝜉)=𝜇31−𝜇^23/2𝜉³−𝜇1−𝜇²𝜉+𝜇2𝐴𝜇, (2.28) where (𝜑(𝜇) is the quantity introduced in (2.27)) 𝐴𝜇=1𝜇𝐹𝜑𝜇,𝜇+12+11−𝜇Π𝜇+12;𝜑𝜇,𝜇+12. (2.29) Now we expand the elliptic integrals of I and III kind with respect to 𝜇, with initial point 𝜇₀=0𝐹𝜑𝜇,𝜇+12=√2𝜇^1/2+3160√2𝜇^5/2+𝑂𝜇^7/2,Π𝜇+12;𝜑𝜇,𝜇+12=√2𝜇^1/2+√23𝜇^3/2+𝑂𝜇^5/2, (2.30) arriving at 𝜇2𝐴𝜇=23𝜇+𝑂𝜇². (2.31) Moreover recalling the Maclaurin first-order expansion with respect to 𝜇𝜇31−𝜇^23/2𝜉³−𝜇1−𝜇²𝜉𝜉=𝜇³3−𝜉+𝑂𝜇², (2.32) adding (2.31) and (2.32), we obtain eventually 𝜂(𝜉)=𝜇32−3𝜉+𝜉³+𝑂𝜇², (2.33) that is, the usual approximation.

Going again to formula (2.25), his peculiarities are: it is nondimensional; it is new, in closed form, exact, and explicit. Moreover, for being implemented, no previous transcendental equation has to be solved; it requires only to compute the dimensionless parameter 𝜇=(𝑅𝐿²)/(2𝐸𝐽) holding all the rod physical data.

On the purpose, it will be observed that in [33], an eigenvalue problem is solved through a double boundary condition of no rotation at the ends of the rod. The authors get the coordinates of elastica parametrized to the external load and marked by the mode number. For our same load, they show four elastica each relevant to a given load, but referring to the static mode 𝑛=1. Other curves are provided referring to the first dynamic mode, 𝑛=2. Furthermore, the module 𝑘 of the Jacobi elliptic functions involved there, is on its own depending on the tip load 𝑃 via an equation to be solved previously. Nevertheless the (essential!) relationship 𝑘⇆𝑃 keeps unrevealed, because the authors at page 740 write:The modulus of elliptic functions 𝑘 and the parameter 𝐹₁ play the role of integration constants and they are related to the force 𝑃 and the angle 𝜑₀ by boundary conditions for each case of rod bending.

In fact, they do not provide the 𝑘 computation at all, marking their plot solutions not in terms of 𝑃 (as it should be done), but of a sequence of 𝑘 values 0.5; 0.6; 0.85; 1−10⁻⁵ see Figure  2, page 742 of [33], so that their (𝑥,𝑦) solution seems to be not easily accessible from 𝑃.

2.4. The Rod Inflected by Hydrostatic Pressure

We are looking for the shape of a (𝐿,𝐸𝐽) slender and thin cantilever clamped at the bottom of a pool of calm water. The rod weight has been deliberately ignored, being partially balanced by the Archimedes buoyancy. The rod is at full length submerged, being obvious the adjustments for a rod leaning out the water or whose tip is under the sheet of water. We assume, as in Figure 6, a reference with the 𝑦 axis downwards and 𝑥 over the water surface. The hydrostatic head, being 𝜌 the liquid mass density and 𝑔 the gravity, will have a linearly deep-growing profile (Stevin), so that, integrating twice such a head, one arrives at the bending moment 𝑀(𝑥)=16𝜌𝑔𝐿𝑥³. (2.34)

Figure 6 

Submerged vertical rod under hydrostatic pressure.

By the usual identification, we obtain the solution to our elastica through the Lauricella function of three variables 𝑥₁=𝜎𝜉⁴−1,𝑥₂=𝜎1−𝜉⁴,𝑥₃=1−𝜉⁴, (2.35) where 𝜎=𝜌𝑔𝐿⁵/24𝐸𝐽, so that𝜂(𝜉)=𝜎81−𝜉⁴²𝐹_𝐷⁽³⁾12;2,12,34;3;𝜎𝜉⁴−1,𝜎1−𝜉⁴,1−𝜉⁴. (2.36) For being compliant with the IRT requirements, it will be |𝑥_𝑘|<1, and then 𝜎<1.

2.5. The Rod Loaded by a Sinusoidal Bending Moment

Assumption A3 implies that our systems are ruled by differential equations affected by a geometrical nonlinearity due to the nonlinear link between the curvature 𝜘 and 𝑦_𝑥. But 𝜘, which is the consequence on the rod, is linked to its cause, namely, the bending moment 𝑀, by a Linear Constitutive Law 𝜘=±𝑀/𝐸𝐽, being 𝐸𝐽 the flexural stiffness. Such a linearity implies that, by superposition, whenever some moment distribution is repeated periodically, we can make its Fourier expansion, so that the rod bending moment is broken out into the addition of an infinitely many harmonics sine/cosine. In such a way, one will get the elastica by superposition, of course retaining those harmonics deemed useful or necessary.

What above provides enough motivation for analyzing once for all, the effect of the sinusoidal load, whose a half wave only is taken, omitting all the possible variants or subcases. Seeing Figure 7, let 𝑞(𝑥) be the continuous transverse load (whose highest value is 𝑄) over the unit of an 𝐿-length flexible cantilever 𝑞(𝑥)𝜋=𝑄sin𝐿𝑥. (2.37) Integrating twice such a load, one arrives at the bending moment affecting the 𝑥-section 𝑀(𝑥)𝐿=−²𝑄𝜋²𝜋sin𝐿𝑥. (2.38) In such a way, introducing 𝜀=𝐿³𝑄/(𝜋³𝐸𝐽)>0, with⁸  𝜀≪1, following (2.9), we find 𝜂(𝜉)=−𝜀^𝜉₁1+cos(𝜋𝜁)1−𝜀²1+cos(𝜋𝜉)²d𝜁. (2.39) If 1+cos(𝜋𝜁)=𝑢/𝜀, so that𝜂(𝜉)=1𝜋₀^{𝜀[1+cos(𝜋𝜉)]}𝑢(1−𝑢)(1+𝑢)(2𝜀−𝑢)d𝑢. (2.40) By the integration Tables  [25], entries 254.12 page 113 and 337.01 page 201, we can establish an integration formula, embedding (2.40) as a particular case: let 𝑎>𝑏≥𝑦>𝑐>𝑑 and define𝑘²=(𝑎−𝑑)(𝑏−𝑐)(𝑎−𝑐)(𝑏−𝑑),𝛼²=𝑏−𝑐𝑏−𝑑,𝜑=arcsin(𝑏−𝑑)𝑦−𝑐(𝑏−𝑐)𝑦−𝑑, (2.41) then^𝑦_𝑐𝑢−𝑐()2()𝑎−𝑢)(𝑏−𝑢)(𝑢−𝑑d𝑢=𝑏−𝑐)(𝑐−𝑑(𝑏−𝑑)𝛼²(𝑎−𝑐)(𝑏−𝑑)Π𝜑,𝛼²,𝑘−𝐹𝜑,𝑘, (2.42) where Π(𝜑,𝛼²,𝑘) is again the elliptic integral of third kind and 𝐹(𝜑,𝑘) of first kind. Being 2𝜀<1, in (2.40), we do the identification (𝑎𝑏𝑦𝑐𝑑↕↕↕↕↕12𝜀𝜀1+cos𝜋𝜉)0−1 (2.43) so that our elastica is𝜂()=2𝜉;𝜀𝜋√Π𝜃(),2𝜀+1𝜉,𝜀2𝜀2𝜀+1,2𝜀𝜃()2𝜀+1−𝐹𝜉,𝜀,2𝜀2𝜀+1, (2.44) where 𝜃()√𝜉,𝜀∶=arcsin2𝜀+1𝜀cos(𝜋𝜉)+𝜀+1cos𝜋𝜉2. (2.45) Figure 8, performed through Mathematica, provides the overlapping of the integral (2.39) numerically computed and the plot consequent to our compact smart formula (2.44).

Figure 7 

Elastica of a cantilever rod bent by a sinusoidal continuous load.

Figure 8 

Integral (2.39) numerical, and by formula (2.44) with 𝜀=0.1. The plot shall be conveniently oriented; the rod's clamped end is at the right side, that is, 𝜉=1, and the load is thought to act upwards.

3. How to Compute the Tip Position after the Strain

Considering the small deflections, the deformed centerline's curvature is approximately equated to the deflection's second derivative, assuming the rod deformation field to proceed along purely vertical displacements. If it is so, the free end 𝑂(𝜉=0,𝜂=0) will be mapped to 𝜉=0, with 𝜂(0)≠0, but having the elastica as horizontal projection, its full unstrained length, then the rod will undergo extension.

With large displacements all, this is not true any more, and a further unknown will appear, namely, the abscissa 𝜉_Ω of the new position Ω, where 𝑂 has been mapped to. Rod's behavior under strain can be modeled in several ways. For instance in [9] after integrating the elastica, the authors write five congruence relationships involving geometric quantities before and after the strain and then leave the rod free of extension.

Following a completely different method, in [34] some possible alternatives like rotational load, load always perpendicular to the tangent at the tip and so on, are taken into account.

Anyway, the simplest criterion is the inextensibility assumption A4; the rod keeps its originary length⁹. For the horizontal weightless cantilever inflected by a tip shear (see Figure 3), the elastica 𝜂(𝜉) has been provided by means of the elliptic integrals, see formula (2.25). Wishing to compute the tip deflection 𝑓, the approximate beam theory would answer 𝑓=𝜂(0), but with large deflection, we have first to detect 𝜉_Ω what is done assuming the unextensibility¹_𝜉Ω1+𝜂²_𝜉d𝜉=1. (3.1) We see that 𝜉_Ω is the sole unknown, having been previously found 𝜂_𝜉, so that, after 𝜉_Ω, we will compute 𝑓=𝜂(𝜉_Ω). let us provide the full treatment for such a case. By plugging 𝜂_𝜉=𝜇1−𝜉²1−𝜇1−𝜉²1+𝜇1−𝜉², (3.2) in (3.1), after the change 1−𝜉²=𝑡, we get 2𝜇=^1−𝜉2Ω₀d𝑡1/𝜇−𝑡(1−𝑡)1/𝜇+𝑡. (3.3) Our unknown appears in the upper integration limit; in order to solve (3.3), we need to establish an existence-uniqueness Lemma.

Lemma 3.1. Consider the transcendental 𝑥-equation where 𝑎>𝑏≥𝑥>𝑦≥𝑐^𝑥_𝑦d𝑡(𝑎−𝑡)(𝑏−𝑡)(𝑡−𝑐)=2𝜇. (3.4) Defining 𝑘²=𝑏−𝑐𝑦𝑎−𝑐,𝜑=arcsin(𝑎−𝑐)𝑏−𝑦(𝑏−𝑐)𝑎−𝑦, (3.5) then the equation (3.4) has a unique real root 𝑥 if and only if 0≤2𝜇≤^𝑏_𝑦d𝑡()=2𝑎−𝑡)(𝑏−𝑡)(𝑡−𝑐√𝐹𝜑𝑦𝑎−𝑐,𝑘, (3.6) and such a root is given by 𝑏()()𝑥=𝑐−𝑎+𝑎𝑏−𝑐sn²𝜇√𝑎−𝑐−𝐹arcsin𝜑𝑦,𝑘,𝑘𝑐−𝑎+(𝑏−𝑐)sn²𝜇√𝑎−𝑐−𝐹arcsin𝜑𝑦,𝑘,𝑘, (3.7) where sn(𝑢,𝑘) is the Jacobi sine amplitude of module 𝑘.

Proof. Let 𝑎>𝑏≥𝑥>𝑦≥𝑐. Observe that the function 𝑥⟼^𝑥_𝑦d𝑡(𝑎−𝑡)(𝑏−𝑡)(𝑡−𝑐), (3.8) is strictly increasing, and that (3.6) ensures existence of solution of equation (3.4). In order to prove the Lemma, we will refer to [25], formula 234.00, page 74, to establish that ^𝑥_𝑦d𝑡()=2𝑎−𝑡)(𝑏−𝑡)(𝑡−𝑐√𝐹𝜑𝑦𝜑(𝑥)𝑎−𝑐,𝑘−𝐹,𝑘, (3.9) with the 𝜑 and 𝑘 meanings previously introduced. The very last step follows from the inversion: 𝐹arcsin𝜓,𝑘=𝛼⟺𝜓=sn(𝛼,𝑘), (3.10) which delivers 𝑥 from (3.9), after equating the right hand side to 2𝜇.

We can then use the Lemma to solve (3.3) via the identification 1𝑎𝑏𝑥𝑦𝑐↕↕↕↕↕𝜇11−𝜉²_Ω0−1𝜇, (3.11) where it has been assumed¹⁰  𝜇<1. Minding that 𝑎=1/𝜇 and 𝑐=−1/𝜇, it should be seen that (3.6) provides an upper limit to the dimensionless load 𝜇 solving numerically a transcendental inequality. To such a limit, a physical meaning can be found as the highest allowable load beyond which the linear constitutive law (curvature versus bending couple) is not true any more.

Solving to 𝜉_Ω, we find out that 𝜉²_Ω=1−𝜇²2𝜇sn2𝜇−𝐹arcsin2𝜇/𝜇+1,𝜇+1/2,𝜇+1/2dn2𝜇−𝐹arcsin2𝜇/𝜇+1,𝜇+1/2,𝜇+1/2². (3.12) Therefore, increasing 𝑅 intensities¹¹, and then 𝜇, the successive positions of the free tip, that is, the components of the displacement 𝑂Ω_𝑘=(𝜉_Ω,|𝜂(𝜉_Ω)|) will be described as functions of 𝜇 at Figures 9 and 10.

Figure 9 

Sheared rod: tip's abscissa as a function of the load.

Figure 10 

Sheared rod: dimensionless tip deflection as a function of the load. Dotted line refers to the approximate theory |𝜂(0)|=2𝜇/3.

First, the (3.12) relationship is plotted as shown in Figure 9.

Plugging (3.12) in (2.25), we obtain for the tip 𝜂 as a function of 𝜇 alone. But such a formula is too long, and then we will omit here, providing its relevant plot.

4. Conclusions

Five elasticae have been analytically solved of a thin cantilever rod bent by different coplanar “dead” loadings¹², namely, whose direction remains unchanged during the deformation of the rod. They are concentrated couple, tip shear, evenly continuous load, hydrostatic load, and sinusoidal load. Each of them describes a single physical situation for which we obtain the large deflections in closed form through elliptic and/or hypergeometric functions. Some care has been put (tip sheared rod) in providing analytically, via the isometric assumption, the free end position after the strain, parametrizing its coordinates as a function of the load. Some of our solutions, specially that to the heavy rod by the Lauricella functions, have enriched the literature on the elastica.

Acknowledgments

The authors are indebted to their friend Aldo Scimone who drew some Figures of this paper; they hereby take opportunity for thanking him warmly. They wish also to thank the anonymous referee whose valuable criticism improved their article.

Endnotes

1. A space tether is a long cable used to connect spacecraft to other orbiting bodies such as space stations, boosters, payload, and so forth in order to transfer energy and momentum.
2. The problem of 3D-elastica is very old, but out of our purpose. We will say only that it first appeared in Lagrange's Mécanique Analytique 2nd edition (1811). J. Binet (Mémoire sur l'intégration des équations de la courbe élastique à double courbure, 1844) and M. Wantzel (Note sur l'intégration des équations de la courbe élastique à double courbure, 1844) arrived at the same differential system, vaguely formulated, and not affronted. A modern treatment of a particular case can be read in [35], with its specific bibliography.
3. We quote three works: (a) De progressionibus transcendentibus, Op. omnia, S.1, vol. 28; (b) De curva hypergeometrica Op. omnia, S.1, vol. 16; (c) Institutiones Calculi integralis, 1769, vol. II.
4. A.M. Legendre, Exercices de calcul intégral, II, quatriéme part, sect. 2, Paris 1811.
5. In [36] putting 𝑅=2/𝑎, he gets the elastica 𝑦(𝑥) through the integral 𝑦(𝑥)𝐶=±¹±𝑥²𝑎⁴−𝐶₂±𝑥²d𝑥. (4.1) “Die Kurve welche diese Gleichung entspricht, heißt die elastische Linie weil ein elastischer Stab der an dem Ende befestigt und an dem andern Ende belastet ist, diese Form annimmt.”
6. Under transversal load at the cantilever tip, this problem in strength of materials is solved in cartesian coordinates following the linear approximation; the solution is found (e.g., see [37, page 305]) to be 𝜂(𝜉=𝜇)32−3𝜉+𝜉³, (4.2) in dimensionless form.
7. Notice that the plots are coming after Mathematica's computation of our explicit solution (2.25), but adding the map 𝜉→1−𝜉 in order to provide the geometrical true configuration of the stressed rod.
8. Inserting in 𝜀 some numerical values, for example, 𝐿=5m,𝑞=6N×m⁻¹,𝐸𝐽=100N×m²,𝑔=9,81m×s⁻², one gets: 𝜀≈3×10⁻², so that, unless of assuming unrealistic lengths, 𝜀 is essentially much less than unity.
9. Only in some special cases such a tip path description is quick. For instance for a rod loaded by a fixed couple (Section 2.1), one gets 𝜉_Ω=1−𝜒,𝑓=𝜒.
10. Having for instance 𝐿=1m,𝑅=10N,𝐸𝐽=100N×m², we get 𝜇=𝑅𝐿²/2𝐸𝐽=0.05.
11. Nothing to say, apart from the more strict link imposed by the transcendental inequality (linearity range), (3.12) requires 𝜇≤1.
12. The contrary case is that of the follower load which is being unchanged in its relationship with the rod, for example, keeping its direction perpendicular—or inclined—to it during and after the strain.