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

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 F 3 D 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.


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 y y x , of the rod ab appenso pondere curvata to be dy

International Journal of Mathematics and Mathematical Sciences
Not succeeding in solving it, he was led to look for its arclength, arriving at the same elliptic integral 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   International Journal of Mathematics and Mathematical Sciences   3 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 1 , 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 3Delastica 2 . 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.

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.

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 for |x| < 1 and real parameters a, b, c. The product of n factors λ n λ λ 1 · · · λ n − 1 , 1.4 called Pochhammer symbol or truncated factorial allows to write 2 F 1 as The first meaningful contributions about various 2 F 1 representations are due to Euler 3 ; nevertheless, he does not seem to have known the integral representation traditionally ascribed to him, but really due to A. M. Legendre 4 . For all this and for the Stirling contributions, see 26 . The above integral relationship is true if c > a > 0 and for |x| < 1, even if this bound can be discarded due to analytic continuation.
International Journal of Mathematics and Mathematical Sciences 5 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 F n D of n ∈ N variables and n 2 parameters , see 27, 28 , defined as F n D a, b 1 , . . . , b n ; c; x 1 , . . . , x n : with the hypergeometric series usual convergence requirements allowing the analytic continuation to C n deprived of the cartesian n-dimensional product of the interval 1, ∞ with itself.

Bent Slender Rod: A Large Deflections 2D-Model
Let us tackle a L-long, thin, and weightless rod whose end A is clamped and B free. We put at B the origin O of a x, y cartesian reference frame with x along the undeformed rod, from B to A; and y downwards, normal at B to x. Our basic assumptions are A 1 the rod is thin, initially straight, homogeneous, with uniform cross-section and uniform flexural stiffness EJ, where E is the Young modulus, and J the crosssection moment of inertia about a "neutral" axis normal to the plane of bending and passing through the central line; A 2 the slender rod is charged by coplanar dead loads; A 3 linear constitutive load: the induced curvature is proportional via 1/ EJ to the bending moment. The rod undergoes large deflections but small strains; A 4 the shear transverse deformation is ignored, and the rod is assumed to keep constant its previous length inextensibility ; A 5 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 EJ and bending moment M to the consequent curvature κ of the deflected centerline ±EJκ M x . With our traditional reference frame x, y , the curvature has always opposite sign to the bending moment EJ y xx to be solved together to the perfect clamping boundary conditions Throughout all the paper, it will be minded that the Cauchy problem can be solved in such a way. Putting In the sequel, we are going to use dimensionless variables ξ x/L, η y/L, so that the previous relationship becomes whose solution is 2.9

Rod Inflected by a Constant Bending Couple at Its Tip
The clamped section is at x L, see Figure 1. In such a case, the equation to be integrated is 2.5 , with M x −M 0 , and 2.9 provides, in nondimensional form, where χ M 0 L/EJ > 0. The elastica is more readable reverting at a dimensional form half circle of centre C L, Lχ and radius Lχ, see Figure 1, as it is well known.

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. Let us consider at Figure 2 a horizontal, flexible, homogeneous, heavy and cantilever; it will be inflected by the uniform evenly continuous load q weight for unity span length . Its elastica, according to the approximate beam theory, would consist of a fourth-order rational function of x. Being we obtain, in nondimensional coordinates, using 2.9 η ξ ν 3 and finally, passing to the variable λ σ/ ξ 3 − 1 : 14 By identification, we obtain the solution to our elastica through a Lauricella function of 3 variables x 1

2.15
For being compliant with the IRT requirements, it shall be | x k | < 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 |η ξ | < h 2 , h ∈ 0, 1 , we get i if 0 < ν < 0.0099, the first-order approximation can be accepted; ii if 0.0099 < ν < 1 by IRT , the hypergeometric function F 3 D 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 F 3 D 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.

A Tip-Sheared Horizontal Rod
We refer to Figure 3, where the force R is transversal in 29 is studied the inclined end load but not a follower one; we have M x Rx, so that, after defining the nondimensional μ RL 2 /2EJ, we obtain being the bending moment always greater than zero, so will be H x , and then it will be η ξ > 0. For this load too, the hypergeometric treatment is possible. Passing from ζ to λ, with

International Journal of Mathematics and Mathematical Sciences
We can see that, recognizing the Lauricella parameters a 2, b i 1/2, c 3 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 ϑ, s they integrate a first time for getting ϑ as a function of s. But dx ds cos ϑ s and dy ds sin ϑ s , so that they get, through a second integration, the elastica's x and y coordinates parametrized on s. Only few develop a different x, y 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 5 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 t 1 − ζ 2 , obtaining For the reader convenience, let us fix some notation. Let a > b > x > y > c. Consider the definite elliptic integral Combining formulae given in 25 : integral 234.16 page 76 and integral 340.01 page 205, introducing we obtain: where: Formula 2.25 provides the dimensionless deflection η to the position ξ and to the load μ: its 3D-plot, via Mathematica, as shown in Figure 4. 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 6 . 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' 7 comparison shows at the low loads an agreement which is becoming worst for increasing loads when the approximate theory is progressively failing.
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 where ϕ μ is the quantity introduced in 2.27 International Journal of Mathematics and Mathematical Sciences 13 Now we expand the elliptic integrals of I and III kind with respect to μ, with initial point

2.30
arriving at Moreover recalling the Maclaurin first-order expansion with respect to μ adding 2.31 and 2.32 , we obtain eventually 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 μ RL 2 / 2EJ 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 n 1. Other curves are provided referring to the first dynamic mode, n 2. Furthermore, the module k of the Jacobi elliptic functions involved there, is on its own depending on the tip load P via an equation to be solved previously. Nevertheless the essential! relationship k P keeps unrevealed, because the authors at page 740 write: The modulus of elliptic functions k and the parameter F 1 play the role of integration constants and they are related to the force P and the angle ϕ 0 by boundary conditions for each case of rod bending.
In fact, they do not provide the k computation at all, marking their plot solutions not in terms of P as it should be done , but of a sequence of k values 0.5; 0.6; 0.85; 1 − 10 −5 see Figure 2, page 742 of 33 , so that their x, y solution seems to be not easily accessible from P .

The Rod Inflected by Hydrostatic Pressure
We are looking for the shape of a L, EJ 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 y axis downwards and x over the water surface. The hydrostatic head, being ρ the liquid mass density and g the gravity, will have a linearly deep-growing profile Stevin , so that, integrating twice such a head, one arrives at the bending moment M x 1 6 ρgLx 3 .

2.34
By the usual identification, we obtain the solution to our elastica through the Lauricella function of three variables where σ ρgL 5 /24EJ, so that

2.36
For being compliant with the IRT requirements, it will be |x k | < 1, and then σ < 1.

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 y x . But κ, which is the consequence on the rod, is linked to its cause, namely, the bending moment M, by a Linear Constitutive Law κ ±M/EJ, being EJ 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 q x be the continuous transverse load whose highest value is Q over the unit of an L-length flexible cantilever q x Q sin π L x .

2.37
Integrating twice such a load, one arrives at the bending moment affecting the x-section In such a way, introducing ε L 3 Q/ π 3 EJ > 0, with 8 ε 1, following 2.9 , we find η ξ −ε

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 O ξ 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 O has been mapped to. Rod's behavior under strain can be modeled in several ways. For instance in 9 after integrating  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. 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 9 . 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 f, the approximate beam theory would answer f η 0 , but with large deflection, we have first to detect ξ Ω what is done assuming the unextensibility 1 ξ Ω 1 η 2 ξ dξ 1.

3.1
We see that ξ Ω is the sole unknown, having been previously found η ξ , so that, after ξ Ω , we will compute f η ξ Ω . let us provide the full treatment for such a case. By plugging

3.3
Our unknown appears in the upper integration limit; in order to solve 3.3 , we need to establish an existence-uniqueness Lemma.

International Journal of Mathematics and Mathematical Sciences
Defining and such a root is given by where sn u, k is the Jacobi sine amplitude of module k.
Proof. Let a > b ≥ x > y ≥ c. Observe that the function 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 with the ϕ and k meanings previously introduced. The very last step follows from the inversion: F arcsin ψ, k α ⇐⇒ ψ sn α, k , 3.10 which delivers x from 3.9 , after equating the right hand side to 2μ. We can then use the Lemma to solve 3.3 via the identification a b x y c

3.11
where it has been assumed 10 μ < 1. Minding that a 1/μ and c −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

3.12
Therefore, increasing R intensities 11 , and then μ, the successive positions of the free tip, that is, the components of the displacement − −−− → OΩ k ξ Ω , |η ξ Ω | will be described as functions of μ at Figures 9 and 10.
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.

Conclusions
Five elasticae have been analytically solved of a thin cantilever rod bent by different coplanar "dead" loadings 12 , 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.